Determination of diesel physical properties at injection pressures and temperatures via all-atom molecular simulations.
The rheological properties of hydrocarbons are of great interest for the petroleum and automobile industries. Optimization of such properties in fuels and lubricants through a judicious combination of compounds and additives has been a long-standing matter of research. One defining characteristic of such fluids is their shear viscosity. A uniform shear viscosity profile across the range of temperature and pressure is desired for a good lubricant and specific patterns of material behavior would be desired for fuels. Although the experimental determination of this quantity for pure fluids and for mixtures has been pursued, its prediction from inter-molecular interactions is sought. Experimental measurements at elevated temperatures and pressures can also incorporate potential uncertainties related to experimental setup and instrumentation. Under these conditions calculations via molecular simulations takes full account of interactions at an atomic level to calculate the properties of interest without involving the cost of samples and the hazards of extreme conditions under experiments. Thus these models enable us to study such properties under extreme conditions that are either not possible or hazardous for experimental investigation. These molecular methods will provide the ability to fine tune the composition for a product optimization as well as to potentially understand the influence of novel additives to fuels and lubricants. Molecular simulation has long been used to gain understanding of phenomena observed in physical systems. In particular, macroscopic transport properties such as viscosity, diffusivity, thermal conductivity, etc., may be computed by using molecular simulation methods such as the molecular dynamics or Monte Carlo methods. There are two types of molecular dynamics approaches to these calculations: equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD). In the EMD approach, transport properties are obtained based on the unforced response of a collection of energy conserving molecules by using the so-called Green-Kubo relations. The basis of such approximation is the observation that the molecular response due to a weak external excitation decays in the same way as a spontaneous fluctuation at equilibrium.
Prediction of alkane viscosity at high pressures up to 2000 bar (typical in a common rail diesel engine condition) is a rather challenging problem. A critical review of the literature reveals that majority of the molecular simulation studies of alkane viscosity use various United Atom (UA) models which typically underestimate viscosity by 20% to 80% . There are relatively fewer reports on simulation of transport properties using all-atom (AA) potentials  . Allen and Rowley optimized the OPLS-AA force field interaction parameters for both short and long alkanes . They observed an increasing positive deviation for viscosity of n-butane with increasing density. For n-nonane, Martin et. al. observed 40% deviation from experimental viscosity using the OPLS-AA force field . In case of perfluoroalkanes, the OPLS-AA force field interaction parameters shows 15% to 35% deviation from experimental viscosity values .
Recently we have developed an EMD based approach using multiple independent trajectories with all-atom interaction potential to reproduce n-hexadecane (representative of GTL diesel fuel) and ndecylbenzene (representative of conventional diesel fuel) viscosity (up to 750 bar pressure and at 423 K temperature) within 12% and 8% accuracy respectively (unpublished data). For lubricant size molecules we found NEMD to be more effective in reproducing viscosity temperature relation (within 6% accuracy for C30 isomers) at ambient pressure as well as magnitude of viscosity even at high pressure up to 2000 bar (within 10% accuracy for C30 isomers) 
In this article we have extended our EMD protocol to study the viscosity profile of two binary mixtures as a function of pressure. One mixture involves two n-alkanes, n-octane ([C.sub.8][H.sub.18];Figure 1a) and n-dodecane ([C.sub.12][H.sub.26];Figure 1b) (case 1) and the other involves two diesel surrogates, n-hexadecane and n-decylbenzene (case 2). n-hexadecane ([C.sub.16][H.sub.34]) represents GTL (Gas to Liquid) like diesel and while n-decylbenzene ([C.sub.16][H.sub.26]) represents regular aromatic diesel as shown in Figure 1c and 1d respectively.
In case 1, we compare our simulation results with experimental data available from literature. In case 2, we compute the viscosities of various binary mixtures (20:80, 40:60, 60:40 and 80:20) of the two representative diesel surrogates (n-hexadecane and n-decylbenzene).
In our previous work with n-hexadecane, we have used the TuTobias-Klein force field  with which we were able to reproduce densities and viscosities within 1% and 10% relative errors vs. experimental values respectively . The Tu-Tobias-Klein force field has been particularly developed for linear molecules and parameters are not available for other topologies e.g. cyclic aromatic/aliphatic molecules. For this reason, we have chosen to adopt a more generic force-field OPLS (Optimized Parameters for Liquid Simulations)  for our simulations. It is known from literature that OPLS all-atom force field fails to reproduce viscosity results for long hydrocarbons in its native format. Alternative approaches have been reported [9,10] with changes in Lennard-Jones parameters, backbone dihedral (C-C-C-C) parameters and partial charges of C and H atoms to reproduce densities, heats of vaporization, self-diffusion and viscosities of long hydrocarbons. We have tested the parameters reported in some of the external publications to see their ability to accurately reproduce density/viscosity of n-hexadecane. Based on our study we found that the parameters developed by Sui et. al  (LOPLS) reproduced densities and the heats of vaporization of alkanes and alkenes of different lengths (both short and long hydrocarbons) and therefore were chosen for n-hexadecane. As the study was limited to alkanes and alkenes, no specific parameters or comments were reported for aromatic components. In case of n-decylbenzene, parameters for aromatic/ring carbons and the respective bonded hydrogen atoms are not reported by the Sui et. al., we have used OPLS parameters. Further in this report, we continue to refer to this force field as LOPLS.
Our protocol includes two steps. The first step is to estimate the density of the system. For this we run isotropic NPT simulation (where the number of particles, pressure and temperature remain fixed) for ~4ns at the concerned temperature and pressure. The damping coefficient for the temperature thermostat and pressure barostat was 1ps. Average convergence density is identified using the running average of density and the box density is then fixed to the average density by adjusting the box volume. Using this box as a starting point, several independent NVT simulations (where the Number of particles, Volume and Temperature remain fixed) are started in a protocol of heat [right arrow] cool [right arrow] production run. In order to have good statistics typically 20 such simulations are run. The heat/cool cycle ensures a different equilibrated system for each independent production run. Typically the heat and cool cycles are of 1ns each while the production run is ~8-16ns. Pressure tensor components are recorded during the production run. We used LAMMPS  package for all our simulations with default thermostat (Nose-Hoover chain) and barostat (implemented by Sinhoda et. al. ).
The Green-Kubo formalism has been used to estimate the viscosity from EMD simulations. The EMD approach considers a collection of energy conserving particles. The equation of motion is given by the Newton's Second Law:
m[??] + [partial derivative]P/[partial derivative]r = 0 (1)
where P is the potential energy. Transport properties such as viscosity, elasticity, diffusion etc., may be approximately obtained through the EMD simulation by using the Green-Kubo relations. The basis of such approximation is the observation that the response due to a weak external excitation decays in the same way as a spontaneous fluctuation at equilibrium. Specifically for the shear viscosity, the Green-Kubo relation is given by:
[[eta].sup.xy] = [v/[k.sub.B]T][[integral].sup.[infinity].sub.0]<[[sigma].sup.xy](0)[[sigma].sup.xy](t)>dt (2)
where <.> donates the ensemble average but computed using time average (based on the ergodicity assumption), V is the volume, [k.sub.B] the Boltzmann constant, T the temperature, and [[sigma].sup.xy] the xy component of the symmetric traceless stress tensor.
RESULTS AND DISCUSSION
Case 1: n-decane, n-dodecane and Their 50:50 Binary Mixture
Densities and viscosities of n-decane, n-dodecane and their 50:50 binary mixture have been reported by Dymond et. al.  up to 100C, and pressures up to 5000 bar using a modified falling body viscometer. The experimental viscosities are shown in Figure 2 (the line is to guide the eye). As expected using mixing rules (see Appendix), the binary mixture viscosity falls roughly between the individual pure components. In this study we predict the densities and viscosities of the pure components and their binary mixture at 100C and up to 2000 bar.
For simulations, we studied the systems using 100 molecules for pure components, i.e., n-octane and n-dodecane and in case of 50:50 mixture, we had 50 molecules of each n-octane and n-dodecane. The results are shown in Tables 1, 2, 3, Figures 2 and Figure 3.
As seen from the tables (Tables 1, 2, 3) and Figure 3, we were able to predict the densities within 2% relative to the experimental values using LOPLS interaction potential. Also the densities follow a linear addition rule (equation 3, where [x.sub.A] and [x.sub.B] are mole fractions of respective two components with densities [[rho].sub.A] and [[rho].sub.B] respectively) that can be used to predict densities of binary mixtures. In case of viscosities, pure component viscosities were predicted within 15% while the binary mixture viscosities were predicted within 25% for temperatures up to 100C and pressures up to 2000 Bar (see Tables 1, 2, 3 and Figure 4).
[[rho].sub.m] = [x.sub.A][[rho].sub.A] + [x.sub.B][[rho].sub.B] (3)
There are many empirical correlations reported in literature quantifying pressure effects on viscosity. One of the most widely used correlations is given by the Barus  which is
[[eta].sub.p] = [[eta].sub.0][e.sup.[alpha]p] (4)
where, [[eta].sub.p] is the viscosity at pressure p, [[eta].sub.0] is the atmospheric viscosity, [alpha] is the coefficient andp is the pressure of concern. Barus relation is valid for pressure below 5000 bar. Thus for diesel fuel case we can rely on it. However, for applications in lubrication sciences where pressure can go beyond 10000 bar (1 GPa) alternative correlations are reported.
In literature, several empirical relations exist for predicting binary mixture viscosities from individual component viscosity and other physical properties. A simple one is Grunberg-Nissan model , which is widely used for predicting binary mixture viscosities. The Grunberg-Nissan model estimates the binary mixture viscosity by
ln [eta] = [x.sub.1] ln [[eta].sub.1] + [x.sub.2] ln [[eta].sub.2] + [x.sub.1][x.sub.2][G.sub.12] (5)
where, x is the mole fraction, subscripts 1 and 2 denote the mixture components. [G.sub.12] is the interaction parameter which is characteristic of the mixture composition. [G.sub.12] is a positive number and typically would lie between 0 and 2 for alkanes. For alkane binary mixtures [G.sub.12] increases with the difference in carbon number of the pure components. The values tend to increase with an increase in pressure and the rate of increase becomes larger with the difference in carbon number. For a given system, [G.sub.12] tends to decrease with increasing temperature, but the effect is small and becomes even less significant as the carbon numbers of the two pure components become closer. Using viscosity prediction models, binary mixtures viscosities of alkanes are generally predicted between the respective pure component viscosities (or in some cases, based on the value of interaction parameter, the binary mixture viscosity can be slightly higher than the respective pure component viscosities). Our predictions as well as the experimental data follow the expected behavior of Arrhenius relationships (see Appendix), or single parameter dependent relationships such as Grunberg-Nissan for viscosity predictions mentioned above and viscosities can be predicted after fixing the interaction parameter (positive number) for the molecules under study.
Case 2: n-hexadecane, n-decylbenzene and Their Binary Mixtures
We performed EMD simulations for all four mole fractions, i.e., 20:80, 40:60, 60:40 and 80:20 using the modified LOPLS force field interaction parameters as reported in Section 2. We performed NPT simulations fixing the desired temperature and pressure and calculated the equilibrated densities as the system evolved and adjusted to the conditions. These equilibrated boxes at the correct densities were then used as the starting point for several NVT simulations as describe in the methods section to record the pressure tensor components. Viscosities were calculated using the Green-Kubo formalism. As previously noted that LOPLS can predict viscosities within 20% relative errors for pressures only up to 750 bar as we observed that the force field interaction parameters were inaccurate for higher pressures. Hence we calculated binary mixture viscosities only up to 750 bar in this case. The number of molecules in the simulations was 100 for all 4 binary mixtures with varying composition of n-hexadecane and n-decylbenzene and the respective and pure components. This provides with a system size which is sufficient for accurate prediction and does not improve by increasing the number of molecules. We tested system size effects for pure hexadecane with 256 molecules at 423K and 1 bar without any significant change in predicted densities (0.7459 gm/cc for 200 mols as compared to 0.7497 gm/cc for 100 mols). Throughout this study a mole fraction of 0 represents pure n-hexadecane and a mole fraction of 1 represents pure n-decylbenzene.
The calculated densities are shown in Figure 5 and Table 4 for all 4 binary mixtures at 150C and pressures of 1,250, 500 and 750 bars. Binary mixture densities can generally be estimated by simple addition mixing rules as expressed by Equation 3 and it can be noted that all calculated binary mixture densities follow such uniform addition mixing rule.
We also calculated the viscosities of the all four mixtures at all four pressure points at 150C using EMD simulations and using the Green-Kubo analysis on the recorded pressure tensor components. The results are summarized in Table 5 and Figure 6.
As seen from Figure 6 the predicted viscosities from EMD simulations follow the Barus relationship (increasing viscosity with pressure) with a (viscosity-pressure coefficient) ranging from 0.001 to 0.002. Further the predicted binary mixture viscosities lie within their respective pure component viscosities at a given pressure (expectedly, viscosity increases as we increase the mole fraction of ndecylbenzene to substitute n-hexadecane) which can be explained by using empirical mixing rules for viscosity such as Grunberg-Nissan relation.
Fuel injection systems work at high pressure and temperatures. The individual components in the fuel influence the behavior of fuel at these conditions and in turn control the injection behavior correlating with the power achieved from the engine. High temperature (up to 150C) and high pressure (up to 2000 bar) are the operating conditions of a CRDi (Common Rail Diesel injection) engine. Experimental measurements are also subject to errors at these conditions.
Molecular modeling is a valuable tool in this aspect as we have previously shown  that molecular modeling can reliably predict the properties of interest (density and viscosity) as a function of pressure and temperature including high temperature and high pressure conditions relevant to CRDi engine operating conditions. In this report we studied the binary mixture viscosities as a function of mole fraction and pressure at high temperature (150C). For modeling the pure components and mixtures, we have identified transferable interaction parameters for the physical property prediction (density and viscosity) for n-alkanes and n-alkyl benzene molecules. Corrected OPLS force field interaction parameters for long hydrocarbons  in conjunction with the OPLS parameters for ring atoms was used for all molecules studied in this report. We also extended our viscosity prediction protocol from single component to two component binary mixtures using EMD simulations. Densities and viscosities for pure components as well as their binary mixtures are predicted within ~4% and ~25% relative error respectively. In the case of binary mixtures of n-octane and n-dodecane, we could compare our viscosity predictions to the published experimental data by Dymond et. al.  and they are within 25% relative error to the respective measurements.
We conclude that all atom molecular dynamics simulations allow a transferable, reliable and predictive methods to calculate densities and viscosities for pure components and their binary mixtures at a wide range of pressures and temperatures including realistic and operating elevated conditions of fuel injection systems and combustion engines.
Abhinav Verma, Roger Cracknell, David Doyle, and Indranil Rudra
[1.] Dysthe, D.K., Fuchs A. H. and Rousseau, B. "Fluid transport properties by equilibrium molecular dynamics. III. Evaluation of united atom interaction potential models for pure alkanes, " J. Chem. Phys. 112: 7581-7590, 2000
[2.] Allen, W. and Rowley, R. L., "Predicting the viscosity of alkanes using nonequilibrium molecular dynamics: Evaluation of intermolecular potential models, " J. Chem. Phys., 106: 10273, 1997
[3.] Martin, M. G. and Thompson, A. P., "Industrial property prediction using Towhee and LAMMPS," Fluid Phase Equilib., 217:105-110, 2004
[4.] Morgado, P, Laginhas, C. M. C., Lewis, J. B., McCabe, C., Martins, L. F. G. and Filipe, E. J. M., "Viscosity of liquid perfluoroalkanes and perfluoroalkylalkane surfactants," J. Phys. Chem. B, 115:9130, 2011
[5.] Shell Oil Company, "Molecular Model to Understand the Effect of Pressure on Base Oil Viscosity," Research Disclosure, 618059, 2015
[6.] Tobias, D. J., Tu, K and Klein, M. L., "Assessment of all atom potential for modelling membranes: Molecular dynamics simulations of solid and liquid alkanes and crystals of phospholipid fragments," J. Chim. Phys., 94:1482-1502, 1997
[7.] Payal, R. S., Balasubramanian, S., Rudra, I., Tandon, K., Mahlke, I., Doyle, D. and Cracknell, R., "Shear viscosity of linear alkanes through molecular simulations: quantitative tests for n-decane and n-hexadecane,"Mol. Sim., 38:1234-1241, 2012
[8.] Kaminski, G. A., Friesner, R. A., Tirado-Rives, J. and Jorgensen, W. L., "Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides," J. Phys. Chem. B, 105:6474-6487, 2001
[9.] Siu, S. W. I., Pluchackova, K. and Bockmann, R. A., "Optimization of the OPLS-AA Force Field for Long Hydrocarbons, " J. Chem. Theory Comput, 8:1459-1470, 2012
[10.] Murzyn, K., Bratek, M. and Pasenkiewicz-Gierula, M. "Refined OPLS All-Atom Force Field Parameters for n-Pentadecane, Methyl Acetate, and Dimethyl Phosphate," J. Phys. Chem. B, 117:16388-16396, 2013
[11.] Plimpton, S., "Fast Parallel Algorithms for Short-Range Molecular Dynamics," J. Comp. Phys., 117:1-19, 1995
[12.] Shinoda, W., Shiga, W. and Mikami M., "Rapid estimation of elastic constants by molecular dynamics simulation under constant stress," Phys. Rev. B., 69:134103, 2004
[13.] Dymond, J. H., Robertson, J. and Isdale, J. D., "Transport Properties of Nonelectrolyte Liquid Mixtures-III. Viscosity Coefficients for n-Octane, n-Dodecane, and Equimolar Mixtures of n-Octane + n-Dodecane and n-Hexane + n-Dodecane from 25 to 100 C at Pressures Up to the Freezing Pressure or 500 MPa, " Int. J. Thermophys., 2:133, 1981
[14.] Barus, C., "Isothermals, isopiestics and isometrics relative to viscosity, " Am. J. Sci, 45:87-96, 1893
[15.] Marczac, W., Adamczyk, N. and Lezniak, M. "Viscosity of Associated Mixtures Approximated by the Grunberg-Nissan Model," Int. J. Thermophys., 33:680-69, 2012
FIE--Fuel Injection Equipment
GTL--Gas to Liquids
OPLS--Optimized Parameters for Liquid Simulations
NPT--No of Particles, Pressure, Temperature
NVT--No of Particles, Volume, Temperature
CRDi--Common Rail Diesel injection
ARRHENIUS MIXING RULES FOR VISCOSITIES
Various mixing rules are available in the literature for predicting the viscosity of binary mixtures. Several of these rules are Arrhenius in nature while some other rules take interaction of the two components into consideration. We looked into the several such rules and plotted using the viscosities of two pure components of our interest (n-hexadecane at 3.5 mPas and n-decylbenzene at 12.30 mPas).
Pure Mixing Rules
We studied the following pure mixing rules and are shown in Figure 7.
1. Arrhenius equation: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]
2. Bingham equation: 1/[mu] = [x.sub.A]/[[mu].sub.A] + [x.sub.B]/[[mu].sub.B]
3. Kendall and Monroe: [[mu].sup.1/3] = [x.sub.A][[mu].sup.1/3.sub.A] + [x.sub.B][[mu].sup.1/3.sub.B]
4. Linear: [mu] = [x.sub.A][[mu].sub.B] + [x.sub.B][[mu].sub.B]
Mixing Rules with Additional Interaction Parameters
We also studied two mixing rules with interaction parameters.
1. Power law: [mu] = [([x.sub.A][[mu].sup.n.sub.A] + [x.sub.B][[mu].sup.n.sub.B]).sup.1/n], where n is the interaction parameter.
2. Grunberg-Nissan Model: ln [mu] = [x.sub.A] ln [[mu].sub.A] + [x.sub.B] ln [[mu].sub.B] + [x.sub.A][x.sub.B][G.sub.AB]
It may be noted that for all mixing rules, the mixture viscosities lie within the neighborhood of the pure component viscosities. There are cases when the mixture viscosities are out of the range defined by the pure components, but these are still not drastically different from the respective pure component viscosities. Arrhenius types of mixing rules are parameter independent, but they force the viscosities to lie within the pure component viscosities. In case of parameter dependent mixing rules, the extra parameter then allows for the deviations to incorporate mixture viscosities that may lie slightly outside the pure component viscosities. However the parameter is not transferable and is dependent on the two components used.
Table 1. Experimental  vs Simulation results (density and viscosity) for n-octane T (C) Press Density (gm/cc) % error Viscosity % error (Bar) (mPas) SIM EXPT SIM EXPT 25 1 0.6991 0.6985 0.1% 0.60 0.51 18% 100 1 0.6248 0.6355 -1.7% 0.29 0.25 16% 100 1036 0.7146 0.7215 -1.0% 0.63 0.57 10% 100 1977 0.7551 0.7594 -0.6% 1.07 0.95 13% Table 2. Experimental  vs Simulation results (density and viscosity) for n-dodecane T (C) Press Density (gm/cc) % error Viscosity % error (Bar) (mPas) SIM EXPT SIM EXPT 25 1 0.7494 0.7453 0.6% 1.46 1.36 8% 100 1 0.6833 0.6904 -1.0% 0.48 0.51 -6% 100 241 0.7067 0.7140 -1.0% 0.71 0.67 6% 100 959 0.7535 0.7597 -0.8% 1.37 1.23 12% Table 3. Experimental  vs Simulation results (density and viscosity) for 50:50 binary mixture of n-octane and n-dodecane. T (C) Press Density (gm/cc) % error Viscosity % error (Bar) (mPas) SIM EXPT SIM EXPT 25 1 0.7284 0.7259 0.3% 0.95 0.87 9% 100 1 0.6593 0.6680 -1.3% 0.42 0.37 14% 100 966 0.7370 0.7410 -0.5% 1.10 0.89 24% 100 2026 0.7798 0.7809 -0.1% 2.01 1.61 25% Table 4. Predicted densities for the different mole fractions at T=150C and different pressures (0 and 1 represent pure n-hexadecane and pure ndecylbenzene). Pressure (Bar) Density (gm/cc) 0.8 0.6 0.4 0.2 1 0.7552 0.731466 0.711181 0.691435 250 0.7793 0.758478 0.736636 0.718854 500 0.7986 0.777401 0.75742 0.738499 750 0.8147 0.793295 0.773305 0.753894 Table 5. Predicted viscosities for the different mole fractions at T=150C and different pressures (0 and 1 represent pure n-hexadecane and pure ndecylbenzene). Pressure (Bar) Viscosity (mPas) 1 (n-dec) 0.8 0.6 0.4 0.2 0(n-hex) 1 0.61 0.63 0.63 0.59 0.56 0.42 250 0.90 0.83 0.80 0.68 0.64 0.50 500 1.30 1.18 1.14 1.01 0.93 0.88 750 1.92 1.75 1.60 1.50 1.25 1.15
|Printer friendly Cite/link Email Feedback|
|Author:||Verma, Abhinav; Cracknell, Roger; Doyle, David; Rudra, Indranil|
|Publication:||SAE International Journal of Fuels and Lubricants|
|Date:||Nov 1, 2016|
|Previous Article:||Formation of intake valve deposits in gasoline direct injection engines.|
|Next Article:||Future specification of automotive LPG fuels for modern turbocharged DI SI engines with today's high pressure fuel pumps.|