Comparison of the Accuracy and the Efficiency among Different Molecular Dynamics Methods for Calculating the Viscosity

The viscosity is one of the important parameters which determine the characteristics of the fluid flow. In this paper, we employed three molecular dynamics methods so-called the periodic poiseuille flow method, the shear flow method and the stress autocorrelation method to calculate the viscosity of the Lennard-Jones fluid. We have found that, among the three, the periodic poiseuille flow method has the greatest convergent speed and the highest accuracy; the shear flow method with a small shear rate has the smallest convergent speed. The accuracy of the stress autocorrelation method evidently depends on the computational time. With 10 time-steps it is as accurate as the shear flow method, but with double time-steps it is more accurate than the shear flow method.


Introduction
With the developments of the micro-and nanotechnologies, the fluid flow in micron or nano scale is paid more and more attentions in many scientific and industrial fields.The viscosity is one of the important parameters which determine the characteristics of the fluid flow and thus greatly affects the accuracy of the theoretical and engineering calculations.However, there is no exact theoretical expression for calculating the viscosity [1].An experimental measurement or a numerical computation is always performed to get the viscosity value.The way of experimental measurements has the drawbacks of high cost and complex implementing procedure.As the developments of physics and the technologies concerning to computing, the numerical methods are put into realities with considering more complex models with high efficiency and high accuracy.So the way of adopting numerical computations has come to be one of the most effective methods for calculating the viscosity.The numerical models mainly include the ones based on the molecular dynamics [2], the dissipative particle dynamics [3] and the Lattice Boltzmann algorithms [4], [5].
Different numerical methods based on the molecular dynamics of calculating the viscosity can be divided into two categories so-called the equilibrium and the nonequilibrium methods.For equilibrium methods, the viscosity is obtained according to the stress-stress autocorrelation function [6] in an equilibrium state without the disturbance of external forces.In nonequilibrium methods, the external forces are exerted to the system and the non-equilibrium steady state of the fluid flow reflects the effects of the viscosity which can be obtained by calculating the shear stress.The common ways to form steady shear flows include the method with a sliding boundary and the one with periodic shear flow [6].J. A. Baker [7] proposed the periodic poiseuille flow method in which the periodic boundary condition was adopted as a substitute of the wall boundary condition [8].The problems of density oscillation and velocity slippage near the walls were solved out in his method, and the viscosity could be obtained by calculating the average flow velocity.This paper discusses the three methods of the periodic poiseuille flow, the shear flow and the stress autocorrelation in more details.Each of them is applied to calculate the viscosity of the simple Lennard-Jones fluid and the three methods are evaluated through comparing their results.

The periodic poiseuille flow method
In molecular dynamics simulation, the poiseuille flow can be obtained by exerting a body force to the system, that is, a corresponding force to every particle.For example, two plates parallel to the xz-plane are located with the distance D and to the fluid between them is applied a body force in x-direction.The resulting flow field has a parabolic velocity pofile Integrating Equation ( 1) and the expression of the average velocity through the yz-plane is given as where ² ¢" denotes the ensemble average.The viscosity can be calculated by Equation (2).
In order to get a zero velocity upon the boundaries, the wall model is used, which yet leads to the density oscillation and the velocity slippage near the walls.J. A. Baker [7] proposed the method of periodic poiseuille flow, in which the length of the fluid region in y-direction is increased to be 2D, and the walls are replaced by the periodic boundary conditions.The system is divided into two domains by the symmetric plane in y-direction and each is exerted a body force in the opposite directions.Thus a zero velocity on the boundaries and on the symmetric plane can be obtained; and two opposite parabolic velocity profiles are formed at the both sides of the symmetric plane.

The shear flow method
The shear flow method is also a kind of non-equilibrium method in which a steady shear rate should be reached by the common method of sliding wall.For example, in a system, two parallel plates of the xz-plane are located in a distance H and the upper wall moves in a velocity V shear while the lower wall is fixed.When the system achieves the steady state a stable Couette flow is formed.The resulting flow field has a linear velocity profile.

( )
In the steady state, the fluid deforms because the velocities in every layer are different, which results in the production of the inner friction.According to Newton's inner friction law, in a laminar flow, the shear stress which can be given by the elements of the pressure tensor should be in direct proportion to the velocity gradient at a certain point [9].The law can be represented as ( ) After putting the Equation (3) into the Equation ( 4), the expression of the viscosity can be formulated as where xy p denotes the xy component of the pressure tensor.The method can be used for the non-Newtonian fluids.However, it would take a long time to finish the simulation to get the viscosity when the shear rate is small.

The stress autocorrelation method
The stress autocorrelation method is a kind of equilibrium method.The viscosity is related to the oscillation of the off-diagonal elements of the pressure tensor.In an equilibrium simulation, the viscosity is obtained by integrating the Green-Kubo equation [6], which can be represented as 0 (0) ( ) where ² ¢" denotes the ensemble average, V is the volume, k B is the Boltzmann constant, T is the absolute temperature, t is the time and p xy is the xy component of the pressure tensor.
The method can be used for calculating the viscosity of non-Newtonian fluid in zero shear rate.

Computational details
In this paper, the three methods are applied to the simple Lennard-Jones fluid-argon.The Lennard-Jones potential function used in the simulation, is shown as Eq.(7).
The system is non-dimensionalized by using the reduced units [2], which are widely used in the molecular dynamics simulation.For all the models, the energy parameter and the length parameter among the atoms in the fluid are all set to be 1.0, i.e., ll V and ll .The neighbour list is updated every time step with a cut-off radius of 2.5.The time step is set to be 0.00233.The size of the simulation box is shown in Table 1.The periodic boundary condition is applied to all the three directions in the periodic poiseuille flow method and the stress autocorrelation method, but only to the xand z-directions in the shear flow method.When the viscosity is calculated by the non-equilibrium methods, the viscous heating would be brought to the system by the external force.Hence, a thermostat should be chosen to keep the temperature of the system constant.The DPD thermostat is used with the advantage of Galilean invariance which does not disturb the hydrodynamic interaction [10].For the stress autocorrelation method, the DPD thermostat should also be applied to the system in equilibrium to make sure that the three methods use the same thermostat.In this way, the results of the three methods can be compared.
For the periodic poiseuille flow method, different external forces are applied to the system and the relevant temperature profiles are shown in Fig. 1.For example, the minimum temperature of the system with the body force 0.12 is about 0.74, which is 2% more than the input value of 0.722.It is clear that the temperature of the system rises with increasing the external forces.For the simple DPD algorithm, the temperature of the system depends strongly on the time step [11].The DPD thermostat is not suitable for a system with the large external force because it means that a very small time step should be taken.Moreover, the computation with very small time step would cost too much time to assure the system reaching the state we want.Correspondingly, for a given time step, the body force should not be too large.In this paper, the body force 0.02 is chosen and the relevant temperature profile is shown in Fig. 2. The deviation from the expected value is very small and the viscous heating effect can be neglected.The similar researches are undertaken when we choose the shear rate in the shear flow method.Because of the artifact of the simple DPD algorithm, the small shear rate should be chosen.
In short, the body force is set as 0.02 for the periodic poiseuille flow method, and the velocity of the upper wall is set as 0.04 for the shear flow method.

Viscosity of a Lennard-Jones fluid
For the periodic poiseuille flow method, there are 2313 fluid atoms which are initially located in the FCC lattice in the cubic box; and the velocities of the atoms are set randomly according to the given temperature.When the system reaches its steady state, the velocity distribution of the particles in the system will become the steady Maxwell distribution.At the beginning of the simulation, the NVE ensemble and the DPD thermostat are applied to the system.The body force is not applied to the system until it reaches the equilibrium state.The velocities of the atoms are corrected every 100 time steps to make sure that the momentum of the system is zero.The velocity profile is shown in Fig. 3.The measurements are accomplished in each part and are fitted with a quadratic function by the least square method with the correlation coefficients larger than 0.97.The velocity profile agrees well with the analytic laminar poiseuille flow pattern from Eq. (1).

Figure 3 Velocity profile of the periodic poiseuille flow
The velocities are output every time step and the results are used to calculate the viscosity every 5  10 time steps.The process is repeated ten times.The dimensionless viscosity is 18 .0 47 .3 r .For the shear flow method, there are 1324 wall atoms and 2321 fluid atoms in the cubic box.The rigid wall has a FCC lattice.Initially, the fluid atoms are located in the FCC lattice and the velocities are set randomly according to the given temperature.At the beginning of the simulation, the NVE ensemble and the DPD thermostat are applied to the system and the upper wall moves at a constant velocity of 0.04.The viscosity is not calculated until the system reaches the steady state.The xy component of the pressure tensor is output every time step and the results are used to calculate the viscosity every 10 7 time steps.The process is repeated ten times.The method has a slow convergence and it takes about 10 7 time steps to get a reasonable result.The reason is that the shear rate we have adopted is very close to zero .The velocity profile in Fig. 4 shows that the steaming velocity of the particles near the fixed wall in the distance range 0 to 7 in y direction is almost zero, deviating from the theoretical values largely.A stable shear rate is formed far away from the fixed wall, which is larger than the theoretical one from Equation (3).But for all that, the method can get a viscosity comparable with that obtained by using the periodic poiseuille flow method.

Comparison of the three methods
The viscosities calculated by the three methods are in agreement, and also accord with that by J. A. Baker [7].All the results are shown in Table 2.

1.2 r
In terms of the convergent speed, the reasonable results can be obtained with a simulation time of 10 5 time steps by the periodic poiseuille flow method and the stress autocorrelation method.However, a simulation time of about 10 7 time steps is required to get a reasonable result for the shear flow method.The shear flow method has the slowest convergent speed, while the periodic poiseuille flow method and the stress autocorrelation method have comparable ones.In terms of the accuracy, the periodic poiseuille method gives a more accurate result which coincides well with the viscosity of 3.4 0.2 r given by J. A. Baker.The result obtained by the stress autocorrelation method with 10 5 time steps is as accurate as that obtained by the shear flow method, while the result with time steps is more accurate than the one by the shear flow method.In short, in terms of convergent speed and accuracy, the periodic poiseuille flow is better than the methods of stress autocorrelation and shear flow.

Conclusion
In this paper, the three methods for calculating the viscosity with molecular dynamics simulation have been compared.The methods differ in applicability, accuracy and computational time.For a Lennard-Jones fluid, the three methods can give comparable results.
For the periodic poiseuille flow method, it can only be used to calculate the viscosity of the Newtonian fluid.However, both of its convergent speed and accuracy are the best among the three.For the shear flow method, it can be used to calculate the viscosity of the non-Newtonian fluid.When the shear rate is very small, its convergent speed is the slowest among the three methods.For the stress autocorrelation method, it can be used to calculate the non-Newtonian fluid with a zero shear rate.With a computational time of 10 5 time steps, it is as accurate as the shear flow method, but while with double time steps it is more accurate than the shear flow method.

H
are both equal to 1.0.For the shear flow method the model includes the rigid walls in which the atoms are located in the structure of face-centered cubic (FCC) with a lattice constant of 2.467.The energy parameter and the length parameter among the atoms in the walls and those in the fluid are set to be 1.0, i.e., wl V and wl H are both equal to 1.0, ensuring that the liquid wets the wall.In all the simulations, the number density of the particles is

Figure 1 .Figure 2
Figure 1.Comparison of the temperature profiles under different body forces

Figure 4
Figure 4 The velocity profile of the shear flow For the stress autocorrelation method, there are 2313 fluid atoms which are initially located in the FCC lattice in the cubic box and the velocities of the atoms are set randomly according to the given temperature.At the beginning of the simulation, the NVE ensemble and the DPD thermostat are applied to the system.The viscosity can be calculated when the system reaches the equilibrium state for it is essentially an equilibrium method.Two simulations with different computation time are performed.One has 10 5 time steps and the other has double.Both simulations are repeated ten times.The viscosities this paper have gotten from the two simulations are 90 .0 97 .3 r and 67 .0 72 .3 r , respectively.

Table 1 .
The size of the simulation box

Table 2 .
Comparison among our results and J. A. Baker's