Investigation of Contact Clusters Between Rough Surfaces

Abstract

We investigate the persistence of micro-contacts between two elastic random rough surfaces by means of a simple model for quasi-static sliding. Contact clusters are calculated with the Boundary Element Method, then surfaces are repeatedly displaced to study the evolution of the original contact area. While the real contact area remains constant due to the rejuvenation of micro-contacts, the original contact clusters are progressively erased and replaced by new ones. We find an approximate exponential decrease of the original real contact area with a characteristic length that is influenced both by statistics of the contact cluster distribution and physical parameters. This study aims to shine light on the microscopic origins of phenomenological rate-and-state friction laws and the memory effects observed in frictional sliding.

1 Introduction

The macroscopic friction force between sliding surfaces originates from the interaction of micro-contacts at the scale of surface asperities [1]. This evidence has been encoded in the Bowden-Tabor law, stating that the friction force is proportional to the real contact area, which is only a small fraction of the apparent one and is, as a first approximation, proportional to the applied load [2]. From this, two strategies can be developed to calculate the friction force between contact surfaces: either by calculating the real contact area as a function of physical and surface topological parameters, or by finding an effective coarse-grained friction law reproducing the emergent macroscopic behavior.

The former approach has been developed since the first work of Greenwood and Williamson [3] and several techniques have been used to calculate the real contact area [4], ranging from analytical calculations [5,6,7,8,9] to finite element [10,11,12] and boundary element methods [13,14,15,16], also extended to plasticity [17]. These studies have allowed to conclude that, for elastic surfaces with a given topography h(x, y), the real contact area is dominated by the root mean square of the slopes $R_{\text{dq}} \equiv \sqrt{\langle |\nabla h|^2\rangle }$, so that $A_{\text{real}}= \kappa N/(E^{\prime } R_{\text{dq}})$, where N is the load, $E^{\prime }$ the effective Young’s modulus and $\kappa$ a proportionality constant.

Although we neglect the dependence on other parameters included in $\kappa$ and the practical difficulty to estimate $R_{\text{dq}}$ [18], the story is further complicated for dynamic frictional sliding, since experimental observations have shown that the real contact area is not fixed: for example, it is decreasing during the transition from static to dynamic friction [19], it grows with time due to aging effects [20], and during rough on rough sliding micro-contacts continuously change. Moreover, given the analogies between the transition to sliding with a fracture mechanics process [21], memory effects due to the history-dependent interface stress concentrations [22, 23] play a crucial role.

In order to take into account these phenomena, effective friction laws can be a viable alternative approach. In particular, since the early works of Dieterich [24] and Ruina [25], rate-and-state friction laws have been successfully used to extend the classical Amontons–Coulomb friction force by introducing, as suggested by the name, a dependence on the velocity and a state variable, which takes into account the memory effects at the interface. These models have found particular relevance in geophysics for the modeling of friction along faults, and the state variable has been recognized as an essential factor to model the macroscopic sliding [26,27,28,29,30,31]. These laws have been used to address problems not only about the frictional sliding between elastic solids [32,33,34,35], but also to glaciers [36], gouges with granular material [37, 38] and boundary lubrication problems [39]. In other studies, the role of the viscoelasticity on the collective behavior of multiple contacts has been investigated during the frictional sliding, establishing a connection between the evolution of the state variable and microscopic features [40, 41].

One of the most common and simplest form of rate-and-state law was proposed by Ruina and relates the instantaneous friction coefficient f(x, t), as a function of the time t and surface point x, to the velocity field v(x, t) and the state variable field $\phi (x,t)$ in a set of two coupled equations [25]:

$$\begin{aligned} f(x,t)&= f_0 + a \ln {\frac{v(x,t)}{v^{*}}} + b\ln {\frac{\phi (x,t)}{\phi ^{*}}} \end{aligned}$$

(1)

$$\begin{aligned} {\dot{\phi }}(x,t)&= 1 - \frac{v(x,t)\phi (x,t)}{D_{\text{c}}} \end{aligned},$$

(2)

where $f_0$, a, b are constants that must be fixed for each system and $v^{*}$ and $\phi ^{*}$ are normalizing factors. The second equation takes into account the memory part and, for constant velocity, it implies an exponential decrease of the state variable over time ruled by a characteristic slip distance $D_{\text{c}}$. This is connected to the sliding distance required for the interface to erase the information of previous slips and nucleation events [24]. In other words, the original real contact area is continuously replaced by new contacts, and the amount of original contacts decreases exponentially with the sliding distance, characterized by a decay length $D_{\text{c}}$ roughly corresponding to the average length of the micro-contacts [28]. This observation can provide the connection between the macroscopic coarse-grained framework of the rate-and-state friction laws with the microscopic approach involving the calculation of the real contact area. For this reason, it is important to investigate the relation between parameters characterizing the surface roughness and those used in the effective laws.

In this paper, we investigate the persistence of the original contact clusters between two elastic random rough surfaces. The contact problem is solved with Tamaas, a boundary element method solver [42, 43] to calculate the real contact area, then surfaces are gradually displaced one grid point at a time, to assess the distance required to erase the original clusters as a function of the model parameters. In particular, we investigate the influence of the Hurst exponent, which characterizes surface roughness, and elastic properties, on the rejuvenation of micro-contacts. We also discuss how the results differ between the fully mechanistic treatment of contact offered by Tamaas, and a purely geometrical approach known as the overlap model.

2 Methods

2.1 Surface Generation

Many natural surfaces can be described as self-affine fractal surfaces [44, 45], that is the same surface features can be observed on several length scales. In this case, it is possible to encode all the topographic information into the language of power spectrum density (PSD) [18]. Thus, the spectral analysis provides a method to generate a self-affine topography starting from the PSD, performing the inverse Fourier transformation. These methods are also implemented in Tamaas [43], so that in our simulations surface generation and BEM solutions are both provided by this software. The idealized PSD of a self-affine randomly rough surface follows the distribution:

$$\begin{aligned} C^{iso}(q) = C_{0}\left\{ \begin{array}{rcl} 0 &{}&{}\text { if } q<q_{\text{r}}\\ q_{\text{L}}^{-2-2H} &{} &{} \text { if } q_{\text{r}}\le q<q_{\text{L}}\\ q^{-2-2H} &{} &{} \text { if } q_{\text{L}}\le q<q_{\text{s}}\\ 0 &{} &{} \text { if } q\ge q_{\text{s}} \end{array} \right. \end{aligned},$$

(3)

where $q_{\text{r}}<q_{\text{L}}<q_{\text{s}}$ are roll-off, large and small wavenumber cutoffs, respectively. H is the Hurst exponent, which is related to the fractal dimension $H = 3 - D$. An example of PSD is shown in Fig. 1a. In the following, we will adopt square surfaces with $N=1024$ points for each side, with periodic boundary conditions and an elementary discretization length l conventionally fixed to $l=1$ $\mu$m. We use the conventions introduced in [46] to define the spectrum parameters, i.e., the adimensional wavenumber $q \equiv k N\cdot l/2\pi$, being k is the wavenumber, so that q is the number of waves per length. By default, in the frequency domain we will use $q_{\text{r}}=8$, $q_{\text{L}}=8$, $q_{\text{s}}=128$, corresponding to wavelengths ($\lambda = Nl/q$) spanning from 8l to 128l, and Hurst exponent $H=0.7$, which have been chosen to have a statistically significant number of contact clusters in the domain and are within the value ranges used in [46]. In [43] we do not explicitly set $C_{0}$, but with these settings the theoretical values of the root mean square of heights and slopes are, respectively, $R_q = 16.8 \; l$ and $R_{\text{dq}}=2.7 \cdot 10^3$.

2.2 BEM Solver

After the generation of the random rough surfaces, we solve the elastic normal contact problem with Tamaas, implementing the Polonsky–Keer algorithm [47, 48] accelerated with the Fast Fourier Transform [49]. An example of the solution for the contact clusters is shown in Fig. 1b with the default settings and a ratio between applied normal pressure and effective Young’s modulus corresponding to $P/E^{\prime } = 0.01$. This ratio is the relevant parameter of the model that will be used in the following sections.

Fig. 1

An example of idealized PSD (left). Example of contact clusters between two self-affine randomly generated surfaces, whose axis length are expressed in units of the elementary length l (right)

2.3 Sliding Algorithm and Contact Clusters

In order to simulate a quasi-static sliding of the model, we shift one of the two surfaces along the horizontal direction by l, wrapping-around the foremost edge due to the periodic boundary conditions. The sliding direction is fixed along the x-axis. Then, the calculation of the contact solution is repeated. After a fixed number n of these sliding steps, whilst the real contact area is constant on average, some of the initial contact clusters have disappeared, others have been modified and new contact clusters have been created. The simulation of the evolution of the real contact area has been provided in the videos attached as supplemental material.

In order to study the persistence of the initial contact points, we must identify and classify the clusters evolving during the quasi-static sliding. The main difficulty is that the contact clusters may merge or split during the evolution. To address this problem, we insert all the n contact maps of size $N\times N$ into a 3-D matrix $n\times N \times N$ and label all its pixels according to a 26-connected pixel connectivity. Then, we use the Hoshen–Kopelman algorithm to solve this connected-component labeling problem [50].

The Hoshen-Kopelman algorithm consists of two rounds of scans over all pixels which needs to be labeled. In the first round, we scan each objective pixel and its surrounding pixels (only half of the neighbors, i.e., 13 pixels, to avoid scanning two times the same pairs). If there is any surrounding pixel that has already a label, then the chosen pixel will be labeled by the smallest index we can have. Also, if there are more than one label that is surrounding the chosen pixel, all these labels will be marked as ”equivalent”. If no label can be found from the surrounding pixels, a new label will be assigned to the chosen pixel. For the second round, all the pixels that have the ”equivalent” labels will be unified under the same label which has the smallest index.

In this way, we can identify all the contact clusters evolving with the sliding and consider only those already present at the beginning as original contact clusters. An example of the evolution of a single contact cluster identified by means of this algorithm is shown in Fig. 2. This allows to calculate the persistence of the initial contact area, as highlighted in Fig. 3. By repeating this procedure for different surfaces obtained with new seeds of the random generator, it is possible to increase the statistics on the clusters and calculate the averages over all results. Movies with an example of the evolution of the real contact area and the persistence of the original clusters for the default spectrum parameter and $H=0.9$ are provided as supplemental material. This formulation does not consider wear or aging effects of the junctions, so that statistically its behavior is invariant under time and space translation, i.e., initial conditions are not influential.

In this model, what matters from the physical point of view is the junction between asperities, which are represented by the contact clusters. We do not interpret results in term of the single pixels, because they can be affected by discretization issues in simulations but also in experiments due to scan resolution. But independently of the resolution, when two rough surfaces are in contact, junctions are formed and evolve during sliding. For example, in [19], the authors study the evolution and deformation of the original contact area as a collection of the contact junctions, both globally and at single cluster level. A reduction of the real contact area, calculated as a collection of clusters, is observed in the static phase before sliding. Our aim is to provide an approximate model of similar processes involving contact clusters and establish connections with the macroscopic effective laws.

Macroscopic rate-and-state friction laws assume an exponential decay of the original contact area with sliding distance. The state variable in rate-and-state friction is interpreted as the age of micro-contacts. Our model explores the decay of micro-contacts through the combination of geometry, described by rough surfaces in relative sliding, and mechanics, since each configuration is obtained with a BEM solver. Although it does not account for aging in terms of creep, it is nonetheless relevant to investigate the decay of the original contact area with respect to the assumption of the rate-and-state friction.

Fig. 2

Evolution of a single contact cluster for the sliding steps $n=0$, $n=4$, $n=8$, respectively, with the default spectrum parameters

Fig. 3

Evolution of another contact cluster, obtained with the default spectrum parameters. In the snapshots following the first one, we have maintained the shadow of the original cluster, highlighting in red the overlap representing the remaining contact area, i.e., the fraction of area persisting after the translation steps

3 Results

3.1 Remaining Contact Area

We define the remaining contact area (RCA) as the fraction $A_{\text{rc}}$ of the initial contact clusters that persists after n translation steps. Although the area of the contact clusters can also increase during the evolution, the remaining contact area under consideration includes only the fraction of the initial one, so that it is decreasing function of the sliding steps, as shown in Fig. 3. It can also be described as the overlapping part between the original contact area and that one after n steps. As observed in previous studies [24], for rough surfaces the RCA decreases approximately as an exponential function of the sliding distance. Our results confirm this behavior, as shown in Fig. 4, where the RCA value on y-axis is normalized to be between 0 and 1. The interpolating function follows $A_{\text{rc}}= \exp {(-x/d)}$, being $x\equiv l\cdot n$ the sliding distance and d the characteristic decay length. This length corresponds to the parameter $D_{\text{c}}$ of the rate-and-state friction law 1.

It is interesting to observe that also single clusters follow this exponential law with good approximation, i.e., the exponential decrease is found both for the total RCA and the RCA calculated by taking into account only clusters whose initial area is included in a small limited range. It is possible to scale the sliding distance with the root square of the mean value of the area of the clusters, namely $L=\sqrt{\langle A \rangle }$, being A the area of single contact clusters. This rescaling aims to remove the influence of the absolute value of the real contact area using only an a-dimensional quantity characterizing the RCA. Thus, all RCA curves collapse approximately on the same exponential function $A_{\text{rc}}= \exp {(-x/(d L))}$, where d is the a-dimensional decay length expressed in units of the average cluster length L. We have estimated d from the RCA curves with a numerical spline cubic interpolation, finding the distance such that $A_{\text{rc}}=e^{-1}$. For the case of Fig. 4, we have found $d=1.08\pm 0.01$ in units of L, obtained with the default spectrum parameters and by averaging over 1300 repetitions.

Fig. 4

Remaining contact area versus the sliding ratio x/L for the total area and various ranges of cluster area, indicated in the legend with the limits expressed in terms of the average cluster area. For each color, the thick solid lines represent the best fitting exponential function

In Fig. 4, the curves of different cluster bins converge approximately to the decay curve of the total average RCA. This illustrates that the behavior and the decay length d is largely independent of size, shape and distribution of contact clusters.

The collapse of the curves is not perfect due to the statistical dispersion of the probability density function (PDF) of the area of contact clusters. This can be understood with the following consideration: supposing that all surface asperities have the same diameter, i.e., their PDF is a delta function, the sliding length required to separate them would be on average their diameter, with a small dispersion due to their random initial relative position. However, if there are asperities of various diameters, according to a PDF with large variance, there could be clusters generated by a small asperity indenting over a large one, so that the sliding length to erase their contact would correspond to the length of the larger one. For this reason, the normalized RCA obtained by averaging over all the PDF is larger than the RCA obtained over a narrow band of the PDF, independently of the average value of the cluster area. This is confirmed by Fig. 4, where the curve corresponding to the total RCA, having the largest variance, is the top limit for the other ones.

3.2 Geometrical Overlap

Fig. 5

Schematic of the geometrical overlap model

From the results of Sect. 3.1, the behavior of RCA is a combined effect of both contact mechanics and statistics of contact clusters. In order to decouple the problems, in this section we calculate the RCA obtained with a simple model for the geometrical overlap of rough surfaces, without any indentation. Similar models have been used to estimate the real contact area [10, 51, 52]. Supposing that there are two random rough surfaces not in contact, as shown in Fig. 5, we can consider their gap g(x, y) as a function of the coordinates. Given the minimum and the maximum gap, $g_{\text{min}}$ and $g_{\text{max}}$, respectively, we define the contact area as the points (x, y) such that:

$$\begin{aligned} \frac{g(x,y)-g_{\text{min}}}{g_{\text{max}}-g_{\text{min}}} \le f \end{aligned},$$

(4)

where f is an arbitrary number between 0 and 1. In this way, for any pair of rough surfaces, we can obtain a contact map similar to that of Fig. 2 but without any indentation. The number f can be fixed such that, on average, the total real contact area matches that obtained with the BEM solver for the same surfaces. In Fig. 6, we compare the contact clusters obtained for the same surfaces with the contact solution and geometrical overlap model. The two cases share many large contact clusters in the same regions, but it is evident that, in the presence of elastic contact, they are smaller, since the elastic resistance prevents a further indentation. Moreover, there are a larger number of small contact spots, which are points with a gap larger than the threshold for the geometrical overlap model, but coming in contact due the surface deformation in the BEM solution.

This difference is summarized by the PDF of the area of contact clusters for both cases, as shown in Fig. 7. Although the average total real contact area is the same, in the geometrical overlap model there is a larger number of big clusters, so that the average cluster length $L=\sqrt{\langle A \rangle }$ is also larger (in the example in Fig. 7, $\langle A \rangle$ for the geometrical overlap model is twice larger). According to the literature [53], the PDF of the cluster area A of random self-affine surfaces is constant for small clusters, whereas for larger clusters it is the power-law $p(A)\sim A^{-2+H/2}$ ruled by the Hurst exponent, with an exponential cut-off for the largest clusters. The area ranges defining these regions are determined by the spectrum parameters. Our results recover this behavior for small clusters, although the range for the power-law is barely visible due to our parameter choice.

Although the PDF for the two cases are different in both the power-law and the cut-off, the RCA curves follow a similar exponential decrease in both cases, as shown in the inset of Fig. 7. For the geometrical overlap model we observe again a non-perfect collapse of the curves by considering smaller ranges, confirming that this effect is due statistics and not a specific feature of the contact solution or the cluster PDF. The estimated value of d, i.e., the length in units of L required to erase the original contacts, is $1.24\pm 0.01$ and $1.08\pm 0.01$ for the geometrical overlap model and the elastic solution, respectively, both expressed in units of their average cluster length L. Thus, for both models it is slightly larger than L, which is a consequence of the statistical dispersion of the PDF as explained in Sect. 3.1.

From these results, we conclude that the exponential decrease is a feature of the shape of the contact clusters obtained with a random rough self-affine surface whose spectrum is characterized by Eq. 3. The persistence of the original feature is completely erased at most in two average cluster length L, whereas the value of the decay length d corresponds approximately to this value. The non-perfect collapse for limited ranges of cluster area is a statistical effect independent of the cluster distribution and the physics of the contact.

Fig. 6

Comparison of the contact maps obtained with the same surfaces for the geometrical overlap model and the BEM elastic solution

Fig. 7

Comparison of the PDF of the area of the contact clusters between the geometrical overlap model and the BEM elastic solution. Surfaces have been generated with $H=0.7$, $q_{\text{L}}=q_{\text{r}}=8$, $q_{\text{s}}=128$ and results have been averaged over 2000 repetitions. In the inset plot, the comparison of the total RCA curve is reported

3.3 Influence of the ratio $E^{\prime }/P$ and Hurst exponent

As shown in Sect. 2.3, The total remaining contact area follows an exponential decay $A_{\text{r}}= e^{-x/d L}$, where $L=\sqrt{\langle A \rangle }$ is the average cluster length, e.g., square root of the average area. In general, the parameter d, i.e., the distance required to erase the original contacts in units of L, is a function of the system parameters, in particular the ratio between Young’s modulus and applied pressure $E^{\prime }/P$ and the Hurst exponent H. By repeating the simulations illustrated in Sect. 2 for various parameter values, we have found the behavior reported in Fig. 8.

From these results, modifications span in a range of $20 \%$ with respect to the average cluster length L, however some trends are observed. In particular, d is an increasing function of up to $E^{\prime }/P=100$. This seems a counter-intuitive result: if we consider a junction between two asperities, we could expect that in an elastic material prone to deformation, it would be more difficult to displace them and erase their contact, i.e., d should be larger. This result can be understood by considering that it is obtained by averaging all the contact clusters, spanning from small and point-like clusters to large and fractal-like shapes. Thus, similarly to the mechanism explained in Sect. 3.1, in the presence of a larger variance in the cluster area distribution, we expect an increase of the decay length.

This is confirmed by considering the distribution of the contact cluster as a function of their area and shape, reported in Fig. 9 for two representative cases. The area is expressed as a fraction of the total real contact area $A/A_{\text{real}}$, so that, for example, a distribution limited to a value $A/A_{\text{real}}<0.1$ implies that the contact map is only comprised clusters large at most $10\%$ of the real contact area. Note that the root mean square of the slopes is not constant for different H. The shape is classified with a compactness parameter c, that is the ratio between the area of the contact cluster and the area of the square encircling it, so that $c=1$ implies a perfect spot-like cluster and c towards zero a fractal-like shape. Although more rigorous mathematical definitions of compactness exist, this definition is simple for calculations and is sufficient for our aims. The cluster shown in the first snapshot of Fig. 2 is an example of smaller c, whereas the cluster in the first snapshot of Fig. 3 has c close to 1.

For $E^{\prime }/P=10$, the PDF is limited to clusters smaller than $5\%$ of the real contact area, but it has a significant tail for small c, implying the presence of small fractal-like clusters. For a $E^{\prime }/P=900$, the distribution is broader and single clusters large $40\%$ of the total real contact area may appear, although most of the clusters are small and spot-like. This large variance explains why the curve of the decay length stabilizes for larger $E^{\prime }/P$.

The trend with H is approximately constant, except for the smaller value of $E^{\prime }/P$, where in any case variations are limited within $10\%$. In other words, the length in units of L required to erase the original contact configuration is independent of the power-law exponent of the surface spectrum. This is relevant for realistic surfaces, which are not perfectly self-affine and can display spectra with several slopes or power-law exponents depending on the length scale under consideration. This implies that a similar universal behavior is expected to be valid for many real surfaces, as in the Dieterich experiments [24], with a characteristic decay length corresponding approximately to the average contact cluster.

Fig. 8

Influence of the Hurst exponent (left) and the ratio between Young’s Modules and applied pressure (right) on the exponential decay length. Corresponding values of the total real contact area are reported in Fig. 10 in the Appendix

Fig. 9

Histograms of the PDF of the area of contact clusters as a function of the compactness parameter c and their area as a fraction of the real contact area. Results have been obtained with default spectrum parameters, $E^{\prime }/P=10$ and $H=0.4$ (left) and $E^{\prime }/P=900$ and $H=0.7$ (right)

Fig. 10

Real contact areas of the datasets of Fig. 8. In this plot, we are reporting the absolute value of the real contact area, thus we have calculated the values with the correction introduced in [46]. For the remaining contact area used in this study, such correction does not influence the results due to the adopted normalization. For information about the behavior of the real contact area as a function of the spectrum parameters, we refer the reader to the cited paper

4 Conclusions

In this paper, we have investigated the persistence of the original contact clusters between elastic random self-affine surfaces. After their random generation from the given power-law spectrum, the contact problem has been solved with the boundary element method implemented in the software Tamaas. We have shown that the original contact area decreases approximately as an exponential function of the sliding distance. This behavior is also found for limited ranges of area. The characteristic decay length scale corresponds closely to the average cluster length, so that all curves can be approximately scaled by this factor. This is also observed using a simple geometrical overlap model instead of the contact solver, demonstrating that the general behavior is a feature of the statistical properties of the random self-affine surface rather than the contact mechanics, which in turn determines the contact cluster distribution and the value of the average cluster length. By exploring the parameter space determined by the Hurst exponent and the ratio between Young’s modulus and pressure, variations limited within a $20\%$ of the decay length have been observed. In particular, the decay length increases for larger $E^{\prime }/P$, which can be explained by a larger distribution of the area of single contact clusters. The Hurst exponent does not have a significant impact, implying that a similar behavior is expected for realistic surfaces whose spectra can display several slopes. These results provide a microscopic confirmation of the state equation for the widely used macroscopic rate-and-state friction laws, linking the characteristic distance to erase an initial surface state to the characteristic length of the micro-contacts.

Appendix

Supplemental material: The videos of the evolution of the total real contact area, with the default spectrum parameters are provided, video_E=1_H=4.mp4 ($H=0.4$, $E^{\prime }/P=100$) and video_E=9_H=7.mp4 ($H=0.7$, $E^{\prime }/P=900$), respectively. For the former case, we have also provided the video of the remaining contact area video_E=1_H=4_RemainingContactArea.mp4.