# Single Level Fast Multipole Method for frictionless rough contact problems

Claudia Stiebritz, Hai-Ping Yin, Julien Cesbron

► To cite this version:

Claudia Stiebritz, Hai-Ping Yin, Julien Cesbron. Single Level Fast Multipole Method for frictionless rough contact problems. International Journal of Mechanical Sciences, 2025, 286, pp.109810. 10.1016/j.ijmecsci.2024.109810. hal-05210230

HAL Id: hal-05210230

https://univ-eiffel.hal.science/hal-05210230v1

Submitted on 14 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Single Level Fast Multipole Method for frictionless rough contact problems

Claudia Stiebritz<sup>a,*</sup>, Hai-Ping YIN<sup>b</sup>, Julien Cesbron<sup>a</sup>

## ARTICLE INFO

### Keywords:

Contact mechanics

Surface roughness

Boundary element method

Boussinesq's contact problem

Fast Multipole Method

## ABSTRACT

A perfectly smooth contact surface does not exist in nature and industrial applications. Even a body, that seems perfectly smooth to the naked eye, will show surface roughness at a higher magnification. Due to the roughness of the surface, there are areas of contact and separation, which increases the complexity of the contact calculation. However, this computational complexity increases further due to the multi-scale nature of road surface texture and large contact area in the case of tyre/road interaction. To reduce the computational complexity of this contact problem, the Single Level Fast Multipole Method (SLFMM) is developed within this paper. The contact problem is based on Boussinesq's contact theory and for the time being, the influence of friction and lateral displacement are neglected. To validate the accuracy and reduction in computational complexity, the SLFMM was applied to rough surfaces of different complexities and compared to a reference method, the so-called Matrix Inversion Method (MIM). Results indicate that the new method computes the pressure distribution and displacement accurately, with a global error of less than 1%. The advantage of the new method compared to the MIM is the multipole expansion, which clusters adjacent contact points to a single center point. As a result, the computational complexity of the contact calculation is reduced. Overall, the Single Level FMM computes the results faster than the reference method. These results demonstrate that the Fast Multipole Method meets the requirements of accuracy and accelerated computation for rough contact problems.

## 1. Introduction

Rough contact problems are of particular interest since contact surfaces are typically not topographically smooth. It is even assumed that surface roughness occurs on different magnification scales [1,2]. Consequently, there is no continuous pressure distribution over the entire contact surface, but rather areas of contact and separation [3,4]. Rough contact problems can be found in numerous engineering applications, such as the contact between bearings [5], electrical contact of relays [6], electrical contact resistance in proton exchange membrane fuel cells [7], the contact between artificial joints [8,9], the rail/wheel contact [10,11] or tyre/road contact [12–15].

Particularly, Hertz and Boussinesq are associated with classical contact mechanics. Hertz (1881) accentuates, that the contact between two elastic spheres induces a slight deformation of both, within their contact area. As a result, the vertical contact pressure is evenly distributed over the contact area [16]. Boussinesq (1885) was the first to study the normal contact between a rigid surface and an elastic half-space, by applying the potential theory [17]. Later, Love (1929) and as well Sneddon (1948) derived from Boussinesq theory the analytical result for a single axisymmetrical asperity such as a cylinder or a cone [18,19]. Nevertheless, the computational complexity for a

numerical method increases dramatically ($O(N^3)$) when studying the tyre/road contact problem, due to the time variation of contact points and the relatively large contact surface.

To solve Boussinesq's contact problem, several methods were applied, like the Matrix Inversion Method (MIM), developed by Johnson [4], the Cattaneo method focusing on plane problems [4] and its enhancements [20–22]. For smooth surfaces the Kalker method [23] and for contact problems with elastic solids of dissimilar materials the Chen method [24] were developed. Another approach is the Two-Scale iterative Method (TIM) for rough contact problems, developed by Cesbron et al. [13,25] and later Dubois et al. [14]. Further methods to compute rough contact problems are the Cholesky factorization used by Webster and Sayles [26], the multigrid method based on Brandt and Lubrecht [27] and Fast Fourier Transform (FTT) [28–32]. Polonsky and Keer [28,29] demonstrated that the FFT is in some cases slower and sometimes even less accurate for rough contact problems. The reason for this is, that the FFT requires for each contact point a minimum of four mesh points.

The motivation and originality of the present study is to use the Fast Multipole Method (FMM) to reduce the complexity and improve the

<footnote><sup>a</sup> Univ Gustave Eiffel, Cerema, UMRAE, Bouguenais, 44344, France</footnote>

<footnote><sup>b</sup> Laboratoire Navier, École Nationale des Ponts et Chaussées, Champs-sur-Marne, 77420, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: claudia.stiebritz@univ-eiffel.fr (C. Stiebritz), honore.yin@enpc.fr (H.-P. YIN), julien.cesbron@univ-eiffel.fr (J. Cesbron).</footnote>
efficiency of the computation of rough contact problems. Indeed, the FMM was developed by Rokhlin and Greengard (1980s) to speed up the calculation of N-body interactions, by reducing their complexity [33–35]. In the early days, the FMM was mainly applied to potential field problems, but its advantages for different physical applications were soon discovered. Its efficiency in solving Matrix-Vector Multipli-cations was applied to BEMs [36,37] and their various applications. Amongst others the FMM was used to solve large-scale 2D elastostatic problems [38,39], 2D elastic wave surface integral equations [40], 3D elasticity and fracture problems [41,42], 3D elastodynamics in the time [43] or frequency [44] domains, time domain partial ele-ment equivalent circuit analysis [45], 3D potential problems [46,47], 3D half-space acoustic wave problems [48], poroelastodynamic prob-lems [49], electromagnetic scattering problems [50,51] and anisotropic material problems [52]. Regarding the Fast Multipole Method, two approaches have to be mentioned, the Single and Multi Level FMM. The Single Level approach, recursively divides the domain into 2 levels, resulting in 16 leaf cells. In comparison, the Multi Level FMM recursively divides the domain further than level 2, which leads to a minimum of 64 leaves.

Therefore, this paper introduces a Single Level Fast Multipole Method (SLFMM) to rough contact problems, focusing mainly on Boussinesq's theory. For this purpose, Section 2 presents Boussinesq's contact problem, the multipole expansion, as a key element of the Fast Multipole Method and the SLFMM itself. Section 3 illustrates and discusses the results obtained for the MIM and Single Level approach for surfaces of different complexity (single half-sphere/hole, 25 half-spheres/holes and a real road surface). The final Section 4 outlines the key findings of the new approach.

## 2. Theoretical background

This chapter explores Boussinesq's contact problem combined with the Fast Multipole Method to solve rough contact problems. Boussinesq theory was chosen to exclude friction at the current state, to focus on an accurate computation of the normal contact pressure distribution.

### 2.1. Boussinesq contact problem

To calculate the pressure distribution between two bodies the states of contact and separation have to be distinguished [3]. This phenomenon is called the “law of unilateral contact” and is described by:

$$ \forall M \in \Sigma_c, h(M) = 0 \quad \text{and} \quad p(M) > 0, \qquad (1) $$

$$ \forall M \in \bar{\Sigma}_c, h(M) > 0 \quad \text{and} \quad p(M) = 0, \qquad (2) $$

$$ (\text{with: } h(M) = u_z(M) - z_r(M) + \delta + z(M)), \qquad (3) $$

where Eq. (1) defines the state of contact ($\Sigma_c$) and Eq. (2) separation ($\Sigma_c$bar). Here, $h(M)$ defines the gap function and $p(M)$ the normal contact pressure. The gap function (Eq. (3)) in general is determined by the distance between the half-space and the rough surface, where $\delta$ is the global penetration between the two bodies, $u_z(M)$ the displacement, $z_r(M)$ and $z(M)$ the heights of the rigid surface and the half-space, respectively [3,4].

Boussinesq contact problem (Fig. 1) describes a concentrated force acting on the surface of an elastic half-space. Accordingly, a normal pressure $p(\xi, \eta)$ acting in a point $S(\xi, \eta)$, within the contact area $\Sigma_c$. Now it is of interest, how this force influences the displacement and stress field in any point $M(x, y)$, at the surface of the half-space $\Sigma$, with $\Sigma = \Sigma_c \cup \bar{\Sigma}_c$ and $\Sigma_c \cap \bar{\Sigma}_c = \emptyset$ [14]. Within Boussinesq contact problem the displacement $u(M)$ is defined as the integral of the influence function $T(M, S)$ multiplied by the normal contact pressure $p(S)$. The influence function calculates the interaction between the points $M(x, y)$ and $S(\xi, \eta)$:

$$ \forall M \in \Sigma, u(M) = \int_{\Sigma_c} T(M, S)p(S)dS, \qquad (4) $$

<figure><img src="image_4.png" /><figcaption>Fig. 1. Boussinesq's contact problem is defined by a concentrated force $p(\xi, \eta)$ acting on the surface of an elastic half-space. This leads to a displacement at the point $M(x, y)$, at the surface of the half-space $\Sigma$.</figcaption></figure>

$$ \text{with: } \forall (M, S) \in \Sigma^2, T(M, S) = \frac{1 - v^2}{\pi Er}, \qquad (5) $$

where $E$ and $\nu$ are the Young modulus and Poisson coefficient of the half-space and $r$ is described by $\sqrt{(x - \xi)^2 + (y - \eta)^2}$. Provided that the global penetration $\delta$ is known, the contact area $\Sigma_c$ and pressure distribution $p(M)$ can be determined. Therefore, the contact area must be discretized before the Boussinesq contact problem can be solved. The surface of the half-space is partitioned into $n$ identical rectangular elements (Fig. 2), with a length $h_x$ and height $h_y$ and its center at $M_i(x_i, y_i)$, where $i \in [1, n]$ [53]. According to the discretization, Boussinesq's contact problem is characterized by:

* $\mathbf{p} = \{p_i\}_{i \in [1, n]}^T$: vector of contact pressure, where $p_i$ is uniform on each element $i$,

* $\mathbf{u} = \{u_i\}_{i \in [1, n]}^T$: vector of displacement.

Consequently, the general relation between the contact pressure distribution $\mathbf{p}$ and the displacement $\mathbf{u}$ is described by $\mathbf{A}\mathbf{p} = \mathbf{u}$, where the influence matrix $\mathbf{A}_{ij}$ is defined by $(\forall (i, j) \in [1, n]^2)$:

$$ A_{ij} = \frac{1 - v^2}{\pi E} \int_{y_j - \frac{h_y}{2}}^{y_j + \frac{h_y}{2}} \int_{x_j - \frac{h_x}{2}}^{x_j + \frac{h_x}{2}} T(M(x, y), S(\xi, \eta)) d\xi d\eta. \qquad (6) $$

Due to the direct Matrix-Vector Multiplication, Boussinesq's contact problem for large contact areas is computationally expensive.

<figure><img src="image_5.png" /><figcaption>Fig. 2. Discretization of the elastic half-space surface in $n$ identical rectangular elements (length $h_x$, height $h_y$ and center at $M_i(x_i, y_i)$, where $i \in [1, n]$).</figcaption></figure>
## 2.2. Multipole expansion

Before the FMM can be applied to Boussinesq's contact problem, its influence function $T(M, S)$ has to be approximated. Therefore, the points $M(x, y)$ and $S(\xi, \eta)$ have to be defined by the spherical coordinates $(R, \theta, \phi)$ and $(\rho, \alpha, \beta)$, respectively [54]. Appendix A provides a detailed definition of the mathematical derivation of the regular and irregular solid harmonics, given by:

$$R_{n,m}(\mathbf{OS}) = \frac{1}{(n+m)!} P_n^m(\cos \alpha) e^{im\beta} \rho^n, \quad (7)$$

$$\overline{S_{n,m}(\mathbf{OM})} = (n-m)! P_n^m(\cos \theta) e^{im\phi} \frac{1}{R^{n+1}}, \quad (8)$$

where the indices $n$ and $m$ describe the degree and order of the solid harmonic. The complex conjugate is denoted by an overline like $\overline{S_{n,m}}$. In addition, the distance between the points and the origin is defined by $(|\mathbf{OS}| < |\mathbf{OM}|)$. Then the inverse of the distance $1/r$ from Boussinesq's contact problem can be approximated by solid harmonics:

$$u(M) = \frac{1 - v^2}{\pi E} \int_{\Sigma} \sum_{n=0}^{nL} \sum_{m=-n}^{n} \overline{S_{n,m}(\mathbf{OM})} R_{n,m}(\mathbf{OS}) \rho(S) d\Sigma, \quad (9)$$

which is called the multipole expansion [55]. To determine Boussinesq contact problem, it is necessary to estimate the recurrence relations of the solid harmonics (cf. Appendix B).

## 2.3. Single level fast multipole method

To begin with, the studied domain has to be divided recursively into 2 levels by a quadtree to allocate the grid points into its leaf cells (Fig. 3). Therefore, the entire domain is defined as root cell zero, level 1 as its four child cells, and level 2 as its 16 grandchildren or leaves. Since the Fast Multipole Method requires for the far-field approximation only cells at level ≥ 2, this approach is called Single Level FMM. In addition, the near-field and far-field need to be defined. The near-field consists of Cell C plus its adjacent cells and is computed directly. In contrast, the influence of the far-field cells has to be approximated by the FMM approach. Accordingly, the advantage of the FMM is the reduction of complexity by limiting the total count of direct computation.

<figure><img src="image_6.png" /><figcaption>Fig. 3. Schematic of the Single Level FMM (16 cells or leaves) consisting of 2 levels with exemplary illustration of Cell C (blue), adjacent (yellow) and far-field (white) cells.</figcaption></figure>

To calculate the far-field contribution, the multipole moments of each of the 16 leaf cells must be computed. Therefore, the multipole moments approximate the influence of points within a leaf to its cell center $S_c$ [39,56]. The multipole moments can be derived from the multipole expansion (Eq. (9)) and are written as:

$$M_{n,m}(S_c) = \int_{\Sigma_l} R_{n,m}(S_c S) \rho(S) d\Sigma. \quad (10)$$

To estimate, how the far field influences the local displacement within Cell C, the moments of the far-field have to be translated to its center $M_c$. This translation is the so-called moment-to-local translation (M2L)

$L_{n,m}$ [39], written as:

$$L_{n,m}(M_c) = (-1)^n \sum_{n'=0}^{nL} \sum_{m'=0}^{n'} M_{n',m'}(S_c) S_{n+n',m+m'}(S_c M_c). \quad (11)$$

During the last step, the result obtained by the M2L (Eq. (11)) is used to calculate the local displacement $u_{k, far}$. Therefore, the far-field influence is projected to every field point $M$ within Cell C [34,54], given by the local expansion:

$$u_{k, far}(M) = \frac{1 - v^2}{\pi E} \sum_{n=0}^{nL} \sum_{m=-n}^{n} R_{n,m}(M_c M) L_{n,m}(M_c). \quad (12)$$

With the estimated far-field contribution, the total local displacement $u_{k,z}(M)$ can be determined by summing up the near-field $u_{k,near}(M)$ and far-field displacement $u_{k, far}(M)$. In the final step, the local pressure distribution $p_k$ is updated by a local matrix inversion based on a prediction/correction algorithm [53]. Based on this procedure, the local pressure distribution for all leaves has to be updated. After the pressure distribution is entirely updated, the algorithm iterates until the solution converges within a prespecified tolerance $\epsilon$ [54], written as:

$$\frac{\| \mathbf{p}^{m_c+1} - \mathbf{p}^{m_c} \|}{\| \mathbf{p}^{m_c} \|} \le \epsilon \quad \text{with} \quad \| \mathbf{p} \| = \sum_{i=1}^{n} |p_i|^2. \quad (13)$$

Fig. 4, shows the iterative approach to estimate the local pressure distribution with a Single Level FMM. At first, the local influence matrix

<figure><img src="image_7.png" /><figcaption>Fig. 4. Flowchart for the iterative Single Level FMM to determine the contact pressure distribution.</figcaption></figure>
$A_{k,loc}$ for each leaf cell is computed and stored by Eq. (6). This process is computationally expensive since each grid point is considered, which results in a dense and symmetric matrix. Furthermore, the program consists of three loops, an external and two internal loops. The external loop controls the iterations ($It$) of the algorithm, and the first internal loop considers each leaf once per iteration as Cell C. With the last internal loop, the program verifies if the updated pressure distribution meets the predefined tolerance.

With the start of the first iteration, the multipole moments (Eq. (10)) for each leaf cell are determined. Since the multipole moments include the integration, they have to be re-estimated for each iteration. Then the first Cell C, its adjacent cells, and far-field are defined within the internal loop. The near-field displacement $u_{k,near}$ for the adjacent cells is calculated directly by matrix-vector multiplication. Afterwards, the far-field contribution has to be determined by the M2L-translation (Eq. (11)) and the local expansion (Eq. (12)). The total displacement $u_{k,z}$ and gap function (Eq. (3)) are determined and finally, the local pressure distribution $p_k^{m_c+1}$ for the first Cell C is updated. After that, the remaining 15 cells are determined, following the same principle. When the pressure distribution for all 16 leaves is updated, the program verifies, if the new pressure distribution converges within the prespecified tolerance $\epsilon$ (Eq. (13)). If it does not meet the defined tolerance, it continues with the next iteration, otherwise, the program stops. Beyond that, the algorithm stops when the maximum number of iterations $It_{max}$ is reached.

### 3. Results and discussion

In this chapter, the SLFMM is used to evaluate surfaces of different complexity, within Cartesian coordinates. To begin with, the multipole expansion regarding its accuracy for a simplified case has to be validated. According to that, the SLFMM is applied to less complex surfaces consisting of axisymmetrical asperities or holes. Finally, the SLFMM is applied to a complex rough surface. The SLFMM is compared to a reference method derived from the Matrix Inversion Method, which is based on a prediction/correction algorithm [53]. The relative error between the two approaches is determined by:

$$ \epsilon_g = 100 \frac{\| \mathbf{P}_{SLFMM} - \mathbf{P}_{MIM} \|}{\| \mathbf{P}_{MIM} \|}, \quad (14) $$

which is called the global error $\epsilon_g$. Both methods are implemented in Python using a computer with 32 GB RAM, 16 cores and Intel Core i9-11950H11th (2.60 GHz) as the CPU.

#### 3.1. Validation of the multipole expansion

To prove that the SLFMM can be used for rough contact problems, the multipole expansion has to be verified. Therefore, the multipole expansion is compared to the analytical results for a single rigid paraboloid of revolution (Fig. 5). The paraboloid is in contact with an elastic half-space by a known penetration $\delta$, resulting in a total load $P$. To verify the accuracy of the multipole expansion (Eq. (9)), it is compared to the analytical displacement outside the contact area. The analytical displacement $u(r)$ is calculated by the following expression,

<figure><img src="image_8.png" /><figcaption>Fig. 5. A simplified axisymmetrical asperity (paraboloid of revolution) in contact with an elastic half-space by a known penetration $\delta$.</figcaption></figure>

referring to Hertz contact theory [16].

$$ \forall r > a \quad u(r) = \frac{u_0}{\pi} \left( \left( 2 - \frac{r^2}{a^2} \right) \arcsin\left(\frac{a}{r}\right) + \frac{r}{a} \sqrt{1 - \frac{a^2}{r^2}} \right), \quad (15) $$

where $a$ defines the radius of the contact area and $r$ the radial distance. The multipole expansion, determined by Eq. (9), requires the recurrence relations for the solid harmonics, seen in Appendix B. To calculate the displacement, the analytical pressure distribution has to be applied for both approaches. To determine the analytical pressure distribution, the following equations are used:

$$ p(r) = \begin{cases} \frac{3p_{nl}}{2} \left(1 - \frac{r^2}{a^2}\right)^{\frac{1}{2}} & \forall r \le a, \\ 0 & \forall r > a, \end{cases} \quad (16) $$

where $p_m$ is the average contact pressure. Fig. 6 illustrates the analytical and the approximated normalized displacement $u/u_0$ (multipole expansion).

<figure><img src="image_9.png" /><figcaption>Fig. 6. Result for an axisymmetrical asperity: normalized displacement $u/u_0$ determined analytically and by the multipole expansion with a given contact pressure distribution $p_r$ (order 0 to 8).</figcaption></figure>

The approximated displacement was determined for order 2 to order 6. Since the symmetry of the revolution of the paraboloid compensates the moments for odd numbers, the result does not improve from order 0 to order 1 for example. This shows how the accuracy improves when the order number increases. Nonetheless, it can be observed that the multipole expansion shows an acceptable result for order 0 with a global error of 0.24% (Fig. 7). For order 8 the multipole expansion approximates the analytical displacement with a global error of 0.0007%. Consequently, the multipole expansion meets the accuracy requirement to compute the displacement within the contact area.

<figure><img src="image_10.png" /><figcaption>Fig. 7. Result for an axisymmetrical asperity: global error $\epsilon_g$ between the analytical and multipole expansion result (order 0 to 8).</figcaption></figure>
## 3.2. Surface generation

To evaluate the Single Level FMM surfaces of different complexity were generated and studied. The first set of surfaces consists either of half-spheres or holes (Fig. 8, b, c and d), which are numerically created.

<figure><img src="image_11.png" /><figcaption>Fig. 8. Numerical generated rigid surfaces: (a) 1 half-sphere (10 mm × 10 mm), (b) 25 half-spheres (50 mm × 50 mm), (c) 1 hole (4 mm × 4 mm) and (d) 25 holes (20 mm × 20 mm). The gray-scaled color bar on the right indicates the depth of the surface.</figcaption></figure>

For the surfaces with one half-sphere (Fig. 8a) and 25 half-spheres (Fig. 8b), the following expression was applied:

$$z_r(x, y) = \sqrt{R^2 - r_d} \quad \text{if } r_d \le a, \qquad (17)$$

$$\text{with: } r_d = (x - x_c)^2 + (y - y_c)^2, \qquad (18)$$

where a is equal to 5 mm and R is 10 mm. The coordinates $x_c$ and $y_c$ describe the position of the sphere center. Furthermore, the surfaces consisting of one hole (Fig. 8c) and 25 holes (Fig. 8d) were defined by:

$$z_r(x, y) = \begin{cases} -1 - \cos(\pi r_d / b) & \text{if } -b < r_d < b, \\ 0 & \text{if } r_d > b \text{ and } -b < x, y < b, \end{cases} \qquad (19)$$

where $r_d$ is equal Eq. (18) and b is 2. To generate the surface with the randomly placed half-spheres with radius of 4, 5 or 6 mm (Fig. 9), Eq. (17) was applied. Solely the information about the position and radius of each half-sphere were extracted from a scanned surface.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Numerical generated rigid surface: 24 randomly placed half-spheres (50 mm × 50 mm). The gray-scaled color bar on the right indicates the depth of the surface.</figcaption></figure>

Finally, a sample of real road surface with a total size of 200 mm × 191 mm was measured. This measurement was obtained by a micro-measurement station using an optical sensor (PH-Line 2100 from Image Optique in France). The horizontal displacement step of the sensor is 0.2 mm and the vertically measured values are between 0 and 30 mm with a resolution of 0.015 mm [25]. From the texture scan

the three-dimensional coordinate ($x_i, y_i, z_i$), the resolution, the length of the surface, and the number of grid points in x and y direction were extracted. Fig. 10 shows a 50 mm by 50 mm sample of this texture scan used for contact calculation in this study.

<figure><img src="image_13.png" /><figcaption>Fig. 10. 2D-measured rigid surface: rough surface (50 mm × 50 mm). The gray-scaled color bar on the right indicates the depth of the surface.</figcaption></figure>

## 3.3. Results single level FMM

To prove that the SLFMM can be applied to rough contact surfaces, four simplified and two more complex surface are evaluated. For this purpose, the main parameters are defined by:

* The elastic half-space (rubber block) is characterized by a Young Modulus of 6 MPa and a Poisson coefficient of 0.5.

* The penetration and the initial pressure distribution are $\delta = -1.0$ mm and $p^{m_c} = \{0, ..., 0\}^T$, if not stated differently.

* The convergence criteria $\epsilon$ is equal to $10^{-3}$, the maximum number of iterations is set to 20 and the order number is 4.

The proposed Single Level approach is applied to simplified surfaces of spherical asperities or holes and validated by comparing it to the MIM. To begin with, the pressure distribution and displacement for one half-sphere (Fig. 8a) and 25 half-spheres (Fig. 8b) are determined. The radius of the half-sphere is 10 mm and the surface consists of 10 000 grid points with a resolution of $h_x = h_y = 0.1$ mm. The pressure distribution was determined by order 0, 2 and 4 to demonstrate the approximation accuracy of the new method (Fig. 11). Order 0 approximated the result with a global error of 37.84%, which was drastically reduced to 1.96% and 0.13% for orders 2 and 4, respectively. Consequently, the SLFMM is congruent to the reference method, with a satisfying result for order 4 within 14 iterations. This statement can be verified by comparing the result for the displacement shown in Fig. 12. Here, the global error between the two methods is 0.005%, which proves the accuracy of the SLFMM.

<figure><img src="image_14.png" /><figcaption>Fig. 11. Result for 1 half-sphere: pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Result for 1 half-sphere: displacement for MIM and SLFMM on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>

For the surface with 25 half-spheres, each of the half-spheres has a radius of 10 mm. The Single Level FMM requires 10 iterations to converge within the pre-defined accuracy. Fig. 13 displays the contact pressure over the total surface calculated by the new approach. It illustrates, that the SLFMM considers the interactions well between the half-spheres. As a result, there is a larger contact area on the outer asperities than on its center asperity. Additionally, the pressure distribution for the center asperity is symmetric and asymmetric for the asperities in the corners.

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Result for 25 half-spheres: pressure distribution over the entire surface (SLFMM). The orange-scaled color bar on the right indicates the contact pressure value.</figcaption></figure>

To prove how accurately the SLFMM approximates the results, Fig. 14 provides a close-up of five asperities. Firstly, it is obvious, that there is an interaction between the asperities proven by the contact pressure differences between the inner and outer asperities.

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Result for 25 half-spheres: pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>

<caption>**Table 1**</caption>

Comparison of the CPU time and global error $\epsilon_g$ for the SLFMM (orders 0, 2 and 4) for 25 half-spheres.

<table><thead><tr><th>Order</th><th>Convergence</th><th>CPU time</th><th>Global error ε<sub>g</sub></th></tr></thead><tbody><tr><td>0</td><td>10<sup>-3</sup></td><td>37 s</td><td>5.92%</td></tr><tr><td rowspan="2">2</td><td>10<sup>-2</sup></td><td>24 s</td><td>0.28%</td></tr><tr><td>10<sup>-3</sup></td><td>38 s</td><td>0.10%</td></tr><tr><td rowspan="2">4</td><td>10<sup>-2</sup></td><td>31 s</td><td>0.15%</td></tr><tr><td>10<sup>-3</sup></td><td>50 s</td><td>0.02%</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Result for 25 half-spheres: displacement for MIM and SLFMM on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>

Foremost, both methods estimate identical results for the pressure distribution with a global error of 0.02% (order 4). Table 1 shows the remaining global errors for orders 0 and 2 and the computation time of the SLFMM. The displacement in Fig. 15 demonstrates, with an error of 0.01%, that the SLFMM determines the result accurately. By changing the order number and/or increasing the convergence criteria, the SLFMM can be accelerated. However, these changes implicate a higher global error for the pressure distribution. The MIM needs 130 s to calculate the result, whereas the new method requires 50 s (order 4). Nevertheless, by increasing the convergence to $10^{-2}$ the SLFMM is faster, with a slight increase of the global error from order 4 to 2. Furthermore, the development of the global error concerning the SLFMM iterations, the deviation between the two methods, and the chosen order number are studied in Fig. 16. The trend of the global error is similar for the single half-sphere, which is why only the development of the 25 half-spheres is shown. During the initial iteration, the global error is infinite, but it decreases rapidly and reaches the defined accuracy after 10 iterations. The global error between the MIM and the SLFMM approach decreases exponentially and changes insignificantly from iterations 7 to 10. As seen before, the result differences between the order numbers are negligible, merely order 0 shows a higher global error than the other order numbers.

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Result for 25 half-spheres: global error $\epsilon_g$ between SLFMM iterations, MIM vs. SLFMM and order numbers 0 to 8.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 17. Results for 1 hole: contact pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis y = 0 and δ = -1.0 mm.</figcaption></figure>

Besides the surface with asperities, two plane surfaces with one (Fig. 8c) and 25 hole(s) (Fig. 8d) are calculated to prove that the SLFMM can determine different types of contact problems. Both surfaces consist of 10 000 grid points with a resolution of $h_x = h_y = 0.1$ mm. The pressure distribution for order 0, 2 and 4 is shown in Fig. 17 and demonstrates how the SLFMM approximates the result of the reference method. For order 0, the global error is 28.76% and for orders 2 and 4, it drops to 1.17% and 0.05%, respectively. Fig. 18 illustrates the displacement for both methods, with a global error of 0.02% (order 4).

<figure><img src="image_21.png" /><figcaption>Fig. 18. Results for 1 hole: displacement for MIM and SLFMM on axis y = 0 and δ = -1.0 mm.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. Result for 25 holes: pressure distribution over the entire surface (SLFMM). The orange-scaled color bar on the right indicates the contact pressure value.</figcaption></figure>

Additionally, the SLFMM was applied to a surface with 25 holes (8), where parameter b is 2. Fig. 19 presents the pressure distribution influenced by the interaction of the contact points for the total surface.

At the outer edges of the surface and the edges of each hole, high pressures are detected. Moreover, the pressure distribution of the exterior holes is asymmetric, which is in particular noticeable at the four corner holes. Again, this difference in the pressure distribution can be seen in Fig. 20.

<figure><img src="image_23.png" /><figcaption>Fig. 20. Results for 25 holes: pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis y = 0 and δ = -1.0 mm.</figcaption></figure>

In total, the SLFMM required 15 iterations to compute the pressure distribution, with an overall error of 31.67% (order 0), 0.80% (order 2) and 0.05% (order 4). For the displacement (Fig. 21) the global error between the two methods is 0.05%. A current limitation for the SLFMM for this case is the computational time to determine the pressure distribution.

<figure><img src="image_24.png" /><figcaption>Fig. 21. Results for 25 holes: displacement for MIM and SLFMM on axis y = 0 and δ = -1.0 mm.</figcaption></figure>

The MIM required 143 s, whereas the SLFMM needed 288 s (order 4) or 248 s (order 2). Even by changing the convergence criteria to $10^{-2}$, the SLFMM still required more computational time than the old method (Table 2). In comparison to the surface with 25 spheres, the SLFMM required more iterations to estimate the result. Furthermore, the elastic half-space is with around 80% of the surface in contact, whereas only 60% of the surface with the 25 half-spheres is in contact. In addition, the fluctuating global error for the iterations and between the two methods (Fig. 22) shows that the new method requires more time to stabilize.

<caption>Table 2</caption>

Comparison of the CPU time and global error $ε_g$ for the SLFMM (order 0, 2 and 4) for 25 holes.

<table><thead><tr><td>Order</td><td>Convergence</td><td>CPU time</td><td>Global error ε<sub>g</sub></td></tr></thead><tbody><tr><td>0</td><td>10<sup>-3</sup></td><td>250 s</td><td>31.67%</td></tr><tr><td rowspan="2">2</td><td>10<sup>-2</sup></td><td>183 s</td><td>0.78%</td></tr><tr><td>10<sup>-3</sup></td><td>248 s</td><td>0.80%</td></tr><tr><td rowspan="2">4</td><td>10<sup>-2</sup></td><td>217 s</td><td>0.78%</td></tr><tr><td>10<sup>-3</sup></td><td>288 s</td><td>0.05%</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>**Fig. 22.** Results for 25 holes: global error $ε_g$ between SLFMM iterations, MIM vs. SLFMM and order numbers 0 to 8.</figcaption></figure>

Next, the surface with randomly distributed half-spheres (Fig. 9) is studied. The radii of the half-spheres vary between 4 and 6 mm, whereby the surface consists of 10 000 grid points with a resolution of $h_x = h_y = 0.5$ mm. To guarantee sufficient contact between the surface and the half-space, the penetration was changed to –4mm for this evaluation only. For a general overview, Fig. 23 illustrates the pressure distribution across the entire surface. Here, the most significant observation is the asymmetric pressure distribution. This asymmetry is especially visible in the transition area between contact and separation, which is caused by the interaction of the asperities. Fig. 24 verifies this asymmetric pressure distribution, where the contact points of each half-sphere are shifted. Apart from this, the SLFMM approximates the result within 10 iterations for order 4 with an error of 0.04%. In Fig. 25 the displacement for both methods is illustrated, with a global error of 0.01% (order 4). To calculate the surface with the randomly distributed half-spheres, the Matrix Inversion Method needed 123 s. Table 3 indicates the computation time of the SLFMM for each order number with its corresponding error.

<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Results for a surface with randomly distributed half-spheres: pressure distribution over the entire surface (SLFMM). The orange-scaled color bar on the right indicates the contact pressure value.</figcaption></figure>

It is observed that the new method is faster than the reference method. For instance, order 4 requires 85 s, which can be reduced to 64 s by increasing the convergence criteria. It can be seen, that the time reduction from order 4 to 2 is often low, with a drastic increase for the global error (cf. Table 3). Consequently, in some cases, it is more beneficial to change the convergence criteria rather than the order number. Nonetheless, the SLFMM offers a good compromise between computational time and accuracy. Finally, the development of

<figure><img src="image_27.png" /><figcaption>**Fig. 24.** Results for a surface with randomly distributed half-spheres: pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis $y = 0$ and $δ = -4.0$ mm.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 25.** Results for a surface with randomly distributed half-spheres: displacement for MIM and SLFMM on axis $y = 0$ and $δ = -4.0$ mm.</figcaption></figure>

<caption>**Table 3**</caption>

Comparison of the CPU time and global error $ε_g$ for the SLFMM (order 2 and 4) for a surface with randomly distributed half-spheres.

<table><thead><tr><td>Order</td><td>Convergence</td><td>CPU time</td><td>Global error ε<sub>g</sub></td></tr></thead><tbody><tr><td>0</td><td>10<sup>-3</sup></td><td>61 s</td><td>16.54%</td></tr><tr><td rowspan="2">2</td><td>10<sup>-2</sup></td><td>47 s</td><td>0.51%</td></tr><tr><td>10<sup>-3</sup></td><td>64 s</td><td>0.46%</td></tr><tr><td rowspan="2">4</td><td>10<sup>-2</sup></td><td>59 s</td><td>0.26%</td></tr><tr><td>10<sup>-3</sup></td><td>85 s</td><td>0.04%</td></tr></tbody></table>

the global error was examined (Fig. 26). The trend of the global error is similar to that of the 25 half-spheres.

<figure><img src="image_29.png" /><figcaption>**Fig. 26.** Results for a surface with randomly distributed half-spheres: global error $ε_g$ between SLFMM iterations, MIM vs. SLFMM and order numbers 0 to 8.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Fig. 27.** Result for a rough surface: pressure distribution over the entire surface (SLFMM). The orange-scaled color bar on the right indicates the contact pressure value.</figcaption></figure>

After demonstrating that the SLFMM can approximate the contact pressure distribution for simple rough surfaces, the method is now applied to a complex rough surface taken from a real road surface texture scan. Afterwards a surface of 50 mm × 50 mm, with a total of 62 500 grid points and a resolution of $h_x = h_y = 0.2$ mm was chosen (Fig. 10). Fig. 27 shows the computed pressure distribution across the rough surface. Since the emerging parts of the rough surface are randomly shaped, the interactions between them are difficult to recognize. To prove that the MIM and the SLFMM are congruent to one another, chosen sections of the contact pressure distribution and of the displacement are illustrated in Figs. 28 and 29. The SLFMM

<figure><img src="image_31.png" /><figcaption>**Fig. 28.** Results for a rough surface: pressure distribution for MIM and SLFMM (order 0, 2 and 4) on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 29.** Results for a rough surface: displacement for MIM and SLFMM on axis $y = 0$ and $\delta = -1.0$ mm.</figcaption></figure>

approximates the contact pressure distribution within 12 iterations, with a global error of 0.03% (order 4), and the displacement, with an error of 0.02%.

It was demonstrated, that the SLFMM meets the requirement of accuracy. However, does it fulfill the condition of accelerated calculation time as well? Fig. 30 shows the CPU time for the MIM and SLFMM for convergences of $10^{-2}$ and $10^{-3}$. The total CPU time for the MIM is 76 175 s (~21:09 h) and 3767 s (~1:22 h) for the Single Level FMM (order 4), which is a reduction of around 95%. Almost 98% of the SLFMM CPU time is related to the direct computation of the near field, the actual SLFMM remaining requires around 2% of the total time. By choosing order 2 instead of 4, the CPU time decreases slightly, with a moderate rise of the global error to 0.24%.

<figure><img src="image_33.png" /><figcaption>**Fig. 30.** Comparison of the CPU time used by the MIM and the SLFMM (order 4) for a convergence $\epsilon$ equal $10^{-2}$ and $10^{-3}$. Subdividing the SLFMM CPU time in the directly computed part $u_{k,near}$ and remaining operations (local matrix and SLFMM calculation).</figcaption></figure>

Table 4 shows the CPU time and the global error for the SLFMM for orders 0, 2, 3, and 4. Since the asperities for the rough surface are not symmetric, order 3 was also evaluated. It can be observed that the global error from order 2 to order 3 decreases significantly, whereas the improvement from order 3 to 4 is comparably low. Moreover, there is no particular acceleration from order 4 to 2, the latter is only 1% faster. However, the computational time decreases by around 22% when the convergence increases, whereby the global error rises slightly.

<table><caption>**Table 4**
Comparison of the CPU time and global error $\epsilon_g$ for the SLFMM (order 2, 3 and 4) for rough surface.</caption><thead><tr><td>Order</td><td>Convergence</td><td>CPU time</td><td>Global error $\epsilon_g$</td></tr></thead><tbody><tr><td>0</td><td>$10^{-3}$</td><td>~ 1:20 h</td><td>13.52%</td></tr><tr><td rowspan="2">2</td><td>$10^{-2}$</td><td>~ 1:02 h</td><td>0.32%</td></tr><tr><td>$10^{-3}$</td><td>~ 1:21 h</td><td>0.24%</td></tr><tr><td rowspan="2">3</td><td>$10^{-2}$</td><td>~ 1:03 h</td><td>0.23%</td></tr><tr><td>$10^{-3}$</td><td>~ 1:21 h</td><td>0.06%</td></tr><tr><td rowspan="2">4</td><td>$10^{-2}$</td><td>~ 1:03 h</td><td>0.22%</td></tr><tr><td>$10^{-3}$</td><td>~ 1:22 h</td><td>0.03%</td></tr></tbody></table>

Furthermore, it is of interest to investigate when the break-even point between the two methods is reached. Fig. 31 shows the time development of the MIM, the SLFMM, and complexity $O(N)$ regarding the surface size. It illustrates, that the SLFMM becomes significantly faster than the Matrix Inversion Method when the domain is greater than 5000 DoFs. In addition, the MIM reaches its limitation at 62 500 grid points, whereas the SLFMM can be used for larger surfaces. Apart from that, the difference between the two convergence criteria tested is comparably low in comparison to the overall improvement between the SLFMM and the MIM.
<figure><img src="image_34.png" /><figcaption>Fig. 31. Comparison of the CPU times used by complexity O(N) the MIM and the SLFMM (order 4) for a convergence ε equal 10⁻² and 10⁻³ per DoF.</figcaption></figure>

Additionally, the development of the total force, global error and ΔCPU were examined regarding the penetration depth for orders 0, 2, 4, and 6. Fig. 32(a) presents the total contact force, which is the sum of

<figure><img src="image_35.png" /><figcaption>Fig. 32. Comparison of the developed SLFMM (order 0, 2, 4 and 6) for (a) the total force, (b) the global error εg and (c) the ΔCPU time for penetrations δ between -3.0 mm to -0.5 mm regarding the MIM.</figcaption></figure>

all contact pressures. For penetrations between -2.0 to -0.5 mm, the SLFMM is equivalent to the reference method with a maximal deviation of 0.20%. When the penetration decreases to -3.0 mm, the variance between the two methods rises to 0.5% (order 0). This deviation is also visible in Fig. 32(b), where the global error is 21.40% (order 0), 0.55% (order 2) and 0.04 (order 4 and 6) for a penetration of -3.0 mm. According to this result, the SLFMM computes accurate results for an order number equal to 2 or bigger. The last Fig. 32(c) illustrates the ΔCPU time between the MIM and SLFMM. When the result is located under zero, the SLFMM is faster than the MIM. Firstly, the new method is faster for all tested penetrations and with each order number. Moreover, there is no significant difference in time between order 2, 4 and 6. To summarize, Fig. 32 shows that the difference between orders 0, 2, 4 and 6 is not substantial for the total force and CPU time, but for the global error. As a consequence, in the case of a road surface, order 4 is a good compromise between the global error and the computational time.

## 4. Conclusion

A Single Level FMM has been developed, to accelerate the calculation between a rough surface and an elastic half-space. Therefore, Boussinesq's contact problem has been approximated by solid harmonics, the so-called multipole expansion. Subsequently, the accuracy of the multipole expansion was validated regarding the analytical displacement. The analytical displacement was computed for a paraboloid of revolution and compared to the multipole expansion. After verifying the accuracy of the multipole expansion, the Single Level FMM was used to compute the pressure distribution for rough surfaces of different complexities.

The numerical validation for the surfaces consisting of half-spheres, holes and randomly distributed half-spheres was performed by the Matrix Inversion Method and the SLFMM. It has been observed that the displacement and pressure distribution were congruent. For all evaluations, the global error was under 0.15% (order 4). Furthermore, the new method estimated the results for the 25 spheres and the surface with randomly distributed half-spheres faster than the reference method. Furthermore, the efficiency was further enhanced by increasing the convergence criteria. It was observed that the CPU reduction from order 4 to 2 is insignificant, comparable to the change of the convergence criteria.

Finally, the new method was used to determine the pressure distribution for a rough surface (road surface sample) and compared to the reference method. It has been demonstrated that both methods are in good agreement, whereas the SLFMM approximated the result within 12 iterations (order 4). In addition, the new method required only 5% of the CPU time of the reference method. In addition, different order numbers for the SLFMM were evaluated to verify the total force, global error, and CPU time. Here, it was proven, that the new approach requires at least order 2 to meet the required accuracy. To summarize, it can be concluded that the Single Level FMM fulfills both the requirement of accuracy and efficiency for solving the contact problem.

Nevertheless, these results can be improved further by a Multi Level FMM approach, which recursively divides the domain further than level 2. With the Multi Level FMM approach, the field of applications could be extended to friction [53,57], fractal surfaces [57–59] and nonlinear contact stiffness [12].

## CRediT authorship contribution statement

**Claudia Stiebritz:** Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis. **Hai-Ping YIN:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Formal analysis. **Julien Cesbron:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Formal analysis.
## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work has been jointly funded by the Région Pays de la Loire and Université Gustave Eiffel.

## Appendix A. Derivation of solid harmonics

Consider two points $M(R, \theta, \phi)$ and $S(\rho, \alpha, \beta)$ that are well-separated by distance $r$ and placed in a reference frame, with its origin at $O(x_0, y_0)$ (Fig. A.33). The point $S$ is within the contact area (orange) and close (distance: $\rho$) to the origin $O(x_0, y_0)$. In contrast, the point $M$ lies outside the contact surface and is separated from the origin ($\rho < R$) by distance $R$. Then, the inverse of the distance $1/r$ can be approximated by the cosine theorem ($r^2 = \rho^2 + R^2 - 2\rho R \cos \gamma$) with $\rho < R < r$.

<figure><img src="image_36.png" /><figcaption>Fig. A.33. Points $S(\xi, \eta)$ and $M(x, y)$ well-separated by distance $r$, placed in a reference frame, with its origin at $O(x_0, y_0)$.</figcaption></figure>

$$ \frac{1}{r} = \frac{1}{\sqrt{\rho^2 + R^2 - 2R\rho \cos \gamma}}, \qquad (A.1) $$

$$ = \frac{1}{R\sqrt{\frac{\rho^2}{R^2} + 1 - 2\frac{\rho}{R}\cos\gamma}}, \qquad (A.2) $$

$$ = \frac{1}{R\sqrt{t^2 + 1 - 2zt}} \quad (t = \frac{\rho}{R}, z = \cos \gamma). \qquad (A.3) $$

The term $\sqrt{t^2 + 1 - 2zt}$ can be rewritten as a generating function of the Legendre polynomials $P_n$.

$$ \frac{1}{\sqrt{t^2 + 1 - 2zt}} = \sum_{n \ge 0} P_n(z) t^n \quad (t < 1), \qquad (A.4) $$

$$ \frac{1}{r} = \frac{1}{\rho} \sum_{n=0}^{\infty} P_n(z) t^n, \qquad (A.5) $$

$$ \frac{1}{r} = \sum_{n=0}^{\infty} P_n(\cos \gamma) \frac{\rho^n}{R^{n+1}} \quad (\rho < R). \qquad (A.6) $$

After combining Eq. (A.6) with the addition theorem for the Legendre polynomials, the regular $R_{n,m}$ and irregular $S_{n,m}$ solid harmonics can be derived [56,60]. The solid harmonics are given by:

$$ R_{n,m}(\text{OS}) = \frac{1}{(n+m)!} P_n^m(\cos \alpha) e^{im\beta} \rho^n, \qquad (A.7) $$

$$ S_{n,m}(\text{OM}) = (n-m)! P_n^m(\cos \theta) e^{im\phi} \frac{1}{R^{n+1}}, \qquad (A.8) $$

where $R_{n,m}$ defines the regular and $S_{n,m}$ the irregular solid harmonics.

## Appendix B. Recurrence relation for solid harmonics

When using solid harmonics, it is important to study their recurrence relations. The first recurrence relations, which are fundamental to calculate higher-order solid harmonics, are written as following:

$$ R_{0,0} = 1, \qquad (B.1) $$

$$ R_{1,0} = x_3, \qquad (B.2) $$

$$ S_{0,0} = \frac{1}{R}, \qquad (B.3) $$

$$ S_{1,0} = \frac{y_3}{R^3}. \qquad (B.4) $$

Afterwards, the higher-order regular and irregular solid harmonics can be determined as following:

$$ R_{n,n} = \frac{x_1 + ix_2}{2n} R_{n-1,n-1}, \qquad (B.5) $$

$$ R_{n,m} = \frac{(2n-1)x_3 R_{n-1,m} - \rho^2 R_{n-2,m}}{n^2 - m^2}, \qquad (B.6) $$

$$ S_{n,n} = \frac{(2n-1)(y_1 + iy_2)}{R^2} S_{n-1,n-1}, \qquad (B.7) $$

$$ S_{n,m} = \frac{(2n-1)y_3 S_{n-1,m} - ((n-1)^2 - m^2)S_{n-2,m}}{R^2}. \qquad (B.8) $$

The vector $x$ can be described by the spherical coordinates $x_1 = (\rho \sin \alpha \cos \beta)$, $x_2 = (\rho \sin \alpha \sin \beta)$, $x_3 = (\rho \cos \alpha)$ and the vector $y$ by $y_1 = (R \sin \theta \cos \phi)$, $y_2 = (R \sin \theta \sin \phi)$, $y_3 = (R \cos \theta)$, where $\rho$ and $R$ describe the distance to the origin.

## Data availability

Data will be made available on request.

## References

[1] Persson B, Albohr O, Tartaglino U, Volokitin A, E. T. On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion. J Phys Condensed Matter 17 2005;1-63.

[2] Persson B. Contact mechanics for randomly rough surfaces. Surf Sci Rep 2006;61:201-27.

[3] Allwood J. Survey and performance assessment of solution methods for elastic rough contact problems. J Tribol 2005;10-2.

[4] Johnson KL. Contact mechanics. Cambridge University Press; 1985.

[5] Zhou W, Zhao Y, Yuan H, Liu J, Wang X. Research on the contact pressure calculation method for the misaligned elastomeric journal bearing. J Mar Sci Eng 2022;10(2):1-19.

[6] Biao W, Yanyan L, Yuanlei Z, Xiaojun Z. Study on detection technology of the contact pressure on the electrical contacts of relays. 26th Int Conf Electr Contacts (ICEC 2012) 2012;161-4.

[7] Zhou P, Lin P, Wu C, Li Z. Effect of nonuniformity of the contact pressure distribution on the electrical contact resistance in proton exchange membrane fuel cells. Int J Hydrog Energy 2011;(36):6039-44.

[8] Askari E, Flores P, Dabirrahmani D, Appleyard R. Study of the friction-induced vibration and contact mechanics of artificial hip joints. Tribol Int 2014;(70):1-10.

[9] Hidayat T, Ammarullah M, Saputra E, Lamura M, Chethan K, Ismail R, et al. A method for estimating the contact area of a dual-mobility total hip prosthesis. AIP Adv 2024;(14):1-9.

[10] Vo K, Tieu A, Zhu H, Kosasih P. A 3D dynamic model to investigate wheel-rail contact under high and low adhesion. Int J Mech Sci 2014;(85):63-75.

[11] Sun Y, Zhai W, Ye Y, Zhu L, Guo Y. A simplified model for solving wheel-rail non-hertzian normal contact problem under the influence of yaw angle. Int J Mech Sci 2020;(174):1-13.

[12] Andersson P, Kropp W. Time domain contact model for tyre/road interaction including nonlinear contact stiffness due to small-scale roughness. J Sound Vib 2008;318:296-312.

[13] Cesbron J, Yin H-P, Anfosso-Lédée F, Duhamel D, Le Houédec D, Feng Z. Numerical and experimental study of multi-contact on an elastic half-space. Int J Mech Sci 2009;51(1):33-40.

[14] Dubois G, Cesbron J, Yin H, Anfosso-Lédée F. Numerical evaluation of tyre/road contact pressures using a multi-asperity approach. Int J Mech Sci 2012;54(1):84-94.

[15] Cesbron J, Anfosso-Lédée F, Duhamel D, Yin H, Le Houedec D. Prediction of contact stresses for tyre-road noise modelling. Euronoise 2006 2018;1-6.
[16] Hertz H. Über die berührung fester elastischer körper. Journal für die reine und angewandte Mathematik 1882;92:156–71.

[17] Boussinesq J. Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques. Gauthier-Villars; 1885.

[18] Love A. The Stress Produced in a Semi-Infinite Solid by Pressure on Part of the Boundary. Trans. R. Soc. Lond. Ser. A, Containing Papers of a Mathematical or Physical Character 1929;228:377–420.

[19] Sneddon IN. The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile. Internat J Engrg Sci 1965;3:47–57.

[20] Ciavarella M. The generalized cattaneo partial slip plane contact problem. i. theory. Int J Solids Struct 1998;35(18):2349–62.

[21] Ciavarella M. The generalized cattaneo partial slip plane contact problem. II. examples. Int J Solids Struct 1998;35(18):2363–78.

[22] Truman C, Sackfield A, Hills D. Contact mechanics of wedge and cone indenters. Int J Mech Sci 1995;37(3):261–75.

[23] Kalker J. Three-dimensional Elastic Bodies in Rolling Contact. Kluwer Academic Publishers, Dordrecht; 1990.

[24] Chen W, Wang Q. A numerical model for the point contact of dissimilar materials considering tangential tractions. Mech. Mater 2008;40(11):936–48.

[25] Cesbron J, Anfosso-Lédée F, Yin H-P, Duhamel D, Houédec DL. Influence of road texture on tyre/road contact in static conditions. Road Mater Pavement Des 2008;9(4):689–710.

[26] Webster MN, Sayles RS. A numerical model for the elastic frictionless contact of real rough surfaces. J Tribol 1986;314–20.

[27] Brandt A, Lubrecht AA. Multilevel matrix multiplication and fast solution of integral equations. J Comput Phys 1990;210–348.

[28] Polonsky IA, Keer LM. Fast methods for solving rough contact problems: A comparative study. J Tribol 2000;36–41.

[29] Allwood J, Ciftci H. An incremental solution method for rough contact problems. Wear 2004;258:1601–15.

[30] Frérot L, Bonnet M, Molinari J-F, Anciaux G. A Fourier-accelerated volume integral method for elastoplastic contact. Comput Methods Appl Mech Engrg 2019;351:951–76.

[31] Röttger M, Sanner A, Thimons LA, Junge T, Gujrati JM, Nöhring WG, et al. Contact engineering—Create, analyze and publish digital surface twins from topography measurements across many scales. Surf Topogr Metrol Prop 2022;10:1–20.

[32] Borodich F, Pepelyshev X. A multiscale statistical analysis of rough surfaces and applications to tribology. Mathematics 2024;1–19.

[33] Rohklín V. Rapid solution of integral equations of classical potential theory. J Comput Phys 1985;60:187–207.

[34] Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comput Phys 1987;315–48.

[35] Greengard LF. The Rapid Evaluation of Potential Fields in Particle Systems. MIT Press, Cambridge; 1988, p. 14–86.

[36] Gomez JE, Power H. Multipole direct and indirect BEM for 2D cavity flow at low Reynolds number. Eng Anal Bound Elem 1997;19:17–31.

[37] Fujiwara H. The fast multipole method for solving integral equations of three-dimensional topography and basin problems. Geophysical Journal International 2000;140:198–210.

[38] Peirce AP, Napier JAL. A spectral multipole method for efficient solution of large-scale boundary element models in elastostatics. Int J Numer Methods Eng 1995;38:4009–34.

[39] Liu Y. A new fast multipole boundary element method for solving large-scale two-dimensional elastostatic problems. Internat J Numer Methods Engrg 2005;65:863–80.

[40] Chen YH, Chew WC, Zeroug S. Fast multipole method as an efficient solver for 2D elastic wave surface integral equations. Springer-Verlag; 1997, p. 495–505.

[41] Fu Y, Klimkowski KJ, Rodin GJ, Berger E, Browne JC, Singer JK, et al. A fast solution method for three-dimensional many-particle problems of linear elasticity. Int J Numer Methods Eng 1998;42:1215–29.

[42] Duc Pham A, Mouhoubi S, Bonnet M, Chazallon C. Fast multipole method applied to symmetric Galerkin boundary element method for 3D elasticity and fracture problems. Eng Bound Elem 2012;36:1838–47.

[43] Takahashi T, Nishimura N, Kobayashi S. A fast boundary integral equation method for three-dimensional elastodynamics in time domain. Eng Anal Bound Elem 2002;27:165–80.

[44] Chaillat S, Bonnet M, Semblat JF. A fast multipole method formulation for 3D elastodynamics in the frequency domain. Comptes Rendus Mécanique 2007;335:1–7.

[45] Antonini G. Fast multipole method for time domain PEEC analysis. IEEE Trans. Mob. Comput. 2003;2:275–87.

[46] Shen L, Liu YJ. An adaptive fast multipole boundary element method for three-dimensional potential problems. 2006, p. 681–91.

[47] Hu B, Li C, Niu Z. A new multi-level strategy of numerical integration in the fast multipole BEM for analyzing 3D potential problems. Comput Math Appl 2024;161:174–89.

[48] Bapat M, Shen L, Liu Y. Adaptive fast multipole boundary element method for three-dimensional half-space acoustic wave problems. Eng Anal Bound Elem 2009;33:1113–23.

[49] Schanz M. Fast multipole method for poroelastodynamics. Eng Anal Bound Elem 2017;89:50–9.

[50] He H-J, Zhang Y, Li K-Y. An efficient method for analysis of electromagnetic scattering from a 3-D coated object half-buried in PEC rough surface. Electromagnetics 2023;43(1):14–23.

[51] Takahashi T. A fast time-domain boundary element method for three-dimensional electromagnetic scattering problems. J Comput Phys 2023;482:1–26.

[52] Mateus D, Dias Jr A, Campos L, dos Santos J, Albuquerque E. The fast multipole boundary element method for anisotropic material problems under centrifugal loads. Eng Anal Bound Elem 2024;(162):75–86.

[53] Cesbron J, Yin H-P. Contact analysis of road aggregate with friction using a direct numerical method. Wear 2010;268(5):686–92.

[54] Liu Y. Fast Multipole Boundary Element Method, Theory and Application in Engineering. Cambridge University Press; 2009, p. 47–83.

[55] Hou J, Chen Y. A new approach to multi-domain fast multipole boundary element method. Eng Anal Bound Elem 2024;166:1–17.

[56] Liu YJ, Nishimura N, Otani Y, Takahashi T, Chen XL, Munakata H. A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model. J Appl Mech 2005;72(1):115–28.

[57] Monti J, Pastewka L, Robbins M. Fractal geometry of contacting patches in rough elastic contacts. J Mech Phys Solids 2022;(160):1–18.

[58] Wei D, Zhai C, Hanaor D, Gan Y. Contact behaviour of simulated rough spheres generated with spherical harmonics. Int J Solids Struct 2020;(193-194):54–68.

[59] Zhai C, Hanaor D, Gan Y. Contact stiffness of multiscale surfaces by truncation analysis. Int J Mech Sci 2017;(131-132):305–16.

[60] Bonnet M. Méthode d'éléments de frontière. ENSTA Paris; 2011, p. 54–91.