1 Introduction

Many lubricated machine components (such as gears, bearings, constant velocity joints, and cam/follower systems) contain non-conformal surfaces that can roll and slide together. In these components, the solid surfaces elastically deform, leading to elastohydrodynamic lubrication (EHL) conditions [1]. Here, the solid surfaces are generally separated by thin lubricant films (< 1 μm) subjected to very high pressures (> 1 GPa) and shear rates (> 105 s−1). The behaviour of lubricants under such extreme EHL conditions remains shrouded in mystery, which complicates the development of macroscale models capable of accurately predicting the frictional response that ultimately determines energy loss and component efficiency.

Numerical methods, analytical solutions, and experimental techniques have all advanced significantly over the last few decades [2]. However, the measurement and prediction of EHL friction remains a significant challenge and continues to attract much controversy [3,4,5]. In addition to the aforementioned methods, molecular dynamics (MD) simulations have emerged as tools capable of giving unique insights into the nanoscale behaviour of lubricants, including under EHL conditions [6]. Much of this work has been pioneered by Mark Robbins alongside his collaborators. In the remainder of this article, we give an overview of the state of the art in the use of MD simulations to help answer outstanding questions in EHL. We pay particular attention to the work of Robbins and co-workers, showing how his work has considerably advanced the field.

2 Shear Thinning

One of the key outstanding questions in EHL is the most suitable equation to describe the non-Newtonian response of lubricants at high shear rates. The two phenomenological models at the centre of this active debate are those due to Eyring [7] and Carreau [8]. The Eyring model [7] assumes that flow occurs by stress-augmented thermal activation (SATA) over a single characteristic energy barrier height. It can be used to describe the change in shear stress, \(\tau\), with shear rate, \(\dot{\gamma }\):

$$\tau ={\tau }_{e}{sinh}^{-1}\left(\frac{{\eta }_{P}\dot{\gamma }}{{\tau }_{e}}\right),$$
(1)

where \({\tau }_{e}\) is the Eyring stress, \({\eta }_{P}\) is the Newtonian shear viscosity at pressure \(P\). The Eyring model predicts that, at high \(\dot{\gamma }\) (when \(\tau\) exceeds \({\tau }_{e}\) by a factor of at least two), \(\tau\) increases linearly with \(\mathrm{log}(\dot{\gamma })\). The model has been used to describe a wide range of lubricant rheology data, mostly from tribometer experiments under EHL conditions rather than viscometers [3]. The related Ree–Eyring [9] model, which allows for multiple flow units (via different barrier heights), has also been applied to lubricant shear thinning under EHL conditions [3].

The Carreau model [8] can be used to describe the variation in shear stress, \(\tau\), with shear rate, \(\dot{\gamma }\):

$$\tau =\frac{{\eta }_{P}\dot{\gamma }}{{(1+{\left(\lambda \dot{\gamma }\right)}^{2})}^{(1-n)/2}}$$
(2)

where \(\lambda\) is a time constant and \(n <1\). The Carreau model predicts that, at high \(\dot{\gamma }\), \(\tau\) increases as a power law with \(\dot{\gamma }\). This equation was derived based on the network theory of Lodge [10], but it is now recognised that power-law behaviour can arise from several different physical pictures of flow, including SATA over a broad range of energy barrier heights. One immediate benefit of the Eyring equation over the Carreau equation is its simplicity, i.e. the fact that it contains fewer tuneable parameters [3]. The Carreau–Yasuda [11] equation contains modifications to improve fitting in the transition region between Newtonian and strong shear thinning regions [3].

The Carreau equation contains many other models as limiting cases [12], including the Eyring model [13]. While the Eyring and Carreau models predict similar shear stress (or viscosity) at low shear rates, they diverge at high shear rates [3]. In fact, the shear stresses from the models can differ by several order of magnitude at shear rates (> 107 s−1) encountered in lubricated components inside modern engines [14].

Testing models for shear thinning at such high shear rates and pressures represents a considerable experimental challenge. High-pressure viscometers can reach shear rates of only ~ 104 s−1 at EHL pressures (> 1 GPa) for low-viscosity fluids before shear heating becomes uncontrollable [15]. At similar pressures, tribometer experiments under thin-film EHL conditions can reach much higher shear rates (106 s−1) before shear heating becomes excessive [16], but it has been debated how useful the rheological data from such experiments are to predict EHL friction [3,4,5]. On the other hand, non-equilibrium molecular dynamics (NEMD) simulations are ideally suited to study the rheology of fluids at high shear rates, since thermostats can be used to suppress shear heating. In fact, NEMD simulations suffer from the opposite problem, in that to obtain a reasonable signal-to-noise ratio and a non-equilibrium steady state, high shear rates (usually > 107 s−1) must be used [6]. Consequently, through a combination of experiment and NEMD simulation, our understanding of shear thinning behaviour under EHL conditions has been substantially improved in recent years.

In 2002, Bair et al. [17, 18] conducted the first comparison between the non-linear rheology predicted by NEMD simulations with that obtained experimentally. They used the inexpensive united-atom force field developed by Siepmann et al. [19] to study the viscosity of squalane (C30H62), a model lubricant base oil. The SLLOD [20] equations of motion were applied to a bulk system of squalane molecules and NEMD simulations were conducted at shear rates between 107 and 1011 s−1 and pressures of up to 0.3 GPa. As in most bulk NEMD simulations, a global Nosé–Hoover thermostat was used to prevent shear heating [21]. The experiments were conducted at shear rates between 102 and 104 s−1 and pressures of up to 1.0 GPa. Although the experimental and NEMD data were separated by several orders of magnitude in shear rate, they collapsed onto the same time–temperature superposition (TTS) master curve. This master curve was successfully fit using the Carreau model [8].

Liu et al. [22] also successfully used the Carreau model [8] to describe the shear thinning behaviour of squalane, as well as several types of poly-α-olefin (PAO) molecules, from bulk NEMD simulations at high shear rate (108–1010 s−1). They correlated the onset of shear thinning with changes to the conformation of the molecules, quantified through the radius of gyration. Santak and Conduit [23] used the Carreau model [8] to estimate Newtonian viscosities from NEMD simulations preformed at high shear rate (109–1012 s−1) for 74 alkanes containing 5–10 carbon atoms with weighted least squares regression.

In 2017, Jadhao and Robbins [24] used bulk NEMD simulations to determine the non-equilibrium shear response of squalane over a much wider range of pressures (0.1–1.0 GPa) and shear rates (105–1010 s−1). They used the united-atom force field due to Mondello and Grest [25], which is an extension of the force field developed by Siepmann et al. [19]. This force field was previously shown to reproduce the experimental Newtonian viscosity of squalane within 15% [26]. The Newtonian viscosities obtained from the NEMD simulations were in good agreement with those from previous experiments performed at low shear rates (102–104 s−1) [15, 17].

At low pressures (and high temperatures), Jadhao and Robbins [24] found that shear thinning could be described using a power law, as in the Carreau model [8]. At higher pressures (and low temperatures), where the viscosity increased above ∼1 Pa.s, squalane began to exhibit glassy flow, and shear thinning was increasingly dominated by SATA processes. In this high-viscosity regime, the shear thinning behaviour was more accurately described using the Eyring model [7], as shown in Fig. 1. These observations have important implications in terms of predicting EHL friction since they imply that the Eyring model [7] is more suitable than the Carreau model [8] to describe the non-linear rheology in high-pressure EHL contacts.

Fig. 1
figure 1

Reproduced with permission from Ref. [24]. Copyright American National Academy of Sciences, 2017

Comparison between shear rate dependence of shear stress (a) and viscosity (b) from NEMD simulations and viscometer experiments. Dashed black lines are fits to the Eyring equation [7]. Data at lower pressure show Carreau-like [8] behaviour (fits not shown).

Full size image

This study by Jadhao and Robbins [24] was a significant breakthrough and was highlighted as a key study in our 2018 reviews on both SATA [27] and NEMD simulations [6] in tribology. Moreover, the authors made their large atomic/molecular massively parallel simulator (LAMMPS) [28] input files publically available in the Supporting Information, ensuring that their results were transparent, reproducible, usable by others, and extensible (TRUE) [29].

In 2019, Jadhao and Robbins [30] undertook a more detailed investigation of the transition between Eyring [7] and Carreau [8] shear thinning behaviour under EHL conditions using NEMD simulations. They showed that, in the high-viscosity regime (> 1 Pa.s), squalane viscosity data could be collapsed onto the Eyring model [7] over almost 30 orders of magnitude in shear rate using TTS. Conversely, the Carreau model [8] (with a power-law exponent of 0.5 as proposed by Bair et al. [17]) severely overestimated the viscosity from the NEMD simulations at high shear rate. The comparison is shown in Fig. 2a. For high-viscosity conditions, the response is dominated by SATA rearrangements over the lowest energy barrier, consistent with Eyring theory. For low-viscosity conditions, many different energy barriers contribute to molecular motion, leading to Carreau-like behaviour. It is noteworthy that modifications of the Eyring equation, such as the Ree–Eyring equation [9], can be used to account for multiple barrier heights. However, this increases the number of tuneable parameters and reduces the predictive power of the model [3, 30].

Fig. 2
figure 2

Reproduced with permission from Ref. [30]. Copyright Springer Nature, 2019

a Comparison of Eyring [7] and Carreau [8] models to viscosity data for squalane collapsed using TTS for values of P and T where η is high (> 1 Pa.s). Experimental and NEMD simulation data are, respectively, indicated by filled and open symbols. b Order parameters Sxx and Syy for the systems shown in a scale by their maximum magnitudes at high rates. Results for the order at all state points collapse when plotted against the same reduced rate that collapses viscosities. Note that alignment saturates at a reduced strain rate of about 10, where the viscosity has dropped by less than an order of magnitude.

Full size image

Jadhao and Robbins [30] also studied the prevalence of different mechanisms that have been proposed to control shear thinning behaviour. One commonly invoked mechanism for power-law (e.g. Carreau) shear thinning is a change in some order parameter, more specifically alignment of chain-like molecules along the flow direction [4, 31]. This could perhaps lower the viscosity by decreasing the rate of collisions between molecules, as well as the associated momentum transfer and dissipation. In contrast, shear thinning in the Eyring model [7] only involves a change in the rate of molecular rearrangements through activated motion over fixed barriers. Jadhao and Robbins [30] showed that the alignment of squalane molecules in the NEMD simulations (Fig. 2b) saturates after the viscosity (Fig. 2a) has dropped by only about a factor of three due to shear thinning. On the other hand, the Eyring model [7] describes a reduction of six or more decades in viscosity resulting from shear thinning.

Using a binary atomic Lennard–Jones fluid [32], Jadhao and Robbins [30] also showed that the transition from Carreau to Eyring is generic, rather than specific to chain-like squalane molecules. Lin and Kedzierski [33] have recently shown that the Carreau to Eyring transition also occurs in a typical polyol ester lubricant, pentaerythritol tetrahexanoate (PEC6), using bulk NEMD simulations. They used the long-chain optimised potential for liquid simulations (L-OPLS) all-atom force field [34]. Prentice et al. [35] successfully employed the Eyring model to estimate the Newtonian viscosities from NEMD simulations of mixtures of hydrogenated 1-decene trimer and 1-decene tetramer as a model for PAO.

It is worth noting that some earlier work of Robbins and collaborators on glassy systems provided an initial guide to the existence of a crossover from power-law (Carreau) to log-rate (Eyring) behaviour driven by increasing Newtonian viscosity. For example, Rottler and Robbins [36] observed that, for bi-disperse Lennard–Jones fluids, a transition from Newtonian behaviour to a logarithmic shear rate dependence of shear stress (consistent with Eyring theory [7]) occurred as the temperature was lowered. Furthermore, Baljon and Robbins [37] studied the Kremer–Grest [38] model of short-chain molecules with 16 beads and observed a change from power-law to Eyring-like behaviour with decreasing temperature. Thompson et al. [39] also found similar changes in the behaviour of confined chains as the Newtonian viscosity was increased, due to decreasing film thickness.

Most of the bulk NEMD simulations performed by Jadhao and Robbins [24, 30] were performed at constant volume. This differs from viscometers and tribometer experiments that are performed at constant pressure. Jadhao and Robbins [24, 30] estimated the viscosity at each target pressure by interpolating between the constant volume results. Results from previous NEMD simulations have shown non-negligible differences between viscosities obtained in the constant-volume and constant-pressure ensembles [40]. Jadhao and Robbins [30] performed a subset of NEMD simulations at constant pressure (in the normal and shear directions) and found that this resulted in similar shear stress (or viscosity) up to shear rates of approximately 106 s−1, but somewhat lower shear stress at high shear rates.

Some potentially significant assumptions have been made in the bulk NEMD simulations that have been performed to date of lubricants under EHL conditions. The first assumption is that linear flow profiles develop between the sliding surfaces, which may not be the case in experiments (as will be discussed in detail in a subsequent section). This could be avoided by using a configurational thermostat [41, 42]; however, these are not yet implemented in open-source MD programs such as LAMMPS [28]. Another potential issue is that the NEMD equations of motion implemented in LAMMPS [28] are the so-called atomic SLLOD [20], rather than the molecular SLLOD [43]. The use of atomic SLLOD [20] for molecular systems results in enhanced orientational ordering and non-zero values for the antisymmetric part of the molecular stress tensor, indicating that a torque is applied to the molecules [44]. For a more detailed discussion of the potential importance of molecular SLLOD [44] in bulk NEMD simulations, see Ref. [45]. Finally, Bair has also recently questioned the ability of NEMD thermostats to maintain isothermal conditions at the high shear rates [46] employed by Jadhao and Robbins [24, 30]. The energy released in these NEMD simulations is certainly considerable; however, the thermostatting methods applied have been shown to maintain isothermal conditions up to even higher strain rates (1012 s−1) [47] than those considered by Jadhao and Robbins [24, 30].

In addition to bulk NEMD simulations, confined NEMD simulations have also shown good agreement with shear stress extrapolations from tribometer experiments conducted at lower shear rates for several molecular fluids, such as model lubricants squalane [48, 49] and diethylhexyl sebacate [48], traction fluids DM2H and Santotrac 2000 [48], as well as benzyl benzoate [49]. These studies were all performed at constant pressure (consistent with experiments) and used the L-OPLS force field [34]. Although, they showed good agreement with experimental trends, these studies did not directly compare shear thinning models.

3 High-Pressure Newtonian Viscosity

Accurate measurements of high-pressure viscosity are crucial for macroscale predictions of EHL film thickness and friction [3,4,5]. The viscosity in the inlet, where pressure is approximately 100 MPa, controls EHL film thickness, while viscosity in the central region of the contact, where pressure often exceeds 1 GPa, governs EHL friction. Through recent advances in computational performance, highly parallelised simulation software, and force field accuracy, MD simulations are now ideally suited to predicting the high-pressure viscosity of lubricant-like molecules [6].

In 2001, McCabe et al. [50] used NEMD simulations to study the high-pressure viscosity of 9-octylheptadecane. They used the united-atom force field developed by Mondello and Grest [25] to study the viscosity up to 1 GPa at high temperature (373 K). The Newtonian viscosity values from the simulations were significantly under-predicted (30–70%) compared to experiment; however, the pressure–viscosity coefficient, α, was within 10% of the experimental value. In calculating α, the authors assumed an exponential increase in viscosity over the entire pressure range studied (0.1–1.0 GPa), as predicted by the Barus equation [51].

More recently, Liu et al. [52] used the same force field to investigate the pressure dependency of viscosity of hydrogenated 1-decene trimer, as a model for PAO4. They found that, as previously observed experimentally, the pressure increase was slower than exponential up to pressures of 1 GPa, which is consistent with the McEwen model [53]. Liu et al. [52] observed a similar viscosity under-prediction as McCabe et al. [50] compared to experiment, and also noted that the pressure–viscosity coefficients they obtained under ambient (α0) and high-pressure (α*) conditions were in close agreement with experiment [15].

The 10th Industrial Fluid Properties Simulation Challenge (IFPSC) [54] involved determining the Newtonian viscosity of 2,2,4-trimethylhexane at 293 K and pressures between 0.1 and 1000 MPa using molecular simulations. New experimental data determined by Bair using a high-pressure viscometer [55] were kept secret from participants until they had all submitted their entries. First place in this competition was awarded to Galvani Cunha and Robbins [56], who used a similar approach to Jadhao and Robbins [24, 30]. At low pressure, the Newtonian viscosity could be obtained directly from the NEMD simulations [24, 30, 56], while at higher pressure this was estimated using extrapolations to lower shear rates with the Eyring equation [7]. Instead of the united-atom Modello and Grest force field [25], Galvani Cunha and Robbins [56] used the all-atom adaptive intermolecular reactive bond-order potential-Morse (AIREBO-M) [57] force field, developed by O’Connor and Robbins. Although they increase the computational cost by approximately an order of magnitude, all-atom force fields have previously been shown to be much more accurate than united-atom ones for the viscosity prediction of long-chain hydrocarbons [58]. Unlike the original AIREBO force field [59], AIREBO-M [57] uses a Morse potential, rather than the standard 12–6 Lennard–Jones potential for the non-bonded interactions. The exponent for the repulsive component in the 12–6 Lennard–Jones potential has never been confirmed from first principles or determined experimentally, and it is used mostly out of convention and convenience [60]. O’Connor and Robbins suggested that the unphysically divergent power-law repulsion causes AIREBO with 12-6 Lennard–Jones to fail when applied to systems at high pressure [57].

Galvani Cunha and Robbins [56] obtained 2,2,4-trimethylhexane viscosity values that deviated by less than 40% from the experimental values over the entire pressure range [55], as shown in Fig. 3a. The agreement was excellent in the 100–700 MPa range, but deviations were larger at lower and higher pressure. The viscosity increased at a slower than exponential rate at low pressures, reached an inflection point, and then increased at a faster-than-exponential rate. The authors fitted their pressure–viscosity data using a hybrid model of the McEwen model [53] and Paluch model [64], as proposed by Bair [62]. The McEwen model [53] is well suited to represent the slower than exponential increase in viscosity with pressure that is frequently observed at low pressure. It is equivalent to the Tait model [63] for ambient reference pressures. On the other hand, the Paluch model [64] can be used to fit the faster-than-exponential increase in viscosity with pressure that is commonly observed at high pressure. The Paluch model [64] assumes the existence of a finite pressure at which relaxation times diverge at a given temperature, which is analogous to the temperature where the viscosity diverges in the popular Vogel–Fulcher–Tammann (VFT) equation [65]. The hybrid model can be used to fit both types of pressure–viscosity response, including the inflection point that is often observed at intermediate pressure [62].

Fig. 3
figure 3

Reproduced with permission from Ref. [56] (a) and Ref. [61] (b). Both copyright Elsevier, 2019

Pressure–viscosity relation for 2,2,4-trimethylhexane from experiment [55] and MD simulations by Galvani Cunha and Robbins [56] (a) and Kondratyuk and Pisarev [61] (b). Black dotted line in (a) is a fit to the MD data with the hybrid model [62] and red dotted line in (b) is a fit to the MD data with a Tait-like model [63].

Full size image

Kondratyuk and Pisarev [61] were awarded second place in the 10th IFPSC. They employed the all-atom condensed-phase optimised molecular potentials for atomistic simulation studies (COMPASS) force field [66]. Instead of using NEMD simulations, they applied the Green–Kubo method [67, 68] in equilibrium molecular dynamics (EMD) simulations. Although EMD is a more direct route to the Newtonian viscosity, it is typically limited to low-viscosity systems (< 50 mPa.s) [69]. Systems with larger viscosities require longer integration times, which can lead to error accumulation. Kondratyuk and Pisarev [61] employed the time decomposition method (TDM) [70], which uses the time-averaged running integral from multiple independent, short MD trajectories to calculate the viscosity. Using the TDM, the authors were able to accurately predict the viscosity of 2,2,4-trimethylhexane up to 600 MPa. At higher pressure, the viscosity was significantly under-predicted compared to experiment [55], by up to around 70% at 1 GPa. Since they did not observe an inflection point or faster-than-exponential viscosity increase with pressure in the simulations, the authors fitted the data with a Tait-like model [63], which is equivalent to the McEwen model [53].

Messerly et al. [71] were awarded third place in the 10th IFPSC, using EMD simulations and Green–Kubo with the TDM [70]. They used a version of the united-atom transferable potentials for phase equilibria (TraPPE) force field [72] for branched alkanes with a 16-6 Mie potential [73]. Heyes [74] showed that varying the steepness of purely repulsive potentials has a dramatic effect on transport properties, while Gordon [75] showed that an n-6 Mie potential (where n = 14–20) could accurately describe the high-pressure viscosity of n-alkanes. The 16-6 Mie TraPPE force field [73] gave near-quantitative agreement with experimental viscosity data [55] (to within the combined uncertainties) from atmospheric pressure to 600 MPa. The force field also reproduced the inflection point and faster-than-exponential viscosity increase at high pressure that was observed experimentally [71]. However, the viscosity at high pressure was dramatically over-predicted (> 100%) compared to experiment at higher pressure. This suggests that the overly repulsive Mie 16-6 potential limits the pressure range over which the force field is reliable [71].

It is interesting to note that none of the entries used force fields with the standard 12-6 form of the Lennard–Jones potential for the non-bonded interactions. This is despite the recent observation that experimental data for a large number of molecular systems at high pressure may be described by a common repulsive potential \(U(r)\propto {r}^{-12}\) [60]. Thus, one area certainly worthy of further study is the identification of the most suitable potential to represent the non-bonded interactions for MD simulations of lubricant-sized hydrocarbons (C10–C40 [76]) at high pressures. This could perhaps be facilitated by re-parameterizing MD force fields using first principles calculations with accurate representations of van der Waals interactions [77].

Another interesting observation from the 10th IFPSC is that the majority of the successful entries used an all-atom force field. While accurate viscosities can be obtained for some hydrocarbons (e.g. squalane that contains multiple short methyl branches [78]) the viscosity of most hydrocarbons are severely under-predicted using united-atom force fields [58]. Although Messerly et al. [71] demonstrated that this can be overcome to some extent by using a more repulsive potential, this leads to large over-predictions of viscosity at high pressure. Thus, allowance for the atomic-scale roughness imposed by hydrogen atoms for systems of interacting molecules in MD simulations seems critical to accurately predict viscosity under high pressures. In most cases, the additional accuracy of all-atom force fields justifies the significantly associated computational cost.

Several subsequent studies have used MD simulations to predict the high-pressure viscosity of other systems. For example, Lu et al. [79] investigated the viscosity of traction fluid molecules containing linked cyclohexane rings using NEMD simulations and the united-atom TraPPE force field [80]. They found that the pressure–viscosity coefficient from NEMD was underestimated for these molecules by between 20 and 50% compared to experiment [79]. Using NEMD and the all-atom L-OPLS force field [34], Lin and Kedzierski [33] found that the McEwen model [53] accurately described the experimental viscosity of PEC6 up to 1 GPa.

Commercial lubricants are complicated mixtures rather than single well-defined molecules. Kondratyuk et al. [81] recently used the all-atom COMPASS [66] force field and Green–Kubo with the TDM [70] to study the high-pressure viscosity of mixtures. They studied the viscosity and density of methylcyclohexane, 1-methylnaphthalene, and their binary mixtures at 323 K and pressures of up to 500 MPa. The simulation results were in quantitative agreement with the previous experiments available up to 100 MPa for both pure compounds (200 MPa for 1-methylnaphthalene) and their mixtures. For 1-methylnaphthalene, the viscosity initially increased at a slower than exponential rate with pressure before it reached an inflection point and then increased faster than exponential. The inflection point (300 MPa) occurred at a pressure slightly below the one at which 1-methylnaphthalene is expected to enter the supercooled phase (400 MPa). For methylcyclohexane, the increase in viscosity with pressure was slower than exponential over the entire pressure range studied. The methylcyclohexane data were fitted using the McEwen model [53], while the 1-methylnaphthalene data were fitted using Bair’s hybrid model [62]. The binary mixtures showed intermediate pressure–viscosity responses between the two pure cases. The mixtures became increasingly non-ideal and the pressure was increased. The mixture viscosities were successfully predicted using the Grunberg–Nissan mixing rules [82] over the entire pressure range.

Despite significant recent progress, MD simulations of the high-pressure viscosity of lubricant-sized molecules (C10–C40 [76]) are still a computationally demanding and time-consuming task. An interesting recent development by Kajita et al. [83] has been to dramatically accelerate (by three orders of magnitude) the evaluation of viscosity using a femtosecond stress-tensor correlation. Using this method, they screened the viscosity–temperature behaviour of 55,000 lubricant-sized molecules. This enabled a Monte-Carlo tree search algorithm to be used to design molecules with a high-viscosity index. Kajita et al. [83] justified the use of femtosecond stress-tensor correlations through the shoving model developed by Dyre et al. [84]. Usually, nanosecond stress-tensor correlations are required in the Green–Kubo method to obtain accurate viscosities [69]. The shoving model was developed for glasses [84], and it is not yet clear whether the model is suitable to describe the viscosity of lubricant-sized molecules in the liquid phase. Some aspects of the model have been experimentally validated for some lubricant-like molecules, such as 5-phenyl-4-ether (5P4E) [85]. While this synthetic polyether has been used as a model lubricant [86], it has a much lower glass transition temperature than conventional lubricants [87]. Therefore, 5P4E is a more natural application of the shoving model than lubricant-sized molecules [84].

4 Flow and Phase Behaviour

A major assumption of most macroscale EHL models is that a linear velocity profile develops in the fluid between the sliding surfaces. The same assumption of planar Couette flow is made in bulk NEMD simulations that use the SLLOD [20] equations of motion alongside standard thermostats (e.g. Nosé-Hoover [21]). However, recent experiments have highlighted that the through-thickness velocity profiles in some systems can show significant deviations from linearity under EHL conditions. Some experiments have shown that the central region can become solid-like and flow more slowly than the liquid-like region close to the surfaces, which is often referred to as plug slip or plug flow [88,89,90]. Experiments for different systems and conditions have shown that the fluid close to the surfaces can become solid-like and flow more slowly than the liquid-like central region, which is known as central localisation or shear banding [91,92,93]. The different types of through-thickness flow profiles possible in experiments under EHL conditions are shown in Fig. 4.

Fig. 4
figure 4

Reproduced from Ref. [6]

a Schematic showing possible lubricant flow behaviour under EHL conditions for two surfaces separated by distance hc, moving at velocities of Ua and Ub. Snapshots showing examples of the flow behaviour from NEMD simulations [48].

Full size image

Confined NEMD simulations in which the wall atoms (rather than the sheared fluid molecules) are thermostatted are more representative of tribometer experiments [94]. In these simulations, no assumptions are made regarding the flow profile that develops in the fluid between the sliding surfaces [95]. This enables deviations from Couette flow behaviour such as boundary slip [96] and shear banding [97] to be accurately captured. Moreover, temperature gradients are allowed to form within the fluid film, as in experiments [94]. The amount of shear heating can be controlled through the coupling time of the thermostat to the wall atoms [98].

In 1990, Thompson and Robbins [99] performed the first systematic study of the structure and flow of fluids near solid boundaries. They performed NEMD simulations of atomic Lennard–Jones fluids confined and sheared between solid surfaces. They observed a wide range in flow conditions, including slip, no-slip, and locking. They showed that the degree of slip was directly related to the amount of structure induced in the fluid by the solid walls. For weak wall–fluid interactions, there was little ordering and slip was observed. For large interactions, substantial epitaxial ordering was induced and the first few fluid layers became locked to the wall, leading to slip within the fluid (i.e. shear banding). Subsequent confined NEMD simulations of Lennard–Jones fluids by Butler and Harrowell [100], binary Lennard–Jones fluids by Varnik et al. [101], and short polymers by Thompson et al. [102] and Rottler and Robbins [36] have also shown various forms of shear banding. Washizu and Ohmori [103] performed the first confined NEMD simulations using a molecular fluid (n-hexane) with film thicknesses (10 nm) and shear rates (106 s−1) that were representative of EHL tribometer experiments.

Since 2012, Dini and co-workers have studied the friction and flow behaviour of Lennard–Jones and binary Lennard–Jones fluids under EHL conditions [104,105,106]. They observed that, above a certain pressure, the flow profiles transitioned from Couette-like to shear banding [104]. Non-equilibrium phase diagrams were constructed for Lennard–Jones fluids at different pressures and shear rates [105]. At low shear rate and high pressure, plug slip occurred where shear was localised close to the walls, while at higher shear rates, central localisation occurred where shear mainly occurred in the middle of the fluid film [104,105,106]. At very high shear rates, the central liquid region widened until eventually a Couette profile was recovered [104,105,106]. Unusual friction behaviour was noted in the central localisation and particularly plug slip regions [106]. Martinie and Vergne [107] verified the feasibility of mapping the lubricant phase diagram derived from the NEMD simulations to the results obtained from previous experiments. From this mapping, they suggested that the physical mechanism that initiates ‘limiting shear stress’ (LSS) behaviour in EHL contacts is shear banding [107]. The LSS describes the shear stress at which the friction coefficient becomes insensitive to pressure and shear rate. The LSS concept is commonly invoked when explaining and predicting EHL friction curves.

In recent years, a number of studies have investigated the potential links between friction and flow behaviour under EHL conditions for molecular fluids. Washizu et al. [108] observed shear banding in NEMD simulations on n-hexane confined between gold slabs at high pressure (0.1–8.0 GPa). They also observed LSS behaviour at the transition from Couette flow to CL. Ewen et al. [48] performed tribometer experiments and NEMD simulations of a range of lubricant (squalane, diethylhexyl sebacate) and traction fluid (Santotrac 2000, DM2H) molecules at pressures (0.5–2.0 GPa) and shear rates (104–1010 s−1) representative of EHL conditions. On increasing pressure, they observed a transition from Couette flow to CL for the lubricant molecules and from CL to PS for the traction fluids. While the lubricants showed a decreased dependence of the friction coefficient on shear rate at high pressure, they did not reach a LSS. Conversely, the traction fluids displayed a LSS at the same pressure as the CL-PS transition occurred. Porras-Vazquez et al. [49] used tribometer experiments and NEMD simulations to study the relationship between the LSS and flow behaviour for squalane and benzyl benzoate. They found similar results to Ewen et al. [48] on atomically smooth surfaces, but also found similar friction coefficients on rough surfaces, where the flow behaviour remained Couette-like even at high pressure. Based on these observations, they suggested that the LSS was caused by a pressure-induced liquid-glass transition in the lubricant [49]. This suggestion is consistent with the previous NEMD observations, since shear banding is well known to occur in glassy materials [101] and the LSS has been linked to the yield stress of glassy lubricants [109]. For n-alkanes, which crystallise more readily than branched lubricant molecules, LSS behaviour and non-linear flow were observed in NEMD simulations of solid-like C60, but not liquid-like C16 [110].

The glass transition in liquid lubricants has also been studied with MD simulations. Thompson et al. [39] showed that the dynamics of simple chain-like molecules became glassy when confined below around six molecular diameters. Recent MD simulations by Gao and Müser [111] of hydrogenated 1-decene trimer molecules between iron walls suggest that while confinement leads to drastically increased layering, load-bearing capacity, and effective viscosity, this does not always indicate crystallisation or a glass transition. However, under EHL conditions, lubricants are subjected not only to strong confinement (< 1 μm), but also very high pressures (> 1 GPa). As suggested by Robbins and Smith [112], there could be connections between lubricant behaviour under confinement and at high pressure.

The effect of the surface properties can also be important in determining flow under EHL conditions [113]. Ewen et al. [98] performed NEMD simulations to study the friction and flow of Santotrac 2000 and hydrogenated 1-decene trimer between surfaces with different thermal conductivities. Reducing the thermal conductivity of one of the surfaces led to large temperature rises close to that surface, resulting in asymmetric flow profiles and significant reductions in friction [98]. The traction fluid (Santotrac 2000) gave significantly higher friction and temperature rises compared to the lubricant (hydrogenated 1-decene trimer). Savio et al. [114] performed NEMD simulations of short n-alkanes (C3–C10) and showed that both regular and random slip patterns on solid surfaces can cause cavitation, particularly at small contact separations.

Although it is commonly observed in highly confined NEMD simulations with atomically smooth substrates [96, 115], slip directly at the solid–liquid interface (‘boundary slip’ in Fig. 4) is unlikely to be significant for lubricants on the rough, oleophilic steel surfaces usually found inside tribological contacts [116]. Gao et al. [117] showed that, for n-hexadecane molecules confined between gold surfaces, slip at the solid–liquid interface was eliminated in the presence of random nanoscale roughness. Furthermore, Washizu et al. [118] showed that, for n-hexane molecules confined between atomically smooth gold surfaces, the slip length decreased with increasing surface separation and was negligible for realistic EHL film thicknesses (400 nm).

Experimental studies have shown that boundary slip can occur inside EHL contacts when an oleophobic coating is applied to one of the surfaces [119]. Organic friction modifier (OFM) additives, such as fatty acids, are known to form oleophobic coatings on metal surfaces [120]. OFMs are generally believed to reduce friction in the boundary regime by forming strong monolayers that can prevent solid–solid asperity contact [121]. However, Choo and Spikes [122] showed experimentally that OFMs also reduce friction in the low-load hydrodynamic lubrication regime. They postulated that this was due to the OFMs forming oleophobic monolayers that promoted boundary slip [122]. Ewen et al. [123] confirmed this hypothesis with NEMD simulations of n-hexadecane lubricant confined between OFM monolayers adsorbed on iron oxide surfaces. The NEMD simulation results also suggested that larger reductions in friction might be possible in the high-pressure EHL regime, due to larger slip lengths [123]. This hypothesis was recently confirmed experimentally by Kalin and Kus [124], who showed that OFMs cause large reductions in friction compared to pure base oils under EHL conditions.

5 Hybrid Simulations

Our understanding of many problems in tribology, including EHL contacts, would benefit from multiscale modelling techniques [125]. The main limitation of MD simulations is the short time and length scales that can be accessed [6]. While coarse graining could help in this respect [126], parameterisation of coarse-grained models that can accurately describe lubricant rheology has proved a considerable challenge [127]. Macroscale techniques, such as continuum fluid dynamics (CFD), have been successfully applied to a range of EHL problems [128,129,130,131]. However, such techniques require several assumptions with respect to the rheology of the lubricant and flow boundary conditions. One approach to improve the accuracy of continuum models is to perform hierarchical coupling to smaller scales [6]. For example, the viscosity [132, 133] and slip length [134, 135] have been calculated in NEMD simulations and then used as inputs in continuum models.

A more significant challenge, which could transform the modelling of EHL contacts, would be to perform concurrent molecular-continuum coupling. The domain decomposition approach proposed by O’Connell and Thompson [136] defines atomic and continuum domains that are coupled through an overlap region in which the solutions in both domains are consistent. This limits the use of MD to regions where the atomistic scales need to be resolved, while employing a continuum-based solver for the remainder of the domain. This approach has been widely adopted to study a range of multiscale flow problems over the last few decades [137].

Robbins and co-workers applied the domain decomposition approach to perform hybrid molecular-continuum simulations of a range of factors important to EHL contacts including nanoscale roughness [138], cavity flow [139], and solid–liquid heat transfer [140]. Aspects linked to asperity interactions and elastic–plastic surface deformation have also been successfully considered [141,142,143], showing that hybrid coupling could be potentially used not only to capture the fluid behaviour and complex solid/fluid interactions, but also the evolution of solid surfaces.

One major limitation of concurrent coupling techniques is that the simulation time step is limited to the shortest timescale that needs to be modelled at the smallest scale [6]. Despite this, considerable performance enhancements can be achieved compared to modelling the entire system at the atomic level [137]. Moreover, Liu et al. [144] have developed a multi-timescale algorithm, which allows the molecular dynamics during one macro time step to be treated as in quasi-steady state, further decreasing the required computer time.

The majority of published modelling simulation research now utilises highly parallelised, open-source codes such as LAMMPS [28] and OpenFOAM [145]. In recent years, significant efforts have been made to enable hybrid molecular-continuum simulations using these software packages [146, 147] while building on the coupling algorithms and methods [137, 148] developed by Robbins and his collaborators [138,139,140,141,142]. The combined gains in computational efficiency from highly parallelised software and hybrid simulations could enable the study of EHL contacts with atomistic detail close to the surfaces in the near future. These multiscale methods could uncover many aspects of EHL behaviour that remain shrouded in mystery.

6 Summary

MD simulations are becoming integral to our fundamental understanding of lubricant behaviour inside EHL contacts. Both bulk and confined NEMD simulations have become promising tools to study the non-Newtonian rheology of lubricants under EHL conditions. Robbins and co-workers have shown that at relatively low pressure, many energy barriers contribute to molecular motion, leading to a power-law increase in shear stress with shear rate, as predicted with the Carreau equation. At high pressure, where the viscosity is much greater, the response is increasingly dominated by SATA rearrangements over the lowest energy barrier, leading to a logarithmic increase in shear stress with shear rate, consistent with Eyring theory.

MD simulations can now show quantitative agreement with experiments used to measure the high-pressure viscosity of lubricant-like molecules. Both EMD and NEMD simulations have been successfully applied to predict viscosities up to pressures of around 1 GPa. The most accurate predictions have used all-atom force fields that do not employ the standard 12-6 form of the Lennard–Jones potential for the non-bonded interactions. Since commercial lubricants are complicated multicomponent fluids, it will be interesting to investigate the predictive capability of MD simulations for the high-pressure viscosity of complex liquid mixtures.

The flow of thin, confined lubricant films under EHL conditions can deviate from the linear profile generally assumed in continuum models. Recent insights from NEMD simulations suggest that pressure-induced glass formation causes a transition in both the EHL friction and flow behaviour. Hybrid molecular-continuum modelling could give further insights into EHL problems by making larger domain sizes and timescales accessible than is possible through conventional MD simulations. We expect MD simulations and related methods, most of which build on Robbins’ pioneering ideas, to play an increasingly central role in both improving our understanding of lubricant behaviour under EHL conditions and improving models to predict macroscale behaviour.