On the effect of confined fluid molecular structure on nonequilibrium phase behaviour and friction

Abstract

A detailed understanding of the behaviour of confined fluids is critical to a range of industrial applications, for example to control friction in engineering components. In this study, a combination of tribological experiments and confined nonequilibrium molecular dynamics simulations has been used to investigate the effect of base fluid molecular structure on nonequilibrium phase behaviour and friction. An extensive parameter study, including several lubricant and traction fluid molecules subjected to pressures (0.5–2.0 GPa) and strain rates (10⁴–10¹⁰ s⁻¹) typical of the elastohydrodynamic lubrication regime, reveals clear relationships between the friction and flow behaviour. Lubricants, which are flexible, broadly linear molecules, give low friction coefficients that increase with strain rate and pressure in both the experiments and the simulations. Conversely, traction fluids, which are based on inflexible cycloaliphatic groups, give high friction coefficients that only weakly depend on strain rate and pressure. The observed differences in friction behaviour can be rationalised through the stronger shear localisation which is observed for the traction fluids in the simulations. Higher pressures lead to more pronounced shear localisation, whilst increased strain rates lead to a widening of the sheared region. The methods utilised in this study have clarified the physical mechanisms of important confined fluid behaviour and show significant potential in both improving the prediction of elastohydrodynamic friction and developing new molecules to control it.

---

Introduction

To reduce energy consumption and thus CO₂ emissions it is necessary to increase the energy efficiency of engineering components and an important way to achieve this is to design lubricants to decrease friction. Many lubricated engineering components include elements that roll and slide together, for example; rolling bearings, gears, constant velocity joints and cam/follower systems.¹ In these components, much of the friction loss is in the elastohydrodynamic lubrication (EHL) regime, where a very thin (<μm) lubricant film is sheared at very high strain rate (10⁴–10¹⁰ s⁻¹) and pressure (0.5–2.0 GPa) between non-conforming contact surfaces.² Generally, low friction is desirable in order to increase efficiency; however, in some components, such as traction drives, high friction is required and can be obtained using bulky ‘traction fluid’ molecules.³ Despite extensive research, the prediction of friction in the EHL regime remains a considerable challenge, primarily because of the combination of extreme pressure and strain rate which need to be reproduced without an uncontrollable temperature rise. As a result, there is still disagreement regarding the most appropriate rheological model to describe the behaviour of lubricant films in rolling-sliding contacts under EHL conditions.¹

Friction in the EHL regime is significantly affected by the molecular structure of the base fluid.³ To predict EHL friction, one must also consider how the extreme conditions affect the behaviour of the lubricant molecules. Increasing the applied hydrostatic pressure generally increases the viscosity of confined lubricants.² High viscosity, in combination with high strain rate, leads to very high shear stress and this can also significantly affect the lubricant rheology through, for example, shear thinning.¹ There is also growing evidence to suggest that very high shear stresses can lead to nonequilibrium phase transitions which result in deviations from Couette flow.^4–6 However, such phase transitions have proved difficult to capture experimentally and neither the physical driving force behind them, nor their influence on the friction behaviour, is yet well understood.^2,6

Confined nonequilibrium molecular dynamics (NEMD) simulations have proved a valuable tool with which to investigate the friction behaviour of thin films of a range of fluids.^7–12 Significant advantages of confined NEMD simulations include the fact that atomic-scale effects (e.g. phase transitions) can be identified and that the friction behaviour emerges naturally from the simulation, rather than being parameterised into the model.² Also, thermostats can be applied to control the large local temperature rises¹³ which complicate experimental investigations.¹ Confined NEMD simulations of Lennard-Jones fluids have revealed that liquids can transform, at steady state, into ordered, solid-like regions in coexistence with disordered, liquid-like regions.^14,15 Such systems exhibit non-affine flow, whereby the solid-like and liquid-like regions have different average strain rates; a phenomenon known as ‘shear localisation’ or ‘shear banding’.²

Comprehensive parameter studies have revealed many forms of shear localisation in single-component Lennard-Jones fluids, with different locations of the liquid-like and solid-like regions depending on the applied conditions.^2,16,17 When the applied conditions yield shear localisation, the friction coefficient of Lennard-Jones fluids no longer shows the usual increase with strain rate; an observation which has also been made for traction fluids above their limiting shear stress.¹⁸ In addition to Lennard-Jones fluids, shear localisation has also been observed in experiments and NEMD simulations of a range of systems,¹⁹ including granular beads,²⁰ glasses,^21,22 colloids,^23,24 and polymer melts.^25,26 In the field of tribology, shear localisation has only been experimentally observed in very viscous model lubricants, such as polybutadiene and polyphenyl ethers (e.g. 5P4E).^27–33 However, there is currently a growing consensus that the phenomenon may also be important for more realistic lubricants under extreme conditions, for example those experienced in the EHL regime.^4–6

In NEMD simulations of binary mixtures of Lennard-Jones particles, crystallisation was suppressed and only liquid or glassy phases were formed.² In these binary systems, the friction coefficient increased linearly with log(sliding velocity),² as is frequently observed for realistic lubricants, such as polyalphaolefin (PAO), under EHL conditions.^1,18 In this case, the friction–strain rate behaviour can be adequately described using stress-promoted thermal activation theory,³⁴ according to the rheological models of Eyring.³⁵ It has also been suggested that the Carreau³⁶ or Carreau–Yasuda models,^36,37 which were developed for polymer shear thinning, can be used to predict the friction–strain rate behaviour of lubricants in the EHL regime.¹

An underlying assumption in current models for EHL friction is that, in the absence of thermal effects, or until a critical strain rate is reached, the rheological properties of the lubricant do not vary through its thickness. This is implicit in applying the Eyring³⁵ or Carreau–Yasuda^36,37 shear thinning models to predict friction in the EHL regime since they generally assume that the film is subject to Couette shear (i.e. the strain rate is the same as the macroscopic velocity gradient). Thus significant modifications of these models may be required where there are deviations from planar Couette flow, such as shear localisation.¹

In this study, tribometer experiments and NEMD simulations have been used to investigate the effect of base fluid molecular structure, pressure and strain rate on nonequilibrium phase behaviour and friction in the EHL regime. NEMD can simultaneously probe both the friction and flow behaviour within the EHL film, but is limited to high strain rates (>10⁷ s⁻¹) to achieve a satisfactory signal to noise ratio within the available computational resources.^10,38 Experimental friction data serves the dual purpose of providing insights into the behaviour at lower strain rates (<10⁷ s⁻¹) as well as some validation of the NEMD results. This study represents an important step towards utilising NEMD simulations, supported by experiments, to provide information which is critical to; (i) improve macroscale models for EHL friction, and (ii) design new molecules to control EHL friction.¹

Methodology

Tribometer experiments

Test fluids. In this study, two model lubricants were compared with two traction fluids. Squalane (2,6,10,15,19,23-hexamethyltetracosane) is a linear C₂₄ alkane with six methyl branches that has commonly been employed as a model lubricant in both experiments and NEMD simulations.^38–40 Bis[2-ethylhexyl]decanedioate or diethylhexyl sebacate (DEHS) is an example of a synthetic ester oil⁴¹ and its alkyl groups can be varied in order to modify its friction behaviour.³ 2,3-Dimethyl-2-[(3-methylbicyclo[2.2.1]hept-2-yl)methyl]bicyclo[2.2.1]heptane (DM2H) and 2-methy-2,4-dicyclohexylpentane (DCMP) are both traction fluids, which have been specifically designed to give high friction even at low speed in traction drives.³ Squalane (Sigma Aldrich), DEHS (Sigma Aldrich), DM2H (Idemitsu Kosan) and DCMP/Santotrac 2000 (SantoLubes LLC) are all well-characterised, monodisperse fluids which are free of any additives; their molecular structures are shown in Fig. 1. In this manuscript, the term ‘lubricants’ is used to refer to fluids chosen to generally give low friction (squalane, DEHS), while ‘traction fluids’ will be used to refer to molecules chosen to give high friction (DM2H, DCMP). It should, however, be noted that traction fluids also serve a lubricating function.

Fig. 1 Molecular structures of the base fluids investigated; lubricants – squalane (a), DEHS (b) and traction fluids – DM2H (c), DCMP (d).

Methodology. Rolling-sliding EHL friction tests were carried out using a minitraction machine (MTM2, PCS Instruments, UK). In this tribometer, a ball is loaded and rotated against the flat surface of a rotating disc immersed in lubricant at a controlled temperature. The ball and disc are driven by independent motors to enable any combination of sliding and rolling. The drive motors, applied load and temperature are all computer-controlled, and friction is measured from a load cell attached to the ball drive shaft.¹⁸

The main test conditions are listed in Table 1. The choices of load, entrainment speed and slide-roll ratio are explained further below. The MTM2 has a maximum load of 75 N and at this load the mean Hertz pressure in steel/steel contacts cannot exceed 0.8 GPa using a ball-on-flat contact with the smallest ball that can be accommodated. In order to obtain higher contact pressures in a ball-on-flat contact, tests were therefore carried out using tungsten carbide (WC) balls and discs whose high elastic modulus gives a significantly higher contact pressure than steel at the same applied load.¹⁸ Details of the materials, radii, applied loads and corresponding mean Hertz pressures for ball-on-flat contact are shown in Table 2. An EHL ball-on-flat contact produces a pressure distribution that is close to ellipsoidal, with a maximum pressure in the centre of the contact that is 1.5 times the mean pressure.

Table 1 Conditions for tribometer experiments

Entrainment speed, U	2.5 m s^−1
Slide-roll ratio, SRR	0–0.5
Film thickness, h	80–500 nm
Strain rate, [small gamma, Greek, dot above]	10^4–10^7 s^−1
Temperature, T	30–120 °C
Mean pressure, FN	0.55–1.75 GPa

---

Table 2 Contact materials, radii, loads and corresponding mean Hertz pressures for ball-on-flat contact

Materials	Ball radius/mm	Applied load/N	Mean pressure/GPa
Steel/steel	19.05	20.0	0.55
Steel/steel	19.05	50.0	0.75
WC/WC	19.05	15.0	1.01
WC/WC	19.05	30.0	1.27
WC/WC	12.70	22.7	1.51
WC/WC	19.05	75.0	1.73

---

In each test, the fluid temperature was increased from 30 °C to 120 °C in 10 °C increments and at each temperature a ‘traction curve’ was measured. This consisted of measuring the friction while progressively increasing the slide-roll ratio (SRR) from 0 to 0.5 in 24 stages. The SRR is the ratio of the sliding speed, v_s, to the entrainment speed, U. If v₁ and v₂ are the speeds of the ball and disc surfaces respectively, v_s is the difference in speed between the two surfaces, (v₁ − v₂) while U, which was held constant at 2.5 m s⁻¹, is the mean speed of the surfaces with respect to the contact, (v₁ + v₂)/2. In EHL, the lubricant film thickness depends on U, but not on v_s, so a traction curve describes how EHL friction increases with sliding speed at a fixed film thickness.¹⁸ The traction curves can be very simply converted to a curve of shear stress versus strain rate by dividing the measured friction by the applied load to obtain the mean shear stress, and dividing v_s by the film thickness, h, to determine the strain rate. In this study, only the data at 80 °C were used for comparison with NEMD results. Measurements at other temperatures were made to enable correction of the influence of shear heating on friction, as described below.

A feature of EHL contacts is that, because of the combination of high strain rate and high pressure, a large amount of energy is dissipated within the lubricant film and this can lead to large rises in local film temperature.^42,43 This should be accounted for when interpreting EHL friction results and was achieved here by calculating the mean fluid film temperature rise due to shear and then correcting the measured friction data back to their bulk test temperature values.¹ Because sets of shear stress versus strain rate data were collected at several test temperatures for each fluid and load, these could be combined to determine the dependence of the shear stress on temperature at a fixed strain rate. From this information, it was possible to correct the measured shear stresses back to the bulk test temperature.^1,44 Further details of the interpolation thermal correction procedure are shown in the ESI.† It is important to note that the correction was not used to determine values of shear stress at the bulk test temperature a priori but only to correct the original shear stress measurement. Moreover, the correction is generally relatively minor for the small temperature rises that are considered in this study (maximum thermal correction accounted for was 10 °C).^1,44

Fig. 2 shows sets of raw and thermally-corrected mean shear stress versus strain rate curves for squalane and DM2H at 80 °C. The minor thermal corrections applied significantly alter the friction behaviour of the lubricants, maintaining a linear increase in the shear stress with log(strain rate), which would otherwise tail off as a result of shear heating.¹ Conversely, the same thermal correction procedure does not significantly change the shear stress in the traction fluids, since their shear stress is almost independent of temperature.^3,18

Fig. 2 Thermally-corrected (filled symbols) and raw (open symbols) mean shear stress versus log₁₀(strain rate) for; a representative lubricant, squalane (a) and a traction fluid, DM2H (b) at 80 °C. The thermal correction was only applied when calculated fluid film temperature rise was <10 °C.

The composite surface roughness was 12 nm for the steel ball and disc and 16 nm for the WC ball and disc. For all four fluids tested, the EHL film thickness was greater than 80 nm at the entrainment speed of 2.5 m s⁻¹ and a temperature of 80 °C (Table 1). Therefore, the film thickness was always at least five times greater than the composite surface roughness, ensuring a full EHL fluid film with negligible asperity contact.^3,18

NEMD simulations

Model details. Classical, all-atom NEMD simulations were performed using LAMMPS.⁴⁵ The MD equations of motion were integrated using the velocity-Verlet algorithm with a time step of 1 fs. All the systems were constructed using the Materials and Processes Simulations Platform (MAPS) from Scienomics SARL. The systems consisted of a lubricant film confined between atomically-smooth iron oxide (α-Fe₂O₃[100]) slabs, a representative model for steel contact surfaces.⁴⁶ It is important to note that the shear localisation behaviour critical to this study can only be captured in NEMD simulations if the shearing process is boundary-driven, rather than imposed directly on the fluid molecules.²⁵ WC surfaces were not considered in the NEMD simulations since the high pressures, which necessitated the use of WC surfaces in the accompanying experiments, could be readily achieved using iron oxide slabs, for which well-tested, accurate MD parameters are available.⁴⁶ The x, y, z dimensions of the slabs were approximately 5.5 × 5.5 × 1.5 nm.

An example system is shown in Fig. 3. The systems contained a fixed number of molecules (600–1040 depending on fluid) to give a film thickness of approximately 15 nm at the lowest pressure considered (0.5 GPa). The total number of atoms in these systems was approximately 60 000. To make this wide parameter study (Table 3) computationally feasible, the simulated film thickness was much lower than in the accompanying experiments (Table 1). However, in this study, the results from the experiments and NEMD simulations are expected to be directly comparable. This is because, apart from for very small film thicknesses,¹¹ and in the absence of boundary slip,^46–51 the shear stress measured in bulk and confined NEMD simulations are virtually indistinguishable.^7,46,52 Thus the difference in film thickness between the simulations and experiments is not expected to influence the shear stress results. To ensure the validity of this approach, some of the systems were also tested with half and double the normal film thickness. At equal strain rate, these systems gave the same shear stress result within the statistical uncertainty, as shown in the ESI.†

Fig. 3 Representative system setup for the confined NEMD simulations of squalane at 80 °C, after compression (1.0 GPa), before sliding. Iron atoms are shown in pink, oxygen in red, carbon in blue and hydrogen atoms are not shown for clarity. Rendered using VMD.⁵⁴

Table 3 Conditions for NEMD simulations

Sliding speed, vs	0.5–100 m s^−1
Film thickness, h	12–16 nm
Strain rate, [small gamma, Greek, dot above]	10^7–10^10 s^−1
Temperature, T	80 °C
Mean pressure, FN	0.50–2.0 GPa

---

The α-Fe₂O₃ slabs were restrained in the corundum crystal structure using harmonic bonds with a spring constant of 25 kcal mol⁻¹.⁴⁶ Forces between the iron oxide and fluid atoms were governed by Lennard-Jones and electrostatic interactions; the Fe and O parameters used were similar to those developed by Savio et al.⁴⁶ and Berro et al.⁵³ to study the behaviour of alkane films confined between α-Fe₂O₃ surfaces, and are given in the ESI.†

Interactions between the fluid molecules were represented with an updated form of the all-atom ‘optimized potential for liquid simulations’ force-field,⁵⁵ L-OPLS-AA.⁵⁶ In L-OPLS-AA, both bonded and non-bonded parameters have been updated for several atom types to provide a more realistic description of long-chain alkanes⁵⁶ and, more recently, alcohols and esters.⁵⁷ More details of the molecular dynamics force-field used, as well as validation for the molecules used here, are shown in the ESI.† Fast-moving bonds involving hydrogen atoms were constrained with the SHAKE algorithm.⁵⁸ Lennard-Jones interactions were cut off at 12 Å⁵⁶ and electrostatic interactions were evaluated using a slab implementation of the particle–particle–particle–mesh (PPPM) algorithm,⁵⁹ with a relative force accuracy of 10⁻⁵. ‘Unlike’ interactions were evaluated using the geometric mean mixing rules.⁵⁵ Although the computational expense of all-atom force-fields is much greater than for simplified united-atom variants, it has been shown that they are critical for obtaining accurate experimental viscosities of long-chain alkanes.⁶⁰ It is therefore expected that the friction behaviour of these molecules will only be accurately reproduced by all-atom force-fields.⁵⁶ Moreover, all-atom force-fields have also been shown to be critical to obtain liquid–solid phase transitions in strongly confined fluids which agree with experimental observations.⁶¹

Simulation procedure. First, fluid molecules (600–1040) were randomly inserted between the slabs to achieve a density comparable to the fluids under ambient conditions using the ‘Amorphous Builder’ tool in MAPS. The entire system was then energy minimised. The pressure was increased by giving a normal force to the outer layer of atoms in the top slab while keeping the outer layer of atoms in the bottom slab fixed in z (purple box in Fig. 3). The pressures, F_N, considered were between 0.5–2.0 GPa, to provide an envelope for the experimental data, with upper and lower limits slightly above and below those covered in the experiments (Table 1). During the compression phase, a global Langevin thermostat⁶² with a damping coefficient of 0.1 ps was applied in all directions and the systems were allowed to equilibrate at 80 °C in order to match the experiments. Initially, the slab separation varied in a damped harmonic manner, so sliding was not applied until a constant average slab separation was obtained and the average hydrostatic pressure within the liquid film reached its target value.¹⁰ These compression simulations were generally around 2 ns in duration.

After compressive oscillation became negligible, a shear velocity gradient was imposed on the system by means of a constant velocity, v_x = ±v_s/2 applied to the outermost layer of atoms in each slab (purple box in Fig. 3) in the x direction.² Sliding simulations were conducted until the mean shear stress reached a steady state. The necessary simulation time was between 2–50 ns, with lower sliding velocities, v_s, requiring longer simulations. For all of the sliding velocities considered, the shear stress reached a steady state after 5–10 nm of sliding distance (see ESI†). The values of v_s applied were between 0.5–100 m s⁻¹,⁵² yielding strain rates of approximately 10⁷–10¹⁰ s⁻¹ for the film thickness simulated (Table 3). Although lower NEMD strain rates are desirable to overlap with those used in the tribometer experiments (<10⁷ s⁻¹), they are not yet practical for extensive parameter studies using accurate all-atom force-fields.¹⁰ Heat generated during the sliding simulations was dissipated using a Langevin thermostat⁶² acting only on the middle 1 nm of the iron oxide slabs (green box in Fig. 3), applied in the spanwise direction (y) i.e. perpendicular to both sliding and compression.⁶³ This approach is known to be more physically meaningful than applying the thermostat directly to the confined fluid, which has been shown to significantly affect its behaviour under sliding conditions.⁶⁴ In the context of temperature rise in experimental EHL conditions, this approach is similar to removing the effect of bounding surface flash temperature rise, but leaving the effect of temperature rise within the lubricant film due to finite thermal conductivity (see ESI†).¹ However, for the thin films simulated here (h ≈ 15 nm), the latter was also relatively small and only became significant at the highest strain rates considered.¹³

Previous NEMD studies have shown that the boundary thermostatting method applied here was effective in controlling the temperature of confined molecular fluids at similar sliding velocities.^13,65 However, the NEMD shear stress data can only be considered as isothermal at up to sliding velocities of 20 m s⁻¹ ([small gamma, Greek, dot above] ≈ 10⁹ s⁻¹), where there is no observable temperature rise within the film above the thermal noise.⁶⁵ The temperature in the centre of the film¹⁷ increased by approximately 10 °C at 40 m s⁻¹ and 30 °C at 100 m s⁻¹; these temperature rises depended only weakly on the pressure and fluid. For the lubricants (squalane, DEHS), this temperature increase lead to a decrease in the shear stress at the highest sliding velocities considered (open circles in Fig. 4). Conversely, for the traction fluids (DM2H, DCMP), the temperature rise had a negligible effect on the shear stress, which is more weakly dependent on temperature.^3,18

Fig. 4 Mean shear stress versus log₁₀(strain rate) for the fluids at: 80 °C and 0.5–2.0 GPa; squalane (a), DEHS (b), DM2H (c), DCMP (d). Thermally-corrected experimental data shown as filled diamonds. Isothermal NEMD data shown as filled circles, NEMD data with a temperature rise shown as open circles. NEMD data time-averaged for the final 10 nm of sliding. Error bars for NEMD data represents maximum variation between independent trajectories for some state points. Error bars for experimental data represent the uncertainty from the thermal correction. Data from higher applied pressures are shown in darker colours.

Results and discussion

Friction behaviour

The variation in mean shear stress with pressure and strain rate was measured in both tribology experiments and in NEMD simulations, allowing a much wider parameter space to be covered than would be possible through either technique independently. The shear stress in the NEMD simulations was monitored through the mean lateral (friction) force acting on the outer layer of atoms in the top and bottom slabs in response to the fluid.^53,66

Fig. 4 shows the variation in the mean shear stress with log(strain rate) for the lubricants (squalane, DEHS) and traction fluids (DM2H, DCMP) at 80 °C, under several applied mean pressures. Note that experimental data have been thermally-corrected to isothermal conditions as described in the Methodology (Fig. 2) and ESI.† NEMD data with filled symbols are isothermal whilst those at the highest strain rates, with a detectable film temperature rise, are shown as open circles.

Fig. 4 shows good agreement between the thermally corrected experimental data ([small gamma, Greek, dot above] < 10⁷ s⁻¹) and the NEMD simulation data ([small gamma, Greek, dot above] > 10⁷ s⁻¹) for most of the fluids and applied conditions studied. In the experiments, there will be significant variation in the conditions within the macroscopic contact,^67,68 which cannot all be captured in a NEMD simulation limited to the nanoscale.⁶ However, the general agreement found between the experimental and NEMD mean shear stress results for the fluids and conditions studied here suggests that NEMD simulations performed under average conditions for the contact are sufficient to reproduce the important trends. Similar conclusions have been drawn from experimental comparisons between microscopic local viscosity measurements and more conventional global average viscosity measurements for contacts under EHL conditions.^67,68

Fig. 4 shows that the mean shear stress varies significantly with the applied load as well as the strain rate for all of the fluids studied. The friction behaviour is markedly different between the lubricants and the traction fluids. For the lubricants (Fig. 4a and b), the thermally-corrected experimental data and NEMD mean shear stress data are in excellent agreement. Generally, the lubricants show a linear increase in shear stress with log(strain rate), as predicted by the Eyring model.³⁵ This behaviour has been consistently observed in previous experiments and NEMD simulations of lubricants.^1–3,66 At the highest pressures considered (>1.5 GPa), there is a sharp increase in experimentally measured shear stress with log(strain rate) at low strain rates (<10⁵ s⁻¹), due to viscoelastic accommodation of strain by the surface in the entry zone of the contact, in response to the very steep friction increase in this region.^3,69 However, at higher strain rates (>10⁵ s⁻¹), the experimental shear stress shows the same linear increase with log(strain rate) as observed at lower pressure.

The traction fluids (Fig. 4c and d) give much higher shear stress than the lubricants, particularly at low strain rates; however, the slope of the shear stress with log(strain rate) is shallower than for the lubricants. In both the simulations and the experiments, the slope of the shear stress with log(strain rate) (i.e. the Eyring stress) generally increases with increasing pressure for the lubricants, and decreases with increasing pressure for the traction fluids, as has been noted previously.⁶⁹ At low strain rates (<10⁵ s⁻¹), the traction fluids also show a sharp increase in shear stress with log(strain rate) due to viscoelastic accommodation of strain,^3,69 but this occurs at lower pressure (0.55 GPa) than for the lubricants. There is good agreement between the experimental and NEMD simulation mean shear stress data for the traction fluids at low pressure (≤0.55 GPa). At higher pressure, the trends in the NEMD and experiments are similar, but the mean shear stress value, as well as the slope of the shear stress with log(strain rate), are overestimated in the NEMD simulations relative to the experiments. Possible reasons for these discrepancies are discussed later.

The critical strain rate transition, [small gamma, Greek, dot above]_c, above which shear thinning is observed, and below which there is a Newtonian plateau (i.e. where the viscosity, η, is independent of [small gamma, Greek, dot above]), typically occurs at [small gamma, Greek, dot above] ≈ τ⁻¹, where τ is the longest relaxation time at equilibrium (in the absence of shear). This is usually τ_rot, the rotational relaxation time.³⁸ The experimental data probes both the Newtonian plateau and the shear thinning regime for the fluids and conditions considered. The high strain rate NEMD data are mostly in the shear thinning regime, but the transition to the Newtonian plateau can be seen for the lubricants at 0.5 GPa at the lowest strain rates considered (≈10⁷ s⁻¹).

The variation in the friction coefficient with strain rate is also shown to make comparisons between different applied pressures easier.² The friction coefficient, μ, was obtained using the extended Amontons–Coulomb law under the high load approximation: F_L/F_N = F₀/F_N + μ ≃ μ. F_L and F_N are respectively the mean lateral (friction) force and normal force acting on the outer layer of atoms in each slab, and F₀ is the load-independent Derjaguin offset representing adhesive surface forces. The validity of this approximation was confirmed by ensuring that a linear fit of F_L as a function of F_N gave an insignificant value of F₀.¹⁰ The trends evident in Fig. 5 are statistically significant; when reproducing some of the NEMD simulations with different starting configurations, but the same parameters, the maximum difference between the friction coefficient between independent trajectories was approximately 10%. Error bars are omitted from Fig. 5 for clarity.

Fig. 5 EHL friction coefficient versus log₁₀(strain rate) for the fluids at: 80 °C and 0.5–2.0 GPa; squalane (a), DEHS (b), DM2H (c), DCMP (d). Thermally-corrected experimental data shown as filled diamonds. Isothermal NEMD data shown as filled circles, NEMD data with a temperature rise shown as open circles. NEMD data time-averaged for the final 10 nm of sliding. Error bars are omitted for clarity. Data from higher applied pressures are shown in darker colours.

The friction coefficients for the lubricants (Fig. 5a and b) generally increase linearly with log(strain rate) and also increases with the applied pressure. This type of behaviour is commonly observed in experiments and NEMD simulations of lubricants under EHL conditions.^1–3,66 Interestingly, the slope of the friction coefficient with log(strain rate) decreases with increasing pressure, as has been observed in NEMD simulations of binary Lennard-Jones mixtures.² Similar observations have also been made in recent NEMD simulations of n-hexadecane confined between α-Fe₂O₃ slabs at high pressure.⁶⁶ Consequently, the friction coefficient at 0.5 GPa exceeds that at 2.0 GPa at very high strain rates (>2 × 10⁹ s⁻¹).

Although the friction coefficient of the traction fluids (Fig. 5c and d) also increases linearly with log(strain rate), it does so with a much lower slope than for the lubricants. This behaviour is similar to that observed in NEMD simulations of single-component Lennard-Jones fluids subjected to high pressures.^2,16,17 The traction fluids generally gave much higher friction coefficients than lubricants, particularly at low strain rate. The friction coefficient at 0.5 GPa exceeds that at 2.0 GPa at much lower strain rates (>10⁸ s⁻¹) for the traction fluids than for the lubricants (>2 × 10⁹ s⁻¹). The high friction from the traction fluid molecules can be attributed to interlocking of the bulky cyclohexane and bicyclo[2.2.1]heptyl groups DCMP and DM2H respectively.⁷⁰ This, coupled with their internal molecular stiffness,⁷¹ increase energy barriers for neighbouring molecules to slide over one another,³⁵ leading to high friction.³

Flow behaviour

The physical state of fluids in the EHL regime is still largely unresolved, owing to experimental difficulties in probing the contact zone under the extreme conditions.² The nonequilibrium phase behaviour of the fluids can be directly probed using velocity and atomic mass density profiles from the NEMD simulations.^2,16,17 Although not part of this current study, experimental velocity profiles for these fluids are also currently under investigation using recently developed techniques.^30–33 The NEMD velocity and atomic mass density profiles are spatially and time averaged to improve their signal to noise ratio. The profiles are both shifted such that a z-coordinate of zero is at the centre of the fluid film. The velocity profiles generally reach a steady state after 5–10 nm of sliding, depending on the fluid and conditions.

Fig. 6 shows how the x-velocity profile in the z-dimension, v_x(z), and the atomic mass density profile in the z-dimension, changes with the applied pressure for a representative lubricant, squalane. At 0.5 GPa (Fig. 6a), a Couette flow profile develops, with a linear velocity gradient between the two slabs. The outer molecular layer of fluid moves at the same velocity as the wetting α-Fe₂O₃ slabs, indicating that no boundary slip occurs (see ESI†). The oscillatory atomic mass density profile indicates strong layering of the squalane molecules close to the surface and weaker layering in the centre of the film. These profiles are similar to those from NEMD simulations of thin n-hexadecane films confined by α-Fe₂O₃ slabs at high pressure.^9,66 At 1.0 GPa (Fig. 6b), there is stronger layering of the fluid extending further from the slabs, suggesting the formation of an ordered, solid-like region close to the slabs, with a liquid-like region in the centre of the film. The flow profiles in Fig. 6 confirm the transition to a ‘central localisation’ (CL) phase, with the solid-like regions of the fluid moving at the same velocity as the slabs and a steeper velocity gradient within the central liquid-like region than for the Couette flow case. CL is one of the nonequilibrium phases identified in previous confined NEMD simulations of Lennard-Jones fluids,^2,16,17 and has also been observed experimentally for 5P4E.^27–29,32 At 2.0 GPa (Fig. 6c), the molecular layering is even more pronounced, and extends further into the centre of the film. Thus, the system shows stronger CL, with thicker solid-like regions moving with the slab and a steeper velocity gradient in the central liquid-like region. It is important to note that the strong layering extends much further into fluid than is directly influenced by the solid slab (Lennard-Jones interactions cut off at 12 Å). It is not possible from these simulations to definitively show the physical state of the solid-like region, whether glassy or fully crystalline; this will be investigated in future studies. The change in the velocity profiles with pressure, moving from Couette-like to increasingly strong CL, is similar to that observed experimentally for viscous 5P4E at lower temperature, strain rate, and pressure.³²

Fig. 6 Atomic mass density (blue) and velocity (orange) profiles for squalane at: 80 °C, 20 m s⁻¹ ([small gamma, Greek, dot above] ≈ 10⁹ s⁻¹), and; 0.50 GPa (a), 1.00 GPa (b), 2.00 GPa (c). Time-averaged for the final 5 nm of sliding.

Fig. 7 shows how the x-velocity profile in the z-dimension, v_x(z), and the atomic mass density profile in the z-dimension, change with the applied pressure for a representative traction fluid, DM2H. The trends for the traction fluids are similar to those for the lubricants, but they show greater divergence from Couette flow under the same conditions. Even at 0.5 GPa (Fig. 7a), DM2H shows weak CL, with layering of the molecules close to the surface resulting in a sinusoidal velocity profile. At 1.0 GPa (Fig. 7b), the traction fluids show much stronger CL, which is similar to the lubricants at 2.0 GPa (Fig. 6c). At 2.0 GPa (Fig. 7c), the traction fluids show ‘plug slip’ (PS), another phase identified in confined NEMD simulations of Lennard-Jones fluids.^2,16,17 Here, the velocity profile shows that shear localises close to the slabs and the density profile indicates strong layering throughout the film.¹⁶ During PS, the first molecular layer still moves at the same speed as the slab and ‘slip planes’ form within the second and third layers.¹⁶ Note that PS is distinct from the boundary slip behaviour, which is commonly observed in NEMD simulations of very thin films (see ESI†).⁹ PS has been observed experimentally for viscous polybutadiene^30,31 and has also been inferred from film thickness measurements in realistic mineral oils.^4,5

Fig. 7 Atomic mass density (blue) and velocity (orange) profiles for DM2H at: 80 °C, 20 m s⁻¹ ([small gamma, Greek, dot above] ≈ 10⁹ s⁻¹), and; 0.50 GPa (a), 1.00 GPa (b), 2.00 GPa (c). Time-averaged for the final 5 nm of sliding.

Fig. 8 shows the change in the velocity profile with sliding velocity for a representative fluid, DM2H. The width of the central liquid-like region increases with sliding velocity from 5–100 ms⁻¹, as has been observed previously for Lennard-Jones fluids undergoing CL.^2,16,17 This is intuitive given that, at higher sliding velocities (strain rates), there is increased shear heating in the centre of the film.¹⁷ At the lowest sliding velocities considered (0.5–2.0 ms⁻¹), the lower signal to noise ratio in the velocity profiles makes it difficult to definitively state the flow behaviour.

Fig. 8 Atomic mass density (blue) and velocity (orange) profiles for DM2H at: 80 °C, 1.0 GPa, and; 5 m s⁻¹ (a), 40 m s⁻¹ (b), 100 m s⁻¹ (c). Time-averaged for final 5 nm of sliding.

In the confined molecular fluids simulated here, the nonequilibrium phase behaviour is broadly similar to that observed in NEMD simulations of Lennard-Jones fluids.^2,16,17 As the pressure is increased, molecules are forced closer together, and eventually localisation of the shear becomes energetically favourable to global shearing of the film. Under the relatively high strain rates accessible through NEMD, this generally manifests as CL, but the transition to PS was also observed for the traction fluids at the highest pressures considered.¹⁷ The results of this study suggest that, precluding excessive shear heating, shear localisation will arise in any confined molecular fluid providing that the pressure is sufficiently high.²

The shear bands form parallel to the confining walls in these and, to our knowledge, all previous NEMD simulations which have investigated this phenomenon.⁶ Conversely, bands have been observed at a 20° angle to the confining surfaces in experiments using 100 μm thick 5P4E films,²⁸ which was later supported by theoretical predictions.²⁹ The bands in the NEMD simulations probably form parallel to the surfaces rather than at an angle because relatively thin films are used with periodic boundary conditions in the x and y directions in order to make the simulations computationally feasible.^6,17

Relationship between flow and friction

Given the fact that shear localisation occurs for a wide range of condensed phase materials (see Introduction), it might be reasonably expected that a common principle applies across the range of different types.¹⁹ The critical question to consider is: why would a fluid film shear inhomogeneously rather than with a linear velocity gradient? The principle of least work or minimisation of heat production may go some way to explain this effect, which is encapsulated in basic form in the following equation, for a fluid undergoing Couette flow:

q̇ = η[small gamma, Greek, dot above]²hS² (1)

where q̇ is the rate of heat production, η is the shear viscosity, [small gamma, Greek, dot above] is the strain rate, and S is a cross-sectional dimension. For Couette flow, [small gamma, Greek, dot above] = v_s/h, where v_s is the difference in velocity between the slabs and h is the measured thickness of the lubricant film. If there is shear localisation over a thickness, Δh, then, assuming the solid-like region is not sheared, the ‘actual’ strain rate is v_s/Δh. The localisation thickness, Δh, can range from near zero to h, the measured film thickness. The product of the actual strain rate and the localisation thickness, [small gamma, Greek, dot above]Δh, must be constant to a first approximation (the exact value of Δh can be difficult to specify). Hence, from eqn (1):

q̇ ∝ η(γ)[small gamma, Greek, dot above] (2)

There is an energetic advantage, therefore, in having localisation of the shear as, if the actual strain rate is larger than the Couette case and the fluid is shear thinning,³⁸ then the rate of heat production (or equivalently work performed) will be reduced in proportion to the decrease in viscosity due to enhanced shear thinning.¹⁶ Moreover, the solid-like region will have an increased thermal conductivity,³ leading to more efficient thermal dissipation when the shear is localised. A combination of these factors mean that, at high pressure shear localisation can reduce the magnitude of the friction coefficient relative to Couette flow.

The slope of the friction coefficient with log(strain rate) is also reduced when fluids undergo shear localisation. This is intuitive since the region in which the shear is localised can be viewed as a strongly shear thinning liquid or, in extreme cases, as a slip plane between two solid-like layers.¹⁶ In the former case, because the thickness of the liquid-like region grows with increasing strain rate (see Fig. 8), the actual strain rate, and thus the friction coefficient, becomes less sensitive to the applied strain rate. In the latter case, invariance of the friction coefficient to the applied strain rate is expected, as it is for sliding between elastic–plastic solids.⁶⁹

Recall that increased pressure generally leads to stronger shear localisation (Fig. 6 and 7), an increase in the magnitude of the friction coefficient, but a reduction in the slope of the friction coefficient with log(strain rate) (Fig. 5). Thus, although the reduction in the slope from stronger shear localisation can clearly be observed in Fig. 5, the decrease in the magnitude of the friction coefficient must be masked by other effects. Specifically, increasing the pressure also increases the fluid viscosity, leading to higher friction.¹ Although the magnitude of the friction coefficient increases in systems at higher pressure, which show more shear localisation, it is lower than would be expected were Couette flow to be enforced under the same conditions.

When the limiting shear stress is reached, the friction coefficient becomes insensitive to both strain rate and pressure.⁶ The results of this study (Fig. 5) suggest that this will only occur for bulky traction fluid molecules at high pressure, which show significant deviation from Couette flow, i.e. strong shear localisation (Fig. 7). Thus, these results support strong shear localisation as the physical origin of the limiting shear stress, one of several proposed explanations in the literature.^6,72 One might also tentatively suggest that the limiting shear stress is reached upon the transition from CL to PS,⁶ although more research is needed to confirm this.

The experimental results suggest that the limiting shear stress was approximately 0.1 GPa for both DM2H and DCMP at 1.0 GPa and 80 °C, which is similar to that reported from earlier experiments for the latter.⁶⁹ The NEMD simulations (Fig. 5c and d) suggest a higher limiting shear stress (approximately 0.25 GPa) and at 1.0 GPa applied pressure, the friction coefficient still increases with log(strain rate). This discrepancy, along with the overestimated friction coefficient for the traction fluids at high pressure (Fig. 5c and d) may have at least two explanations. Firstly, the phase transitions critical to the friction behaviour, may occur under different conditions in the simulations (dependent on the MD force-field)⁵⁶ than in the experiments. Specifically, for the traction fluids, the friction results suggest that shear localisation, which reduces both the magnitude of the friction coefficient and its slope with log(strain rate), occurs under milder conditions, and to a greater degree, in the experiments compared to the NEMD simulations. Secondly, the thinner NEMD film thickness (Table 3) relative to the experiments (Table 1), could lead to differences in the relative sizes of the solid-like and liquid-like regions under equivalent conditions as well as variations in thermal dissipation. It should be noted, however, that there is generally good agreement between the results of the tribology experiments and NEMD simulations, particularly for the lubricants.

The differences in the friction behaviour between the two lubricants are clear in both the experiments and the NEMD simulations. Squalane (Fig. 4a) gives a consistently higher friction than DEHS (Fig. 4b) under all the conditions studied. Molecular factors that promote low EHL friction include linear chains and flexible groups, such as ester linkages, that reduce energy barriers for neighbouring molecules to slide over one another.³ Both molecules include linear backbones of the same length (24 atoms); however, squalane includes several methyl branches, while DEHS has two ethyl branches as well as two flexible ester linkages (Fig. 1). This additional flexibility and molecular ‘smoothness’ leads to a significantly lower friction for DEHS than squalane. In future studies, more molecular structures could be compared and similar methodologies to those used here could be used to design new molecules with specific levels of EHL friction.

The differences between the two traction fluids are far more subtle. At low pressure (≤0.55 GPa), the shear stress of DCMP (Fig. 4d) has a slightly steeper slope with log(strain rate) than for DM2H (Fig. 4c). This difference can be explained through the NEMD simulations, where DCMP has a more Couette-like profile at 0.50 GPa (see ESI†) while DM2H shows some central localisation (Fig. 7a). At higher pressures (>0.55 GPa), DM2H and DCMP have virtually identical friction and flow behaviour in both the experiments and NEMD simulations. Both traction fluids have similar molecular structures (Fig. 1), with cyclohexane groups in their chair confirmation meaning that they have steric bulk in three dimensions.

An increased knowledge of the friction and flow behaviour of confined fluids from tribology experiments and NEMD simulations, as used in this study, could help to accelerate the design of both new lubricants which give lower friction, and new traction fluids that give higher friction. For example, lubricants which are more susceptible to shear localisation under EHL conditions could lead to lower friction. Shear localisation could be perhaps be promoted by using strongly layered lubricants, such as ionic liquids.^12,73,74 Conversely, including a mixture of different traction fluid molecules to suppress the formation of solid-like regions and thus shear localisation should increase EHL friction, particularly at high strain rates.

Conclusions

A combination of tribometer experiments and confined NEMD simulations has been used to investigate the effect of molecular structure on friction and nonequilibrium phase behaviour. An extensive parameter study, including several lubricant and traction fluid molecules subjected to pressures (0.5–2.0 GPa) and strain rates (10⁴–10¹⁰ s⁻¹) typical of the EHL regime, suggests clear relationships between the friction and flow behaviour.

In both the tribometer experiments and NEMD simulations, lubricant molecules generally give a low friction coefficient which increases with both strain rate and pressure. Conversely, traction fluid molecules give a high friction coefficient which has a much weaker dependence on strain rate and pressure. This behaviour can be rationalised through analysis of flow profiles from the NEMD simulations which suggest significant shear localisation, particularly in the traction fluids.

Nonequilibrium phase changes, which were first identified in NEMD simulations of simple Lennard-Jones fluids, have also been observed here, for the first time, in well-characterised molecular fluids. All of the fluids investigated display some shear localisation at high pressure. Both lubricants and traction fluids show central localisation, where the regions of the fluid close to the slab become solid-like and only the central, liquid-like region is sheared. At very high pressure, the traction fluids also show plug slip, where the centre of the film becomes solid-like and the shear occurs close to, but not at, the confining walls. The traction fluid molecules, which have bulky cycloaliphatic groups, show a higher propensity for shear localisation, and only exhibit Couette flow at low pressure or very high strain rate. Higher pressure induced more pronounced shear localisation, whilst increasing the strain rate resulted in a widening of the sheared, liquid-like region. Shear localisation generally leads to a decrease in both the magnitude of the friction coefficient and the slope of its increase with log(strain rate) compared to Couette flow.

For most of the fluids and conditions considered, the variation in the shear stress with log(strain rate) in the NEMD simulations is in good agreement with thermally corrected values from tribometer experiments conducted at lower strain rates. This study has important implications with respect to utilising tribometer experiments and NEMD simulations to predict and control EHL friction. Specifically, they suggest that a range of molecules show significant deviations from planar Couette flow under EHL conditions, meaning that models currently used to predict friction in this regime may require modification to account for this. Moreover, the techniques applied here show significant potential to accelerate the design of new lubricant and traction fluid molecular structures.

Acknowledgements

J. P. E acknowledges the financial support of the Engineering and Physical Sciences Research Council (EPSRC) via a Case Conversion studentship and Shell Global Solutions via the University Technology Centre for Fuels and Lubricants at Imperial College London. D. D. thanks the EPSRC for an Established Career Fellowship EP/N025954/1. The authors acknowledge the use of Imperial College HPC service, URL: http://www.imperial.ac.uk/admin-services/ict/self-service/research-support/hpc/. We would like to thank Janet Wong and Stephen Jeffreys (Imperial College) and Neal Morgan (Shell Global Solutions) for useful discussions. All data reported in the manuscript can be obtained by emailing the corresponding author or http://tribology@imperial.ac.uk.

References

1. H. Spikes and J. Zhang, Tribol. Lett., 2014, 56, 1–25.
2. C. Gattinoni, D. M. Heyes, C. D. Lorenz and D. Dini, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 52406.
3. J. Zhang, A. Tan and H. Spikes, Tribol. Lett., 2017, 65, 13.
4. P. Sperka, I. Krupka and M. Hartl, Tribol. Lett., 2014, 54, 151–160.
5. P. Sperka, I. Krupka and M. Hartl, Friction, 2016, 4, 380–390.
6. L. Martinie and P. Vergne, Tribol. Lett., 2016, 63, 21.
7. S. A. Gupta, H. D. Cochran and P. T. Cummings, J. Chem. Phys., 1997, 107, 10335–10343.
8. A. Jabbarzadeh, P. Harrowell and R. I. Tanner, J. Phys. Chem. B, 2007, 111, 11354–11365.
9. D. T. Ta, A. K. Tieu, H. T. Zhu and B. Kosasih, J. Chem. Phys., 2015, 143, 164702.
10. J. P. Ewen, C. Gattinoni, N. Morgan, H. Spikes and D. Dini, Langmuir, 2016, 32, 4450–4463.
11. S. T. Cui, P. T. Cummings and H. D. Cochran, J. Chem. Phys., 2001, 114, 7189–7195.
12. N. Voeltzel, A. Giuliani, N. Fillot, P. Vergne and L. Joly, Phys. Chem. Chem. Phys., 2015, 17, 23226–23235.
13. D. Toton, C. D. Lorenz, N. Rompotis, N. Martsinovich and L. Kantorovich, J. Phys.: Condens. Matter, 2010, 22, 74205.
14. S. Butler and P. Harrowell, J. Chem. Phys., 2003, 118, 4115–4126.
15. S. Butler and P. Harrowell, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 51503.
16. S. Maćkowiak, D. M. Heyes, D. Dini and A. C. Brańka, J. Chem. Phys., 2016, 145, 164704.
17. D. M. Heyes, E. R. Smith, D. Dini, H. A. Spikes and T. A. Zaki, J. Chem. Phys., 2012, 136, 134705.
18. J. Zhang and H. Spikes, Tribol. Lett., 2016, 63, 24.
19. G. Ovarlez, S. Rodts, X. Chateau and P. Coussot, Rheol. Acta, 2009, 48, 831–844.
20. P. Schall and M. van Hecke, Annu. Rev. Fluid Mech., 2010, 42, 67–88.
21. Y. Shi and M. L. Falk, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 214201.
22. F. Varnik, L. Bocquet, J.-L. Barrat and L. Berthier, Phys. Rev. Lett., 2003, 90, 95702.
23. K. Yeo and M. R. Maxey, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 51502.
24. W. M. Holmes, P. T. Callaghan, D. Vlassopoulos and J. Roovers, J. Rheol., 2004, 48, 1085.
25. J. Cao and A. E. Likhtman, Phys. Rev. Lett., 2012, 108, 28302.
26. M. Mohagheghi and B. Khomami, ACS Macro Lett., 2015, 4, 684–688.
27. S. Bair and W. O. Winer, J. Tribol., 1992, 114, 1–9.
28. S. Bair, F. Qureshi and W. O. Winer, J. Tribol., 1993, 115, 507–514.
29. S. Bair and C. McCabe, Tribol. Int., 2004, 37, 783–789.
30. A. Ponjavic, M. Chennaoui and J. S. S. Wong, Tribol. Lett., 2013, 50, 261–277.
31. A. Ponjavic, L. di Mare and J. S. S. Wong, J. Polym. Sci., Part B: Polym. Phys., 2014, 52, 708–715.
32. B. Galmiche, A. Ponjavic and J. S. S. Wong, J. Phys.: Condens. Matter, 2016, 28, 134005.
33. A. Ponjavic and J. S. S. Wong, RSC Adv., 2014, 4, 20821–20829.
34. H. Spikes and W. Tysoe, Tribol. Lett., 2015, 59, 14.
35. H. Eyring, J. Chem. Phys., 1936, 4, 283–291.
36. P. J. Carreau, J. Rheol., 1972, 16, 99.
37. K. Yasuda, R. C. Armstrong and R. E. Cohen, Rheol. Acta, 1981, 20, 163–178.
38. S. Bair, C. McCabe and P. T. Cummings, Phys. Rev. Lett., 2002, 88, 58302.
39. S. Bair, C. McCabe and P. T. Cummings, Tribol. Lett., 2002, 13, 251–254.
40. S. Bair, Proc. Inst. Mech. Eng., Part J, 2002, 216, 139–149.
41. Y. He, T. J. Zolper, P. Liu, Y. Zhao, X. He, X. Shen, H. Sun, Q. Duan and Q. Wang, Friction, 2015, 3, 243–255.
42. J. F. Archard, Wear, 1959, 2, 438–455.
43. K. L. Johnson and J. A. Greenwood, Wear, 1980, 61, 353–374.
44. H. Spikes and J. Zhang, Tribol. Lett., 2015, 58, 17.
45. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
46. D. Savio, N. Fillot, P. Vergne and M. Zaccheddu, Tribol. Lett., 2012, 46, 11–22.
47. N. Fillot, H. Berro and P. Vergne, Tribol. Lett., 2011, 43, 257–266.
48. H. Washizu, S. Hyodo, T. Ohmori, N. Nishino and A. Suzuki, Tribol. Online, 2014, 9, 45–50.
49. A. Jabbarzadeh, J. D. Atkinson and R. I. Tanner, J. Chem. Phys., 1999, 110, 2612–2620.
50. A. Martini, H. Y. Hsu, N. A. Patankar and S. Lichter, Phys. Rev. Lett., 2008, 100, 4.
51. J. P. Gao, W. D. Luedtke and U. Landman, Tribol. Lett., 2000, 9, 3–13.
52. M. J. Stevens, M. Mondello, G. S. Grest, S. T. Cui, H. D. Cochran and P. T. Cummings, J. Chem. Phys., 1997, 106, 7303–7314.
53. H. Berro, N. Fillot and P. Vergne, Tribol. Int., 2010, 43, 1811–1822.
54. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics Modell., 1996, 14, 33–38.
55. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
56. S. W. I. Siu, K. Pluhackova and R. A. Bockmann, J. Chem. Theory Comput., 2012, 8, 1459–1470.
57. K. Pluhackoya, H. Morhenn, L. Lautner, W. Lohstroh, K. S. Nemkovski, T. Unruh and R. A. Bockmann, J. Phys. Chem. B, 2015, 119, 15287–15299.
58. J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
59. I. C. Yeh and M. L. Berkowitz, J. Chem. Phys., 1999, 111, 3155–3162.
60. J. P. Ewen, C. Gattinoni, F. M. Thakkar, N. Morgan, H. Spikes and D. Dini, Materials, 2016, 9, 651.
61. H. Docherty and P. T. Cummings, Soft Matter, 2010, 6, 1640–1643.
62. T. Schneider and E. Stoll, Phys. Rev. B: Condens. Matter Mater. Phys., 1978, 17, 1302–1322.
63. S. Y. Liem, D. Brown and J. H. R. Clarke, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 3706–3713.
64. S. Bernardi, B. D. Todd and D. J. Searles, J. Chem. Phys., 2010, 132, 244706.
65. H. Berro, N. Fillot, P. Vergne, T. Tokumasu, T. Ohara and G. Kikugawa, J. Chem. Phys., 2011, 135, 10.
66. T. D. Ta, A. K. Tieu, H. Zhu, B. Kosasih, Q. Zhu and H. T. Phan, Tribol. Int., 2017, 113, 22–35.
67. A. Ponjavic, J. Dench, N. Morgan and J. S. S. Wong, RSC Adv., 2015, 5, 99585–99593.
68. J. Dench, N. Morgan and J. S. S. Wong, Tribol. Lett., 2017, 65, 25.
69. C. R. Evans and K. L. Johnson, Proc. Inst. Mech. Eng., Part C, 1986, 200, 303–312.
70. H. Yamano, K. Shiota, R. Miura, M. Katagiri, M. Kubo, A. Stirling, E. Broclawik, A. Miyamoto and T. Tsubouchi, Thin Solid Films, 1996, 281, 598–601.
71. H. Washizu and T. Ohmori, Lubr. Sci., 2010, 22, 323–340.
72. H. Kobayashi and Y. Fujita, J. Appl. Phys., 2014, 115, 223509.
73. A. C. F. Mendonca, A. A. H. Padua and P. Malfreyt, J. Chem. Theory Comput., 2013, 9, 1600–1610.
74. F. F. Canova, H. Matsubara, M. Mizukami, K. Kurihara and A. L. Shluger, Phys. Chem. Chem. Phys., 2014, 16, 8247–8256.

---

Footnote

† Electronic supplementary information (ESI) available: Discussion of the thermal correction procedure, the molecular dynamics force-field, boundary slip and film thickness as well as additional flow profiles. See DOI: 10.1039/c7cp01895a

---