Contact area of rough spheres: Large scale simulations and simple scaling laws

Contact, friction, and adhesion control the behavior of mechanical devices from cars¹ to artificial joints² and determine the properties of granular matter such as colloids or sand piles.³ The canonical geometry used to measure these quantities is a sphere of radius R pushed into a flat substrate with a normal force N. Recent metrology devices with this geometry have the ability to measure relative motion of surfaces with nanometer resolution, including the surface force apparatus (⁠$R∼10 mm$⁠),⁴ colloidal probes (⁠$R∼10 μm$⁠),⁵ and atomic force microscopes (⁠$R∼10 nm$⁠).⁶ The surfaces in these devices generally have small scale roughness superimposed on the macroscopic curvature that can dramatically alter the contact area, friction, tangential stiffness, adhesion, and transport properties. There has been significant recent progress in treating roughness on nominally flat surfaces but not on the interaction between roughness and macroscopic curvature.

In this paper, we use molecular simulations to study contact of rough spheres with R from 30 nm to $30 μm$⁠. By combining recent results for rough flat surfaces^7–14 with classical results for smooth spheres,^15–20 we develop parameter free equations that describe the full variation of the area of molecular contact A with N. A is difficult to obtain from experiments, but it plays a central role in determining contact properties. The variation of A with N determines whether there is adhesion, and the magnitude and geometry of contact area are widely used to calculate stiffness, transport, and the friction force F. A simple approximation that appears to be valid in many cases, even down to sliding AFM tips,^21–23 is that of a constant interfacial shear stress τ,²⁴ giving a friction force $F=τA$⁠.

More than a century ago Hertz showed¹⁵ that contact of a rigid smooth sphere with radius R on an elastic flat of contact modulus $E*$ occurs in a circle with radius

The work of Hertz has been subsequently extended to many other situations, such as shear loading, interfacial friction, and adhesion.^16–20 

Roughness reduces the area in intimate repulsive contact since only the highest asperities are able to touch. We therefore expect that the true contact area A is less than $AHertz≡πaHertz2$ in most situations. The prevalent theory of contact of a rough sphere on a flat was presented by Greenwood and Tripp.²⁵ Their analysis is based on the traditional Greenwood-Williamson (GW) model²⁶ for nominally flat surfaces that approximates the rough surface by a collection of noninteracting Hertzian contacts of identical radius. In the last decade, analytical⁷ and numerical^8–11 works have revealed important limitations of GW results for total contact area¹³ and contact geometry^8,27,28 motivating us to revisit the contact of rough spheres.

Recent studies of nonadhesive contact between nominally flat surfaces of area A₀ all find a linear relation between A and N up to fractional contact areas $A/A0∼10%$⁠.^7–14 This implies a constant normal pressure $prough=N/A$ within the contacting regions. Roughness enters only through the dimensionless root mean square (rms) slope of the surface $hrms′≡〈|∇h|2〉1/2$⁠, where $h(r→)$ is the height as a function of position $r→$ in the plane. Numerical studies find $prough≡hrms′E*/κ$⁠, where the dimensionless constant $1/κ≈1/2$ in the continuum, hard-wall limit.^8–13 The constant $1/κ$ becomes smaller for finite range soft repulsion and adhesive interactions. Negative values are found for strongly adhesive interactions and correspond to macroscopic adhesion or stickiness.¹⁴ At high loads, the surfaces are pushed into compliance and A = A₀. Over the whole range of loads, the fractional area of contact $f=A/A0$ approximately follows an error function:^7,12,29,30

where $p=N/A0$ is the nominal pressure pushing the two surfaces together.

Here, we extend these calculations to contact of a rough rigid sphere of radius R and a semi-infinite elastic substrate of contact modulus $E*$⁠. The results can be mapped to contact of two elastic rough spheres using a standard mapping.¹⁸ The rigid surface has height $hR+h$⁠, where the first term is the average curvature and the second is the random roughness. As in flat surface studies,^8–13 h has a self-affine fractal form from lower wavelength λ_s to upper wavelength λ_L with Hurst exponent H.³¹ 

The calculations use an efficient Green's function method^32–35 supplemented by a padding region that prevents any effect from repeating images.³⁶ We use the Greens function for an isotropic solid with Poisson ratio 1/2.¹⁸ Atoms on the flat elastic substrate are arranged on a square grid with spacing a₀. This can be regarded as the (100) surface of an fcc lattice or a convenient discretization of the continuum equations of elasticity. We show results for different interfacial interactions with $λs=4a0, λL=512a0$ and H = 0.8, but simulations with other parameters were carried out with no fundamental differences in the results. Contact area is obtained from atoms that feel a repulsive load.¹⁴ Contact radii up to 1024 atoms were studied, corresponding to $∼0.3 μm$⁠.

Fig. 1 shows the result of a nonadhesive calculation with hard-wall interactions between the rough sphere and elastic flat.³³ Initial contact occurs at the highest protrusion which is not necessarily located at the bottom of the equivalent smooth sphere. In the particular case shown in Fig. 1(a), contact occurs outside of the Hertz contact radius $aHertz$ corresponding to the given load (solid black line).

One can estimate the radius $a′$ of the first asperity to contact by setting the nominal separation between a smooth sphere and substrate equal to the rms height $hrms$ of the roughness: $a′2=2Rhrms$⁠. For $a′<λL, hrms$ depends on the scale $a′$ over which it is measured.³⁷ Solving for $a′$ gives the analytical expression

This radius is shown by the dashed line in Fig. 1(a).

The bottom row of panels in Fig. 1 shows two measures for the radially averaged pressure in the contact. The black line shows the full average pressure, $p¯(r)$⁠, that includes noncontacting regions with zero pressure in the average. As the normal force increases [going from panel (a) to panel (d)], $p¯(r)$ increases but is close to Hertz's expression (thin solid line) in all cases. The largest deviation from Hertz's solution is near the contact edge because part of the contacting area sits outside $aHertz$⁠.

Also shown in Fig. 1 is the average contact pressure, $p¯c(r)$⁠, where in contrast to $p¯(r)$ non-contacting regions are excluded from the average. For low loads, $p¯c$ is always close to $prough$ and is independent of radius and load (see Figs. 1(a) and 1(b)). At larger loads, $p¯c$ rises above $prough$⁠. The upper panels of Fig. 1 show that this increase in $p¯c$ occurs when the surfaces are being pushed into nearly complete contact. In this limit, $p¯c$ becomes identical to $p¯$ and both follow the Hertzian form.

The numerical results for the variation of total contact area with load are accurately described by simply assuming that the fraction of the Hertzian area that contacts is given by the result for a nominally flat rough surface at the same mean pressure $N/AHertz$

While the pressure profile predicted by Hertz varies with r (see Fig. 1), the associated distribution of local pressures is peaked around the mean pressure $N/AHertz$⁠. Given this assumption, there should be a crossover from randomly rough to Hertz scaling when $N/AHertz=prough$ or at a load

Figure 2(a) shows that this simple scaling assumption (thin solid line) captures numerical results for hard wall interactions over a wide range of sphere radii. For each load, the area $〈A〉$ is averaged over five realizations of a rough sphere and normalized by $AHertz(N)$⁠. When the load is normalized by N_c, data for different R collapse onto Eq. (4) at intermediate to large loads. Note that the collapse covers up to 10 orders of magnitude in both $〈A〉$ and N, since $AHertz∝N2/3$⁠.

Fig. 2(a) reveals yet another scaling regime at the lowest loads where $A/AHertz=const.≪1$⁠. Here, contact is dominated by the first rough peak to make contact. The peak radius ρ can be estimated from the rms curvature $〈|∇2h|2〉$ of the rough surface as³¹ 

The gray horizontal lines in Fig. 2(a) show that the result of this calculation gives an excellent agreement with the numerical data at low loads. Note that Eq. (5) with R replaced by ρ also describes the crossover from this first contacting asperity to multiasperity contact. The multiasperity regime where $A∝N$ is suppressed for small spheres (⁠$ρ∼R$⁠) because there is never a statistical number of asperities.

Real interatomic interactions have a finite stiffness and an adhesive tail.³⁸ A common model for the force between surfaces is a 9–3 Lennard-Jones interaction that is truncated at its minimum, so it is purely repulsive. Fig. 2(b) shows results for this model with the strength of the potential chosen so that the interaction has the same stiffness as the bulk (see supplementary section S-II). We find that the area is generally larger than for the hard-wall case. At intermediate loads, the increase is captured by Eq. (3) with $1/κ$ decreased to $1/κrep$⁠, consistent with past studies of nominally flat surfaces with soft repulsion.^10,12,39 The finite compliance of the interfacial repulsion leads to variations in surface separation, and supplementary section S-2³¹ provides an analytic expression for the resulting changes in $1/κrep$⁠. Fig. 2(b) shows that the renormalized value of $κrep=2.9$ (thin black line labeled $A∝N$⁠) obtained from this expression agrees well with the numerical data in the linear scaling regime. The contact of the first asperity is approximately Hertzian. Compliance increases the effective radius to $ρeff≈3ρ$ (gray horizontal lines in Fig. 2(b)).

Finally, we consider longer-range attractive interactions that favor adhesion between sphere and substrate.³¹ This introduces a new quantity, the work of adhesion w which is the energy gained per unit contact area.^19,38 The competition between adhesion and elastic deformation is quantified by the elastocapillary length $w/E*$⁠, which is normally much less than an atomic spacing.^14,40 The main effect of adhesion in the intermediate $A∝N$ scaling regime¹⁴ is to approximately renormalize the value of $1/κ$ to $1/κadh=1/κrep−1/κatt$⁠. The value of $1/κatt$ increases with $w/E*$⁠, and there is macroscopic adhesion when $κadh$ becomes negative. Using the analytical theory of Ref. 14, we find $1/κatt≈0.04$ and $1/κatt≈0.25$ for $w/E*=0.0005a0$ and $w/E*=0.0050a0$⁠, respectively. These values correspond to van der Waals adhesion (⁠$w∼50 mJ m−2$⁠) on stiff ceramic substrates (⁠$E*∼300 GPa$⁠) and metals (⁠$E*∼30 GPa$⁠), respectively.

Renormalizing κ works well for the weakly adhesive case where $1/κatt≲1/κrep$⁠. For the strongly adhesive case, we replace the expression for the fractional area, Eq. (2), with tabulated values obtained from a reference calculation on nominally flat surfaces and an identical interaction potential. The thin black lines in Fig. 2(c) show that the result of this procedure is in excellent agreement with our calculations of adhesive spheres, even in strongly adhesive cases. The influence of the curvature of the macroscopic sphere, and hence, the change in nominal contact area, is still captured well by Hertz's expression as long as $κadh$ remains positive.

The numerical data for adhesive interactions show that the contact of the first asperity leads to different asymptotic behavior than in the nonadhesive case and much larger variations between different realizations of the random surfaces. For an analytical description of the asymptotic behavior at low loads, Hertz's expression needs to be replaced by its adhesive counterpart and the Johnson-Kendall-Roberts (JKR)¹⁶ or Derjaguin-Muller-Toporov (DMT)¹⁷ models are commonly used. These models present bounding limits of continuum contact for soft substrates with large w and stiff substrates with small w, respectively. The general intermediate case can be described by the Maugis-Dugdale solution.²⁰ Our numerical data are well-described by simply assuming a constant contact area $Aadh$ for the first asperity, which all of these models predict at low positive load. The results of fitting $Aadh$ individually to our ensemble of five rough surfaces are shown by the gray area in Fig. 2(c).

For all of the surface interactions in Fig. 2, the range of linear scaling (⁠$A∝N$⁠) decreases as $R/λs$ and $h′rms$ decrease. Linear scaling requires a statistical number of asperities,⁸ and this requires a large separation between R and the size of the first asperity to contact. The smallest $R=100a0∼30nm$ corresponds to a typical atomic force microscope tip. In contrast to some earlier results,⁴¹ we find no significant linear region for tips that are this small, even when the smallest scale of roughness is an atomic diameter and an artificial hard sphere interaction is used on this scale. In contrast, there is a substantial linear region for $R=105a0∼30 μm$⁠, which is typical of colloidal probe experiments or characteristic asperity radii inferred from macroscopic friction experiments.⁴² 

Our results have important consequences for fits to experimental data on a wide range of scales. At low loads and contact areas, recent results for rough flat surfaces can be applied to rough spheres.^7–14 Smooth sphere theories^15–20 can be used to fit experimental data in the limit of full contact, i.e., when the contact radius is larger than the value corresponding to N_c, $a/R=(3π/4κ)hrms′$⁠. Fit values of w will, in general, be reduced from the intrinsic interaction energy due to the energy cost of flattening the roughness.^43,44 Excellent fits to JKR theory have been obtained for atomically flat mica and very smooth elastomer surfaces with R from 1 to $10 mm$ and a/R from 0.001⁴⁵ to 0.1.^46,47 Many AFM experiments with $R∼30$ nm are also consistent with DMT or JKR theory.^21–23 This is consistent with our finding that there is no intermediate $A∝N$ region for these small tips (Fig. 2). Of course, the roughness may affect the effective radius that should be used in extracting w.⁴⁸ 

In summary, we have developed simple analytic expressions for the contact of rough spheres and tested them with large scale simulations with tip radii from $30 nm$ to $30 μm$ over 10 orders of magnitude in load. Our equations combine established continuum theories for smooth spheres^15–20 and recent results for nominally flat but rough surfaces.^7–14 The latter are extended from the ideal limit of hard wall interactions to include finite compliance and adhesion which become important at nanometer scales. The results provide detailed predictions for the functional form of the contact area as well as simple equations for estimating which regime describes an experiment or simulation. Above a crossover load $Nc$ (Eq. (5)), the contact area follows continuum theories for smooth spheres using the macroscopic radius of curvature R. The same theories can be used at very small loads, but with R replaced by the radius ρ of the first asperity to make contact (Eq. (6)). When $R≫ρ$⁠, there is an intermediate range of pressures where the area is proportional to load and is consistent with results for rough flat surfaces. A simple criterion describes whether surfaces adhere in this limit.¹⁴ 

This material is based upon work supported by the Deutsche Forschungsgemeinschaft under Grant No. PA 2023/2 and U.S. National Science Foundation Grants Nos. DMR-1411144, OCI-0963185, and CMMI-0923018. We additionally thank the Jülich Supercomputing Center for providing computer time (Project No. hfr13) and the Large Scale Data Facility at KIT for providing storage space.