Onset of Sliding of Elastomer Multicontacts: Failure of a Model of Independent Asperities to Match Experiments

Modeling of rough frictional interfaces is often based on asperity models, in which the behavior of individual microjunctions is assumed. In the absence of local measurements at the microjunction scale, quantitative comparison of such models with experiments is usually based only on macroscopic quantities, like the total tangential load resisted by the interface. Recently however, a new experimental dataset was presented on the onset of sliding of rough elastomeric interfaces, which includes local measurements of the contact area of the individual microjunctions. Here, we use this more comprehensive dataset to test the possibility of quantitatively matching the measurements with a model of independent asperities, enriched with experimental information about the area of microjunctions and its evolution under shear. We show that, despite using parameter values and behavior laws constrained and inspired by experiments, our model does not quantitatively match the macroscopic measurements. We discuss the possible origins of this failure.

1. Introduction

The mechanical behavior of contact interfaces between rough solids is crucial to understand their tribological properties. The rough contact mechanics community has been developing models in two main directions (see Vakis et al., 2018 for a recent review). First, asperity models in which the contact interface is divided into well-defined microjunctions actually carrying the normal and tangential loads applied to the contacting solids (Braun and Röder, 2002; Ciavarella et al., 2006; Violano and Afferrante, 2019). Each microjunction is ascribed a set of individual properties (e.g., its height, radius of curvature or friction coefficient) necessary to apply some assumed behavior laws [e.g., any contact (Johnson, 1987) or friction law (Le Bot et al., 2019)] when submitted to an external stimulus. The macroscopic behavior of the interface is then the emerging, collective response of the population of microjunctions (Trømborg et al., 2014; Braun and Peyrard, 2018; Costagliola et al., 2018). Second, continuum models in which the input quantity is the full topography of the rough surfaces, and an exact solution of the unilateral contact and friction problem is seeked (Pastewka and Robbins, 2014; Yastrebov et al., 2017; Ponthus et al., 2019), again under some assumptions on the interfacial behavior, concerning, e.g., elasticity, friction, and adhesion.

Each approach can be used to produce two types of results, either deterministic or statistical. Deterministic results are obtained for a given topography (for continuum models) or a given set of model parameters (for asperity models, including the properties of each microjunction) and are thus specific to those input data. They are relevant for quantitative comparison with a particular experiment. In contrast, statistical results are the expected results of a large number of deterministic calculations performed on statistically similar random surfaces. In asperity models, statistical results are obtained when using probability density functions (pdfs) of the microjunction properties (Greenwood and Williamson, 1966; Braun and Peyrard, 2008; Thøgersen et al., 2014). In continuum models, they are usually obtained using the power spectrum density (psd) of the topographies under study (Persson, 2001). In the following, we aim at finding a quantitative match with a specific set of measurements, so we will consider deterministic results.

Both asperity and continuum models have been widely explored in the context of rough contacts under purely normal load, with a recent study explicitly comparing the relative merits of the two approaches (Müser et al., 2017). Several studies aimed at a quantitative comparison between deterministic model results and local, microjunction level measurements (see e.g., McGhee et al., 2017; Acito et al., 2019). In contrast, to our best knowledge, such comparisons have not been reported in the case of sheared multicontacts. Here we will attempt to build an asperity model able to quantitatively match recent measurements performed on the incipient tangential loading and onset of sliding of a rough elastomer slab in contact with a smooth glass plate (Sahli et al., 2018, 2019) (Figure 1A). Those measurements (see a typical example in Figures 1C,D) are particularly interesting and constraining for models because, in addition to the macroscopic loads on the interface, they include the evolution under shear of the individual contact areas and shapes of the many microjunctions forming the interface (Figure 1B).

FIGURE 1

Figure 1. Experiments that our model attempts to reproduce. (A) Sketch of the experimental setup. (B) Typical segmented image of the interface showing individual microcontacts in white, for P = 6.40 N. (C) Concurrent time evolutions of the tangential force Q (red) and the area of real contact (blue), for P = 3.10 N. (D) Area of real contact as a function of the tangential force, for the same data as in (C).

The philosophy of this work is to start with a model of independent asperities like the earthquake model of Braun and Peyrard (2008), enrich it with the recently identified phenomenology of shear-induced area reduction, and genuinely ask the question whether such a model is sufficient to quantitatively match a particular experimental dataset. In other words, we do not aim at a definitive model of the incipient tangential loading and onset of sliding of rough elastomer contacts. Rather, we make one single step ahead compared to the models in the literature, and try to conclude whether this step (including shear-induced area reduction) is sufficient or not. Such an approach can only be fruitful if the values of the model parameters are sufficiently constrained by the experimental dataset, so that one avoids fortuitous agreement. This can be achieved (i) by limiting to the strict minimum the number of parameters that cannot be directly measured experimentally, and (ii) by performing a thorough exploration of the parameter space for those remaining, unconstrained parameters. In this work, we did our best to apply this strategy, which in our case leads to an unsatisfactory agreement. This result is nevertheless a progress in the sense that it clarifies the range of assumptions that remain to be questioned and improved in future studies.

In section 2.1, we describe the asperity model and provide the experimental constraints on the model parameters in section 2.2. Quantitative comparisons between the model and measurements are given in section 3, while in section 4 we discuss the possible reasons for the absence of good matching between the two.

2. Materials and Methods

2.1. Model Description

We consider the frictional interface between a slider of mass M and a track. The tangential displacement of the slider, X(t), is assumed to obey the following equation of motion:

$\begin{array}{ll}M\ddot{X}+M\eta \dot{X}=k_L(vt-X)-F,& (1)\end{array}$

where M is the slider's mass, k_L is the stiffness of the loading spring through which the slider is pulled at constant velocity v, η is an effective viscous parameter accounting for dissipation, e.g., in the air or in the loading spring, the dot indicates the time derivative, and F is the resistive force due to interfacial/adhesive friction.

We assume that the interface is a multicontact made of N independent individual microjunctions, each resisting a force f_i, so that $F=\sum _{i=1}^N{f}_i$. Each microjunction can be in either of two states. First, a pinned state during which the junction acts like a (time-evolving) elastic spring of stiffness k_i, so that f_i = k_i(X(t) − x_i), with x_i the slip displacement of the junction with respect to the track (e.g., x_i = 0 as long as junction i has never been slipping). When a threshold force f_si is reached, the junction enters a slipping state, during which f_i = ϵf_si. Note that ϵ <1, so that f_si and ϵf_si are the analogs, at the junction level, of a static and a dynamic friction force, respectively.

The mechanical behavior of individual junctions is inspired by experimental observations made on the same setup and materials in contact as in Figure 1A, but when the rough slab is replaced by a single smooth sphere (Sahli et al., 2018, 2019; Mergel et al., 2019). The resulting sphere/plane contact is assumed to be representative of an individual microjunction within a multicontact like that of Figure 1B. Those experiments, carried out both for large normal loads (Sahli et al., 2019) and for small (even negative) normal loads (Mergel et al., 2019), have shown that, under increasing shear, the initially circular contact shrinks anisotropically and becomes increasingly ellipse-like. As shown in Sahli et al. (2019) the shrinking minor axis of the ellipse is parallel to the shear loading direction, while the variations of the major axis (in the direction orthogonal to shear) can be neglected.

Defining ℓ_∥i and ℓ_⊥i the sizes of an elliptic microjunction along and orthogonal to shear, respectively, we can define its area as $A_i=\frac{\pi }{4}{\ell }_{\parallel i}{\ell }_{\perp i}$. Following Mindlin (1949), the stiffness of such an elliptic contact along the shear direction is, assuming no-slip contact conditions:

$\begin{array}{ll}{k}_i=\frac{\frac{\pi }{2}{\ell }_{\perp i}E}{(1+\nu )\left[\text{K}(e)-\frac{\nu }{e^2}(\text{K}(e)-\text{E}(e)\right)]},& (2)\end{array}$

with E and ν being the Young's modulus and Poisson's ratio of the material that constitutes the microjunctions, $e=\sqrt{1-\frac{{\ell }_{\parallel i}^2}{{\ell }_{\perp i}^2}}$ is the excentricity of the junction, K and E are the elliptic integrals of the first and second type, respectively. Note that assuming that microjunctions are elliptic is the simplest increment of realism compared to a circular assumption, in order to account for the complex shapes observed for microjunctions in the experiments.

Assuming that each microjunction is initially circular, we can define the common initial value, ℓ_0i, of ℓ_⊥i and ℓ_∥i from its initial individual area A_0i as ${\ell }_{0i}=\sqrt{\frac{4A_{0i}}{\pi }}$. As already mentioned, ℓ_⊥i varies negligibly under shear, so we will consider that ℓ_⊥i = ℓ_0i at all times. The evolution of each ℓ_∥i is then deduced from the shear-induced area reduction reported in Sahli et al. (2018):

$\begin{array}{ll}{A}_i=A_{0i}-{\alpha }_b\frac{1}{A_{0i}^p}{f}_i^2,& (3)\end{array}$

with α_b and the exponent p two constant parameters of the model. The size of junction i along the shear direction is thus simply ${\ell }_{\parallel i}=\frac{4A_i}{\pi {\ell }_{0i}}$. Note that there is currently no rigourous contact mechanics theory for the evolution of the shear stiffness of a sheared sphere/plane contact that would incorporate anisotropic contact area reduction. Here, such a behavior is approximated at all times by the combination of Equation (2), which is valid under no-slip assumption, and of Equation (3), which was empirically found at macroscale. Doing so, we assume that Equation (3) also applies at microscale, as suggested by the existence of common values of α_b and p for both the macro- and micro-scales (Sahli et al., 2018).

For each microjunction, Equation (3) is used from the beginning of the experiment, when f_i assumed to be 0, up to when the junction first starts to slip (when f_i = f_si). At that instant, A_i takes the value $A_{si}=A_{0i}-{\alpha }_b\frac{1}{A_{0i}^p}{f}_{si}^2$. For later times, based on the observation of the typical behavior of A_i during the experiments of Sahli et al. (2018) (see Figure 2), we assume that A_i always remains equal to A_si.

FIGURE 2

Figure 2. Left axis: time evolution of the contact area A_i of five typical microjunctions in the experiment at P = 4.01 N. Right axis: concurrent evolution of the macroscopic tangential force Q.

In contrast, the force resisted by a microjunction can vary with time after the first onset of slipping. When the slider's velocity, Ẋ(t), gets smaller than a minimum value Ẋ_min = c_min × v, with c_min a scalar parameter, we assume that all the slipping contacts will repin, with a position x_i = X−ϵf_si/k_i. Doing so, at repinning, there is no force discontinuity, as the repinning force k_i(X−x_i) is equal to the slipping force ϵf_si.

Following Sahli et al. (2018), the threshold force f_si at which the microjunction starts to slip is assumed to be proportional to its area at the same instant, i.e., f_si = σA_si, with σ the frictional shear strength of the contact.

The algorithm used to solve the model numerically is provided in Appendix 1.

2.2. Experimental Constraints

In order to quantitatively compare the model results with the multicontact experiments reported in Sahli et al. (2018, 2019), we need to feed the model with parameter values based on the measurements. In the following list, we first provide all constant parameter values that are directly accessible experimentally, with the error bars when relevant.

• M = M₀ + M₁, where M₀ = 100 g is the mass of the slider and M₁ = 0, 55, 111, 215, 308, 552 g an additional mass, for the six experiments performed. All masses are given at ±1g.

• k_L = 9,200 ± 200 N/m (Sahli et al., 2018).

• v = 0.1 mm/s.

• E = 1.6 ± 0.1 MPa (Sahli et al., 2018). This value and the associated error bar are the mean value and standard deviation over 32 estimates using 5 different spherical PDMS samples prepared in the same conditions.

• ν is assumed to be equal to 0.5, as is classically done for elastomers.

• the individual values of the initial areas of all microjunctions, A_0i, are extracted from the initial image (for Q = 0), segmented as described in Sahli et al. (2018). The fact that all microjunctions have different areas is the result of the random nature of the elastomer surface topography (see a typical Power Spectrum Density in Supplemental Material of reference Sahli et al., 2019) and of its elastic contact interaction with the rigid glass plate.

• N is also extracted form the same segmented image.

• σ_exp = 0.23 ± 0.02 MPa (Sahli et al., 2018), is the experimental value of the frictional shear strength of the interface, determined from a linear fit of (A_s, Q_s) for the 6 experiments. Q_s is the macroscopic static friction (peak) force and A_s is the total area of real contact at the same instant. We will discuss below how the value of σ in the model is related to σ_exp.

There are three model parameters which cannot be directly measured in experiments: η, ϵ, and c_min.

η is introduced to enable energy dissipation in the system, thus avoiding spurious oscillations of the slider. However, the value of η should not be too large, because it would prevent the possibility of stick-slip in the model. We found that stick-slip exists up to η between 180 and 200, but for those large values, the initial stick-slip cycles are significantly different from the experiments. In practice, we found that

$\begin{array}{ll}\eta =100& (4)\end{array}$

is a good compromise between oscillation reduction and a reasonable reproduction of the stick-slip sequence. The results are rather insensitive to the precise value of η, since η = 50 gives virtually identical results.

ϵ has a leading order control on the amplitude and period of the tangential force fluctuations during stick-slip. Systematic tests of the model for various values of ϵ led us to choose

$\begin{array}{ll}ϵ=0.90.& (5)\end{array}$

In particular, this value is sufficiently small to enable stick-slip for all six normal loads (as observed experimentally), while reproducing reasonably well the amplitude and period of the stick-slip sequences in all cases.

Note that in the model, if there was no stick-slip, the steady-state sliding friction force would be $\sum _{i=1}^Nϵ\sigma A_{si}$ (all microjunctions are in their slipping state). In order for this value to match the macroscopically measured value Q_s = σ_expA_s, one has to impose that

$\begin{array}{ll}\sigma =\frac{{\sigma }_{exp}}{ϵ},& (6)\end{array}$

and this is what we do in the following.

Our tests showed that the value of c_min has no impact on the results as far as it is sufficiently small. For instance, simulations with c_min=10⁻⁵ are essentially undistinguishable from those using 0.01. The reason is that, when |Ẋ(t)| crosses the value c_min × v, the velocity drop is so fast that the time at which the crossing occurs is almost independent on the value of c_min. In our calculations, we will use

$\begin{array}{ll}{c}_{min}=0.01.& (7)\end{array}$

Extracting values for p and α_b in Equation (3) requires fitting the power law relationship between the individual area reduction parameters, α_i, and the initial areas, A_0i, presented as purple squares in Figure 3 of Sahli et al. (2018). Such a fitting is actually difficult due to the large dispersion of the data, as can be inferred from the large difference in total area decay of microjunctions 1–3 in Figure 2, although they start with almost identical areas. A fit letting both α_b and p as fitting parameters gives 95% confidence error bars as large as 600% for the optimum value for α_b, which is not a viable option. We then tried to fix the value of p and fit the data with α_b being the only fitting parameter. We found that the quality of the fit (quantified by its R² value) was essentially independent of p (as long as it is not too different from the value 3/2 proposed in Sahli et al., 2018), preventing any objective choice of p.

Based on those observations regarding the determination of p and α_b from experimental data, in our model studies we decided to fix p and, for each value of p, we determined the value of α_b that gives the best agreement between the area decay predicted by the model and that measured in the experiments. To do that, we fitted both the experimental and model version of the curve A(Q) by a quadratic function of the form $A=A_0-\alpha Q^2$. A₀ being the same in the model as in the experiment (because $A_0=\sum _{i=1}^N{A}_{0i}$), the fitting procedure enables identification of an α_b which provides an exact match between the two quadratic decays. Importantly, we found that, for all tested values of p close to 3/2 (the value suggested in Sahli et al., 2018), the model results (when using the corresponding fitted α_b) were almost undistinguishable. So, in practice, we chose p = 3/2, for which the model studies give an optimal α_b = 0.45 10⁻¹⁵ m⁵/N² for the experiment with the smallest normal load, and α_b = 1.00 10⁻¹⁵ m⁵/N² for the experiment with the largest normal load. We then adopted the average value between both, α_b = 0.725 10⁻¹⁵ m⁵/N², as a constant to be used for all experiments.

3. Results: Quantitative Comparison

We run the model of section 2.1 with the parameter values described in section 2.2, for the six different PDMS/glass multicontact experiments reported in Sahli et al. (2018). Figure 3 compares, for two different normal loads, the measured time evolution of the area of real contact and tangential force to their corresponding model predictions. Note that the initial real contact area is essentially proportional to the normal load, as widely discussed in the contact mechanics modeling literature [see e.g., the reviews (Persson et al., 2005; Vakis et al., 2018)] and confirmed in the experiments discussed here (see Figure S2 of Sahli et al., 2018). To facilitate comparison between model prediction and measurement, the time origin of the experimental data has been offset by the amount necessary to superimpose the measured and predicted force curves in the central portion of their initial increase. Note that the initial non-linear increase of the measured force is due to the non-vanishing bending stiffness of the steel wire used to pull the slider, when it first bends around a pulley before a significant tension arises along the wire. The apparent difference between the measured and predicted values of the initial area of real contact is due to the above mentioned time offset: the initial predicted value exactly corresponds to the measured value from the first image, but the latter image now corresponds to a negative time and is thus not shown in the figure. The observed difference is of the order of the area measurement noise, presumably due to temporal fluctuations in the illumination and noise in the camera's sensor.

FIGURE 3

Figure 3. Direct comparison between measurements and model predictions, for two typical experiments with either P = 1.53 N (A) or P = 3.10 N (B). Time evolutions of the measured (dashed, black) and predicted (solid, blue) area of real contact, and of the measured (dashed, magenta) and predicted (solid, red) tangential force.

Figure 4 then shows, for all normal loads, the evolution of the area of real contact as a function of the tangential force, for both the measured and predicted data. This figure is similar to Figure 2A in Sahli et al. (2018), but shows all measurements points rather than just 1 of 130. Note that stick-slip is responsible for the accumulation of nearly horizontal cycles close to the minimum area/maximum force point of each curve. Also note that the model forces can transiently exceed the value σ_expA, but always remain smaller than $\sigma A=\frac{{\sigma }_{exp}}{ϵ}A$, as expected.

FIGURE 4

Figure 4. Real area of contact vs. tangential force for the six experiments. Black (blue) curves show the measured (predicted) data. The red line has slope $\frac{1}{{\sigma }_{exp}}$ and passes through the origin.

4. Discussion

Although other combinations of model ingredients may have been proposed, we believe that our model incorporates all of the currently available knowledge on the system that we tried to reproduce. As such, it can be seen as the most comprehensive independent asperities model of shear multicontacts so far, to be used for deterministic comparison with the experiments of Sahli et al. (2018, 2019).

Most of the model parameters (M, k_L, v, E, ν, A_0i, N, σ_exp) take their value directly from the measurements. Three adjustable parameters have been systematically varied to choose the most relevant value: c_min has no effect on the results, while η and ϵ have been adjusted to reproduce at best the stick-slip regime. Ideally, p and α_b should not be adjustable, but the dispersion in the experimental estimates of α_i is such that their values were not sufficiently constrained. In practice, the value of p was chosen equal to the one suggested from experiments incorporating not only microjunctions within multicontacts, but also millimetric smooth sphere/plane individual contacts (Sahli et al., 2018). The value of α_b was then adjusted to best match the overall decay of real contact area during the incipient loading of the interface.

With those values, the time evolution of the tangential load Q is quite well-reproduced (see Figure 3). In particular, the slope of the incipient loading is correct, which suggests that the stiffness used for the individual microjunctions is also correct. In contrast, the time evolution of the real contact area is not satisfactory. Of course, the total amplitude of the real area decay, from the initial contact to macroscopic sliding, is correct, because we start from the measured initial value ($\sum _{i=1}^N{A}_{i0}$), and we adjusted α_b to get the correct final value. So we argue that the quality of the comparison between the model and experimental results can only be assessed through the shape of the real area decay. And as can be seen from Figure 4, while the shape of the experimental curves A(Q) is essentially quadratic, that of the model curves is much more linear [except from the very beginning, when all microjunctions are pinned and thus decay quadratically according to Equation (3)]. We emphasize that this quasi-linear shape is a very robust feature of our model, because we found that the predictions are essentially unaffected by changes in the model assumptions (elliptic vs. circular microjunctions, Equation (3) applied at all times or only before the first depinning event) and in the parameter values (for values of η and ϵ enabling stick-slip or not).

The shape of the curve A(Q) results from a sum of a large number (N) of complex individual behaviors (non-linear area decay while pinned, constant area while slipping) with distributed parameters (initial area, stiffness, threshold), and is therefore unlikely amenable to a simple explanation. We can however mention an instructive particular case where all microjunctions would have the same initial area. In those (unrealistic) conditions, all microjunctions would behave identically when submitted to a common displacement X and thus depin at the same instant. The total area decay would be the sum of N identical quadratic decays, and thus be itself quadratic with the total shear load, until macroscopic sliding. With those specific (but wrong) initial conditions, we would recover a macroscopic area curve with the correct quadratic shape and a simple adjustment of the value of α_b would allow us to provide a good matching with the measurements. This example illustrates the major influence of the distribution of initial areas on the final shape of A(Q). It also clearly shows that the fact that we did not succeed in reproducing a quadratic area decay is not a generic problem of our model, but partly relates to the initial conditions (through the A_0i) imposed by the experimental dataset.

Could there be other reasons for the failure of our model to reproduce the evolution of the real contact area? The main model ingredient responsible for this evolution is Equation (3). The first possibility is that, despite the evidence brought in Sahli et al. (2018, 2019), the anisotropic area reduction under shear would not follow a single behavior law at all junction scales, from millimeter- to micrometer-sized junctions. This possibility is indeed suggested by a recent adhesion-based model of sheared sphere/plane junctions (Papangelo et al., 2019), where the authors find that the exponent p varies systematically, for a given sphere, with the normal load applied, and hence the initial area. Here we did not try to apply the model of Papangelo et al. (2019), because it would require the knowledge of the characteristic radius of curvature of, and normal load on, each individual microjunction. In contrast, experimental measurements only provide a combination of both quantities, through the area of the microjunction.

We now argue that the solution to the failure of our model will presumably be much more complex than a mere improvement of the form of Equation (3). The problem may very well be that the predicted individual force f_i is significantly different from the one that really applies on the microjunction. This is substantiated by Figure 2 which shows the time evolution of the contact area of various microjunctions. Two of them (4 and 5) were selected to show that the time window over which the area decay occurs can be very different: microjunction 4 has not started to shrink yet when the decay of microjunction 5 is already complete. This observation suggests that the individual tangential forces f₄ and f₅ evolve very differently during the experiment, although they have very similar initial areas and are submitted to the same tangential displacement by the glass substrate. We speculate that such a difference may be the result of elastic interactions between microjunctions, with junctions in a crowed neighborhood evolving differently from those far from neighboring junctions¹.

Those interactions are ignored in our model of independent microjunctions. We thus believe that, in order for an asperity model to have a chance to quantitatively match experiments like those considered here, tangential elastic interactions must be accounted for to describe the shear behavior of individual microjunctions. Such improved models may incorporate those tangential interactions in ways similar to models already developed for the normal interactions during normal loading of rough surfaces (see e.g., Ciavarella et al., 2006; Afferrante et al., 2012) or for the friction of multicontacts (see e.g., Braun and Scheibert, 2014; Trømborg et al., 2014; Braun and Peyrard, 2018).

5. Conclusion

We developed an independent asperity model for the incipient shear loading and onset of sliding of dry multicontact interfaces between a rough elastic solid and a smooth rigid surface. We used it to attempt the first deterministic comparison with experiments which, in addition to the macroscopic loads and displacements, also considers the individual areas of the many microjunctions forming the interface.

The main outcome is that, although we did our best to incorporate experimentally-based behavior laws, parameter values and initial conditions into the model, it fails to quantitatively reproduce the measurements of Sahli et al. (2018, 2019). Based on observations at the microjunction scale, we suggest that an interesting starting point for future attempts to improve the quantitative deterministic comparison between asperity models and experiments, may be to incorporate a description of the tangential elastic interactions between microjunctions.

Nevertheless, we anticipate that asperity-based friction models, although accounting for tangential elastic interactions, may suffer from the same limitations as their counterparts for purely normal contact (in particular the difficulty to define asperities when a continuum of length scales is involved in the topography, see e.g., Müser et al., 2017; Vakis et al., 2018), and may still be unsuccessful to quantitatively match friction experiments. We thus urge for the concurrent development of continuum models suitable to reproduce friction experiments like those of Sahli et al. (2018, 2019).

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

JS and MP designed and analyzed the results of the model. MP developed the numerical implementation of the model, and produced the model data shown in Figures 3, 4. RS performed the experiments and produced the experimental data shown in Figures 1–4. JS wrote the manuscript. All authors commented on the manuscript.

Funding

This work was supported by LABEX MANUTECH-SISE (ANR-10-LABX-0075) of Université de Lyon, within the program Investissements d'Avenir (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). It has been supported by CNRS-Ukraine PICS Grant No. 7422.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Oleg M. Braun, who initiated the work, but passed away before its completion.

Footnotes

1. ^The slight initial increase in area of junction 4 in Figure 2 may be due not only to elastic interactions but also to a slight aging due to viscous creep.

References

Acito, V., Ciavarella, M., Prevost, A. M., and Chateauminois, A. (2019). Adhesive contact of model randomly rough rubber surfaces. Tribol. Lett. 67:54. doi: 10.1007/s11249-019-1164-9

CrossRef Full Text | Google Scholar

Afferrante, L., Carbone, G., and Demelio, G. (2012). Interacting and coalescing Hertzian asperities: a new multiasperity contact model. Wear 278–279, 28–33. doi: 10.1016/j.wear.2011.12.013

CrossRef Full Text | Google Scholar

Braun, O. M., and Peyrard, M. (2008). Modeling friction on a mesoscale: master equation for the earthquakelike model. Phys. Rev. Lett. 100:125501. doi: 10.1103/PhysRevLett.100.125501

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, O. M., and Peyrard, M. (2018). Seismic quiescence in a frictional earthquake model. Geophys. J. Int. 213, 676–683. doi: 10.1093/gji/ggy008

CrossRef Full Text | Google Scholar

Braun, O. M., and Röder, J. (2002). Transition from stick-slip to smooth sliding: an earthquakelike model. Phys. Rev. Lett. 88:096102. doi: 10.1103/PhysRevLett.88.096102

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, O. M., and Scheibert, J. (2014). Propagation length of self-healing slip pulses at the onset of sliding: a toy model. Tribol. Lett. 56, 553–562. doi: 10.1007/s11249-014-0432-y

CrossRef Full Text | Google Scholar

Carnahan, B., Luther, H. A., and Wilkes, J. O. (1969). Applied Numerical Methods. New-York, NY: Wiley.

Google Scholar

Ciavarella, M., Delfine, V., and Demelio, G. (2006). A “re-vitalized” Greenwood and Williamson model of elastic contact between fractal surfaces. J. Mech. Phys. Solids 54, 2569–2591. doi: 10.1016/j.jmps.2006.05.006

CrossRef Full Text | Google Scholar

Costagliola, G., Bosia, F., and Pugno, N. M. (2018). A 2-D model for friction of complex anisotropic surfaces. J. Mech. Phys. Solids 112, 50–65. doi: 10.1016/j.jmps.2017.11.015

CrossRef Full Text | Google Scholar

Greenwood, J. A., and Williamson, J. B. P. (1966). Contact of nominally flat surfaces. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 295, 300–319.

Google Scholar

Johnson, K. L. (1987). Contact Mechanics. Cambridge: Cambridge University Press.

Google Scholar

Le Bot, A., Scheibert, J., Vasko, A. A., and Braun, O. M. (2019). Relaxation tribometry: a generic method to identify the nature of contact forces. Tribol. Lett. 67:53. doi: 10.1007/s11249-019-1168-5

CrossRef Full Text | Google Scholar

McGhee, A. J., Pitenis, A. A., Bennett, A. I., Harris, K. L., Schulze, K. D., Urueña, J. M., et al. (2017). Contact and deformation of randomly rough surfaces with varying root-mean-square gradient. Tribol. Lett. 65:157. doi: 10.1007/s11249-017-0942-5

CrossRef Full Text | Google Scholar

Mergel, J. C., Sahli, R., Scheibert, J., and Sauer, R. A. (2019). Continuum contact models for coupled adhesion and friction. J. Adhes. 95, 1101–1133. doi: 10.1080/00218464.2018.1479258

CrossRef Full Text | Google Scholar

Mindlin, R. D. (1949). Compliance of elastic bodies in contact. J. Appl. Mech. ASME 16, 259–268.

Google Scholar

Müser, M. H., Dapp, W. B., Bugnicourt, R., Sainsot, P., Lesaffre, N., Lubrecht, T. A., et al. (2017). Meeting the contact-mechanics challenge. Tribol. Lett. 65:118. doi: 10.1007/s11249-017-0900-2

CrossRef Full Text | Google Scholar

Papangelo, A., Scheibert, J., Sahli, R., Pallares, G., and Ciavarella, M. (2019). Shear-induced contact area anisotropy explained by a fracture mechanics model. Phys. Rev. E 99:053005. doi: 10.1103/PhysRevE.99.053005

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastewka, L., and Robbins, M. O. (2014). Contact between rough surfaces and a criterion for macroscopic adhesion. Proc. Natl. Acad. Sci. U.S.A. 111, 3298–3303. doi: 10.1073/pnas.1320846111

PubMed Abstract | CrossRef Full Text | Google Scholar

Persson, B. N. J. (2001). Theory of rubber friction and contact mechanics. J. Chem. Phys. 115, 3840–3861. doi: 10.1063/1.1388626

CrossRef Full Text | Google Scholar

Persson, B. N. J., Albohr, O., Tartaglino, U., Volokitin, A. I., and Tosatti, E. (2005). On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion. J. Phys. Condens. Matter 17, R1–R62. doi: 10.1088/0953-8984/17/1/R01

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponthus, N., Scheibert, J., Thøgersen, K., Malthe-Sørenssen, A., and Perret-Liaudet, J. (2019). Statistics of the separation between sliding rigid rough surfaces: simulations and extreme value theory approach. Phys. Rev. E 99:023004. doi: 10.1103/PhysRevE.99.023004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahli, R., Pallares, G., Ducottet, C., Ben Ali, I. E., Al Akhrass, S., Guibert, M., et al. (2018). Evolution of real contact area under shear and the value of static friction of soft materials. Proc. Natl. Acad. Sci. U.S.A. 115, 471–476. doi: 10.1073/pnas.1706434115

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahli, R., Pallares, G., Papangelo, A., Ciavarella, M., Ducottet, C., Ponthus, N., et al. (2019). Shear-induced anisotropy in rough elastomer contact. Phys. Rev. Lett. 122:214301. doi: 10.1103/PhysRevLett.122.214301

PubMed Abstract | CrossRef Full Text | Google Scholar

Thøgersen, K., Trømborg, J. K., Sveinsson, H. A., Malthe-Sørenssen, A., and Scheibert, J. (2014). History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework. Phys. Rev. E 89:052401. doi: 10.1103/PhysRevE.89.052401

PubMed Abstract | CrossRef Full Text | Google Scholar

Trømborg, J. K., Sveinsson, H. A., Scheibert, J., Thøgersen, K., Amundsen, D. S., and Malthe-Sørenssen, A. (2014). Slow slip and the transition from fast to slow fronts in the rupture of frictional interfaces. Proc. Natl. Acad. Sci. U.S.A. 111, 8764–8769. doi: 10.1073/pnas.1321752111

PubMed Abstract | CrossRef Full Text | Google Scholar

Vakis, A. I., Yastrebov, V. A., Scheibert, J., Nicola, L., Dini, D., Minfray, C., et al. (2018). Modeling and simulation in tribology across scales: an overview. Tribol. Int. 125, 169–199. doi: 10.1016/j.triboint.2018.02.005

CrossRef Full Text | Google Scholar

Violano, G., and Afferrante, L. (2019). Contact of rough surfaces: modeling adhesion in advanced multiasperity models. Proc. Inst. Mech. Eng. J. Eng. Tribol. 233, 1585–1593. doi: 10.1177/1350650119838669

CrossRef Full Text | Google Scholar

Yastrebov, V. A., Anciaux, G., and Molinari, J.-F. (2017). The role of the roughness spectral breadth in elastic contact of rough surfaces. J. Mech. Phys. Solids 107, 469–493. doi: 10.1016/j.jmps.2017.07.016

CrossRef Full Text | Google Scholar

Appendix 1: Integration of the Equation of Motion of the Slider

Integration Scheme for the Differential Equation Giving X(t)

Equation (1) is integrated by a second order Leapfrog algorithm.

$\begin{array}{ll}\gamma (t)=\ddot{X}(t)=\frac{1}{M}\left[k_L(vt-X)-{\displaystyle \sum _{i=1}^N}{f}_i\right]-\eta \dot{X}& (\mathrm{A1})\end{array}$

is the acceleration of the slider at the beginning of an integration step. γ(t = 0) has to be calculated at the start of the simulation. In our case, we start with X(t = 0) = 0 and Ẋ(t = 0) = 0 so that γ(t = 0) = 0.

Let dt be the time step. The algorithm first updates X(t) according to:

$\begin{array}{ll}X(t+dt)=X(t)+dt\dot{X}(t)+\frac{1}{2}d{t}^2\gamma (t).& (\mathrm{A2})\end{array}$

With this updated value of X we calculate the new acceleration γ(t + dt). In the second stage of the algorithm Ẋ is updated according to

$\begin{array}{ll}\dot{X}(t+dt)=\dot{X}(t)+\frac{1}{2}dt(\gamma (t)+\gamma (t+dt)),& (\mathrm{A3})\end{array}$

and then everything is ready for the next step starting at t+dt.

To properly select dt we calculate the elasticity constant k_i for each contact at the start of the simulation. The angular frequency of oscillation of the slider of mass M under the total force of the contacts is ${\Omega }_0=\sqrt{(\sum _{i=1}^N{k}_i)/M}$ and the corresponding period is T₀ = 2π/Ω₀. The calculation uses dt = T₀/C_t with $C_t=10^4$. The exact value of dt depends on the experiment, but a typical value is dt = 1μs. This value is sufficiently small with respect to all the frequencies entering in the dynamics of the slider so that even a second order algorithm gives a very good numerical accuracy. We have however run some calculations with a 4th order Runge-Kutta method (Carnahan et al., 1969), which is significantly slower, but has errors that decay as dt⁴, to check the accuracy of our calculations.

Algorithm of the Subroutine Which Calculates γ

To compute γ(t+dt), X(t+dt) and x_i are known. The main point is to compute all the forces f_i on the junctions. The state of each junction is recorded with two flags: θ_i records its instantaneous state, θ_i = 1 for a pinned junction, θ_i = 0 for a slipping junction, and h_i keeps track of its history, h_i = 1 for a junction which has never been slipping switches to h_i = 0 the first time the junction starts to slip, when f_i ≥ f_si.

The program scans all the junctions and performs the following steps:

• Compute f_i for each junction

⋆ if h_i = 1 the area of the junction depends on f_i according to Equation (3), which determines ℓ_∥i and then k_iA_i(f_i) according to Equation (2). Thus

$\begin{array}{ll}{f}_i=k_i{[A}_i(f_i)](X(t+dt)-x_i)& (\mathrm{A4})\end{array}$

gives an equation for f_i. It is too complex to be solved analytically. We solve it by an iterative process, using a dichotomy method starting from the value of f_i from the previous step. Once f_i is known we update A_i(f_i), k_i(A_i) and f_si = σA_i for further steps.

⋆ if h_i = 0

- if θ_i = 1 the junction is pinned but A_i is fixed, as well as k_i(A_i), and they are known from previous iterations so that f_i = k_iX(t)−x_i.

- if θ_i = 0 the junction is slipping. f_i = ϵf_si.

• Check for transitions in the junction state

⋆ if θ_i = 1 (pinned junction) then if f_i ≥ f_si the junction starts to slip so that θ_i switches to 0, f_i = ϵf_si. h_i switches to 0 if it was still equal to 1.

⋆ if θ_i = 0 (slipping junction), if |Ẋ(t)| < c_min × v the junction repins, θ_i switches to 1 and we set x_i = X−ϵf_si/k_i so that the junction starts in the pinned state with f_i = ϵf_si.

Once all junctions have been scanned and all f_i are determined, we can compute γ from Equation (A1).