Molecular dynamics simulations of sliding silica layers passivated with hydroxyl groups and n-pentanol chains that act as lubricants at the contact surface under an external load (normal pressure) and shear forces (shear stress) were performed. The systems could evolve from unequilibrated conformations to stationary states at constant speed. As the system evolved from the initial unequilibrated conformation, the hydroxyl groups of the n-pentanol chains moved to positions near the silica surfaces and formed hydrogen bonds with the surface hydroxyl groups. The number of n-pentanols placed in the contact zone matched the number of hydroxyl groups on the two contacting surfaces; therefore, the silica surfaces were fully saturated with a monolayer of adsorbed n-pentanol molecules. Due to the lower frictional forces between the terminal methyl groups of the n-pentanol chains (compared to the hydrogen-bond forces between the surface and n-pentanol hydroxyl groups), the contact surface where friction and sliding occurs is formed by the terminal methyl groups. A snapshot of the two silica surfaces exposed to external forces in the normal and lateral directions at a constant speed, several ns after the initial equilibration period (i.e., stationary state), is shown in Fig. 3.
An external force of 200 kcal/(mol Å) was applied in the normal direction, which is equivalent to a force of 13,984 pN or PN = 552.56 MPa (A = 53.9 × 46.7 Å2), and as explained in the methodology section, this magnitude of pressure was applied in the normal direction to both surfaces in opposite directions to simulate a compressive load in the system. We applied two PN of 552.56 MPa to the system; however, compared with experiments, these pressures do not add or cancel because they are applied on different parts of the system. The first represents the load, and the second represents the response to the applied load of the surroundings of the system. The applied shear stress on each layer was PN/4 = 138.14 MPa in opposite directions.
To compare with experiments, the magnitude of the shear stress should be the sum of the shear stresses applied to the two layers in the simulations; therefore, τ = PN/2 = 276.28 MPa. Once the system reached a steady state where all the n-pentanols had been adsorbed to the hydroxyl groups of the silica surface, the layers periodically underwent stick and slip events or continuous sliding (depending on the magnitude of τ). The n-pentanol chains were not completely perpendicular to the surfaces; they had an angle of inclination because of the load (Fig. 3, left). The silica layers were free to move in the lateral direction perpendicular to the sliding direction, but the layers did not move in this direction (Fig. 3, right).
As the system evolved, the silica layers began to slide in opposite directions, as shown in the time profile of the position of the center of mass of each layer (Fig. 4). Since the layers were under periodic conditions in the lateral directions, the layers can slide infinitely in the sliding direction. In the first moments, the layers moved very quickly because the layers under load were not yet in contact. The layers contacted after ~ 0.03 ns and, for a period that lasts ~ 0.26 ns, the layers moved at a low speed. Then, probably after a critical level of adsorption of the lubricant molecules was reached, the layers accelerated their displacement until they reached a larger constant speed in a stationary state at times greater than 0.97 ns. The initial point of displacement is not zero; it corresponds to the position of the initial center of mass of the layers, which is half the size of the simulation cell in the direction of displacement.
Three regimes can be observed in the temporal profile of the absolute displacement of the layers (Fig. 5), obtained from the difference between the displacement of the layers. These regimes can be distinguished by the slopes of the profile. These slopes (speeds) can be associated with the level of adsorption and ordering of the n-pentanol molecules. In the first regime, at the beginning of the simulation, none of the n-pentanols are adsorbed on the silica surfaces, and at long times (> 0.97 ns), all the n-pentanols are likely adsorbed. The three temporal regimes are denoted in Fig. 5 as periods I, II, and III. Periods I and II are characterized by having logarithmic and exponential growth, respectively, while period III is characterized by a constant speed, resembling a stationary state. The time limits between periods II and III are clearly delimited by the slope change of the absolute displacement profile, as shown in the inset plot of Fig. 5. However, the time limit between periods I and II is not clearly defined, only the change in the behavior of the profile indicates the existence of two periods.
The distribution of the oxygen atoms of the n-pentanol chains as a function of the position between the two silica layers is shown in Fig. 6a for the first and last 0.05 ns of period I. The corresponding results for period II are shown in Fig. 6b. In Fig. 6c, the distribution of the first 0.05 ns of period III is shown. The probability plots show normal distributions of the position of the hydroxyl groups of the silica layers, which became narrow and taller as the simulation progressed and reached the stationary state in period III. An analysis of the integration of the probability profiles showed how many lubricant molecules were not adsorbed at the silica surfaces. Most n-pentanol molecules were adsorbed on the silica layers in a very short period. From the first 0.05 ns of the simulation in period I, a large fraction (81%) of the lubricant molecules were adsorbed on the silica surfaces through hydrogen bonds between the hydroxyls of the n-pentanol chains and the surface hydroxyls of the silica layers. Through period I, characterized by a logarithmic growth of the absolute displacement of the layers, more n-pentanol chains were adsorbed, reaching 86% of the total chains in the last 0.05 ns of period I. During the first 0.05 ns of period II, characterized by a constant acceleration of the absolute displacement, the fraction of adsorbed molecules reached 88%, while at the end of period II (last 0.05 ns), almost 100% of the n-pentanol chains were adsorbed on the silica layers. In period III, characterized by the constant velocity of the absolute displacement, the fraction of adsorbed chains did not vary.
In period I, the hydroxyls of the n-pentanols that are not adsorbed at the silica layers likely produced some hydrogen-bond forces between those lubricant molecules away from the silica surfaces, increasing the frictional forces between the layers. These layers moved at a very low average velocity of 35.21 nm/ns, which is six times lower than that obtained at times greater than 0.97 ns (~ 205 nm/ns), as shown in Fig. 5 (period III). In period II, a very large fraction of n-pentanols was adsorbed on the passivated silica surfaces (Fig. 6, bottom) and the absolute displacement of the layers accelerated until it reached a stationary state of sliding in period III. The maximum adsorption obtained through period II is probably the result of the accelerated displacement of the layers, which exponentially increased the velocity of the displacements, acting as a homogenizing force.
At the silica surfaces, the network of hydrogen bonds between the hydroxyl groups of the silica layers and the hydroxyl groups of the lubricant molecules developed quickly. In Fig. 7a, a snapshot of the network of hydrogen bonds at 0.12 ns shows that most of the hydroxyl groups of the lubricant molecules are already present between two hydroxyl groups, even at this early stage of the simulation. The hydroxyl groups from the lubricant molecules can be located between contiguous hydroxyl groups from the silica layer in the sliding direction or between laterally contiguous groups, forming a 30° angle with the sliding directions. Further, a few of the hydroxyl groups from the n-pentanols formed small clusters in period I. When the simulation entered the constant acceleration period II (Fig. 7b), all the small clusters disappeared, and all the hydrogen bonds were homogeneously distributed at the surfaces. In the stationary state with constant velocity (period III), the homogeneous distribution of contiguous and laterally contiguous hydrogen bonds was maintained (Fig. 7c). The hydroxyl groups of the lubricant molecules did not stay between the same hydroxyl groups for the entire simulation: they slowly moved. When instability events occur during the sliding, large and local shear forces induce some hydroxyl groups from the lubricant molecules to move from contiguous to laterally contiguous hydrogen bonds and vice versa, which is an energetically less expensive route than a hydroxyl group directly moving in the sliding direction of the surface layer.
A more detailed analysis of the acceleration period (period II), which occurs between the logarithmic growth of period I and the long-term stationary behavior of period III, is also shown in Fig. 5, which plots a quadratic expression between the absolute displacement and the time, t:
$$\text{s}\left(\text{t}\right)\text{ = 16.81 nm + }\left(\text{36.95 nm/ns}\right)\text{ }\left(\text{t - 0.29 ns}\right)\text{ + }\left(\text{163.75 nm/}{\text{ns}}^{\text{2}}\right){\text{ }\left(\text{t - 0.29 ns}\right)}^{\text{2}}$$
1
.
We compared the previous expression with the expression corresponding to a movement subject to constant acceleration:
$$\text{s}\left(\text{t}\right)\text{ = }{\text{s}}_{\text{0}}\text{ + }{\text{v}}_{\text{0 }}\left(\text{t -}{\text{ t}}_{\text{0}}\right)\text{ + }\frac{\text{1}}{\text{2}}\text{ }{\text{a}}_{\text{c}}\text{ }{\left(\text{t - }{\text{t}}_{\text{0}}\right)}^{\text{2}}$$
2
,
where s0, v0, and t0 are the initial position, initial velocity, and time when period II begins, respectively, and ac is the constant acceleration. To obtain the best fit to Eq. (2), several subsets of data, starting at different initial times and all ending at the same time of 0.97 ns (the time frontier between periods II and III, located as shown in the inset plot of Fig. 5) were analyzed. The best match between the end speed in period I and v0 in Eq. 2 was chosen. Comparing equations (1) and (2), we obtained the constant acceleration for this period (~ 327.50 nm/ns2), which is ~ 33.4 times the acceleration due to gravity. The rapid acceleration in the transition period (II) to the final stationary state at constant velocity (period III) is the result of the high shear stress used (half the normal load).
Various additional shear stresses were studied: 82.88, 110.51, 138.14, 165.77, 193.40, 221.02, and 248.65 MPa, while maintaining the normal load of 552.56 MPa and a constant temperature of 300 K. At τ = 248.65 MPa (0.45 PN), a behavior similar to that observed at τ = 276.28 MPa was observed (Fig. 8). Three characteristic periods were present: the initial transition state characterized by the presence of some hydroxyl groups from the n-pentanol chains in the contact surface (period I) was two times longer than the period observed at τ = 248.65 MPa. Once the critical fraction of adsorbed molecules of the lubricant was achieved, the constant-acceleration period developed (period II). Under these conditions, it took longer to reach the final stationary state of period III: for this system, it took about 2 times longer (~ 2 ns) than that estimated using τ = 248.65 MPa. The constant speed that was reached at long times (period III) was 147.05 nm/ns. A 10% reduction in the shear stress produced a 29% reduction in the speed achieved over long times. Analysis using Eq. 2 produced values of v0 = 1.21 nm/ns and ac = 76.28 nm/ns2, with time limits on the transition from period I to period II at 0.55 ns and an upper limit for period II of 2 ns.
Further reductions in τ increased the duration of period I and reduced the magnitude of the constant velocity of period III. For τ = 221.02 MPa, a 10% reduction from the original value, period I was extended to 2.6 ns (Fig. 9), and the constant speed for period III was reduced to 52.97 nm/ns. Under this τ condition, the first stick and slip events (inset plots of Fig. 9) occurred in period III, lasting ~ 0.02 ns. At higher τ (> 40% PN), the layers displaced with continuous sliding, and it was expected that displacements at lower velocities will induce discontinuous displacements due to the presence of conformations that temporally exhibit higher cohesivity [46–48]. It is expected that in simulations under a forced constant velocity (at the same PN and equivalent average τ), the stick and slip behavior will be more periodic due to periodic lapses where the conformations can reorganize and increase cohesivity [29, 43, 49, 50].
As the shear stress was reduced, the adhesion periods became more frequent. The lowest shear stress at a net displacement between the two layers corresponded to τ = 110.51 MPa (0.20 PN). In Fig. 10, the typical displacements that occurred at this low τ are shown, which are characteristic of the single slip regime, where the layers did not move for large periods. Instead, they remained adhered and slid over each other for very short periods. From the probability distributions of the absolute displacements (inset of Fig. 10), the slides had a constant displacement (~ 0.27 nm), equal to half the separation of the surface hydroxyl groups in the sliding direction. This indicates that two types of slip occurred, one that brings the terminal methyl groups from the lubricant molecules into contact, followed by a second type of slip that separates them. The relatively long time (few ns) between slips is necessary to allow the end groups of the chains to pass and laterally surround each other. Similar behavior has been observed in alkyl chains attached to silica layers [51].
The profile for τ = 110.51 MPa does not start at time zero because for shear stresses less than 40% of the normal load, the time for the simulation to reach equilibrium is very long and computationally prohibited, so initial conformations correspond to stationary states at 40% of the normal load. At this τ, the adhesion periods were longer than the slip periods (single slip regime), in contrast to higher τ, where multiple slips (and single slips) were present, and for higher τ, adhesion periods disappeared (smooth sliding) [17]. A linear regression for the τ = 110.51 MPa profile showed a very low velocity, ~ 0.128 nm/ns. Due to the procedure used, it is unknown what happened in the previous regime where the layers were displaced at constant acceleration.
At the shortest τ used in this work (82.88 MPa, 0.15 PN), in 10 ns of simulation, the layers moved slowly, but remained stuck together, moving in the same direction, unlike at higher τ, where the layers moved in the direction of the applied τ. In this way, the minimum τ to overcome the static frictional forces between the layers was found between 15% and 20% of PN (552.56 MPa), which corresponds to the static coefficient of friction, µs, which is between 0.15 and 0.20. The same behavior has been found experimentally: the µs for silica with n-pentanol as lubricant has been measured as between 0.15 [2] and 0.20 [52] at room temperature. Asay et al. [2] have found this coefficient by studying the rotation of quartz balls on silica surfaces treated with perfluoro-alkylated chains, with n-pentanol molecules adsorbed on the surface (95% partial pressure) at room temperature, PN = 115 MPa, and speed of 1.5 mm/s. The coefficient obtained by Yau et al. [52] was obtained by studying the rotation of borosilicate balls with adsorbed n-pentanol, PN = 199 MPa, at room temperature and speeds between 0.11 and 0.33 m/s (nm/ns). Lubrication with n-pentanol molecules considerably reduced the frictional forces naturally found between silica surfaces, which can be passivated with hydroxyl groups or hydrogens.
The results found in this work are summarized in Fig. 11: we show the dependence between the applied τ and the stationary speed obtained from the absolute displacements between the layers. A similar behavior between the shear stress and the speed has been obtained in experiments of sliding contact using steel at 540 oC, and vapor phase lubrication using acetylene [53]. The results follow a logarithmic dependence between the variables across the range of speeds studied in this work; this is consistent with the behavior predicted by models for certain regimes, such as the thermal Prandtl–Tomlinson (PTT) model [54–56], which is specific to single slip displacements. The magnitude of the average velocities obtained in this work by fixing the shear stresses is much higher than those studied experimentally in tribological devices and atomic force microscopy experiments. These experimental studies probably are in the single slip regime [57–59] but the high velocities found in this work are similar to the sliding conditions of the moving parts of NEMs [1]. At very low speeds, such as those employed in experiments, slip events are less frequent, and the corresponding simulations require prohibitive amounts of computation time. Even so, the single slip regime was observed in this work at the lowest shear stress studied (110.51 MPa), which produced an average speed of 0.128 nm/ns.
In the range of velocities between 0.128 nm/ns and the following higher speed studied in this work (0.630 nm/ns), there was a transition between the single slip and the multiple-slip regimes. This transition occurred near the characteristic speed, va, found in simulations for the sliding of Pt tips on gold surfaces [59]; va is the limit velocity where the thermal energy assists single slip events. Above va, the frictional forces approached a plateau in the PTT model [55]. The value for va in the Pt/Au simulations at constant speed is higher than that in the corresponding AFM experiments. The difference was attributed to the failure of the simulations to accurately reproduce the mechanical response of the experimental cantilever [59]. The same may also be true for our work due to the small thickness and small area of the simulated layers. Experimental studies at high speeds are difficult to perform, partly because achieving high speeds requires equipment that generates resonant frequencies due to mutual alignment of moving parts within, which impacts the accuracy of the measured tribological properties [60]. The frequency of the oscillations that occur naturally in the shear stresses and normal loads of experiments performed at constant speed is particularly affected by the resonant frequencies of the equipment.
In the multiple-slip regime shown in Fig. 12 for τ = 238.1 MPa (0.25 PN) at 300 K, the multiple slips could be the result of various single slips occurring after very short stick periods (~ 0.02 ns), or they could truly be multiple slips, occurring after single slips and long stick events. Multiple slips are not continuous; they are rare events that occur after a series of single slip events.
For systems at different τ conditions, the distribution of the atomic positions along the normal axis to the silica surfaces was used to compare them. To highlight the differences, the systems using the smallest (20% PN) and largest (50% PN) τ conditions were compared, they correspond to regimes of single slip sliding and continuous sliding, respectively. We first focused on the contact created by the hydrogen bonds between the hydroxyl groups of the lubricant and the layer. Figure 13 shows the temporal profiles of the average position of the oxygen atoms in the hydroxyl groups belonging to the silica surface (averaged for all the oxygen atoms in each layer).
The system under the highest τ had a slightly expanded volume in the normal direction occupied by the lubricant molecules of ~ 0.58 Å, while maintained the area occupied by the lubricant layers. This expansion in the normal direction corresponds to half the van der Waals radius of a hydrogen atom (1.20 Å) [61]. This slight increase translated into an increased average volume occupied by each n-pentanol molecule of ~ 7.30 Å3, large enough to accommodate the van der Waals volume of a single hydrogen atom. However, due to the compact localization and orientation of the n-pentanol chains, that space could only be accessed by the hydrogen of the hydroxyl group at the end of the n-pentanol chains or the hydrogens of the methyl group at the other end. The additional space gives those hydroxyl or methyl groups more freedom to move, not only in the normal and sliding directions but also in the lateral direction, facilitating the sliding of the layers [20].
The temporal profiles of the average positions also showed the breathing of the layers (Fig. 13a), which manifested as small oscillations of at most 0.2 Å; they are not due to slipping events because the periods of oscillation are very short, ~ 3.2 ps. The probability distributions of these oscillations are normal (Figs. 13b and 13c); the distribution widened for the profile corresponding to the single slip regime using the shortest τ. The normal distribution of the separations indicates that using constant PN as a barostat in the simulations will produce reliable results [10], which may not be obtained using constant-separation simulations. These oscillations represent the expansions and contractions of the layers that are naturally present in any solid. Their normal distribution probably indicates that the behavior of the hydroxyl groups is not affected by the external forces imposed on the outermost layer of silicon atoms.
The probability distributions of the positions of the oxygen atoms of the hydroxyl groups in the silica surface and the hydroxyl groups in the n-pentanol chains are shown in Figs. 14a–14d using the smallest (20% PN) and largest (50% PN) τ conditions. The mean value of the distribution of the oxygen atoms of the silica layers is the same as the mean value of the average position of the entire layer of hydroxyl groups (Figs. 13b and 13c). The distributions of the oxygen atoms of the n-pentanol chains were shorter and wider than the distributions of the oxygens of the hydroxyl groups of the surface layers, showing how free are these hydroxyl groups from the lubricant to move in the normal direction. The average separation between the mean values of both distributions (surface and chains) did not change under the smallest and largest τ conditions, which indicates that the hydrogen bond forces dominated the interactions between the hydroxyl groups and the shear forces at the contact surface that were transmitted to the hydroxyl groups were not large enough to affect the hydrogen bond network.
Once the lubricant chains adsorb on the silica surfaces, they formed upper and lower monolayers, and sliding occurred at the contact surface formed by the methyl groups of the upper and lower monolayers of lubricant. At the other end of the n-pentanol chains, which form the active contact surface where sliding occurs, the probability distribution of the methyl groups of the n-pentanol chains (Fig. 15) showed normal distributions. The separation of the mean values of the upper and lower monolayers of the lubricant chains under the smallest (20% PN) τ condition was ~ 2.54 Å, while at the largest (50% PN) τ condition, it was ~ 2.70 Å. The largest (50% PN) τ condition increased the separation between the lubricant layers by only 0.14 Å.
The mean values for the rest of the carbons, using the largest (50% PN) τ conditions, are also plotted in Fig. 15. The mean values of the distributions for the methyl carbon and the bonded methylene carbon were very close (difference of 0.55 Å). A typical separation between carbons is ~ 1.54 Å. For the next two methylene groups in the chain, the carbon separations were similar due to the zigzag conformation of the n-pentanol chains. This conformation is preferred by normal alkanes, which probably determined the level of tilting of the chains with respect to the silica surfaces. A chain with an orientation perpendicular to the surface will show regular separations between its constituent carbons.