CC BY
22
Molecular dynamics simulation of fluid viscosity in nanochannels
V. Rudyak, A. Belkin
Novosibirsk State University of Architecture and Civil Engineering, Leningradskaya 113, Novosibirsk, 630008, Russia
valery.rudyak@mail.ru, a_belkin@ngs.ru PACS 47.61.-k, 47.10.-g DOI 10.17586/2220-8054-2018-9-3-349-355
The viscosity of fluids in a plane nanochannel has been studied by molecular dynamics method. The effective viscosity coefficient was determined using the fluctuation-dissipation theorem derived previously by the authors from the nonequilibrium statistical theory of fluid transport in confined conditions. It has been found that the fluid viscosity in a nanochannel is strongly dependent on the interaction potential between the fluid and channel wall molecules. In particular, increasing the depth of the potential well of this interaction leads to an increase in the viscosity. At the same time, if the depth of the potential well is small, the fluid viscosity in a nanochannel may be even lower than its viscosity in an unconfined (bulk) system. Thus, the fluid viscosity in a nanochannel and hence the channel flow resistance can be varied by changing the material of the nanochannel walls.
Keywords: nanochannel, confined system, viscosity, molecular dynamics, correlation functions, nonequilibrium statistical mechanics.
Received: 11 April 2018
1. Introduction
Interest in fluid-based transport processes (i.e. liquids and gases) in nanochannels is due to the rapid development of nanotechnologies related to power systems, medicine, the development of advanced filter elements, etc. At present, it is obvious that the properties of fluids in nanochannels and nanometer-sized pores differ fundamentally from their bulk properties. In these channels, the fluid density becomes nonuniform [1], the molecular diffusion is anisotropic [2], and other transport properties should also change. However, so far, there is no consensus about what relations should be used to describe the fluid transport coefficients, in particular viscosity, under confined conditions.
Simple models have been developed [3,4] which consider a fluid layer with altered properties, in particular viscosity, near the channel walls. However, these models do not provide an answer to the question of how to determine the size of this layer and how the viscosity in it changes. The need to use constants that cannot be determined within the proposed theoretical model is a serious problem (see, e.g., [5, 6]). Of course, the unknown constants can be determined using experimental data; however, the available data for nanochannels are contradictory [7,8]. Moreover, these experimental data are often unreliable as there are no direct methods for measuring physical characteristics in nanochannels and they have to be extracted from certain measured integral characteristics and interpreted using conventional macroscopic theories.
A statistical theory of transport processes under confined conditions was developed in [9,10]. It has been established that the transport properties of fluids in nanochannels and nanopores are determined not only by the characteristics of the fluid, but also by the properties of the entire fluid - channel (pore) wall system. These transport properties are described by effective transport coefficients, which in [9,10] are determined using fluctuation-dissipation relations between the transport coefficients and correlation functions of the corresponding microscopic fluxes. These relations extend the well-known Green-Kubo formula for fluids in bulk.
Analytically, the correlation functions can be calculated only in the simplest case of a rarefied gas (see [11-13]). In all other cases, it is necessary to use methods of molecular modeling. The purpose of the present study was to perform a molecular dynamics (MD) method to calculate the viscosity of fluids in a plane nanochannel (nanoslit) and study the factors that determine the fluid viscosity. We modeled the viscosity of argon and benzene in channels of different heights and simultaneously studied the fluid structure in these channels.
2. Method for determining the effective viscosity coefficient
In [9,10], explicit expressions for the effective viscosity coefficient n of fluids in confined conditions were constructed using methods of nonequilibrium statistical mechanics. This coefficient is expressed in terms of time correlation functions of microscopic stress tensors. Along with the contribution from the stress tensor of the fluid,
it includes cross correlations characteristic of the fluid-wall surface system and due to the interaction between the fluid molecules and the molecules of the channel walls. The viscosity is given by the relation:
n = Vff + Vfb,
(1)
nff
5kTf
Jf (0) : Jf (t)) dt,
1
nfb =
Jb(0) : Jf (t) + Jfb(tn + Jf (0) : Jfb(t)) dt +
1
Jb(0): Jfb(t )) dt, (2)
5kTf J \ " \ 1 w • >) • ^ v 5kTb
0 0
where the integrands in formulas (2) are classical correlation functions [13, 14], the angular brackets denote averaging over the equilibrium distribution function, T is the temperature, k is the Boltzmann constant, and t is the time when the viscosity reaches the so-called plateau value [15] that corresponds to the time of decay of the correlation functions. Here and below, the subscript f corresponds to the fluid molecules, and b to the molecules of the channel walls.
The time correlation functions in formulas (2) are determined by the corresponding microscopic fluxes:
Na
Ja = £pp + - rj) Fij,
1
Na fbb Nv
2
fb
=1 V j= 1
1 Nf Nb -1 ^^ (r¿rj)Fij.
i=1 j=1
Here r and p¿ are the radius vector of the center of mass and momentum of the i-th molecule of phase a and Fj is the intermolecular interaction force.
In this paper, we present an investigation of fluid viscosity in a plane channel. The channel walls were modeled by two square plates, each consisting of two rows of molecules. The size of the plates was chosen so that the results remained unchanged as it increased. The molecules were located at the sites of a face-centered cubic lattice. On the channel boundaries along the walls, periodic boundary conditions were used.
The viscosity coefficient (1), (2) was calculated by the standard molecular dynamics method (see e.g., [16]). This was performed with the original SibMD software package which has previously been used to solve various nanofluid transport problems [1,15,17,18]. The Newton equations were integrated by the Schofield scheme [19]. The integration step is equal to 1.09 fs. Due to the local instability and mixing of phase trajectories of the system in the molecular dynamics calculation [1,3,20,21], the results should be averaged over an ensemble of independently obtained phase trajectories. In the calculations presented in this paper, averaging was carried out over five thousand independent phase trajectories.
All intermolecular interactions were described by the cutoff Lennard-Jones potential:
$(r)
4e
aft
r
12
&a¡3
r
0,
- $0, r < Rc;
r > Rc ,
(3)
where a is the effective diameter, e is the depth of the potential well, RC is the cutoff radius of the potential, and r = |r — Tj | is the distance between the centers of molecules, a, 3 = f, b. The cutoff radius of the potential was set equal to RC = 2.5aff, the potential shift was determined from the condition $(RC) = 0. Interaction parameters between the fluid molecules and the walls were calculated from the interaction constants of individual substances using the following combination relations afb = y/affabb and efb = y/effebb.
The effective viscosity coefficient of a fluid in the nanochannel was compared with the viscosity coefficient of the same fluid in a bulk system, which was also simulated by the molecular dynamics method. In this case, for the simulated cubic cell filled with molecules, we used periodic boundary conditions in all directions. The comparison was carried out at the same pressure and fluid temperature. The pressure in the channel was determined from the force per unit area of the walls, and the pressure in bulk was determined by the virial theorem [14].
T
1
T
T
6
3. Simulation results
In formulas (2), we can distinguish three types of correlation functions ^ which define the effective viscosity coefficient: the kinetic correlation functions, ~ (pixpiy(0) • PixPiy(i)), the potential correlation functions, ^p ~ ((xi — xj)Fijy(0) • (xi — xj)Fijy(t)), and the cross-correlation functions, ~ ((xi — xj)Fijy(0) • pixpiy(t)). The corresponding contributions also appear when calculating the viscosity coefficient (1). The kinetic contributions are due to momentum transfer during molecular motion. It is these contributions that determine the momentum
transfer in a rarefied gas and the viscosity of the gas. The above-mentioned potential contributions are due to momentum transfer in the interaction of molecules.
Figure 1 shows the time dependences of these three components and the total correlation function for the viscosity of argon in the nanochannel considering both fluid - fluid and fluid - wall interactions. The spherically symmetric potential (3) adequately describes inert argon molecules (off = 3.405 A, eff/k = 119.8 K). In an unconfined system, viscosity can be determined with good accuracy by the MD method, at least up to the saturation line [22]. The temperature in our MD simulation was 160 K, and the reduced concentration of molecules no3 = 0.4 (here and below, a = off). All the correlation functions are normalized to (p^Py), so that the initial value of the kinetic component is unity. The density of the system is high enough, and the potential contribution is the major one.
Fig. 1. Viscosity correlation function and its components versus time (in picoseconds) for argon in a nanochannel of height h = 27.2 A (solid line). The dashed line corresponds to the potential component, the dotted line to the kinetic component, and the dash-dotted line to the cross-correlation. The parameters of the potential of the wall molecules correspond to argon
Let us consider in more detail the relaxation of the components of the correlation functions in bulk and in the nanochannel. A comparison of the two cases is shown in Fig. 2 (the data in the figure are given on a logarithmic scale).
Two important points should be noted. First, the values of the potential component in the nanochannel far exceed those in bulk throughout its time evolution. Second, there is anisotropy in momentum transfer during collisions in the nanochannel (see curves 4 and 5 in Fig. 2). The correlation function increases more strongly in the plane perpendicular to the walls of the channel.
In contrast, the kinetic component in the nanochannel is practically isotropic and decays faster than in bulk (see curves 1 and 2 in Fig. 2), which is due, in particular, to the faster velocity relaxation of the molecules in collisions with the channel walls.
The viscosity coefficient (1), (2) is a function of the time. Its evolution for the system described above (Figs. 1 and 2) is shown in Fig. 3. The value of the viscosity coefficient is obtained only when this function reaches a plateau value in about 3 psec. The main contribution to the viscosity coefficient comes from the potential component. The viscosity coefficient in nanochannels is approximately 40 % higher than that in bulk.
To answer the question of how and which properties of the channel walls affect the viscosity, we performed a systematic simulation of fluid viscosity in channels where the interaction parameters abb and ebb of the wall molecules (3) were varied.
It was found that the variation in the effective size of the wall molecules abb had little effect on the effective viscosity in the channel. At the same time, the viscosity changed greatly with a change in ebb. Fig. 4 shows the relative viscosity coefficient r/* = n/no in a channel of height h = 27.2 A normalized by the viscosity in bulk no versus the parameter ebb. Increasing this parameter leads to a rapid increase in the effective viscosity coefficient.
On the logarithmic scale of Fig. 4, this increase is almost linear. For a copper wall (eff /k = 1247 K), and for a zinc wall (eff/k = 1040 K), the viscosity coefficient of argon can be expected to increase by a factor of about four. In contrast, for small values of this parameter, ebb < eff, the viscosity coefficient decreases and becomes even less than no.
0 0.1 0.2 t
Fig. 2. Components of the correlation function of the argon viscosity coefficient versus time (in picoseconds) in bulk (the solid line 1 corresponds to the kinetic component, and the solid line 2 to the potential component) and in a nanochannel of height h = 27.2 A (the dotted line 3 corresponds to the kinetic component, the dashed line 4 to the potential component in the plane perpendicular to the channel walls, and the dash-dotted line 5 to the potential component in the plane of the walls)
*
7
1.4
1
0.6
0.2
I 2 3 t
Fig. 3. Viscosity coefficient (solid line) and its components versus integration time (in picoseconds) for argon in a nanochannel of height h = 27.2 A (the dotted line corresponds to the kinetic component, the dashed line to the potential component, and the dot-dashed line to the cross-correlation component)
This conclusion is supported by the data in Fig. 5, which shows curves of the normalized viscosity coefficient of benzene (abb = 5.27 A, ebb/k = 440 K) at room temperature versus height of channels with different wall materials. The lines in the figure correspond to the approximations of the results by the function r/* = 1 ± B/h (plus for higher line and minus for lower line) in which the constant B depends on the properties of the fluid and the walls. In the nanochannel with walls of carbon, whose molecules have a small value of e (abb = 3.4 A, ebb/k = 28 K), the viscosity of argon indeed decreases.
On the other hand, in the channel with aluminum walls (abb = 2.56 A, ebb/k = 857.6 K), there is an increase in viscosity. It can be argued that it is the great depth of the potential well of aluminum molecules that leads to an appreciable increase in viscosity. In all cases, the differences between the viscosity coefficient in the nanochannel and the no value increase monotonically with decreasing channel height.
Fig. 4. Normalized viscosity coefficient of argon in a nanochannel versus potential well depth £bb(K), h = 27.2 A
Fig. 5. Viscosity coefficient versus channel height (in angstrom units) for the following types of fluid and channel walls: C6H6-C (•) and C6H6-Al (■)
4. Conclusions
Strictly speaking, the only consistent method for nanochannel flow simulation is the molecular dynamics method [23]. However, experience has shown that conventional hydrodynamic approach is usually applicable well beyond its range of validity. Therefore, it is possible to attempt to model flows in nanochannels and small microchannels using conventional hydrodynamic methods (see, e.g., [24-26] and the reference therein). For this, however, it is necessary to have adequate data on fluid viscosity in such channels. The present study has shown that this viscosity is largely determined by the walls of the channel. This is due to the fact that the momentum redistribution process in the system depends significantly on the interaction between the fluid molecules and the channel walls, and it is this process that determines the viscosity.
The active interaction of fluid molecules with the channel walls leads to a significant change in the fluid structure near the wall. Fig. 6 shows the MD simulation data for the radial distribution function g2 of benzene near the wall of a channel with carbon walls. Its comparison with the radial distribution function of benzene in bulk shows not only a several-fold increase in the maximum values of the radial function, but also the occurrence of a quasi-long-range order. This increase in the order of the fluid near the surface will lead to an increase in viscosity compared to its value in bulk.
Si
6-
4
2
0
1
2
3
4 r
Fig. 6. Radial distribution function of benzene molecules in bulk (solid line) and in a channel with h = 40 A (dotted line). Intermolecular distances are in the units of a for benzene (5.27 A). Values are normalized by the average concentration of the molecules
However, the contribution to the viscosity coefficient due to the interactions between the fluid molecules and the walls also plays a key role. This should be taken into account, and it should be realized that the fluid viscosity in nanochannels is not only determined by the fluid properties. In this case, it is necessary to introduce the effective viscosity (and thermal conductivity) of the entire fluid-wall system. This system is a special two-phase medium in which transport processes are in a sense similar to those in two-phase suspensions, where it is also necessary to introduce effective transport coefficients (as in the well-known Einstein's and Maxwell's formulas) which are determined by the interaction between base fluid molecules and nanoparticles [18].
Many experiments have demonstrated a significant reduction in flow resistance in microchannels [26,27]. This is usually associated with the slip effect in such channels [28-31]. However, as indicated above, such a reduction in flow resistance may also be due to a decrease in the effective viscosity of the fluid compared to its viscosity in bulk. In this connection, it should be noted that in this paper, we considered flow in a plane channel, in fact, in a nanoslit. However, in closed nanochannels with a circular or square cross-section, these effects will be more pronounced.
Finally, another important point should be noted. The MD simulation results show that transport processes in nanochannels are inhomogeneous, so that, generally speaking, the fluid viscosity also turns out to be inho-mogeneous. Therefore, a consistent approach to the study of nanofluids' transport properties requires the use of generalized nonlocal transport equations and nonlocal transport constitutive relations. Such approaches are well known and have been successfully used [13,14,32-35].
Acknowledgements
This work was supported by the Russian Science Foundation (Grants # 17-01-00040, # 17-58-45023).
[1] Rudyak V.Ya., Belkin A.A., Egorov V.V., Ivanov D.A. Modeling fluid flows in nanochannels by molecular dynamics method. Nanosystems: Physics, Chemistry, Mathematics, 2011, 2(4), P. 100-112.
[2] Andryushchenko V., Rudyak V. Self-diffusion coefficient of molecular fluid in porous media. Defect and Diffusion Forum, 2011, 312-315, P. 417-422.
[3] Blake T. Slip between a liquid and a solid: D.M. Tolstoi's (1952) theory reconsidered. Colloids and Surfaces, 1990, 47, P. 135-145.
[4] Myers T. Why are slip lengths so large in carbon nanotubes? Microfluidics and Nanofluidics, 2010, 10(5), P. 1141-1145.
[5] Popov I.Yu. Statistical derivation of modified hydrodynamic equations for nanotube flows. Phys. Scr., 2011, 83, P. 045601.
[6] Calabro F., Lee K.P., Mattia D. Modelling flow enhancement in nanochannels: Viscosity and slippage. Applied Mathematics Letters, 2013, 26, P. 991-994.
[7] Judy J., Maynes D., Webb B.W. Characterization of frictional pressure drop for liquid flows through microchannels. Int. J. Heat Mass Transfer, 2002, 45, P. 3477-3489.
[8] Pfahler J., Harley J., Bau H., Zemel J.N. Gas and liquid flow in small channels. Micromech. Sensors, Actuators, Syst., 1991, 32, P. 49-58.
[9] Rudyak V.Ya, Belkin A.A. Fluid viscosity under confined conditions. Doklady Physics, 2014, 59(12), P. 604-606.
References
[10] Rudyak V.Ya, Belkin A.A. Statistical mechanics of transport processes of fluids under confined conditions. Nanosystems: Physics, Chemistry, Mathematics, 2015, 6(3), P. 366-377.
[11] Khonkin A.D. Equations for spatial-temporal and temporal correlation functions and proof of the equivalence of the results of the Chapman-Enskog method and temporal correlation functions method. Theoretical and Mathematical Physics, 1970, 5(1), P. 125-135 (in Russian).
[12] Ernst M.H. Formal theory of transport coefficients to general order in the density. Physica, 1966, 32(2), P. 209-243.
[13] Rudyak V.Ya. Statistical aerohydromechanics of homogeneous and heterogeneous media. Vol. 2. Hydromechanics. NSUACE, Novosibirsk, 2005, 468 p. (In Russian).
[14] Zubarev D.N. Nonequilibrium statistical thermodynamics. Consultants Bureau, New York , 1974, 489 p.
[15] Rudyak V.Ya., Belkin A.A., Ivanov D.A., Egorov V.V. The simulation of transport processes using the method of molecular dynamics. Self-diffusion coefficient. High Temperature, 2008, 46(1), P. 30-39.
[16] Rapaport D.C. The Art of Molecular Dynamics Simulation. Cambridge University Press, 1995, 549 p.
[17] Rudyak V.Ya., Krasnolutskii S.L., Ivanov D.A. Molecular dynamics simulation of nanoparticle diffusion in dense fluids. Microfluidics and Nanofluidics, 2011, 11(4), P. 501-506.
[18] Rudyak V.Ya., Krasnolutskii S.L. Dependence of the viscosity of nanofluids on nanoparticle size and material. Physics Letters A, 2014, 378, P. 1845-1849.
[19] Schofield P. Computer simulation studies of the liquid state. Comput. Phys. Comm., 1973, 5(1), P. 17-23.
[20] Norman G.E., Stegailov V.V. The classical molecular dynamics method: purpose and reality. Nanostructures. Mathematical physics and simulation, 2011, 4(1), P. 31-59 (in Russian).
[21] Norman G.E., Stegailov V.V. Stochastic Theory of the Classical Molecular Dynamics Method. Mathematical Models and Computer Simulations, 2013, 5(4), P. 305-333.
[22] Ashurst W.T., Hoover W.G. Argon shear viscosity via a Lennard-Jones Potential with equilibriu m and nonequilibrium molecular dynamics. Phys. Rev. Lett., 1973, 31, P. 206-208.
[23] Rudyak V.Ya., Minakov A.V. Modern problems of micro- and nanofluidics. Nauka, Novosibirsk, 2016, 298 p. (in Russian).
[24] Belonenko M.B., Chivilikhin S.A., Gusarov V.V., Popov I.Yu., Rodygina O.A. Soliton-induced flow in carbon nanotubes. Europhysics Letters, 101(6), P. 66001.
[25] Golovina D.S., Chivilikhin S.A. Reproduction of the evolution of the liquid front profile in inhomogeneous nanoporous media. Nanosys-tems: Physics, Chemistry, Mathematics, 2017, 8(5), P. 567-571.
[26] Liu C., Li Z. On the validity of the Navier-Stokes equations for nanoscale liquid flows: The role of channel size. AIP Advances, 2011, 1, P. 032108.
[27] Liakopoulos A., Sofos F., Karakasidis T.E. Friction factor in nanochannel flows. Microfluidics and Nanofluidics, 2016, 20(1), P. 24.
[28] Ou J., Perot B., Rothstein P. Laminar drag reduction in microchannels using ultrahydrophobic surfaces. Phys. Fluids, 2004, 16, P. 46354643.
[29] Lauga E., Brenner M.P., Stone H.A. Microfluidics: the no-slip boundary condition. Chapter 19 in: Handbook of Experimental Fluid Dynamics, C. Tropea, A. Yarin, J. F. Foss (Eds.), Springer, 2007.
[30] Andrienko D., Diinweg B., Vinogradova O.I. Boundary Slip as a Result of a Prewetting Transition. J. Chem. Phys., 2003, 119, P. 1310613112.
[31] Laturkar S.V., Mahanwa P.A. Superhydrophobic coatings using nanomaterials for anti-frodt applications (review). Nanosystems: Physics, Chemistry, Mathematics, 2016, 7(4), P. 650-656.
[32] Huang Z. Formulations of nonlocal continuum mechanics based on a new definition of stress tensor. Acta Mechanica, 2006, 187, P. 1127.
[33] Eringen A.C. Nonlocal continuum mechanics based on distributions. International J. Engineering Sci., 2006, 44, P. 141147.
[34] Meshcheryakov Yu.I., Khantuleva T.A. Nonequilibrium processes in condensed media. Part 1. Experimental studies in light of nonlocal transport theory. Physical Mesomechanics, 2015, 18(3), P. 228-243.
[35] Khantuleva T.A., Shalymov D.S. Modelling nonequilibrium thermodynamic systems from the speed-gradient principle. Phil. Trans. R. Soc. A, 2017, 375, P. 20160220.
