A non-empirical free volume viscosity model for alkane lubricants under severe pressures

TL;DR

This paper develops a non-empirical, microscopic free volume viscosity model for alkane lubricants under severe pressures, combining molecular dynamics data with theoretical concepts to accurately predict viscosity without empirical parameters.

Contribution

It introduces a parameter-free viscosity model based on equilibrium molecular dynamics and microscopic ensemble averages, applicable to high-pressure alkane lubrication.

Findings

Accurately predicts viscosity and diffusion coefficients at high pressures and temperatures.

Provides a microscopic basis for viscosity modeling without empirical fitting.

Enhances understanding of alkane lubricant behavior under severe conditions.

Abstract

Viscosities $\eta$ and diffusion coefficients $D_s$ of linear and branched alkanes at high pressures $P$$<$0.7 GPa and temperatures $T$=500-600 K are calculated by equilibrium molecular dynamics. Combining Stokes-Einstein, free volume and random walk concepts results in an accurate viscosity model $\eta(D_s(P,T))$ for the considered P and T. All model parameters (hydrodynamic radius, random walk step size and attempt frequency) are defined as microscopic ensemble averages and extracted from EMD simulations rendering $\eta(D_s(P,T))$ a parameter-free predictor for lubrication simulations.

Equations

$$D_s=k_{B}T/(n\pi \eta R_h)$$

$$R_h=\sqrt{⟨a⟩/\pi }$$

$$D_s=D_0\exp (-v_c/v_f).$$

$$D_s(a,r_c)=r_c^2/(6\Delta t)=r_c^2/(6{\tau }_0)\exp (-r_{c}a/v_f).$$

$$D_s\approx ⟨r_c{⟩}^2/(6{\tau }_0)\exp (-⟨r_c⟩⟨a⟩/v_f).$$

$$⟨\mathbf{e}(0)\cdot \mathbf{e}(t)⟩=C{\text{e}}^{-t/{\tau }_{\alpha }}+(1-C){\text{e}}^{-(t/{\tau }_{\beta }{)}^b}$$

$$\log \frac{\eta (T,\rho )}{{\eta }_0(T)}=\frac{⟨r_c⟩⟨a⟩}{v_{mol}}({\rho }_m/\rho -1{)}^{-1}$$

Full text

Knowlegde-based design and optimization of liquid lubricants require a quantitative modelling of their rheological properties under relevant tribological conditions Bair (2007). For instance, lubricants in roller element bearing and gear applications are subject to pressures of the order of GPa Spikes (1994). Traditional empirical viscosity models (such as Barus or Roelands equation) fail to describe $\eta(P)$ over the relevant pressure range Bair and Kottke (2003) indicating that improved viscosity models for the extreme pressure regime Bair et al. (2013); De la Porte and Kossack (2014) should be based on physical insights Vergne and Bair (2014). A promising approach employs the Stokes-Einstein relation Einstein (1905)

[TABLE]

that connects the viscosity with the self-diffusion coefficient $D_{s}$. Here, $R_{h}$ denotes the hydrodynamic radius and $n$ lies between the slip and no-slip hydrodynamics limits, 4 and 6. The applicability of Eq. (1) on a microscopic level has been theoretically motivated and is well established under normal conditions Zwanzig and Bixon (1970); Hansen and McDonald (2013), also for non-spherical molecules if $R_{h}$ is considered a free parameter. However, a breakdown of Stokes law has been observed in various dense liquids, including molecular glass formers Bordat et al. (2003); Brillo et al. (2011). The mechanism of this breakdown is still subject to extensive research, mostly focused on densification obtained by supercooling Charbonneau et al. (2014); Henritzi et al. (2015); Kawasaki and Kim (2017). This raises the question whether Eq. (1) remains valid for liquids which are densified by pressurisation instead of cooling.

In this letter, EMD simulations are utilized to validate Eq. (1) for linear and branched alkanes (constituents of ordinary lubricant base stock) for a pressure and temperature range representing typical tribological high load applications. We suggest a microscopic defintion of

[TABLE]

by introducing an EMD averaged molecular cross section $\langle a\rangle$. This represents an important step towards a parameter-free structure-property relationship (SPR) relateing molecular structure with macroscale viscosities.

Employing Eq. (1) necessitates an additional SPR for $D_{s}$. Here, we utilize the free volume (FV) concept Doolittle (1951); Cohen and Turnbull (1959)

[TABLE]

with the mean FV per molecule $v_{f}$ and the critical volume $v_{c}$. Note, that since $v_{f}$ is determined by the lubricant density $\rho$ an equation of state $\rho(P)$ is required to arrive at a pressure dependent viscosity law.

Although widely used for soft matter systems Almeida et al. (1992); Javanainen et al. (2010); Vrentas and Vrentas (2012), the FV concept is being challenged due to concerns about the relevance of FV compared to energetic effects and about the physical interpretation of the free parameters $D_{0}$ and $v_{c}$ Vrentas and Vrentas (2003); Falck (2005); Betancourt et al. (2015); Berthier and Biroli (2011); Jadhao and Robbins (2017). We show in the second part of this letter that Eq. (3) can be applied to our alkane lubricants and suggest microscopic definitions of $D_{0}$ and $v_{c}$. The latter relies on the observation that self-diffusion can be considered a random walk of a molecule’s center of mass (COM) with step length and attempt frequency determined by EMD simulations.

Two linear alkanes and three poly-$\alpha$-olefines (PAO) (for structures see Fig.1a) are modelled with the all-atom optimized potentials for liquid simulations (AA-OPLS) Jorgensen et al. (1996). The EMD of these lubricants is simulated within a constant volume subject to periodic boundary conditions for densities ranging roughly from $470$ to $850$ kg/m3. Time integration is performed employing the LAMMPS software suite Plimpton (1995); http://lammps.sandia.gov with timestep $0.5$ fs, and a Nosé-Hoover thermostat with relaxation time $0.1$ ps Allen and Tildesley (1989); Frenkel and Smit (2001). Viscosities $\eta$ are determined via the Green-Kubo formalism and self-diffusion coefficients $D_{s}$ via the mean squared displacement (MSD)Hansen and McDonald (2013) (see supp. Figs. S1/2).

As shown in Fig. 1b, results for $\eta$ and $D_{s}$ vary over 3 orders of magnitude and are fully compatible with Eq. (1) (assuming slip boundary conditions $n$=4 Hansen and McDonald (2013); Zwanzig and Bixon (1970)). A parameter-free quantitative agreement is achieved by introducing the hydrodynamic radius via Eq. (2) as follows. Since the Stokes drag on macroscopic solid objects with slippery surfaces scales with the object’s cross section in the direction of a displacement Leith (1987), we calculate a directional molecular cross section $a$ as indicated in Fig. 1c. To each molecule a volume $v_{mol}$ is assigned using a coarse grained hard sphere approach based on a CHX ($X=1,2,3$) united atom representation Martin and Siepmann (1998) (see details in supp. Fig. S3). Then the molecule is displaced over a short distance $\varepsilon$ and its effective cross section area is defined by $a=\Delta V/\varepsilon$, where $\Delta V$ is the newly occupied volume. Finally, the mean cross section is obtained by $\langle a\rangle=\int aq(a)\textnormal{d}a$, where $q(a)$ is the probability for the molecule to have a configuration with cross section $a$. Interestingly, $\langle a\rangle$ scales linearly with the molecule size, despite the different morphologies (i.e. number of branches in alkanes - see inset in Fig. 1b). This scaling can be rationalized by a cylindrical shape estimate of long alkane chains, neglecting the contribution of chain ends and knots (supp. Fig. S4).

After having established a parameter-free relation between $\eta$ and $D_{s}$, we now focus on the diffusive motion of the alkanes. Fig. 2a shows part of a C10-trimer COM trajectory with a behaviour which is characteristic of a caging effect. The COM position oscillates within a compact volume due to confinement by the neighboring molecules (snapshots 1,3,5 of Fig. 2a). Elementary diffusion steps (EDS) take place via occasional irreversible translations (indicated by red arrows in snapshots 2,4,6 of Fig. 2a). The FV ansatz leading to Eq. (3) assumes that the probability for an EDS is given by $p(v_{c})=\exp(-v_{c}/v_{f})$. Here, $v_{c}$ is the critical void size in the cage formed by a molecule’s neighbours allowing for an irreversible COM jump. Following the simple argument that this void has to accommodate the molecule, the critical volume $v_{c}$ is expected to be of the order of the hard core molecule volume $v_{mol}$. Note, that $p(v_{c})$ depends parametrically on the mean free volume per molecule $v_{f}=v-v_{mol}$, where $v$ denotes the molar volume Cohen and Turnbull (1959).

Indeed, the self diffusion coefficient follows the form Eq. (3), as shown in Fig.2b, where lines are best fits for the $600$ K data. Fig. 2c displays the dependence of the fit parameters $\tilde{v}_{c}$ and $\tilde{D}_{0}$ on the size of the molecules. Surprisingly, the critical volume is about $3$ times larger than $v_{mol}$ in contradiction to the simple argument stated above. An alternative interpretation of $v_{c}$ is based on the following consideration. To perform an EDS, a molecule needs to move from its cage center to a void in the cage wall. The necessary critical volume for this displacement over the cage size $r_{c}$ is then $r_{c}\cdot a$ with $a$ the molecule’s cross section. For constant $a$ and $r_{c}$, the diffusion process could then be pictured as a random walk with stepsize $r_{c}$ and step frequency $1/\Delta t$. The latter is the product of an attempt frequency $1/\tau_{0}$ and the success probability $\exp(-r_{c}a/v_{f})$ resulting in

[TABLE]

However, both $a$ and $r_{c}$ depend a priori on the molecules’ configurations with respect to the direction of the EDS. Fluctuating shapes and distances in molecular fluids result in a probability distribution $q(a,r_{c})$ for $a$ and $r_{c}$ that determine the diffusion coefficent $D_{s}(a,r_{c})$ in a certain direction and thus the total diffusion coefficient is $D_{s}=\int D_{s}(a,r_{c})q(a,r_{c})\,\textnormal{d}a\,\textnormal{d}r_{c}$. In the following, we demonstrate by sampling $q(a,r_{c})$ over all possible configurations and orientations that the first moments $\langle a\rangle=\int aq(a,r_{c})\,\textnormal{d}a\,\textnormal{d}r_{c}$ and $\langle r_{c}\rangle=\int r_{c}q(a,r_{c})\,\textnormal{d}a\,\textnormal{d}r_{c}$ dominate the diffusion process:

[TABLE]

For a given configuration, Eq.(4) implies a direction dependent stepsize and success probability for an individual EDS. To study this anisotropy we consider an auxiliary system of preferentially oriented alkanes immersed in a bath of unconstrained molecules. This artificial test situation is realized for n-hexadecane at two different densities by applying opposing external forces $\pm 0.05$ eV/Å to the head and tail carbon atoms of $2.6\%$ of the molecules (Fig.3a). The resulting preferential orientation leads to a permanently anisotropic cross section $\langle a(\theta)\rangle$ (Fig.3b) and cage radius $\langle r_{c}(\theta)\rangle$ (red dots in Fig.3c), where $\theta$ denotes the angle between the applied forces and the EDS direction. Here, $\langle a(\theta)\rangle$ was calculated as previously defined (see Fig. 1c) for a given direction $\theta$. Lacking an unambiguous definition of $\langle r_{c}(\theta)\rangle$, a pragmatic estimate was based on the direction dependent radial distribution function $g_{COM}(\theta,r)$ of the molecules’ COM (Fig. 3c) via $g_{COM}(\theta,\langle r_{c}(\theta)\rangle)=1$. In the same spirit, the isotropic cage size $\langle r_{c}\rangle$ (circle in Fig. 3c) was determined from the isotropic radial distribution function $g_{COM}(r)$ (bold line in Fig. 3c). Note, that this value is close to the result obtained via $\langle r_{c}\rangle=\int\langle r_{c}(\theta)\rangle d\,\cos{\theta}$.

As expected the preferentially oriented molecules exhibit a pronounced anisotropy in the MSD $\langle r^{2}(\theta)\rangle$ (Fig.3d) and consequently in the diffusion coefficient $D_{s}(\theta)$. Interestingly, the mean diffusion coefficient $D_{s}=(D_{s}(\theta=0)+2D_{s}(\theta=\pi/2))/3$ is equal to the isotropic diffusion coefficient of the unperturbed molecules, confirming that diffusion is given by an average over all directions with respect to the main axis of a hexadecane molecule. Most importantly however, $D_{s}(\theta)$ scales with the respective static structure properties as predicted by Eq.(5) validating the applicability of a FV ansatz on the microscopic level of an EDS (Fig. 3e).

Comparing Eq.(3) with (5) leads to $D_{0}=\langle r_{c}\rangle^{2}/(6\tau_{0})$ and $v_{c}=\langle r_{c}\rangle\langle a\rangle$. Indeed, applying the above structure evaluation for $\langle a\rangle$ (inset in Fig. 1b) and $\langle r_{c}\rangle$ (Fig. 4a) to the unperturbed systems of all 5 fluid types reveals that the product $v_{c}=\langle r_{c}\rangle\langle a\rangle$ agrees well with the fitted critical volumes $\tilde{v}_{c}$ (Fig. 4a inset). Both, $\langle a\rangle$ and $\langle r_{c}\rangle$ are only weakly dependent on density and temperature (supp. Figs. S5/6) and can be conveniently estimated from an EMD simulation for a single $\rho$ and $T$.

The remaining free parameter $\tau_{0}$ (time between random walk attempts) can be interpreted as the time scale for structure decorrelation in the molecule/cage system. The connection between structural relaxation, diffusion and viscosity is subject of ongoing research Ma et al. (2019), and fully unravelling the underlying mechanisms goes far beyond the scope of this work. Nevertheless, as a starting point we consider the time autocorrelation function $\langle{\bf e}(0)\cdot{\bf e}(t)\rangle$ of the molecules’ end-to-end vector orientation ${\bf e}(t)$ (see Fig. 4b) to quantify the intramolecular structure decorrelation. For long hydrocarbon chains, $\langle{\bf e}(0)\cdot{\bf e}(t)\rangle$ is well described by a double exponential function Morhenn et al. (2012)

[TABLE]

with two separate characteristic decay times $\tau_{\alpha}$ and $\tau_{\beta}$. On the one hand, a long time decay is observed on the time scale of the diffusion process $\tau_{\alpha}\approx\Delta t=r_{c}^{2}/(6D_{s})$, ranging from $10$ ps - $1$ ns. On the other hand, $\tau_{\beta}$ is of the order of $1-10$ ps and is insensitive to the fluid density (within statistical uncertainties, see supp. Fig. S7). This $\beta$-relaxation time fits well with the expected attempt frequency for the random walk $\tau_{\beta}\approx\tau_{0}=r_{c}^{2}/(6D_{0})$ as illustrated in Fig. 4c, which suggests it as a good measure for the relevant structure decorrelation on short times.

The presented results are further validated in a series of scaling tests with modified model parameters for both intra- and intermolecular interactions (supp. Fig. S8). In particular, the strength of nonbonded interactions has little influence, but a scaling of the atomic radii $\sigma_{LJ}$ results in strong variations of the diffusivity caused by the exponential term in Eq. (5). Moreover, the prefactor $D_{0}$ is sensitive to variations of the energy barrier for bond rotation, which influences intramolecular relaxations. All scaling tests are also in quantitative agreement with the Stokes-Einstein relation, supporting our definition of a hydrodynamic radius in Eq. (2).

Finally, combining Eqs.(1), (2) and (5) the viscosity can be expressed as a function of density

[TABLE]

with $\rho_{m}=M/v_{mol}$ the maximum hypothetical density for zero free volume ($1-v_{f}/v=\rho/\rho_{m}$; supp. Fig. S4). Apart from the density, the r.h.s. contains only equilibrium structure properties, namely the molecules’ volume $v_{mol}$, mean cross section $\langle a\rangle$ and mean next neighbor distance $\langle r_{c}\rangle$. The temperature dependence of $\eta$ enters via $\eta_{0}(T)=1.5k_{B}T\tau_{0}(T)/(r_{c}^{2}\sqrt{\pi\langle a\rangle})$. By employing Eq.(7) and identifying $\tau_{0}$ with $\tau_{\beta}$ a parameter free rescaling of the simulated viscosities can be established (Fig. 4d). This scaling law can also be applied to experimental high $T$ and $P$ viscosity data for *n-*dodecane Caudwell et al. (2004). We find a good agreement of the experimental data with our parameter-free viscosity model Eq. (7) (inset of Fig. 4d).

To conclude, a combination of basic random walk and FV theory fully describes the self diffusion mechanism of long alkane chains, linear and branched alike, in the high $T$ and $P$ regime. A crucial part of the presented work is the introduction of the mean cross section $\langle a\rangle$ and mean cage size $\langle r_{c}\rangle$ as novel molecule shape parameters. While $\langle a\rangle$ establishes a quantitative link between viscosity and self diffusion via the Stokes-Einstein relation, $\langle a\rangle$ and $\langle r_{c}\rangle$ allow for a parameter-free density scaling of both transport coefficients. The viscosity model can be directly implemented in density-based Reynolds-solvers Elrod (1981) and will contribute to a new cutting-edge simulation tool for tribological applications. The proposed shape parameters also opens new possibilities to quantify the role of molecular structure on rheology in anisotropic situations, such as shear thinning Liu et al. (2017); Ingebrigtsen and Tanaka (2018) or in nanometer-thin boundary lubrication films Rosenhek-Goldian et al. (2015). Our approach might also be useful for other soft materials, such as self assembled membranes Almeida et al. (1992); Javanainen et al. (2010), polymer–solvent systems Vrentas and Vrentas (2012) or adsorbates in nanopores Falk et al. (2015).

The authors gratefully acknowledge funding by the industrial partners of the MikroTribologie Centrum $\mu$TC (Karlsruhe, Germany), computing time within project HFR14 at NIC Jülich and useful discussions with L. Joly, S. Kapfer and L. Bocquet.