In the schematic of Figure 1 we illustrate the concept of the proposed hybrid multiscale contact approach. We start from the typical roughness power spectral density (Figure 1C) as obtained by three-dimensional profilometry (B) of a real rough surface, such as a road asphalt surface (A). After splitting of the roughness spectral content in low (red in (C)) and high (blue in (C)) partitions, the contact mechanics problem is solved with the kinematically non-linear model on the high-slope scale range for a range of pressures and sliding velocities, and the results are averaged in terms of an effective contact interaction law (D) to be applied to the VHS-based linear model (E). Since the latter typically includes most of the roughness spectral content, the overall cost of the numerical computations is comparable to the cost of the VHS-based model computations.

More in detail, on the smallest length scales we simulate frictionless contact using the kinematically non-linear model, described in Section 2.4, for a number of applied average pressures, from which the average separation, the normalized contact area and the viscoelastic shear stress are computed. The effective contact interaction law giving the average pressure as a function of the average separation (Figure 1D) is then applied to the surface profile containing only the low-slope roughness components (E). A typical repulsive law, obtained for the adhesionless case, is shown in Figure 2 for a (two-dimensional) sinusoidal rigid roughness in sliding contact against a viscoelastic slab. The kinematically non-linear model results, reported for two different sliding speeds, show a monotonic increase of the pressure at decreasing average interface separations.

We consider now a rough rigid surface in steady sliding contact against a smooth rubber sample. The rigid surface is assumed periodically rough, with periodicity L ₀ in the x- and y-directions (coordinates in the average contact plane, with x = x , y ). The apparent projected interaction area is A ₀. The surface roughness z r x is decomposed, as described in Section 2.2, in the low-slope and high-slope scales, z ls x and z hs x , respectively, with z r x = z ls x + z hs x . First, the kinematically non-linear simulation (Section 2.4) is carried out with the surface roughness z hs x , at varying sliding speed v ₀ (directed along the x-direction), leading to the homogenized solution fields σ hs u ̄ , v 0 , α hs u ̄ , v 0 and Δ τ hs u ̄ , v 0 depending on the sliding speed v ₀ and on the average separation u ̄ . Here σ _hs is the average contact pressure σ hs u ̄ , v 0 = A 0 , h s − 1 ∫ A 0 , h s d 2 x p hs x , (1) where p hs x and A _0,hs are respectively the contact pressure field and the nominal in-plane area at the high-slope roughness scale. α _hs is the normalized projected contact area α hs u ̄ , v 0 = A 0 , h s − 1 ∫ A 0 , h s d 2 x H p hs x , (2) where H denotes the Heaviside step function, and Δτ _hs is the effective average shear contact stress resulting from the asymmetric deformation of the rubber sample induced by viscoelasticity Δ τ hs u ̄ , v 0 = A 0 , h s − 1 ∫ A 0 , h s d 2 x p hs x , u ̄ ∂ z hs ∂ x . (3)

The VHS-based frictionless simulation (Section 2.3) is now carried out on the surface z l s x , at a given externally applied average pressure σ ₀ and at sliding speed v ₀, by making use of the effective interaction law computed above, σ hs u ls x , v 0 , where u ls x is the local interface separation in the BEM contact geometry. The VHS-based model provides the contact pressure σ ls x and the interface separation u ls x . The total viscoelastic friction force F _T can thus be calculated as follows F T ≈ ∫ A 0 , l s d 2 x σ ls x ∂ z ls ∂ x + Δ τ hs  u ls x , v 0 , (4) where A _0,ls is the nominal in-plane area of the low-slope scale, leading to the total viscoelastic friction coefficient μ r = μ ls + μ hs , (5) with μ hs = 1 σ 0 A 0 , l s ∫ A 0 , l s d 2 x Δ τ hs u ls x , v 0 (6) μ ls = 1 σ 0 A 0 , l s ∫ A 0 , l s d 2 x σ ls x ∂ z ls ∂ x . (7) The total projected contact area, from which the adhesive contribution to friction can be calculated, is A c = ∫ A 0 , l s d 2 x H σ ls x α hs u ls  x . (8)

It is well established that most surfaces show roughness over a broad range of wavelengths (Persson, 2006), with few exceptions belonging to nature-optimized surfaces (showing few-scale hierarchical topography arrangement). A typical example are asphalt or concrete road surfaces, which are rough over many length scales and exhibit fractality (Klüppel and Heinrich, 2000; Persson, 2001). The characterization of roughness for contact mechanics applications involves the adoption of at least two-scales profilometry, e.g. with atomic force microscopy and white light interferometry to cover the full range of roughness wavelengths. Once the roughness data are available, the statistics are fully contained in the surface power spectral density (PSD) in case of Gaussian roughness (Persson, 2006). An equivalent alternative to the PSD is the height distance correlation (HDC) function R zz λ , introduced by Heinrich and Klüppel (Klüppel and Heinrich, 2000). Assuming a homogeneous roughness z r x , the HDC is R zz λ = z r x + λ − z r x 2 = 2 z rms 2 − R λ , (9) where 〈〉 denotes the ensemble average, and where the autocorrelation (AC) function R λ (for homogeneous roughness) is given by R λ = z r x + λ z r x . (10) The PSD C q is given by the Fourier transform of the AC C q = 2 π − 2 ∫ d 2 λ R λ e − i q ⋅ λ . (11) The simplest procedure to calculate the PSD is from the Fourier transform of the surface roughness z r q ; in particular we observe that z r q ′ z r q ′ ′ = 2 π − 4 ∬ d 2 x d 2 y z r x z r y e i q ′ ⋅ x e i q ′ ′ ⋅ y = C q ′ δ q ′ + q ′ ′ . (12) The Fourier transform C zz q of the HDC function can be easily linked to the Fourier transform C q of the AC function with C zz q = 2 z rms 2 δ q − C q . (13)

We now approximate the measured surface by a finite number of sinusoidal waves, following the approach suggested by Reinelt and Wriggers (Wriggers and Reinelt, 1996) and adopted by other researchers (Falk et al., 2016). For simplicity, we focus on a two-dimensional surface roughness, however, the approach can be equally applied to the more general case of three-dimensional roughness. In particular, the single sine wave function (characterized by amplitude Δ and generic wavelength 2π/q), has the following HDC function R zz λ = q 2 π ∫ 0 2 π / q d x Δ ⁡ sin q x + λ − Δ ⁡ sin q x 2 = 2 Δ 2 ⁡ sin 2 q λ 2 , (14) which we write in the following more convenient form R zz λ > 0 = 2 Δ 2 ⁡ sin 2 q λ 2 , λ ≤ π / q 2 Δ 2 , elsewhere . (15) We note that the slope in the double logarithmic diagram of R zz λ ≈ 0 is 2 for any wavelength, whereas road profiles usually have slopes in range of 2.1–2.5 (Schramm, 2003). Thus, as expected, multiple sine waves with different amplitudes and wavelengths must be superimposed in order to reproduce the power law behaviour of real rough surface HDCs. The root mean square slope is defined as m 2 = ∇ z r ⋅ ∇ z r , (16) which for the sinus specializes to m 2 q = q Δ / 2 .

Finally, the measured track is approximated by a sum over a discrete number M of wavelengths z r x = ∑ i = 1 M Δ i ⋅ sin q i x + ϕ i , (17) with a uniformly distributed random shift angle ϕ _i , and q _i = 2 ⁱ ⋅ q ₀, where the periodicity length is 2π/q ₀.

The low-slope scale contact model is based on the adoption of the VHS assumption. The interaction occurs between a rigid rough surface (z _ls) and a smooth VHS, in relative steady sliding. This contact model has been first reported in (Persson and Scaraggi, 2014), and is summarized in the following with the addition of the interface contact law contribution coming from the high-slope scale roughness.

In Figure 3 we show the schematic of the low-slope scale contact interface. The local interface separation, corresponding to the local distance between the mating surfaces, measured in a reference moving with the rigid rough substrate, is u ls x = u ̄ + w x − z ls x , (18) where u ̄ is the average interface separation, w x the surface viscoelastic normal displacement, and z ls x is the low-slope surface roughness, with w x = z ls x = 0 . We define the following Fourier transforms (q is the transformed variable of x) w q = 2 π − 2 ∫ d 2 x w x e − i q ⋅ x (19) and σ ls q = 2 π − 2 ∫ d 2 x σ ls x e − i q ⋅ x , (20) where σ ls x is the distribution of contact pressures ( σ 0 = σ ls x is the average contact pressure). Following the discussion reported in (Persson and Scaraggi, 2014), w x can be related to σ ls x through a simple equation in the Fourier space w q = M z z q σ ls q , (21) where M z z q is the surface response of the VHS in the frequency domain M z z q , ω = 2 q E r ω . (22) Here, E r ω is the frequency-dependent complex reduced Young’s modulus, given by E r ω = E ω 1 − ν 2 , (23) with ν as the Poisson’s ratio. The Poisson’s ratio undergoes only a minor variation with frequency and thus can be assumed to be constant (Persson, 2001). The linear viscoelastic modulus E ω can be measured using standard techniques, and its real and imaginary part are typically fitted by a generalized Maxwell model (Lorenz et al., 2014; Scaraggi and Persson, 2015). The relation between separation u _ls(x) and interaction pressure σ _ls(x) within the Derjaguin’s approximation (Derjaguin, 1934) can be written in terms of a generic interaction law (Persson and Scaraggi, 2014) σ = f u . (24) Derjaguin’s approach allows to describe the force acting between close elementary portions of interacting surfaces in terms of the force per unit area acting between two planar semi-infinite walls. This approximation is typically adopted when i) the solids are separated by a distance that is much smaller than the local radius of curvature of the surfaces, and ii) the interaction decays sufficiently rapidly with distance, so that interactions between far enough elementary surface portions are negligible.

In this work we integrate the repulsive term of the Lennard-Jones potential to obtain the traction due to the adhesionless interaction in (24), whereas the repulsive contribution coming from the high-slope scale is taken into account by adding to the integrated Lennard-Jones the high-slope scale traction σ hs u ̄ , calculated as described in the following.

The kinematically non-linear model considers a deformable body undergoing adhesionless and frictionless contact with a superposition of rigid sinusoidal surfaces. The generic sinusoidal surface is supposed to represent one of the high-slope wavelengths belonging to the roughness spectrum of the surface. In order to compute the effective interface properties, a micromechanical numerical test is conducted, see Figure 4 (for simplicity we show only one wavelength). Because of the periodic geometry of the surface, the width of the rubber slab is taken equal to the largest wavelength of the sinusoidal spectrum; the specimen is considered as a periodic cell extracted from a long boundary layer. The specimen must have sufficient height in order to incorporate the mostly stressed material region, so that the whole amount of energy dissipation taking place at the micro-level is accounted for. In (De Lorenzis and Wriggers, 2013), a ratio H/L ₀ = 0.75 was found sufficient for all analyzed dragging velocities. Thus, this ratio is also adopted in this work. From the numerical standpoint, a sufficiently fine mesh must be used in order to achieve a satisfactory resolution of the contact area and hence an accurate estimation of the contact stresses. Again exploiting results in (De Lorenzis and Wriggers, 2013), we adopt a discretization with 120 × 90 elements. At the beginning, the specimen is pressed onto the sinusoidal surface with a predetermined normal pressure of absolute value p ̄ , which stays constant during the sliding process. Then, the upper side of the specimen is moved horizontally with a specified constant velocity v ̄ . Since the specimen is considered as a periodic cell, periodic boundary conditions are applied on its lateral surfaces. The finite linear model for viscoelasticity proposed by Reese and Govindjee (1998) is used. The model is able to capture the large deformations arising in the contact zone of the rubber sliding over a rough surface. The model is still “linear” with respect to the evolution equation (Lion, 1997a; Lion, 1997b), i.e. it assumes small deviations from thermodynamic equilibrium. Full details can be retrieved in (Lion, 1997a) or (Reese and Govindjee, 1998), see also (De Lorenzis and Wriggers, 2013). In order to extract a single value for the macroscopic friction coefficient from each micromechanical test, time averaging is applied (Wriggers and Reinelt, 1996; De Lorenzis and Wriggers, 2013; Wagner et al., 2015).

A confocal optical profilometer is used in this work to measure the surface roughness of asphalt specimens. In order to obtain representative data, 30 line scans in different regions of the specimen are performed using the optical profilometer in different directions. Each line scan has a measurement length of 16 mm that is consistent with the typical grain size of asphalt. The lateral resolution of the profilometer (0.05 mm) gives an acceptable degree of accuracy (within the purpose of this work). The lower cut-off is selected as 0.0625 mm. The upper cut-off is selected as the typical largest grain size of asphalt, and ranges between 5 and 20 mm. The asphalt mixture is AC-8, the measurement speed is 0.1 mm/s.

The superposed sinusoidal functions are fitted to the averaged measured surface based on their HDC functions, using the least square algorithm on a double logarithmic scale, which identifies the amplitudes Δ _i and the spatial frequencies q _i such that ∑ j = 1 N ∑ i = 1 M log 2 Δ i 2 ⁡ sin 2 q i λ j 2 − log R zz λ j 2 1 2 → min λ j ≤ π q i (29) ∑ j = 1 N ∑ i = 1 M log 2 Δ i 2 − log R zz λ j 2 1 2 → min λ j > π q i , (30) where N is the number of discrete sample points where the profile height is measured, and M is the number of superposed sine waves (each with its amplitude and wavelength) in the approximation. The approximation is carried out with different values of M. The slopes of the measured HDC and of the HDC approximated by more than eight sine waves and their cut-off-lengths were found to be identical. The HDC functions of the sine waves and of their superposition are depicted in Figure 6, where the red curve denotes the averaged measured asphalt road track. The agreement is very good. Table 2 lists the wavelengths and amplitudes of the 8 waves along with the corresponding root mean square slopes. It is important to note that only one scale has a root mean square slope larger than 1.

As follows, we illustrate the results obtained for the smallest scale of the approximated rough surface. These results are obtained applying the kinematically non-linear model of Section 2.4 (and by comparison the VHS-based model of Section 2.3 to a contact problem with a sinusoidal surface with a wavelength of 0.0625 mm and an amplitude of 0.01867 mm. The root mean square slope roughness is 1.3275, see Table 2.

Figure 7A shows the obtained rolling friction coefficient as a function of the applied pressure for a sliding velocity of 10 mm/s belonging to the rubbery region. For small applied pressures, the VHS-based model only slightly underestimates the friction coefficient compared to the kinematically non-linear model. This can be justified by the fact that at very small applied pressures, the contact occurs only at the top of the asperity, so the geometric non-linearity can be neglected. A decreasing friction coefficient for high pressure values is observed for both approaches. One can note that for applied pressures larger than 2 MPa, the difference between the approaches becomes larger, as the deformation regime induced by large root mean square slope roughness significantly deviates from the small displacement assumption and thus the effect of the kinematic non-linearity comes into play. Figure 7B gives the same curve for a sliding velocity of 50 mm/s belonging to the maximum dissipation region. Similar trends as in Figure 7A are observed, with a higher maximum friction value. The deviation between the kinematically non-linear and the linear VHS-based approaches is very large in this case. These results are in agreement with the observations in (Scaraggi et al., 2016).

The friction coefficient at small-slope scales (wavelength 2–8) from Table 2 is now calculated using the repulsive law deduced from the high-slope scale, which is shown in Figure 2. Figure 8 illustrates examples of the contact pressure field as a function of the dimensionless contact position for one realization of the rough surface, with dragging velocities of 10 mm/s and 50 mm/s, respectively. It should be noted that according to these contact pressures, the corresponding friction coefficients at the contact points are approximated using the results of the friction coefficients at the smallest scale, which are shown in Figures 7A,B.

Figure 9 illustrates the friction coefficient predicted using the VHS-based model and the hybrid multiscale approach, for two different velocities belonging to the rubbery and the maximum dissipation regions. The discrepancy in the friction coefficient is consistent with the results in Figure 7. For a low dragging velocity of 10 mm/s, the difference between the VHS-based prediction and the results of the hybrid multiscale approach is relatively small. This can be justified by the fact that for low sliding velocities, the local contact pressures are more evenly distributed and they are small, see Figure 8A. Moreover, for small pressures, the results of the kinematically non-linear model and of the VHS-based model are comparable, see Figure 7A. However, for a dragging velocity of 50 mm/s, the difference is large, which can be ascribed to the irregular distribution of the contact pressures at larger dragging velocities and to the resulting large contact pressures, see Figure 8B. The effect of the non-linearity is remarkable, see Figure 7B. From comparison of the results in Figure 5 to those in Figure 9, one can note that in the VHS the smallest wavelengths of the roughness, ranging from 0.0625 to 0.25 mm, give the largest contribution to the friction coefficient, which seems unrealistic.