## Abstract

We analyse contact between a rigid cylindrical indenter and an elastic–perfectly plastic solid with a fractal surface roughness, which is idealized as a Weierstrass profile. A hierarchical scheme is used to calculate quantities such as the number of contact spots; their size; the true contact pressure acting on the contact spots, and the total true area of contact. In addition, we calculate critical conditions that will initiate plastic flow in the asperities, and determine the extent of this plastic flow as a function of the applied loading, the material properties of the deforming solid, and the surface roughness. An important conclusion of our study is that a perfectly fractal description of surface roughness appears to lead to unphysical predictions of the true contact size and number of contact spots, for both elastic and elastic–plastic solids. Possible approaches to resolving these difficulties are discussed.

## 1. Introduction

When two rough surfaces are pressed into contact, they meet only at the highest points on the two surfaces. To predict the size of the resulting contact spots, the true contact pressure, and the true area of contact is a venerable problem in contact mechanics, and is clearly of considerable practical importance, since phenomena such as friction, wear, and contact fatigue are controlled by the nature of asperity contacts. Nevertheless, a fully satisfactory solution remains elusive.

Difficulties in modelling contact between rough surfaces arise because of the fractal nature of surface roughness. Early models, such as the classic work of Greenwood & Williamson (1966) and subsequent refinements by Whitehouse & Archard (1970), Nayak (1971), Bush *et al*. (1975), Greenwood (1984) and McCool (1986) among others, idealized a rough surface as a collection of asperities with spherical tips of equal radius, but with a Gaussian distribution of heights. These models predict a well-defined number of total contact spots, contact size, contact pressure and contact area, which depend on the curvature of the asperity tips, the asperity height distribution, and the elastic and plastic properties of the two contacting surfaces. The Greenwood–Williamson model has proved difficult to apply to real surfaces, however. Surface profiles are measured by sampling the surface at discrete intervals. Asperity tips are identified by finding points that are higher than their neighbours. The curvature of each asperity tip is to be determined by fitting a curve (in two-dimensions) or surface (in three-dimensions) through the points surrounding the peak. For a differentiable surface, this procedure converges to identify a well-defined set of asperities as the sampling interval is reduced. When this procedure is applied to experimentally measured surfaces, however, the number of asperities and their curvature appear to increase without limit as the sampling interval is reduced to extremely small values.

This is, of course, a consequence of the fractal nature of surface roughness. There has been some debate as to the most accurate description of realistic surfaces, (Whitehouse & Archard 1970; Sayles & Thomas 1978; Majumdar & Bhushan 1990; Majumdar & Tien 1990, 1991) but all models agree that surfaces have an approximately self-affine fractal character over a wide range of length-scales, and usually continue to exhibit fractal behaviour down to the limits of the resolution of standard surface profilometres. Indeed, Majumdar and Bhushan (1990, 1991) show data that suggests that roughness spectra for hard disk drives remain fractal down to wavelengths of a few hundred nanometres. For a perfectly fractal surface, Greenwood and Williamson's description of the surface topography breaks down completely, since such a surface does not contain well-defined summits.

Many recent studies of rough surface contact have, therefore, abandoned the Greenwood–Williamson description of surface geometry in favour of one that captures in some way the fractal character of surface roughness. Various approaches have been proposed, including procedures for adapting Greenwood and Williamson's calculations to fractal surfaces by introducing generalized definitions of asperities (Majumdar & Bhushan 1990, 1991), direct numerical simulations of experimentally measured profiles (e.g. Sayles 1996; Borri-Brunetto *et al*. 1999; Polonsky & Keer 2000*a*,*b*; Hyun *et al*. 2004; Pei *et al*. in press), or hybrid schemes which combine numerical simulations with Greenwood and Williamson's approach (see Mihailidis *et al*. 2001 for a review). Attempts to adapt Greenwood–Williamson to a fractal surface appear to suffer from the same difficulties as the original Greenwood–Williamson analysis: there is no unambiguous way to characterize asperities on a fractal surface. Numerical simulations provide considerable insight into the deformation fields associated with rough surface contact, but there is a limit to the spatial resolution of any numerical scheme, which artificially truncates the roughness spectrum. This again introduces problems very similar to those associated with the Greenwood–Williamson model: for a fractal surface, predictions of the contact size, number of contacts, and true contact area will not converge as the resolution of the numerical scheme is refined. Great care must be taken to ensure that this truncation does not lead to spurious predictions (Borri-Brunetto *et al*. 1999).

A significant step towards a satisfactory model of the elastic contact between rough fractal surfaces was made recently by Ciavarella *et al*. (2000) and by recent work of Persson (2001, 2005). The approximate approach of Ciavarella *et al*. (2000) will be followed closely in our work and so will be described in more detail in §§2 and 3 of this paper. Briefly, they chose to approximate surface roughness using the two-dimensional Weierstrass function (given in equation 2.1 below). The Weierstrass function consists of a series of discrete harmonic functions, with progressively decreasing wavelength and amplitude. Ciavarella *et al*. devised a hierarchical scheme to predict the apparent contact pressure and area on each successive scale of roughness. For the limit of a perfectly fractal surface, prediction that the true contact area consists of an infinite number of contact spots, with zero size, infinite pressure, and vanishing total area. Interestingly, they also showed that the same conclusion can be derived from Archard's (1957) early models of rough surface contact (Ciavarella *et al*. 2000; Ciavarella & Demelio 2001). This is a purely formal mathematical result, of course: in practice phenomena such as plastic flow or adhesion might alter these predictions, and in any event no real surface is perfectly fractal. There are some applications, such as modelling friction, where the true contact size and area are of great interest, and predicting these quantities remains an unresolved problem. For modelling the compliance of the surface, or its electrical contact resistance, and possibly physical phenomena such as plastic flow, wear, and the initiation and growth of contact fatigue cracks, the behaviour of very short-wavelength roughness is less important, since damage events occur at sub-micron length-scales rather than nanometre length-scales. The main advantage of the approach devised by Ciavarella *et al*. is that the behaviour of long-wavelength roughness scales can be predicted in terms of measurable surface parameters without needing to identify real asperities on the surface.

The work of Ciavarella *et al*. (2000) was restricted to elastic solids. In this paper, our goal is to extend their calculations to account for plastic deformation. The problem to be solved is illustrated in figure 1. We consider an isotropic, elastic–perfectly plastic solid with Young's modulus *E*, Poisson's ratio *ν* and yield stress *σ*_{Y}. The solid has a rough surface, which is characterized by the two-dimensional Weierstrass profile given in equation (2.1). The solid is indented by a smooth cylinder with radius *R,* which is subjected to a force *P* per unit length out of the plane. The goal of our analysis is to estimate quantities such as the number of contact spots, their size (or more precisely, their size distribution); the true contact pressure acting on the contact spots, and the total true area of contact. It is of particular interest to calculate critical conditions that will initiate plastic flow in the asperities, and to determine the extent of this plastic flow as a function of the applied loading, the material properties of the deforming solid, and the surface roughness.

The remainder of this paper is organized as follows. In §2, we give a more careful statement of the problem to be solved, the assumptions underlying our calculations, and the results that will be calculated. Section 3 gives a brief review of the indentation response of an elastic–perfectly plastic solid with sinusoidal roughness, as the first step towards analysing the Weierstrass profile. In §4, we review the procedure proposed by Ciavarella *et al*. to analyse contact of the Weierstrass profile, and devise a simplified, approximate procedure that can be used to derive analytical results for both elastic and plastic contacts. Section 5 presents detailed results and §6 contains conclusions.

## 2. Rough surface model

Majumdar & Bhushan (1991) suggest that the two-dimensional Weierstrass function(2.1)with *γ*≈1.5 and *D*≈1.4 provides a good qualitative description of many practical surface profiles. The power spectrum for the ideal Weierstrass profile is a series of Dirac delta functions, but Berry & Lewis (1980) argue that the Weierstrass distribution can be used to approximate a continuous power spectrum by ensuring that the integral of the continuous power spectral density matches that of the Weierstrass function. The resulting continuous approximation to the power spectral density is given by(2.2)where the circular frequency *ω* has dimension of L^{−1}. Although, the Weierstrass function captures the multi-scale nature of surface roughness, it is not a completely realistic description of real surfaces. The actual power spectrum for measured surfaces is continuous, not a series of Dirac delta functions. In addition, roughness wavelengths in the range 50 μm–1 mm are strongly sensitive to processing and in general are not well characterized by even the continuous approximation to the Weierstrass function. Similarly, the Weierstrass profile becomes unphysical at very short wavelengths: as *m*→∞, the surface becomes needle-like, with roughness amplitude much greater than its wavelength. Any departure from ideal behaviour at long wavelength (low *m*) will lead to quantitative, but not qualitative changes in our predictions and could readily be incorporated in the computations. The departure from ideal behaviour at short wavelengths (large *m*) has a critical influence on the true size of the contact spots, and for ideally elastic surfaces, on the total true area of contact and contact pressure. It will not, however, have a significant influence on the behaviour of longer-wavelength scales of roughness.

Analysing contact of a Weierstrass surface with a realistic value of *γ* is difficult, because roughness scales with successive values of *m* do not clearly decouple. Ciavarella *et al*. (2000) show that calculations are greatly simplified by approximating a fractal surface with a larger value of *γ*>5. The surface then consists of a series of harmonic profiles, with widely separated wavelength, each of which can be analysed independently. Although, the resulting surface does not resemble any actual surface, it does capture the self-affine fractal character of real surfaces, and so provides an invaluable qualitative picture of the behaviour of two contacting fractal surfaces.

We now consider an elastic–perfectly plastic solid, with Young's modulus *E*, Poisson's ratio *ν* and yield stress *σ*_{Y}, which has a surface roughness defined by equation (2.1). The solid is indented by a smooth cylinder with radius *R* as illustrated in figure 1. The cylinder is loaded by a total force (per unit out of plane length) *P*, to induce a contact pressure distribution *p*(*x*). Of course, the two surfaces actually meet only at the tips of asperities. We shall find that for an elastic–perfectly plastic solid with an ideal Weierstrass profile, the true contact consists of an infinite number of contact spots of zero size, which are all subjected to the same, uniform contact pressure. This is a purely formal result, however, and does not characterize the loading or deformation of the solid in a useful manner. We obtain a better picture by regarding the loading and deformation to consist of a series of separate scales, , as illustrated in figure 1. The lowest scale, *n*=−1 corresponds to the nominal contact. Successive terms characterize the behaviour of the solid at length-scales comparable to successive wavelengths in the Weierstrass profile. The fields associated with each term interact weakly through the requirement that the pressure at scale *n* must be statically equivalent to the pressure at scale *n*+1.

The nominal contact area and pressure (scale *n*=−1) can be computed using results for smooth contacts (see e.g. Johnson 1985). It is helpful to review these results briefly. The severity of the nominal load is conveniently characterized by a dimensionless load parameter . For , the contact is elastic, and the contact width and nominal contact pressure distribution follow as(2.3)where is the average nominal contact pressure. For , the contact is fully plastic, and the contact width and stress distribution can be approximated using the rigid plastic slip-line field solution, which gives(2.4)where is the ratio of the plane strain hardness of the solid to its tensile yield stress. In the intermediate regime , the contact width and pressure distribution would need to be computed using numerical simulations. This will not be attempted here: we shall focus instead on the limiting cases of perfectly elastic and fully plastic bulk contacts.

It is important to note that the nominal contact geometry introduces a natural and well-defined long-wavelength cut-off for the roughness spectrum. Roughness wavelengths greater than *a*_{−1} clearly do not influence the contact pressure. The roughness component with wavelength closest to *a*_{−1} should be regarded as part of the nominal contact geometry (its effect is equivalent to modifying the shape of the indenter). The first component of the profile to influence the contact pressure distribution has index , and has wavelength and amplitude:(2.5)Henceforth, we will characterize the surface roughness within the area of contact as(2.6)where *g*_{0} and *λ*_{0} represent the amplitude and wavelength of the first roughness component that has wavelength less than the nominal contact width. It will be important to keep in mind that while *γ*, *D*, *G*_{0} and *Λ*_{0} are properties of the surface geometry only, *g*_{0} and *λ*_{0} depend also on the nominal contact geometry as shown in equation (2.5). We shall return to this issue when presenting our results in §5.

If we now examine a small region of the nominal contact, with length and centred at some point *x* within the contact area, we observe a sinusoidal surface with wavelength *λ*_{0} in contact with a nominally flat indenter. The indenter is subjected to a uniform nominal contact pressure , which induces a contact width 2*a*_{0} and pressure *p*_{0}(*x*) on the zeroth roughness scale. We can then examine more closely a region with size within one of these contacts. Here, we observe a sinusoidal surface with wavelength *λ*_{1}, subjected to nominal contact pressure , inducing a contact size *a*_{1}, and contact pressure *p*_{1}(*x*). This process can evidently be continued indefinitely.

The goal of our analysis is to determine, for each roughness scale *n*, the total number of contact spots *N*_{n}, the apparent size of the contact spots *a*_{n}, and the apparent contact pressure distribution *p*_{n}(*x*). Calculations are complicated by the fact that the quantities of interest vary with position on the surface. It is possible (but exceedingly cumbersome) to calculate this spatial variation. In practice, it is preferable to quantify the variation by computing distribution functions, which, for example, specify the fraction of asperity contacts with size between *a*_{n} and *a*_{n}+d*a*_{n}, or the fraction of the nominal contact area that is subjected to pressure between *p*_{n} and *p*_{n}+d*p*_{n}, for each scale .

The distribution function for the pressure plays a particularly important role in subsequent calculations. For each scale of roughness, we define a dimensionless function such that is the fraction of the nominal contact (which has total width 2*a*_{−1}) that is subjected to a dimensionless pressure between *p*_{n}/σ_{Y} and . As an example, the distribution function characterizing the nominal contact can easily be computed. The cumulative distribution *Q*(*p*) (the proportion of the contact over which the pressure exceeds *p*) can be computed from the contact pressure distribution as . The distribution function then follows as *q*_{−1}(*p*)=−∂*Q*_{−1}/*p*. For example, if the nominal contact is elastic, we have(2.7)The distribution function follows as(2.8)where *ξ*=*p*/*σ*_{Y} . Similarly, for fully plastic contact, the distribution function is evidently(2.9)where *δ* denotes the Dirac delta distribution and H_{0}≈2.9 is the ratio of the hardness to the tensile yield strength of the solid, as before. Similar distributions will characterize the pressure distribution for each higher roughness scale. Although, it is not possible to write down a simple expression relating the contact pressure *p*_{n}(*x*) to *q*_{n}(*p*) for *n*≥0, Ciavarella *et al*. (2000) show that this is not necessary. Their procedure for computing *q*_{n} will be reviewed in detail in §4.

The distribution function *q*_{n}(*p*) provides a convenient description of the contact pressure acting on the *n*th roughness scale, but is difficult to interpret physically. To provide a more direct measure of the deformation of each roughness scale, we introduce the contact size distribution function *S*_{n}(*a*) defined such that(2.10)with *N*_{n} the total number of contact spots at the *n*th scale, gives the number of contacts with size between *a* and *a*+d*a*. Similarly, as a more direct measure of the loading applied to the contacts in the *n*th scale, we define the mean contact pressure for an individual contact spot as(2.11)We then introduce the contact pressure distribution *P*_{n}(*p*^{mean}) such that(2.12)gives the number of contacts in the *n*th scale that are subjected to mean pressure between *p*^{mean} and *p*^{mean}+d*p*^{mean}. Naturally, since there is a direct relationship between the contact size and mean contact pressure for individual contact spots, there will be a similar direct relationship between *S*_{n}(*a*) and *P*_{n}(*p*^{mean}).

Although, we shall formally compute expressions for *S*_{n}(*a*) and *P*_{n}(*p*^{mean}) for our model surface, we readily concede that in view of the rather crude description of surface topography that forms the basis for our calculations, the exact distribution functions may have little practical significance. In many cases, the *average* values of contact size and pressure for each scale of contact are sufficient. These are formally related to *S*_{n}(*a*) and *P*_{n}(*p*^{mean}) through(2.13)In practice, however, it is simpler to compute and directly, rather than via the distribution functions. The total contact area and nominal contact pressure:(2.14)

We shall calculate these quantities using the elegant procedure developed by Ciavarella *et al*. (2000) to model contact between elastic solids with a Weierstrass profile. A key assumption in their work is that the behaviour of a particular roughness scale *n* at a point on the surface can be calculated from the solution to a sinusoidal surface indented by a flat, frictionless punch. Successive roughness scales interact only through the requirement that the nominal pressure that acts on scale *n* is equal to the true pressure for the preceding scale *n*−1. This assumes that no interaction occurs between the deformation fields associated with successive roughness scales.

For an elastic solid, this assumption can easily be justified using Saint Venant's principle. For elastic–plastic materials, no rigorous proof exists. We speculate, however, that provided the bulk contact (at scale *n*−1) consists of a single, isolated contact area, with a size much smaller than any characteristic structural dimension, the same decoupling can be applied for elastic–plastic solids. As an extreme example, consider the limiting case of a rigid-perfectly plastic solid. The slip-line field solutions for contact between a semi-infinite solid and a single indenter invariably show a rigid region under the contact area (Hill 1950; also reproduced in Johnson 1985). Consequently, the deformation associated with the first roughness scale in a Weierstrass surface (scale 0) may be embedded within the rigid region associated with the nominal contact. The behaviour of the first scale of roughness is, therefore, equivalent to that of a sinusoidal surface, with rigid boundary conditions at infinity. Our numerical solution to this problem shows that a rigid region forms under each contact area. Consequently, the second roughness scale can similarly be embedded within this region, and this process can be continued for all higher roughness scales. Since, the deformation field constructed in this manner forms a kinematically admissible collapse mechanism, our solution should strictly be regarded as an upper bound to the collapse load, and a lower bound to the contact size and contact area. We present some further support for our calculations in electronic supplementary material, where the analytical calculations are compared with the predictions of full-field finite element computations for up to three roughness scales. It is important to recognize, however, that if the deformation field for the bulk contact (scale −1) is not contained, a rigid region may not form below the bulk contact in the fully plastic limit. In this case, the first (longest-wavelength) scale of roughness is not contained, and so cannot be idealized as a sinusoidal surface with rigid boundary conditions at infinity. Similarly, there is no reason to suppose that rigid regions would develop to contain shorter-wavelength scales. Such unconstrained flow is particularly likely to occur in metal forming operations, and existing slip-line field solutions suggest that asperities are easily flattened in the presence of unconstrained far-field deformation (Sutcliffe 1988).

## 3. Indentation response of a sinusoidal surface

Clearly, establishing the response of a simple sinusoidal surface to indentation loading is a key step towards solving the contact problem outlined in §1. The indentation response of an elastic–perfectly plastic solid with a harmonic roughness was analysed in detail by Gao *et al*. (in press). Only a few of their results are needed here: these are reviewed briefly below.

We focus attention on a single roughness scale with wavelength and amplitude(3.1)as illustrated in figure 1. The substrate is an elastic–perfectly plastic solid with Young's modulus *E*, Poisson's ratio *ν* and yield stress *σ*_{Y}. Gao *et al*. (in press) show that the strength of the surface is conveniently characterized by a material and geometric parameter given by(3.2)with . In the discussion to follow, the value of *ψ*_{n} for the zeroth roughness scale (*n=*0) plays a central role in characterizing the resistance of the rough surface to plastic flow. In terms of this parameter, we have(3.3)

Note that the resistance to plastic flow varies in inverse proportion to *ψ*: a rough surface with a low yield stress has a high value of *ψ*, while a smooth, hard surface would have a low *ψ*. In addition, note that *ψ*_{n} increases with *n*, so we anticipate that roughness scales with long-wavelength will be elastic, but all scales beyond a critical value of *n* will deform plastically.

The solid is indented by a flat, rigid surface, under the action of a nominal pressure. The loading induces a contact width 2*a*_{n}, pressure distribution *p*_{n}(*x*) and mean pressure given by the functional relationships(3.4)

Here, we have shown that formally that *A*, *Φ* and depend on both *ψ*_{n} and *g*_{n}/λ_{n}. In fact, if the surface remains elastic, exact calculations show that the solution depends only on *ψ*_{n} and is independent of *g*_{n}/λ_{n}. If the surface deforms plastically, numerical simulations show that the functions depend only weakly on *g*_{n}/λ_{n}. It would not be difficult to include this variation in our computations, but for simplicity we have chosen to neglect it.

Finally, we must also compute the pressure distribution function *q*_{n}(*p*/σ_{Y}), which is defined so that is the fraction of the sinusoidal surface that is subjected to normalized pressure between *p*_{n}/σ_{Y} and . This distribution function can be determined from the contact pressure, following the procedure given for the nominal contact in equations (2.7) and (2.8). The result will be specified by the functional relationship(3.5)

We have calculated approximations to the functions *A*, *Φ*, and *q* by fitting curves to the numerical simulations reported in Gao *et al*. (in press), which are reproduced in figure 2. The resulting expressions are somewhat lengthy, and are given in detail in the electronic supplementary material A. Briefly, we observe four general regimes of behaviour, depending on the normalized contact pressure .

For , the contact remains elastic. Expressions for the contact size and pressure distributions in the elastic regime are given in Westergaard (1939) and Ciavarella

*et al*. (2000), and are reproduced in electronic supplementary material.For , where

*H*_{0}≈2.9 is the ratio of the (plane strain) hardness to the yield stress of the solid, the contacts are in the elastic–plastic regime. In this regime, we estimate the contact size and pressure using a smooth interpolation between elastic and fully plastic behaviour. Full expressions are given in electronic supplementary material.For , the contacts behave like isolated (plane strain) hardness indentations, with a uniform contact pressure with magnitude .

For , the contact fraction satisfies

*A*>2/3, so that the plastic zones under neighbouring asperities begin to interact. In this regime, the contact pressure under each contact remains approximately uniform, but the magnitude of the pressure increases monotonically as the contact fraction approaches one. The limiting value of the contact pressure as*A*→1 plays an important role in subsequent discussions. The limiting value of is a function of*ψ*_{n}, and as approaches : approximately twice the material hardness.

Our approximate curves are compared with the results of finite element simulations in figure 2, which shows a parametric plot of the nominal and mean contact pressures as a function of the contact size as the load increases, for various values of *ψ*.

## 4. Contact mechanics for a general Weierstrass profile

We proceed to analyse contact of the Weierstrass profile. Our analysis follows exactly the work of Ciavarella *et al*. (2000) and so will be reviewed only briefly. We will first present a general procedure for computing distribution functions for contact size and pressure for each scale. Subsequently, we suggest a simple approximation that can be used to obtain closed-form estimates for the number of contact spots, their size, and mean pressure for each scale of loading.

### (a) Results for general surface contact

Although, it is difficult to compute the spatial contact pressure distribution for the Weierstrass surface, it is straightforward to determine the probability distribution function defined in §3. The calculation begins by computing the distribution function , which specifies the proportion of the nominal contact area that is subjected to a pressure between *p*/σ_{Y} and (*p*+d*p*)/σ_{Y}. The result is given in equation (2.8).

The distribution *q*_{0}(*p*) for the zeroth roughness scale may, therefore, be computed as follows. Recall that we introduced the dimensionless function *q*, such that specifies the proportion of the sinusoidal surface that is subjected to pressure between *p*/σ_{Y} and (*p*+d*p*)/σ_{Y} when loaded by a normalized nominal pressure *ξ*. The zeroth roughness scale is subjected to a distribution of nominal pressure, as quantified by . Summing over all possible values of the contact pressure for the nominal contact, we find that(4.1)The upper limit on the integral here is purely formal, of course: the nominal contact pressure does not exceed , so for . This result can be applied recursively, to see that(4.2)

Once the hierarchy of distribution functions has been computed, various quantities of interest can be deduced. For example, to calculate the total (apparent) area of contact for the *n*th scale, note that the area of contact for one wavelength of the *n*th scale subjected to pressure is . In addition, there are such wavelengths within the nominal area of contact. The total contact area, therefore, follows as(4.3)The total number of contact spots at the *n*th scale within the nominal contact is(4.4)whereupon the average size of the contacts at the *n*th scale follows as(4.5)The average contact pressure can be found by a similar approach, with the result(4.6)Finally, the nominal contact pressure for the *n*th scale follows as(4.7)

We can also find expressions for the contact size and mean pressure distribution functions *S*_{n}(*a*) and defined in equations (2.10) and (2.12). To compute the distribution function for the contact size, denote the pressure required to induce a contact size *a* is given by the functional relationship , where *A*^{−1} denotes the inverse of the function *A* defined in equation (3.4). Given the distribution function for the contact pressure at (*n*−1)th scale, we can easily compute the probability distribution function of the contact fraction for scale *n* as follows. For a single contact, , so that let . The probability of finding a contact fraction in between and is(4.8)

### (b) Simplified analysis—the uniform contact pressure assumption

The results listed in the preceding section can be simplified dramatically if the contact pressure distribution for each scale in the contact (including the nominal contact area) is *uniform*. This is an excellent approximation when each scale of contact is fully plastic. At lighter loads, the nominal contact, and some long-wavelength roughness scales, will remain elastic, and for these scales the distribution of contact pressure is not uniform. Surprisingly, however, we have found that averaged quantities such as the total area of contact, the average contact size, and the average contact pressure can be computed remarkably accurately even for elastic surfaces by replacing this non-uniform contact pressure distribution with a statically equivalent, uniform pressure.

Proceeding along these lines, therefore, we suppose that the entire nominal contact area is subjected to the same contact pressure. For this approximation, the distribution function , which specifies the fraction of the nominal contact area (scale −1) subjected to pressure between *p*/σ_{Y} and (*p*+d*p*)/σ_{Y}, must, therefore, be given by(4.9)where *δ*(*x*) denotes the Dirac delta distribution and . The distribution function for the zeroth scale of roughness follows from equation (4.2) as(4.10)

Recall now that the function specifies the fraction of a sinusoidal surface that is subjected to pressure between *p*/σ_{Y} and (*p*+d*p*)/σ_{Y}, when indented by a nominal contact pressure . If the pressure is approximated by a statically equivalent uniform distribution, with magnitude *p*^{mean}, it follows that(4.11)where the functions *A* and *Φ*^{mean} are defined in equations (3.4), and given specific forms in equations (A2)–(A18) of the electronic supplementary material. Repeating this procedure for subsequent roughness scales gives the following result for the distribution functions *q*_{n}(*p*):(4.12)Here, the mean contact pressure at each scale can be determined from(4.13)where the sequence is started with given in equation (4.9).

Given for each scale, and the distribution functions , we can determine the remaining quantities of interest using equations (4.3)–(4.7). The total area of contact at the *n*th scale follows as(4.14)where the function *A* is defined in equation (3.4). The number of contacts follows as(4.15)and the average contact size is given by(4.16)Finally, equation (4.7) provides an additional way to compute as(4.17)from which we also recover the expected relationship between the mean pressure at successive scales:(4.18)

### (c) Approximate analytical expressions for average contact size, pressure, number of contacts, and total contact area for elastic contact

If the contacts remain elastic, the procedure outlined in the preceding section can be used to derive simple analytical expressions for the average contact size, average contact pressure.

Under these conditions, we may approximate(4.19)whereupon equation (4.13) reduces to(4.20)This recursion relation is satisfied by(4.21)The average contact size may be estimated using equation (4.16) as(4.22)Similarly, the total area of contact at the *n*th scale follows as(4.23)The total number of contacts at the *n* th scale can be estimated as(4.24)For large *n*, these results may be approximated by their asymptotes(4.25)

(4.26)

(4.27)

(4.28)

As expected, equation (4.27) shows that the total contact area at the *n*th scale (for large *n*) is proportional to the contact load *P* and independent of the nominal contact size *a*_{−1} or the indenter geometry (other than through its effect of setting the wavelength of the zeroth roughness scale). This result also provides an opportunity to check the predictions obtained using the constant pressure approximation (§4*b*) with the exact solution. Ciavarella *et al*. (2000) find that , giving a coefficient that is within 4% of the approximate value in equation (4.27). A more detailed comparison of the approximate and exact calculations in §5 shows similar agreement under all conditions.

Note also that equation (4.25) shows that the mean contact pressure acting on each roughness scale becomes progressively more severe as *n* increases. Consequently, for an elastic–plastic solid, the short-wavelength roughness scales will always deform plastically. We also see that the formal limit *n*→∞ for a perfectly elastic solid consists of an infinite number of contact spots, with zero size, which are subjected to infinite contact pressure. The total true contact area vanishes.

### (d) Deformation and loading for the most heavily loaded region at each scale (elastic contact)

To calculate the conditions necessary to initiate yield at each successive scale of roughness, or to identify conditions that will just begin to bring a roughness scale into complete contact with the indenter, it is of interest to determine the contact size and contact pressure for the most heavily loaded asperities at each successive scale of roughness. We assume that all roughness scales, as well as the nominal contact, remain elastic; and assume also that . Equation (2.3) shows that the maximum contact pressure exerted by the nominal contact is(4.29)

The most severely loaded asperities on the zeroth roughness scale are, therefore, subjected to a nominal contact pressure . Equation (A3) in electronic supplementary material gives the maximum contact pressure for a sinusoidal surface that is subjected to a normalized nominal pressure as(4.30)whence the maximum contact pressure on the zeroth roughness scale is . Similarly, the maximum contact pressure on successive scales of roughness follows as . The recursion relation is satisfied by(4.31)Other than a numerical factor, this result is identical to equation (4.21), which estimates the average contact pressure on each roughness scale. However, equation (4.31) does not assume isolated Hertz-type contacts and so is valid for all values of *ψ*_{0}, provided that no scale reaches full contact.

### (e) A graphical approach to calculating contact size, pressure and area for roughness scales in the fully plastic regime

If the bulk contact is fully plastic (which requires ), and the surface roughness is such that , then the procedure outlined in §4*b* gives exact results, since all roughness scales are in the fully plastic regime, so the contact pressure distribution at each scale is uniform. In addition, under these conditions, the method can be illustrated using a simple graphical procedure, originally suggested by Childs (1973, 1977) in a discussion of the persistence of asperities under contact loading.

We begin by simplifying the behaviour of the sinusoidal surface for the fully plastic regime. Consider a single scale of roughness with wavelength *λ*, which is indented by a rigid platen. The surface is subjected to a nominal contact pressure (force per wavelength) , which induces a contact size *a*. Each contact spot is subjected to a uniform contact pressure with magnitude *p*^{mean}. For values of *ψ*_{0}>16, the relationship between the mean normalized contact pressure and contact size for the fully plastic regimes given in equations (A14)–(A18) in electronic supplementary material can be approximated by(4.32)where , and . In addition, the normalized nominal contact pressure follows as . Both these functions are plotted in figure 3.

Next, we apply the iterative scheme outlined in §4*b* for this simplified situation. Since the nominal contact (scale −1) is in the fully plastic regime, the nominal contact pressure is uniform, and has magnitude . This implies that the nominal pressure exerted on asperities in the zeroth roughness scale is . Figure 3 can be used to deduce the normalized contact size corresponding to this nominal contact pressure; whence the true contact pressure can also be found. The same procedure may be used to progress to higher scales of roughness: the condition determines the normalized contact size at each successive roughness scale, whereupon can be deduced. In addition, the apparent total contact area at each scale follows as(4.33)In particular, the total contact area in the limit *n*→∞ is(4.34)Thus, we find that the contact area is again proportional to the normal load, but instead of converging to zero as *n*→∞, the total contact area converges to a finite value.

A similar conclusion was reached by Willner (2004) using a more sophisticated series of numerical simulations. His analyses did not account for interactions between neighbouring contact spots, however, and consequently gave a limiting area of . More recently, Pei *et al*. (in press) have conducted an extensive study of the indentation response of an elastic–plastic solid with a nominally fractal surface roughness using three-dimensional finite element computations. They found that the true area of contact is a function of the ratio of yield stress to Young's modulus for the solid: this behaviour is most likely due to an artificial truncation of the roughness spectrum resulting from the resolution of the finite element mesh. For sufficiently low values of *σ*_{Y}/*E*, however, they found that the true contact area was related to the load by , in excellent agreement with our approximate estimate.

The convergence of contact size and contact area to their limiting values are illustrated in figure 4. Since the normalized contact size as *n*→∞, while , we have plotted and as functions of *n*. Evidently, both contact size and total contact area converge to their asymptotic values exponentially: indeed, it is straightforward to show that for sufficiently large *n*, the contact size and total contact area may be approximated by(4.35)

These results can also be used to deduce the asymptotic behaviour of fine-scale roughness on an elastic–plastic solid. We saw in §4*c* that even if the bulk contact and the first few roughness scales remain below yield, the contact pressure on an elastic surface increases for progressively higher scales of roughness. Consequently, there will always be a critical scale of roughness that reaches yield. By the same argument, there will always be a critical roughness scale which reaches the fully plastic regime. For higher (shorter-wavelength) scales of roughness the graphical method illustrated in figure 3 can again be used to deduce the contact pressure and contact fraction. For sufficiently large *n*, the true contact area and contact fraction will again approach their asymptotes as shown in figure 4.

It is interesting to note the important role played by interactions between neighbouring asperities in determining the asymptotic behaviour of fine scale roughness. Without this interaction, the contact pressure on each scale of roughness is bounded by the material hardness (i.e. the limiting pressure under an isolated indenter). For a rigid plastic contact, therefore, all roughness scales (including *n*=0) must reach full contact to support the load. For a more general elastic–plastic contact, the total true contact area and the true contact spot size would be determined by the first roughness scale to reach the fully plastic limit, since all shorter-wavelength scales of roughness must be at full contact to support the contact pressure. The true contact spot size could then be computed unambiguously. This, in fact, is the simple picture of plastic contact between two rough surfaces proposed by Tabor (1951), although his focus was on computing the total area of contact, and did not address the question of the contact spot size.

## 5. Results

We now discuss in detail the predictions of the calculations outlined in the preceding sections. We consider a surface with Weierstrass profile defined in equation (2.1), which is indented by a cylinder of radius *R* subjected to a load *P* per unit length out of plane, as illustrated in figure 1. It is convenient to characterize the applied load, surface geometry and material properties using three dimensionless groups:(5.1)

Here, *Σ* is a dimensionless load parameter, and is proportional to the mean nominal contact pressure. The parameters *Ψ* and *K* together characterize the surface roughness and its resistance to plastic flow.

As a preliminary step, we express several parameters of interest in terms of *Σ*, *Ψ* and *K*. If the nominal contact is elastic () the contact size is , while for a fully plastic nominal contact () . The wavelength of the zeroth roughness scale (see figure 1) follows as , so we find that the plasticity index (defined in equation (3.3)) for the zeroth roughness scale is given by(5.2)

### (a) Conditions for inducing yield

Controlling plastic flow under contact loading is a central issue in designing failure resistant surfaces for tribological applications. With this motivation, we begin by estimating the conditions necessary to initiate yield in our idealized model of a surface with fractal roughness. With reference to figure 1, we focus attention on each successive scale of roughness under the contact in turn. Strictly, our calculations show that short-wavelength roughness scales will *always* deform plastically, so the surface does not have an elastic limit. Instead, our goal is to calculate the number of roughness scales that remain elastic, as a function of the load (parameterized by *Σ*) and surface roughness and material properties, parameterized by *Ψ* and *K*.

We assume for simplicity that the nominal contact remains elastic, which requires *Σ*≤3.2 (Johnson 1985). The most heavily loaded region of the zeroth scale of roughness occurs at the centre of the nominal contact, and is subjected to a normalized contact pressure . The elastic solution for indentation of a sinusoidal surface (§3) shows that the zeroth roughness scale will behave in one of two ways, depending on the value of *ψ*_{0}. For *ψ*_{0}<0.2, the zeroth roughness scale will remain elastic up to full contact, which occurs at a critical load . For *ψ*_{0}>0.2, the zeroth scale of roughness will reach yield before full contact. The critical nominal contact pressure at yield is given by .

Similar arguments can be applied to each successive scale of roughness . The maximum nominal contact pressure acting on each roughness scale can be computed from equation (4.31); the criteria and give the conditions for full contact and yielding on the *n*th scale of roughness, respectively.

The results of these calculations are summarized in figure 5, which shows the critical normalized load *Σ* required to: (i) just cause yield in each roughness scale (shown as solid lines) and (ii) just bring each roughness scale to full contact (shown as dashed lines). The results depend (weakly) on *D* and *γ*, and are shown for the particular case *D*=1.5, *γ*=5. To understand the significance of figure 5, it is helpful to visualize the behaviour of a contact under progressively increasing load. At extremely low values of applied load, the nominal contact width is small, and so all roughness scales under the contact have short wavelength. Since, short wavelengths in the Weierstrass series have a high amplitude/wavelength ratio, these roughness scales easily deform plastically, and so only the zeroth scale of roughness (with longest wavelength) under the contact remains elastic. As the load is increased, the nominal contact width increases, and so additional roughness scales are brought within the contact. These roughness scales have a longer wavelength, and hence a lower amplitude/wavelength ratio, making them more resistant to plastic flow. Consequently, for a sufficiently large load, and sufficiently smooth surface, the new roughness scales may remain elastic. For example, when a surface with is subjected to a load *Σ*=0.1, scales *n*=0, 1, 2 are all elastic, while all shorter-wavelength roughness scales exceed yield. If the load is increased still further, the roughness scales with longest wavelength under the contact may reach full contact. For the surface with , scale *n*=0 reaches full contact at , scales *n*=0,1 reach full contact shortly thereafter. If a roughness scale reaches full contact without yielding, it remains elastic until the plastic zone associated with the nominal contact reaches the surface.

The parameter evidently plays the role of a plasticity index for the surface: a contact with a low value of has a large number of roughness scales that deform elastically up to the load where the nominal contact yields, while a contact with high will experience extensive plastic flow even at modest loads.

For practical purposes, it is of particular interest to estimate the depth of the plastically deformed region under the contact. As a rough guide, the depth of the plastically deformed region can be taken to be one half the wavelength of the first roughness scale to reach yield, since longer wavelengths remain below yield, while all shorter wavelengths will exceed yield. This approximation gives , where *k* is the number of roughness scales that remain in the elastic regime, which may be estimated from figure 5. The normalized depth of the plastically deformed region as a function of the load parameter *Σ* is shown in figure 6, for various values of the plasticity index , for *D*=1.5, *γ*=5. We observe three general regimes of behaviour: for large values of , only the zeroth roughness scale remains elastic at light loads, so the depth of the deformed layer is proportional to *Σ*. For , the zeroth roughness scale yields at large loads, which results in a sudden increase in the depth of the deformed layer. For , the zeroth scale remains elastic up to full contact, while all remaining scales yield—in this case there depth of the deformed layer is proportional to *Σ* over the full range. Finally, for , the one or more roughness scales transition to elastic behaviour at sufficiently large loads. Under these conditions the depth of the deformed layer is proportional to the load parameter *Σ* at light loads, but is independent of *Σ* above a critical load.

### (b) Contact pressure distribution functions

We proceed next to examine the contact pressure distributions in more detail. The pressure acting on each scale of contact is completely characterized by the distribution functions *q*_{n}(ξ) specifying the fraction of the contact. For an elastic surface, *q*_{n}(ξ) are smooth functions (see e.g. figure 5 of Ciavarella *et al*. 2000). For plastic contacts, *q*_{n}(ξ) contains singularities (this is a consequence of the fact that contact pressure under fully plastic contacts is uniform—see electronic supplementary material for a more detailed discussion). Consequently, it is more convenient to plot the cumulative contact pressure distribution functions , which specify the proportion of each scale of contact over which the pressure exceeds *p*/σ_{Y}. The cumulative contact pressure distributions can be computed from the definition . Results are shown in figure 7 for a dimensionless nominal pressure *Σ*=0.001 and two values of surface roughness, with *γ*=5 and *D*=1.5. figure 7*a* shows results for *ψ*_{0}=0.0002, in which case a large number of scales under the contact are elastic. Under these conditions, the pressure distributions for roughness scales *n*=0, …, 9, which deform elastically, tend towards a self-similar distribution. (Ciavarella *et al*. 2000). The most heavily loaded region of the 10th roughness scale begins to deform plastically, while the 11th reaches the fully plastic limit (which occurs when the contact pressure reaches 2.9*σ*_{Y}), and begins to transition towards a uniform distribution. The pressure acting on subsequent scales becomes increasingly uniform, and gradually increases towards the limiting value 5.8*σ*_{Y}. A contrasting result for a rougher surface with *ψ*_{0}=0.1 is shown in figure 7*b*. In this case, only scales *n*=0, …, 1 are elastic; scale *n*=2 exceeds yield, and a small fraction of scale *n*=3 reaches the fully plastic limit. Scale *n*=4 is dominated by contacts in the fully plastic regime. The pressure acting on subsequent scales of roughness again becomes progressively more uniform, and increases towards the limit 5.8*σ*_{Y}. Pei *et al*. (in press) have noted the tendency of plasticity to drive contact size and pressure distributions to uniform distributions in numerical simulations (they refer to plasticity as a ‘grand equalizer’). Our model suggests an explanation for this phenomenon. The behaviour of a sinusoidal surface is shown in figure 2, and is approximated by equations (A2)–(A18) in the electronic supplementary material. Evidently, nominal pressures in the range acting on the sinusoid all induce the same true contact pressure. The first roughness scale to reach fully plastic behaviour, therefore, exerts a nearly uniform pressure on the next scale of roughness.

### (c) Average contact size, contact pressure, total contact area, and number of contacts for each scale

The overall behaviour of the contacts at each scale may be computed using the distribution functions given in §5*b*, and may also be estimated using the approximate procedure outlined in §4*b*. Figures 8 and 9 show the results of these calculations. Figures 8*a* and 9*a* show the contact fraction and the ratio of true to nominal contact area for successive roughness scales *n* under a contact; while figures 8*b* and 9*b* show the total number of contacts on each successive scale, normalized by the number of contacts in the zeroth roughness scale, *N*_{n}/*N*_{0}, and the dimensionless nominal contact pressure *P*/σ_{Y}*A*_{n} on each scale. All variables have been plotted as a function of to illustrate their asymptotic behaviour clearly. Results are shown in figure 8 for a nominal contact pressure *Σ*=0.001 and for several values of surface roughness (parameterized by ); while figure 9 shows results for a fixed dimensionless surface roughness and progressively increasing load *Σ*.

To interpret the figures, note that each curve labelled (a)–(d) corresponds to a different surface. The circles on each curve show the values of variables on successive scales of roughness *n*=0, 1, 2, …, obtained using the procedure outlined in §4*b*, while the solid lines show the exact procedure described in §4*a* and electronic supplementary material. As a representative example, consider case (*a*) in figure 8. This is a very smooth surface, so the zeroth roughness scale (the left-most point on the curve) has a relatively high contact fraction . Note also that the contact area and number of contact spots for this scale must always satisfy and *N*/*N*_{0}=1. Finally, since the surface is relatively smooth, and the value *Σ*=0.001 corresponds to a light load, the nominal contact pressure on the zeroth scale for case (*a*) has a low value . Successive circles on the curve marked (a) show the behaviour of subsequent roughness scales in the same way. For this surface, a large number of roughness scales (10) remain elastic, so the contact fraction and area for these scales approach the elastic asymptotic limits and , respectively, while the number of contact spots and nominal pressure approach and . The transition from elastic to elastic–plastic behaviour occurs at , marked as a chain line in figure 8 (this approximation is strictly valid only if the contact size reaches its elastic asymptote). For case (*a*) the 11th and 12th roughness scales exceed yield, and are loaded in the elastic–plastic regime. Finally, the transition to fully plastic behaviour occurs near *ψ*_{n}≈2, and for case (*a*) scales 13 and higher are in the fully plastic regime. For these scales, the contact size and contact area approach the plastic asymptotes: ; , while the number of contact spots and the normalized nominal pressure approach and . Contrasting behaviour for a rough surface under the same load is shown in curve (d) in figure 8. In this case, only a scale of roughness (*n*=0) is elastic; scale *n*=1 reaches yield, so the elastic asymptote is not reached. Scales *n*=4 and higher exhibit fully plastic behaviour.

The behaviour of roughness scales under a surface with fixed roughness, which is indented by a cylinder under progressively increasing load is illustrated in figure 9. Results are shown for roughness and five different loads. At the lightest load (*Σ*=10^{−5}) only the zeroth scale of roughness is in the elastic regime; there are 2–3 scales in the elastic–plastic regime, but the majority of roughness scales lie in the fully plastic regime and the contact size, area, pressure and the number of contact spots quickly approach their fully plastic asymptotes. As load is increased, additional roughness scales with a lower ratio of amplitude to wavelength are brought within the contact, so the number of scales in the elastic regime increases with load. At *Σ*=0.03, there are six roughness scales in the elastic regime: a sufficient number to allow the contact size, area, pressure and the number of contact spots approach their elastic asymptotes. Shorter wavelengths again transition to the fully plastic regime.

## 6. Discussion

Ciavarella *et al*. (2000) showed that the contact between two elastic, perfectly fractal surfaces consists of an infinite number of contact spots, with zero size, subjected to infinite contact pressure. This paper was motivated by the hope that accounting for plastic flow would resolve these difficulties. We have found, however, that the total contact area and contact pressure between elastic–plastic solids are well defined, but it remains impossible to predict the actual number of contacts or their size. The ambiguous nature of the true contact could be resolved in either of two ways: (i) it is possible that some physical process, which was not included in our contact model, can regularize the ill posed nature of contact even for perfectly fractal surfaces. For example, Persson *et al*. (2005) shows that accounting for adhesion has this effect, for surfaces with *D*<1.5. (ii) No real surface can be perfectly fractal, and provided that the surface geometry can be measured, and characterized, at extremely short wavelengths where the fractal description breaks down, standard techniques—including perhaps the Greenwood–Williamson method—can be used to model the contact. Available experimental studies (Majumdar & Bhushan 1990; Bora *et al*. 2004) provide strong evidence that surfaces cease to be fractal at nanometre length scales.

Although, it may ultimately be necessary to abandon perfectly fractal descriptions of surface roughness, it is clear that a fractal description of a surface has several advantages over the conventional approach; and that some features of fractal behaviour must be retained in a complete analysis. In particular, the fractal description captures the important fact that asperity contacts are spatially correlated. In the Greenwood–Williamson description, the focus is entirely on the smallest scale of contacts, which are assumed to be uniformly distributed spatially, and to exert a nominally uniform pressure. In contrast, the fractal description, together with the Ciavarella model, allows one to calculate the collective effect of the contact clusters, without needing to describe the contacts themselves in detail. There are also some applications where the behaviour of the surfaces is determined by long-range surface roughness, and is insensitive to the detailed nature of the contact spots themselves. For example, contact compliance, and contact resistance, are determined by long-range roughness, and may converge even for elastic surfaces with fractal roughness (Ciavarella *et al*. 2004; Manners & Gholami 2005)

Our analytical calculations also help to elucidate the behaviour that one would expect to observe in a numerical simulation of contact between rough, ideally fractal, surfaces. The calculations of Ciavarella *et al*. (2000), as well as those reported here, raise some questions regarding the interpretation of such simulations. It appears to be an inescapable conclusion that the true contact area for elastic contact between two perfectly fractal surfaces consists of an infinite number of contact spots of zero size, which are subjected to infinite contact pressure, and have a vanishing total area of contact. Numerical simulations that predict contact sizes, contact pressures, or the true area of contact between two elastic solids are, therefore, valid only if: (i) the actual surface geometry ceases to be fractal at some short wavelength (e.g. the power-spectral density could be truncated at some short wavelength; or else decay sufficiently rapidly to give a differentiable autocorrelation function at the origin); and (ii) the resolutions of both the surface profilometre used to measure the three-dimensional surface, as well as the numerical simulations, are sufficient to detect and to characterize surface geometry at these short wavelengths. Numerical simulations of elastic–plastic contact are somewhat less sensitive to artificial truncation of the roughness spectrum, but the predicted total number of contacts, and their size, would not converge for a perfectly fractal surface.

The qualitative effects of introducing an artificial short-wavelength cut-off into numerical simulations of contact between fractal surfaces can be estimated using our simple model. As a representative example, we consider a Weierstrass surface (equation (2.6)) containing only terms *n*=0,1,2,…, 5, which is indented by a perfectly flat, rigid indenter under a uniform nominal pressure *p*_{−1}. We restrict attention to modest nominal pressures, so that the true area of contact is much less than the nominal contact area. Under these conditions, the true contact consists of a well-defined set of contact spots. The average size of these spots, and the true contact pressure, are essentially independent of nominal pressure, while the number of contact spots and the total true contact area increase in proportion to the nominal pressure. The contact size and mean contact pressure, as well as the relationship between true contact area and the number of contact spots, is determined by the elastic and plastic properties of the surface and the surface roughness, which may be characterized together by *ψ*_{0}. We find that for *ψ*_{0}<0.006, the surface deforms elastically: under these conditions the behaviour of the contacts can be estimated using the expressions given in §4*c*, with *n*=5. For *ψ*_{0}>0.006, the contact spots exceed yield. The effect of plastic flow is to increase the size and number of contacts, to increase the true area of contact, and reduce the mean pressure. As *ψ*_{0} increases, the true contact progressively transitions to the limiting conditions for a rigid plastic surface: the contact fraction for the fifth scale approaches *a*_{5}/*λ*→1, the true contact pressure approaches 〈*p*_{5}〉=5.8*σ*_{Y} and the ratio of true to nominal contact area approaches . The transition with *ψ*_{0} is gradual, however. To illustrate this transition from elastic to plastic behaviour, we have plotted the variation of the true contact pressure *P*/*A*_{n}*σ*_{Y}*ψ*_{0}=*P**λ*_{0}/*A*_{n}*E*^{*}*y*_{0} as a function of 1/*ψ*_{0}=*σ*_{Y}*λ*_{0}/*E*^{*}*y*_{0} in figure 10 for a surface with *D*=1.5 and *γ*=5. Results are shown for surfaces truncated at scales *n*=5, 8 and 11. We have chosen to plot the results in this form to compare our predictions with the detailed finite element simulations of indentation of a nominally fractal surface reported by Pei *et al*. (in press). In particular, their fig. 4 shows the predicted true contact area as a function of the ratio of yield stress to elastic modulus *σ*_{Y}/*E*. The numerical results are in excellent qualitative agreement with our estimate. They observe a clear transition from elastic behaviour at high *σ*_{Y}/*E* to plastic behaviour at low *σ*_{Y}/*E*; moreover, they predict in the limit , in excellent agreement with our estimate of . In our calculations, the elastic–plastic transition is a direct consequence of truncating the Weierstrass series, however. As the series is truncated at progressively shorter wavelengths, the elastic regime disappears. These results illustrate the challenges associated with numerical simulations of contact between fractal surfaces.

The simple model we have developed here could be extended in several ways. For example, it would be straightforward to extend the calculations described here to three dimensions. In addition, it would not be difficult to extend our calculations to account for the scale dependent nature of plastic flow: due to indentation size effects (Nix & Gao 1998), we might expect short-wavelength roughness to return to elastic behaviour. It would be of great interest to determine where this transition occurs, because then it would be possible to design surfaces that remain elastic. Finally, the general procedure employed here is clearly not restricted to fractal surfaces, and could be applied to any discrete approximation to a power spectral density, provided that the discrete scales remain widely separated.

A more important issue, however, is the need to adopt a more realistic model of the surface roughness. The Weierstrass profile gives a convenient way to characterize the multi-scale nature of surface roughness, but is clearly not a realistic model of surface roughness. Furthermore, the scale decoupling introduced by Ciavarella *et al*. (2000), which forms the basis of our analysis, strictly holds only in the limit *γ*→∞. A number of questions, therefore, need to be addressed. Firstly, it is of interest to determine the range of *γ* for which the scale decoupling can be applied. For elastic surfaces, this is straightforward. The scales will decouple as long as there are sufficient contacts of scale *n*+1 within each contact of scale *n*. Recall that (for sufficiently large *n*), the contact spot size for the *n*th scale can be estimated as . The number of contact spots of scale *n*+1 within each contact of scale *n* follows as . Consequently, the critical value of *γ* to ensure *m* contact spots is . Taking *m*=2 as an optimistic estimate requires *γ*≥4; in practice, an estimate of *γ*≥10 is more reasonable. It is difficult to obtain similar analytical estimates for plastic contacts, so instead, we have investigated this issue using finite element computations with a finite number of roughness scales, some of which are described in electronic supplementary material. Simulations with two and three scales of roughness suggest that scales of roughness completely decouple for *γ*≥10. For smaller *γ*, we find that a few asperities in the *n*th scale that are located near the peaks of scale *n*−1 may coalesce, in contrast to the large *γ* approximation. This is because at light loads, the *n*th scale contains only a few contacts, and insufficient lateral constraint develops to increase the contact pressure for the finest scale. However, asperities on scale *n* that are located in the valleys of scale *n*−1 are strongly constrained, and so exhibit behaviour that is in good agreement with the behaviour of a surface with purely sinusoidal roughness. Consequently, the predictions of the scale decoupling analysis continue to predict qualitatively the *average* behaviour of each roughness scale even for low values of *γ*. Nevertheless, further work is required to properly characterize the behaviour of random surfaces, for both elastic and plastic solids.

## 7. Conclusions

We have outlined an approximate model that is intended to estimate the behaviour of a rough elastic–perfectly plastic solid with (two-dimensional) self-affine fractal surface roughness subjected to contact loading. The surface was idealized by a Weierstrass profile (equation (2.1)) with widely separated scales of roughness (*γ*≫1). Under these conditions, each scale of roughness can be regarded as a sinusoidal surface, indented by a flat surface. Successive scales are coupled through the condition that the nominal pressure acting on scale *n* is equal to the true pressure acting on scale *n*−1. On this basis, a recursive expression can be found for the contact pressure distribution, the contact size, and the total area of contact for each scale of roughness, in terms of the response of a sinusoidal surface (Ciavarella *et al*. 2000). In our calculations, the response of an elastic–plastic solid with sinusoidal roughness was approximated by fitting curves to finite element computations, as outlined in §3. We found that simple but highly accurate estimates for the average contact fraction, pressure and area of contact on successive roughness scales can be obtained by replacing the pressure distribution acting on each scale by a statically equivalent uniform distribution.

We have used this approach to analyse the behaviour of an elastic–perfectly plastic Weierstrass surface when indented by a rigid cylinder. We found that the behaviour of the solid is governed by three dimensionless functions, defined in equation (5.1). The parameter *Σ* quantifies the magnitude of the applied load; while *Ψ* and *K* together characterizes the surface roughness and its resistance to plastic flow. The conditions necessary to induce plastic flow in successive roughness scales under the cylinder are illustrated in figure 5. For light loads, the nominal contact area is small, so roughness scales under the indenter have a large amplitude compared with their wavelength, and consequently a low resistance to plastic flow. As a result, at light loads only the zeroth scale of roughness deforms elastically; all other roughness scales are fully plastic. As the applied load is increased, the nominal contact area spreads, and new roughness scales that come into the contact area may remain elastic, and for sufficiently high loads will reach full contact. The parameter was identified as a plasticity index for a surface with fractal roughness. Surfaces with low values of have a large number of scales of roughness in the elastic regime, and consequently plastic flow is restricted to a thin layer at the surface, figure 6. As increases, the number of plastic roughness scales, and the corresponding depth of the plastically deformed layer, increase.

We also calculated the average contact fraction, total contact area, number of contacts, and average contact pressure for successive scales of roughness. The results are illustrated in figures 8 and 9. The behaviour of the surface depends primarily on the dimensionless parameter *ψ*_{0} defined in equation (3.3), which characterizes the resistance to plastic flow of the longest-wavelength scale of roughness under the nominal contact. This parameter can be computed in terms of *Ψ* and *K* using equations (5.2). For low values of *ψ*_{0}, a large number of long-wavelength roughness scales under the contact remain elastic: their behaviour can be calculated using the closed-form expressions given in §4*c*. Short-wavelength roughness scales always deform plastically, and for large *n* approach the asymptotes given in §4*e*. For an elastic–plastic solid with a perfectly fractal surface, the true area of contact consists of an infinite number of contact spots of zero size, with a true contact pressure . The true contact area is bounded, and proportional to the load.

The implications of these findings were discussed in detail in §6. Our calculations suggest that idealizing surface roughness as a perfect fractal leads to unphysical predictions for the contact size and number of contact spots. It is likely that these difficulties can only be resolved by taking into account the departure of actual surfaces from perfectly fractal behaviour at very short wavelengths.

## Acknowledgements

This work was performed as part of the General Motors Collaborative Research Laboratory on Computational Materials Research at Brown University. The authors are grateful to M. Ciavarella, R. W. Carpick and J. F. Molinari for helpful discussions and for providing pre-prints of their recent publications.

## Footnotes

The electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2005.1563 or via http://www.journals.royalsoc.ac.uk.

- Received June 1, 2005.
- Accepted August 3, 2005.

- © 2005 The Royal Society