Elastic Rough Surface Contact and the Root Mean Square Slope of Measured Surfaces over Multiple Scales

: This study investigates the predictions of the real contact area for perfectly elastic rough surfaces using a boundary element method (BEM). Sample surface measurements were used in the BEM to predict the real contact area as a function of load. The surfaces were normalized by the root-mean-square (RMS) slope to evaluate if contact area measurements would collapse onto one master curve. If so, this would conﬁrm that the contact areas of manufactured, real measured surfaces are directly proportional to the root mean square slope and the applied load, which is predicted by fractal diffusion-based rough surface contact theory. The data predicts a complex response that deviates from this behavior. The variation in the RMS slope and the spectrum of the system related to the features in contact are further evaluated to illuminate why this property is seen in some types of surfaces and not others.


Introduction
Practically, all surfaces are rough over many scales, and when brought together, what is actually touching is a fraction of the surface area. The reduced real contact area due to roughness also governs the friction and wear between the surfaces [1][2][3][4][5][6][7][8][9][10][11][12]. It is now well accepted that roughness results in the real contact area being a proportional load, which results in the macroscale phenomena of the friction force being proportional to the normal force (i.e., the friction coefficient). This real contact area and small contact spots also result in the phenomena of contact or spreading electrical and thermal resistance. The contact spots occurring between the peaks or asperities act as bottlenecks of the electrical current and heat conducted between the surfaces. This results in spreading resistance to occur between the contacting surfaces. Many rough surface contact models exist to predict the spreading resistance [13][14][15][16][17][18][19][20]. Although this work focuses on the total real contact area, the area of the individual asperities can be estimated by dividing by the number of asperities. The number of asperities can be statistically calculated using several ways [14,[21][22][23] or counted individually from a deterministic rough surface contact model.
The real area of contact between rough surfaces is often assumed to only be a function of the average pressure normalized by the equivalent elastic modulus and the root-mean-square (RMS) slope [24,25]. This general relationship has been derived by several statistical and fractal rough surface contact methods, and although the relationship between the contact area and load is shown to be related solely by the RMS slope, the equations differ [26][27][28][29]. In these works, the surfaces are often modeled as a set of asperities of different scales and sizes based on the geometric shape (e.g., Archards' bumps on bumps) or as self-affine fractals [25,30]. When using the same idealized surface the majority of contact models developed from these assumptions converge, more or less, on to a single contact area vs. applied pressure plot [31]. A similar master curve is replicated experimentally by producing several self-affine fractal surfaces (via 3-D-printing) of varying roughness and performing in situ contact measurements [32]. However, not all models and experimental measurements produce the same results. Campana and Muser [33] noted a 20% variation between the contact area as a function of the pressure normalized by the RMS slope of real measured and random rough surfaces. The RMS slope is calculated by squaring the slope at each point on a rough surface, taking the mean value of these, and then taking the square root of the mean. These models based around the RMS slope appear to have become the most popular rough surface contact model, having been extensively cited, and are being applied to a wide variety of applications, including friction, wear, leakage, dynamic impact, and contact resistance.
Putignano et al. used deterministic modeling of elastic rough surface contacts to find the proportion between the normalized contact area and average pressure to be approximately 2 [34]. Additionally, a recent analysis of generated fractal surfaces suggests that the contact area is not solely dependent on the root mean squared gradient but must also factor in the bandwidth or Nayak parameter, α [35]. These discrepancies may tie back to the original surfaces many of the current models are based upon. Greenwood and Williamson [21], Whitehouse and Archard [36], McCool [37], and others have observed that some real surfaces have nearly Gaussian distributions in surface heights and are highly random, but they are never perfectly Gaussian nor perfectly fractal. Sayles and Thomas also showed that many different surfaces appear to follow a self-affine fractal structure but just approximately [38].
Recent work by Zhang and Jackson has found discrepancies between different mathematical theories of fractals and the parameter extraction methods employed [39,40]. Green also recently derived corrected equations to analytically relate statistical surface parameters to Weierstrass-Mandelbrot fractal surfaces [41]. Even more recently, researchers have used tunneling electron microscopy and atomic force microscopy to measure surface roughness and spectrums all the way to the atomic scale [42,43]. These measurements appear to suggest that surfaces do not usually follow a fractal description over all scales. Many other researchers have questioned the applicability of assuming a fractal geometry for rough surface contact since it began usage in the 1990s [44][45][46][47][48]. Nonetheless, the main focus of this work is on the usage of the RMS slope to characterize and predict rough surface contact behavior, whether the surfaces are fractal or not.
In this work, we examine the predicted relationship between contact area and force with real (measured) imperfect surfaces using a boundary element method. If the fractal and statistical-based theories mentioned above are correct, by using surfaces that each have the same RMS slope, these theories predict the full collapse of the contact area over the pressure data. If they do not, it suggests that additional parameters should at least be considered to refine contact results and emulate manufactured surfaces. As a preliminary step toward identifying such parameters, an examination of the hierarchy of the multiscale surface features on the surface will be explored and compared to ideal spectrums.

Methods
The same measured surfaces used in Jackson and Green [49] are again employed here to consider the effectiveness of the RMS surface slope on normalizing the contact area to pressure function. Their measured surface properties are listed Table 1. The properties of the original surfaces are given in Table 1. σ is the root mean square roughness and m n are the spectral moments of the surface. B is the ratio of the amplitude to the wavelength ratio taken from a Fourier transform of the surface, and B ave is the average value and B max is the minimum value. Finally, g is the root mean square slope or gradient, which is of course related to m 2 by: The 1001 nodex1001 node surfaces in [49] were reduced to 512 × 512 in the current work so that an existing boundary element method [50,51] and code can be used to solve the contact problem. Two additional experimentally measured surfaces are added and also considered (Surfaces 4 and 5). They were also scaled to be 100 micrometers in width on each side. These surfaces can be described as (1. mirror finished silicon surface, 2. unfinished machined steel surface, 3. a plasma sprayed coating, 4. the 4 L surface of a s-22 micro-finish comparator, and 5. a sand finished aluminum surface). Care was taken such that all of the surfaces' measurements were made in the range of resolution of the profilometry equipment.
Next, the surface heights of each were scaled by a single factor such that all the surfaces had an RMS slope of unity. This same factor also adjusts the roughness by the same proportion. The resulting surfaces have the parameters listed in Table 2. Note that the values of m 2 and g are both now unity due to the normalization. In addition, many of the other parameters are now within an order of magnitude. The effective elastic modulus (E*) between the surfaces is also set to unity. The surfaces are considered to be perfectly elastic in nature. As mentioned previously, a boundary element model is used to predict the real contact area of the surfaces as they are effectively brought into contact with a frictionless rigid flat surface.

Results and Discussion
The predictions of the boundary element method (BEM) method for the contact area as a function of the average contact pressure of each surface are then compared. The real contact area, A r (what is actually touching), is normalized by the nominal or apparent area, A n . In this case, the nominal or apparent area, A n , is the entire surface. The applied normal load is also normalized by the nominal or apparent area, A n , and the effective elastic modulus, E*. If the results follow the theory mentioned previously, the predictions should collapse onto one curve when plotted this way. As shown in Figure 1, the real area of contact curves as a function of the normalized pressure do not collapse onto one curve. The five surfaces immediately start to diverge after the surfaces are put into contact and do not begin to converge until they reach a dimensionless contact pressure of 2, which is associated with over 90% of the surfaces in contact. Surface 2 had the strongest response as a function of the applied pressure and surface 1 and 5 had the weakest responses over different ranges. Although the differences may visually appear small, they are as large as 50% in some locations on the curves. Note also that surfaces 1 and 2 have relatively little agreement, even though they have nearly identical values of the RMS slope as well as the Nayak parameter, α, which in a recent work [52] was stated to be the two most important parameters governing rough surface contact. The differences in the results might be due to the anisotropic and random nature of real measured surfaces that preclude perfectly statistically Gaussian or fractal self-affine surfaces.
ct. 2021, 5, x FOR PEER REVIEW 4 of 13 50% in some locations on the curves. Note also that surfaces 1 and 2 have relatively little agreement, even though they have nearly identical values of the RMS slope as well as the Nayak parameter, α, which in a recent work [52] was stated to be the two most important parameters governing rough surface contact. The differences in the results might be due to the anisotropic and random nature of real measured surfaces that preclude perfectly statistically Gaussian or fractal self-affine surfaces. The differences might therefore be explained by how the structure of the surface changes as the scales are traversed. The spectrum of the surfaces is plotted in Figures 2  and 3. A spectrum is a description of a signal's (here a surface) amplitude at different harmonic scales, or wavelengths (frequency is the inverse of wavelength and defines the number of peaks per length). In this work, the spectrum was calculated using a fast Fourier transform. Here, the spectrums are plotted in two ways: first, the spectrums are plotted simply as the amplitude as a function of the wavelength (see Figure 2); next, they are plotted by normalizing the amplitudes of the Fourier series by the wavelength at each scale (see Figure 3). In previous literature, the aspect ratio (amplitude to wavelength ratio) has been assigned the symbol, B [17,53]. If B were constant across all scales, this indicates that a surface's asperities (the peaks on a rough surface) have the same aspect ratio or shape across all scales. This structure would indicate a self-similar fractal surface. Figure  3 demonstrates that the normalized amplitude for each surface is divergent from each other for wavelengths exceeding 10 −6 m and non-constant at all length scales. For most of the surfaces, the ratio B appears to increase as the scale (wavelength) decreases. Surface 5 does have a shoulder in the data spanning between the 300-nm and 10-µm wavelength, which you would expect in a self-similar structure.
In addition, the power spectral density of each surface is calculated and plotted as shown in Figure 4. The complete original PSDs of each surface is plotted in Appendix A as well. Overall, all the PSDs are not isotropic (axisymmetric), but their radially average The differences might therefore be explained by how the structure of the surface changes as the scales are traversed. The spectrum of the surfaces is plotted in Figures 2 and 3. A spectrum is a description of a signal's (here a surface) amplitude at different harmonic scales, or wavelengths (frequency is the inverse of wavelength and defines the number of peaks per length). In this work, the spectrum was calculated using a fast Fourier transform. Here, the spectrums are plotted in two ways: first, the spectrums are plotted simply as the amplitude as a function of the wavelength (see Figure 2); next, they are plotted by normalizing the amplitudes of the Fourier series by the wavelength at each scale (see Figure 3). In previous literature, the aspect ratio (amplitude to wavelength ratio) has been assigned the symbol, B [17,53]. If B were constant across all scales, this indicates that a surface's asperities (the peaks on a rough surface) have the same aspect ratio or shape across all scales. This structure would indicate a self-similar fractal surface. Figure 3 demonstrates that the normalized amplitude for each surface is divergent from each other for wavelengths exceeding 10 −6 m and non-constant at all length scales. For most of the surfaces, the ratio B appears to increase as the scale (wavelength) decreases. Surface 5 does have a shoulder in the data spanning between the 300-nm and 10-µm wavelength, which you would expect in a self-similar structure. However, a completely self-affine fractal surface would not experience the drop off seen at higher wavelengths, and would be represented by a straight line without deviation. The deviation seen here for the normalized amplitude for the five RMS-equal surfaces indicates a likely vector for the divergence between traditional models and real surfaces.  In addition, the power spectral density of each surface is calculated and plotted as shown in Figure 4. The complete original PSDs of each surface is plotted in Appendix A as well. Overall, all the PSDs are not isotropic (axisymmetric), but their radially average PSD is quite similar to the previous spectrum plots. From Figure 4 (also Figure 3), it is clear that surfaces 3-5 might be considered representative of self-affine rough surfaces, which follow a piece wise power law. Surfaces 1 and 2 are different, their high-frequency PSD is flattened, which is not as often observed. This might be the reason why the contact results of surfaces 1 and 2 significantly deviate from the rest of the surfaces. Additionally, the gap between the real contact area of surface 1 and 2 could be due to the difference of their PSD curves.   However, a completely self-affine fractal surface would not experience the drop off seen at higher wavelengths, and would be represented by a straight line without deviation. The deviation seen here for the normalized amplitude for the five RMS-equal surfaces indicates a likely vector for the divergence between traditional models and real surfaces.
To further explore the potential link between variable RMS roughness and surface contact response, the RMS slope (g) is plotted as a function of the surface height in Figure 5. For a certain truncation height (h), only the portion of the surface taller than h is used to calculate the RMS slope (g). As shown in Figure 5, the RMS slope (g) varies based on the truncation height (h). Since only the taller peaks on the surface will be in contact at lighter loads and the rest of the surface will only move into contact as the load increases, it stands to reason that the RMS surface roughness will only be dependent upon the surface area currently in intimate contact. In this way, the RMS surface roughness will vary throughout the loading of a real surface. We can call the RMS slope of the contacting surface the effective RMS slope. The data in Figure 5 shows that the effective RMS slope will start at a higher value for lightly loaded contacts and decrease gradually until approaching the RMS slope of the entire surface for heavier loads as expected. The effective RMS slope appears to approach the RMS slope of the entire slope at an h value of approximately -σ or −2 σ, Fractal Fract. 2021, 5, 44 7 of 12 which is just after passing the mean height of the surface. The trends for each surface are complex, containing plateaus but are roughly monotonic. It is likely that this variation in the slope during loading is a direct cause of the disparity between the contact response of real imperfect surfaces and computational models, and is a result of these natural surfaces' lack of perfectly self-affine features.
loads and the rest of the surface will only move into contact as the load increases, it stands to reason that the RMS surface roughness will only be dependent upon the surface area currently in intimate contact. In this way, the RMS surface roughness will vary throughout the loading of a real surface. We can call the RMS slope of the contacting surface the effective RMS slope. The data in Figure 5 shows that the effective RMS slope will start at a higher value for lightly loaded contacts and decrease gradually until approaching the RMS slope of the entire surface for heavier loads as expected. The effective RMS slope appears to approach the RMS slope of the entire slope at an h value of approximately -σ or −2 σ, which is just after passing the mean height of the surface. The trends for each surface are complex, containing plateaus but are roughly monotonic. It is likely that this variation in the slope during loading is a direct cause of the disparity between the contact response of real imperfect surfaces and computational models, and is a result of these natural surfaces' lack of perfectly self-affine features. Contact theories incorporating self-affine features remain one of the most critically important tools we use to understand how surfaces will generally behave when creating an interface. It is not surprising that the data set that is most nearly self-affine (surface 5) sat in between all the other data sets at appreciable loads. However, there is always room for refinements. Moving forward, it would be interesting to compare and evaluate the many existing rough surface contact mechanics models when real surface measurements are employed. Contact theories incorporating self-affine features remain one of the most critically important tools we use to understand how surfaces will generally behave when creating an interface. It is not surprising that the data set that is most nearly self-affine (surface 5) sat in between all the other data sets at appreciable loads. However, there is always room for refinements. Moving forward, it would be interesting to compare and evaluate the many existing rough surface contact mechanics models when real surface measurements are employed.

Conclusions
This work considered five different measured rough surfaces that are scaled such that they all have an equal RMS slope. A BEM method was used to predict the real area of contact for these surfaces as a function of load. Based on existing fractal-based contact theories, the predicted real area of the contact versus load curve of these surfaces should collapse onto one unified curve. However, they do not. This work also provides several possible reasons for this deviation. Measured surfaces are not perfect self-affine fractals that many of these theories are based on, and therefore this is probably why the theories do not hold for measured surfaces. The spectrums of the measured surface clearly show that they are not perfect fractals because the curves are not linear. The real contact area as a function of the load curves for surfaces 1 and 2 do not agree, although both their RMS slope and Nayak parameter, α. In addition, it is also shown by using a truncation approximation that the RMS slope of the surface portion actually in contact (the real area of contact) can actually change non-linearly with the load and area of contact. This too provides an explanation as to why the real contact area vs. load curves do not collapse. Data Availability Statement: Some data is available in a previously published work. Additional data will be made available in a data depository.
Acknowledgments: Thank you to Auburn University and Hefei University for providing the resources to support this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Power Spectral Densities of Measured Surfaces.