Transient molecular dynamics simulations of viscosity for simple fluids
A transient molecular dynamics (TMD) method has been developed for simulation of fluid viscosity. In this method a sinusoidal velocity profile is instantaneously overlaid onto equilibrated molecular velocities, and the subsequent decay of that velocity profile is observed. The viscosity is obtained by matching in a least-squares sense the analytical solution of the corresponding momentum transport boundary-value problem to the simulated decay of the initial velocity profile. The method was benchmarked by comparing results obtained from the TMD method for a Lennard-Jones fluid with those previously obtained using equilibrium molecular dynamics (EMD) simulations. Two different constitutive models were used in the macroscopic equations to relate the shear rate to the stress. Results using a Newtonian fluid model agree with EMD results at moderate densities but exhibit an increasingly positive error with increasing density at high densities. With the initial velocity profiles used in this study, simulated transient velocities displayed clear viscoelastic behavior at dimensionless densities above 0.7. However, the use of a linear viscoelastic model reproduces the simulated transient velocity behavior well and removes the high-density bias observed in the results obtained under the assumption of Newtonian behavior. The viscosity values obtained using the viscoelastic model are in excellent agreement with the EMD results over virtually the entire fluid domain. For simplicity, the Newtonian fluid model can be used at lower densities and the viscoelastic model at higher densities; the two models give equivalent results at intermediate densities.
Transient molecular dynamics simulations of viscosity for simple fluids Jason C. Thomas and Richard L. Rowleya Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, USA Received 22 June 2007; accepted 21 August 2007; published online 6 November 2007 A transient molecular dynamics TMD method has been developed for simulation of fluid viscosity. In this method a sinusoidal velocity profile is instantaneously overlaid onto equilibrated molecular velocities, and the subsequent decay of that velocity profile is observed. The viscosity is obtained by matching in a least-squares sense the analytical solution of the corresponding momentum transport boundary-value problem to the simulated decay of the initial velocity profile. The method was benchmarked by comparing results obtained from the TMD method for a Lennard-Jones fluid with those previously obtained using equilibrium molecular dynamics EMD simulations. Two different constitutive models were used in the macroscopic equations to relate the shear rate to the stress. Results using a Newtonian fluid model agree with EMD results at moderate densities but exhibit an increasingly positive error with increasing density at high densities. With the initial velocity profiles used in this study, simulated transient velocities displayed clear viscoelastic behavior at dimensionless densities above 0.7. However, the use of a linear viscoelastic model reproduces the simulated transient velocity behavior well and removes the high-density bias observed in the results obtained under the assumption of Newtonian behavior. The viscosity values obtained using the viscoelastic model are in excellent agreement with the EMD results over virtually the entire fluid domain. For simplicity, the Newtonian fluid model can be used at lower densities and the viscoelastic model at higher densities; the two models give equivalent results at intermediate densities. 2007 American Institute of Physics. DOI: 10.1063/1.2784117 I. INTRODUCTION Viscosity is an important physical property of a fluid. While viscosity can be measured at many temperatures and pressures, it is not feasible to do so at all conditions required for process design and many of today s applications. For example, the viscosity of lubricants at engine temperatures can be difficult or expensive to obtain. Likewise measurement of viscosity at the variety of process conditions of likely importance is economically unrealistic. A reliable and efficient method for prediction of transport properties over a wide range of conditions is desirable. Molecular dynamics MD simulations is a convenient approach for accurate prediction of vapor and liquid viscosities, because in theory a complete and consistent set of thermophysical properties can be obtained over the entire fluid domain for a substance if an efficacious intermolecular potential model is available. A number of different MD methods have been developed in an effort to predict accurate viscosity values. In the Green- Kubo formalism of equilibrium molecular dynamics EMD simulations, the viscosity is related to the decay of natural fluctuations of mechanical properties occurring in an equilibrium simulation through a time correlation function.1,2 Obtaining an accurate integral over the time correlation function usually requires long simulation times and many averages of the correlation function, and one must be concerned about when the correlation function has actually decayed to zero the so-called long-time-tail issue . Nonequilibrium MD NEMD methods3 6 enhance the signal-to-noise ratio over equilibrium simulations by imposing a velocity gradient, usually linear corresponding to Couette flow, in the simulation cell. This decreases the averaging requirements at any one shear rate, but the gradient is generally extreme relative to experimental possibilities so that multiple simulations must often be used to identify the zero-shear value.6,7 As an alternative to imposing Couette flow, another NEMD method uses an oscillating field and pseudo-steady-state simulation data are used to obtain the viscosity.8 Standard NEMD methods typically obtain the viscosity from the steady-state response to an imposed velocity gradient. Transient molecular dynamics TMD simulations, on the other hand, obtain the viscosity from the transient decay of a system that is instantaneously perturbed from equilibrium. Like NEMD simulations, an imposed gradient gives a stronger signal-to-noise ratio than the decay of natural fluctuations used in EMD simulations. Unlike NEMD simulations, one need only be interested in the initial decay without concern about when the steady state is reached. Furthermore, simulations at multiple shear rates to obtain a low-shear value can be avoided. However, multiple TMD simulations must generally be averaged to obtain the desired statistical accuracy. While we have not made an efficiency comparison for comparable accuracy of results between the various methods, Maginn and co-workers,9,10 found a marked decrease in computational time compared to EMD or NEMD when applying an analogous TMD method. In this paper, we develop a TMD method for viscosity prediction and test the method for Lennard-Jones LJ fluids. a Author to whom correspondence should be addressed. Electronic mail: rowley@byu.edu THE JOURNAL OF CHEMICAL PHYSICS 127, 174510 2007 0021-9606/2007/127 17 /174510/8/$23.00 127, 174510-1 2007 American Institute of Physics Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp While the idea and preliminary implementation of transient methods for transport property predictions are not new,9 12 the approach developed here 1 is compatible with standard simulation algorithms including periodic boundary conditions and Ewald sums, 2 is formulated in terms of a factorable analytic solution to the corresponding boundary-value problem, and 3 uses an initial velocity profile that involves momentum transfer over the entire cell domain to enhance the statistics of the viscosity obtained. II. THEORY In the TMD method, the simulation configuration conforms to that of a physically meaningful solvable boundaryvalue problem with an initial nonequilibrium profile of a property that relaxes toward equilibrium. The transport property governing the profile s rate of decay toward equilibrium is obtained by adjusting it to obtain the best fit usually in a least-squares sense of the solution of the macroscopic boundary-value problem to the corresponding response observed in the simulation. For example, in a TMD method developed by Hulse et al.,12 molecules within a small spherical region of a previously equilibrated simulation cell are instantaneously heated by scaling their velocities. The subsequent transient temperature decay of this spherical region is then matched, by adjusting the value of the thermal conductivity in the analytical solution, to the analytical solution of the macroscopic boundary-value problem corresponding to the temperature decay of a spherical fluid with an initial temperature higher than the infinite, initially isothermal fluid in which it is embedded. Similarly, the TMD method developed by Maginn and co-workers9 11 matches the transient decay of an initial composition or velocity profile to the analytical solutions of the corresponding boundary-value problems for mass and momentum transfer, respectively, by adjusting the diffusion coefficient or viscosity. The method developed in this study is analogous to these previous methods in requiring the transient simulation to conform to a well-defined macroscopic boundary-value problem for which an analytical solution is available. Here we use the simulation cell and coordinates shown in Fig. 1. An initial sinusoidal velocity profile is overlaid on the molecular velocity vectors in a previously equilibrated simulation cell as sketched in Fig. 1. This two-dimensional configuration in Cartesian coordinates with vx y and vy=0=vz meshes well with standard MD simulations performed in Cartesian coordinates. Assuming the fluid is relatively incompressible and that there are neither pressure gradients nor external forces, one can write the equations of continuity and momentum for this configuration as vx t = yx y , 1 where is density, vx is the velocity in the x direction, t is time, and yx is the shear stress. A constitutive equation relating the shear stress to the viscosity is required to complete the formulation of the boundary-value problem. We consider in development of this TMD method only fluids that can be adequately modeled with either a Newtonian or a linear viscoelastic constitutive equation. Restriction of the method to these two classes of fluid guarantees that the solution of the boundary-value problem is separable, but does not significantly restrict application of the simulation method for the model fluids of interest in this and related work. A. Newtonian fluid constitutive equation For the geometry depicted in Fig. 1, the relation between the shear stress and shear rate for a Newtonian fluid is yx = vx t , 2 where is viscosity. Insertion of Eq. 2 into Eq. 1 gives the working partial differential equation that relates the transient behavior of the velocity to the viscosity, vx t = 2vx y2 . 3 In the TMD method developed here, we use an initial velocity profile of the form vx 0,y = vmax cos 2 L y , 4 which exhibits maxima in the velocity vmax at the cell boundaries in the y direction as illustrated in Fig. 1 and conforms to the periodic boundary conditions of the simulation and Neumann-Neumann boundary conditions that require zero shear at y=0 and y=L, namely, vx y y=0 = 0 = vx y y=L . 5 The solution of this boundary-value problem is vx t,y = vmax cos 2 L y exp 2 L 2 t = vx 0,y t , 6 where FIG. 1. Schematic of the simulation cell showing the x and y coordinates with an overlay representation of the initial cosine velocity profile for vx. 174510-2 J. C. Thomas and R. L. Rowley J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp t = exp 2 L 2 t . 7 The initial streaming velocity distribution is thus seen to decay exponentially with time t by the transient function t . The corresponding MD simulation is performed in a box of length L in the y direction. The streaming velocity given by Eq. 4 is overlaid on the instantaneous molecular velocities of an initially equilibrated simulation. The exponential decay in time of the velocity profile is tracked, and the viscosity is obtained by matching in a least-squares sense the analytical solution given in Eq. 6 to the simulated transient data. A convenient form for t that optimizes Eq. 6 to the simulation values can be derived that uses individual particle positions and velocities instead of a continuous profile, thus avoiding binning of molecular velocities and fitting to a velocity profile. The sum of squared residuals is first written as S = i=1 N vxi t vmax cos 2 L yi 2 . 8 To minimize S, the derivative of S with respect to is set equal to zero, and the resultant equation is solved to find the optimum value of at each t. This gives t = i=1 N vxi t /vmax cos 2 /L yi t i=1 N cos2 2 /L yi t . 9 For large enough values of N and under the assumption that particles are evenly distributed along the y axis, the denominator can be treated as a Riemann sum which can be taken as a definite integral in the limit of small particle spacing y. Thus, i =1 N cos2 2 L yi t = 1 y 0 L cos2 2 L y dy = 1 y L 2 = N 2 , 10 where the orthogonality of the trigonometric cosine function has been used to perform the integration. This gives the final expression for the optimized t calculated in terms of individual molecular yi and vxi values at each time step, t = 2 N i=1 N vxi vmax cos 2 L yi . 11 B. Viscoelastic constitutive equation Keshavarzi et al.13 suggest that LJ fluids exhibit viscoelastic behavior at higher densities, for example, when * 0.7, where *= / .3 Here is the number density and is the LJ size parameter. Linear viscoelastic fluids can be represented by a number of different models, but all can be rearranged into the form = t G t t t dt , 12 where G t t is the relaxation modulus. For simplicity the Maxwell equation14 which includes both elastic and viscous effects is used here. For a fluid with flow in the x direction and the velocity gradients in the y direction as in Fig. 1, the chosen constitutive equation is yx + yx t = 0 yx. 13 The variable is the relaxation time and 0 is the viscosity at zero shear. If is zero, the equation for a Newtonian fluid is returned. This differential equation can be solved to give yx = t 0 exp t t t dt . 14 The quantity in brackets is the relaxation modulus. Substitution of Eq. 14 into Eq. 1 gives vx t,y t = t 0 exp t t 2vx t ,y y2 dt . 15 This equation simplifies to t t = a t exp t t t dt , 16 where a = 0 2 L 2 , 17 for the case of a factorable velocity profile of the form v t,y = v t cos 2 L y , 18 as used in the TMD simulations. Equation 16 can be solved using the Laplace transforms for the initial sinusoidal pulse represented by Eq. 4 . The transform of Eq. 16 is s s 0 = a 1 1 + s s , 19 where s is the Laplace transform of in terms of the Laplace variable s. The inverse transform of Eq. 19 gives the solution for t , which can be written in the convenient form t = 1 2 A exp B+t/2 + A+ exp B t/2 , 20 where A = 1 1 4a 1/2, B = 1 1 4a 1/2, 21 a = 0 2 L 2 . As was the case for the Newtonian constitutive equation, Eq. 11 can be used with the transient velocity data obtained 174510-3 Simulations of viscosity J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp from the simulation to find the optimized values for t , in the least-squares sense. These values for t are then used to regress the best values for 0 and using Eqs. 20 and 21 . Although not straightforward, it can be shown that Eq. 20 is equivalent to Eq. 7 in the limit as approaches zero. III. SIMULATION DETAILS Molecular dynamics simulations were performed for LJ fluids using periodic boundary conditions and the minimum image convention. Pair-wise additivity of the LJ pair potentials was assumed, and the LJ potential was truncated at 3.5 . A velocity Verlet algorithm was used with a dimensionless time step of t*=0.003 t*=t /m 2 1/2, where m is the molecular mass and is the LJ well-depth parameter for all but the two highest temperatures studied; a shorter time step of t*=0.001 was used for these two temperatures. Although 256 particles were used in the standard simulations, additional simulations were performed for higher densities with 512 particles by doubling the simulation cell length in the y direction. These larger simulations were performed to help decouple the two parameters 0 and in Eqs. 20 and 21 when analyzing the data using the viscoelastic model. Simulations were made at a large number of state points within the domain 0.8 T* 4.0 and 0.05 * 1.0 as shown in Table I. Here T*=kT/ , where k is Boltzmann s constant and T is the temperature. Standard NVT MD simulations were used to create the initially equilibrated starting configurations at each desired temperature and density. These NVT simulations were allowed to equilibrate for 20 000 time steps after which equilibrated configurations were saved at 1,000 time step intervals to create 80 different equilibrated configurations from which to initiate the transient simulations. An instantaneous cosine velocity profile conforming to Eq. 4 was added to the specific particle velocities of the equilibrated configurations, and standard NVE simulations were then performed allowing the velocity profile to relax. While NVE simulations were used to avoid thermostat effects on velocities, viscous heating of the fluid does occur and is not included in the model used to analyze the transient velocities. Parametric studies of the shear rate suggested that velocity gradients of dvx * /dy* =0.056 for 256-particle simulations and dvx * /dy*=0.028 for 512-particle simulations were optimal for generating a clear transient response without producing shear-rate dependent viscosities that can occur at higher shear rates. Additional simulations were performed for dimensionless densities less than 0.5 using an initial maximum velocity of vx *=0.42 lower initial velocity gradient to examine the effects of viscous heating on the results. Values of vx * t , y generated by the simulation, such as those shown in Fig. 2, were collected as a function of time and converted to optimum t values using Eq. 11 . A heuristic was used to dynamically determine when to terminate collection of the transient data. The time required for t to TABLE I. Lennard-Jones dimensionless viscosity values * and calculated standard deviations in parenthesis obtained from TMD simulations using the Newtonian fluid assumption and Eq. 7 . The open areas correspond to the liquid-vapor two-phase region and the solid-phase region. * T* 0.80 1.00 1.10 1.20 1.25 1.30 1.50 1.80 2.10 2.50 3.00 3.50 4.00 0.05 0.104 0.009 0.110 0.005 0.124 0.010 0.126 0.008 0.137 0.008 0.155 0.011 0.182 0.013 0.200 0.032 0.231 0.033 0.273 0.067 0.241 0.024 0.10 0.133 0.008 0.140 0.012 0.143 0.014 0.164 0.018 0.180 0.006 0.204 0.013 0.234 0.027 0.276 0.027 0.313 0.044 0.323 0.025 0.15 0.174 0.011 0.194 0.015 0.222 0.012 0.244 0.018 0.276 0.027 0.321 0.026 0.325 0.025 0.356 0.050 0.20 0.203 0.013 0.220 0.006 0.262 0.016 0.280 0.022 0.311 0.022 0.368 0.036 0.400 0.048 0.399 0.075 0.30 0.297 0.020 0.309 0.014 0.355 0.041 0.372 0.046 0.445 0.038 0.485 0.056 0.460 0.072 0.40 0.410 0.039 0.445 0.054 0.484 0.032 0.489 0.041 0.515 0.066 0.614 0.091 0.622 0.082 0.50 0.559 0.034 0.585 0.053 0.579 0.037 0.628 0.085 0.636 0.072 0.728 0.106 0.770 0.101 0.734 0.126 0.60 0.805 0.060 0.803 0.066 0.804 0.079 0.840 0.087 0.851 0.097 0.878 0.097 0.876 0.103 0.881 0.101 0.945 0.144 1.013 0.124 0.70 1.167 0.115 1.186 0.122 1.237 0.084 1.291 0.084 1.212 0.098 1.200 0.133 1.191 0.190 1.276 0.142 1.341 0.199 1.330 0.121 1.347 0.193 0.80 2.516 0.288 2.436 0.275 2.231 0.300 2.282 0.228 2.272 0.164 2.107 0.179 1.893 0.274 1.843 0.313 1.966 0.208 2.029 0.333 1.780 0.255 1.949 0.274 1.681 0.373 0.85 4.359 0.366 3.639 0.351 3.241 0.261 3.310 0.297 3.142 0.270 2.989 0.346 2.693 0.271 2.550 0.203 2.744 0.483 2.505 0.365 2.256 0.313 2.151 0.433 2.194 0.427 0.90 5.556 0.619 5.252 0.378 4.611 0.337 4.395 0.325 4.541 0.487 3.953 0.583 3.510 0.473 3.184 0.296 3.248 0.531 2.777 0.444 2.959 0.504 2.931 0.740 0.95 5.951 0.593 5.155 0.551 4.566 0.459 4.117 0.598 3.458 0.618 3.309 0.726 3.630 0.640 1.00 7.023 0.876 6.057 0.717 5.451 0.809 5.087 0.987 4.290 1.051 4.137 0.963 174510-4 J. C. Thomas and R. L. Rowley J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp decay from 0.5 to 0.25 was monitored, and the simulation was continued beyond the t =0.25 time for nine times this time interval. Stored values for each of the 80 t transient decays were averaged to obtain a smoother data set, and ten replicates of these smoother data sets were generated at each state point. A single regression was performed to obtain the predicted viscosity at each state point from the ten replicate data sets 800 transient decays , but separate regressions on the individual replicates were also performed to estimate the standard deviation of the viscosity value. Three examples of fitting the transient t decay to the macroscopic analytical equations for the two different constitutive models, Eqs. 7 and 20 , are shown in Fig. 3. At lower densities, roughly * 0.3, t is well represented with the Newtonian model. At high densities, * 0.7, t can exhibit a negative region which the viscoelastic model can reproduce but the Newtonian model cannot. At intermediate densities, 0.3 * 0.7, either model can be used to fit the transient response data. A chi-squared test was used to determine the quality of the regressions. In all cases, the observed decay retained the cosine shape conformal to the original applied velocity distribution. As can be seen from Fig. 3, t for low densities and very short times t* 0.5 is inconsistent with the purely exponential decay of the Newtonian model given in Eq. 7 . This short-time inconsistency, prior to the purely exponential decay of t , correlates roughly with the inverse of the collision frequency that one would calculate for hard spheres with the same value of , suggesting that multiple molecular collisions may be required before the continuum equations apply. For hard spheres, simple kinetic theory1 predicts v* = 8T* / for the expectation of the velocity and l* =1/ 2 * for the mean free path. A straightforward estimate of the collision time is therefore tcoll * = l* v* = 1 4 * T* . 22 In regressing low-density data with the Newtonian model, we have not used the early simulation data before 1.4tcoll * to accommodate this molecular relaxation period. IV. RESULTS Table I shows the regressed dimensionless viscosity values * *= 2 / m 1/2 and the calculated standard deviations obtained over the whole T domain using the Newtonian constitutive equation and the 256-particle simulation results. Figure 4 compares these results to the correlation FIG. 2. Typical decay of the initial velocity profile of the molecules in the simulation cell. This particular response is for *=0.95 and T*=1.5. FIG. 3. Transient t data from the simulation solid black line and the fitted macroscopic equation based on the Newtonian model and on the viscoelastic model solid gray line at low density upper figure: * =0.05, T*=1.5 , moderate density middle figure: *=0.4, T*=1.5 , and high density lower figure: *=0.95, T*=1.5 . The two sets of curves for the moderate- and high-density plots are for the two simulation sizes of 256 lower curves and 512 particles upper curves . FIG. 4. Percent difference between simulated LJ viscosity values obtained in this study using the TMD method based on the Newtonian fluid assumption and a previous correlation Ref. 15 at various temperatures: T*=0.8 , gray fill , 1.0 , gray fill , 1.1 , 1.2 , 1.25 , 1.3 , 1.5 , 1.8 , 2.1 , 2.5 , 3.0 , 3.5 , and 4.0 . 174510-5 Simulations of viscosity J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp based on EMD values reported by Rowley and Painter.15 We have also made a similar comparison not shown to EMD viscosity values reported by Woodcock16 for larger simulation cells. That comparison gives similar results except at very large densities near the solid-liquid transition line where the simulated viscosity appears to have a size dependence. Extensive and very accurate EMD simulation results 3% for liquid densities, increasing to 10% for low-density gases were recently reported by Meier et al.17 The assumption of Newtonian behavior yields good agreement of the predicted viscosity values with the EMD results over a large range of fluid conditions 0.1 * 0.7 , but the systematic deviation between the two methods is evident at dimensionless densities above 0.8. This increasingly positive deviation from the EMD values with increasing density corresponds with the transient behavior of t that we observed above a dimensionless density of 0.7, namely, that t exhibited increasingly larger negative regions with increasing density for * 0.7. This behavior is attributed to viscoelastic behavior. Table II provides the * values and their standard deviations obtained when the simulation results are analyzed in terms of the viscoelastic model, the positive systematic bias observed at the higher densities is removed. Figure 5 compares the viscosity values obtained using the viscoelastic model to the same EMD results used in Fig. 4. Figure 6 shows a comparison of the TMD results to the more recent Meier EMD values.17 Both Figs. 5 and 6 suggest a slightly negative bias of the TMD results relative to the EMD results, but the agreement is well within the uncertainty of the two data sets. Values of the viscosity in the region 0.3 * 0.7 obtained by both the Newtonian and viscoelastic models were the same within the uncertainty of the simulated results. As mentioned above, the use of the viscoelastic model provided very good fits of all of the transient data and is seen to remove the high-density bias. Viscosity predictions show a negative systematic bias in Fig. 4 at very low densities * 0.1 . The TMD method is less well suited for prediction of gas viscosities at low density. This is likely because the low relatively speaking molecular collision frequency in dilute gases makes momentum transport less easily modeled with the continuous equations. TABLE II. Lennard-Jones dimensionless viscosity values * and calculated standard deviations in parenthesis obtained from TMD simulations using the viscoelastic model and Eq. 20 . The open areas correspond to the liquid-vapor two-phase region and the solid-phase region. * T* 0.80 1.00 1.10 1.20 1.25 1.30 1.50 1.80 2.10 2.50 3.00 3.50 4.00 0.05 0.10 0.15 0.172 0.006 0.193 0.009 0.220 0.010 0.243 0.011 0.277 0.013 0.313 0.018 0.343 0.021 0.370 0.048 0.20 0.201 0.010 0.219 0.010 0.255 0.009 0.280 0.035 0.299 0.013 0.335 0.016 0.374 0.014 0.392 0.022 0.30 0.301 0.040 0.332 0.044 0.352 0.021 0.393 0.047 0.422 0.029 0.471 0.039 0.494 0.051 0.40 0.401 0.013 0.436 0.010 0.457 0.030 0.484 0.055 0.511 0.062 0.571 0.128 0.586 0.079 0.50 0.548 0.077 0.590 0.142 0.588 0.032 0.615 0.038 0.637 0.150 0.654 0.086 0.732 0.175 0.749 0.106 0.60 0.788 0.089 0.786 0.092 0.790 0.021 0.801 0.033 0.829 0.102 0.831 0.051 0.854 0.043 0.944 0.243 0.925 0.260 1.008 0.122 0.70 1.203 0.063 1.214 0.082 1.183 0.137 1.180 0.048 1.164 0.047 1.166 0.142 1.235 0.299 1.221 0.156 1.262 0.080 1.293 0.181 1.314 0.186 0.80 2.164 0.064 2.025 0.071 1.961 0.077 1.930 0.079 1.948 0.059 1.947 0.085 1.867 0.097 1.833 0.213 1.813 0.124 1.876 0.275 1.802 0.053 1.761 0.432 1.871 0.264 0.85 3.007 0.106 2.719 0.114 2.598 0.167 2.497 0.105 2.564 0.111 2.480 0.098 2.391 0.103 2.279 0.128 2.260 0.128 2.203 0.159 2.248 0.205 2.127 0.207 2.099 0.549 0.90 3.801 0.081 3.602 0.146 3.399 0.157 3.482 0.118 3.243 0.191 3.119 0.187 2.868 0.167 2.851 0.175 2.667 0.179 2.602 0.373 2.751 0.199 2.622 0.275 0.95 4.318 0.244 3.981 0.233 3.705 0.175 3.541 0.240 3.295 0.381 3.165 0.408 3.146 0.746 1.00 5.178 0.401 4.898 0.282 4.273 0.355 4.236 0.330 3.963 0.264 3.822 0.453 FIG. 5. Percent difference between simulated LJ viscosity values obtained in this study using the TMD method based on the viscoelastic assumption and a previous correlation Ref. 15 at various temperatures: T*=0.8 , gray fill , 1.0 , gray fill , 1.1 , 1.2 , 1.25 , 1.3 , 1.5 , 1.8 , 2.1 , 2.5 , 3.0 , 3.5 , and 4.0 . 174510-6 J. C. Thomas and R. L. Rowley J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp Additional simulations with longer times and larger system sizes suggest that accurate results can be obtained at dimensionless densities lower than 0.1, but we have not pursued that avenue further because of the decreased efficiency of the method in that region and because the Chapman-Enskog theoretical equations1 are accurate and easier to use in this density range. Potential errors in the results due to viscous heating are largest for low densities and low temperatures. The final temperature of the transient data permits calculation of an upper bound on the error due to the viscous-heating temperature rise. Using the temperature dependence of * given in Ref. 15, this upper bound on the error is less than 1% for the vast majority of the results. Even for points at lower T* and/or *, the maximum error is less than 4% for the viscoelastic model. Simulations repeated at lower T* and * conditions with smaller initial gradients produced equivalent * values to those in Table II within the standard deviations shown, but the uncertainties were larger. An upper-bound estimate of the possible effect of viscous heating on the Newtonian-model results suggested that larger errors, up to perhaps 12%, might be possible at the extreme low-density limit. However, again repeated simulations using the lower initial velocity gradients were nearly indistinguishable from those with initial gradients, suggesting that the upper-bound analysis significantly overestimates viscous heating errors. It should be noted that there are three missing entries in Tables I III at very high densities near the liquid-solid phase boundary. These values were omitted from these tables because the chi-squared statistics indicated an inability to adequately fit the simulated decay of t even with the viscoelastic equations. While the response from either the 256- or 512-particle systems could be fitted individually, both responses could not be simultaneously regressed suggesting an inadequacy in the linear viscoelastic model near the solidliquid phase boundary. The two sizes of simulation cells were used both to ensure size consistency of the model and to provide good decoupling of the two parameters in the FIG. 6. Percent difference between simulated LJ viscosity values obtained in this study using the TMD method based on the viscoelastic assumption and EMD simulated results Ref. 17 at various temperatures: T*=0.8 , gray fill , 1.0 , gray fill , 1.1 , 1.2 , 1.25 , 1.3 , 1.5 , 1.8 , 2.1 , 2.5 , 3.0 , and 4.0 . TABLE III. Values of *, the viscoelastic relaxation time, and its standard deviation in parentheses as determined from the TMD simulations using the viscoelastic assumption and Eq. 20 to simultaneously fit the transient response in the two different sized cells. * T* 0.80 1.00 1.10 1.20 1.25 1.30 1.50 1.80 2.10 2.50 3.00 3.50 4.00 0.05 0.10 0.15 0.501 0.204 0.519 0.199 0.537 0.271 0.426 0.194 0.355 0.091 0.435 0.162 0.343 0.098 0.361 0.221 0.20 0.459 0.159 0.390 0.243 0.323 0.167 0.337 0.225 0.344 0.122 0.364 0.161 0.369 0.156 0.254 0.138 0.30 0.232 0.180 0.179 0.170 0.176 0.120 0.125 0.182 0.197 0.089 0.184 0.136 0.162 0.125 0.40 0.176 0.073 0.154 0.157 0.159 0.117 0.152 0.112 0.119 0.127 0.109 0.147 0.195 0.207 0.50 0.110 0.114 0.058 0.091 0.080 0.125 0.137 0.071 0.105 0.131 0.070 0.105 0.021 0.101 0.075 0.104 0.60 0.030 0.064 0.049 0.046 0.092 0.039 0.081 0.051 0.033 0.046 0.068 0.066 0.141 0.095 0.009 0.070 0.033 0.114 0.013 0.106 0.70 0.039 0.029 0.035 0.034 0.060 0.049 0.062 0.032 0.047 0.041 0.038 0.041 0.044 0.053 0.042 0.044 0.066 0.052 0.068 0.050 0.021 0.072 0.80 0.120 0.016 0.103 0.032 0.092 0.028 0.074 0.027 0.063 0.036 0.072 0.032 0.073 0.042 0.037 0.038 0.067 0.041 0.022 0.031 0.037 0.029 0.022 0.055 0.036 0.047 0.85 0.170 0.012 0.131 0.031 0.111 0.018 0.097 0.028 0.080 0.019 0.081 0.018 0.074 0.035 0.066 0.034 0.055 0.021 0.049 0.033 0.037 0.033 0.040 0.051 0.020 0.061 0.90 0.157 0.022 0.148 0.024 0.126 0.027 0.111 0.023 0.115 0.019 0.091 0.040 0.061 0.024 0.054 0.032 0.058 0.041 0.037 0.042 0.052 0.026 0.040 0.022 0.95 0.129 0.018 0.115 0.022 0.071 0.028 0.067 0.023 0.037 0.028 0.064 0.028 0.049 0.052 1.00 0.112 0.025 0.095 0.029 0.085 0.019 0.053 0.024 0.053 0.035 0.045 0.011 174510-7 Simulations of viscosity J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp viscoelastic model. Values of the dimensionless viscoelastic relaxation time * obtained and their standard deviations are reported in Table III. V. CONCLUSIONS A TMD method has been developed for prediction of fluid viscosity. The simulation produces transient velocities from an initial cosine velocity profile. These transient velocities are fitted to a theoretical analytical solution of the boundary-value problem describing the macroscopic momentum transfer for the geometry and initial condition used in the simulation. The viscosity is adjusted in this procedure to give the best least-squares fit of the simulation data to the macroscopic equation. To obtain the macroscopic equation explicit in the viscosity, a constitutive equation is used to relate the stress tensor to the velocity gradient. In this work we have used two different constitutive relations. If the fluid is considered to be Newtonian, the analysis equation is capable of fitting the transient response observed from the simulation only for dimensionless densities less than 0.7. At higher densities, the velocity decay shows a region of velocity reversal indicative of viscoelastic behavior. Correspondingly, the resultant viscosity values obtained from the method agree with previous EMD results in the density range 0.1 * 0.7. The TMD method as developed here is not recommended for dimensionless densities below 0.1, where the frequency of molecular collisions is low. Increased simulation cell sizes can be used to ameliorate this problem. This does not limit the method because simulations are generally not required in the low-density region where theoretical equations can be used to predict dilute gas viscosities. At higher densities, the transient velocity decay is well represented by the equations obtained using a linear viscoelastic model to relate the shear stress to the velocity gradient. This model adequately represents the regions of reversed velocity flow observed in the transient velocity decay for * 0.7. We have used two cell sizes in our simulations and required the model to adequately fit both transient responses to ensure a decoupling of the two parameters in the viscoelastic model. The resultant viscosity values agree well with EMD results over the density domain * 0.3 and also with the TMD results obtained using the Newtonian fluid assumption in the region 0.3 * 0.7. However, the transient responses in the two different size systems could not be fitted simultaneously using the viscoelastic model for three densities along the liquid-solid phase transition line. The TMD method developed here can be used to simulate fluid viscosity for dimensionless densities above 0.1. The viscoelastic model can be used for densities above 0.3, but must be used to obtain accurate results at dimensionless densities greater than 0.7. The velocity relaxation in the LJ fluid is very fast and data are collected over a dimensionless time of 10 20, corresponding to approximately 20 40 ps for Ar. An average of 800 of these velocity relaxations produces a viscosity estimate with an uncertainty of less that 10%. ACKNOWLEDGMENTS Partial funding of this project through Design Institute for Physical Properties DIPPR Project 801 is gratefully acknowledged. 1 D. A. McQuarrie, Statistical Mechanics Harper & Row, New York, 1976 . 2 J. M. Haile, Molecular Dynamics Simulation: Elementary Methods Wiley, New York, 1992 . 3M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids Clarendon, Oxford, 1987 . 4 D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids Academic, New York, 1990 . 5 D. J. Evans, W. G. Hoover, B. H. Failor, B. Moran, and A. J. C. Ladd, Phys. Rev. A 28, 1016 1983 . 6 R. Edberg, D. J. Evans, and G. P. Morriss, J. Chem. Phys. 84, 6933 1986 . 7 S. T. Cui, S. A. Gupta, P. T. Cummings, and H. D. Cochran, J. Chem. Phys. 105, 1214 1996 . 8F. Muller-Plathe, J. Chem. Phys. 106, 6082 1997 . 9 G. Arya, E. J. Maginn, and H. C. Chang, J. Phys. Chem. 113, 2079 2000 . 10M. S. Kelkar and E. J. Maginn, J. Chem. Phys. 123, 224904 2005 . 11 E. J. Maginn, A. T. Bell, and D. N. Theodorou, J. Phys. Chem. 97, 4173 1993 . 12 R. J. Hulse, R. L. Rowley, and W. V. Wilding, Int. J. Thermophys. 26, 1 2005 . 13 E. Keshavarzi, M. Vahedpour, S. Alavi, and B. Najafi, Int. J. Thermophys. 25, 1747 2004 . 14 R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Dynamics of Polymeric Liquids, 2nd ed. Wiley, New York, 2002 . 15 R. L. Rowley and M. M. Painter, Int. J. Thermophys. 18, 1109 1997 . 16 L. V. Woodcock, AIChE J. 52, 438 2005 . 17 K. Meier, A. Laesecke, and S. Kabelac, J. Chem. Phys. 121, 3671 2004 . 174510-8 J. C. Thomas and R. L. Rowley J. Chem. Phys. 127, 174510 2007 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp