The article information
 Ruoyu Dong, Bingyang Cao, Heming Yun, Baoming Chen
 董若宇, 曹炳阳, 云和明, 陈宝明
 Study on NonNewtonian Behaviors of LennardJones Fluids via Molecular Dynamics Simulations
 LJ流体非牛顿现象的分子动力学模拟
 Chinese Journal of Chemical Physics, 2016, 29(6): 754760
 化学物理学报, 2016, 29(6): 754760
 http://dx.doi.org/10.1063/16740068/29/cjcp1606129

Article history
 Received on: June 20, 2016
 Accepted on: August 11, 2016
b. Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China
Complex liquids like colloidal solutions and polymer melts can exhibit distinct nonNewtonian behaviors, namely the change of viscosity and other physical properties with shear rate [13]. This is the common sense because the rheological properties can be observed, measured and explained in the laboratory. In the meantime, with the wide application of computer simulations, people realize that simple fluids like argon at high shear rates can also shear thin, shear thicken or display viscoelasticity. Usually, the Weissenberg number (
Molecular dynamics (MD) serves as an efficient tool to generate fluid flows and enables analyses on the related phenomena [7, 8]. Some studies were already carried out to generalize the rheological laws of simple fluids under various conditions with MD simulations. Erpenbeck [9] reported a disorderorder transition with a hardsphere model, where the alignment of particles occurs along the flow direction, called the "string phase". This phenomenon was further explained as the occurrence of hydrodynamic instability [10]. Unfortunately, researchers found that this "string phase" is actually an unphysical phenomenon originated from the improper thermostat applied [11, 12], and then the configurational thermostat [13, 14] was developed to deal with this issue. In MD simulations, Delhommelle et al. [15] applied the configurational thermostat to both steady and oscillatory states under extensive shear rates and frequencies, and the existence of shearinduced ordering under oscillatory shear was discovered. One of the famnewpage theoretical analyses for nonNewtonian behavior of simple fluids was conducted by Gan and Eu [16, 17]. They derived a heierarchy of nonlinear integral equations for the nonequilibrium fluctuations based on the generalized Boltzmann equation. The shear effects appear in the nonequilibrum potential of fluids at the equilibrium states. Using MD simulation methods, the GE theory was tested to predict the distortion of fluid structure under steady shear [18, 19] and only a qualitative agreement was observed.
From the above brief introduction, the nonNewtonian phenomenon of simple fluids is indeed a very popular research area. Many open questions still remain, especially about the mechanisms of the interesting similarities of rheological behaviors between simple and complex fluids. In this work, we use nonequilibrium molecular dynamics simulations to investigate the rheological behaviors of monoatomic fluids governed by the LennardJones (LJ) potential for steady Couette flows and oscillatory flows. For the steady Couette flow, the shear thinning behavior is observed, which is quite similar to several previous literature reports [13, 15]. We specifically calculate the potential energy for different shear rates and find it has a close relationship with the pair distribution function. For oscillatory flows, we observe the change of the phase difference between shear stress and shear rate with the frequency, which was seldom studied in the literature. Also, the proportion of elastic and viscous modulus, corresponding to the imaginary and real parts of viscosity respectively, varies with the frequency. These findings can also serve as a guide to understand the viscoelasticity of LennardJones fluids.
Ⅱ. SIMULATION METHOD AND DETAILSOur simulation is carried out on a monoatomic fluid governed by the LennardJones (LJ) potential in the form of
$\begin{eqnarray} \varphi (r) = 4\varepsilon \left[{\left( {\frac{\sigma }{r}} \right)^{12}\left( {\frac{\sigma }{r}} \right)^6 } \right] \end{eqnarray}$  (1) 
Here r is the intermolecular distance,
${\dot r_i}{\rm{ = }}\frac{{{p_i}}}{m}$  (2) 
${\dot p_i} = {F_i}  \xi {p_i}$  (3) 
$\begin{eqnarray} \dot \xi = \frac{1}{Q}\left( {\sum {\frac{{p^2 }}{m}}  gk_{\rm{B}} T} \right) \end{eqnarray}$  (4) 
Here subscript i refers to y, z directions, r, m, p, F are the distance, mass, momentum, and force respectively,
Figure 1 exhibits the schematic diagram of our simulation method. At every time step, we select the atoms in slab 1 to give them an additional positive velocity
$\begin{eqnarray} \tau _{xy} = \frac{{\displaystyle\sum\limits_t {\displaystyle\sum\limits_i {v_x } } }}{{2L_x L_y t}} \end{eqnarray}$  (5) 
where i denotes the atoms in slab 1 with an additional
${v_x} = {v_x}_0\cos \left( {\omega t} \right)$  (6) 
Here
In the rest of the paper, the LJ reduced units are introduced to exhibit the simulation results indicated by the superscript "
First, we plot some typical velocity and temperature profiles under different shear rates in Fig. 2. In Fig. 2(a), the velocity profiles are all Vshaped, which are the results of the periodic boundary condition and the current method of generating linear shear flows, i.e. positive and negative additional velocities imposed on the atoms in slab 1 and slab n/2+1 (n=20). Then, in Fig. 2(b), when the shear rates are relatively low (
The shear viscosity as a function of the shear rate is depicted in Fig. 3. The viscosity is evaluated as the ratio of shear stress to shear rate
$\begin{eqnarray} \eta =  \left\langle {\tau _{xy} /\dot \gamma } \right\rangle \end{eqnarray}$  (7) 
From Fig. 3(a), we can easily see a Newtonian plateau before
$\begin{eqnarray} \eta \left( {\dot \gamma } \right) = \frac{{\eta _0 }}{{\left[{1 + \left( {\lambda \dot \gamma } \right)^2 } \right]^q }} \end{eqnarray}$  (8) 
which was developed empirically from the behavior of complex fluids. In Eq.(8),
To discuss the rheological behavior, we need to take a look at the stress tensor T of the fluids. It can be written as
$\begin{eqnarray} {\bf{T}} = \left( {\begin{array}{*{20}c} {t_{xx} }&{t_{xy} }&{t_{xz} } \\ {t_{yx} }&{t_{yy} }&{t_{yz} } \\ {t_{zx} }&{t_{zy} }&{t_{zz} } \\ \end{array}} \right) = \frac{1}{3}\left[\textrm{tr} ({{\bf{T}}}) \right]{\bf{I}} + {\boldsymbol{\sigma }} \end{eqnarray}$  (9) 
where tr (T) is the trace of the stress tensor, I is the unit tensor and
$\begin{eqnarray} p =  \frac{1}{3}[\textrm{tr}({\bf{T}})] \end{eqnarray}$  (10) 
and the stress tensor is rewritten as
$\begin{eqnarray} {\bf{T}} =  p{\bf{I}} + {\boldsymbol{\sigma }} \end{eqnarray}$  (11) 
or equivalently,
$\begin{eqnarray} T_{ij} =  p\delta _{ij} + \sigma _{ij} \end{eqnarray}$  (12) 
When a simple Couette flow is constructed in the system, for purely viscous fluids, Eq.(9) is in the form of
$\begin{eqnarray} {\bf{T}} =  p\left( {\begin{array}{*{20}c} 1&0&0 \\ 0&1&0 \\ 0&0&1 \\ \end{array}} \right) + \left( {\begin{array}{*{20}c} 0&\tau&0 \\ \tau&0&0 \\ 0&0&0 \\ \end{array}} \right) \end{eqnarray}$  (13) 
However, for the viscoelastic fluids
$\begin{eqnarray} {\bf{T}} =  p\left( {\begin{array}{*{20}c} 1&0&0 \\ 0&1&0 \\ 0&0&1 \\ \end{array}} \right) + \left( {\begin{array}{*{20}c} {\sigma _{xx} }&\tau&0 \\ \tau&{\sigma _{yy} }&0 \\ 0&0&{\sigma _{zz} } \\ \end{array}} \right) \end{eqnarray}$  (14) 
Actually, Eq.(9) and Eq.(10) give just one possible way of defining the pressure which is not necessarily the real hydrostatic pressure of the fluids. This means that the absolute values of
To calculate normal stress differences along with isotropic pressure, we first distinguish n planes perpendicular to x direction. Then for every time step, the
$\begin{eqnarray} t_{xx} = \frac{{J_x }}{{tL_y L_z }} \end{eqnarray}$  (15) 
where
Figure 4 shows the dependences of the normal stress differences and isotropic pressure on the 3/2 power of shear rate. At low shear rates, s_{1} and s_{2} are both close to zero, which indicates the equilibrium state of the fluid. While with the increase of the shear rate, s_{1} increases greatly to large positive values and s_{2} has a slight decrease, leading to negative values. And the absolute values of s_{1} are greater than that of s_{2} at all time. Also, linear variance can be observed which agrees with the theory in Ref.[28]. In terms of the macroscopic hydrodynamics of nonNewtonian fluids, the normal stress effects discussed above are responsible for lots of interesting phenomena, like Weissenberg effect, of liquid polymers.
To explain the nonNewtonian behaviors under steady shear, we use the pair radial distribution function
$\begin{eqnarray} g\left( r \right) = \rho ^{  2} \left\langle {\sum\limits_i {\sum\limits_{j \ne i} {\delta \left( {{\bf{r}}_i } \right)\delta \left( {{\bf{r}}_j  {\bf{r}}} \right)} } } \right\rangle\nonumber\\ = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j \ne i} {\delta \left( {{\bf{r}}  {\bf{r}}_{ij} } \right)} } } \right\rangle \end{eqnarray}$  (16) 
Figure 5 shows
As shown in Fig. 1, we choose a fixed value of the velocity amplitude
$\begin{eqnarray} \tau = \tau _0 \cos (\omega t) \end{eqnarray}$  (17) 
which can be rewritten in the complex form
$\begin{eqnarray} \tau = \tau _0 \textrm{e}^{i\omega t} \end{eqnarray}$  (18) 
Correspondingly, the response of the velocity field at the boundary gives the shear rate
$\begin{eqnarray} \dot \gamma = \dot \gamma _0 \textrm{e}^{i(\omega t + 180^\circ + \phi )} \end{eqnarray}$  (19) 
Then, the shear viscosity is in the form of
$\begin{eqnarray} \eta =  \frac{\tau }{{\dot \gamma }} = \frac{{\tau _0 }}{{\dot \gamma _0 }}\cos \phi  i\frac{{\tau _0 }}{{\dot \gamma _0 }}\sin \phi \end{eqnarray}$  (20) 
And the real and imaginary parts can be defined, respectively, as
$\begin{eqnarray} \eta _\textrm{r} = \frac{{\tau _0 }}{{\dot \gamma _0 }}\cos \phi \end{eqnarray}$  (21) 
$\begin{eqnarray} \eta _\textrm{i} = \frac{{\tau _0 }}{{\dot \gamma _0 }}\sin \phi \end{eqnarray}$  (22) 
Finally, we can obtain the shear strain and modulus
$\begin{eqnarray} \gamma = \gamma _0 \textrm{e}^{i(\omega t + 90^\circ + \phi )} = \frac{1}{\omega }\dot \gamma _0 \textrm{e}^{i(\omega t + 90^\circ + \phi )} \end{eqnarray}$  (23) 
$\begin{eqnarray} G = \frac{\tau }{\gamma }=\frac{{\omega \tau _0 }}{{\dot \gamma _0 }}\sin \phi i\frac{{\omega \tau _0 }}{{\dot \gamma _0 }}\cos \phi \end{eqnarray}$  (24) 
where the elastic modulus and viscous modulus are
$\begin{eqnarray} G' = \frac{{\omega \tau _0 }}{{\dot \gamma _0 }}\sin \phi = \omega \eta _\textrm{i} \end{eqnarray}$  (25) 
$\begin{eqnarray} G" = \frac{{\omega \tau _0 }}{{\dot \gamma _0 }}\cos \phi = \omega \eta _\textrm{r} \end{eqnarray}$  (26) 
We assume that the fluid exhibits linear viscoelastic response to the imposed shear stress, which is a sinusoidal function of time. This means the shear rate amplitude should be small enough to avoid the nonlinear effects. Then, under this circumstance, Eqs.(17)(26) are applicable. Otherwise, to derive the complex viscosity, the calculation of a Fourier transform of the equilibrium autocorrelation function of the shear stress will be needed [15].
In Fig. 7, we plot the velocity profiles along the x direction for three different angular frequencies. The velocity data points are averaged from different time steps t=
The time dependent shear stress and shear rate under
Figure 10 shows the real and imaginary parts of complex viscosity. The real part first increases with the frequency and then decreases. On the other hand, the imaginary part tends to increase monotonically with the frequency. We can conclude that at relatively high frequency, the imaginary part holds a bigger proportion than the real part, which denotes a more important role of elasticity. In fact, from Eq.(25) and Eq.(26) the imaginary part of viscosity has similar physical meaning as the elastic modulus and they both indicate the strength of the elasticity. The elastic modulus and viscous modulus are presented in Fig. 11. At low frequency, the viscous modulus is larger than the elastic one, while the case reverses at high frequencies. All these findings manifest the same point that under oscillatory shear the LJ fluid does exhibit viscoelasticity and at high frequencies the elasticity dominates over viscosity.
Ⅳ. CONCLUSIONThe nonNewtonian behaviors of LennardJones fluids undergoing steady and oscillatory shear flows are studied. For steady Couette flow, we observe that the shear viscosity has a plateau at low shear rates known as the Newtonian zone and then experiences a sharp drop at high shear rates. This phenomenon is widely known as shear thinning. And the normal stress differences are also calculated. The results reveal that the first normal stress difference increases from zero at equilibrium state to relatively large positive values with increasing shear rates, while the second normal stress difference decreases to negative values. This effect illustrates the viscoelastic nature of the simulated LJ monatomic fluids. We have analyzed the radial distribution function and potential energy to illustrate that the repulsion between molecules is more powerful while the attraction is weaker when the shear rate is higher. This may be interpreted as it is the change of microscopic arrangements of the molecules that leads to these nonNewtonian phenomena. In the study of the oscillatory shear flow, we demonstrate that LennardJones fluids are actually viscoelastic fluids and the elasticity dominates over viscosity at relatively high frequencies. The elaboration is carried out through the following perspectives: The phase difference between shear stress and shear rate increases with the frequency; the imaginary part of viscosity dominates over the real part at high frequencies; the elastic modulus dominates over viscous modulus at high frequencies. These three findings are in great coincidence with each other and support the idea of viscoelasticity in LennardJones fluids. There are actually no onehundredpercent viscous fluids and the elasticity should be put forward to explain the nonNewtonian phenomena. From our simulations, it can be summarized that LennardJones fluids do have very similar rheological behaviors to complex fluids.
Ⅴ. ACKNOWLEDGMENTSThis work is supported by the National Natural Science Foundation of China (No.51322603, No.51136001, No.51356001), the Science Fund for Creative Research Groups (No.51321002), the Program for New Century Excellent Talents in University, the Tsinghua National Laboratory for Information Science and Technology of China (TNList), the Foundation of Key Laboratoty of Renewable Energy Utilization Technologies in Buildings of the National Education Ministry in Shandong Jianzhu University (No.KF201301).
[1]  H. L. Yang, J. M. Ruan, J. P. Zou, Q. M. Wu, Z. C. Zhou, and Y. Y. Xie, Chin. J. Chem. Phys. 22 , 46 (2009). DOI:10.1088/16740068/22/01/4650 
[2]  Q. Wang, H. Z. Li, Y. J. Xie, H. Y. Li, and H. Y. Yang, Chin. J. Chem. Phys. 25 , 448 (2012). DOI:10.1088/16740068/25/04/448456 
[3]  V. Citro, F. Giannetti, and J. O. Pralits, Fluid Dyn. Res. 47 , 015503 (2015). DOI:10.1088/01695983/47/1/015503 
[4]  D. M. Heyes, J. NonNewtonian Fluid Mech. 21 , 137 (1986). DOI:10.1016/03770257(86)800325 
[5]  R. Laghaei, A. E. Nasrabad, and B. C. Eu, J. Chem. Phys. 123 , 234507 (2005). DOI:10.1063/1.2138700 
[6]  D. J. Evans, Phys. Rev. A 23 , 1988 (1981). DOI:10.1103/PhysRevA.23.1988 
[7]  R. F. Cracknell, D. Nicholson, and N. Quirke, Phys. Rev. Lett. 74 , 2463 (1995). DOI:10.1103/PhysRevLett.74.2463 
[8]  X. B. Nie, S. Y. Chen, and M. O. Robbins, J. Fluid Mech. 500 , 55 (2004). DOI:10.1017/S0022112003007225 
[9]  J. J. Erpenbeck, Phys. Rev. Lett. 52 , 1333 (1984). DOI:10.1103/PhysRevLett.52.1333 
[10]  P. F. Lutsko, and J. W. Dufty, Phys. Rev. Lett. 57 , 2775 (1986). DOI:10.1103/PhysRevLett.57.2775 
[11]  D. J. Evans, and G. P. Morriss, Phys. Rev. Lett. 56 , 2172 (1986). DOI:10.1103/PhysRevLett.56.2172 
[12]  D. J. Evans, S. T. Cui, H. J. M. Hanley, and G. C. Straty, Phys. Rev. A 46 , 6731 (1992). DOI:10.1103/PhysRevA.46.6731 
[13]  O. G. Jepps, G. Ayton, and D. J. Evans, Phys. Rev. E 62 , 4757 (2000). DOI:10.1103/PhysRevE.62.4757 
[14]  L. Lue, O. G. Jepps, J. Delhommelle, and D. J. Evans, Mol. Phys. 100 , 2387 (2002). DOI:10.1080/00268970210122145 
[15]  J. Delhommelle, J. Petravic, and D. J. Evans, J. Chem. Phys. 120 , 6117 (2004). DOI:10.1063/1.1652014 
[16]  H. H. Gan, and B. C. Eu, Phys. Rev. A 45 , 3670 (1992). DOI:10.1103/PhysRevA.45.3670 
[17]  H. H. Gan, and B. C. Eu, Phys. Rev. A 46 , 6344 (1992). DOI:10.1103/PhysRevA.46.6344 
[18]  Kalyuzhnyi. V. Yu, S. T. Cui, P. T. Cummings, and H. D. Cochran, Phys. Rev. E 60 , 1716 (1999). 
[19]  Kalyuzhnyi. V. Yu, S. T. Cui, and H. D. Cochran, Phys. Rev. E 63 , 011209 (2000). DOI:10.1103/PhysRevE.63.011209 
[20]  G. S. Grest, and K. Kremer, Phys. Rev. A 33 , 3628 (1986). DOI:10.1103/PhysRevA.33.3628 
[21]  B. Y. Cao, M. Chen, and Z. Y. Guo, Appl. Phys. Lett. 86 , 091905 (2005). DOI:10.1063/1.1871363 
[22]  B. Y. Cao, M. Chen, and Z. Y. Guo, Phys. Rev. E 74 , 066311 (2006). DOI:10.1103/PhysRevE.74.066311 
[23]  M. P. Allen, and D. J. Tildesley, Computer Simulation of Liquids. New York: Oxford University (1989). 
[24]  H. Eslami, F. Mller, and  Plathe, J. Phys. Chem. B 114 , 387 (2010). 
[25]  I. Borzs, P. T. Cummings, and D. J. Evens, Mol. Phys. 100 , 2735 (2002). DOI:10.1080/00268970210137275 
[26]  P. J. Carreau, PhD Thesis, Madison:University of Wisconsin (1968). 
[27]  C. M. Tenney, and E. J. Maginn, J. Chem. Phys. 132 , 014103 (2010). DOI:10.1063/1.3276454 
[28]  K. Kawasaki, and J. D. Gunton, Phys. Rev. A 8 , 2048 (1973). DOI:10.1103/PhysRevA.8.2048 
b. 清华大学航天航空学院, 热科学与动力工程教育部重点实验室, 北京 100084