Boundary lubrication of heterogeneous surfaces and the onset of cavitation in frictional contacts

See allHide authors and affiliations

Science Advances  25 Mar 2016:
Vol. 2, no. 3, e1501585
DOI: 10.1126/sciadv.1501585


Surfaces can be slippery or sticky depending on surface chemistry and roughness. We demonstrate in atomistic simulations that regular and random slip patterns on a surface lead to pressure excursions within a lubricated contact that increase quadratically with decreasing contact separation. This is captured well by a simple hydrodynamic model including wall slip. We predict with this model that pressure changes for larger length scales and realistic frictional conditions can easily reach cavitation thresholds and significantly change the load-bearing capacity of a contact. Cavitation may therefore be the norm, not the exception, under boundary lubrication conditions.

  • Boundary lubrication
  • confined fluid
  • wall slip
  • surface heterogeneity
  • cavitation
  • reynolds equation
  • molecular Dynamics
  • tribology


Many surfaces are heterogeneous. They exhibit natural local variations in surface chemistry (13) and roughness (4), and they can be artificially patterned to the desired interfacial chemistry (5, 6) and topography (7, 8). These factors affect wettability (9) and the amount of friction between solids (3) or between wall and fluid during flow (10). Wall friction is usually quantified by the Navier slip length b (11), the fictive distance from the wall where the flow profile becomes zero when extrapolated linearly. Wall slip can be neglected when the channel height h is much larger than b. This is the case in most macroscopic flow geometries where no slip is assumed to take place at the wall (12). In confined geometries, such as those encountered in nanometer-thin lubricants in sliding contacts or in boundary lubrication (13), h ~ b and slip can no longer be ignored. Slip lengths on the order of nanometers and longer (10, 14, 15) are observed on surfaces with weak interactions and low commensurability with the fluid (1620). Slip reduces friction in lubricated devices (21, 22) and can be used to control the flow of fluid mixtures (23). It thus can be exploited in engineering applications (24) by tailoring the surface into slipping and nonslipping domains (25, 26).

Here, we use molecular dynamics (MD) simulations (see Materials and Methods) and analytical models to show that such local variations in surface slip can give rise to large pressure excursions within a lubricated contact under shearing. We present simple scaling expressions that capture the magnitude of the pressure oscillations as a function of geometry, viscosity, and slip length for regular and random slip patterns. The amplitude of these pressure variations can be orders of magnitude larger than the applied pressure and therefore have a profound impact on flow.


We first simulate the simple shearing of a confined lubricant over alternating patches with small and large slip lengths. Our MD system is periodic and has a length L of 95.94 nm and a width of 6.92 nm in the x and y directions, respectively, with two surfaces separated at h = 5 nm along z (see Fig. 1). We use linear pentane (C5H12) as the reference fluid and change the chain length to vary fluid viscosity (see Materials and Methods). The upper surface has no slip, whereas the lower one contains one slipping and one nonslipping domain of equal length λ = L/2 oriented transversely to the shearing direction. Shearing is performed under conditions representative of a highly loaded lubricated contact. An external pressure Pext = 250 MPa is applied to the upper wall, which is displaced at a constant velocity of u2 = 30 m/s.

Fig. 1 MD simulation setup.

The fluid is confined between two surfaces. There is no slip on the top surface and the bottom surface is patterned into slipping and sticking domains.

We also carry out simulations of smaller homogeneous systems without local variations in wall-fluid interaction to assess slip behavior on homogeneous surfaces and compute the slip length. In this case, the system is 6.0 nm long and its lower surface is either slipping or sticking. In agreement with many previous studies on this type of geometry (16, 18), we obtain linear (Couette) velocity profiles u(z) across the film thickness z (red points in Fig. 2, A and B). Slip causes partial plug flow (Fig. 2B), which increases the flow rate relative to the configuration with two nonslipping walls (Fig. 2A).

Fig. 2 Profiles of dimensionless velocity Embedded Image across the normalized film thickness Embedded Image for the heterogeneous system.

(A to C) Sticking domain (A), slipping domain (B), and corresponding normalized slip length Embedded Image (C) along the shearing direction Embedded Image .

The flow profile in the heterogeneous system is not linear despite the simple shear geometry. The blue points in Fig. 2 (A and B) show that the fluid indeed sticks and slips in the respective domains, but the velocity profile curves toward higher velocity in the sticking region and toward lower velocity in the slipping region. Because mass is conserved, the global flow rate through the channel needs to be constant. Slip increases mass transport near the interface, which must be compensated by the lower mass transport in the bulk of the fluid flowing over the slipping region and the enhanced mass transport in the fluid flowing over the sticking region.

To test whether Navier-type slip boundary conditions apply to the heterogeneous system, we quantify the local slip length for each of the 192 bins from the local velocity profiles as shown in Fig. 2 (A and B). The local slip length is given by b(x)=u(x,z)/(∂u(x,z)/∂z)|z=0 , where u is the fluid velocity and z is the position within the channel. Local variations in slip are shown in Fig. 2C and are large. This can be attributed to the large error introduced by extrapolating the velocity profile. Nevertheless, we find an approximately constant value on each domain: in the nonslip region, we find b = 0, whereas on the slipping side, b ≈ 2.3h, consistent with calculations from the reference system with homogeneous walls (shown by the red line in Fig. 2C).

Slip is hence a local property. This is because the geometrical expression for the slip length stems from the fluid viscosity η and the tangential momentum transfer at the wall-fluid interface (11, 27), which are local properties independent of the neighboring domains. The transition zone from no slip to slip is located at the domain boundaries within the slipping domain (Fig. 2C). Its length is ~ 4 nm, which is approximately 10 times the chain length of pentane. It does not change significantly for other short n-alkanes or for other domain lengths λ used in this work and may thus be considered negligible compared to typical pattern sizes (on the order of micrometers) that can be realized experimentally.

The velocity profiles shown in Fig. 2 (A and B) are fully captured by a simple hydrodynamic model that assumes a constant slip length on each of the surface domains. The starting point for this model is the classical Navier-Stokes equation for fluid flow along the shearing (x) direction. Because the estimated Reynolds number for the chosen system Re = 0.22 is low and the thin-film approximation h/λ ≪ 1 applies, the expression can be simplified for laminar flow and negligible inertia effects (12). In steady state, this givesEmbedded Image(1)where we have normalized the fluid velocity by the driving velocity u2, Embedded Image; all lengths by λ, Embedded Image; and all heights by h, Embedded Image. This leads to a normalized pressure Embedded Image, where P0 = λτ0/h and τ0 = ηu2/h is the shear stress necessary to drive the Couette flow of a Newtonian fluid of shear viscosity η. Integrating Eq. (1) twice in the z direction with a no-slip condition and velocity u2 on the upper wall and a Navier slip boundary with slip length Embedded Image on the lower surface gives the velocity profileEmbedded Image(2)The first summand in Eq. (2) is of linear order and hence Couette-like, whereas the second summand is of quadratic order and Poiseuille-like.

The pressure gradient Embedded Image in Eq. (2), which is the origin of the Poiseuille contribution, is given by the well-known Reynolds equation, obtained from combining Eqs. (1) and (2) with the continuity equation (12)Embedded Image(3)Equation (3) accounts for slip on the lower surface and is valid for incompressible, isoviscous fluids under isothermal and steady-state conditions. It can be solved directly by integration (see discussion S1). For the geometry used in our MD calculations, Embedded Image for Embedded Image and Embedded Image for Embedded Image. Accounting for the periodicity of the system along the shearing direction Embedded Image, the pressure gradient isEmbedded Image(4)with Embedded Image. This factor, valid for two grains of the same length λ, varies between κ = 0 for Embedded Image (homogeneous limit) and κ = 1 for Embedded Image.

The black solid lines in Figs. 2 and 3 show that this simple hydrodynamic model perfectly reproduces MD results despite the severe confinement of the fluid. The curved velocity profiles in Fig. 2 (A and B) are captured well by Eq. (2), confirming that fluid flow even down to a channel height of a few nanometers for molecules that are ~ 1 nm in size can be described by continuum hydrodynamics (28). The pressure along the channel is obtained trivially by integrating Eq. (4) once. This gives a pressure drop of δP = ±6κηu2 λ/(5h2) over each of the domains. Figure 3 shows that the total pressure Embedded Image indeed drops linearly on the sticking domain, followed by an increase in the slipping region. Pressure excursions become larger with the applied shear rate through either higher wall velocity u2 or smaller film thickness h, with the viscosity η of the fluid, and with the system length L. For example, longer n-alkane chains lead to higher fluid viscosity η and factor κ, which gives larger pressure gradients.

Fig. 3 Dimensionless pressure distribution along the shearing direction.

Points indicate MD results; the line is the result from the hydrodynamic model. Reference conditions are u2 = 30 m/s , h = 5 nm , L = 95.91 nm , and pentane viscosity η = 0.50 mPa•s at an external applied pressure Pext = 250 MPa . The legend indicates variations of these values. (Inset) Probability distribution of the dimensionless pressure in the presence of N randomly distributed slipping and sticking domains along the shearing direction. Points are obtained from MD simulations with L = 191.8 nm and for six realizations each of N = 24 (λ = 8 nm), N = 48 (λ = 4 nm) , and N = 192 (λ = 1 nm) . Scaling with Embedded Image collapses MD data to a single Gaussian distribution with standard deviation Embedded Image.

Changing the length of the system shows that dP/dx is independent on the domain length; thus, the overall pressure variations along the shearing directions scale linearly with the periodicity length L. This leads to vanishing pressure excursions when L is very small, that is, when it is comparable to the chain length of the fluid molecules. Additionally, the heterogeneous slip/no-slip character of the lower surface disappears if the chain length becomes larger than the domain width for regular patterns (28), which also contributes to vanishing δP . Conversely, large pressure variations, of the order of magnitude of the external pressure, can be reached when L is increased. Figure 3 shows that results of a wide variety of MD simulations with varying u2 , h, η, and L all collapse on the same curve when the pressure is expressed as Embedded Image. The same scaling relationship is obtained in numerical solutions of the Reynolds equation for two-dimensional patterns (see fig. S1). The pressure decreases and increases almost linearly along the shearing direction on the sticking and slipping regions, respectively, with maximum excursions amounting to 66 to 80% of the analytical expression for the pressure drop δP.

Many surfaces may not be artificially patterned in regular slipping and nonslipping domains but exhibit intrinsic, random variations in slip length. Examples are a metal surface that exposes grains of random orientation to the fluid or a diamond surface with mixed hydrogen and oxygen terminations. A simple model for a realization of such a configuration is a channel with a random sequence of domains with either b = 0 or b = bs. Assuming that the total length of the slipping and nonslipping domains is the same, the normalized pressure gradient within each of the domains (of length λ) is then given by Eq. (4) with a random sign, Embedded Image. The resulting pressure Embedded Image in our periodic MD simulations is a random walk in channel position Embedded Image. For a periodic system with N = L/λ domains, the pressure distribution is Gaussian and has a standard deviation of Embedded Image (see discussion S2). Pressure excursions thus scale as Embedded Image, which is in agreement with the MD results shown in the inset of Fig. 3. This estimate works remarkably well even for very short patterns of atomic-scale length, λ ~ 1 nm, where the thin-film approximation h ≪ λ must break down. Such patterns could be realized, for instance, by patches of oleophilic and oleophobic terminating groups.

We therefore obtain similar scaling relationships for regular and random patches. ΔP increases with pattern size, fluid viscosity, and surface speed, or by lowering film thickness. In some cases, the pressure drop may exceed the externally applied pressure Pext, leading to localized spatial regions where P decreases below ambient. This favors the occurrence of dewetting, bubble formation, and dynamic cavitation. We test this in MD calculations on a system that features regular slipping and sticking domains at low external pressure. Figure 4 shows that the pressure excursions are indeed larger than Pext. On the nonslipping domain, a bubble forms and the fluid detaches from the lower surface in a region of negative pressure (see movie S1). This creates hydrodynamic lift under shearing, and the surface separation increases by 10% compared to the stationary case with zero wall velocity. We also solve the Reynolds equation (Eq. (3)) on the noncavitated domain where the pressure calculated from MD is positive. Specifically, we use a no-slip boundary condition on λ1, a slip boundary condition on λ2, and a zero-pressure condition at the crossover from both regions to λcav (see Fig. 4 and discussion S1). This simple cavitation model almost exactly reproduces the pressure in the noncavitating regions.

Fig. 4 Snapshot of a cavitation bubble under shearing and corresponding dimensionless pressure distribution along the shearing direction.

The hydrodynamic model is applied to the noncavitated domains λ1 and λ2, whose extents are obtained from MD results with the condition P > 0. System parameters are u2 = 10 m/s, h = 5.5 nm, L = 143.9 nm , and Pext = 1 MPa. At this pressure, the fluid (n-decane) has viscosity η = 0.30 mPa•s .


The pressure variations described here can be obtained experimentally when surfaces have very low roughness or are treated with oleophobic coatings in only limited regions of the whole contact area. Data for slipping walls at ambient pressure are available for a hexadecane lubricant on photoresist-coated glass surfaces, where b = 25 nm (14), and for glycerol confined between gold surfaces with chemisorbed thiol layers, where b = 40 nm is found (15). For a film thickness of 50 nm and a low shearing speed of 10–4 m/s, the continuum formulation predicts pressure gradients of the order of 108 to 1010 Pa/m. In the presence of surface patterns with a length of ~ 10 μm, which can be easily realized in experiments, this leads to pressure variations of the order of 10 to 100 kPa.

In addition to the previous low-pressure examples, slip has also been reported in highly loaded contacts. An estimate of slip length Embedded Image for diamond-like carbon (DLC) lubricated by polyalphaolefin with h = 100 nm and u2 = 0.4 m/s was given by Kalin et al. (29). In the presence of bare steel (sticking) and DLC-coated (slipping) domains, and using an oil viscosity η = 70 Pa·s at a peak contact pressure of 1 GPa (29), one obtains dP/dx ~ 1015 Pa/m, which is several orders of magnitude larger than the above low-viscosity examples.

A pressure below ambient should be reached for patterns of ~ 1 μm in all aforementioned cases. Thus, cavitation phenomena such as dewetting or bubble formation should then occur, which can significantly enhance the load-bearing capacity of such contacts (30). This opens exciting perspectives for exploiting the properties of slipping surfaces in engineering applications and has important implications for surfaces with natural heterogeneity in slip. The scaling relationships for pressure presented here can be used to estimate the onset of surface-induced cavitation. We note that even on surfaces with naturally occurring random variations in slip, cavitation of the form shown in Fig. 4 can ensue. During squeeze-out of a contact, the channel height h drops continuously and our scaling relation predicts a divergence ∝ h–2 of the amplitude of pressure fluctuations as the channel height approaches zero. When transitioning from hydrodynamic to boundary lubrication during sliding, cavitation may therefore be unavoidable in many situations and could control friction and wear of the contact. This picture is fundamentally different from our current understanding of squeeze-out of a contact at rest that involves fluid layering (31) and liquid-solid transitions (32).

In summary, we found a Couette-like and a pressure-driven, Poiseuille-like contribution to shear-driven flow between two surfaces with heterogeneous slip lengths. The Poiseuille contribution emerges because a constant flow rate along the shearing direction must be sustained. The amplitude of the overall pressure excursion increases with the applied shear rate, fluid viscosity, and system length, and decreases with channel height for regular and random arrangements of slipping and sticking domains. This can lead to dramatic pressure changes under normal tribological conditions and the onset of cavitation in the contact area.


All lubricants considered in this study are short n-alkanes. We varied their chain lengths from 3 to 10 carbon groups to study the influence of fluid viscosity. Intramolecular and intermolecular interactions are described by a United Atom TraPPE force field (33) with a cutoff distance of 1 nm.

The walls are composed of single rigid layers of a (111) surface of a face-centered cubic crystal with a gold lattice constant of 4.08 Å. Wall-fluid interactions are of Lennard-Jones nature with σ = 2.65 Å for the wall atoms. Different slip behaviors can be obtained by simply changing the well depths of the wall-fluid potential energy (16). We thus set this parameter for the slipping and sticking regions to ϵwf,slipping = 0.1ϵf and ϵwf,sticking = 5ϵf, respectively, with ϵf being the average interaction energy within the fluid.

MD simulations were performed with a time step of 1 fs and ran for 100 ns in the case of the heterogeneous system. After reaching steady-state flow and a stable film thickness in approximately 1 ns, the local shear response was obtained through statistical averaging over 192 bins and 99 ns of simulation time in the shearing direction. In the case of smaller systems with homogeneous surfaces, simulations ran for 10 ns. The fluid velocity u(z) was obtained by averaging over the whole domain and 9 ns of simulation time. Slip length b and fluid viscosity η = τ/(∂u/∂z) can then be calculated directly from this velocity profile. Here, τ is the shear stress measured at the walls.

In all simulations, the temperature of the fluid was kept constant at 303 K by means of a Nosé-Hoover thermostat applied in the y direction only, that is, perpendicular to the shearing direction.


Supplementary material for this article is available at

Discussion S1. Analytical solution of the Reynolds equation with slip and a cavitated zone.

Discussion S2. Pressure excursion for random alternations of slip.

Fig. S1. Pressure excursion for two-dimensional slip patterns.

Movie S1. Formation of a cavitation bubble under shearing.

Reference (34)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Computer time at the Jülich Supercomputing Centre (project hfr13) is gratefully acknowledged. We thank the Large-Scale Data Facility at the Karlsruhe Institute of Technology for providing storage space. Funding: This work is based on research partially funded by the Deutsche Forschungsgemeinschaft (grants PA 2023/2 and GU 367/30). Author contributions: All authors designed the study, discussed results, and contributed to writing the manuscript. D.S. carried out MD simulations and derived the hydrodynamic model. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data are available from the authors upon request.

Stay Connected to Science Advances

Navigate This Article