Chinese Journal of Chemical Physics  2016, Vol. 29 Issue (6): 754-760

#### The article information

Ruo-yu Dong, Bing-yang Cao, He-ming Yun, Bao-ming Chen

Study on Non-Newtonian Behaviors of Lennard-Jones Fluids via Molecular Dynamics Simulations
LJ流体非牛顿现象的分子动力学模拟
Chinese Journal of Chemical Physics, 2016, 29(6): 754-760

http://dx.doi.org/10.1063/1674-0068/29/cjcp1606129

### Article history

Accepted on: August 11, 2016
Study on Non-Newtonian Behaviors of Lennard-Jones Fluids via Molecular Dynamics Simulations
Ruo-yu Donga,b, Bing-yang Caob, He-ming Yuna, Bao-ming Chena
Dated: Received on June 20, 2016; Accepted on August 11, 2016
a. Key Laboratory of Renewable Energy Utilization Technologies in Building of Ministry of Education, Shandong Jianzhu University, Jinan 250101, China;
b. Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China
*Author to whom correspondence should be addressed. Bing-yang Cao,E-mail:caoby@tsinghua.edu.cn,Tel./FAX:+86-10-6279-4531
Abstract: Using nonequilibrium molecular dynamics simulations, we study the non-Newtonian rheological behaviors of a monoatomic fluid governed by the Lennard-Jones potential. Both steady Couette and oscillatory shear flows are investigated. Shear thinning and normal stress effects are observed in the steady Couette flow simulations. The radial distribution function is calculated at different shear rates to exhibit the change of the microscopic structure of molecules due to shear. We observe that for a larger shear rate the repulsion between molecules is more powerful while the attraction is weaker, and the above phenomena can also be confirmed by the analyses of the potential energy. By applying an oscillatory shear to the system, several findings are worth mentioning here:First, the phase difference between the shear stress and shear rate increases with the frequency. Second, the real part of complex viscosity first increases and then decreases while the imaginary part tends to increase monotonically, which results in the increase of the proportion of the imaginary part to the real part with the increasing frequency. Third, the ratio of the elastic modulus to the viscous modulus also increases with the frequency. These phenomena all indicate the appearance of viscoelasticity and the domination of elasticity over viscosity at high oscillation frequency for Lennard-Jones fluids.
Key words: Non-Newtonian    Viscoelasticity    Lennard-Jones fluids    Molecular dynamics
Ⅰ. INTRODUCTION

Complex liquids like colloidal solutions and polymer melts can exhibit distinct non-Newtonian behaviors, namely the change of viscosity and other physical properties with shear rate [1-3]. 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 ($W_{\rm{i}}$), which is the structural relaxation time and shear rate, can be used to relate rheological behaviors of simple liquids to complex liquids [4]. It's believed that under the same $W_{\rm{i}}$, fluids show similar flow behaviors. This was also illustrated through the universal scaling obtained from the statistical-mechanical analysis [5]. For the complex liquids in industry or labs, the value of $W_{\rm{i}}$ is of order unity [4]. Considering the relaxation time of simple fluids is of order picosecond, the very much high shear rates of THz are unrealizable experimentally because the macroscopic turbulence prevents the attachment of high shear rates [6]. That's why people still do not totally accept the standpoint that all fluids can be expected to be non-Newtonian.

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 disorder-order transition with a hard-sphere 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 shear-induced ordering under oscillatory shear was discovered. One of the fam-newpage theoretical analyses for non-Newtonian 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 non-Newtonian 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 Lennard-Jones (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 Lennard-Jones fluids.

Ⅱ. SIMULATION METHOD AND DETAILS

Our simulation is carried out on a monoatomic fluid governed by the Lennard-Jones (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, $\sigma$ is the molecular diameter of 3.405×10-10 m and $\varepsilon$ is the potential well depth of 1.67×10-21 J. Orthorhombic cells are used in the system with periodic boundary conditions (PBC) applied in three directions. The systems all consist of N=1000 particles and are divided into n slabs (labeled from 1 to n) along the z direction. A constant-T simulation also known as the canonical ensemble (NVT) with a Nose-Hoover thermostat is used. The thermostat is applied only in two directions (y, z) except the flow direction (x) [20-22] and the equations of motion are

 ${\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, $k_{\rm{B}}$ is the Boltzmann constant, g=$2N$, Q=$2Nk_{\rm{B}}T\tau_{\rm{r}}$2 and $\tau_{\rm{r}}$=0.05 is the relaxation time. Applying the thermostat in only two directions will guarantee that the flows are not affected by the thermostat. To reduce the time-consuming calculations of the inter-particle interactions, a cutoff distance of $r_{\rm{cut}}$=2.0$\sigma$ and a cell-linked list method [23] is applied in the simulations. A leap-frog Verlet algorithm with a time step of dt=0.005$\tau$ is used to integrate the equations of motion.

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 $v_x$, while the atoms in slab n/2+1 are with a negative velocity $-v_x$. By using the periodic boundary conditions, Couette flow can be established at both halves of the system. The numbers of atoms from the two slabs involved in this process are selected to be the same. Then after time t, the momentum flux, i.e. shear stress, can be expressed as

 $\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$. The factor 2 arises from the periodicity of the system. Besides the common steady flow, we are also concerned with oscillatory shear flow and the imposed velocity is in the form of

 ${v_x} = {v_x}_0\cos \left( {\omega t} \right)$ (6)
 FIG. 1 Schematic diagram of the simulation method.

Here $v_x_0$ is the amplitude of the oscillatory velocity and $\omega$ is the frequency.

In the rest of the paper, the LJ reduced units are introduced to exhibit the simulation results indicated by the superscript "$^*$", e.g. the temperature $T^*$=$Tk_{\rm{B}}/\varepsilon$, number density $\rho^*$=$\rho\sigma^3$, velocity $v^*$=$v\sqrt{m/\varepsilon}$, time $t^*$=$t\sqrt{\varepsilon/m}/\sigma$, diameter $r^*$=$r/\sigma$, shear viscosity $\eta^*$=$\eta\sigma^2/\sqrt{\varepsilon m}$, shear rate $\dot \gamma^*$=$\dot \gamma\sigma\sqrt{m/\varepsilon}$, shear stress $s^*$=$s\sigma^3/\varepsilon$ and angular frequency $\omega^*$=$\omega\sigma\sqrt{m/\varepsilon}$. The state points of $T^*$=0.722 and $\rho^*$=0.8442 near the triple point of the LJ fluid are selected as the case study in the following work.

Ⅲ. SIMULATION RESULTS AND DISCUSSION A. Steady shear flow

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 V-shaped, 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 ($\dot\gamma^*$=0.006, 0.12), the temperature profiles are almost horizontal straight lines with a value very close to 0.722, which indicates a negligible viscous heating behavior. While at high shear rate ($\dot\gamma^*$=1.61), the profile is composed of two distinct parabolic shapes, interconnected at the center of the system. Parabolic temperature profiles for fluids under linear shear have also been observed in some previous study [24]. The non-negligible viscous heating causes the remarkable temperature differences, and the two connected parabolic shapes correspond with the V-shaped velocity profiles. More interestingly, temperature jumps can be observed at the regions of slab 1 and slab n/2+1 (n=20).

 FIG. 2 (a) Velocity profiles and (b) temperature profiles of different shear rates.

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)
 FIG. 3 (a) Dependence of the shear viscosity on the shear rate, a Carreau model; (b) dependence of the shear viscosity on the 0.5 power of shear rate with the linear fit (intercept: 3.41, slope: -1.33).

From Fig. 3(a), we can easily see a Newtonian plateau before $\dot\gamma^*$=0.01 where the viscosity is independent of the shear rate, which agrees with the reference data reported by Borzsák et al. [25]. From around $\dot\gamma^*$=0.1, the shear viscosity decreases markedly with the increase of the shear rate, which indicates the appearance of shear thinning. A comparison is made between our results and the famous Carreau model [26],

 $\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), $\eta_0$ is the Newtonian viscosity and the parameter $\lambda$ is a characteristic time constant, approximately equal to the inverse of the shear rate when shear thinning behavior begins [27]. The fitting parameters are estimated to be $\eta_0$=3.31, $\lambda$=7.08 and q=0.11, consistent with the results in Ref.[27]. Also, just as the mode coupling theory stated, we find in Fig. 3(b) the shear viscosity does vary linearly with the square root of the shear rate. In fact, this shear thinning phenomenon is not only common to complex fluids, but also can be observed in the rheology of Lennard-Jones fluids.

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 ${\boldsymbol{\sigma}}$ is the Deviatoric stress tensor. We can define the isotropic pressure as

 $\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 $\sigma_{ii}$ may vary, but the differences between them do not. Then the first and second normal stress differences are defined as s1=$\sigma_{xx}$-$\sigma_{yy}$ and s2=$\sigma_{yy}$-$\sigma_{zz}$. In this case, three variables are needed to describe the flow characteristics of viscoelastic fluids, namely shear stress $\tau$, s1, and s2.

To calculate normal stress differences along with isotropic pressure, we first distinguish n planes perpendicular to x direction. Then for every time step, the $v_x$ of atoms moving across plane i and the $\Delta v_x$ generated by pairs of atoms interacting across plane i are recorded as the momentum flowing past this plane. We can derive the total momentum flux i.e. normal stress $t_{xx}$ of plane i:

 $\begin{eqnarray} t_{xx} = \frac{{J_x }}{{tL_y L_z }} \end{eqnarray}$ (15)

where $J_x$ is the total x component of momentum passing through Plane i over time. Finally, for n planes, the average of $t_{xx}$ can be calculated, which corresponds to the normal stress component of the x direction of the stress tensor T. Apparently, $t_{yy}$ and $t_{zz}$ can be derived in the same way, and following Eqs.(9)-(14), we can obtain the isotropic pressure and normal stress differences.

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, s1 and s2 are both close to zero, which indicates the equilibrium state of the fluid. While with the increase of the shear rate, s1 increases greatly to large positive values and s2 has a slight decrease, leading to negative values. And the absolute values of s1 are greater than that of s2 at all time. Also, linear variance can be observed which agrees with the theory in Ref.[28]. In terms of the macroscopic hydrodynamics of non-Newtonian fluids, the normal stress effects discussed above are responsible for lots of interesting phenomena, like Weissenberg effect, of liquid polymers.

 FIG. 4 Dependences of the normal stress differences and isotropic pressure on the 1.5 power of shear rate with the linear fit. (For s1, intercept: 0.31, slope: 7.31; for s2, intercept: -0.14, slope: -0.4; for p, intercept: 2.46, slope: 4.07).

To explain the non-Newtonian behaviors under steady shear, we use the pair radial distribution function $g(r)$ as the monitor of the shear induced structural changes. It gives the probability of finding a pair of atoms a distance r apart, relative to the probability of a random distribution at the same density [23]. To directly compute $g(r)$, we can use the following equation [23]

 $\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 $g(r)$ in two shear rates. The first peak of $g(r)$ arises just around $r/\sigma$=1.0, where the attraction and repulsion have a balance. We find that $g(r)$ is larger for a stronger shear flow when $r/\sigma <$1.0, i.e. at the ascent stage of the first peak, and this indicates the molecules are getting closer to each other for a higher shear rate. Oppositely, $g(r)$ becomes smaller for a stronger shear at the descent stage of the first peak. The above observation may suggest that for a larger shear rate the repulsion is more powerful while the attraction is weaker. Through the calculation of the repulsive and attractive parts of the potential energy, we confirm this conjecture: In Fig. 6, the attractive $U_{\rm{a}}$, repulsive $U_{\rm{r}}$ and the total potential energies U under different shear rates are plotted. The absolute values of $U_{\rm{a}}$ and U decrease monotonically with the shear rate, while the trend for $U_{\rm{r}}$ reverses. Apparently, the attractive potential energy dominates over the repulsive one and the trends agree with the rule obtained from the analyses of the radial distribution function.

 FIG. 5 Pair radial distribution function for two shear rates.
 FIG. 6 Dependence of the total (U), attractive ($U_{\rm{a}}$) and repulsive ($U_{\rm{r}}$) parts of the potential energy on the shear rate.
B. Oscillatory shear flow

As shown in Fig. 1, we choose a fixed value of the velocity amplitude $v_{x0}$ and by controlling the number of atoms applied with $\pm v_{x0}$ in slab 1 and slab n/2+1, the stress amplitude is also fixed. We gather the statistics from slab 1 under different angular frequencies. The shear stress at slab 1 can be expressed as

 $\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=$T_{\rm{p}}/4$+$kT_{\rm{p}}$ (k=0, 1, 2$\ldots$) at the same phase, where $T_{\rm{p}}$ stands for the period. It should be noted that at a relatively small frequency ($\omega^*$=0.01$\pi$) the profile is almost linear, while at large frequencies ($\omega^*$=0.25$\pi$, $\pi$) the profiles exhibit distinct nonlinear features. For the latter case, even some approximate still flow fields in the middle can be observed. The above phenomena are mainly due to the fact that the oscillation is only applied at the boundary and at the center of the system. When the angular frequency is too large, the oscillation does not have enough time to transfer deep into the system under the present selected time step, and thus some intermediate regions are not affected.

 FIG. 7 Instantaneous velocity profiles along the x direction for different angular frequencies obtained from the average at t=$T_{\rm{p}}$/4+$kT_{\rm{p}}$ (k=0, 1, 2, $\ldots$).

The time dependent shear stress and shear rate under $\omega^*$=5$\pi$ are shown in Fig. 8. Two periods are selected to exhibit the results. Both the shear stress and shear rate are sinusoidal functions of time, and a phase difference can be seen between them. Figure 9 shows the phase differences between the shear stress and the inverse of shear rate under a wide range of angular frequencies. A monotonic increase with the increase of frequency can be observed. For the viscous fluid, the phase difference between shear stress and rate should equal zero. And for the elastic solid it should equal 90°, which means the shear stress and shear strain are in phase. Our results suggest that the elasticity dominates at high frequency even for Lennard-Jones fluids.

 FIG. 8 Time dependent shear stress and shear rate under $\omega^*$=5$\pi$.
 FIG. 9 Dependence of the phase difference between the shear stress and the inverse of shear rate on frequency.

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.

 FIG. 10 Dependence of the real ($\eta_{\rm{r}}$) and imaginary ($\eta_{\rm{i}}$) parts of the complex viscosity on frequency.
 FIG. 11 Dependence of the elastic modulus ($G'$), and viscous modulus ($G"$) on frequency.
Ⅳ. CONCLUSION

The non-Newtonian behaviors of Lennard-Jones 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 non-Newtonian phenomena. In the study of the oscillatory shear flow, we demonstrate that Lennard-Jones 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 Lennard-Jones fluids. There are actually no one-hundred-percent viscous fluids and the elasticity should be put forward to explain the non-Newtonian phenomena. From our simulations, it can be summarized that Lennard-Jones fluids do have very similar rheological behaviors to complex fluids.

Ⅴ. ACKNOWLEDGMENTS

This 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).

Reference
 [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/1674-0068/22/01/46-50 [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/1674-0068/25/04/448-456 [3] V. Citro, F. Giannetti, and J. O. Pralits, Fluid Dyn. Res. 47 , 015503 (2015). DOI:10.1088/0169-5983/47/1/015503 [4] D. M. Heyes, J. Non-Newtonian Fluid Mech. 21 , 137 (1986). DOI:10.1016/0377-0257(86)80032-5 [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
LJ流体非牛顿现象的分子动力学模拟

a. 山东建筑大学, 可再生能源建筑利用技术教育部重点实验室, 济南 250101;
b. 清华大学航天航空学院, 热科学与动力工程教育部重点实验室, 北京 100084