Asperity-based modification on theory of contact mechanics and rubber friction for self-affine fractal surfaces

Modeling the real contact area plays a key role in every tribological process, such as friction, adhesion, and wear. Contact between two solids does not necessarily occur everywhere within the apparent contact area. Considering the multiscale nature of roughness, Persson proposed a theory of contact mechanics for a soft and smooth solid in contact with a rigid rough surface. In this theory, he assumed that the vertical displacement on the soft surface could be approximated by the height profile of the substrate surface. Although this assumption gives an accurate pressure distribution at the interface for complete contact, when no gap exists between two surfaces, it results in an overestimation of elastic energy stored in the material for partial contact, which typically occurs in many practical applications. This issue was later addressed by Persson by including a correction factor obtained from the comparison of the theoretical results with molecular dynamics simulation. This paper proposes a different approach to correct the overestimation of vertical displacement in Persson’s contact theory for rough surfaces with self-affine fractal properties. The results are compared with the correction factor proposed by Persson. The main advantage of the proposed method is that it uses physical parameters such as the surface roughness characteristics, material properties, sliding velocity, and normal load to correct the model. This method is also implemented in the theory of rubber friction. The results of the corrected friction model are compared with experiments. The results confirm that the modified model predicts the friction coefficient as a function of sliding velocity more accurately than the original model.


Introduction
Despite numerous efforts in the field of contact mechanics, modeling of the real area of contact between two solids squeezing together is an ongoing challenge. A summary of recent attempts to solve a contact mechanics problem can be found in Ref. [1]. Accurate modeling of the real contact area is fundamental to understanding and modeling of many physical processes such as friction, adhesion, wear, contact resistivity, frictional heating, etc. The study on contact mechanics was pioneered by Hertz [2], who modeled the contact between two elastic spheres with smooth surfaces. Following his approach, several classical contact theories were developed to model the contact between a rough surface and a flat rigid plane. One of the well-known models is Greenwood and Williamson (GW) contact theory [3], which models the rough surface as a plane covered with spherical-shape bumps with the same radius of curvature and Gaussian height distribution. The GW model considers a single length scale roughness and neglects the interactions between the asperities. Therefore, this contact model, in its original form, is unable to model the real contact area accurately. Since the original publication, several revisions to this model have been published in the  Refs. [4][5][6][7][8]. Alternatively, Bush et al. [9] developed a more detailed contact model, known as the BGT model, based on the Hertz contact theory. They modeled multiscale roughness using elliptical paraboloid asperities with random aspect ratio and orientation. However, they also neglect the interaction between asperities, which results in similar errors as the GW model. In reality, the contact areas are neither circular nor elliptical. Greenwood [5] argued that using the geometric mean of the summit curvatures in GW theory, the results of the BGT model can be reproduced more conveniently, and perhaps more accurately. In addition, the effect of interactions between asperities was later considered (approximately) in an improved version of the GW theory [6], and this theory remains popular in the contact mechanics community. A review of other contact theories based on similar concepts can be found in Ref. [10].
As a pioneer in the hierarchical modeling of the contact area, Archard [11] developed a contact model, in which surface roughness was described as hierarchical uniform spherical asperities on top of larger asperities. Using Hertz contact theory in his model, he showed that by increasing the number of length scales in the model, a more realistic representation of the real contact area versus normal load would be achieved. Although Archard's model can represent the multiscale nature of roughness and its contribution to the real area of contact, it is not useful in practical applications, as it idealizes the roughness as uniformly distributed spherical asperities. Several decades later, Persson [12] used the idea of multiscale modeling of rough surfaces to introduce a completely different approach to model the real contact area. He considered the probability distribution of normal stress (contact pressure) at different length scales to develop a theory of contact mechanics. He implemented the surface roughness in his theory using its power spectrum, calculated from the surface height profile. He assumed that the smooth, soft solid could deform and penetrate in large substrate valleys; therefore, the displacement on its contact area could be approximated by the height profile of the rough rigid surface. This assumption leads to notable errors in the modeling of the real contact area under small normal loads and partial contacts [13][14][15]. Therefore, he later proposed a correction factor to match the theoretical results to the molecular dynamics (MD) simulations [16]. However, the physical interpretation of this factor is not clear as it is referred to as a fudge parameter by Dapp et al. [17], who systematically analyzed Persson's contact mechanics theory.
In this paper, an asperity-based approach is proposed to correct the displacement assumption for contact mechanics on self-affine fractal surfaces using an affine transformation. A self-affine fractal surface has the property that if it is scaled by different factors along parallel and perpendicular directions, the scaled surface has the same topography as the unscaled one [18]. Many surfaces in nature have the properties of a self-affine fractal topography, and fractal descriptions can be used to describe their roughness [19][20][21]. Several studies [22][23][24][25][26][27] have shown that sandpaper and asphalt have approximately self-affine fractal surfaces over several orders of magnitude in spatial scale. Therefore, this approach can be applied in the modeling of contact and friction of a tire tread compound on sandpaper and asphalt surfaces. Two different methods are considered to find the affine transformation parameter. The first method is based on the local peak distribution of the surface profile for nonsymmetric height distribution, such as abraded sandpaper. The second method is based on an improved Greenwood-Williamson (IGW) contact theory for normal height distribution. This paper is organized as follows: Section 2 presents a brief review of Persson's multiscale contact theory and its main assumptions. In Section 3, the asperitybased approach to correct the displacement assumption in the Persson theory is presented, and two methods to find the correction factor are explained. Section 4 summarizes the rubber friction theory modified using the correction method introduced in this paper. Then, Section 5 presents the input parameters of the theoretical models. In Section 6, the numerical results of the modified theory are presented and discussed, and the asperity-based factors are compared with Persson's correction factor. In addition, the results of the original and modified friction models are compared with some experimental results. Finally, Section 7 presents the summary and conclusions of this work.

Summary of Persson's multiscale contact theory
Attempting to model rubber friction on a rough rigid surface, Persson [12] developed a multiscale contact theory for an isotropic, linear viscoelastic (and elastic), semi-infinite half-body based on a novel approach. He used the surface roughness power spectrum to model the contact area of a multiscale rough surface with a small-slope height profile. He started from the assumption of a complete contact between two solids, and then he introduced the partial contact by imposing a detachment boundary condition. In the original theory, the detachment occurs when the normal pressure becomes zero, assuming the adhesive force between two surfaces is negligible. Later, he also considered the effect of adhesion in the detachment boundary condition [28,29] with application to smooth contact surfaces and very soft materials such as gelatin.
Persson introduced a pressure probability distribution on the contact area as a function of magnification (length scale), at which the contact area is studied. At the lowest magnification (the largest length scale), it appears that there is complete contact between two surfaces, and the apparent contact area is the same as the nominal contact area. Similarly, the pressure distribution at the lowest magnification is equal to the nominal pressure distribution. However, partial contact between two surfaces can be observed by studying the contact area at higher magnifications, e.g., under an optical microscope. Complete contact regions can be detected on the top of the asperities, while surface separations and gaps exist in the valleys between asperities. Therefore, the apparent contact area under an optical microscope, which is the sum of all the contact patches, is smaller than the nominal contact area. If each contact patch at the macroscale is studied further at much higher magnifications, e.g., under an electronic microscope, smaller asperities, and accordingly smaller contact patches can be detected on each macro-asperity as shown in Fig. 1. Therefore, the apparent contact area decreases and approaches the real contact area as the magnification increases.
If A0 and σ0 denote the nominal contact area and the pressure measured at the lowest magnification, respectively, then the apparent contact area and the average pressure at magnification  must satisfy the following equation: where     and  A denote the average pressure and the apparent contact area at magnification  .
Using the stress probability distribution, the average contact pressure as a function of magnification can be calculated using: where   ( , ) P represents the pressure probability distribution. Therefore, the ratio of the apparent contact area at magnification  to the nominal contact area can be calculated using: As the apparent contact area decreases with the increase of magnification, the contact pressure distribution becomes broader [30] as shown in Fig. 1, and it can be described by a diffusion type equation as: where the symbol <…> stands for ensemble average, where  δ( ) stands for the Dirac delta function. In addition to these boundary conditions, the diffusion coefficient which is a function of multiscale deformation on the contact surface. The stress-displacement relation can be found using the constitutive equation, and the displacement field must satisfy the linear momentum balance law in terms of displacement (Navier equation). In this theory, it is assumed that the vertical displacements can be approximated by the height profile of the rough surface while the lateral displacements are neglected. Then, the mean square of the height profile is calculated from the surface roughness power spectrum, therefore: where z u and h are displacement and surface height profile, respectively, and ( ) C q is the surface roughness power spectrum. They are calculated as a function of wavevector   ( , ) x y q q q of the surface roughness. The surface roughness power spectrum is defined as: . The details of the derivation of the diffusion coefficient as well as the solution for Eq. (4) can be found in the original paper by Persson [12]. The final results provide the ratio of the apparent contact area to the nominal contact area, which for a linear viscoelastic material can be calculated using the following equation: where E is the viscoelastic modulus, which is a function of the frequency of loading applied by the multiscale asperities along the sliding direction  ( cos ) vq , and  is the Poisson ratio whose dependency on the frequency of loading can be neglected. The parameter 0 q is the short cutoff wavenumber corresponding to the largest wavelength that contributes to the contact area.
The main shortcoming of Persson's multiscale contact theory is that the vertical displacement field is assumed for complete contact between two solids while deriving the diffusion coefficient. However, this assumption is only valid when there is complete contact with no gap anywhere between two contact surfaces, which is generally not true and has been criticized in Refs. [17,31]. This overestimation of displacement results in a larger diffusion coefficient and accordingly broadening of pressure distribution. When the pressure distribution becomes broader, the ratio of the apparent contact area to the nominal contact area decreases based on Eq. (3). Consequently, the theory underestimates the real contact area as it overestimates the elastic energy due to large displacements. This error was later addressed by Persson as he introduced a correction factor based on the comparison of his theoretical results for stored elastic energy with calculations obtained from MD simulation [16]. This factor can be implemented in the original theory using the following equation: where   0.5 [32] and This correction factor can be multiplied by Eq. (7) to correct the relationship between the elastic energy and vertical displacement for partial contact. However, the physical interpretation of this factor is not clear. As Dapp et al. [17] systematically analyzed Persson's contact mechanics theory, they referred to it as a fudge factor that is parameterized to match the numerical results. In the next section, an asperity-based approach is introduced to obtain a correction factor based on roughness parameters of a self-affine surface and the effective elastic modulus.

Asperity-based correction factor for Persson's contact theory
Several studies [22][23][24][25][26][27]33] showed that the concept of self-affinity could be applied to many rough surfaces of interest, such as asphalt and sandpapers. For selfaffine fractal surfaces, an affined transformed height profile can be introduced to substitute the original height profile. This transformed height profile can represent a more accurate fraction of the surface that contributes to the contact deformation. When an affine transformation with a fixed maximum height of max z is used for a self-affine fractal surface, the fractal dimension and the shape of its height distribution do not change. In this transformation, the original height profile z = h(x) is shifted towards the upper regions of the macro-asperities that have apparent full contact with the soft smooth solid at the macroscopic scale. Using this approach, the large valleys between those asperities that have no contact with the deformed material and do not contribute to the displacement field can be eliminated. Consequently, the errors due to the large displacement assumption corresponding to these valleys could vanish. Considering  max z max{ ( )} h x as the upper fixed boundary, the surface profile is shifted towards the upper regions using the following affine transformation: where T p is an affine shift parameter and ( ) T h x is the transformed height profile. After this transformation, the mean height profile shifts from zero to    and this becomes the new nominal contact plane and coordinate system. This transformation from full contact to partial contact for a sinusoidal roughness is shown schematically in Fig. 2.
After eliminating the offset of the transformed profile by zeroing the new mean   T h , the mean square of the transformed height profile can be calculated as: Accordingly, the corrected mean square of displacement in the vertical direction is equal to the mean square of the transformed height profile of the rough surface as: The surface roughness power spectrum can be obtained as an analytical function for a self-affine fractal surface and written as: where H is Hurst exponent, 0 h is the root mean square (RMS), and 0 q is the smallest wavenumber of the surface that contributes to the contact mechanics.
Two different simple asperity-based methods are proposed to determine the shift factor T p . The first method is based on surface characterization and distribution of the local maximum height of the surface profile. Using the first method, the errors due to the asymmetry of the surface profile can be reduced to some extent. However, in this method, the effect of material properties is neglected in determining the fraction of the height profile involved in the surface deformation. On the other hand, the second method is based on macroscale deformation and the separation distance between two surfaces using an improved version of the classical asperity model of GW contact theory.

Local peak distribution method
The simplest approach to find the shift factor is based on the analysis of the surface roughness profile. In this method, it is assumed that the fraction of the substrate profile that is in contact with the deformable surface can be found using the local peak (maxima) distribution of height profile. Thus, the local peaks within the range of macroscale wavelengths must be found from numerical analysis of the surface profile as shown in Fig. 3. The average of local peak distribution, as shown by the red dashed line in Fig. 3, can be used to shift the original height distribution to the upper regions of the rough surface using the affine transformation   | https://mc03.manuscriptcentral.com/friction presented in Eq. (9). It is assumed that the average of the transformed profile with respect to the original midplane is equal to the average of local peak distribution   Peaks h therefore the shift parameter becomes: The shortcoming of this method is that the shift factor is independent of nominal pressure and sliding velocity, and only considers the surface roughness characteristics to provide a better approximation of displacement field. However, when the surface profile is asymmetric and has remarkable skewness, a combination of this approach and the macro-scale deformation analysis explained in the following section can be used to obtain more accurate results.

GW theory and the improved version
The original GW [3] theory is based on contact between a rough surface and a flat plane. The rough surface is modeled as a plane covered with spherical shape asperities (bumps) with the same radius of curvature and a Gaussian height distribution. The original GW model considers a single-scale roughness and neglects the effect of elastic interaction between the asperities. Under a small normal load and at large length scales when the contact patches are sufficiently separated, it might be reasonable to assume that the deformation fields induced by the macro-asperities are independent. Therefore, the penetration depth of each rigid asperity into the elastic surface can be calculated using the Hertz contact theory [2]. The penetration depth, as shown in Fig. 4, for a single asperity indenting a semi-infinite flat surface is: where N F is the nominal normal load of a single asperity, R is the radius of curvature of asperity, and * E is the effective modulus calculated as: In the case of rubber in contact with sandpaper or asphalt, the elastic modulus of rubber is about 2 to 3 orders of magnitude less than the substrate surface, and it is reasonable to assume that the asperity is rigid and    In modeling multiple macro-asperities, it is assumed that the normal load distributes over all the macro contact regions, and it can be calculated as the sum of the normal forces applied over all the asperities on the nominal contact area. Assuming ( ) z is the height distribution of the surface profile and defining    d z as the separation distance between the rubber surface and the centerline of the rough substrate (  0 z ), the rubber is in contact with macro asperity if  z d . Then, based on the original GW theory, the nominal pressure is equal to: where a n is the number of asperities per unit area. In Eq. (16), the magnitude of the dynamic modulus at the minimum loading frequency (   min 0 q v ), which is applied by the macro-asperities with wavenumber 0 q , must be used. If the surface height profile has a normal distribution, then where   is the standard deviation. Typically, the asperities form a disordered hexagonal-like distribution [34] with a lattice constant   0 0 2π q , then: Assuming the largest length scale (macro-scale) roughness can be described using a sinusoidal function as h is the RMS of the surface height profile. Then, the average radius of macro-asperities can be calculated as: Using Eqs. (17) and (18) in Eq. (16) and knowing the nominal pressure, the only unknown variable is the separation distance "d" that has to be found using a numerical method. However, the original GW model might underestimate the separation distance since it does not consider the effect of interaction between the asperities on the penetration distance. When the normal load on the nominal contact area is not sufficiently small, the interaction between the asperities cannot be neglected. Ciavarella et al. [7] also showed that there are significant differences between the results of the original GW theory and the numerical results of a contact simulation between a rough surface and a flat plane at intermediate loads. Therefore, Ciavarella et al. [6] proposed an improved version of the GW theory to include the interaction between the asperities. In this improved version, it is assumed that the asperities are uniformly distributed over the contact area and hence, the deformation of the contact area is uniform and can be approximated as  * 0 a / A E over a compact area a A according to Timoshenko and Goodier [35]. Therefore, the interaction between asperities results in an increase of the effective separation distance between the mean planes of two surfaces, so that the effective separation distance becomes: where  0 2π a A q and using Eqs. (17)- (19) in Eq. (16) gives: (20) which is an implicit function of nominal pressure (  0 ) and must be solved iteratively.
The effective separation distance at the largest length scale represents the fraction of the rigid height profile, which does not contribute to the vertical displacement of the deformable surface. Shifting the height profile to the higher regions so that the minimum height of the transformed profile equals the effective separation distance eliminates the valleys where no contact occurs. The shift factor is then calculated using the following equation: Using this method, the shift factor is not only a function of roughness parameters but also a function of nominal pressure and elastic properties of the materials. When the height distribution of the rough surface is non-gaussian, e.g., worn or polished surfaces, it is suggested that instead of the original height distribution, the local peak distribution as discussed in the previous section is used in Eq. (20).
Since in rubber materials, the viscoelastic modulus highly depends on sliding velocity, the shift factor is also a function of velocity. At high sliding velocities, the rubber exhibits glassy behavior due to the high frequency of oscillating loads applied by substrate asperities, and its stiffness increases. Therefore, it cannot deform sufficiently to fill the gap between the asperities and the effective separation distance becomes larger. On the other hand, at low sliding velocity, the rubber undergoes large deformation, and the effective separation distance decreases, consequently; it is expected that the shift factor approaches unity at very Since the IGW method is based on the Hertz contact theory with the assumption of small deformation, it gives more accurate results for the rubber in the glassy state rather than the rubbery state. Therefore, the results are more reliable at high velocities when the rubber has higher stiffness. Hongyan et al. [36] also studied the applicability of the Hertz contact theory to large deformation using finite element methods. They found that the Hertz contact theory can predict the maximum contact stress under large deformation with about 2% error and overestimates the contact length with about 33%, but they did not validate the results with experiments.

Theory of rubber friction
Persson used his multiscale contact mechanics theory to develop a theory of hysteresis friction in viscoelastic materials, assuming that the energy dissipation due to friction is equal to the energy loss due to internal friction (hysteresis effect) in the rubber. Based on this theory, the hysteresis friction coefficient can be calculated using the following equation: (22) where max ζ is the magnification corresponding to the maximum wavenumber (upper cutoff) that contributes to the hysteresis friction. It is proposed by Lorenz et al. [37] that the upper cutoff wavenumber is determined by the condition that including all the roughness components up to the cutoff results in the (cumulative) root mean square slope of 1.3 or reciprocal of the road contamination particle size (whichever is smaller). Note that the correction factor ( ) S q did not exist in the original model, and was added later to this model to match the MD numerical results, resulting in the reduction of hysteresis friction. Although adding this correction factor is necessary to improve the model to match the numerical simulations [38], but as Dapp et al. [17] stated: "it has so far been implemented only heuristically".
Persson developed a friction theory assuming the dominant friction mechanism at high velocities is the hysteresis friction; however, many experimental results [39][40][41][42][43][44][45] suggest that the hysteresis friction could not accurately explain the rubber friction over a wide range of velocities and surface conditions. Rubber friction is generally due to more complex processes occurring on the contact area. In addition to hysteresis loss due to oscillating loadings applied by multiscale asperities, contact processes include adhesion, interfacial crack propagation, shearing on possible lubricant and contamination particles at the interface, and different wear processes [44,45]. Therefore, Lorenz et al. [44] proposed a semi-empirical model of friction coefficient for low sliding velocities. This semi-empirical model includes the contribution from the area of real contact to the friction, while the real contact area is calculated from Persson's contact theory. This semi-empirical model is formulated as where  This friction model, without the correction parameter ( ) S q , is used along with the asperity-based correction factor proposed in this paper to calculate the friction coefficient of a rubber sample sliding on an asphalt track at different velocities. In order to evaluate the modified model, the theoretical results are compared with some experimental results obtained using a sliding friction test set-up designed and developed in the Center for Tire Research (CenTiRe) [46]. The input parameters for theoretical models are explained in the next section.

Rubber dynamic modulus
One of the input parameters of the theoretical model is the viscoelastic master curve of the rubber sample for a wide range of frequencies. Dynamic mechanical analysis (DMA) is a standard method for the measurement of viscoelastic modulus, which uses an oscillating load to deform the sample with constant strain amplitude at a limited range of frequencies.
Since there are always some limitations to the testing frequency, the tests must be repeated at different temperatures. Temperature-frequency superposition is typically used to obtain the master curve for a broader range of frequencies. Since the variation of the viscoelastic modulus is very high near the glass transition temperature, in this work, a smaller temperature step was used around this temperature so that the viscoelastic modulus curves could overlap one other. The overlapping of curve segments is necessary to obtain the shift factor. The amplitude of strain was kept constant at 0.5% to avoid the effect of softening in the filled rubber and to remain in the linear region. The horizontal shift factor was obtained from shifting the curves with continuous overlaps to have a smooth master curve. The dynamic modulus of two Styrene-Butadiene Rubber (SBR) compounds used in this study is shown in Fig. 5. Compound A has a glass temperature of Tg= -32 °C and compound B has a glass temperature of Tg= -21.7 °C.
When a rubber block slides on a rough surface, its surface undergoes large strains of the order of 100% or even more as observed by Iwai and Uchiyama [47]. Therefore, the material properties at large strain provide a more accurate representation of the material properties close to the contact surface. To obtain a viscoelastic master curve for large deformation, the rubber samples were tested at a fixed frequency of 1 Hz and the strain amplitude is changed from a small value (about 0.5%) to the maximum limit of the sample. This procedure was repeated for different temperatures ranging from -25 to 80 °C. Then, using the frequency-temperature shift factor aT, the master curve for large strain values was approximated over a wide range of frequencies and used in the theoretical model.

Surface roughness power spectrum
An optical profilometer (Nanovea JR25) was used to measure the surface height profile. This profilometer offers high accuracy and precision for both lateral and height positions with submicron resolution. The profile measurement is considered valid only if the drop-out rate for the height signal is less than 10% [48]. Linear interpolation was used to replace the invalid samples. Before the analysis of height distribution, any curvature, slope, or offset that might exist in the measured profile must be eliminated. A sample of line-scan of asphalt profile before and after slope and offset suppression is shown in Fig. 6.
Before the calculation of the power spectrum of the surface profile, it is necessary to eliminate the | https://mc03.manuscriptcentral.com/friction spectral leakage of the signal using an appropriate windowing technique. The spectral leakage occurs when the discrete Fourier transform (DFT) is used for surface profile analysis. In DFT analysis, it is assumed that the input signal (surface profile) is repeated periodically with a period equal to the length of a measured profile in horizontal (x or y) direction. Therefore, at the edges of the height signal, a jump in the composite signal might occur. In order to prevent spectral leakage, a window must be applied to reduce the signal to zero at the edges. There are various windowing techniques available for this purpose, and two of the common windows widely used for calculation of the power spectrum of the surface are Hann and Split Cosine Bell (also known as Cosine Digital Tapering) windows. The drawback of using a Hann window is that it reduces the effective length of the signal, which influences the low-frequency content of the height signal. On the other hand, the Split Cosine Bell Window (SCBW) keeps 80% of the original length of the signal unaffected because it goes from zero to unity and vice versa much faster than the Hann window. Therefore, SCBW is very suitable for short length signals [33], and was used in this study, which is defined mathematically as:  (26) where is the discretized wavevector. The fast Fourier transform (FFT) algorithm was used to calculate    h q , and then the surface roughness power spectrum was calculated using the following equation [30]: The calculated power spectra of the 120-grit sandpaper and asphalt surface used in numerical calculations are shown in Figs. 7(a) and 7(b), respectively.
The short cutoff wavenumbers of both surfaces are shown in Fig. 7, both are about 2,000 m -1 . The Hurst exponent of the fractal surfaces can be calculated using the slope of the power spectrum as: The Hurst exponent of the sandpaper is 0.7, and for the asphalt surface, it is 0.87. A smaller value of Hurst exponent indicates that more irregularities exist in the topography of the surface. The RMS values of sandpaper and asphalt are 76.572 and 265.7 μm, respectively. The maximum magnification for the large cutoff wavenumber is chosen as 500 for both surfaces.

Modified contact area
The asperity-based approach for modification of the contact theory was initially evaluated by comparing the results with the factor that Persson introduced in his model to fit the results with MD numerical simulations. In the asperity-based approach introduced in this paper, the surface roughness power spectrum has been transformed in the model, which affects the real contact area as presented in Eq. (7). On the other hand, in Persson's approach, the real contact area remains the same but a correction parameter is multiplied by the original contact model. Therefore, to compare the results of these two different approaches, the correction factor (CF) of the asperity-based approach is expressed as m

A q C q
A q C q is also approach similar value as ( ) S q but to provide a better demonstration of the com- Next, the modified multiscale contact area is compared with the original theory to demonstrate how it affects the calculation of the real contact area. As mentioned previously, the validity of the Persson original contact model for partial contact is questionable. It is expected that the original model underestimates the real contact area [49] due to the overestimation of displacement and elastic energy. However, using the methods proposed in this paper, some of the errors due to this issue are eliminated.

Modification of contact area based on local peak distribution
The local peak distribution method was used to find the shift factor for 120-grit sandpaper. The material properties of compound A were used in Eq. (7). When the local peaks were selected within a fixed minimum distance of   0 0 2π q , few asperities affecting the contact area were not captured using the peak finding algorithm and resulted in an overestimation of the shift factor,  . Figure 8 shows the local peaks used in the calculation of the shift factor and the comparison between the local peak correction factor (LP-CF) and the Persson fitted correction factor (PF-CF).
On the other hand, when a more flexible range of distance between the local peaks, e.g., within  0 to  0 2 was considered to find the effective asperities, the shift factor reduced to  1 2 T p , which results in a perfect match between LP-CF and PF-CF as shown in Fig. 9. The comparison between the results suggests that the local peak distribution can give a similar correction factor as found using MD simulation.
Additionally, the results of the original contact theory (without correction factor) compared with the modified version using local peak distribution for sliding velocity of 0.1 m/s under three different nominal pressures are shown in Fig. 10. As expected, the modified theory predicts a larger real contact area.

Modification of contact area based on improved GW method
In the second method, the effective height profile of the rough surface was calculated using the improved GW (IGW) contact theory. The advantage of this method compared to the local peak distribution method is that the material properties and normal load are also considered to determine the effective height profile of the rough surface. Therefore, the shift factor does not only depend on the surface roughness parameters, but also on the nominal pressure and the sliding velocity (in the case of viscoelastic materials). Figure 11 shows the comparison between the IGW-based correction factor (IGW-CF) and Persson fitted correction factor (PF-CF) for compound A with the sliding velocity of 0.1 m/s and three different nominal pressures on a 120-grit sandpaper surface. The shift factors for nominal pressures of 0.3 and 0.5 MPa found to be 2.02 and 1.83, respectively, which are close to the constant shift    factor found from the local peak distribution method. Therefore, both asperity-based methods give results similar to Persson's correction factor. However, the shift factor for a nominal pressure of 0.1 MPa was found to be 2.46, which predicts a smaller penetration depth compared to the other methods. The comparison between the original contact model and the modified version using the IGW method for sliding velocity of 0.1 m/s and different nominal pressures is shown in Fig. 12.
As expected, the IGW based modified contact model predicts a much larger real contact area at low normal pressures. The separation distances calculated using the IGW method and the corresponding shift factors as functions of nominal normal pressure and sliding velocity are shown in Figs. 13 and 14, respectively. When the normal pressure increases, the separation distance between the nominal contact surfaces decreases, and consequently, the shift factor decreases. On the other hand, when the sliding velocity increases, the elastic modulus of rubber also increases as the rubber becomes stiffer. Therefore, the separation distance between the two surfaces increases, and a larger shift factor must be implemented in the contact model. The numerical results suggest that for a wide range of velocities and intermediate normal pressures, the average shift factor is about two, and it gives similar results as the correction factor based on MD numerical simulation.

Modified friction model: Theory and experiments
The modified displacement assumption and contact theory, obtained from the IGW method, were implemented in the theory of rubber friction as expressed in Eqs. (22) to (24) excluding Persson's correction factor S(q). Then, the numerical results were compared with some experimental results obtained from a sliding friction tester. The details of the test set-up and the experimental procedure can be found elsewhere [46]. In this experiment, a small rubber block of compound B was tested on the asphalt surface with the power spectrum shown in Fig. 7(b). All the parameters used in the theoretical model were summarized in Table 1.
In the velocity range used in the experiment, it is expected that the hysteresis friction is the dominant | https://mc03.manuscriptcentral.com/friction mechanics of friction. Therefore, the displacement assumption has a significant effect on the theoretical results. The friction test was repeated several times for each sliding velocity for at least 3 m sliding distance. The experimental friction coefficients (blue dots) and the average friction coefficient (red marks) for each velocity are shown in Fig. 15 along with the theoretical model calculated based on (a) original multiscale contact theory and (b) the modified version of the theory. Red and orange dashed lines show the contributions of hysteresis friction using Eq. (22) and adhesive friction using Eq. (23), respectively.
As expected, the original model overestimates the contribution of hysteresis friction and predicts a higher value of the friction coefficient. On the other hand, the corrected model gives a better prediction of the real contact area and friction coefficient. Therefore, the numerical results of the modified friction model are in good agreement with the experimental results.

Conclusions
An asperity-based approach was introduced to modify the displacement approximation in Persson's multiscale theory of contact mechanics and rubber friction. This approach is applicable to contact mechanics of self-affine fractal surfaces, including various types of sandpapers and asphalt textures. In this approach, an affine transformation is used to shift the height profile of the rigid rough surface to   the upper regions of the asperities, which effectively contribute to the displacement field of the deformable body at the macro scale. This affine transformation reduces the effective height profile without altering the fractal properties of the surface. Two asperity-based methods were studied to find the shift factor for the affine transformation. The first method uses the height profile distribution to find the local peaks of asperities at the largest length scale of the surface roughness. The basic assumption in this method is that the height distribution of local peaks mainly contributes to the displacement field at the contact area. The main drawback of this method is that it only considers the height profile of the surface and neglects the material properties to find the shift factor. On the other hand, the second method uses material properties in addition to surface roughness parameters to find the shift factor for the affine transformation. It finds the effective penetration depth of the deformable body at the largest length scale using the improved Greenwood-Williamson (IGW) contact theory, which includes the interaction between asperities in an approximate way. Then, this penetration depth is used to transform the height profile so that the maximum vertical displacement of the smooth deformable surface becomes equal to the large-scale penetration depth. These two methods can also be combined to obtain a better approximation of the vertical displacement field on the contact area when the surface profile has a non-Gaussian height distribution or has significant skewness. In this condition, at first, the original nonsymmetric height profile must be shifted using local peak distribution to eliminate the deep valleys, which do not contribute to the contact mechanics. Then, the shifted profile can be used within the IGW contact theory to find the penetration depth at the largest length scale.
The numerical results showed that the asperitybased correction factors are in good agreement with the factor obtained based on fitting the theoretical results to MD numerical simulation. The local peak distribution method can perfectly match the fitting factor when the local peaks are found appropriately. On the other hand, the results of the IGW method deviate slightly from the MD-based correction factor at very low normal pressure as it predicts smaller penetration depth. However, it is also not clear in what conditions and to what extent the MD based factor can effectively correct the model. Besides, different values of γ for Eq. (8), though in the limited range of 0.4 to 0.5, were reported in Refs. [16,32,44,50], which can give slightly different results. The effect of nominal pressure and sliding velocity on the penetration depth and the shift factor in the IGW method were also studied. The results suggest that, on average, the vertical displacement on rough surfaces, such as sandpaper and asphalt, is approximately half of the height profile for intermediate pressures and a wide range of velocities. Therefore,  2 z h u is a better approximation than  z u h in the contact model, at lease for the case study done in this paper.
Finally, the modified displacement and contact models were implemented in the theoretical model of rubber friction. The results of both modified and original theoretical models were compared with the experimental results of friction between rubber and Friction 9(6): 1707-1725 (2021) | https://mc03.manuscriptcentral.com/friction asphalt at different sliding velocities. The comparison of the results indicates that the modified theoretical friction model is in good agreement with the experiment, while the original theory (without the parameter S(q)) overestimates the hysteresis friction component.

Open Access This article is licensed under a Creative
Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Anahita EMAMI. She received her Ph.D. degree in engineering mechanics in 2018 from Department of Biomedical Engineering and Mechanics at Virginia Polytechnic Institute and State University. She is currently an assistant professor at Texas State University. Her research interests include contact, friction, and wear of flexible materials, stretchable sensors, and experimental and numerical stress analysis.