History-Dependent Stress Relaxation of Liquids under High-Confinement: A Molecular Dynamics Study

Abstract

When liquids are confined into a nanometer-scale slit, the induced layering-like film structure allows the liquid to sustain non-isotropic stresses and thus be load-bearing. Such anisotropic characteristics of liquid under confinement arise naturally from the liquids’ wavenumber dependent compressibility, which does not need solidification to take place as a prerequisite. In other words, liquids under confinement can still retain fluidity with molecules being (sub-)diffusive. However, the extensively prolonged structural relaxation times can cause hysteresis of stress relaxation of confined molecules in response to the motions of confining walls and thereby rendering the quasi-static stress tensors history-dependent. In this work, by means of molecular dynamics, stress tensors of a highly confined key base-oil component, i.e., 1-decene trimer, are calculated after its relaxation from being compressed and decompressed. A maximum of 77.1 MPa normal stress discrepancy has been detected within a triple-layer boundary film. Analyses with respect to molecular morphology indicate that among the effects (e.g., confinement, molecular structure, and film density) that can potentially affect confined stresses, the ordering status of the confined molecules plays a predominant role.

1. Introduction

Liquid under confinement, as a subject, has been studied extensively in diverse fields, such as tribology and lubrication, porous media and biomembranes, and electrochemistry, over the past several decades [1,2,3,4]. When liquids are confined into a nanometer-scale slit, they often exhibit layered structures [5,6], significantly enhanced shear viscosities [7,8,9], and oscillatory solvation forces while being continuously compressed [5,10,11]. Those characteristics, deviating from the bulk, are typically due to perturbation (i.e., the presence of confining surfaces) induced molecular ordering, which is inherently a natural consequence of liquid wavenumber dependent compressibility [12,13]. In the framework of a density-functional theory-based description [12,13,14,15], one may discern that the one-particle density of a confined liquid holds a similar functional form as that of the two-particle density of a bulk liquid. Such a similarity has been confirmed by recent simulation works [16,17].

Compared to the Brownian motion of molecules in bulk liquids, as described by the Rouse model [18], the mobility of a liquid under high-confinement (with only a few molecular layers) can be largely impeded. Instead of random walks, motions of confined molecules become collective (known as molecular ordering or orientational correlation), which leads to prolonged structural relaxation times up to several orders of magnitudes [19,20]. In the case of large surface area to film thickness ratio, e.g., when using surface force apparatus (SFA) [21], the confining/shear rate typically, if not significantly, exceeds the inverse structural relaxation time with an increasing discrepancy as the number of confined monolayers decreases. Effective shear viscosity (${\eta }_{\mathrm{eff}}$) increases correspondingly given such an inverse dependence on the relaxation time (${\tau }_{\mathrm{relax}}$), i.e., ${\eta }_{\mathrm{eff}}\propto {\tau }_{\mathrm{relax}}^{-1}$, beyond which liquid flow obeys Newton’s law. Reversible phase transformation can take place when a liquid is close to a liquid–solid or liquid–glass transformation [4,22,23,24], for which, however, one should be aware that confinement alone [24] cannot be served as the sole explanation since, for example, the effects of thermodynamic anomalies are omitted [25]. In fact, the behavior of a liquid confined within a finite space shall be affected by a series of cross effects including molecular structures and chemical compositions [8,24,26], confining/shear rates [27], the nature of confining walls [8,11,28,29] and dimensions [30], and contaminations [31].

Solidification is not a necessity for a confined liquid to manifest solid-like behavior (e.g., a finite yield point and critical shear stress [32]), particularly when the liquid has complex molecular structures. It has been reported [17,19,33,34] that highly confined liquids can remain fluidity in a quasi-static state, with molecules being sub-diffusive. The heterogeneity of liquids under confinement, in light of the oscillatory film density, allows them to possess non-isotropic stress tensors and thereby to sustain shear and load irrelevant to the surface tension [35]. Our previous work [17], for example, shows that a triple-layer boundary film (made of 1-decene trimer) can sustain normal loads on the order of 100 MPa over sufficiently long times, which requires neither crystallization nor glass transition to occur. Note, merely ‘anisotropy’ cannot guarantee a sustainable load-bearing ability as flow may still occur in response to anisotropic stresses. In fact, for many liquids, such a load-bearing ability builds up rather continuously with a decreasing number of the confined monolayers [5]. Interestingly, we also found that the mobility of confined molecules can rise (against an overall trend) with increasing confinement right after a monolayer molecule is expelled from the central confined region [17]. This observation suggests that, in addition to confinement, confining pressure can matter as well [7]. It also explains why a less dense film is able to endure higher loads than a denser film [34], which is far-fetched to be simply understood as a competition of disorderly statistics of a polymer chain and specific surface-segment interaction.

Drainage dynamics of liquids under confinement stems from the change of confining geometry induced molecular restructuring, which can be history-dependent owing to the prolonged structural relaxation time [27]. As an extension of our previous work [17], the focus herein is to investigate the hysteresis of confined stress tensors after the liquid relaxed from being compressed (from $n+1$ to n monolayers) or decompressed (from $n-1$ to n monolayers) by a blunt tip using molecular dynamics (MD). In the two cases studied, the initial pressure differences between the central confined region and the ambient drive the liquid to flow in opposite directions, of which the aim is to clarify to what extent the stress anisotropy persists and whether the relaxed stresses end up comparable. Connections between the stress status and molecular morphology are established to clarify the origin of confinement-dependent stress discrepancy.

2. Methodology

The simulation model is shown in Figure 1, where the liquid is sandwiched between a blunt ridge and a flat substrate. Here we continue to focus on a branched hydrocarbon, namely 1-decene trimer (C${}_{30}$H${}_{62}$, a snapshot of a single molecule is shown as the inset of Figure 1), as in our previous work [17] since it is an important base-oil component and less likely to order compared to linear alkanes. Both the tip and the substrate are made of iron with a (100) surface and are atomically flat. Using a flat punch instead of a curved tip simplifies the problem and improves statistics, but it can also lead to a slight increase in the free energy barrier when a monolayer is expelled. From sideways along the x-direction, molecules can be squeezed into (or sucked from) a reservoir with a nominally constant chemical potential [36]. The in-plane dimensions of the substrate and the bottom of the ridge are 20.0 × 5.1 nm${}^2$ and 9.9 × 5.1 nm${}^2$, respectively. The region underneath the latter is, by default, defined as the ‘central confined region’. There is a total of 450 liquid molecules in the model system, giving rise to an adsorbing liquid film of an approximately six-monolayer thickness in the absence of a counter surface. The substrate and the bottom of the ridge (in mauve) are strongly adhesive to the liquid molecules, whereas other solid walls (in gray) are repulsive with a cutoff of $\sqrt[6]{2}\sigma$ ($\sigma$ is the distance at which the pairwise potential energy is zero). Atoms from the solid walls remain fixed relative to their center-of-mass velocities throughout the simulations provided that deformation of the solids is negligible.

The inter- and intra-molecular interactions of the liquid are described by a united-atom variant—the transferable potentials for a phase equilibria (TraPPE) [37,38] force-field, in which monomers, i.e., CH${}_3$, CH${}_2$, and CH, are represented by single interaction sites (pseudo-atoms). This force-field can faithfully reproduce experimentally measured thermophysical properties [17,39] and can significantly lower the computational cost compared to all-atom variants [40,41]. The intra-molecular (bonded) interactions contain harmonic energy expressions for bond stretching, bond angle bending, and torsional motions. Meanwhile, the inter-molecular (non-bonded) interactions are governed by 12–6 Lennard–Jones (LJ) potential with a cutoff of 1.0 nm. Long-range Coulomb interactions are not taken into account because the pseudo-atoms are electronically neutral. The Lorentz–Berthelot (L-B) mixing rules [42,43] are applied for differing atom pairs, including both the liquid–liquid and liquid–solid interactions [44]. Detailed force-field parameters and functional forms can be found in References [37,38,44]. Note, the purpose here is prone to construct a boundary condition to mimic macroscopic observations rather than to accurately simulate iron-hydrocarbon interactions at an interface, otherwise the iron surfaces would have been passivated with oxygen [45].

There are two sets of simulations carried out: (1) compress/decompress the liquid at a constant rate, and (2) relax the liquid at fixed separations. Regarding (1), the displacing rate of the indenter is 0.2 m/s (along z-axis), which is the lowest rate used in our previous work [17]. In the case of decompression runs, unloading of the indenter takes place either right after it reaches the lowest position (at which the separation between the tip and the substrate is ∼0.8 nm) or after halting the indenter for 100 ns. Regarding (2), the liquid is relaxed at a fixed separation d (d = 1.68, 1.48, and 1.38 nm) for a time-span of 105 ns after a compression or a decompression run. At each separation, a boundary film with a total number of n = 3 monolayers persists. The relaxation duration is chosen based on the fact that all the confining stresses reach ostensibly stable (quasi-static) states at the end of the simulation. Results reported in Section 3.1 correspond to (1), while those in Section 3.2, Section 3.3, Section 3.4 and Section 3.5 correspond to (2). In both sets of simulations, periodic boundary conditions are applied in the $xy$-plane; in the z-direction, atom movements are confined by the solid walls. The system is maintained at 300 K using a Langevin thermostat with a damping time constant of 100 fs. The thermostat is only applied in parallel to the y-component of the velocities of liquid molecules so as not to interfere with the flow dynamics and molecular diffusion. The simulation time-step is 2 fs. Results reported in this work are each averaged over three independent simulations in terms of different seeds. All the simulations are carried out using the open-source code LAMMPS [46].

3. Results and Discussion

3.1. Solvation Force

Solvation forces ($F_z$), measured as the normal forces on the indenter, are shown as functions of wall–wall separations (d) in Figure 2, where results from the two cases appear contrasting. First, forces from prompt decompression (red squares) are predominantly attractive as opposed to the solely repulsive forces during compression (black circles). Given the opposite flow directions and viscous nature of the liquid, this is mainly because the inertia forces from liquid molecules resist tip motion while the molecules are being either squeezed or sucked. Second, undulation, as an indicator of transition of monolayers from ($n\pm 1$) to n, appears only when the molecules are expelled (during compression). Such collective behavior of molecules in response to the external forces, often described as a ‘layer-by-layer’ transition [16,47], is, however, not detected in the case of decompression even after the system is pre-equilibrated for 100 ns before unloading (blue dashed curve). This simply suggests that the current displacing rate (0.2 m/s) is much too fast for molecules to flow and restructure into analogous arrangements within the confining geometry. Also irrelevant to the lowest tip position, discrepancies between the curves in Figure 2 will remain existing as long as the deforming rate surpasses the inverse of the structural relaxation time of molecules confined. In that sense, the transient solvation forces will be dependent on the moving history of the tip.

3.2. Stress Relaxation

To study stress anisotropy and its sustainability, evolution of the stress tensors ($\sigma$) of liquid at the central confined region is monitored at fixed separations ($d_{1–3}$). Using the virial estimator [48], stress tensor components are calculated right after the indenter halted from either a compression or a decompression run with the number of monolayers evolved from $n\pm 1$ to n (n = 3). The virial estimator includes kinetic energy contribution and virial contribution, in which the latter consists of both the intra- and inter-molecular interactions. Although the virial theorem is known to be rigorous for homogeneous systems in thermal equilibrium [49], it predicts, according to our tests, fairly close mean stresses of a multi-layer film as those from the method of planes (MOP) [50]. Furthermore, the normal stresses predicted using the virial estimator agree reasonably well (the errors are less than 5%) with the mechanical definition of the normal pressure: $P_z=F_z/A$, where $F_z$ denotes the total normal forces on the indenter and A the contact area.

Stress anisotropy, detected on both simple [51] and complex [17] molecules, can be attributed to the heterogeneity of liquids under confinement. Compared to being compressed (as demonstrated in our previous work [17]), the anisotropic characteristics of the liquid after decompression appear less pronounced due to the lack of molecular ordering (discuss next). Figure 3 shows a comparison of normal stresses (${\sigma }_{zz}$) evolution after compression (solid symbols) and decompression (open symbols) runs, where stresses from each pair of cases first rapidly approach each other from the opposite ends, then gradually evolve into nominally stable states after approximately 50 ns. Notably, the discrepancy between the quasi-static pair stresses (at ∼100 ns) rises from trivial to somewhat significant as the confining gap decreases. Specifically, according to the fitting curves, the stress differences read 6.1, 28.6, and 77.1 MPa for the gaps of $d_1$, $d_2$, and $d_3$, respectively. Those discrepancies, as one can imagine, would slowly decay with time until finally disappearing once the corresponding structural relaxation times are exceeded, provided that fluidity is remained within the confined zone [17]. However, the corresponding time-span can easily extend beyond the typical time-scale of MD [8], which is assuredly out of the scope of the current study. Nevertheless, before reaching such a threshold, the ability of a confined liquid to sustain non-isotropic stresses persists and ought to be dependent on both the confinement and the moving history of the tip.

3.3. Radius of Gyration

The aforementioned quasi-static history dependence of stress tensors and earlier claimed stress anisotropy [17] are inherently determined by the detailed confined molecular arrangement, which is not readily accessible by the experimental techniques. To elaborate such correlations, the radius of gyration ($R_{\mathrm{g}}$) of central confined molecules is calculated based on time-dependent MD configurations via $R_{\alpha \alpha }^2=\frac{1}{N_{\mathrm{mol}}}{\sum }_{i=1}^{N_{\mathrm{mol}}}\frac{1}{N_{\mathrm{mon}}}{\sum }_{j=1}^{N_{\mathrm{mon}}}\langle \delta r_{\alpha }^{i,j}·\delta r_{\alpha }^{i,j}\rangle$. In these calculations, only molecules located within the mid-monolayer are taken into consideration for reasons to exclude boundary effects from the confining walls. In the expression, $\alpha$ denotes Cartesian coordinates, $\delta r$ is the relative position between monomer j and the center-of-mass of molecule i it belongs to, $N_{\mathrm{mol}}$ and $N_{\mathrm{mon}}$ are the total numbers of embraced molecules and monomers in every single molecule, respectively, and angle brackets denote an ensemble average. The principle moments of $R_{\mathrm{g}}^2$ at $d_{1–3}$ during the last 25 ns of relaxation are shown in Figure 4.

Within the confining space, liquid molecules are highly deformed and parallel to the confining walls, as opposed to those in a bulk phase of which the three principle components of $R_{\mathrm{g}}^2$ are comparable ($R_{\alpha \alpha }^2\approx 10.8$ Å${}^2$ for 1-decene trimer) [17]. As shown in Figure 4, the in-plane components ($R_{xx}^2$ and $R_{yy}^2$) are at least two-fold as that in the normal direction ($R_{zz}^2$), indicating a sheet-like film structure that complements the observed oscillatory density profile. In the case of compression, the extension of molecules along the flow direction (x-axis) appears more sensitive to the confinement change given a rapidly increased pressure gap between the initial and the final states (see Figure 3). High pressure inducing a strong tendency to flow can ease molecular alignment along the streaming direction and thereby yield stretched molecular appearances. In contrast, the bottleneck effect, particularly at high confinement, would hamper the traveling of disordered molecules from an open space into a narrow channel (decompression case), resulting in a relatively low film density in the central void (discuss next). Under low confining pressure, the inert response of molecular conformation to the confinement change suggests that a much longer time is required for molecules to replenish the central void and become orderly arranged thereafter. This explains why there is no sign of collective behavior (e.g., undulations on the solvation force curve) in the course of continuous decompression, as shown in Figure 2.

3.4. Mean-Squared Displacement

Diffusion of central confined molecules at $d_{1–3}$ is characterized via measuring the in-plane mean-squared displacement (MSD) of monomers in the mid-confined-layer: MSD = $\frac{1}{N_{\mathrm{mol}}}{\sum }_{i=1}^{N_{\mathrm{mol}}}\frac{1}{N_{\mathrm{mon}}}{\sum }_{j=1}^{N_{\mathrm{mon}}}{|{\mathit{r}}_{i,j}\left(t\right)-{\mathit{r}}_{i,j}\left(t_0\right)|}^2$, where ${\mathit{r}}_{i,j}\left(t\right)$ is the position of the j-th monomer from the i-th molecule at time t ($t\ge t_0$ = 80 ns). A comparison of the results between compression and decompression runs is shown in Figure 5. For all the cases studied, the confined liquid remains fluidity with molecules being sub-diffusive given that the smallest mean displacement (∼4 Å) of monomers over a 25 ns time window exceeds the typical diffusion rates of polymers adsorbed on solid surfaces [52] as well as atoms within crystalline structures, even those located along the high-diffusivity paths (dislocations, grain boundaries, and free surfaces) [53]. Clearly, at a given separation, molecules relaxed from decompression are more diffusive compared to those after compression, particularly evident in the case of high-confinement (see the blue curves in Figure 5). It is the discrepancies between film densities and thereby between molecular morphologies that render the mobility of confined molecules history-dependent. It ought to clarify that higher molecular mobility in the central confined region does not necessarily refer to an easier incoming flow at the entrance of the confined channel (i.e., the edges of the tip) but the opposite instead. Moreover, as pointed out in our previous work [17], diffusivity of confined molecules may not necessarily hold a monotonic function of confinement given the fact that it is the confining pressure that plays an essential role. Hereby, this claim is partially endorsed by the results shown in Figure 5, where molecules at a higher confinement (blue dashed curve) can be more diffusive than those at a lower confinement (red curves for example). Note, at $d_3$, neither of the cases are in equilibrium states (at least the one from decompression is apparently not); otherwise, the two curves will be comparable with molecules presumably being less diffusive compared to those at $d_2$. However, after relaxing the system for an additional 100 ns, diffusion of molecules at $d_3$ (from decompression) remains momentarily higher than those at $d_2$.

3.5. Density Profile

Figure 6 shows confined liquid density profiles ($\rho$) along the normal direction (z-axis) with a bin size of 0.5 Å. At each separation, three molecular layers, evolved from t = 0 to 100 ns, can be identified (according to the oscillation periods). Evidently, wider wall–wall separation ($d_1$) inhibits molecular transfer to a lesser degree given the analogous density profiles ending up after 100 ns relaxation (Figure 6d). While, at high confinement ($d_3$), the notable discrepancy between the density peaks of the mid-layer (Figure 6f) shall be the origin of the aforementioned dissimilarities in molecular morphology and stress tensors at rest. That being said, the lower the density peak, the less pronounced molecular ordering and thereby the lower confining stresses. Note, the correlation between confining stresses and film density is not ubiquitous. As an example, the average film density in the case of d = 1.68 nm at 100 ns (∼0.97 g/cm${}^3$) is 24% higher than that of the bulk liquid density (0.78 g/cm${}^3$) [17] at ambient conditions, whereas the normal stresses appear attractive (the negative stress at 100 ns reads –27.9 MPa). Essentially, what determines the stress tensors of a liquid confined within a limited space is not merely whether the molecules are densely packed (mean film density) but rather to what extent they are orderly arranged (oscillatory peak height).

4. Conclusions

The heterogeneity of liquids under confinement enables them to sustain non-isotropic stresses that are irrelevant to the surface tension. In this work, stress relaxation of a key base-oil component, 1-decene trimer, under nanoscale confinement was studied by means of molecular dynamics. In terms of a triple-layer boundary film within a fixed slit, a stress discrepancy was identified while the liquid relaxed from tip compression and decompression in nominally quasi-static states. Such a discrepancy, increasing with confinement, stems from the differing ordering status of molecules in the central confined region that is dominated by the confining pressure. Although highly relevant, factors such as confinement and film density are rather extrinsic. The comparably lower normal stresses in the case of decompression can be attributed to the impeded flow of disordered molecules from an open space into a narrow channel, and thus, a much longer time is required for molecules to evolve into an orderly arranged film structure. Such prolonged structural relaxation times can induce hysteresis of molecular restructuring in response to the external forces and thereby yield the stress tensors history-dependent. Findings from this work, along with our previous study [17], may complement the understanding of dynamic and mechanical drainage responses of liquids with non-granular molecular structures under small-scale confinement.