Unraveling the pressure-viscosity behavior and shear thinning in glycerol using atomic scale molecular dynamics simulations

In order to increase the usage and explore new applications of glycerol as a replacement for fossil-based lubri-cants its properties needs to be known at the fundamental level. In this study, the viscosity of pure glycerol at high pressures and strain rates has been investigated using of molecular dynamics (MD) simulations, utilizing both the Green-Kubo (GK) formalism and the SLLOD algorithm. Although the viscosity acquired by the GK method is in agreement with the corresponding experimental values at low pressure, a significant distinction was identified between the viscosity obtained by the GK method and the experimental values at higher pressures (P > 0.5 GPa). This results in a clear difference between the viscosity-pressure coefficient attained by the GK method and the corresponding experimental value. The SLLOD method using a non-equilibrium MD (NEMD) platform was exploited to take into account the simultaneous effects of strain rate and pressure on viscosity. As a result, the pressure-viscosity coefficient acquired by the SLLOD algorithm approaches the experimental value. By combining the experimental outputs for viscosity at low strain rates ( ˙ γ < 10 4 s (cid:0) 1 ) with the SLLOD outputs at higher rates ( ˙ γ > 10 5 s (cid:0) 1 ), the evolutions of glycerol viscosity with pressure and strain rate were ultimately achieved. Implementing this computational platform depicts the shear thinning process in pure glycerol in a wide range of pressures and strain rates.


Introduction
The release of oil-based lubricants into the environment causes many environmental hazards every year, including impeding the growth of plants, endangering the aquatic life cycle, and polluting groundwater [1,2].Industrial lubricants can contaminate a wide area and affect vast domains of the environment by leaking into waterways [3].Regarding fossil oils, considering that 1% of the total consumption of oil is used for lubricants, we must terminate the usage of these oils and replace a vast amount of lubricant base oils with green alternatives.Hence, there is an ever-increasing need and desire to switch to environmentally friendly lubricants, including vegetable oils [4,5].However, weak thermal stability and narrow viscosity range are among the drawbacks of vegetable oils [6,7].Glycerol aqueous mixtures, in addition to more favorable physical properties at low temperatures, are biocompatible, which make them a suitable alternative for vegetable oils [8,9].Since glycerol is one of the by-products of conventional biodiesel production, which is a fast growing industry, and a way to increase its profitability, would be to effectively use glycerol in innovative applications [10,11].Additionally, glycerol has exhibited promising tribological properties, including a very low friction coefficient in the elastohydrodynamic lubrication (EHL) regime [12][13][14][15][16][17].It is important to keep in mind that glycerol exhibits hygroscopic qualities, making it impractical for direct use as a lubricant in its pure state.Moreover, the viscosity of pure glycerol tends to be excessively high for the majority of lubrication applications, thus making a 15 wt% aqueous solution a more suitable choice [18].Shi et al. demonstrated that under certain conditions, the friction coefficient of glycerol aqueous mixtures decreases with the water content up to 15 wt % [19,20].During a wear test, it became evident that this specific water concentration marked the threshold between low and high wear rate in hydrodynamic and boundary lubrication regimes [21].Nonetheless, the threshold is significantly influenced by various conditions, including load, speed, surface roughness, and temperature.
Viscosity is the determinant property of a lubricant from the perspective of load-bearing capacity.Depending on the surface mechanical load, the relative sliding velocity of the substrate, and the viscosity, lubrication can be performed by different mechanisms, including boundary, mixed, and full film (elasto)hydrodynamic.When the applied mechanical load to the lubricant is low enough, the viscosity practically does not experience any alterations, while under heavy mechanical loading, it is highly dependent on pressure and shear strain rate.Moreover, with the addition of 15 wt% water, the viscosity of glycerol experiences a notable decrease, shifting from 1410 cP to 109 cP at 20 • C [18].Figures S1 and S2 in supplementary information illustrate the graphs depicting variations in the viscosity of glycerol aqueous solutions with temperature and the corresponding changes in the viscosity-temperature coefficient with concentration, respectively.Consequently, glycerol presents not only a temperature-viscosity coefficient to consider but also a humidity-viscosity coefficient that warrants attention.Significant viscosity reduction at high strain rates, defined as shear thinning, occurs in polymer-enhanced lubricants due to both polymer chain alignment and molecular bond dissociation [22][23][24].The high-pressure viscometers used to diagnose the shear thinning behavior have their own challenges, among these is the considerable effect of the heat generated at high shear rates on the viscosity.In experimental methods, the thickness of the lubricant layer is measured in the EHL regime, and the average viscosity value is evaluated using the regeneration of the shear thinning results together with the deterministic relation in the EHL regime [25][26][27].However, high pressure falling body viscometers have been utilized as a more accurate method than measuring film thickness at higher pressures [28].
Molecular dynamics (MD) simulation has proven to be a reliable and affordable computational framework that can be exploited to calculate the properties of fluids, including lubricants.This method can be utilized to calculate the physical properties of the lubricant in extreme conditions when it is not technically possible or economically viable to implement an experimental measurement procedure [27,29].Generally, for a successful MD simulation framework, a molecular model with a suitable size that can clarify the atomic-scale alterations taking place in an elastohydrodynamic configuration is needed [30][31][32][33].However, the typical limitations related to computational methods still exist for MD simulations.In a limited simulation time (t < 100 ns), thermal noises of the system can affect the physical properties outputs.Therefore, for shorter timespans in MD simulations, high strain rates (γ > 10 8 s − 1 ) are usually used to suppress the thermal noise effects of the system.In contrast, lower shear rates (γ < 10 5 s − 1 ) are usually used in experimental studies.It is noteworthy that the effects of thermal noises on the outputs also depend on the dimensions of the system of study.As expected, the effects of thermal noises will be intensified by decreasing the size of the supercell describing the lubricant using period boundary conditions.Typical supercell widths in MD simulations (<100 Å) can be significantly thinner than an actual elastohydrodynamic film while increasing the dimensions of the system will impractically boost the amount of required computational resources.Therefore, it is necessary to make a trade-off between designing a small enough system to reduce the computational costs and overcome thermal noises and large enough to accurately describe what is happening in a real elastohydrodynamic film [34,35].On the other hand, under extreme conditions (high pressure/ temperature/shear strain rate), a physical parameter can be determined through MD computation while other parameters of the material can be obtained using empirical equations that establish a correlation between the total physical characteristics.Muraki et al. demonstrated that the pressure-viscosity coefficient (α) is a determinant parameter of the lubrication performance in the EHL regime [36].Afterward, Bair et al. performed viscosity evaluations up to 1.4 GPa to calculate the viscositypressure coefficient for several traction base stocks and sorted the coefficients based on traction coefficients [37].Wang et al. presented a computational model based on non-equilibrium MD and pressureviscosity coefficient evaluation up to 0.6 GPa to design base stock molecules prior to synthesis [37].Limited studies have been conducted to determine the viscosity of lubricants in the range of high pressures and strain rates (P > 1.5 GPa and γ > 10 6 s − 1 ) [38].But as pointed out, MD simulations can provide a more accurate understanding of the rheological behavior of the lubricant by calculating viscosity values at extreme pressure and strain rate ranges.
This study aims to introduce a computational platform based on NEMD simulation that paves the way for understanding how the viscosity of a typical lubricant is affected by high pressures and strain rates and depicts the shear thinning process in glycerol as a green lubricant.In this study, a molecular model is built to calculate the physical properties of pure glycerol, and the CHARMM36 force field is used to extract the atomic interaction parameters of the system.As the first step, the diffusion coefficient of pure glycerol is calculated and compared with the values reported in previous studies to verify the validity of the force field parameters.Then, we evaluate the viscosity of glycerol in a wide range of pressures and strain rates using the Green-Kubo and SLLOD methods.Calculating the evolutions of viscosity with strain rate and pressure provides a consistent insight into the shear thinning process for a typical lubricant.

Computational methods
This section presents the details associated with the system configuration and the applied MD simulation methods.

Remarks on system setup and equilibrium simulation
MD simulations is a reliable computational framework for the evaluation of the physical properties of lubricants for properties such as viscosity under special conditions such as high pressure and high shear rate.Since the interaction parameters of glycerol are not available among the existing fragments of proteins [39], nucleic acids [40,41], lipids [42], and carbohydrates [43] in the CHARMM36 force field, the generalized expression of the CHARMM36 force field (CGenFF) [44][45][46][47][48] was utilized to compute the atomic interactions of the system including atomic species, partial charge, and bond parameters.The general form of the CHARMM36 force field is given by the following expression [49]: where (k b , k θ , k ω , k ub ) and (b 0 , θ 0 , u 0 , ω 0 ) represent the force constants and reference values for the contribution of covalently bonded pairs, angular bending, Urey-Bradley terms, and improper dihedrals in the energy of the system, respectively.K φ , n, and δ denote the force constant, dihedral multiplicity, and reference phase for the dihedral contribution to the energy of the system.The non-bonded terms include Coulomb and van der Waals (vdW) interactions, where q i and q j are the partial atomic charges of atoms i and j, respectively.ε ij and R min,ij correspond to the Lennard-Jones (LJ) well depth and the distance at the LJ minimum, respectively, while r ij is the distance between i and j particles.Based on the CGenFF platform, a deterministic decision tree V. Fadaei Naeini et al.
was used to assign the atomic species of the CHARMM36 to the atoms in the system.Afterwards, the interaction parameters of the atom types in glycerol are specified by replacing the defined bond parameters of the corresponding atomic species in CGenFF.Moreover, the partial charge assignment was performed by the bond-charge approach developed by Vanommeslaeghe, et al. [46].
A cubic simulation box containing 875 glycerol molecules (12250 atoms) was built with the size of 53.5 Å in X, Y, and Z directions using Packmol [50].The number of glycerol molecules in the periodic cell and the initial system volume were adjusted so that the density of glycerol is as consistent as possible with the reported values [51] in the desired simulation temperature.
The equilibrium stages were carried out under the NPT ensemble, and the Nosé-Hoover barostat was used to maintain the system temperature at T = 300 K and the system pressure to a constant level for each set of simulations [52,53].The relaxation times for the Nosé-Hoover thermostat and barostat were adjusted to 100 fs and 1000 fs, respectively.All the simulation steps were carried out using the LAMMPS package [54,55] with a graphics processing unit (GPU), and the simulation timestep was set to 1 fs.The particle-particle particle-mesh solver [56] was exploited to compute long-range Columbic interactions in K-space.
Three series of simulations, including a equilibrium simulation, a production run to implement the Green-Kubo formalism, and a production run using the SLLOD algorithm, were performed (see below).The final configuration of the system after running the equilibrium stage was used to initiate the simulations to calculate the viscosity based on the Green-Kubo and SLLOD methods.The simulation time for the production run of each step is presented in Table 1.
The following sections provide the implementation details of the Green-Kubo and SLLOD methods to calculate viscosity.

Green-Kubo method
The Green-Kubo (GK) formalism is a fundamental approach exploiting the time autocorrelation function to calculate the most important properties during transport processes in a liquid [57][58][59][60].The shear viscosity can be computed through ensemble averaging of autocorrelation of the off-diagonal components of the pressure tensor using the GK formula [61]: where η, V, k B , and T indicate the dynamic viscosity, the system volume, the Boltzmann constant, and the system temperature, respectively.P αβ is the off-diagonal component of the pressure tensor, which is

Table 1
Simulation steps together with simulation time.computed by the virial equation as follows: where r iα and v iα represent the α-components of the coordinate and velocity of the i th particle, respectively.f iα is the α-component of the force that acts on the i th particle, and N is the number of particles.The final configuration of the system after the equilibrium stage was used to initiate the viscosity calculations with the GK method.The outputs of viscosity for pure glycerol are given in Section 3.2.

SLLOD algorithm
The SLLOD equations of motion, primarily proposed by Hoover and Ladd [62], were later proven to be applicable to investigate fluids undergoing planar Couette flow [63].The mechanical response of the system to an external perturbation is exploited to calculate the viscosity, mainly based on the simplified form of the Navier-Stokes equations.The modified SLLOD equations of motion proposed by Tuckerman et al. are generally written in the form of two differential equations as follows [64]: where m i , r i , and p i denote the mass, position, and momentum of particle I, respectively, F φ i is the total interaction force on particle i, and ∇ u represents the velocity gradient tensor.The integration procedure in LAMMPS updates positions and velocities in each timestep based on the SLLOD equations.In this method, viscosity is calculated based on different shear strain rates (γ) through the following equation: where P ij represents the off-diagonal component of the stress tensor, index "i" shows the direction of flow caused by the imposed shear, and index "j" indicates the normal direction to the flow.
To initiate viscosity calculations using the SLLOD method, a segment comprising approximately one-third of the size of the system final configuration after the equilibrium stage was truncated and selected for further analysis in subsequent stages.In this regard, the number of atoms in the system was reduced to less than a third in order to prolong the simulation time span and reach convergence in the viscosity curve (Fig. 1, Panel a).At the beginning of the simulation, as revealed in Fig. 1. b, a pre-equilibrium stage, including a series of successive simulations with different time steps under the NPT ensemble was carried out to prepare the system for the next stage.Prior to the production run, the system undergoes an equilibrium stage under the NPT ensemble for 30 ns.In addition, to apply shear strain to the simulation box, the box type is changed from orthogonal to triclinic (Fig. 1, Panel a).The viscosity of glycerol in different shear strains, γ = 5.92 × 10 5 , 10 6 , 10 7 , 10 8 , 10 9 s − 1 are calculated and displayed in Section 3.4.Fig. 1 demonstrates the preparation steps of the initial system together with the successive steps for carrying out the MD simulation using the SLLOD algorithm.

Results and discussion
In this section, the results associated with evaluating the viscosity using the computational methods in the previous section are presented.The results include the outputs of the equilibrium simulation, viscosity calculation by the GK and SLLOD methods, as well as the viscosity of pure glycerol in terms of strain rate.

Equilibrium simulation
The root-mean-square deviation (RMSD) of glycerol was computed and depicted in the supplementary information, Figure S1, to examine the structural stability of glycerol molecules during equilibrium simulation.We have found that the glycerol molecules show overall structural stability under the specified thermodynamic conditions during the equilibrium stage.
After carrying out the equilibrium simulation, the MSD parameter for glycerol was calculated to ensure the reliability of the outputs.As shown in Eq. 7, MSD indicates the deviation of the positions of particles with respect to their reference positions over time.
r → (0) and r → (t) indicate the reference and current positions of a particle, respectively, and E demonstrates the dimensionality of the system.According to Eq. 7, in sufficiently high lag time, the diffusion coefficient of the desired particle, D(τ), can be obtained by linear curve fitting of the MSD diagram.Fig. 2 displays the MSD diagram of pure glycerol in terms of lag time together with the linear fitting to MSD at higher lag times.In addition, the MSD parameter for glycerol aqueous solutions with different concentrations is illustrated in Figure S2.
According to Eq. 7 and using a linear fitting on the MSD diagram at high ranges of lag time, the diffusion coefficient for pure glycerol can be evaluated.The diffusion coefficient of pure glycerol was calculated to be 6.037 ± 0.2607 Å 2 /ns, which is in agreement with the results previously reported by Fisher et al. [51].Similarly, the MSD parameter for glycerol in aqueous solutions of 10%, 30%, 50%, 70%, and 80% mol. is presented in the supplementary information (Figure S2).As the concentration of glycerol in aqueous solutions increases, we have observed a rise in the gradient of the MSD diagram, particularly at longer lag times.Therefore, it can be deduced that the diffusion coefficient of glycerol increases in aqueous solutions by adding glycerol to the solution.The outputs of RMSD and MSD in the equilibrium stage demonstrate that the molecular configuration of the system is properly prepared for the following simulation stages.

Viscosity: Green-Kubo formalism
In this section, the time evolution of the viscosity of pure glycerol based on the GK method is presented.The assessments for viscosity were made for 40 ns MD simulation using the method described in Section 2.2.Fig. 3 reveals the time variation of viscosity for pure glycerol.
Each curve in Fig. 3 shows the time evolution of viscosity at a specified pressure maintained at a constant level by the Nosé-Hoover barostat [52,53].The rate of change in the viscosity with pressure is determined from the relative disparity between the different curves in Fig. 3.For different ranges of pressure, the time variation of the viscosity curves converges towards a constant value after 25 ns.Nevertheless, there is still a considerable discrepancy between the calculated values of the viscosity by the Green-Kubo method in Fig. 3 and the experimental values previously reported by Herbst et al. [65] (See Section 3.3).The possible reasons for this difference will be discussed in the following sections.

Viscosity: SLLOD algorithm, pressure-viscosity coefficient
In this stage, the time evolution of the viscosity was calculated by applying the SLLOD algorithm.The viscosity values in this section were obtained based on the method described in Section 2.3.Fig. 4 illustrates the viscosity of pure glycerol at different pressures by implementing the SLLOD algorithm in the molecular dynamics simulation.
As can be interpreted from Fig. 4, the viscosity values fluctuate around a constant value after 100 ns.In order to calculate the viscosity values, the shear strain rate and temperature in the SLLOD algorithm were adjusted to 5.92 × 10 5 s − 1 and 300 K, respectively.It can be inferred from Fig. 4 that the higher the pressure the smoother the viscosity curve.The initial oscillations of the viscosity disappear as the pressure increases.Moreover, the higher the pressure exerted the higher 0.5  the growth rate of the viscosity, as expected from Eq. 6.The calculated viscosity values using the SLLOD algorithm are higher than the viscosity acquired by the GK method; however, there is still a remarkable difference between the viscosity of pure glycerol at high pressures compared to the ones reported experimentally [65].Fig. 5 displays the semi-logarithmic graphs of glycerol viscosity in terms of pressure using the reported experimental values [65], the acquired viscosity by the Green-Kubo method, and the outputs of the SLLOD algorithm.
From Fig. 5, we conclude that at lower pressures (P < 0.1 GPa), there is less disparity between the computed viscosity by the GK and SLLOD methods and the experimental values, while at increased pressure (P > 0.1 GPa), an apparent divergence between the computational and experimental values is observed.In all three cases, the semi-logarithmic viscosity-pressure diagram follows a linear trend, suggesting an exponential expression to describe the relationship between the viscosity and the pressure (η = η 0 e αP ).In Table 2, the evaluated viscosity-pressure coefficient by the experimental method is compared with the outputs of MD simulation for the coefficient.
According to Table 2, the exponential evolution rate of viscosity with pressure from experiment is significantly higher than the corresponding values acquired by the computational methods.We deduce from the viscosity-pressure diagram (Fig. 5) that the SLLOD algorithm gives better agreement with the experimental coefficient than the GK method.Accordingly, the coefficient obtained by the GK method is recognizably different from the experimental value.This difference occurs due to the equilibrating nature of the GK formalism in which the effects of the shear strain rate on the viscosity are practically ignored, while in the SLLOD algorithm, the applied strain rate to the simulation box can be adjusted.It is noteworthy that reducing the shear rate in MD simulation results in a significantly longer simulation time for the viscosity outputs to be converged, which is not computationally affordable.Accordingly, the viscosity of glycerol in higher ranges of strain rates is calculated by the SLLOD algorithm in the next step.

Shear thinning
This section presents the viscosity of glycerol in terms of strain rate at different pressure levels.Variations of the viscosity versus the shear rate at several pressure levels were determined using the SLLOD algorithm, and the results are displayed in Fig. 6.A similar assessment was previously made for squalane by Robbins et al. [38].
As displayed in Fig. 6, at a specific pressure level, the viscosity of glycerol decays by increasing the shear strain rate.Such a rheological response indicates that shear thinning occurs for pure glycerol over a relatively wide range of shear strains.Experimental values of viscosity in lower ranges of shear strain (γ < 10 3 s − 1 ) are indicated based on the reported experimental evaluations of Herbst et al. [65], while in higher ranges (γ > 10 5 s − 1 ), the viscosity was determined by the application of the SLLOD algorithm in NEMD simulations.As can be concluded from Fig. 6, the effect of pressure on viscosity becomes more evident as the shear strain rate is reduced.Besides, the shear thinning procedure in glycerol aqueous solutions was formerly identified in the shear strain range of [20000-80000] s − 1 [66].By combining the low strain rate outputs for viscosity from experiment and high strain rate by computational evaluations from MD simulations, a reliable platform to determine the viscosity over a wide range of strain rates and pressures can be achieved.

Concluding remarks
Viscosity, as the internal resistance of lubricants to flow, plays a crucial role in specifying the tribofilm thickness at contact spots.In the first part of this study, by calculating the MSD parameter and the diffusion coefficient of pure glycerol and comparing them with the reported experimental values, the capability of CHARMM36 force field parameters in calculating the physical properties of glycerol was validated.In the next part, the viscosity of glycerol was calculated using Green-Kubo and SLLOD methods.At pressures<0.5 GPa, it seems that the viscosity values with the Green-Kubo and SLLOD methods are in accordance with the experimental values.While at higher pressures, the

Table 2
The comparison between pressure-viscosity coefficient of glycerol by experimental and MD simulation.

Experimental Data
Fig. 6.Viscosity (η) of glycerol as a function of strain rate (γ) at different pressure levels.
V. Fadaei Naeini et al. viscosity outputs acquired by the Green-Kubo approach reveal a considerable deviation.This deviation occurs due to the equilibrium nature of the Green-Kubo approach, which practically makes it impossible to consider the effects of strain rate on viscosity.At the same time, the SLLOD method, considering the effect of strain rate on viscosity, evaluates the viscosity of glycerol at high pressures significantly closer to the experimental values.Therefore, exploiting the SLLOD algorithm leads to computing the viscosity-pressure coefficient notably closer to the experimental values than the Green-Kubo method.
It is worth mentioning that calculating the viscosity with the SLLOD method at shear rates close to experimental values practically confronts computational challenges since the longer simulation time is required to overcome the thermal noises of the system and reach convergence.Hence, by consolidating the viscosity outputs from experimental procedures and non-equilibrium MD simulation, we can detect the influence of shear thinning on the evolution of glycerol viscosity in a wide range of strain rates.It is crucial to remember that glycerol is highly hygroscopic, which means it can not be practically used as a lubricant in its pure form.Additionally, the viscosity of pure glycerol is too high for most lubrication needs.Therefore, opting for an aqueous solution with optimal concentration is a more suitable choice in such cases.However, the central emphasis of this study primarily revolves around utilizing NEMD simulation as a dependable computational technique for calculating lubricant viscosity.The computational platform introduced in this study facilitates the determination of viscosity for a typical lubricant under extreme pressure and shear rate conditions, whose experimental implementation may be challenging.

Fig. 1 .
Fig. 1.The MD simulation using SLLOD algorithm: a.The preparation steps of the initial configuration of the system, b.Successive steps for carrying out the MD simulation.

Fig. 2 .
Fig.2.Mean squared displacement of pure glycerol with respect to the lag time together with the linear fitting in higher lag times.

Fig. 3 .Fig. 4 .
Fig. 3.The time evolution of viscosity for pure glycerol obtained by the Green-Kubo method.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5.A comparison between the outputs for viscosity of glycerol in terms of pressure by experimental procedure[65], Green-Kubo method, and SLLOD algorithm.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Method η-P Expression Viscosity-Press.Coeff.(α) [GPa − 1 ]Experimental by Herbst et al.[65] η exp.= 1314e5