Boundary element method (BEM) applied to the rough surface contact vs. BEM in computational mechanics

In the numerical study of rough surfaces in contact problem, the flexible body beneath the roughness is commonly assumed as a half-space or a half-plane. The surface displacement on the boundary, the displacement components and state of stress inside the half-space can be determined through the convolution of the traction and the corresponding influence function in a closed-form. The influence function is often represented by the Boussinesq-Cerruti solution and the Flamant solution for three-dimensional elasticity and plane strain/stress, respectively. In this study, we rigorously show that any numerical model using the above mentioned half-space solution is a special form of the boundary element method (BEM). The boundary integral equations (BIEs) in the BEM is simplified to the Flamant solution when the domain is strictly a half-plane for the plane strain/stress condition. Similarly, the BIE is degraded to the Boussinesq-Cerruti solution if the domain is strictly a half-space. Therefore, the numerical models utilizing these closed-form influence functions are the special BEM where the domain is a half-space (or a half-plane). This analytical work sheds some light on how to accurately simulate the non-half-space contact problem using the BEM.


Introduction
Boundary element method (BEM) is a numerical method used to approximate the solutions of mechanics problems using the boundary integral equation (BIE). It was first used in a paper by Brebbia and Dominguez [1] in 1977. For a linear elastostatic problem, it leads to an integral equation of the tractions and surface displacements over the boundary. Therefore, this method is also known as the boundary integral equation (BIE) method. The development of the computer technology makes BEM a powerful numerical methods and a general-purpose solver for different problems with arbitrary boundary. In the elastostatic problem, the BEM can be described as the Betti's reciprocal theorem (or the Somigliana identity) with the Kelvin's solution (or Mindlin's solution) as the auxiliary solution. By discretizing the boundary, performing the numerical integration over each boundary element and assembling the solution matrix, the unknown (either the surface displacement or traction per node) can be solved numerically from the system of linear equations. The main advantage of BEM for the elastostatic problem is that only the boundary needs to be discretized. The entire domain does not need to be discretized regardless of its size, if the body load is neglected. Additionally, the evaluation of the internal values inside the domain is exact and does not rely on the nodal density of the domain.
The BEM is first applied to the contact mechanics problem by Anderson et al. [2] and later on widely used in many other contact mechanics problems [3−6].
Since the conventional BEM 1 results in a dense and unsymmetric matrix which requires 2 ( ) O N memory and 3 ( ) O N operations for solving the system of linear equations using the Gauss elimination [7] where N is the number of nodes. For a reasonable computational time, a conventional BEM may be applied to the problem with 10 5 nodes or less. This is the main reason why the contact interfaces in those work are represented by (piece-wise) simple curves and discretized with few points. In those work, most of the contact domains are finite and two-dimensional. The spatial (three-dimensional) contact needs more boundary nodes and thus longer computational time. Therefore, the axisymmetric domain is commonly studied for the sake of low computational time. This disadvantage may be overcome by the fast multipole BEM [7].
Long before the birth of this numerical method, its fundamental, i.e., the boundary integral equation, has been intensively used by Muskhelishvili [8], Galin [9], Timoshenko and Goodier [10], Gladwell [11], Johnson [12], etc., in the contact mechanics to find the analytical solutions. This is due to the nature of the contact mechanics where most of the interesting solutions (e.g., contact pressure, contact area, surface displacement, etc.,) are on the contact interface (i.e., the boundary). In nearly all those work, the domain is assumed to be a half-space or a half-plane. This is because that the solid contact between bodies is often a local phenomenon so that the dimension of the contact area is negligibly small compared to that of the contacting bodies. There are attempts working on the non-halfspace domains (e.g., quarter space [13−16] and thin layer [17]). For an arbitrary boundary, the finite element method (FEM) may be used to calculate the corresponding influence coefficients [18].
Ever since the pioneering work of Conry and Seireg [19] and Kalker and Van Randen [20], a series of similar numerical models [21−49] have been proposed and are dominant in the simulation of rough surface contact. Those models all adopt the same assumption that the flexible bodies of the contact pair are half-spaces or half-planes. For a non-periodic problem, the surface displacement components can be modeled as a convolution of the tractions and the corresponding influence functions. For a non-periodic continuum half-space, the influence function is either the Boussinesq-Cerruti solution (three dimensional elasticity) [42] or Flamant (plane strain/stress condition) [21] solution. For a periodic continuum problem, the influence function is either the Westergarrd solution (plane strain/stress condition) or the Tripp solution [50] (three dimensional elasticity). If the micro-structure of the half-space is considered (e.g., the half-space is composed of lattices), the lattice Green's function can be used as the influence function [41]. In the discretized domain, the surface displacement vector can be calculated by the traction vector multiplied by a discretized influence matrix [42]. The normal boundary condition on the interface follows the Signori inequality (also known as the Kuhn-Tucker condition). Tangential boundary conditions varies with problems and the common ones are full stick, perfect slip and partial slip. These types of numerical models are wildly applied to the purely normal contact [33], sliding contact [40], partial slip [42], rolling contact [47] and adhesive contact [44,48]. The unknown tractions (e.g., contact pressure and shear stress) are solved iteratively through enforcing the corresponding boundary conditions. This may be done by transforming the contact problem to an equivalent optimization problem and solved by the classic solvers (e.g., the conjugate gradient (CG) method) [37,46,47].
In the early days, this type of numerical model is associated with different names, e.g., finite surface element model [21], conventional deformation matrix method [28], moving grid method [28], FFT-based method [33], etc. Starting from the work of Peng and Bhushan [40], this type of numerical methods is commonly accepted as the boundary element method (BEM). The reason lies in the fact that the surface displacement is calculated by the integral equation over the boundary which has a similar manner as the BIE in BEM. Therefore, the numerical models applied to the rough surface contact mentioned above may all be treated as the special BEM where the domain is a half-space. However, a rigorous proof is missing. To distinguish these two models, the so called BEM applied to the rough surface contact is referred to as the special BEM. The other one is referred to as the general BEM. The logic behind this proof is that the boundary integral equation using the Flamant solution and Boussinesq-Cerruti solution as the influence function can be recovered from the general BEM [51] when the domain is a half-space and a half-plane, respectively. The general BEM has recently been applied to the smooth adhesive contact [6] and mixed lubrication where the sub-surface stress is evaluated [52]. Li [52] showed the difference between the sub-surface stress with and without the half-space assumption. As far as we know, the general BEM has rarely been applied to the rough surface contact [52]. This proof can be served as a complimentary material for the BEM community that the half-space problem can be further simplified using more straight-forward integral equation. In the mean time, this study introduces the general BEM to the tribologists to improve the accuracy of a non-half-space problem.
In Section 2, the general BEM of an elastostatic problem with a finite body is briefly discussed. In Section 3, the general BEM is shown to degrade to the special BEM with the Flamant solution as the influence function if the domain is a half-plane. Similarly, in Section 4, the general BEM is shown to degrade to the special BEM where the Boussinesq-Cerruti solution is the influence function if the domain is a half-space.

Boundary element method for a finite domain
BEM is a powerful numerical technique used in continuum mechanics, especially for the linear problems. For the elastostatic problem, BEM relies on the Somigliana identity and the fundamental solutions (e.g., the Kelvin and Mindlin solution). If the body load is neglected, only the boundary of the domain needs to be discretized. The internal values of displacement, stress and strain can be accurately determined by the corresponding integral equation. This is a major difference from the finite element method (FEM) which is another popular methodology in the rough surface contact problems [53,54] , inside the domain  (the boundary  is excluded) is expressed by the following boundary integral equation (BIE) in tensor notation [51,55]: which is also known as the Somigliana identity. The source point is . The kernels, ij u  and ij p  , are the auxiliary (fundamental) solutions of thedisplacement and the traction components at    in an elastic infinite domain due to an point load acting at x (it is commonly known as the Kelvin solution). The explicit forms of the kernels, ij u  and ij p  , for three-dimensional elasticity can be found in Eqs. (where the boundary  is excluded). Based on the strain-displacement relation and Hooke's law, the internal state of stress ij  inside  can be expressed in a similar manner by the following BIE [51,55]: where the kernels, ijk D  and ijk S  , are given in Eqs. (A11) and (A13).
For a discretized form of Eq. (2) associated with higher order elements (e.g., linear, quadratic, etc.), the nearly singular behavior of sub-surface stress occurs when x is close to the boundary. This is caused by the numerical error introduced by the conventional Gauss quadrature [56]. The distance transformation method [52,56] may be used to improve the accuracy of sub-surface stresses for these cases. Note that the stress in the left hand side of Eq. analytically. For a discretized form of Eq. (2), nonsingular behavior may also be held if the boundary is discretized only by constant elements. This is because integrals over constant elements have closed-form solutions [7]. When (1) can be deduced to the following BIE for the surface displacement: if Γ is of class 1 C (continuously or piece-wisely differentiable up to the first order). Otherwise, ( ) ij C x can be determined based on the fact that the rigid body displacement is excluded when the boundary is traction-free. After discretizing the boundary  , the above BIE, together with the boundary conditions, can be formed into a system of linear equations. For more details on numerically solving BIE in Eq. (3), readers should refer to the text books [7,51,55].
The kernel ij p  has a strong singularity Thus, the boundary integrals with kernels, ij p  and ij u  , in Eq. (3) are in the sense of the Cauchy principal value, i.e., in Ref. [51], Consider a discretized boundry consisting of m line segments (elements): In the rough surface contact model, the boundary values are usually assumed to be uniform within each element. In BEM, this type of element is referred to as the constant element. Only one node needs to be placed within each element (commonly at the centroid). Therefore, the group of source nodes and field nodes are: be the traction and surface displacement vectors associated with each field point. Assuming  is of class 1 C , the discretized form of Eq. (3) is After combining the common terms in Eq. (4) and matrix inversion, we have where kl K is the influence sub-matrices (3 3)  representing the surface displacement vector at source point i x due to the unit traction vector acting on l  . When the domain is strictly smooth half-space and problem is purely normal contact, kl K should be the same as the wildly used Love's influence function [12].
Plane strain condition. A finite 3D domain is degraded to a 2D domain over the xy plane where the field and source points are x y x , respectively. Boundary integral equations in Eqs. (1), (2) and (3)

Boundary element method for a halfplane problem
In the plane strain condition, consider a half-plane with a flat boundary: , 0} x y x y , see Fig. 1. Equstions (1) and (3) [57]. This condition requires the vanishing of the state of stress at the far end which is exactly the criterion shown above. In order to guarantee the stress-free condition at the far end, the average shear and normal tractions over the entire  should be strictly zero. Defining a half-plane with the finite boundary: [ , ]; 0} x y x L L y , then the criterion of zero mean traction is: For the non-periodic contact where the traction is constrained in one or several finite strip, Eq. (7) is satisfied. If the average shear and normal tractions are not zero (it happens in the periodic contact problem), then the contribution of the average terms on the state of stress should be studied individually. Due to the simple geometry of the boundary  , the following simplification is available for the terms in the Kelvin solution (see Appendix A for more details): x r r y r r n y r Substituting the above simplification into Eqs. (A1) and (A2), we have the following simplified Kelvin solution: It is clear that all the kernels above are functions of x   and y.
Substituting the simplified Kelvin solution listed in Eq. (9) and 0 y  into Eq. (3), then it can be splitted into two BIEs 2 for the surface displacement components u and v : The above set of coupled integral equations can be decoupled using one-dimensional Fourier transform. The explicit form of the one-dimensional Fourier transform pair can be found in Eqs. (B1) and (B2) in Appendix B. Using the convolution law, Eqs. (10) and (11) in the frequency domain are 3 : is one-dimensional Dirac function. 3 Note that the constant term in 11 ( 0), i.e.,   ( ,0), ( ,0) P k U k and  21 ( ,0) P k are determined similarly. Finally, the closed-forms of U and V are obtained by solving the system of linear equations in Eqs. (12) and (13): After the inverse Fourier transform of Eqs. (14) and (15), u and v in the Cartesian coordinates are: Equations (16) and (17) are exactly the same as the integral forms of the surface displacement using the Flamant solution as the influence functions [12], except for the missing constant terms regarding the rigid body displacement.
Substituting Eqs. (14) and (15) into Eq. (1) in the frequency domain, the displacement components ( , ) U k y and ( , ) V k y in the frequency domain are: Sneddon [58] revisited the half-plane problem using the Airy's function and the Fourier transforms. Equations (18) and (19) are the same as those derived by Sneddon [58] if the consistent form of Fourier transform pair is used. Setting 0 y  in Eqs. (18) and (19), Eqs. (14) and (15) are recovered. After the inverse Fourier transform, the integral forms of the surface displacement using the Flamant solution as the influence function in the Cartesian coordinates are recovered [12,58].
Applying the Fourier transform to the Hooke's law for the plane strain condition, we have [58] ( , ) 2 ( , ) xy k y G U k y i kV k y y (22) Substituting the forms of ( , ) U k y and ( , ) V k y given by Eqs. (18) and (19) into Eqs. (20), (21) and (22) x y z x L L y L L z , the zero traction criterion can be stated as: A half-space assumption results in the following simplifications of the terms in the Kelvin solutions (see Appendix A for more details): The information on the two-dimensional Fourier transform pair used in this study can be found in Appendix B. A Fourier transform table is given in Table B1 in Appendix B for the help of deriving Eqs. (28), (29) and (30). Solving the system of linear equations above, we can get the surface displacement components due to the shear traction, x q , in the frequency domain: Using the inverse Fourier transform with the help of Table B1 in Appendix B and the convolution law, the well-known solutions of i u of the Cerruti problem are recovered [12]: The results due to the shear traction, y q , are the same as above if u and v are interchanged and x q is replaced by y q . Let 0 p  and 0 x y q q   , then the corresponding boundary value problem is referred to as the Boussinesq problem. The Fourier transform of Eq. (3) results in Solving the above system of linear equations, we have the closed form solution of i u due to the normal traction, p in the frequency domain: Using the inverse Fourier transform with the help of Table B1 in Appendix B and the convolution law, the well-known solutions of i u of the Boussinesq problem are recovered [12]: Following the procedures shown in Section 3, the displacement components and the state of stresses inside the domain can be also obtained in the same manner.
In summary, the integral equations using the Boussinesq (Cerruti) solution as the influence functions are a special form of the BEMs in Eqs. (1−3) where the domain is strictly a half-space.

Discussions
In this paper, we do not tend to judge how accurate is the special BEM using the Boussinesq-Cerruti/Flamant solutions. The ultimate goal of this theoretical study is to show the relationship between the boundary element method in the community of computational mechanics and the one used in rough surface contact (or more broadly, Tribology). In the previous tribology literatures, this "common sense" has been taken so granted that it was never proved rigorously. The researchers from the community of boundary element method seldom applied BEM to the half-space/ quarter or other infinite domains. This history causes a misunderstanding of BEM in many tribology literatures that it can only be applied to the halfspace/quarter-space since the corresponding Green's function is available. This paper shows that BEM can be applied not only to the half-space problem, but also to other domains with arbitrary boundary.
In Sections 3 and 4, it is rigorously proved that the so-called BEMs using the Boussinesq-Cerruti (or Flamant) solutions as the influence functions are a special form of the BEM where the domain is strictly a half-space (or half-plane). For the 2D plane strain/ stress condition, the Flamant solution is commonly obtained using either the complex variable method [8,59], the stress function in the polar coordinate [10] or the Airy's function [58]. For the 3D elasticity, the methods include the Boussinesq potentials [12] or the Papkovich-Neuber potentials [60]. This study offers an alternative approach to find the Flamant and Boussinesq-Cerruti solution using the boundary element method.
Therefore, the application of those special BEMs should be restricted to the contact scheme where half-spaces/half-planes are in contact, e.g., the half-space indented by a rigid indenter. However, except for this special contact scenario, there are many schemes where the flexible contact bodies are not exactly half-spaces (e.g., a half-space with a nominally flat rough boundary indented by a rigid flat).
For those cases, the integral equations using the Boussinesq-Cerruti (or Flamant) solutions should be replaced by the BIEs shown in Eqs. (1−3) in order to account for the non-half-space domain. The contact problem is still equivalent to an optimization problem where the classic solver (e.g., Conjugate Gradient method) can still be applied to solve the traction, i p , and the surface displacement, i u . Then, using Eqs. (1) and (2), the displacement and the state of stress inside the domain can be accurately determined. One reason why the Boussinesq-Cerruti (or Flamant) solutions are wildly used on the rough surface contact problems is that the fast algorithm (e.g., the Fast Fourier transform (FFT) [33] and Multi-Level-Multi-Integration (MLMI) [26]) can be applied to accelerate the numerical integration. This benefit no longer exists for Eqs. (1−3). Since the integrals in Eqs. (1−3) are not strictly convolutions, the Fourier transform can no longer be applied, along with the FFT. MLMI [61] significantly reduces the computational time by performing the numerical integration only on the coarse level (with few nodes on it). The final results are approximated through the interpolation between the coarser and finer level. However, it is assumed that the influence coeficient matrix is invariant on all levels. This is not strictly true for a non-half-space (plane) domain, especially when the rough surface exists on the boundary. Therefore, MLMI may not be valid for a non-half-space (plane) problem. This difficulty may be solved by the fast multi-pole method [7].
Recently, BEM is thoroughly discussed in two important papers [62,63]. The former one is by Müser et al. [62] and it is a summary of the results of a rough surface contact problem challenge initiated in 2016. Among all models, four numerical models are built under the framework of BEM. The Boussinesq solutions are adopted in all four models. Their differences mainly lie in (1) how adhesion is introduced and (2) how to achieve a stable iterative process towards a convergent solution. All the models have good agreement with the experimental results and other methods. This conclusion is valid within the scope that the roughness is relatively smooth.
The latter one [63] is a review paper on the simulation technique in tribology and part of Section 2.2 is dedicated to a brief review of the development of BEM. From the authors' perspective, the BEM mentioned in this review paper is all about the special BEM using the Boussinesq/Cerruti solutions (or other equivalent solutions in the plane strain/stress). The limitations (e.g., large deformation, large sliding, plastic deformation, half-space assumption, etc.) of BEM have been discussed. One of limitations is stated as follows: "Such solutions exist for a limited number of cases and mainly under the assumption that the solid can be locally considered as a flat half-space. These limitations imply a more restrictive field of application for the BEM compared to the FEM, which is a versatile numerical method." In this study, we show theoretically that, at least in the linear elastic problem, BEM is also a versatile numerical method for non half-space problem.

Conclusions
In this study, the so-called "BEM" applied to the rough surface contact is rigorously proved to be the special case of the general BEM where the domain is a half-space (or half-plane for plane strain/stress condition). For the plane strain condition, the BIEs in the general BEM results in the well-known Flamant solution. For the three-dimensional elasticity, the same BIE for the relation between the surface displacement and the traction becomes the well-known Bousssinesq (Cerruti) solution if the normal (shear) traction is applied on the boundary alone. The general BEM offers an alternative way to simulate those contact scenarios with a better accuracy where the domain is a nonhalf-space.
A general way of evaluating the Fourier transform and its inverse is using the residual theorem [64]. One-dimensional Fourier transforms frequently used in Section 3 are shown in Table B1.
The two-dimensional Fourier transform pair used in this study is The modulus of ( , ) x y k k is 2 2 x y k k k   . Some Fourier transform is evaluated based on the fact that the Fourier transform of the axisymmetric function is the same as the Hankel transform of zero order [58].
Others are evaluated based the some known solution and the derivative law. Some two-dimensional Fourier transforms frequently used in Section 4 are prepared in Table B1.