Statistics of the separation between sliding rigid rough surfaces: Simulations and extreme value theory approach
Nicolas Ponthus, Julien Scheibert, Kjetil Thøgersen, Anders Malthe-Sørenssen, Joël Perret-Liaudet
aa r X i v : . [ phy s i c s . c l a ss - ph ] A p r Statistics of the separation between sliding rigid rough surfaces:Simulations and extreme value theory approach
Nicolas Ponthus, Julien Scheibert, Kjetil Thøgersen, Anders Malthe-Sørenssen, and Jo¨el Perret-Liaudet Univ Lyon, Ecole Centrale de Lyon, ENISE, ENTPE, CNRS,Laboratoire de Tribologie et Dynamique des Syst`emes LTDS UMR5513, F-69134, Ecully, France Physics of Geological Processes, The NJORD Centre,Department of Geosciences, University of Oslo, Norway Department of Physics, University of Oslo, 0316 Oslo, Norway (Dated: April 30, 2019)When a rigid rough solid slides on a rigid rough surface, it experiences a random motion inthe direction normal to the average contact plane. Here, through simulations of the separation atsingle-point contact between self-affine topographies, we characterize the statistical and spectralproperties of this normal motion. In particular, its rms amplitude is much smaller than that ofthe equivalent roughness of the two topographies, and depends on the ratio of the slider’s lateralsize over a characteristic wavelength of the topography. In addition, due to the non-linearity ofthe sliding contact process, the normal motion’s spectrum contains wavelengths smaller than thesmallest wavelength present in the underlying topographies. We show that the statistical propertiesof the normal motion’s amplitude are well captured by a simple analytic model based on the extremevalue theory framework, extending its applicability to sliding-contact-related topics.
I. INTRODUCTION
The interfacial separation d between the surfaces of twosolids brought close to one another is central to many in-terfacial processes. Those include attractive forces whenthe distance is small but finite (Van der Waals [1], elec-trostatic [2], Casimir forces [3]), repulsive elastic forceswhen the distance vanishes [4, 5], heat transfer and non-contact friction [6], electric conductivity [7] and perme-ability [7, 8]. The evaluation of d becomes difficult whenthe typical separation becomes of the order of the surfaceroughness, because the separation is now a random vari-able of the position along the interface. In this case, d often refers to the average separation, between the meanplanes of the two rough surfaces.In the particular case when the two rough surfacescome into contact, most of the literature has treated theirnormal approach (see e.g. [9, 10], and [11] for shear load-ing). For elastic bodies under sufficient compressive pres-sure, a so-called multi-contact is formed, made of myr-iad individual micro-contacts where mainly the highestantagonist asperities are involved in the actual contact.This situation is typical of elastomer contacts [11]. Theaverage separation between the two bodies is found to de-pend in particular on the ratio p/E ⋆ of the applied pres-sure p to composite elastic modulus E ⋆ and the spectralproperties of the topography [5, 12]. When p/E ⋆ tendstowards zero, i.e. when the pressure becomes very lowcompared to the material stiffness, and when the two sur-faces are brought in contact through a pure normal trans-lation, those two surfaces touch on only one point, whichis the first to come into contact. Such single-point con-tact situations, which are the focus of the present study,have previously been investigated in the context of theprecise measurement of dispersion forces [13, 14] or of thecontact of metallic surfaces under light load [15]. In such cases the measurable quantity is the separation of thetwo mean planes for single-point contact, d (see Fig. 1).Note that if the two solids are shifted one with respectto the other parallel to the contact plane, the measuredvalue of d will likely vary, because the single-point con-tact will involve a different couple of antagonist asper-ities. Such a sensitivity to details of the measurementprocedure is responsible for significant uncertainties inthe evaluation of d [13, 14]. d is also expected to varyas soon as the solids are slid one on another, because thepoint of contact will continuously change, and this is thephenomenon of interest in the following. FIG. 1. (a): Illustration of single-point contact, on the ex-ample of two 1D centered Gaussian white noises, Z and Z .The separation at single-point contact, d , is measured be-tween the mean heights of the two processes. (b) Probabilitydensity functions (pdf) of both processes. From now on, we will consider that, in such a weakly-loaded, single-point contact situation, one body (theslider, with a finite-sized area) is set to slide on the other(the track, having a larger area). Due to the randomnature of the antagonist topographies, the slider will ex-perience a roughness-induced motion in the direction nor-mal to the average contact plane ( z -displacement), d ( u ),with u the tangential displacement ( x -displacement) ofthe slider, as shown on Fig. 1. In the following, d ( u )is referred to as the normal motion. If the sliding veloc-ity was high enough, the slider could loose contact withthe track and enter a bouncing regime ([15–18]). In thefollowing, we only consider slow sliding, in which suchinertia effects can be neglected. In those quasi-staticconditions, the time dependence of the normal motionis irrelevant and the quantity of interest is d ( u ). Tocharacterize this quantity, we perform direct numericalsimulations of the single-point contact between slidingrigid rough surfaces, as described in section II.From the illustration of Fig. 1, it is natural to inter-pret the single-point contact sliding process as a geomet-rical filtering in which the input signals are the two an-tagonist topographies, Z and Z , and the output sig-nal is the roughness-induced normal motion, d ( u ). Inparticular, one expects that the broader the probabilitydensity function of the topographies, the larger the aver-age single-point contact separation. One also expects thespectral properties of d ( u ) to be dependent on in-planefeatures of the topographies, like their spectral contents.This is why our simulations explore a variety of powerspectrum densities (PSD) of the contacting topographies(section II). In section III, we characterize in details therelationship between the properties of the topographiesand that of the resulting normal motion d ( u ).The height of single-point contact being directly re-lated to the altitude of the highest asperities of the an-tagonist surfaces, it is tempting to use the concepts ofextreme value theory (EVT) [19–21] to estimate the sta-tistical properties of d . EVT has been extensively usedin various fields [22], including rupture in disordered me-dia [23], risk in finance or insurance [20], or catastrophicnatural events (preface of [21]). EVT predicts the proba-bility distribution of rare events, and is used in section IVto predict the distribution of the maximum height of thetopographies and thus of d . Those predictions are quan-titatively compared to the simulation results and usedto discuss the applicability of EVT to sliding-contact-related topics. II. DIRECT SIMULATIONS: METHODSA. Properties of the topographies
To characterize the properties of the separation atsingle-point contact between sliding surfaces, d ( u ), weperformed direct simulations of a rough square slider(surface L × L ) moving quasi-statically along a roughtrack (surface L × L , L > L ), and touching it in asingle point for each of the successive positions u of theslider. The two rotations of the slider around the in-plane axis are forbidden and its translation along thetrack is imposed. Its only free motion is that along theout-of-plane, z -axis. The slider and track have the samestatistical roughness properties. Two-dimensional (2D) Gaussian topographies, z , withvarious spectral properties have been generated, fromtheir 2D power spectrum density (PSD). Assuming thatthe topographies are isotropic, they are fully character-ized by the radial profile of their PSD, S zz ( k r ), with k r the radial wave number. Knowledge of the PSD allowsone to calculate a variety of useful estimators of the to-pographies’ properties, among which its rms roughness, R q , from: R q = M , (1)and its central wavelength, λ , from: λ = 12 r M M , (2)with the radial spectral moments M i defined by [24] M i =2 π R + ∞ k i +1 r S zz ( k r )d k r . FIG. 2. Sketch of the radial profiles S zz ( k r ) of the 2D PSDsconsidered for the antagonist topographies and definition ofthe corresponding parameters. We used realistic PSDs corresponding to self-similartopographies such as the one shown on Fig. 2. Such PSDcan be fully described using 4 parameters: S sets the am-plitude of the surface; k l , the low cut-off wave number,and k s , the high cut-off wave number, set the wave num-ber range over which the topographies are self-similar; − α is the slope of the self-similar part and is linked tothe fractal dimension. Indeed, α relates to the Hurst ex-ponent H through [25–27]: α = 2( H + 1). Note that,for a slider of size L , the lowest accessible wave numberis k L = πL . With this choice of PSD profile, the threefirst radial spectral moments M i can be calculated an-alytically (appendix A). Once injected in Eqs. 1 and 2,appendix A provides the explicit expressions of R q and λ .In summary, a given simulation corresponds to a givenset of 5 parameters: S , k l , k s , α and L . In practice,we will use the following equivalent set of 5 parameterswith a more intuitive physical meaning: R q , λ and L as three characteristic length scales and α and b = k l k s astwo shape descriptors of the PSD radial profile. B. Numerical topography generation
The surfaces are represented numerically by a (2Θ +1) × (2Ψ+1) height matrix z . The location along a surfaceis identified by the vector x ij = ( x i , y j ) with i varyingfrom − Θ to Θ and j from − Ψ to Ψ such that x i = i ∆ x and y j = j ∆ y . Thus z ( x ij ) defines the altitude of thetopography at each point. The Fourier transform of z is F [ z ]( k θψ ) = A ( k θψ ) exp(i φ ( k θψ )) with A , the amplitude,and φ , the phase, two real functions. The wave vector is k θψ = ( k θ = θ (2Θ+1)∆ x , k ψ = ψ (2Ψ+1)∆ y ), with θ varyingfrom − Θ to Θ and ψ from − Ψ to Ψ. By setting: A ( n, m ) = A ( − n, − m ) , (3) φ ( n, m ) = − φ ( − n, − m ) , (4)we ensure that z is real, and reads (after inverse Fouriertransform of A ( k θψ ) exp(i φ ( k θψ ))): z ( x ij ) = 12Θ + 1 12Ψ + 1 X θ X ψ A ( k θψ ) cos( φ ( k θψ ) + k θψ · x ij ) (5)The amplitude A ( k θψ ) can be expressed as a function ofthe continuous PSD profile, S zz ( k r ), as: A ( k θψ ) = s (2Θ + 1)(2Ψ + 1) S zz ( | k θψ | )∆ x ∆ y . (6)In order to produce numerical topographies obeyingthe PSDs described in section II A, we use Eq. 5 in whichwe insert both Eq. 6 and phases φ randomly drawn from auniform law over [0 2 π [, yielding a Gaussian distributionof heights.Figure 3 represents four typical topographies obtainedfor various values of α and b , the lateral length of all pan-els corresponding to the same number of central wave-length, λ . One can see that the smaller b and α , thericher the spectral contents of the topography, with b having the strongest effect. The spectral bandwidth canbe quantified by a spreading parameter δ z = q − M M M inspired by [28–32]. δ z can vary between 0 and 1, with δ z being close to 0 for a narrow-band topography. Figure 4,which shows the evolution of δ as a function of α , forvarious b , confirms the trends illustrated in Fig. 3. C. Simulation parameters
The in-plane discretization ∆ x = ∆ y of the surfaceis chosen such that π x = 6 k s . This ensures that the -5-4-3-2-1012345 FIG. 3. Typical topographies generated, for various valuesof α and b . In-plane size in units of λ . Out-of-plane size(grayscale bar) in units of R q .FIG. 4. Spectral bandwidth δ z as a function of the shape-descriptors of the PSD considered here (Fig. 2), α and b . sinus corresponding to the largest wave number is well-resolved, with twelve points per wavelength in the spatialdomain. Sliding motion is simulated by moving the slideralong x with respect to the track, by one grid size at eachstep.For each dimension of the topographies (out-of- and in-plane) a reference length is chosen. For the out-of-planedimension, the rms roughness, R q , is chosen, while forthe in-plane dimension, we chose the central wavelength, λ . Note that, in our case of normal approach of rigidbodies, the in- and out-of-plane dimensions are uncou-pled. In particular, dilating only the in-plane dimensiondoes not affect the value of the normal separation, whiledilating only the out-of-plane dimension does not affectthe index of the topography points that are involved inthe single-contact. Simulations are thus defined by threedimensionless parameters: two for the PSD shape, b and α , and one for the slider size ˜ L = L/λ . The outputquantity is thus the dimensionless separation at single-point contact, ˜ d = d / √ R q , as a function of the di-mensionless sliding distance, ˜ u = u/λ . Note that, forthe contact between two statistically identical topogra-phies, each with an rms roughness R q , as is the casein the present study, the normalizing quantity for d , √ R q = R ∗ q , represents the equivalent rms roughness ofthe sum topography.The value of b is varied from 0 .
05 to 1. Notice that b = 1 is the case of a rectangular-shaped radial PSD,while for b =0, the topography would be purely self-similar. α is varied from 3 to 4, corresponding to a Hurstexponent varying from 0 . L is varied from 43 to 760. The tracksize is then ˜ L × ˜ L . Note that we have limited the rangeof variations of ˜ L to values such that (i) k L < k l so thatthere is a white noise part in the PSD and (ii) L = 5 L and L/ ∆ x is smaller than 17000, to keep topographymatrices computationally tractable. Our computationalresources allowed simulation of topographies with anycombination of parameters within the above-mentionedranges. An additional set of simulations has been per-formed with α = 4 and b = 0 .
46, which allowed us tovary ˜ L from 1.5 to 1389, and thus to explore more widelythe effect of the slider size on the roughness-induced nor-mal motion.For each random draw of phases, a different topogra-phy is generated, with the specified PSD. For each set ofparameters ( α , b and ˜ L ), several draws of topographiesare performed in order to get converged statistical resultsfor the separation at single-point contact, ˜ d . Tests haveshown that with 15 draws, the expected value of each ofthe three first statistical moments (mean, standard devi-ation and skewness) of ˜ d is measured to better than 5%accuracy.Finally, let us define the parameter N , which will beuseful in the following sections: N = 4 π ˜ L . (7) N represents the number of circular patches of diameter λ along the slider’s surface. For a narrow band process,it is close to the number of asperities on the slider’s sur-face. Here and in the following, the term asperity refersto any convex portion of the topography. III. DIRECT SIMULATIONS: RESULTS
On Fig. 5, an example of simulated separation atsingle-point contact, ˜ d (˜ u ), is plotted. On Fig. 6, typ-ical probability density functions (pdf) of ˜ d are shown.One can notice that h ˜ d i , the mean value of ˜ d , is largerthan 0 by several R ∗ q (typically 2 to 5, depending of thesimulated topogaphies). Note that d ( u ) (the distancebetween the two mean planes (see Fig. 1)) can be equal FIG. 5. Typical example of separation at single-point contact,˜ d (˜ u ). α = 4, b = 0 .
46, ˜ L = 327. Inset: zoom showing cusp-like features. to 0 only if one topography would be the exact comple-mentary of the other at position u . Also, the standard de-viation of ˜ d , σ ˜ d , is always found smaller than 1. Finally,the skewness of the distribution, sk ˜ d , is positive, due tothe fatter right tail of the pdf, implying that, unlike theunderlying topographies, the separation at single-pointcontact, d , is not a Gaussian process. FIG. 6. Typical pdf of the separation at single-point contactfor two different slider sizes, ˜ L = 327 ( N = 136118) and˜ L = 44 ( N = 2423). α = 3 . b = 0 .
46. Dashed lines: EVTmodel discussed in section IV.
More quantitatively, Fig. 7 shows the normalized meanvalue, standard deviation, and skewness, respectively h ˜ d i , σ ˜ d and sk ˜ d , for all simulations parameters used, asa function of the slider area, represented by the number N . Figure 7 clearly shows that the statistical moments ofthe dimensionless separation ˜ d only depend on N . Both h ˜ d i and sk ˜ d are increasing functions of N whereas σ ˜ d FIG. 7. (a) Mean value, (b) standard deviation and (c) skew-ness of the pdfs of ˜ d , for all simulations performed (circles).Error bars represent the standard deviation over 15 statisti-cally identical simulations. For clarity, not all error bars areplotted. Lines: various models discussed in section IV. Themarker size increases when b decreases. Lighter gray is forlower α . • , N , (cid:7) , (cid:4) are for L = 100, 250, 500 and 750, respec-tively. The two values of δ z used in Figs. 7a and 7b are themaximum and minimum values used in our simulations is a decreasing function. Note that the N -axis is logarith-mic, indicating that the variations with N are relativelyslow.We now describe the spectral contents of theroughness-induced normal motion of the slider. Theyare described by the power spectrum density of ˜ d (˜ u ), S ˜ d ˜ d (see Fig. 8 for a typical example). We find that theresulting PSDs are of the self-affine type, with a whitenoise part at low frequencies, and a power law decay ofexponent − α ⋆ at high frequencies. The crossover wavenumber is denoted ˜ k ⋆l . FIG. 8. Typical PSD of the separation at single-point contactand its empirical approximation (Eq. 8). α = 3 . b = 0 . L = 337. As can be seen on Fig. 8, ˜ d has non-vanishing spectralcontents for wave numbers higher than ˜ k s , the topogra-phies’ largest wave number. We measured the exponent − α ⋆ of the power-law decay of the PSD, for wave num-bers larger than ˜ k s , and found that it is always roughlyequal to -4. Then, to estimate the value of ˜ k ⋆l , we proposethe following empirical model for the PSD (see black linein Fig. 8): S ⋆ if ˜ k < ˜ k ⋆l , S ⋆ (cid:16) ˜ k ⋆l ˜ k (cid:17) if ˜ k > ˜ k ⋆l . (8)We then fit the value ˜ k ⋆l with the constraint that themoment of order 0 of the model PSD (which is the rms value of the signal) is equal to that of the simulated one,which amounts to impose that S ⋆ = q σ d + h ˜ d i ˜ k ⋆l . Anal-ysis of the dependence of ˜ k ⋆l with the simulation param-eters, α , b and ˜ L , for all the simulations performed, al-lowed us to find the following empirical expression:˜ k ⋆l ≈ f ( b, ˜ L ) = (cid:18) . · − b . + 0 . (cid:19) ˜ L . . (9)Figure 9 shows that Eq. 9 nicely predicts the value of ˜ k ⋆l obtained from the simulations. FIG. 9. Simulated ˜ k ⋆l versus its approximated expression, f ( b, ˜ L ) (Eq. 9). Solid line: equality line. The marker sizeincreases when b decreases. Lighter gray is for lower α . • , N , (cid:7) , (cid:4) are for L = 100, 250, 500 and 750, respectively. It is interesting to reformulate those results in terms ofthe central wavelength and the spectral bandwidth of theprocess ˜ d , which are generic estimators of a PSD, alsovalid in particular for non-self-affine PSDs. They are de-fined as ˜ λ ⋆ = q ˜ m ˜ m and δ ⋆ = q − ˜ m ˜ m ˜ m , respectively,with moments ˜ m i = R + ∞−∞ | ˜ k | i S ˜ d ˜ d (˜ k )d˜ k . Note that withsuch a definition [28–32], odd moments do not vanish.We investigated the relationship between the spectralparameters of ˜ d and those of the contacting topogra-phies and found the results shown in Fig. 10. First, δ ⋆ isa function of δ z only, through (Fig. 10a): δ ⋆ ≈ . δ z + 0 . . (10)Second, ˜ λ ⋆ depends on both δ z and ˜ L , through:˜ λ ⋆ ≈ f ( δ z , ˜ L ) = − . δ z + 1 . L . . (11)Figure 10b shows that this expression is a good approxi-mation of ˜ λ ⋆ , for all the simulations performed. Findingexplanations for Eqs. 9, 10 and 11 would be the subjectof an interesting future work. IV. DISCUSSIONA. Single-point contact as a geometrical filtering
The shape observed for the probability density func-tion of the interfacial separation (Fig. 6) can be under-stood as a geometrical filtering of the antagonist topogra-phies. This filtering process is expected to strongly de-pend on the size of the slider. In the limit of a point-like slider ( L vanishes), it will be able to follow exactlythe track’s topography, so that its normal motion will be FIG. 10. Spectral parameters of the normal motion ˜ d (˜ u ) as afunction of those of the contacting surfaces. (a) δ ⋆ vs δ z . (b)˜ λ ⋆ vs f ( δ z , ˜ L ) (Eq. 11). Solid line: equality line. The markersize increases when b decreases. Lighter gray denotes lower α . • , N , (cid:7) , (cid:4) are for L = 100, 250, 500 and 750, respectively. equal to that topography. In the case of a finite-sizedslider, the slider will not be able to penetrate into thevalleys of the track’s topography but will mainly slideon the highest asperities. Hence, the slider is expectedto successively explore the shape of different asperities:the ones with the smallest distance to the slider’s to-pography. The switching between asperities in contactis expected to be abrupt because the slider will intanta-neously stop following the shape of the previous asperityand start following the new one. This scenario is in per-fect agreement with the typical normal motion shown inFig. 5 (inset), in which one can identify cusp-like featuresat the local minima (when the slider switches asperities)and smooth maxima (when the slider follows the summitof one asperity).Those features of the normal motion are fully consis-tent with the observations made on the pdf and PSD of˜ d . The asymmetry between minima and maxima in thenormal motion explains the fatter tail of the pdf for largealtitudes, which are more probable than the small ampli-tudes (at the cusps), and thus explains why sk ˜ d > | sin ( x ) | is a Dirac comb withamplitude decreasing asymptotically as 1 /k , i.e. withan exponent close to the measured − α ⋆ ≃ −
4. We sug-gest that the presence of the cusps is the origin of the ob-served spectral enrichment (beyond k s ) of the roughness-induced normal motion. The observation that the stan-dard deviation of d is always smaller than that of thesum topography (Fig. 7b) is related to the fact that theslider cannot explore the lowest parts of the track’s to-pography, due to its finite size. Finally, the non-vanishingvalues of the mean of ˜ d are fully consistent with the factthat the slider only touches the highest asperities of thetrack (Fig. 7a).As noted above, the slider’s normal motion only dif-fers from the track’s topography if the slider has a non-vanishing size, which indicates that the geometrical fil-tering is intrinsically a finite-size effect. And indeed, thevarious statistical properties of the normalized normalmotion’s pdf only depend on N (Fig. 7). Larger slid-ers have a larger probability to touch a high asperity, sothat h ˜ d i is larger (Fig. 7a). Similarly, larger sliders pen-etrate less into the tracks’ valleys, so that σ ˜ d is smaller(Fig. 7b). B. Extreme value theory approach
Due to the large number of points used to represent thetopographies, the numerical simulations presented aboveare computationally expensive and require large RandomAccess Memory for the reverse Fourier transform opera-tion (for the largest simulations, we used 512Gb of RAMfor 3h30 on Bi-Xeon E5-2640v3 (16 core 2.6GHz)). Thusbeing able to predict the statistical properties of the sepa-ration at single-point contact, d , directly from the prop-erties of the topographies is highly desirable. Remember-ing that the slider, due to its finite size, can only get intocontact with the highest asperities of the track’s topog-raphy, it is tempting to investigate how the frameworkof extreme value theory (EVT, see e.g. [19–21]) can beapplied to the present single-point contact problem.We represent rough surfaces through N points whichare independent realizations of a centered Gaussian pro-cess. Consider two antagonist such topographies, z and z , with identical rms roughness R q , as sketched in Fig. 1.Their separation at single-point contact, d , is given by d = − min ( x i ,y j ) ( z ( x i , y j ) − z ( x i , y j )). z and z hav-ing symmetric distributions, (i) the distribution of theirdifference is then statistically equal to the distribution oftheir sum and (ii) the opposite of the minimum is statis-tically equivalent to the maximum. We can thus work onthe sum of the two topographies, z = z + z , which has acentered Gaussian distribution with a standard deviationequal to R ∗ q = √ R q , and examine d = max( z ).Let P be the cumulative density function (cdf) of z and p the associated probability density function (pdf).In our case: p ( z ) = 1 R ∗ q √ π exp − z R ∗ q ! , (12) P ( z ) = 12 zR ∗ q √ ! . (13)We will now follow the EVT approach described e.g. in [19]. The probability that the altitude at one pointof z is smaller than a value Y is given by P ( Y ). Theprobability that all N altitudes of z are smaller than Y ,meaning that Y is greater or equal to the highest point of z , is the cdf G ( Y ) = P ( Y ) N . Thus, the pdf correspond-ing to the fact that Y is the largest of the N altitudes,which is precisely the pdf of the sliders’ height, reads: g ( Y ) = G ′ ( Y ) = N p ( Y ) P ( Y ) N − . (14)Typical distributions g are shown in dashed lines onFig. 6 for two values of N , which is analogous to twodifferent sizes of the slider. It appears that those EVT-predicted distributions present similar qualitative fea-tures as the simulated distributions. In particular, themean height of the slider also increases with N , while thestandard deviation of the slider’s height also decreaseswith N . As in simulations, g is not a symmetric distri-bution, but has a fatter tail for large values of Y (positiveskewness), which is typical of extreme value statistics. Asshown in [19, 33], when N is increased, the distributionof the maximum of a variable tends toward either theFrechet, Gumbel or Weibull distribution depending onthe pdf of this variable. In particular, for a pdf with aright tail decaying faster than a power law, as is the casefor Gaussian distributions, the pdf of the maximum willtend toward the Gumbel distribution. As a consequence,we expect the Gumbel distribution to be the limiting caseof g for very large sliders ( N ≫ C. Comparison between methods
Once the topography’s pdf has been chosen to be Gaus-sian, the EVT prediction (Eq. 14) only depends on theparameter N . So, quantitative comparison between thepredictions of EVT and the numerical simulations onlyrelies on a relevant choice of N . Remembering that inEVT, we represent the topographies as a collection of N discrete, statistically independent values, one looks fora number related to the number of asperities present onthe simulated slider’s surface. For a 1D process, thisnumber can be given [29] by dividing the length of ob-servation, L , by the central wavelength, λ (Eq. 2), sothat N = Lλ . For the two-dimensional processes ob-served here, the same path of thought can be followedwith areas of diameter λ : N = L π λ = 4 π (cid:18) Lλ (cid:19) = 4 π ˜ L , (15)which justifies the prefactor used in Eq. 7. Note that in[13], a correlation length is used instead of λ to define N .Using this choice of N in the EVT approach, we over-plotted the analytical results of Eq. 14 on all panels ofFig. 7. This comparison shows a rather good quantita-tive agreement with our simulation results, confirmingthat the good match observed on Fig. 6 is actually truefor all the explored simulation parameters. Yet, an offsetexists on h ˜ d i between the EVT prediction and the nu-merical results. We interpret this offset as a side effectof the approximation of a continuous topography by aset of discrete points: while the mean plane of a single,continuous asperity always lies below its summit, in thecase of a point-like asperity, the mean plane has the ex-act same altitude as the summit itself. Hence, the meanplanes of two continuous topographies in contact are al-ways separated by a larger distance than those of twosets of discrete asperity summits. The exact offset be-tween the two situations depends on both the amplitudeand shape of the asperity. In our case, an empirical cor-rection of R ∗ q seems to correctly capture our simulationdata. Note that in Fig. 6, the analytical pdfs shown in-clude this correction, while the solid line in Fig. 7a doesnot.We emphasize that such an agreement is a priori non-trivial. First, the agreement quality significantly dependson the definition of N , suggesting that the arbitrary def-inition used (Eq. 15) is adequate. Second, the predic-tion is based on the EVT framework, which considerstopographies made of independent realizations of a Gaus-sian process. In constrast, the topographies used in thesimulations incorporate a finite correlation length, dueto the shape of the PSD used to generate them. We be-lieve that this difference is the main reason for the slightdiscrepancies observed in Fig. 7 between simulations andEVT predictions. To improve the agreement, one wouldneed to account for the deviations from EVT induced bya finite correlation length.This is what Preumont attempted in [29], on the prob-lem of finding the maximum value reached during a cer-tain time window by a correlated 1D Gaussian signal.Assuming that the successive extrema of the signal forma Markovian process, he was able to find an exact, butintricate expression for the pdf of this maximum value.By fitting this pdf with a Gumbel distribution, he wasable to identify semi-empirical expressions of its meanvalue and standard deviation, as a function of N and the spectral bandwidth δ of the process: h ˜ d i P reumont = p κ u N + γ √ κ α N , (16) σ ˜ d P reumont = π √ √ κ α N , (17) κ u = ( . − e − . δ ) if δ < . .
94 if δ ≥ . , (18) κ α = ( δ if δ < . .
05 if δ ≥ . , (19)with γ = 0 . √ ζ (3) π ∼ .
14 ( ζ is Riemann’s zeta function), indepen-dently of N (see Fig. 7c).Those semi-empirical expressions are overplotted onFig. 7 using the value of δ z for δ in Eqs. 18 - 19. Those ex-pressions appear to provide an excellent agreement withour simulation data, in particular they capture the cor-rect amplitude of h ˜ d i . Such an improvement of theagreement confirms that the discrepancies observed be-tween EVT and simulations are mainly due to the fi-nite correlation of the simulated topographies. Yet, hereagain, such a good agreement was not expected, sinceEqs. 16- 19 were obtained for 1D processes, while our sim-ulations use correlated 2D processes (the topographies). D. Relation to experiments
There are very few experimental works in the literaturereporting measurements of the roughness-induced normalmotion of macroscopic sliding solids. A notable exceptioncan however be found in [35, 36], where the authors mon-itor the normal acceleration of mild-steel slider-buttonsof centimetric radius of curvature, during sliding on arough mild steel disk. In particular, they report largewave number tails of the normal displacement PSDs ofthe type k − , in close agreement with our numerical find-ings (see Eq. 8).An interesting comparison can also be made with theliterature about stylus measurements of rough surfaces.All wavelengths of the topography that are smaller thanthe tip size will be filtered-out through a geometrical fil-tering process analogous to the one studied here, leadingto erroneous topography measurements (see e.g. [37]).Indeed, in [38], it is shown that while the amplitude oflarge wavelengths is accurately measured, that of smallwavelengths is underestimated. Thus, the rms value ofthe measurement is smaller than that of the topography,consistently with our results of Fig.7b. They also showthat the crossover wavelength separating both regimesscales as R / (2 − H ) , with R the curvature radius of aparabolic tip and H the Hurst exponent. This result in-dicates a size-dependence of the filtering process, whichis analogous to the L - (or N -) dependence that we ob-served. In [39], the authors further showed that geomet-rical filtering induces cusps in stylus measurement, andthat those cusps are responsible for a k − behaviour ofthe large wave number tail of the PSD. Again, this isfully consistent with our results (see Eq. 8). V. CONCLUSION
We addressed the question of the roughness-inducednormal motion during sliding of two solids in the limit ofvanishing normal load, i.e. when the contacting asper-ities do not deform. We considered the simplified caseof the quasi-static evolution of the separation at single-point contact, when the slider has no rotational degreeof freedom. Systematic numerical simulations assum-ing Gaussian self-affine topographies with various powerspectrum densities, and sliders with various sizes havebeen performed. We found that the normal motion re-lates to the topographies through a geometrical filteringprocess which depends on the size of the slider. We alsofound that the resulting normal motion (i) has enrichedspectral contents in the high wave number range, (ii)is non-Gaussian and (iii) has standard deviation muchsmaller than that of the sum topography. We providedempirical expressions relating the characteristics of thetopography to that of the roughness-induced normal mo-tion. We demonstrated that the distribution of the am-plitude of the normal motion can be well predicted withinthe framework of extreme value theory (EVT) as soon asthe number of points representing the topography of theslider is taken equal (to a prefactor close to 1) to thesurface of the slider divided by the square of the centralwavelength of the topography.These results are relevant whenever rough surfaces arebrought into light contact, that is, when there is no signif-icant deformation of the bodies. They can be useful not only for sliding surfaces, but also to assess the variabilityof static measurements made on statistically equivalentcontacts [13, 14]. In particular, such a variability is ex-pected to be much smaller than the characteristic am-plitude of the two antagonist topographies. Our resultsare limited to single-point contacts, when the two solidsare brought into contact through normal translation. Inthe case where a slider would be free to tilt, it would,under gravity, settle on three contact points to satisfyisostatic equilibrium. Accounting for such an effect is aninteresting topic for a future work.The fact that EVT nicely predicts the simulation re-sults indicate that computationally expensive simulationslike those decribed here may not be necessary in the fu-ture. Indeed, simple analytical formula (Eqs. 14- 15) orsemi-empirical expressions (Eqs. 16- 19) are sufficient toevaluate most of the relevant statistical descriptors of theroughness-induced normal motion. Our results thus fur-ther extend the already large range of applicability ofEVT to rough contact situations.
ACKNOWLEDGMENTS
This work was supported by LABEX MANUTECH-SISE (ANR-10-LABX-0075) of Universit´e de Lyon,within the program Investissements d’Avenir (ANR-11-IDEX-0007) operated by the French National ResearchAgency (ANR). It received funding from the People Pro-gram (Marie Curie Actions) of the European Union’s Sev-enth Framework Program (FP7/2007-2013) under Re-search Executive Agency Grant Agreement PCIG-GA-2011-303871. J.P.-L. is member of the Labex CeLyAof Universit´e de Lyon, operated by the French NationalResearch Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). K.T. acknowledges support from EarthFlows - Astrategic research initiative by The Faculty of Mathemat-ics and Natural Sciences at the University of Oslo. [1] H. Hamaker, Physica , 1058 (1937).[2] J. D. Jackson, Classical Electrodynamics (third edition)(John Wiley & Sons, Inc., 1999).[3] H. B. G. Casimir, Proc. Kon. Ned. Akad. Wet. , 793–795 (1948).[4] K. L. Johnson, Contact Mechanics (Cambridge Univer-sity Press, 1987).[5] B. N. J. Persson, Physical Review Letters , 125502 (2007).[6] A. I. Volokitin and B. N. J. Persson,Reviews of Modern Physics , 1291 (2007).[7] F. Plourabou´e, P. Kurowski, J.-M.Boffa, J.-P. Hulin, and S. Roux,Journal of Contaminant Hydrology , 295 (2000).[8] L. Talon, H. Auradou, and A. Hansen,Phys. Rev. E , 046108 (2010).[9] A. I. Vakis, V. A. Yastrebov, J. Scheibert, L. Nicola,D. Dini, C. Minfray, A. Almqvist, M. Paggi, S. Lee,G. Limbert, J. F. Molinari, G. Anciaux, R. Aghababaei,S. E. Restrepo, A. Papangelo, A. Cammarata, P. Nicol-ini, C. Putignano, G. Carbone, S. Stupkiewicz, J. Lengiewicz, G. Costagliola, F. Bosia, R. Guarino,N. M. Pugno, M. H. M¨user, and M. Ciavarella,Tribology International , 169 (2018).[10] M. H. M¨user, W. B. Dapp, R. Bugnicourt, P. Sainsot,N. Lesaffre, T. A. Lubrecht, B. N. J. Persson, K. Harris,A. Bennett, K. Schulze, S. Rohde, P. Ifju, W. G. Sawyer,T. Angelini, H. Ashtari Esfahani, M. Kadkhodaei, S. Ak-barzadeh, J.-J. Wu, G. Vorlaufer, A. Vernes, S. Solhjoo,A. I. Vakis, R. L. Jackson, Y. Xu, J. Streator, A. Ros-tami, D. Dini, S. Medina, G. Carbone, F. Bottiglione,L. Afferrante, J. Monti, L. Pastewka, M. O. Robbins,and J. A. Greenwood, Tribology Letters , 118 (2017).[11] R. Sahli, G. Pallares, C. Ducottet, I. E. Ben Ali,S. Al Akhrass, M. Guibert, and J. Scheibert,Proceedings of the National Academy of Sciences of the USA , 471 (2018).[12] V. A. Yastrebov, G. Anciaux, and J.-F. Molinari,Journal of the Mechanics and Physics of Solids , 469 (2017).[13] P. J. van Zwol, V. B. Svetovoy, and G. Palasantzas,Physical Review B , 235401 (2009). [14] W. Broer, G. Palasantzas, J. Knoester, and V. B. Sve-tovoy, Physical Review B , 155410 (2012).[15] C. Zouabi, Dynamique d’un contact glissant rugueux- rugueux sous faible charge: exp´eriences et mod´elisation,PhD thesis, Ecole Centrale de Lyon (2016).[16] J. Slaviˇc, M. D. Bryant, and M. Boltezar,Journal of Sound and Vibration , 732 (2007).[17] V. Dang, J. Perret-Liaudet, J. Scheibert, and A. Le Bot,Computational Mechanics , 1169 (2013).[18] C. Zouabi, J. Scheibert, and J. Perret-Liaudet,EPL , 50006 (2016).[19] D. Sornette, Critical Phenomena in Natural Sciences,Springer Series in Synergetics (Springer-Verlag,Berlin/Heidelberg, 2006).[20] P. Embrechts, C. Kl¨uppelberg, and T. Mikosch,Modelling Extremal Events (Springer Berlin Heidelberg,Berlin, Heidelberg, 1997).[21] L. de Haan and A. Ferreira, Extreme Value Theory,Springer Series in Operations Research and Financial En-gineering (Springer New York, New York, NY, 2006).[22] J.-Y. Fortin and M. Clusel,Journal of Physics A: Mathematical and Theoretical , 183001 (2015).[23] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi,Advances in Physics , 349 (2006).[24] M. S. Longuet-Higgins,Phil. Trans. R. Soc. Lond. A , 157 (1957).[25] S. Plaszczynski, Fluctuation and Noise Letters , R1 (2007).[26] G. Palasantzas, Physical Review B , 14472 (1993).[27] B. N. J. Persson, O. Albohr, U. Tartaglino,A. I. Volokitin, and E. Tosatti,Journal of Physics: Condensed Matter , R1 (2004).[28] E. H. Vanmarcke, Journal of the Engineering Mechanics Division , 425 (1972).[29] A. Preumont, Journal of Sound and Vibration , 15 (1985).[30] A. Preumont, Random Vibration and Spectral Analysis(Springer Science & Business Media, 1994).[31] A. Der Kiureghian, Journal of the Engineering Mechanics Division , 1195 (1980).[32] D. Benasciutti and R. Tovo,International Journal of Fatigue , 867 (2005).[33] M. R. Leadbetter, G. Lindgren, and H. Rootz´en,Extremes and Related Properties of Random Sequences and Processes,Springer Series in Statistics (Springer New York, NewYork, NY, 1983).[34] S. N. Majumdar and A. Pal, arXiv:1406.6768 (2014).[35] A. Soom and C. Kim,Journal of Lubrication Technology , 221 (1983).[36] A. Soom and C. Kim,Journal of Lubrication Technology , 514 (1983).[37] T. D. B. Jacobs, T. Junge, and L. Pastewka,Surface Topography: Metrology and Properties , 013001 (2017).[38] F. Lechenault, G. Pallares, M. George,C. Rountree, E. Bouchaud, and M. Ciccotti,Physical Review Letters , 025502 (2010).[39] E. L. Church and P. Z. Takacs, inProc. SPIE 1332, Optical Testing and Metrology III: Recent Advances in Industrial Optical Inspection(1991) pp. 504–515. Appendix A: Surfaces parameter
For surfaces described with radial power spectrum den-sities as shown in Fig. 2, the three first radial moments M , M , and M have the following expressions: M = 2 π S (cid:18) k l − k L − k αl ( k − αs − k − αl ) α − (cid:19) (A1) M = π S (cid:16) k l − k L − k αl ( k − αs − k − αl ) α − (cid:17) if α = 32 π S (cid:16) k l − k L + k l ( ln ( k s ) − ln ( k l )) (cid:17) if α = 3(A2) M = π S (cid:16) k l − k L − k αl ( k − αs − k − αl ) α − (cid:17) if α = 42 π S (cid:16) k l − k L + k l ( ln ( k s ) − ln ( k l )) (cid:17) if αα