Elucidating the Properties of Surrogate Fuel Mixtures Using Molecular

Jan 15, 2016 - Calculated results were compared to experimentally determined values and to values obtained with the nonreactive, all-atom version of t...
0 downloads 10 Views 4MB Size
Download PDF
Article pubs.acs.org/EF

Elucidating the Properties of Surrogate Fuel Mixtures Using Molecular Dynamics Barbara L. Mooney,†,‡ Brian H. Morrow,†,‡ Keith Van Nostrand,∥ Dianne Luning Prak,† Paul C. Trulove,† Robert E. Morris,# J. David Schall,⊥ Judith A. Harrison,*,† and M. Todd Knippenberg¶,‡ †

Department of Chemistry, United States Naval Academy, Annapolis, Maryland 21402, United States Department of Chemistry, High Point University, High Point, North Carolina 27262, United States ⊥ Mechanical Engineering Department, Oakland University, Rochester, Michigan 48309, United States ∥ School of Physics and Astronomy, Rochester Institute of Technology, Rochester, New York 14623, United States # Chemistry Division, U.S. Naval Research Laboratory, Washington, DC 20375, United States ¶

ABSTRACT: The wide compositional differences between conventional and alternative fuels have resulted in much research aimed at determining which alternative fuels can be used, and in what proportions, in conventional engines. Atomic-scale modeling is uniquely positioned to lend insight into this question without extensive large-scale tests. The predictive power such modeling affords could narrow the phase space that must be examined experimentally. This study utilizes molecular dynamics (MD) simulations to predict the properties of a set of pure hydrocarbons, as well as binary and multicomponent surrogate fuel mixtures for alternative fuels created from these pure components. The accuracy and transferability of the modified LennardJones adaptive intermolecular reactive empirical bond-order potential (mod-LJ AIREBO) [Liu, A.; Stuart, S. J. J. Comput. Chem. 2008, 29, 601−611] was assessed by calculating densities, heats of vaporization, and bulk moduli of pure hydrocarbons and the mixtures of these hydrocarbons, i.e., surrogate fuels. Calculated results were compared to experimentally determined values and to values obtained with the nonreactive, all-atom version of the optimized potential for liquid simulations (OPLS-AA) [Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem. Soc. 1996, 118, 11225−11236]. The mod-LJ AIREBO potential quantitatively predicts the densities of the pure hydrocarbons and binary mixtures of n-dodecane and 2,2,4,4,6,8,8heptamethylnonane (isocetane). It is interesting to note, that despite doing an excellent job predicting the densities of the pure hydrocarbons, the performance of the mod-LJ AIREBO potential degrades when predicting the densities of the multicomponent surrogates and mixtures of n-hexadecane and isocetane, implying that it is not straightforward to extend potentials fit with pure compounds to mixtures. The OPLS-AA potential also has difficulty quantitatively predicting the densities of mixtures, although a new parameter set for long-chain hydrocarbons (L-OPLS) [Siu, S. W. I.; Pluhackova, K.; Bockmann, R. A. J. Chem. Theory Comput. 2012, 8, 1459−1470] yields some improvement for binary mixtures. Heat of vaporization predictions using both potentials also agree reasonably well with experiment. Bulk moduli predictions using the mod-LJ AIREBO potential are consistently higher than, and do not quantitatively agree with, the experimental values. In contrast, bulk moduli predictions using the OPLS-AA potential are generally in good agreement with experimental values. Despite the success of the OPLS-AA potential predicting the bulk moduli of individual hydrocarbons, it is unable to quantitatively predict the bulk moduli of the multicomponent surrogates. Interestingly, the use of the L-OPLS parameter set improves density predictions but not predicted bulk moduli values.

■

INTRODUCTION Conventional hydrocarbon-based fuels for diesel and jet engines are typically mixtures of straight-chain alkanes, branched alkanes, cyclic alkanes, aromatics, and alkenes, with various impurities. Diesel fuels and their proposed replacements are highly complex and variable mixtures of mostly hydrocarbons,1 and this variability affects a broad range of properties, most of which are critical to performance in engines. In contrast, alternative fuels contain components from a small number of hydrocarbon classes. For example, both Fischer− Tropsch fuels and hydrotreated renewable diesel derived from algae are principally comprised of only straight-chain and branched alkanes.2−4 Because the composition, and thus performance, of alternative fuels can vary widely depending upon the source, much research has focused on developing an understanding of © 2016 American Chemical Society

how composition differences impact properties and performance. One approach to studying fuel composition, property, and performance relationships is to create surrogate fuels with known compositions. Combustion properties and physical properties of the surrogate fuels can then be measured with the hope of determining how each component impacts overall fuel performance. Elucidating molecular-level rules governing the composition-property-performance relationships should aid in the development of fuels tailored to specific needs. Additionally, having the ability to predict properties from composition alone could reduce or even eliminate the need for expensive experimental testing of unsuitable candidate fuels.5 Received: June 29, 2015 Revised: January 13, 2016 Published: January 15, 2016 784

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

Thus, the long-range interactions are dependent on the chemical environment of each carbon atom. The parametrization of the mod-LJ AIREBO potential was accomplished by fitting to data from a collection of low molecularweight hydrocarbons. These hydrocarbons were linear with one to six carbon atoms with benzene being the only cyclic hydrocarbon used in the fitting. To date, the liquid properties of higher molecular-weight, pure hydrocarbons, such as those that would be found in liquid fuels, or mixtures of these hydrocarbons have not been calculated using the mod-LJ AIREBO potential. The published mod-LJ AIREBO intermolecular parameters are used with one modification: the value of σC (4,0), which corresponds to a carbon atom with four carbon neighbors and no hydrogen neighbors, was modified from 2.35 to 2.70 Å. Because this parameter was not fit originally to experimental data, but was chosen by extrapolating from other values,11 this change does not impact quantities reported as part of the original publication. This change was made after initial simulations of highly branched aliphatic molecules, such as isocetane, had densities deviating from experimental values. For comparison purposes, MD simulations were also carried out using a nonreactive force field, known as the optimized potential for liquid simulations (OPLS).22,23 Recent MD simulations have shown that, as the molecular chain length increases, the reliability of united atom (UA) potentials decreases.29 Because fuels are composed of long-chain hydrocarbons, the all atom (AA) version of OPLS, as implemented in the LAMMPS MD code,24 is adopted here with the parameters distributed by Tinker.30 To generate a simulation box, the Cartesian coordinates of a single molecule were generated using chemical drawing software.31 The Packmol open-source molecular-packing tool,32 which randomly places molecules within the simulation box, was then used to replicate this coordinate file many times to fill the box. Packmol uses the GENCAN packing algorithm,33 which allows for rotation of the molecules as well as backtracking if an earlier molecule is placed poorly with respect to distance between neighboring molecules. For simulations involving pure hydrocarbons, 256 molecules were placed in a cubic periodic box, with the length of each side chosen so that the initial box size would be 5 Å larger on each side than expected at the experimental density. For the binary surrogates, 256 molecules were placed in the simulation box with dimensions based upon the target experimental mole fractions for simulations containing two molecule types.20,34 Selected simulations containing 512 and 1024 molecules yield results that are within the error bars of the results using 256 molecules. A larger sample size than that used previously for mod-LJ AIREBO simulations11 was selected here to ensure that adequate numbers of molecules were present for each component in the hydrocarbon surrogate mixtures. Once this initial packing of molecules was complete, NVT simulations, where the number of atoms (N), the volume (V), and the temperature (T) are all constant, for 50 ps were carried out. A thermostat was used to maintain the temperature at 298 K and to allow the positions, bond angles, and torsions of each molecule to randomize relative to the initial placement. The mod-LJ AIREBO is not yet available as part of the LAMMPS distribution; therefore, an inhouse serial MD code was used for all simulations with this potential. A Berendsen thermostat35 has been implemented in that code and was used to maintain the temperature at 298 K. With the number of atoms (N), pressure (P), and temperature (T) remaining constant, simulations were then carried out with a target P of 1 atm until the total box volume reached equilibrium, which is reached when the system pressure oscillates about the target pressure. The Berendsen barostat uses a scaling factor, μ, which is given by κ Δt μ = 1 − τ (Psystem − Ptarget). The instantaneous pressure of the

While various methods exist to predict the physical properties of hydrocarbons, they do not provide molecular insight into mixing, degradation, or other dynamic behavior. Many properties of fuels can be calculated using all-atom molecular dynamics (MD) simulations.6−8 In addition to thermodynamic properties,9 MD yields atomic-scale insight into mixing, other dynamic behavior, and temperature-dependent behavior and can yield insights into fuel degradation. In MD, a potential energy function (PEF) is used to calculate the bonded (intramolecular) and nonbonded (intermolecular) interactions between atoms and is therefore responsible for molecular motion. Because we are interested in fuel degradation and dynamic behavior in addition to thermodynamic properties, we benchmark the performance of a modified version of the adaptive intermolecular reactive empirical bondorder (AIREBO) potential,10 namely, the mod-LJ AIREBO11 in reproducing selected fit-for-purpose (FFP) properties.12 The room-temperature density, enthalpy of vaporization, and bulk modulus at one atm of several pure hydrocarbons found in fuels, as well as two-component and multicomponent fuel surrogates, are calculated. Because many alternative fuels are composed of linear and branched hydrocarbons,13−15 the systems studied here are samples of simple, pure hydrocarbons, binary mixtures of n-dodecane and n-hexadecane with isocetane (2,2,4,4,6,8,8-heptamethylnonane) (i.e., two-component surrogates), and multicomponent surrogates, composed of n-alkanes, branched alkanes, cyclic alkanes, alkenes, and aromatic hydrocarbons. Herein, the performance of the mod-LJ AIREBO in reproducing the selected FFP properties is compared to experimental data, which have been measured recently16−20 or were measured as part of this work, and to quantities calculated using a second, widely available PEF, the all-atom version of the optimized potential for liquid simulations (OPLS-AA).21−23 This nonreactive potential is distributed as part of the open-source MD code LAMMPS.24

■

METHOD

Bond-order potentials were developed originally to model solid-state covalent materials (e.g., diamond).25,26 The hydrocarbon REBO potential27 was extended to include nonbonded interactions to model liquid systems and self-assembled monolayers while maintaining its ability to model chemical reactions. The potential, known as the AIREBO potential,10 uses a Lennard-Jones (LJ) 12-6 equation, with the typical two parameters σ and ε, to model intermolecular interactions between two atoms. This potential has one set of intermolecular parameters for all carbon−carbon interactions, a second set for all carbon−hydrogen interactions, and a third set for all hydrogen−hydrogen interactions. Thus, a molecule could have two carbon atoms with different hybridizations, but both atoms would use the same values of σ and ε. To better account for the long-range interactions based on different hybridizations of carbon, an adaptive treatment of intermolecular interactions was later added to the AIREBO potential.11 In this mod-LJ AIREBO model, different parameters are used for carbon atoms depending upon their local bonding environments. Hydrogen retains single values for σ and ε because it has a single hybridization in covalent compounds. Because this potential retains the ability to model chemical reactions present in the underlying REBO potential, it must be able to account for atoms that are “between” bonded and nonbonded states: that is, during a chemical reaction, an atom may have a fractional number of neighbors. A bicubic-spline function28 is used to interpolate between (integer) values of the number of neighbors that are carbon and hydrogen, NC and NH, respectively.

P

system (Psystem) is calculated at every time step and compared to the target pressure, Ptarget. The simulation time step of 0.25 fs, Δt, is multiplied by the isothermal compressibility, κ, which should be set to allow the system to adequately explore phase space, and divided by a relaxation time of the pressure τP, which is 10 ps. The quantity μ is used to scale each periodic boundary (PB) dimension, as well as the 785

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels Cartesian coordinates of each atom at every time step. The PB lengths will be reduced if pressure needs to be increased, or grow in size if pressure needs to be reduced, scaling the dimensions of the simulation system. Time integration is carried out by solving Newton’s equations of motion using the Verlet algorithm.28 The OPLS-AA potential is distributed as part of LAMMPS,24 which is an open-source, parallel processing program for MD simulations. NPT simulations using this potential within LAMMPS are performed using the Nosé-Hoover method. The Nosé-Hoover algorithm also modifies the dimensions of the simulation box to reach a target pressure. Atomic coordinates are scaled to achieve a target pressure, while performing time integration using the Nosé-Hoover equations of motion to integrate time as well as thermostat to a target temperature.36 Because two different algorithms to adjust the simulation pressure were used, the following test was performed: Two identical MD simulations of benzene using the AIREBO potential in LAMMPS were carried out, one with the Nosé-Hoover algorithm and the second using the Berendsen algorithms for pressure control. After equilibration, both simulations produce the same average volume and bulk modulus within one standard deviation. This indicates that a choice of barostats is not critical to the results reported here as long as the simulation length is sufficient to reach an average volume. The simulations typically reached equilibrium during NPT dynamics during the first 250 ps of simulation time. This was confirmed by calculating the velocity autocorrelation functions of several systems over the last 75 ps of this equilibration phase (Figure 1). It is clear from analysis of

property calculation. Results presented here were calculated from the final 200 ps of this segment of the simulation. For all systems studied, five independent starting configurations were examined with the reported data being averages of these five separate simulations. As a check, longer simulations of 2 ns were also run on selected systems and the results are within the error bars reported here. To reduce the computational time associated with the multicomponent surrogate systems, the simulation-box size for each surrogate was determined by calculating an estimated volume of each single molecule contained in the surrogate and then calculating the total number of molecules required to fill a cubic box of approximately 40 Å on each side. This distance corresponds to approximately 4.5 times the interaction range of the potential, larger than the “rule of thumb” distance of 3 times the interaction range. It also allows for an increased number of molecules of each type, particularly when the mole fraction of a minor component is small (Table 1). Once the total number of molecules was calculated, they were placed in a starting box of 45 Å on each side and equilibrated as described above. The enthalpy of vaporization ΔHvap is calculated at 298 K using the equation ΔHvap = Egas − (1/N)Eliq + RT, where Eliq is the potential energy of simulation box of the pure liquid hydrocarbon at temperature T, R is the gas constant, N is the number of molecules in the simulation, and Egas is the potential energy of an isolated hydrocarbon molecule in the gas phase. (See Wang and Hou for a derivation of this equation.37) Bulk modulus (B) is a measure of compressibility that is an important determinant of fuel suitability in engines, particularly for diesel engines. Within the context of MD simulations, there are several ways to calculate B.38−40 The equation B =

(

kBT ⟨V ⟩ ⟨V 2⟩ − ⟨V ⟩2

) , where B is P

bulk modulus, kB is the Boltzmann constant, T is the target temperature of the thermostat, ⟨V⟩ is the simulation average volume, and ⟨V2⟩ − ⟨V⟩2 is the variance of the volume, was used here. The subscript P denotes calculations using NPT dynamics. This formula was chosen over alternative formulations based on pressure or total energy due to the variations in system pressure for simulation systems of the size used here and the relative insensitivity in total energy of liquid samples to small changes to volume. For systems where literature values of density or bulk modulus are not available, density and speed of sound were measured using methods described previously.16,18 The isentropic bulk modulus (Ev) in Pa can then be calculated using Ev = c2ρ, where c and ρ are the speed of sound and the density in SI units, respectively.18 Briefly, the instruments used were an Anton Paar DSA 5000 density and sound analyzer and an Anton Paar SVM 3000 Stabinger viscometer. Each

Figure 1. Velocity autocorrelation function of the molecular center of mass calculated over the last 75 ps of the 250 ps equilibration phase for the systems indicated in the legend. autocorrelation functions that the motion of the molecules becomes uncorrelated after approximately 2 ps. Once equilibrium was reached, simulations were run an additional 200 ps and data collected for FFP

Table 1. Multicomponent Surrogate Table of Compositiona surrogate number 2 molecule count dodecane hexadecane octadecane isooctane isocetane methylcyclohexane decalin dodecene toluene trimethylbenzene tetralin methylnaphthalene total a

32 4 3 12 4 61 42 7 29 18 18 17 247

4 vol % 25 5 5 5 5 15 15 5 5 5 5 5

molecule count 43 0 0 0 0 82 54 0 43 0 33 6 261

11 vol % 37 0 0 0 0 22 21 0 8 0 10 2

molecule count 53 0 0 25 4 63 43 7 0 18 18 0 231

16 vol % 40 0 0 10 5 15 15 5 0 5 5 0

molecule count 0 8 6 24 8 40 28 13 58 0 35 34 254

21 vol % 0 10 10 10 10 10 10 10 10 0 10 10

molecule count 29 0 0 0 4 18 0 12 104 16 48 45 276

28 vol % 25 0 0 0 5 5 0 10 20 5 15 15

molecule count 0 0 3 10 3 17 12 11 98 73 60 0 287

31 vol % 0 0 5 5 5 5 5 10 20 25 20 0

molecule count 11 0 0 0 0 0 0 0 0 88 89 85 273

vol % 10 0 0 0 0 0 0 0 0 30 30 30

Approximate side length (Å): 40. Box volume (Å3): 64 000. 786

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

Figure 2. Structures of the single component hydrocarbons examined here. These hydrocarbons were used to construct the surrogate systems examined here. instrument was cleaned and calibrated using NIST-certified standards. All chemicals were purchased from Sigma-Aldrich with a purity of greater than 94% as measured by gas chromatography. The simulated values of density, heat of vaporization, and bulk modulus are compared to experimentally determined values herein.

■

RESULTS

The accuracy and transferability of the mod-LJ AIREBO potential was examined for higher molecular weight hydrocarbons using samples of pure liquid hydrocarbons, binary mixtures of these hydrocarbons, and multicomponent mixtures of these hydrocarbons termed “multicomponent surrogates”, which see practical use as models of more complex alternative fuels. The structures of the pure hydrocarbons examined here are shown in Figure 2, and the compositions of the multicomponent surrogates are listed in Table 1. Densities calculated for the pure hydrocarbons at 25.0 °C and 1.0 atm with the mod-LJ potential are compared to the experimentally determined values and those obtained using the OPLS-AA potential in Figure 3. When the error bars are considered, the mod-LJ AIREBO potential is able to accurately predict the densities of the straight-chain and branched alkanes. However, it slightly under predicts the densities of the aromatic and cyclic hydrocarbons. The OPLS-AA potential slightly over predicts the densities for the straight-chain and branched alkanes while under predicting the densities of the cyclic alkanes, such as decalin and tetralin, and aromatic hydrocarbons. The root-mean-square deviation (RMSD) of the densities from the experimental values is 0.023 and 0.010 g/cm3 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. The correlation between the calculated densities and the experimentally determined values is shown in Figure 3b. Both potentials do a reasonable job predicting the densities of these hydrocarbons, as indicated by the R2 values of 0.92 and 0.99 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. For comparison, the values calculated by Caleman et al.41 using the OPLS-AA potential for toluene and 1,2,4trimethylbenzene are included. The values calculated here are in reasonable agreement with those calculated by Caleman et al.; however, it should be noted that a different OPLS-AA parameter set was used in that work. The inability of the OPLSAA potential to accurately predict the densities of long-chain hydrocarbons was reported previously.42 Siu et al. optimized

Figure 3. (a) Densities of pure hydrocarbons calculated using molecular dynamics with the OPLS-AA and the mod-LJ AIREBO potentials. The all atom (AA) version of OPLS, as implemented in the LAMMPS MD code,24 with the parameters distributed by Tinker30 was used here. All simulation conditions and the experimental procedure are described in the text. Error bars represent one standard deviation. (b) Densities of pure hydrocarbons calculated using molecular dynamics with the OPLS-AA and the mod-LJ AIREBO potentials versus the experimentally determined densities. For comparison, the values for toluene and 1,2,4-trimethylbenzene computed by Caleman et al.41 are also included in both (a) and (b).

torsional parameters, partial charges, and LJ parameters for methylene hydrogen atoms to improve the properties of long787

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

ΔHvap for the larger hydrocarbons examined here. The experimental ΔHvap values also fall outside the error bars of the predicted OPLS-AA values for the majority of the hydrocarbons examined. The RMSD of the ΔHvap from the experimental values is 12.8 and 7.6 kJ/mol for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. The values of ΔHvap calculated using the OPLS-AA potential for toluene are in reasonable agreement with those calculated by Caleman et al.41 in contrast to trimethylbenzene. Because fossil fuels contain complex mixtures of many hydrocarbons and sometimes have unknown compositions,44,45 it is common for researchers to create surrogate fuel systems to use in laboratory tests20,34,46 and in test engines.13 To better understand the behavior of liquids containing more than one type of hydrocarbon, two separate binary surrogates composed of two types of hydrocarbons and for which thermophysical data are available in the literature were examined. Densities of mixtures of two structural isomers of C16H34, n-hexadecane, and isocetane (or 2,2,4,4,6,8,8-heptamethylnonane) calculated using the mod-LJ potential are compared to the experimentally determined values and the values obtained using the OPLS-AA potential in Figure 5. The values calculated with the mod-LJ AIREBO potential increase as the mole fraction of isocetane increases similar to the experimental values. The calculated density of pure isocetane is in quantitative agreement with the experiment in contrast to the value for pure n-hexadecane. As a result, the calculated densities of the binary mixtures are generally farther from the experimental value as the amount of n-hexadecane is increased. In contrast, the predicted densities calculated using the OPLS-AA potential do not increase with increasing mole fraction of isocetane (Figure 5). Changes brought about by mixing result in better agreement with the experiment. The correlation between the calculated densities and the experimentally determined values yields R2 of 0.51 and 0.88 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. The worst agreement between the experiment and the predicted OPLS-AA densities is for both pure compounds. The poor performance of the OPLS-AA potential reproducing

chain hydrocarbons so that lipid molecules in biological membranes could be studied. Optimized parameters were not developed for other hydrocarbon types, however. Values obtained using the L-OPLS parameter set are discussed below. Predicted values of ΔHvap calculated using the mod-LJ AIREBO potential (Figure 4) show greater deviations from the

Figure 4. Heats of vaporization of pure hydrocarbons calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions are described in the text. The correlation between the calculated ΔHvap and the experimentally determined values yields R2 values of 0.86 and 0.84 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. Error bars represent one standard deviation. Experimental values were taken from the NIST Web book.43 For comparison, the values for toluene and 1,2,4trimethylbenzene computed by Caleman et al.41 are also included.

experimentally determined values43 than the predicted densities, with no clear trend manifesting itself as a function of hydrocarbon type. Interestingly, accurate prediction of hydrocarbon densities or ΔHvap for small hydrocarbon molecules11 does not translate into accurate predictions of

Figure 5. Densities of mixtures of n-hexadecane and isocetane calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. The all atom (AA) version of OPLS, as implemented in the LAMMPS MD code,24 with the parameters distributed by Tinker30 and with the alterations proposed by Siu et al.42 for long-chain hydrocarbons. Other simulation conditions are described in the text. Error bars represent one standard deviation. Experimental values were interpolated to 25 °C from the measurements of Luning Prak et al.34 788

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

Figure 6. Densities of mixtures of n-dodecane and isocetane calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions are described in the text. Error bars represent one standard deviation. Experimental values were interpolated to 25 °C from the measurements of Luning Prak et al.20

Figure 7. Radius of gyration Rgyr as a function of simulation time of samples of liquid isocetane in (a) and n-hexadecane in (b) from simulations using the OPLS-AA potential with the parameters of Tinker. For comparison, Rgyr as a function of simulation time for liquid n-hexadecane calculated using the L-OPLS parameters in shown in (c). Data shown in the red lines were calculated for that component within the binary surrogate with a mole fraction of 0.504 and black lines represent pure liquids.

properties of long alkanes, such as n-hexadecane, was pointed out by Sui et al.42 In that work, an optimized parameter set (LOPLS) was developed for long alkanes. For comparison, densities for the pure components and a 50−50 mixture of nhexadecane and isocetane calculated using this parameter set are also shown in Figure 5. Using this parameter set, the densities are improved but are still not in quantitative agreement with the experimentally determined value. Not unexpectedly, similar trends in the densities of binary surrogates composed of various mixtures of n-dodecane and isocetane calculated using the mod-LJ potential and the OPLSAA potential are apparent in Figure 6. The predicted densities obtained using the mod-LJ AIREBO are in good agreement with the experimentally determined values20 while the values calculated using the OPLS-AA potential are not, with ndodecane demonstrating particularly large deviations from the experiment, which would be lessened by using the L-OPLS parameters.42 Calculation of the radius of gyration Rgyr of the hydrocarbon molecules during a simulation lends insight into the conformation of the molecules in the liquid at equilibrium. Rgyr is how far the average carbon-containing unit is from the

center of mass (COM) of the molecule. It is calculated between all −CH2− or CH3− units within the molecule using the relation Rgyr2 = (1/N)⟨∑kN= 1(rk − rmean)2⟩, where the index k is over N carbon-containing units, rk represents the distance of the kth unit from the COM, rmean is the average distance of the carbon-containing units from the COM of all the simulation molecules, and the brackets ⟨⟩ denote a time average. Rgyr is calculated at 25 fs intervals for every molecule within the liquid simulations of pure n-hexadecane and pure isocetane and for each of these molecules within the 0.504 mol fraction binary mixture of n-hexadecane and isocetane. Histograms of Rgyr are shown in Figure 7. It is clear from the histograms of Rgyr (Figure 7) that there is a distribution of lengths that each molecule adopts in a liquid. In the case of isocetane, histograms of the isocetane molecule are identical whether or not the molecule is in a 50/50 mixture or a pure sample (Figure 7a). In contrast, histograms of nhexadecane in a 50/50 mixture compared to a pure sample differ when calculated using the OPLS parameter set (Figure 7b). In the 50/50 mixture, the distribution is broader and the most probable Rgyr is shifted to lower values. Thus, the molecule is more stretched out on average in the pure 789

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

Figure 8. Density of carbon atoms in the z direction of the simulation cell projected along the xy plane calculated from MD simulations of the binary surrogate with a mole fraction of 0.504 using the OPLS-AA potential. In (a) and (b), data shown were calculated from simulations using the parameters of Tinker, and in (c), the L-OPLS parameters were used.

of the long-chain and branched hydrocarbons. For example, the pure isocetane peaks around 3.3 Å have a higher probability than the peaks for the pure n-hexadecane. At distances above 4 Å, the pure n-hexadecane system is more structured than the pure isocetane (sharper peaks with larger intensity). In other words, when comparing the pure branched versus straight-chain alkane, the branched alkane is less ordered at longer distances. The same pattern in the gCC(r) is evident in the mixtures of ndodecane and isocetane equilibrated with both potentials. The predicted densities of the multicomponent fuel surrogates in Table 1 are shown in Figure 10. There is good

hydrocarbon compared to the mixture. This difference is largely absent when using the L-OPLS parameter set (Figure 7c). As expected, the most probable distance is shorter and the distribution is narrower for the branched alkane than for the long-chain alkane. Additional information about the mixing of n-hexadecane and isocetane can be gleaned from an examination of the density of carbon atoms in the simulation. The carbon-atom density in only the z direction for isocetane in the 50/50 mixture is projected on the xy plane in Figure 8a and for nhexadecane in the 50/50 mixture in Figure 8b is calculated using the OPLS parameters. Comparison of these figures revealed that the areas of high density for isocetane correspond to voids in the n-hexadecane map and vice versa. In other words, the n-hexadecane and isocetane are not uniformly mixed but segregated when using the OPLS parameters for nhexadecane. This segregation is lessened but does not completely disappear when the optimized L-OPLS parameters are used (Figure 8c). The carbon−carbon radial distribution functions (gCC(r)) calculated using VMD and the OPLS-AA potential for nhexadecane, isocetane, and the binary mixtures of these two hydrocarbons are shown in Figure 9. The gCC(r) is an important quantity because it is a measure of the structure of the equilibrated liquids and is accessible experimentally via helium or neutron scattering.47−49 Comparison of distribution functions shown in Figure 9 reveals differences in the structure

Figure 10. Densities of the multicomponent surrogate systems in Table 1 calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions and the experimental procedure are described in the text. Error bars represent one standard deviation. The correlation between the calculated densities and the experimentally determined values yields R2 of 0.99 and 0.85 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively.

qualitative agreement between the experimental and simulated densities for most multicomponent surrogates using both potentials. Importantly, the relative ordering of the densities is reproduced by both potentials; further, the potentials are reasonably accurate over the broad range of densities spanned by these surrogate fuels. It should be noted, however, that the predicted densities are all higher and lower (and outside the

Figure 9. Unnormalized carbon−carbon radial distribution functions (gcc) for the binary mixtures of n-hexadecane and isocetane shown in the legend. These were calculated using VMD. 790

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

0.82 GPa for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. The mod-LJ AIREBO potential performs slightly better when predicting the bulk modulus of the isocetane−dodecane mixtures and the multicomponent surrogates than it did for the pure hydrocarbons with RMSD values of 0.40 and 0.44 GPa, respectively, with the predicted values within one standard deviation of the experimental value for all multicomponent surrogates except surrogate 4. Interestingly, performance of the mod-LJ AIREBO potential degrades when predicting the values of bulk modulus for the isocetane−hexadecane binary with an RMSD of the bulk moduli from the experimental values of 0.96 GPa. Generally, OPLS-AA potential yields better values of bulk modulus than the mod-LJ AIREBO potential. The RMSD values for the bulk moduli computed with the OPLS-AA potential for isocetane− dodecane, isocetane−hexadecane, and the multicomponent surrogates are 0.064, 0.057, and 0.35 GPa, respectively. Surprisingly, the L-OPLS parameter set does not improve the bulk modulus of n-hexadecane as it does for the density values (Figure 13).

error bars) than the experimentally determined densities when calculated with the mod-LJ AIREBO and the OPLS-AA potentials, respectively. The RMSD of the densities from the experimental values is 0.021 and 0.032 g/cm3 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively. The correlation between the calculated densities and the experimentally determined values yields R2 of 0.99 and 0.85 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively, demonstrating the better performance of the OPLS-AA potential predicting densities of the multicomponent surrogates. The bulk modulus was calculated using a technique known as window averaging over the last 200 ps of the simulation.50 In this technique, an average value is calculated in the window of a certain width, tw. The window is then shifted in time by ts and a new average calculated. The reported value of bulk modulus is the average of the values in each segment tw. The value of 1.25 ps for ts was selected because the velocity autocorrelation function has decayed significantly and is close to zero after 1.25 ps (see Figure 1). The selection of the window size tw can impact the calculated values of B, with values of tw that are too small leading to an overestimation of B. The average bulk modulus should be invariant to the selection of tw. After checking several values of tw, it was determined that a value of tw of 75 ps leads to converged values of B in all cases. The bulk moduli of the pure hydrocarbons, the binary systems, and the multicomponent surrogate fuels were calculated, and the results are shown in Figures 11, 12, 13,

■

DISCUSSION The ability to predict the properties of hydrocarbon fuels from their composition via MD simulation depends directly upon the accuracy of the underlying potential model. Here, the accuracy and transferability of the mod-LJ AIREBO potential was examined. The mod-LJ AIREBO potential is a bond-order potential (BOP) that is designed to model chemical reactions in condensed-phase hydrocarbons using a small number of intermolecular parameters. These intermolecular molecular parameters were derived from low molecular-weight hydrocarbon molecules. The functional form of BOPs has been shown to be based on the second-moment approximation.51 According to the moments theorem, the n-th moment of the local density of states on an atom i is determined by the sum of all paths between n neighboring atoms that start and end at atom i. A good estimate of the bond energy can be obtained by knowing only the second moment, which is related to the energy beginning on atom i and summing over the nearest neighbors.52 Thus, the local electronic bond energy for each atom from molecular orbitals is proportional to the square root of the number of neighbors. This approximation was used to derive an expression for the energy of a solid in terms of pairadditive interactions between atoms, and the resulting expression can be rearranged into an expression that is similar to the bond-order potentials. Thus, the bond-order potentials capture the essence of quantum-mechanical bonding. Values calculated using the mod-LJ AIREBO potential were compared to experimental values and values calculated using the OPLS-AA potential, a nonreactive empirical force-field potential. While the form of the OPLS-AA potential cannot be traced to quantum mechanics, the nonreactive nature of the potential makes it computationally faster than the mod-LJ AIREBO potential. In addition, it contains a large number of parameters that could be adjusted to reproduce desired properties, unlike the mod-LJ AIREBO. The mod-LJ AIREBO potential quantitatively reproduces the densities of a wide variety of hydrocarbons and hydrocarbon mixtures, showing transferability to systems for which it was not specifically parametrized. For the set of pure hydrocarbons examined here, the density predictions of the OPLS-AA potential are slightly worse than those of the mod-LJ AIREBO potential when using the Tinker parameters. It should be noted

Figure 11. Bulk moduli of pure hydrocarbons calculated using MD with the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions and the experimental procedure are described in the text. Error bars represent one standard deviation. The correlation between the calculated bulk moduli and the experimentally determined values yields R2 of 0.72 and 0.19 for the OPLS-AA and the mod-LJ AIREBO potentials, respectively.

and 14. In general, the mod-LJ AIREBO potential is unable to make quantitative predictions about modulus and over predicts the values in all cases but one. In contrast, values calculated using the OPLS-AA potential for pure hydrocarbons are within one standard deviation of the experimental values (Figure 11). Note that the error bars shown for the simulated bulk moduli are given as the square root of the average variance of the five unique starting conditions. The RMSD of the bulk moduli from the experimental values for the pure hydrocarbons is 0.26 and 791

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

Figure 12. Bulk moduli of mixtures of n-dodecane and isocetane calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions are described in the text. Error bars represent one standard deviation. Experimental values were interpolated to 25 °C from the measurements of Luning Prak et al.20

Figure 13. Bulk moduli of mixtures of n-hexadecane and isocetane calculated using molecular dynamics with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions are described in the text. Error bars represent one standard deviation. Experimental values were interpolated to 25 °C from the measurements of Luning Prak et al.34

Figure 14. Bulk moduli of the multicomponent surrogates in Table 1 calculated at 1 bar using MD with both the OPLS-AA and the mod-LJ AIREBO potentials. Simulation conditions are described in the text. Error bars represent one standard deviation. Experimental values were measured and simulations carried out at 25.0 °C.

that improvements in density can be obtained by adopting the L-OPLS parameters for long-chain alkanes. For the two binary surrogates, again the densities predicted by the mod-LJ AIREBO potential are slightly better than those from the

OPLS-AA potential using the Tinker parameters. Using the parameters optimized for long alkanes (L-OPLS) results in predictions that are largely comparable to those from the modLJ AIREBO potential. When the multicomponent surrogates 792

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

components for simulations using the L-OPLS parameters, which is surprising since both materials are hydrocarbons. While many potentials for liquid hydrocarbons exist, e.g., OPLS-AA, many of these potentials require a specified molecular topology that does not change during the course of the simulation. That is, they do not permit the making and breaking of bonds. A potential that can offer insight into chemical reactions, such as the mod-LJ AIREBO, as well as reproduce the thermophysical and structural properties of hydrocarbons may offer a more holistic description of fuel behavior because it has the additional benefit of being able to model fuel degradation.

are considered, the densities predicted by both potentials are not quantitative, although the OPLS-AA potential performs statistically better than the mod-LJ AIREBO potential, which is interesting given the discrepancies in the torsional parameters noted by Siu et al.42 It is interesting to note that, despite doing an excellent job predicting the densities of the pure hydrocarbons, the performance of the mod-LJ AIREBO potential degrades when predicting the densities of the multicomponent surrogates, implying that it is not straightforward to extend potentials fit to pure compounds to mixtures. Generally, the predictions of both potentials are slightly less accurate for ΔHvap, which is likely due to an underestimation of the strength of the intermolecular forces in some compounds. However, both potentials still predict the correct qualitative trends. When considering the mod-LJ AIREBO potential, the inability to produce quantitative ΔHvap values for all pure hydrocarbons is not surprising given the fact that all sp3hybridized carbon atoms, regardless of the molecule in which they appear, share the same well-depth parameter εC for the intermolecular LJ terms. In the same way, all sp2-hybridized (and sp-hybridized) carbon atoms utilize the same εC values regardless of the molecule. Moreover, because the intermolecular parameters for all hydrogen atoms are identical, allowing additional variation in σH and εH could improve accuracy. What is surprising is that the mod-LJ AIREBO potential performs as well as the OPLS-AA potential, which has a far greater number of intermolecular parameters. The calculation of the bulk modulus of liquids using MD is challenging and typically accomplished by examining the volume fluctuations, or conversely pressure fluctuations, over long times. The magnitude of volume fluctuations is impacted not only by the potential but also by the simulation size. In principle, MD simulations containing large numbers of molecules would be carried out for nanoseconds or longer to obtain values of bulk moduli. In practice, computational constraints limit system sizes and simulation times with more complicated potential energy functions requiring longer times to evaluate the forces. Thus, to calculate the bulk modulus using MD, a balance must be struck between compositional efficiency and magnitude of the volume fluctuations (i.e., system size). For most of the systems examined here, the values predicted with the OPLS-AA potential are closer to the experimentally determined values than the mod-LJ AIREBO values, with the mod-LJ AIREBO potential consistently over predicting the bulk modulus. Given the fact that both potentials perform similarly when calculating densities and ΔHvap, the differing performance predicting the bulk modulus is the subject of further investigation to determine what parameters in the mod-LJ AIREBO potential influence the predicted bulk moduli. Interestingly, the values of bulk modulus for the multicomponent surrogates predicted by the OPLS-AA potential are in excellent agreement with the experimentally measured values. This improved performance could be due to a cancellation of errors given that many parameter types are present in these systems. This work also makes the prediction that the intensity of the peaks in the radial distribution function can be used to identify changing composition of mixtures of hydrocarbons. Because the radial distribution function can be measured experimentally, this prediction could be validated. In addition, the binary mixtures demonstrate some slight segregation of the

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions ‡

B.L.M., B.H.M., and M.T.K. contributed to the paper equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Office of Naval Research (ONR) under contracts N0001414WX00681 and N0001413WX20347 (J.A.H., K.V.N.), N0001415WX01853 (B.H.M., J.A.H., P.C.T., D.L.P.), N0001413WX20383 (D.L.P., P.C.T.), and N0001413WX20382 (R.E.M.).

■ ■

ABBREVIATIONS MD = molecular dynamics mod-LJ AIREBO = modified Lennard-Jones Adaptive Intermolecular Reactive Empirical Bond-Order OPLS-AA = optimized potential for liquid simulations allatom L-OPLS = long-chain hydrocarbon optimized potential for liquid simulations REFERENCES

(1) Challon, B.; Baranescu, R. Diesel Engine Reference Book; SAE Publishing: Warrendale, PA, 1999. (2) Carr, M. A.; Caton, P. A.; Hamilton, L. J.; Cowart, J. S.; Mehl, M.; Pitz, W. J. In An Experimental and Modeling-Based Study into the Ignition Delay Characteristics of Diesel Surrogate Binary Blend Fuels; ICEF2011; Morgantown, West Virginia, 2011. (3) Edwards, T.; Maurice, L. Q. Surrogate mixtures to represent complex aviation and rocket fuels. J. Propul. Power 2001, 17 (2), 461− 466. (4) Blakey, S.; Rye, L.; Wilson, C. W. Aviation gas turbine alternative fuels: A review. Proc. Combust. Inst. 2011, 33, 2863−2885. (5) Maginn, E. J. From Discovery to Data: What Must Happen for Molecular Simulation to Become a Mainstream Chemical Engineering Tool. AIChE J. 2009, 55 (6), 1304−1310. (6) Payal, R. S.; Balasubramanian, S.; Rudra, I.; Tandon, K.; Mahlke, I.; Doyle, D.; Cracknell, R. Shear viscosity of linear alkanes through molecular simulations: quantitative tests for n-decane and nhexadecane. Mol. Simul. 2012, 38 (14−15), 1234−1241. (7) Aquing, M.; Ciotta, F.; Creton, B.; Fejean, C.; Pina, A.; Dartiguelongue, C.; Trusler, J. P. M.; Vignais, R.; Lugo, R.; Ungerer, P.; Nieto-Draghi, C. Composition Analysis and Viscosity Prediction of Complex Fuel Mixtures Using a Molecular-Based Approach. Energy Fuels 2012, 26 (4), 2220−2230. (8) Fern, J. T.; Keffer, D. J.; Steele, W. V. Vapor-liquid equilibrium of ethanol by molecular dynamics simulation and voronoi tessellation. J. Phys. Chem. B 2007, 111 (46), 13278−13286. 793

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels

order (REBO) potential energy expression for hydrocarbons. J. Phys.: Condens. Matter 2002, 14 (4), 783−802. (28) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, England, 1992; p 118− 122. (29) Mueller, E. A.; Mejia, A. Comparison of United-Atom Potentials for the Simulation of Vapor-Liquid Equilibria and Interfacial Properties of Long-Chain n-Alkanes up to n-C-100. J. Phys. Chem. B 2011, 115 (44), 12822−12834. (30) Tinker Molecular Modeling Package. http://dasher.wustl.edu/ tinker/ (accessed March 1, 2015). (31) Marvin Sketch, Version 16.1.25.0, https://www.chemaxon.com/ download/marvin-suite/, (accessed March 1, 2015). (32) Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30 (13), 2157−2164. (33) Martinez, J. M.; Martinez, L. Packing optimization for automated generation of complex system’s initial configurations for molecular dynamics and docking. J. Comput. Chem. 2003, 24 (7), 819− 825. (34) Luning Prak, D. J.; Trulove, P. C.; Cowart, J. S. Density, Viscosity, Speed of Sound, Surface Tension, and Flash Point of Binary Mixtures of n-Hexadecane and 2,2,4,4,6,8,8-Heptamethylnonane and of Algal-Based Hydrotreated Renewable Diesel. J. Chem. Eng. Data 2013, 58 (4), 920−926. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (36) Shinoda, W.; Shiga, M.; Mikami, M. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69 (13), 134103. (37) Wang, J.; Hou, T. Application of Molecular Dynamics Simulations in Molecular Property Prediction. 1. Density and Heat of Vaporization. J. Chem. Theory Comput. 2011, 7 (7), 2151−2165. (38) Halac, E. B.; Burgos, E. Bulk modulus of amorphous carbon and silicon-carbon multilayered films grown by molecular dynamics simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80 (4), 045426. (39) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987; p 408. (40) Haile, J. M. Molecular Dynamics Simulation: Elementary Methods; John Wiley & Sons, Inc.: New York, 1992; p 489. (41) Caleman, C.; van Maaren, P. J.; Hong, M. Y.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8 (1), 61−74. (42) Siu, S. W. I.; Pluhackova, K.; Bockmann, R. A. Optimization of the OPLS-AA Force Field for Long Hydrocarbons. J. Chem. Theory Comput. 2012, 8 (4), 1459−1470. (43) NIST Web book. http://webbook.nist.gov/cgi/cbook.cgi?ID= C544763&Mask=4. (44) Edwards, T. In Kerosene Fuels for Aerospace Propulsion Composition and Properties, 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Indianapolis, IN, 2002. (45) Maurice, L. Q.; Lander, H.; Edwards, T.; Harrison, W. E. Advanced aviation fuels: a look ahead via a historical perspective. Fuel 2001, 80 (5), 747−756. (46) Luning Prak, D. J.; Cowart, J. S.; Hamilton, L. J.; Hoang, D. T.; Brown, E. K.; Trulove, P. C. Development of a Surrogate Mixture for Algal-Based Hydrotreated Renewable Diesel. Energy Fuels 2013, 27 (2), 954−961. (47) Zimm, B. H. The Scattering of Light and the Radial Distribution Function of High Polymer Solutions. J. Chem. Phys. 1948, 16, 1093− 1099.

(9) Lustig, R. Microcanonical Monte Carlo simulation of thermodynamic properties. J. Chem. Phys. 1998, 109 (20), 8816−8828. (10) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 2000, 112 (14), 6472−6486. (11) Liu, A.; Stuart, S. J. Empirical bond-order potential for hydrocarbons: Adaptive treatment of van der Waals interactions. J. Comput. Chem. 2008, 29 (4), 601−611. (12) NAVAIR. Standard Work Package, Naval fuels & Lubricants CFT Shipboard Aviation Fuel, JP-5, Qualification Protocol for Alternative Fuel/ Fuel Sources; NAVAIR: Patuxent River, MD, 2011. (13) Mathes, A.; Ries, J.; Caton, P. A.; Cowart, J. S.; Luning Prak, D.; Hamilton, L. J. Binary Mixtures of Branched and Aromatic Pure Component Fuels as Surrogates for Future Diesel Fuels. SAE International Journal of Fuels and Lubricants 2010, 3 (2), 794−809. (14) Bruno, T. J.; Baibourine, E. Comparison of Biomass-Derived Turbine Fuels with the Composition-Explicit Distillation Curve Method. Energy Fuels 2011, 25 (4), 1847−1858. (15) Cowart, J. S.; Sink, E.; Slye, P.; Caton, P. A.; Hamilton, L. J. Performance, Efficiency, and Emissions Comparison of Diesel Fuel and Fischer−Tropsch Synthetic Fuel in a CFR Single Cylinder Diesel Engine during High Load Operation. SAE Tech. Pap. Ser. 2008, 200801-2382. (16) Luning Prak, D. J.; Cowart, J. S.; McDaniel, A. M.; Trulove, P. C. Density, Viscosity, Speed of Sound, Bulk Modulus, Surface Tension, and Flash Point of Binary Mixtures of n-Hexadecane plus Ethylbenzene or plus Toluene at (293.15 to 373.15) K and 0.1 MPa. J. Chem. Eng. Data 2014, 59 (11), 3571−3585. (17) Luning Prak, D. J.; Cowart, J. S.; Trulove, P. C. Density, Viscosity, Speed of Sound, Bulk Modulus, and Surface Tension of Binary Mixtures of n-Heptane+2,2,4-Trimethylpentane at (293.15 to 338.15) K and 0.1 MPa. J. Chem. Eng. Data 2014, 59 (11), 3842− 3851. (18) Luning Prak, D. J.; Jones, M. H.; Cowart, J. S.; Trulove, P. C. Density, Viscosity, Speed of Sound, Bulk Modulus, Surface Tension, and Flash Point of Binary Mixtures of 2,2,4,6,6-Pentamethylheptane and 2,2,4,4,6,8,8-Heptamethylnonane at (293.15 to 373.15) K and 0.1 MPa and Comparisons with Alcohol-to-Jet Fuel. J. Chem. Eng. Data 2015, 60 (4), 1157−1165. (19) Luning Prak, D. J.; Morris, R. E.; Cowart, J. S.; Hamilton, L. J.; Trulove, P. C. Density, Viscosity, Speed of Sound, Bulk Modulus, Surface Tension, and Flash Point of Direct Sugar to Hydrocarbon Diesel (DSH-76) and Binary Mixtures of N-Hexadecane and 2,2,4,6,6Pentamethylheptane. J. Chem. Eng. Data 2013, 58 (12), 3536−3544. (20) Luning Prak, D. J.; Alexandre, S. M.; Cowart, J. S.; Trulove, P. C. Density, Viscosity, Speed of Sound, Bulk Modulus, Surface Tension, and Flash Point of Binary Mixtures of n-Dodecane with 2,2,4,6,6Pentamethylheptane or 2,2,4,4,6,8,8-Heptamethylnonane. J. Chem. Eng. Data 2014, 59 (4), 1334−1346. (21) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (22) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225−11236. (23) Jorgensen, W. L.; McDonald, N. A. Development of an all-atom force field for heterocycles. Properties of liquid pyridine and diazenes. J. Mol. Struct.: THEOCHEM 1998, 424 (1−2), 145−155. (24) LAMMPS Molecular Dynamics Simulator. http://lammps. sandia.gov/index.html. (25) Tersoff, J. Modeling Solid-State Chemistry - Interatomic Potentials for Multicomponent Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 39 (8), 5566−5568. (26) Tersoff, J. New Empirical-Approach for the Structure and Energy of Covalent Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37 (12), 6991−7000. (27) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A second-generation reactive empirical bond 794

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795

Article

Energy & Fuels (48) Soper, A. K. The radial distribution functions of water and ice from 220 to 673 K and at pressures up to 400 MPa. Chem. Phys. 2000, 258 (2−3), 121−137. (49) Soper, A. K.; Ricci, M. A. Structures of high-density and lowdensity water. Phys. Rev. Lett. 2000, 84 (13), 2881−2884. (50) Nevins, D.; Spera, F. J. Accurate computation of shear viscosity from equilibrium molecular dynamics simulations. Mol. Simul. 2007, 33 (15), 1261−1266. (51) Brenner, D. W. The art and science of an analytic potential. Phys. Status Solidi B 2000, 217 (1), 23−40. (52) Finnis, M. W.; Sinclair, J. E. A simple empirical N-body potential for transition metals. Philos. Mag. A 1984, 50, 44−55.

795

DOI: 10.1021/acs.energyfuels.5b01468 Energy Fuels 2016, 30, 784−795