Rough Surface Contact Modelling—A Review

: It has been shown experimentally that boundary friction is proportional to load (commonly known as Amontons’ law) for more than 500 years, and the fact that it holds true over many scales (from microns to kilometres, and from nano-Newtons to Mega-Newtons) and for materials which deform both elastically and plastically has been the subject of much research, in order to more fully understand its wide applicability (and also to ﬁnd any deviations from the law). Attempts to explain and understand Amontons’ law recognise that real surfaces are rough; as such, many researchers have studied the contact of rough surfaces under both elastic and plastic deformation conditions. As the focus on energy efﬁciency is ever increasing, machines are now being used with lower-viscosity lubricants, operating at higher loads and temperatures, such that the oil ﬁlms separating the moving surfaces are becoming thinner, and there is a greater chance of mixed/boundary lubrication occur-ring. Because mixed/boundary lubrication occurs when the two moving rough surfaces come into contact, it is thought timely to review this topic and the current state of the theoretical and experimental understanding of rough-surface contact for the prediction of friction in the mixed/boundary lubrication regime.


Introduction
Leonardo da Vinci's work on tribology was discussed by Dowson [1], and more recently, Hutchings [2] and Betancourt-Parra [3] have reviewed the many friction experiments carried out by him between 1493 and 1515. In addition, Sawyer [4] reported on da Vinci's work on wear. According to Hutchings [2], da Vinci first stated the "laws" of sliding friction in 1493 (based on his experiments), and these were: (1) the force of friction between two sliding surfaces is proportional to the load pressing the surfaces together, and (2) the force of friction is independent of the apparent area of contact between the two surfaces. The first law suggests that the force of friction, F, can be written as F = µW, where W is the applied load, and µ is a constant of proportionality which is known as the friction coefficient. Leonardo da Vinci found that for unlubricated sliding bodies, µ was approximately 1 /4. Leonardo da Vinci's work was written up in his personal notebooks, but it is not clear how widely known his work was (at the time) in the broader scientific community. Amontons [5] was the first to publish the "laws" of friction in a scientific journal, in 1699, and these became known to a much wider audience. Amontons [5] found a value for µ of 1 / 3 . Following on from Amontons, Coulomb [6] also worked extensively on friction in machines. A useful recent historical discussion of the work of Amontons and Coulomb can be found in reference [7].
There has been much discussion as to the origin of Amontons' empirical law and why it applies so widely over many length scales, and to systems in which there is plastic deformation, and also in systems where the deformation is elastic [8][9][10][11][12].
It is, of course, important to realize that real surfaces are rough. Asperities are randomly distributed (up and down) relative to the "centre-line average" (CLA) position. When two such surfaces touch, it will be the highest asperities on each surface that will be in contact, and so the actual area of contact is very much less than the apparent (geometric) area, and for a given load, the contact pressure at the highest asperities will be the largest. Initial attempts to explain Amontons' first law [13] suggested that the force needed to "push" the surface over an asperity with angle θ would be equal to the normal load (W) multiplied by tan (θ), so that if F is the friction force, then F/W = µ = tan(θ). However, it was later realized that the energy expended on dragging an asperity to the top of another is recovered when it falls down the other side [14], so no energy is lost, and the two surfaces should continue to move (once they are set in motion) without the need for any other driving force. Therefore, some other "energy dissipating mechanism" is needed to explain friction.
Bowden and Tabor [15], in their 1930s work on metal-metal contact, assumed that the contacting asperities would deform and flow plastically such that the real contact area would increase until the contact pressure drops below the material yield stress. In this model, the highest asperities are "wiped away", the real contact area is proportional to load, and the contact pressure is constant and equal to the material yield stress. This type of model assumes that material hardness is proportional to the material yield stress. The friction force was calculated by multiplying the real contact area by the shear stress, and so Amontons' law of friction was explained by the fact that the real contact area was proportional to the applied load. Although this explanation seems reasonable for plastic deformation, many machines are known to operate successfully for thousands of hours, and it is difficult to see how these lifetimes can be maintained if the rough surface asperity contact is always plastic.
At first sight, if rough surfaces deform elastically, this seems to disagree with Amontons' first law of friction, as Hertz [16] has shown that for a single elastic contact of radius R contacting an infinitely hard flat surface, the area of contact would be proportional to W 2/3 , where W is the applied load. However, in practice, a rough surface has a distribution of asperity heights, and a more detailed analysis by Greenwood and Williamson [17] of rough, non-adhering surfaces as a collection of asperities of varying height, with spherical tips, deforming elastically, found a linear relationship between the real contact area and load. This initial analysis by Greenwood and Williamson has since been refined and extended by many other researchers [18][19][20][21][22][23][24][25][26][27]. The analysis of Greenwood and Williamson [17] was substantially aided by earlier useful insights from Archard and others [28][29][30][31][32][33]. In particular, in reference [30], Archard showed that if there is one level of "protuberances", of radius R 2 , superimposed on a spherical tip of radius R 1 (with R 2 << R 1 ), then elastic deformation leads to the real area of contact being proportional to W 8/9 (and that if a further level of "protuberances" are added, of radius R 3 , with R 3 << R 2 , then the real contact area is proportional to W 26/27 ). Archard's work [30] showed that the real area of contact can be almost proportional to the load even when there is elastic contact only, for multi-asperity contact. It should also be mentioned that Archard's concept of protuberances on protuberances on protuberances clearly has a link to fractals, where multiple length scales are important, but at the time of Archard's work, fractals were not yet known about, as they were only introduced in the 1970s [34]. The consensus of these various studies is that, in the elastic contact case, the real contact area is approximately proportional to the applied load, and therefore Amontons' first law of friction follows. Persson [8] has reported that the real contact area is also approximately proportional to the applied load even when there is adhesive interaction between surfaces.
The current picture then seems to be as follows: the rough surfaces initially deform plastically on contact, and the highest asperities are "wiped away". This then leads to a "run-in" surface, which has a different asperity height distribution compared to the original surfaces, which leads to lower contact pressures, and which can sustain elastic deformation for many cycles (which explains the longevity of many machine elements). A useful reference on "running-in" may be found in [35]. For both plastic and elastic deformation, then, the real area of contact is approximately proportional to the load, and so Amontons' first law of friction is obeyed (although the coefficient of friction could be different in the plastic and elastic stages).
Despite the simplistic picture above, in a useful review of Amontons' law, Gao et al. [9] commented that because of the various assumptions made in the Bowden Tabor (BT) [15] and Greenwood Williamson (GW) [17] theories, they doubted that both theories could account for all of the situations where Amontons' law has been found to hold, and fundamental insights into why the law holds, which are less model specific, are still needed, in the opinions of those authors. In that paper, Gao et al. [9] also found, using surface force apparatus, atomic force microscopy, and molecular dynamics simulations, that there is not a linear relationship between the friction force and the normal load at the local level (i.e., at the scale of individual asperities), but that there is a linear relationship between these quantities for their time-and space-averaged values (i.e., when averaged over many asperities). Gao et al. [9] also stated that although mechanical energy is converted to heat at a steady rate, when averaged over the whole macroscopic junction, locally, at any instant, the properties are far from their mean steady-state values. Thus, even in the steady state, the local tribological parameters fluctuate widely, are largely uncorrelated, and may not even average out to the system average. This review is organised as follows. In the next section, an overview will be given of how the friction coefficient of lubricated contact is determined for a rough surface, followed by a section considering the statistics of rough surfaces, and then a section on a more detailed discussion of rough-surface contact models. The importance of lubricant additives is discussed, and then there is a summary of recent experimental data on rough-surface contact friction, followed by a section on how the models can be used by engineers and tribologists to estimate mixed/boundary friction. The review will finish with discussion and conclusions sections.
It should be mentioned that this topic has been the subject of previous reviews, for example by McCool [36] in 1986, Bhushan [37] in 1998, Liu et al. [38] in 1999, Adams et al. [39] in 2000, and more recently by Vakis et al. [40] in 2018, although the latter is a more general review of simulation in tribology. This review differs from earlier ones in that it reports more recent work on rough-surface contact modelling (up to 2021), includes recent experimental data on rough-surface contact friction, and is more focussed on how engineers and tribologists can apply such models in order to calculate mixed/boundary friction forces and power losses.

The Friction and Load Carrying Capacity of Lubricated Rough-Surface Contact
Assume that a rough surface is separated by a fluid film. If the centre-line averages of the rough surfaces are separated by a film thickness, h, and if the combined root mean square (RMS) roughness of the surfaces is σ, then a useful dimensionless tribological parameter is the lambda ratio, λ, which is defined to be h/σ. In tribology textbooks (see for example [41][42][43], it is generally considered that the surfaces are completely separated if λ > 3 (hydrodynamic lubrication), whereas if λ < 1, the contact is in the boundary lubrication regime. If 1 < λ < 3, then the contact is often assumed to be in the mixed lubrication regime.
If the total load carried by the contact is W TOTAL , this can be split between the load carried by the asperities W asp and the load carried by the fluid film W fluid as: The friction force, F TOTAL , of the contact can be written as where µ asp is the friction coefficient of the asperities (i.e., when λ = 0) and µ fluid is the friction coefficient of the fluid (i.e., when λ > 3). The friction coefficient, f, of the contact can then be written as f = F TOTAL /W TOTAL , which is equal to where X = W asp /W TOTAL , which can be regarded as the proportion of mixed/boundary lubrication. Clearly, it is expected that X = 1 when λ = 0 and X is approximately zero when λ > 3. It is of great interest, for the modelling of mixed/boundary friction, to know how X varies with λ, and it is also of experimental interest to measure this variation. Equation (3), above, has previously been reported by Spikes and Olver [44], and in that work they proposed that X was simply equal to (1 + λ) −2 .

The Statistics of Rough Surfaces
The fact that real surfaces are rough has been known for many years. However, it is only since the 1960s that digital methods have been used to quantitatively measure the roughness and other statistical properties of surfaces. Useful textbooks on surface roughness include those by Halling [43], Thomas [45] and Whitehouse [46].
Nayak published three impactful papers [47][48][49] on the statistics of surface roughness which were described by Greenwood [50] as being "widely recognized as the foundations of the theoretical study of surface roughness". Greenwood himself later published further papers on surface roughness (including "A Unified Theory of Surface Roughness" [51], and "Problems with Surface Roughness", [52], and also made useful comments [50] on Nayak's third paper [49]).
Nayak [47] focussed on random, isotropic, Gaussian surfaces, and treated the height of a rough surface as a two-dimensional random variable. Nayak's work built on the previous work of Longuet-Higgins [53,54], who studied random ocean surfaces. Nayak [47] found that there were three important parameters that fully characterized a random, isotropic, Gaussian surface, and these were m o , m 2 , and m 4 , where m o is related to the RMS roughness of the surface (σ), m 2 is related to the RMS surface gradient, and m 4 is related to the RMS radius of the curvature of the rough peaks.
Nayak's second paper [48] detailed the difference between profile and surface properties, and gave useful advice on how to reliably extract surface properties from surface profile measurements.
Nayak's third paper [49] considered how the number of contacts varies as the surfaces approach each other. Simple approaches to this problem assume that a contact is formed whenever a summit meets the opposing surface, and that the contact then grows independently of others. What is thought to happen, in practice, as rough surfaces approach is that contacts can grow and also merge with other nearby contacts, such that although the real contact area increases, the actual number of contacts can decrease (when the merging of neighbouring contacts occurs). In the third paper, Nayak looks at the number of closed contours on a rough surface at a given height. If the number of contact patches at a given height is significantly less than the number of summits above that height, then clearly some of the initial contacts have merged into composite contacts, such the assumption that contacts behave independently is likely to be wrong. Greenwood [50] later commented that any asperity contact model applied to the realistic range of separations (z < 3σ), where z is the separation of the surfaces and σ is the RMS surface roughness, must take into account the merging of asperities, the consequent reduction in the number of contacts, and the corresponding distortion of the contact shapes.
A useful summary of the statistics of rough surfaces was given by Greenwood in 2013 [55], in which he commented that measurements of surface roughness show that a nominally flat surface will always have a distribution of height variations, possibly superimposed on a longer wavelength "form" or regular waviness (often related to the way in which the surface was manufactured). Clearly, the height distribution could take many mathematical forms, but often a Gaussian distribution is followed: where h is the height above (or below) the nominal zero level, σ is the RMS surface roughness, and ϕ(h) is the probability density function. Greenwood also notes that although the surface roughness as a whole may not follow a Gaussian distribution function, the heights of the summits on a surface usually do (and it is, of course, the summits on a surface where the first contacts form). However, it is possible that the standard deviation of the Gaussian distribution of the summit heights could be different that of the surface as a whole.
In addition to the surface roughness summit heights, the asperity radius of the curvature is also of major importance. Numerous researchers have studied typical asperity radii of curvature and how these are statistically distributed [43,[56][57][58], and Greenwood [52] commented that the distribution of the asperity radius of the curvature follows a Rayleigh distribution, such that P(r) = r· exp − r 2 2 . One of the issues with the measurement of surface roughness is that the surface properties depend on the size of the measuring probe (or stylus). As mentioned previously, Nayak found there were three key parameters, m o (the RMS roughness), m 2 (the RMS surface slope) and m 4 (the RMS of the radius of the curvature of the asperities). As the size of the measuring stylus decreases, the curvature of the asperities, and their slope, increase (i.e., the asperity radius of the curvature decreases) although the RMS roughness stays relatively constant [57].
An alternative approach to the study of surface roughness is to focus on autocorrelation functions and power spectral densities. The autocorrelation function is defined by where z(x) is the height of the surface at point x and β is the shift in position at which the autocorrelation function is to be evaluated. This kind of analysis can also be extended to surfaces, where the height would be z(x,y), and Equation (5) would be extended to also include an integration in the y direction. It is often assumed that the autocorrelation function takes the form below: where, in the above equation, β o is a constant. If β = 0, then C(β) = σ 2 as required. The Power Spectral Density (PSD), as defined in Halling's textbook [43], is Using the autocorrelation function above (Equation (6)), it is found that For large values of k In an important contribution published in Nature in 1955, Sayles [58] found that G(k) was proportional to 1/k 2 for a wide range of surfaces spanning length scales differing by eight orders of magnitude. The factor β o that appears in the equations above can be regarded as the length scale of the measurement stylus, and so the values of m o , m 2 and m 4 will vary depending on the value of β o . It could be argued that in calculating the moments of the spectral density, rather than integrating between 0 and ∞, one should instead integrate between 2π/L 1 and 2π/L 2 , where L 1 is the size of the sample (typically a few mm) and L 2 is the size of the stylus (typically a few µm). It should be pointed out that the higher moments of the spectral density function do not exist for the spectral density function of Equation (9) (if the integral limits are 0 and ∞). The relationships between m o , m 2 , m 4 and the spectral density function are Additional useful details on the measurement of rough surface properties can be found in references [59][60][61]. It should be mentioned that research in this area is ongoing, and recently additional measurements were reported on surface roughness spanning eight orders of magnitude, with the lowest scale being a few Ångströms [62].
The main summary of this section is that a lot of information is known, or can be measured, on the properties of rough surfaces. Much of the published work has assumed that rough surfaces are generally isotropic, random, and described by a Gaussian probability distribution, although these assumptions can be relaxed if needed.

Bowden and Tabor's Model
Bowden and Tabor predominantly worked on metal-metal contacts [15]. The adhesion concept of friction for dry metal surfaces, put forward in 1734 by Desaguliers [63], was used successfully by Bowden and Tabor for metal-metal interfaces. However, at first sight, the idea that adhesion is involved in the friction process seemed to contradict experiments which found that friction was independent of the contact area. This contradiction disappeared when it was realized that it was the real contact area that was important in friction, not the apparent (geometric) area. Bowden and Tabor's experiments [15] showed that the friction between two metal surfaces was strongly dependent on the real contact area.
Bowden and Tabor [15] also distinguished between the abrasive wear of surfaces (which involves a harder surface sliding against a softer one) and adhesive wear, in which individual asperities adhere to each other, and will only plastically deform if a critical shear stress is applied.
Bowden and Tabor [15] assumed that the points at which contact occurred were like hardness indentations, and so the total area of actual contact would be W/H, where W is the normal load, and H is the hardness of the softer material (note that hardness has units of N/m 2 ). If a critical shear stress, τ s , was needed to shear the asperity junctions, then the frictional force would be (τ s W/H), and so friction is proportional to normal load. A useful retrospective look-back at this work was published by Tabor in 1996 [64]. Barber, in 2013 [26], commented that there were criticisms of this model, as although plastic deformation would be expected on first loading, the effects of work-hardening and residual stresses would reduce the likelihood of plastic deformation during subsequent loading, as for example, during sliding. Barber also commented that this argument was supported by wear measurements on sliding surfaces that suggested that the detachment of wear particles was a rare consequence of asperity-asperity interaction. Criticisms of the plastic deformation approach of Bowden and Tabor [15] motivated the research, by Archard [30], Greenwood and Williamson [17], and others, that attempted to derive an explanation of Amontons' law that did not depend on the assumption of plastic deformation.

Archard's Work
Archard published a number of influential papers on rough-surface contact in the period from the 1950s to the 1980s [28,30,31,[65][66][67][68]. The most influential of these papers, regarding the realization that the elastic deformation of rough surfaces could lead to a real contact area that was approximately proportional to load, was his 1957 paper [30] "Elastic Deformation and the Laws of Friction". This paper considered what happened when asperities at multiple length scales contacted flat surfaces. Archard explained that for a single asperity of radius R 1 , contacting a smooth surface under a load W, there will be a single circular area, A 1 , of radius b, where Archard then said that the load δW supported by an annulus or radius r with width δr is given by Archard then considered that the initial asperity (of radius R 1 ) was covered in smaller asperities, the radii of which were R 2 (where R 2 << R 1 ), which were evenly distributed over the larger asperity with a number m per unit area. The annulus load above was assumed to be unchanged, with the load δW supported by q "protuberances" (with q = 2πrmδr) and that the load w 2 on each one being given by The area a 2 of each of these contacts is where k 2 is the constant of proportionality appropriate to a sphere of radius R 2 , and the area of contact within the annulus is qa 2 , such that the total contact area, A 2 , is The addition of this second layer of asperities results in a total contact area proportional to W 8/9 , whilst the area of a single asperity is proportional to W 2/3 . Archard [30] showed in the same paper, using a similar analysis to that above, that the addition of a third level of asperities (with radius R 3 , where R 3 << R 2 << R 1 ) results in a total contact area that is proportional to W 26/27 .
The essence of Archard's argument [30] is that once it is accepted that there are asperities on multiple different length scales, then the real area of contact is approximately proportional to area, even if the contact is elastic.

The Greenwood-Williamson and Greenwood-Tripp Models
Greenwood and Williamson's highly cited 1966 [17] paper titled "Contact of Nominally Flat Surfaces" built upon Archard's insights [30], but instead of assuming a distribution of asperity radii, they considered asperities with the same radius of the curvature but a differing height distribution. A flat surface was assumed to be covered with asperities of which the density was η per unit area, all the same size, except that their heights followed a statistical distribution function ϕ(z) for which the standard deviation was σ (with the asperities assumed to have the same radius of the curvature, R). When the surface asperities contacted a flat surface, the nominal surface area of which is A, and when the distance between the flat surface and the mean plane of the asperity laden surface is h, contact will only occur at those asperities of which the height z exceeds h. The number of such asperities is given by where such contact does occur, the compression of the asperity, δ, will be equal to (z − h), and so according to Hertzian elastic contact theory, the Hertzian radius, a, is given by a = (Rδ) 1/2 , and the Hertzian load, P, is P = 4E*R 1/2 δ 3/2 /3 (where E* is the effective combined elastic modulus of the surfaces). Therefore, the total (real) contact area, A r , is given by: The total load P will be: If σ is the RMS surface roughness, then the above relationships can be written as: ·ϕ * (s)·ds (22) where in the equations above, ϕ*(s) is the standardized height distribution. Two cases were considered by Greenwood and Williamson [17]. In the first case, they assumed an exponential probability distribution, such that ϕ*(s) = exp(-s). It was possible to solve the above integrals exactly with this assumption, with the results below: In this particular example, the real contact area (A r ) is directly proportional to the load (P), and in addition, the number of contact points (N) is also directly proportional to the load. Greenwood and Williamson [17] pointed out that an exponential distribution of summit heights is a good approximation to a Gaussian distribution for the important range from z = σ to z = 3σ.
The second example considered was when the standardized height distribution was In this case, the integrals in Equations (21) and (22) were not solved in a closed form, but were reported as Table 1 lists the values of these functions for various values of n, taken from reference [21]. Clearly, these functions drop off very rapidly with λ (where λ = h/σ) and for λ > 4 the functions are essentially zero, as expected, as for λ > 3 it is generally assumed that surfaces are completely separated. Recently, Jedynak [69] published the exact solution below for F n (λ) for the case where 2n is an integer (which applies to the interesting cases of n = 3/2 and n = 5/2): where Г(x) is the Gamma function, and U(m,x) is the parabolic cylinder function [70]. A subsequent paper by Greenwood and Tripp, in 1970 [21], extended the analysis to allow for the contact of two rough surfaces (the Greenwood and Williamson paper [17] only considered the contact of one rough surface against a flat surface). Greenwood and Tripp [21] found that the load P carried by the asperities was where, in the above equation, β is the asperity radius of the curvature, d is the separation between the mean level of the surfaces, and A is the geometrical area of the contact. This equation is still widely used today to estimate the load carried by asperities when a contact is in mixed or boundary lubrication. An exact solution to the Greenwood and Williamson model, assuming a Gaussian surface roughness distribution, was published in 2011 by Jackson and Green [71].
It should be noted that in the Greenwood and Williamson [17] and Greenwood and Tripp [21] type models, the fluid between the moving surfaces is effectively just a "spacer layer". However, it is known that in elastohydrodynamically lubricated contacts, the high pressures elastically distort the surfaces, and significantly increase the viscosity of the fluid. In addition, the high pressures can reduce the amplitude of the roughness peaks (see for example references [72,73]). As far as the author of this review is aware, this effect is not accounted for in the majority of rough-surface contact models published to date.
Greenwood and Tripp [21] concluded that if the asperity peak height distribution is Gaussian, then the mode of deformation, the asperity shape, and whether the asperities are on one or both surfaces are all relatively unimportant, as in all cases the real area of contact is proportional to the load.
In a useful paper published in 2001, Greenwood and Wu [25] took a critical look back at some of the assumptions used in the Greenwood and Williamson model [17]. It was concluded that the theory contained the essential feature needed-asperities have a height distribution that approximates to a simple exponential in the relevant range of heights. They did not believe that the assumption of a single radius of the curvature was a major issue for the theory. One aspect of the original model that was, with the benefit of hindsight, considered wrong was the assumption that "peaks" on a measured surface profile (points higher than their immediate neighbours at the sampling interval used) correspond to asperities. This assumption gives a false impression of both the number and radius of the curvature of the asperities.
Another couple of comments are worth highlighting. It was later pointed out, by Whitehouse and Archard [18], that the asperity parameters used in the Greenwood and Williamson [17] and Greenwood and Tripp [21] papers are not independent, and that the quantity (ηβσ) that appears in the models should be treated as a constant, with typical values in the range 0.03 to 0.05.
Greenwood and Williamson [17] also introduced a plasticity index, Ψ, defined as where E* is the combined plane stress modulus for the two surfaces (Pa), H is the hardness (Pa), σ is the combined RMS surface roughness (m), and R is the mean asperity radius of the curvature. Greenwood and Williamson [17] studied numerous different surface topographies, and found that Ψ could vary from around 0.25 up to 30, and for Ψ > 1, the surfaces would primarily deform plastically, whilst for Ψ < 0.6, the surfaces would deform elastically, such that it is only in the region 0.6 < Ψ < 1 where the contact type (elastic or plastic) is uncertain. Greenwood and Wu [25] later commented that the condition for elastic contact is lκ < 2H/E*, where l is the asperity length, κ is its curvature, H is the hardness, and E* is the plane-stress modulus. Small asperities have much smaller radii of curvature, such that their curvature κ will be greater, and so it seems much more likely that small asperities will be deformed plastically, whereas the larger contacts will be deformed elastically (as was essentially pointed out by Whitehouse and Archard [18] in 1970).

Bush, Gibson, Thomas Model
Bush, Gibson and Thomas [19] extended the Greenwood and Williamson analysis [17] by allowing not only for a statistical distribution of asperity heights, but also a statistical distribution of asperity radii. In addition, they assumed that the asperity radius could be different (in the x-and y-directions), and so their 2D surface was effectively covered with ellipsoids. They used the Nayak surface roughness parameters m o , m 2 and m 4 in their work, and found a linear relationship between the load and the real contact area. They also found that the real contact area, A r , was related to λ by

Models Due to Persson
From 2001 onwards, Persson published models [22,23] that approached the problem from a different direction. The works of Archard [30], Greenwood and Williamson [17], Greenwood and Tripp [21], and Bush, Gibson and Thomas [19] all assumed that loads were relatively low, asperity contacts were elastic, and the real area of contact was a relatively small fraction of the apparent (geometric) area. Persson considered the problem of rough-surface contact from the opposite limit, in which he assumed the real area of contact was large (his work was motivated by contact between rubber surfaces). In fact, Persson commented that his theory [22,23] is exact in the limit of full-contact conditions, although this has been debated [74]. A parabolic partial differential equation was used to describe the interfacial pressure, and the diffusivity term in the differential equation was approximated by assuming that the Power Spectral Density (PSD) of the elastically deformed surface was equal to that of the undeformed rough surface. However, the stored elastic energy is calculated on the assumption that the PSD of the deformed surface is equal to the product of the PSD of the undeformed surface times the fraction of contact area. It was found that there was a linear relationship between the real contact area and load, and this finding is more robust than earlier results from Greenwood and Williamson [17] etc., as those earlier findings were only valid for a relatively small real area of contact (5-10% or less). Persson [75] also reported that the pressure p(u), as a function of surface separation u, took the form where u o is a constant. Persson [75] commented that earlier theories, based on elastic contact, predicted that p(u)~u −a exp(−bu 2 ), where a = 1 for the Bush, Gibson and Thomas model [19], and a = 5/2 for the Greenwood and Williamson model [17]. Persson [75] also stated that the functional variation of pressure with the separation of Equation (34) was in better agreement with experimental data than the variation of pressure with separation from the Greenwood and Williamson [17] or Bush, Gibson and Thomas [19] models. It should be mentioned that Persson's models [22,23,75] are not thought to be widely used by the average tribologist for the calculation of mixed/boundary friction, as they are generally seen as being difficult to apply, and no simple equations are available for the prediction of the variation of the real contact area and load carried as a function of surface separation (such a relationship is most likely the reason why the Greenwood and Williamson [17] and Greenwood and Tripp [21] models are still so widely used).
A useful challenge was proposed by the editors of Tribology Letters to researchers to use a variety of models/numerical analyses to analyse a specific rough contact [115]. The challenge was proposed in 2015, and the results were published in Tribology Letters in a special issue in 2017 [116][117][118][119][120][121][122][123]. In total, 12 groups of researchers submitted 13 different approaches to the challenge, including one experimental approach (in which the researchers 3D printed a larger version of the specific contact and carried out experimental measurements on it). In a paper summarising the results of the challenge, Müser et al. [120] concluded that many of the methods used in the challenge identified the correct trends (such as how the real contact area depended on the load, for example, and how the mean gap between surfaces depended on the pressure), although it was pointed out that asperitybased models were expected to become less reliable when the roughness extends to more than 2.5 decades, and in addition, the predicted contact-patch-size distribution was not correct for such models (as bearing-area approaches neglect the effect of long-range elastic deformation). Some of the numerical approaches used would only take a few minutes to run for a 1000 × 1000 grid on a modern laptop.

The Impact of Lubricant Additives
Many of the papers referred to above simply refer to the fluid that separates the rough surface as having a viscosity, which is used to calculate the oil film thickness between the surfaces, and hence the λ ratio. This approach would suggest that all fluids should have friction curves of the same shape, when plotted against λ. In reality, commercially available lubricants contain numerous lubricant additives [124], some of which have a major effect on friction. The two major classes of lubricant additives that significantly affect friction are (i) anti-wear additives, of which zinc dialkyl dithiosphosphate (commonly abbreviated as ZDDP) is the best-known example, and (ii) friction modifiers, of which MoS 2 , graphite, glycerol mono-oleate, and oleyl amide are well-known examples.
Anti-wear additives, such as ZDDP, are thought to operate by forming a thick, highshear tribo-film. In the case of ZDDP, typical tribo-film thicknesses are of the order of 100-200 nm. A thorough review of ZDDP can be found in reference [125]. Recently, the research group at Imperial College reported [126] that the surface roughnesses of ZDDP films were significantly higher than those of the metal surfaces they were formed on (the work was performed using a Mini Traction Machine in which the ball and disk had very low surface roughnesses).
On the other hand, friction modifier additives are thought to form molecularly thin films (only a few molecular layers) which are easily shearable (i.e., they have a low shear stress). The reason for such behaviour is that friction modifiers tend to form surface layers in which there is only weak bonding between each layer. If such films form on the top of asperities, when sliding occurs, the layers of the friction modifiers will shear much more easily than the metal surfaces would if such tribo-films were not present. A review of friction modifiers can be found in reference [127]. Some research has been carried out to adapt rough-contact models to allow for such additives [80][81][82][83][84][85]. Most often, the tribo-film is modelled as a layer with different mechanical and roughness properties from that of the underlying substrate. Data on the mechanical properties of ZDDP tribo-films are available (for example, see reference [128]). If a Greenwood-Williamson-type model [17] is applied to a thick tribo-layer, then the shape of the friction versus λ curve would be influenced entirely by the roughness of the tribolayer (as the pre-factor in these models, which depends on the roughness and mechanical properties of the layer, will not affect the shape of the curve). When the correct tribo-layer roughness is used in the calculation of λ, Dawczyk et al. [126] recently reported that the friction versus λ curve is essentially the same as that for lubricants without additives (although the value of friction at λ = 0 may be higher, with a typical friction coefficient for ZDDP containing lubricants being 0.12 when λ = 0).
The explanation by Dawczyk et al. [126] cannot account for the change in shape of the friction versus λ curve when friction modifiers are used, as these do not form thick tribo-films and so cannot affect the value of λ. As mentioned above, the molecular structure of friction modifiers is such that they are sheared very easily. Measurements of the shear stress of friction modifier films have been reported by numerous researchers [129][130][131], both for organic films and for those formed from molybdenum, and the typical values reported are of the order of 25 MPa, which should be compared with typical shear stresses on ZDDP anti-wear films, which are thought to be in the range 200-300 MPa [132]. Understanding how friction modifiers influence the shape of the friction versus λ curve is still an active area of research.
It should also be mentioned that there is an increasing amount of published research on the topic of nanoparticles as alternative lubricant additives [133][134][135]. It would certainly be of interest to understand how such nanoparticles form tribo-layers, and whether or not their friction is affected in the same way as that of conventional anti-wear or friction modifier additives. However, the author is not aware of any commercially available lubricants in current large-scale production that contain such additives, and a review of the behaviour of nanoparticle tribo-films would be more timely once such additives are in more widespread use.

Experimental Data on Rough-Surface Friction
As discussed in Section 2, the proportion of mixed/boundary friction, X, can be considered to be the load carried by the asperities divided by the total load. In recent published work from Imperial College [126], the fraction of mixed/boundary lubrication, X, versus the lambda ratio, λ, was reported. In that work, the friction coefficient of a range of lubricants was measured using a Mini-Traction Machine (MTM). If f o was the friction coefficient measured at λ = 0, and f EHD was the friction coefficient measured at large values of λ, then the value of X was defined to be (f − f EHD )/(f o − f EHD ), where f was the measured friction coefficient at a particular value of λ. Figure 1 shows measured data from reference [126]. It was pointed out that if the correct value of surface roughness was used to calculate λ, then all of the lubricants would be found to fit onto a universal curve. It should be pointed out that although the RMS roughness of the MTM ball and disk were only about 4 nm each, when lubricants containing ZDDP anti-wear additives were used, relatively thick tribo-films were formed (between 100-200 nm), and the RMS surface roughness at the top of the tribo-films was significantly higher, at around 40 nm. λ was estimated by calculating the elastohydrodynamic oil film thickness in the contact (the radius of the curvature of the ball, the elastic moduli, the loads, speeds, temperatures were known, and the viscosity and pressure-viscosity coefficient of the lubricant were calculated). The authors of this paper checked lubricant friction data from other published MTM experiments [124,136], and found that base oils (lubricants that do not contain any additives) also fit well on the above curve (although in the case of base oils, the value of λ would be calculated using the original MTM ball and disk roughness values). The authors of this paper checked lubricant friction data from other published MTM experiments [124,136], and found that base oils (lubricants that do not contain any additives) also fit well on the above curve (although in the case of base oils, the value of λ would be calculated using the original MTM ball and disk roughness values).
Separately, other experimental work that can be used to calculate X versus λ has been reported by He et al. [137] and Cui et al. [138] on experimental test rigs that are significantly different from the MTM (in both the operating conditions and surface roughnesses of the test specimens). Figure 2 shows a fit through the data of Figure 1, with experimental data from references [137,138] overlaid. The fit through the data used the equation below [139]: where the values of the adjustable parameters k and a were determined using Excel's Solver function, and were found to be approximately k ≈ 3/2 and a ≈ 4/3 [139].
Lubricants 2022, 10, x FOR PEER REVIEW 15 of 23 Figure 2. Graph comparing the measured value of X (the proportion of mixed/boundary friction) versus λ for data from Imperial College [126], from He et al. [137], and from Cui et al. [138].

Application of Rough-Surface Models to the Prediction of Mixed/Boundary Friction
In Section 2, it was reported that if the proportion of mixed/boundary lubrication, X, is defined to be the load carried by the asperities divided by the total load; as such, it is of great practical interest, for the calculation of mixed/boundary friction, to know how X varies with λ.
Olver and Spikes [44] proposed the relationship below: In the previous section, a good fit to experimental data was found with the expression Figure 2. Graph comparing the measured value of X (the proportion of mixed/boundary friction) versus λ for data from Imperial College [126], from He et al. [137], and from Cui et al. [138].

Application of Rough-Surface Models to the Prediction of Mixed/Boundary Friction
In Section 2, it was reported that if the proportion of mixed/boundary lubrication, X, is defined to be the load carried by the asperities divided by the total load; as such, it is of great practical interest, for the calculation of mixed/boundary friction, to know how X varies with λ.
Olver and Spikes [44] proposed the relationship below: In the previous section, a good fit to experimental data was found with the expression below: with k ≈ 3/2 and a ≈ 4/3 [139]. Coy [155] proposed a linear relationship for the calculation of mixed/boundary friction: It was assumed by Coy [155] that X was zero for λ > 3.
The various asperity models reviewed in Section 4 can be used to find the variation of X with λ. For example, in the Greenwood and Williamson [17] and Greenwood and Tripp [21] models, the load carried by the asperities, W(λ), was calculated. If it is assumed that if X is equal to W(λ)/W(0), then the variation of X with λ can be found straightforwardly.
For the Greenwood and Williamson model [17], with an exponential probability distribution of asperity peaks, it is found that In the above equation, a factor of c was introduced, as it is possible that the RMS roughness of surface peaks may not be the same as the RMS roughness of the overall surface. However, in the absence of any further information, most tribologists would assume that c = 1.
For the Greenwood and Tripp model [21], it is found that Similarly, for the Bush, Gibson and Thomas model [19], we find that  [22,23] were not included in the comparison, as it was not clear to the author how to calculate X versus λ from those models. Given that Equation (37) was derived from a good fit to experimental data, it is clear that the contact asperity models which give the closest agreement with experimental data are (1) the Greenwood and Williamson model [17] when the probability distribution function is assumed to be an exponential, and (2) the Bush, Gibson and Thomas model [19]. Simply looking at the amount of mixed/boundary friction at λ = 1, it can be seen that the fit to experimental data gives X ≈ 40%, whilst the Greenwood-Tripp model [21], although it is still widely used today, only gives an estimate for X of around 13%. Although X has been defined in this paper to be equal to the load carried by the asperities divided by the total load on the contact, the majority of rough-surface contact models also calculate the real area of contact. Numerical calculations (by the author of this review) have found that X can also be equated to the real contact area (divided by the real contact area when λ = 0).
It should be mentioned that many researchers, in the calculation of the value of λ, estimate the oil film thickness in the contact by assuming that the surfaces are smooth, and then calculating the oil film thickness using standard hydrodynamic or elastohydrodynamic equations. However, it is known that, for rough surfaces, the standard Reynolds' equation should be modified using flow factors [156,157] to account for surface roughness. Therefore, ideally, oil film thicknesses should be calculated using a modified Reynolds' equation that includes flow factors, but this does not normally take place.

Discussion
Early models developed to explain Amontons' law tended to be analytical, based on plastic [15] or (more usually after 1966) elastic asperity deformation [17,19,21,30]. There has been a trend in the last 10-20 years to move to complete numerical solutions based on a variety of solution techniques [120] (Fast Fourier Transform (FFT), finite elements, the biconjugate gradient stabilized method, boundary element methods, and various types of molecular dynamics simulation). These methods can also be used on realistically sized grids. Most surface roughness profilometers will supply measured data on grids of size 1000 × 1000, and numerical solvers can solve such rough-surface contact problems in just a few minutes on a modern laptop. Despite the impressive progress reported in numerical solutions, assumptions still need to be made in the numerical models as to how the asperities deform. Many numerical models of rough-surface contact assume the elastic deformation of asperities (see for example references [105,107]) whereas other models allow for the elastic and plastic deformation of asperities (which will clearly change the surface shape and roughness parameters as contact progresses) [109,110]. A useful review of elastic-plastic contact mechanics may be found in reference [111]. Some numerical models also allow for adhesive effects [117,123]. It should be mentioned, in this context, that a consensus view seems to be that, for typical machine elements, most plastic deformation occurs during the "running-in stage", and in an insightful paper, researchers from Heriot-Watt University [158] combined numerical simulations of roughsurface contacts during "running-in" that used an elastic-plastic model, with experimental data, and found evidence for plastic deformation during "running-in", but that the volume of plastically deformed materially reduced as the "running-in" progressed. It was also found that the plasticity index (defined in Equation (32)) decreased as the "running-in" progressed, from an initial "plastic" value to one in the "elastic" regime once "running-in" had finished. Therefore, the assumption that asperities deform elastically may be quite reasonable for fully "run-in" surfaces, as originally suggested by Archard [30], Greenwood and Williamson [17] and Whitehouse et al. [18].
The inclusion of plastic contact in numerical rough-surface contact simulations is clearly important. However, the authors are not aware of many rough-surface contact simulations that allow for wear particles to be produced from such contact which can then act as hard contaminants (a process that clearly takes place in realistic contacts). Two papers by Zhang, Moslehy and Rice [159,160] did attempt such an analysis, as reported in 1991, but this work does not seem to have been followed up.
It is known that repeated plastic contacts can result in changes in hardness and material composition, and ideally these effects also need to be included [161][162][163]. The ambition of many current leading tribologists is the development of a "tribology digital twin" in which all of the various effects are accounted for, such that a realistic tribology experiment could essentially be reproduced on a computer. At present, the development of a "digital twin" in which all of these effects are accounted for is still likely many years away, although recently there has been increasing use of artificial intelligence (AI) and machine learning (ML) methods within the tribology community (two recent useful reviews can be found in [164,165]). It should be mentioned that physics-based insights into the transition from mild to severe wear have recently been published by Aghababaei et al. [166], based on atomistic simulations, which could help in the development of tribological digital twins. Another study of a similar type was reported by Han and Zou [167], who studied the evolution of contact characteristics during a scuffing process, using contact mechanics and plasto-elastohydrodynamic lubrication models. In parallel with these simulation studies, there are also experimental studies that provide evidence of surfaces "roughening up" just prior to failure, with examples of such studies being those reported Roper and Bell [168] and by Bell and Willemse [169].
Despite the impressive advances in numerical and analytical rough-surface contact models over the last 50-70 years, there is still a need for "quick and simple" methods by which tribologists and engineers can make reasonable estimates for the proportion of mixed/boundary friction, and the friction losses that arise in this regime. Therefore, simple models, in which the proportion of mixed/boundary friction can be estimated as a function of a simple parameter, such as the lambda ratio, as described in the previous section, are still important. If such simple models can be validated, and/or improved by comparison with detailed analytical/numerical or experimental data, this would be even better.

Conclusions
The study of rough-surface contact has a long history, with the modern understanding of such contacts beginning around the 1950s, when it was generally assumed that contacts behaved plastically, with a realization in the late 1950s and 1960s that-once surfaces had "run-in"-many contacts behaved elastically. In both plastic and elastic contact conditions, Amontons' law of friction (that friction was proportional to load) could be derived approximately, although in the case of elastic contact conditions, this only became apparent when a distribution of asperity heights and shapes was included in the models. Many relatively simple models were developed, and are still widely used today by tribologists and engineers. Perhaps the best known of these models are those of Greenwood and Williamson [17], Greenwood and Tripp [21], and Bush, Gibson and Thomas [19], all of which assume relatively low loads and low real contact areas (such that the assumption of elastic contact is valid). Later, Persson [22,23] developed models which were more accurate for high real areas of contacts (such as those that occur in rubber contacts). However, Persson's models are not thought to be that widely used within the tribology community due to their relative complexity and the lack of simple equations that can be easily used in engineering calculations. Recently, there has been a trend to treat rough-surface contacts with full numerical simulations, using measured surface profilometer data (usually generated on a 1000 × 1000 sized grid); alternatively, typical surface profiles can be generated artificially [170]. Such models are not generally available to typical tribologists and engineers, such that there is still a need for simple methods that give "rough and ready" answers to how much friction occurs in mixed/boundary lubricated contacts. Experimental data are becoming available [126,137,138] in which the proportion of mixed/boundary friction has been measured as a function of the lambda ratio. The availability of such data enables a comparison to be made of the various rough-surface contact models. Experimental data from [126] found that the proportion of mixed/boundary friction, X, was between 30 and 50% when λ = 1. The Bush, Gibson and Thomas model [19] predicts a value for X of about 32% when λ = 1, whilst the Greenwood and Williamson model [17], with an exponential distribution of asperity heights, predicts that X ≈ 37% when λ = 1. On the other hand, the Greenwood and Tripp model [21], which is still widely used today for the prediction of mixed/boundary friction, predicts a value of X of only around 13% at λ = 1, suggesting that mixed/boundary friction is significantly underestimated with that model, which is a conclusion in agreement with the recent work of Leighton et al. [171].
Looking to the future, as energy efficiency becomes a matter of greater and greater importance, there is likely to be more mixed and boundary lubrication in machines, as lubricant viscosities are tending to become lower, operating loads and temperatures are tending to become higher, and machines and the amount of lubricant within them are tending to become smaller. The validation of rough-surface contact models with experimental data will then become a matter of greater importance. The recent use of numerical simulations for the prediction of the rough-surface real area of contact and friction would then be useful in the generation of data to which simpler equations could be fitted, which could be used by the wider tribological community to better predict mixed/boundary friction and wear.
It is anticipated that research in rough-surface contacts will continue for many years to come. Data Availability Statement: The data in this paper can be found in the various references listed.