Molecular dynamics simulation of nanoscale liquid flows

Microfluid Nanofluid (2010) 9:1011–1031 DOI 10.1007/s10404-010-0612-5 REVIEW Molecular dynamics simulation of nanoscale liquid flows Yuxiu Li • Jinliang Xu • Dongqing Li Received: 8 February 2010 / Accepted: 25 March 2010 / Published online: 14 April 2010  Springer-Verlag 2010 Abstract Molecular dynamics (MD) simulation is a powerful tool to investigate the nanoscale fluid flow. In this article, we review the methods and the applications of MD simulation in liquid flows in nanochannels. For pressuredriven flows, we focus on the fundamental research and the rationality of the model hypotheses. For electrokineticdriven flows and the thermal-driven flows, we concentrate on the principle of generating liquid motion. The slip boundary condition is one of the marked differences between the macro- and micro-scale flows and the nanoscale flows. In this article, we review the parameters controlling the degree of boundary slip and the new findings. MD simulation is based on the Newton’s second law to simulate the particles’ interactions and consists of several important processing methods, such as the thermal wall model, the cut-off radius, and the initial condition. Therefore, we also reviewed the recent improvement in these key methods to make the MD simulation more rational and efficient. Finally, we summarized the important discoveries in this research field and proposed some worthwhile future research directions. Y. Li Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences, Guangzhou 510640, China J. Xu The Beijing Key Laboratory of New and Renewable Energy, North China Electric Power University, Beijing 102206, China D. Li (&) Department of Mechanical and Mechatronics Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada e-mail: [email protected] Keywords Molecular dynamics simulation  Liquid flow  Nanochannels  Nanofluidics Understanding the behavior of liquid flow in nano-environments is of important relevance to various applications such as separation and identification of biological and chemical species, drug delivery, programmable catalysis, selective absorption, and bionic energy (Wheeler and Stroock 2008). Extensive research has been conducted to understand nanoscale flow physics. Experimentally, nanochannels of various sizes made from different materials have been fabricated and tested (Stein et al. 2004; Pennathur and Santiago 2005; Raviv et al. 2001; Fan et al. 2003; Tas et al. 2004; Gogotsi et al. 2001; Stern et al. 1997). Computationally, Monte Carle and molecular dynamics (MD) simulations have become the most popular methods to mimic the nanoscale flow. In this article, we try to review the studies of nanoscale liquid flow using MD simulations. Fluids confined in nanoscale chambers/channels possess significant different properties from fluids in the macroscale and the microscale. In nanoscale, the competitive balance of the forces involved in the fluid system often does not exist. The effects of volume forces which dominate the fluid properties in the large scale are weakened, and the impacts of surface forces are enhanced. These result in many special flow phenomena in nanoscale, for example, the slip boundary condition effect. The hypothesis that the fluid velocity is equal to the wall velocity at the interface of fluid and solid, the so-called non-slip boundary condition, is widely accepted because it has been validated by several experimental and computational researches in large-scale systems. However, the non-slip boundary condition is not valid in nanoscale systems. That 123 1012 Microfluid Nanofluid (2010) 9:1011–1031 is, the velocity of fluid is not equal to that of the solid wall in nanoscale systems, which is called the slip boundary effect. For studying the fluid flow, the continuum theory usually is applied, based on the hypothesis of continuum medium. However, this hypothesis may no longer be valid in nanoscale systems where the molecular free path is comparable to the characteristic size of the fluid-confined space, especially for gas. Therefore, the theories based on the particle movement were developed. The Molecular dynamics simulation (MDS) is a computer simulation method that provides a view of the particles’ motion by computing the interactions of atoms and molecules over a period of time based on known physics. With the development of nanofluidics, the MDS has become a promising simulation method for fluid flow in nanoscale channels. This article is focused on MDS of liquid flows in nanochannels, and is organized as follows. First, the method of MDS will be reviewed briefly. In the following sections, we shall discuss several aspects of MDS of liquid flows in nanochannels: pressure-driven flow, electric fielddriven flow, thermal-driven flow, liquid–solid interaction, slip boundary conditions, and improved MD schemes. Finally, a summary of the past and current research in this area and comments on several promising research directions will be provided. 1 Introduction of MD The MDS starts from the integration of the Newton’s second law. Applying appropriate integral methods, such as the Gear finite-difference algorithm or Verlet method, the basic dynamics parameters such as position, velocity, and interaction force can be determined. The macroscopic physical properties such as pressure, mean velocity, temperature, particle number density etc. can subsequently be determined via statistical mechanics. On the basis of the MD method, the Newtonian equation for each fluid molecule is written as mi Nw N X X ~2 dr ~ ~ F ijw þ ~ F sou~i F þ ¼ ij dt2 j6¼i;j¼1 j 6¼i;j¼1 ð1Þ potential model is applied, the interaction force between a pair of molecules comes from the following relation: o/ij Fij ¼  orij ð2Þ The LJ potential /ij is given by:    r 12 r6 /ðrÞ ¼ 4e  r r ð3Þ where r is the distance between two molecules, r is a molecular length scale, and e is an interaction strength parameter. For example, r = 3.4 9 10-10 m and e = 1.67 9 10-21 J are for liquid argon molecules. The cut-off distance (beyond this distance, it is assumed that the molecules do not interact) is usually set to be 2.5r. The values of the parameters r and e change with the type of wall and fluid molecules. They are named as rwf and ewf when Eq. 3 is applied for the potential interaction between fluid particles and solid wall particles. The cut-off length is changed to be 2.5rwf, correspondingly. Finally, combining Eqs. 3 and 2 with Eq. 1 and applying the shifted potential and shifted force concepts of Haile (1992), we obtain the basic equation for each fluid particle. The shifted force technique ensures that the molecular force between fluid particle i and other fluid particles or other solid wall particles is reduced to zero once the distance between such fluid particle and other particles reaches or is beyond the cut-off distance. The main advantage of doing so is to save a lot of computational time for the force integration.   N r13 r7  r 13  r 7 ~rij P e 24 r 2 r  r 2 rc þ rc jrij j j6¼i;j¼1     N w 13  7 ~ rwf 13 rwf 7 rijw ewf P rwf rwf þ 24 rwf 2 r  r 2 rcw þ rcw jrijw j jw ¼1 ai þ~ F sou~i ¼ mi~ ð4Þ The movement of the fluid particles will be predicted by these basic physics laws. However, solving these equations with special conditions, such as the boundary conditions and initial conditions requires specific methods as explained in the following. w where subscript i represents particle i, and ~i is the unit vector in x-coordinate. The first term of the right-hand side of Eq. 1 is the molecular force due to Lennard–Jones (LJ) potential between particle i and other fluid molecules in the calculation domain, and the second term is the molecular force between particle i and all the solid wall particles j. The last term in Eq. 1 represents the external force, which makes the fluid deviate from the equilibrium, such as gravity, electric force, or pressure. When a two-body 123 1.1 Thermal wall model When the fluid flows inside the solid boundary, if a fluid molecular is quite close to the wall surface, the repulsive force between the fluid particles and solid molecules will push the fluid particle away from the wall surfaces. Physically such a force interaction will never let the fluid particles penetrate the wall surfaces. However, in the real computation, one cannot guarantee that any of the fluid Microfluid Nanofluid (2010) 9:1011–1031 particles does not travel beyond the wall surfaces. This is because it is impossible to use an infinitely small time step. Owing to the above reason, a thermal wall model is introduced as an accessorial boundary conditions. During the computational process through the time series, for most time steps, none of the particles travels beyond the wall surface. Incidentally, there are some time steps that have a couple of particles striking the wall surfaces. In this situation, the thermal wall model is used to ensure that all the fluid particles are confined in the two parallel walls. In other words, the thermal wall model is to reallocate velocity vector of the fluid molecules infinitely closed to the solid wall based on the real physical principle. Take the fluid confined between two parallel solid walls as an example. When a fluid molecular strikes a thermal wall at temperature Tw, all the three components of velocities are reset to a biased Maxwellian distribution. In the Cartesian coordinate system, the three velocity components, after the fluid particles strike the wall surfaces, are given by Alexander and Garcia (1997): qffiffiffiffiffiffiffiffi vx ¼ kBmTw wG ; qffiffiffiffiffiffiffiffi ð5Þ vy ¼ kBmTw w0G ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vz ¼  2kBmTw ln w The positive sign of vz is for the fluid particles striking the bottom wall surface, while the negative sign is for the fluid particles striking the upper wall surface. w is an uniformly distributed random number in (0,1) and wG, w0G are Gaussian-distributed random numbers with zero mean and unit variance. The detailed reentry mechanism at the thermal wall is described by Tenenbaum et al. (1982). 1.2 Initial conditions The initial conditions provide the information of the simulation regions to activate the numerical computation. Taking the crystal fluid confined between two solid walls as an example, the procedure is summarized as follows: • • • The solid wall molecules are arranged in terms of FCC structure, and such an arrangement is not changing during the whole calculation process, because the two solid walls are rigid. The initial positions of fluid molecules are structured according to FCC. However, the fluid molecules will be melting once the applied force acts on the fluid molecules. Randomly distributed velocities in the fluid region are assumed to activate the calculation. However, in order to make sure that the steady state can be easily reached after some reasonable time, the randomly distributed 1013 velocities should satisfy that the total kinetic temperature is equal to the kinetic energy in the fluid calculation domain: N   3 1 X NkB Tfluid ¼ m ð6Þ v2i;x þ v2i;y þ v2i;z 2 2 i¼1 • Values for the initial acceleration ai(0) are determined in terms of the positions ri(0) by computing the force on each atom, and applying the Newton’s second law. If the Gear finite-difference algorithm is used as the integral method, the higher derivatives are also needed. These higher derivatives are assumed to be zero for the initial time step to activate the calculation. iv v x000 i ð0Þ ¼ xi ð0Þ ¼ xi ð0Þ ¼ 0 000 iv v ð7Þ yi ð0Þ ¼ yi ð0Þ ¼ yi ð0Þ ¼ 0 iv v z000 ð0Þ ¼ z ð0Þ ¼ z ð0Þ ¼ 0 i i i For other types of fluid and the solid walls, the solidstate molecular structure should be used as the initial position. 1.3 Periodic boundary conditions Because of the limitation of the computational ability, it is impossible to simulate all of the molecules in the flow system. In order to solve this problem, a feasible method involves the use of the periodic boundary conditions. In this method, a unit cell or computational domain far from the system’s edge is considered. The flow field in the whole system is the periodic recurrence of the flow field of this unit cell. When a molecule passes through one face of the unit cell, it reappears on the opposite face with the same velocity. Take the confined fluid flow between two parallel plates as an example. In the Cartesian coordinate system, the periodic boundary conditions are applied in x and y directions, and the length of the unit cell is Lx in x direction. When a fluid particle is leaving the boundary at x = 0 or x = Lx, another fluid molecule will enter the unit cell at x = Lx or x = 0. The particle entering the unit cell has the same velocity as the fluid molecule that is leaving the calculation domain. This principle ensures that the total fluid molecule number in the unit cell is constant. As mentioned above, the unit cell is a small part of the real fluid flow system. When we use the periodic boundary condition, we assume that the real fluid flow system is the periodical repetition of the unit cell. Therefore, when we compute the molecular interaction force, the image force shall be included. For instance, fluid particle i has the molecular interaction with particle j in the unit cell. However, particle j located at (xj, yj, zj) has a series of image particles outside of the unit cell located at (xj ± nLx, 123 1014 Microfluid Nanofluid (2010) 9:1011–1031 yj ± mLy, Lz) where n and m are integers ranging from 1 to infinity. In order to reduce the computational time for the force calculation of particle i with any particle j in the unit cell and its particle images, Eq. 1 is corrected in terms of Haile (1992): ~ F ij ¼  aX x ¼1 aX y ¼1 ~ij  aLÞ d/ij ðr drij a ¼1 a ¼1 x represents the direction normal to the solid walls. We define the function Hn(zi) such that   Hn zi;j ¼ 1 if ðn  1ÞDz\zi \nDz; ð9Þ otherwise Hn ðzi;j Þ ¼ 0 ð8Þ y where L is Lx in x direction, and Ly in y direction. a is the cell translation vector. Eq. 7 accumulates forces for molecule i interacting with molecule j inside the fluid calculation domain and its images outside of the unit cell. Note that Eq. 7 is also valid for the force integration between particle i with any solid wall particle jw and its images outside of two solid wall regions. A minimum-image criterion has been developed to the calculation of Eqs. 7 and 1. For a cubic container in three dimensions, the minimum image criterion applies separately to each Cartesian component of the pair separation vector. For the x component, we use xij ? xij - axLx where ax = 0 if -0.5Lx B xij B 0.5Lx, ax = -1 if xij B -0.5Lx, and ax = 1 if xij C 0.5Lx (Haile 1992). Similar treatment should be used for yij. where the subscript j represents the j time step. Then, the average non-dimensional number density in nth bin from JN time step to JM time step is qn r3 ¼ JM X N X   r3 Hn zi;j ADzðJM  JN þ 1Þ j¼JN i¼1 123 ð11Þ where zi is the coordinate of the midpoint of the nth bin, and A is the area expressed as Lx 9 Ly. JN and JM are, respectively, the starting and the ending time steps of the parameter averaging. The integrated time interval from JN time step to JM time step is larger than 150s. Longer time interval also results in no change in the results of the parameter averaging. The bin average velocity from JN time step to JM time step is given by uðzi Þ ¼ 1 ðJ M  J N þ 1 Þ 1.4 Parameter averaging By integrating the Newton’s second law equation temporally and spatially, the position, velocity, and acceleration of the molecules are obtained. All of these values are for individual molecules. However, what we are interested in are the macro quantities of the bulk fluid, such as the density, velocity, and pressure. Therefore, statistical thermodynamics must be used to complete the transformation from the micro quantities to the macro quantities. The following summarizes the process of the parameter averaging. Temporally, the appropriate integral methods, such as the Gear finite-difference algorithm, are employed for the dynamics simulation of the molecular system. The time step is chosen from Dt = 0.005s to Dt = 0.0005s depending on the amplitudepofffiffiffiffiffiffiffiffi the applied force, where s is the LJ reduced time s ¼ r m=e: For a large driven force or a large shear rate, the fluid particle has a small characteristic time response to the large applied force. Spatially, the unit cell is divided into many bins. The number of bins is dependent on the fluid molecule numbers. For larger number system, more bins are divided. The macroscopic parameters, such as the number density, mean velocity, temperatures, are defined in the center of each bin. These parameters can be obtained via statistical mechanics, and the detailed procedures are summarized as follows. Taking the confined flow in two parallel plates as an example, we calculate these parameters as a function of Z, where Z ð10Þ PN JM X N X i¼1 Hn ðzi;j Þ j¼JN Hn ðzi;j Þvxi;j i¼1 ð12Þ where vxij is the velocity in x component of particle i at j time step. Using Eq. 13, we compute the bin temperature as Tðzi Þ ¼ 1 P 3kB ðJM  JN þ 1Þ Ni¼1 Hn ðzi;j Þ JM X N h i2 X  Hn ðzi;j Þmi vai;j  ua ðzi Þ ð13Þ j¼JN i¼1 where the superscript a represents x, y or z. At steady state, uy(zi) and uz(zi) are zero for all the bins. 2 Pressure-driven flows Pressurized flow is a general mode of flow in nanochannels. Two fundamental physical processes are involved here, namely, the filtration of fluids into the nanochannel and the transport of fluids inside nanochannel. In the filtration process, the flow resistance and the structure change of the liquid are key factors. In the transport process, the bulk average velocity and molecular structure of liquid are focal points. There are many reported research studies of modeling and numerical simulation of the pressure-driven transport processes in nanochannels. The validity of the continuum theory in nanoscale is examined by comparing the result of MDS with those using Navier–Stokes (NS) equations. Microfluid Nanofluid (2010) 9:1011–1031 2.1 Infiltration flow The fundamental infiltration and defiltration mechanisms of nanofluidics are distinct from that of bulk phases. For smaller channel sizes, such as several nanometers, the relationship between molecular structure and the flow velocity is concerned. Thomas and McGauthey (2009) simulated pressure-driven water flow through 0.83–1.66-nm diameter armchair single-walled carbon nanotubes (SWCNTs) in the NVT ensemble. TIP5P potential was used to model the interactions between water molecules (Mahoney and Jorgensen 2000). The interaction of wall molecules and carbon atoms was modeled using LJ potential of Werder et al. (2003). In order to maintain the temperature of the system, a Berendsen thermostat in the velocity components transverse to the flow direction was employed. The authors found that employing the thermostat to all the three directions did not affect the predictions about the flow. The authors drew two conclusions. First, unlike the continuum theory, the flow enhancement in subcontinuum systems may not increase monotonically with the decrease in flow area. On the contrary, when the flow cross-sectional area is comparable to the size of the liquid molecules, confinement-induced changes to the liquid structure may reduce the flow enhancement. The authors identified the critical diameter of the flow cross section, which is between 1.25 and 1.39 nm, just as shown in Fig. 1. Second, the liquid structure and liquid flow are dependent when the characteristic flow time scale is comparable to the molecular structure relaxation time. As Fig. 1 a Relationship between average flow velocity and applied pressure gradient for the 75-nm-long CNTs. b Variation in flow enhancement factor with CNT diameter. The dashed line between D = 1.25 nm and D = 1.39 nm delineates continuum and subcontinuum flow regimes (Thomas and  McGauthey 2009). Note: By Darcy Dp law, the flow velocity v ¼ c Dp L ; where L means the pressure gradient, and c means the hydraulic conductivity. There exist three factors causing the actual hydraulic conductivity to exceed that calculated from the Poiseuille relation, which is the liquid slip at the solid–liquid boundary, confinement-induced reductions in the liquid viscosity, and subcontinuum changes to the liquid structure. Flow enhancement factor means that the ratio of actual c and the no-slip Poiseuille c 1015 shown in Fig. 2, the water molecular structure varies according to the CNT diameter, from single-file molecular chains (0.83 nm), tilted pentagonal rings (0.97 nm), stacked pentagonal rings (1.10 nm), stacked hexagonal rings (1.25 nm), and to disordered bulklike water (1.39 and 1.66 nm). The authors established the relationship between the molecular structure and the flow enhancement. In the subcontinuum scale (D B 1.25 nm), the transport of molecules in these structured subcontinuum systems is contingent upon the movement of nearby layers. A layer of water molecules that fills an energetically stable CNT surface site can impede the movement of other layers. Such localized limitations on transport will reduce velocity and lower the overall flow enhancement. In most infiltration studies, the liquid molecules were placed in vacuum nanotubes or nanochannels (Hummer et al. 2001; Holt et al. 2006). The absolute vacuum condition is difficult to achieve in real experiments. In a bulk liquid, the gas nanophase that typically consists of a few to several hundred gas molecules often has secondary influence on the liquid motion, in part due to the low gas molecule density that is usually a few orders of magnitude smaller than that of liquids. However, in a nanochannel, due to the radial confinement, the effective gas molecular density significantly increases. If the nanotube diameter is smaller than a few nanometers, then a single gas molecule could make it difficult for the infiltrated liquid molecules to bypass it. Qiao et al. (2007) showed that the gas molecules in relatively large nanochannels could be dissolved in the liquid during pressure-induced infiltration. However, the important role of gas phase was ignored in the previous studies. The authors found that the gas molecules tend to form clusters in relatively small nanochannels, which triggers liquid defiltration at a reduced pressure. The authors studied a cluster of 12 carbon dioxide molecules in water confined by CNTs of four different radii. The CO2 cluster was placed in the middle of a long open-ended tube and surrounded by water molecules, and the MD simulations were carried out at a constant temperature with a constant pressure, using a condensed-phased optimized molecular potential for atomistic simulation studies (Sun et al. 1998). The MDS in this article not only shows the dependence of the defiltration behavior on the pore size and the pore surface properties, but also reveals the inherent interactions between gas and liquid phases, where the effective gas solubility highly depends on the characteristic length scale. In a nanochannel, even a single gas molecule can significantly promote the liquid infiltration and may make the confined liquid unstable. The low effective gas solubility in a small nanochannel leads to the formation of trapped gas clusters during infiltration and the outflow behavior during defiltration, whereas the liquid–solid interaction is dominant after the gas phase is dissolved in 123 1016 Microfluid Nanofluid (2010) 9:1011–1031 Fig. 2 a Water structures inside the 0.83–1.66-nm-diameter CNTs. Carbon atoms are removed for clarity, and the chirality vector for each CNT is provided. b Axial distribution function (ADF), structure relaxation time s, and extent of the ADF, Lt, for each CNT. c P(n, tc), the probability that n molecules cross the system midplane over the characteristic flow time tc, and the CNT permeability a is provided for pressure gradient with 4 9 1014 Pa/m (Thomas and McGauthey 2009) larger nanopores or after the gas phase is diffused out of small nanopores. Qiao et al. (2009) performed MDS of water flowing in silicon dioxide nanochannel by using the COMPASS force field. They considered a straight nanochannel with a diameter of 0.83 nm. One end of the nanochannel was closed. The effective pressure was computed from the state function of water, and the redistribution of water molecules with pressure variation reflects the adiabatic pressure-driven infiltration. The authors analyzed the resistance as and when the water molecules enter the nanochannel. The 123 energy barrier that the water molecules need to overcome in order to enter the nanochannel includes two parts. First part is the change of the free energy of the infiltrated molecules due to each water molecule losing two hydrogen bonds when it enters the nanochannel. Second part is the van der Waals and electrostatic interactions between water molecules and nanochannel atoms. After the infiltration occurred, the driving pressure is still increased almost linearly with the volume difference, which implied that the continued infiltration had to overcome not only the initial energy barrier but also a certain resistance that varied with Microfluid Nanofluid (2010) 9:1011–1031 the infiltration volume. It is referred to as column resistance. In a nanopore with the diameter comparable with the liquid molecule size, most of the infiltrated water molecules are in direct contact with the solid surface, which caused column resistance. Liu et al. (2009) use MD to explore insights into the complex behaviors of liquid infiltration into cone-shaped nanopores. The new contribution of this article is the study of flow in nanochannels with a varying cross-section, while most previous studies of nanofluidics were limited to nanochannels with an invariant cross-section. The liquid and the solid are water and carbon nanocone, respectively. The authors used a series of existed molecular model, potential parameters, and experimental value for building CNTS and CNCS model. The authors used the LAMMPS (Plimpton 1995) (large-scale atomic/molecular massively parallel simulator) package, with the NVT ensemble and the temperature fixed at 300 K. Water molecules are modeled by the rigid extended simple point charge potential SPC/E (Berendsen et al. 1987), while the carbon– oxygen LJ parameters are extracted from the experimental low-coverage isotherm data of oxygen adsorption on graphite. The following are the findings of this article. Under ambient conditions, the applied pressure needs to be higher than a critical infiltration pressure so that water molecules may flow into a nanotube, and such an infiltration pressure is higher for smaller nanotubes. Spontaneous infiltration initiates when the apex angle of nanocone is large despite the hydrophobic nature of the solid phase. However, water cannot wet the entire nanocone. Furthermore, with the same blunt cone angle, a long nanocone and a short nanocone may exhibit distinct wetting properties. Enlarging the apex angle always leads to a lower pressure for liquid to reach the same pore size, and hence the easier infiltration. 2.2 Special statistical method for low speed flows In Sect. 1.4, we introduced the transformation method from the micro quantities to the macro quantities. In the MDS method, the total velocity of molecules consists of the thermal motions velocity and the liquid flow velocity caused by the externally applied forces. The total liquid flow velocity is calculated by summing up the thermal motion velocity and the external force-driven flow velocity, temporally and spatially. In general, the thermal motion velocity of molecules is on the order of several hundred meters per second. Theoretically, if the time averaging is performed for infinitely long time, then the summation of the random thermal motion velocity will be zero. However, in practice, the time averaging is conducted over a limited period of time, and hence the summation of the random thermal motion velocity may not be zero. Therefore, when 1017 the external force-driven flow velocity is \1 m/s, the thermal motion velocity will affect the average value dramatically so that the MDS method cannot extract the true bulk liquid flow velocity. In these cases, the bulk liquid flow velocity is overestimated. This has been verified by the error theory and simulation analyses (Zhang and Li 2007). In order to extract the true liquid bulk flow velocity caused by the externally applied forces, Zhang and Ely (2004) proposed a new linearized algorithm, where the total flow velocity is the linear accumulation of the flow velocity increment of each time step. In each time step, using the linear development of the particle interaction force function in MDS method, the force increment associated with the external force can be separated from the total increment of the particle interaction force. Then, the increment of the flow velocity can be calculated. The total flow velocity is the accumulation of the increments over all time steps. When the interval of time step is smaller than 10-16 s, the algorithm is convergent. This interval of time step is smaller than that applied in the conventional MDS method. 2.3 Special applications Conventional injection involved in spray combustion devices has been widely studied. In recent years, numerous studies have attempted to predict nanoscale jet breakup by using non-equilibrium MD (Branam 2005; Ludwig 2004; Micci et al. 2000; Moseler and Landman 2000; Shin 2006). Kim (2007) developed a non-equilibrium MD code and applied it to simulate nanojetting and electrowetting phenomena. The liquid and solid are argon and platinum, respectively. The sum of the LJ potential and the Coulomb potential was adopted. For liquid–solid interaction, the LJ potential was modified by the intermolecular parameters, a and b, between liquid and solid. In this study, the nanojetting mechanism from droplet ejection, breakup, wetting, and the drying process was identified by the simulation, with a = 0.5 and b = 0.5. Three types of typical wetting phenomena were simulated to understand the effects of intermolecular forces on the surface tension and contact angle. Three different wettabilities (full wetting, intermediate wetting, and no wetting) are obtained under three potential parameters condition of b = 1.0, b = 0.25, and b = 0.1. In addition, a conceptual nanopumping system was simulated using the electrowetting phenomenon. For the homogeneous flow in nanotube, the author found that those wall effects become more serious as the diameter decreases from 12 to 3 nm. For flow in a charged nanotube, the negative surface charges attract the counterions to the wall surface. On the other hand, repulsive forces are induced between fluid and wall particles, causing the slip flow at wall surface. Therefore, the authors 123 1018 drew a conclusion that the sign of wall surface charge determines the flow resistance that is inversely proportional to flow rate while applying the same driving force. 3 Electric field-driven flows Pressure-driven flows in nanochannels need a large pressure gradient, which is difficult to achieve experimentally. Therefore, other driving forces, such as a surface tension difference (de Gennes et al. 2003), a chemical gradient (Chaudhury and Whitesides 1992), a temperature gradient (Linke et al. 2006), and electric field (Kuo et al. 2001; Mitchell et al. 2000; Daiguji et al. 2004; Ermakov et al. 1998; Freund 2002; Qiao and Aluru 2003) have been considered. 3.1 Fundamental research Electroosmotic flow is an important fluid transport mechanism in nanofluidic system. In general, there is an electrical double layer (EDL) near the solid–liquid interface when an ionic solution is in contact with a charged solid surface. If an external electrical field is applied in the direction tangential to the surface, the fluid will be dragged by the moving ions in the EDL and an electroosmotic flow is generated. The key parameters of the electroosmotic flow in nanochannels are the zeta potential, the density distribution of the solution and ions, and the velocity profile in the channels. There are several other factors affecting the above physical parameters, such as the surface roughness, the hydrophobicity of the solid wall, the external electric field strength, the charge density, and the polarity of the surface charges. In the MDS of electroosmotic flow, the Coulomb potential should be added to the LJ potential for the charged particle pair. The general expression is " # r 12 r 6 qi qj Vðri;j Þ ¼ e þ ð14Þ þ e0 rij rij rij where qi, qj are the charge of the charged particle pair. In general, the electrostatic interactions were computed using the particle-mesh Ewald (PME) method (Darden et al. 1993) for the Coulomb interactions. For electroosmostic flow, the charged atoms were driven by the electric filed force. Therefore, the Newton’s Second Law should be written as follows:   N X o/ rij mai ¼ þ Eqi ð15Þ orij j¼1;j6¼i Modeling and simulation of electroosmotic flow in nanometer scale channels can address many fundamental 123 Microfluid Nanofluid (2010) 9:1011–1031 issues currently encountered during the design of nanofluidic systems (Kuo et al. 2001). By far, the most popular method to simulate the electro-osmotic flow is based on the classical continuum theories, such as the combination of Poisson–Boltzmann (PB) equation for the ion distribution and NS equations for the fluid transport (Mitchell et al. 2000; Daiguji et al. 2004; Ermakov et al. 1998). However, some important information dominated in nanoscale would be ignored while using the continuum method. For instance, the ions are modeled as point charges, and the water is modeled as a dielectric continuum in the PB equation. The molecular nature of the ions and the discreteness of the water molecules are not accounted for in the PB equation (Darden et al. 1993). Several articles have shown the discrepancies between the MD results and the continuum ones. The MDS can provide a qualitative understanding of the various physical processes without the assumptions built into the continuum theory. Freund (2002) found the profile differences of the ion density and the water velocity near the walls derived from the two methods. Qiao and Aluru (2003, 2005) published a series of articles on this topic. They modified the PB equation by adding MDS and simulated more realistic model systems. They found that the direction of electro-osmotic flows is opposite to that predicted by the continuum equation under the condition that ions have the same sign of charge as the walls. They also found that the ionic distribution is different in positively and negatively charged channel walls due to different hydrogen bond patterns. Moreover, this phenomenon cannot be observed by using continuum method. Joly et al. (2004) reported that the zeta potential increased when the slip near hydrophobic walls existed. Therefore, the continuum theory may be invalid at smaller nanoscale systems. In any case, MDS is a powerful tool to investigate the flow and ion distribution in nanoscale. Qiao and Aluru (2005) investigated the influence of the finite size of the ions and the water molecules on the ion distribution, and the effect of surface–fluid interactions on the electroosmotic transport. The simulation system consists of a slab of KCl solution sandwiched between two oppositely charged channel walls. The models in this article are approximately the same as that used in Darden’s study (Darden et al. 1993). The ions distribution in nanochannels is not coincident with the Boltzmann distribution. Figure 3 shows the distinct results derived from the continuum method and the MD simulation. For Cl- ions, the first concentration peak is about 81% higher compared to the PB theory prediction, and the concentration does not decrease monotonically to its value in the channel center as PB equation predicts. For K? ions, the concentration is substantially different from that of the Cl- ion concentration. Particularly, the K? ion concentration appears to have Microfluid Nanofluid (2010) 9:1011–1031 Fig. 3 K? ion and water concentration profiles across the channel (Qiao and Aluru 2005) two distinct peaks. The authors explained that the first peak is related to the interaction of K? ion and the walls. The second peak is mainly caused by the interaction between the ion and its hydration water molecules. In an electrolyte solution, the water molecules near an ion are loosely bonded to the ion due to the strong charge–dipole interactions between the ion and the water molecules. These water molecules are called hydration water. Because of this interaction, the water distribution will affect the ion distribution. It is well known that the water density oscillates near the wall, and so there is the second peak of ions. For water and ion transport, the MD results showed that the electro-osmotic velocity in a negatively charged channel is much higher than that in a positively charged channel. The authors calculated the water viscosity and found that the viscosity is much higher near a positively charged surface compared to that near a negatively charged surface. The authors established the relationship between the water molecular structure, especially the hydrogen bonding, and the various properties. The hydrogen bonding is referred as the H–H bonding between two water molecules. In the positively charged surface, the number of hydrogen bonding is much larger than that in the negatively charged surface. The larger the number of the hydrogen bondings, the higher the water viscosity. Kim and Darve (2006) studied the effects of solid surface roughness on the density profile, the zeta potential, the ionic distribution, and the velocity profile. One of the key issues in electroosmotic flow related to the roughness was proposed in this article: The zeta potential is influenced by the roughness. It should be noted that the surface potential and the zeta potential are different. The zeta potential is defined as the electrostatic potential at the slip plane. In electrokinetic theory, the surface potential is approximated 1019 by the zeta potential because it is impossible to measure the surface potential. The zeta potential would deviate significantly from the surface potential when the influence of surface roughness increases so that the location of the slip plane is changed. The structure of electric double layers on charged rough walls was also addressed in several books (Bikerman 1970; Dukhin and Derjaguin 1974). In the computational model (Darden et al. 1993), the authors applied several simplified processes, which may influence the reliability of the conclusion. Their key assumptions are listed as follows. The walls were modeled as a simple cubic crystal structure and fixed at each lattice site. When constructing rough wall channels, the extra partial layers of wall atoms were added only in two-dimensional direction. Only counterions were included for simplicity and computational efficiency. The authors combined several existed models together to simulate the effect of the roughness to electro-osmotic flow. The sodium ions were modeled as point charges with the LJ potentials. For water molecules, the authors used the extended simple point charge model [SPC/E (Berendsen et al. 1984)]. The GROMACS code (Spoel et al. 2004) with GROMOS-96 force field was used except for wall atoms. For the interactions of wall atoms with other types of atoms, the authors applied the LorentzBerthelot combination rule. For the calculation of electrostatic potentials, the PME method (Darden et al. 1993) was used. In order to constrain water molecules, the SETTLE algorithm (Miyamoto and Kollman 1992) was used. Very high electric field strength (0.1 V/nm) was applied for avoiding relatively large thermal noises which result in large statistical errors. Because of the layering of water molecules and the finite-size effect of ions and walls, the ionic distribution deviates from the Boltzmann distribution. The polarization density, expression of the structure of water molecules and ions, was computed. The polarization density component in the direction normal to the wall dominates over the other components. This result indicated that there existed the preferred orientation of water molecules. For rough walls, the preferred orientation is more evident inside the groove. With the use of the velocity autocorrelation function, the residence time and the diffusion coefficient showed that the diffusion of water and ions in the direction normal to the wall is reduced near the wall and further reduces inside the groove. The zeta potential would decrease as the period of surface roughness decreases or the amplitude of surface roughness increases. 3.2 Pumping application Fluid can be transported by electric or magnetic fields if there are ions or magnetic quanta in the fluid. However, pure water is charge-neutral and has no magnetic moment. Understanding and controlling the water transportation 123 1020 Microfluid Nanofluid (2010) 9:1011–1031 Fig. 4 Introduction to the main system. The three circles on the right side of the cylindrical channel represent the charges that drive the water flow through the nanochannel. a Side view of the main system. b Top view of the same arrangement (Gong et al. 2007) through nanochannels is of great importance for designing novel molecular devices and has wide applications (Whitesides 2006; Whitby and Quirk 2007; Regan et al. 2004; Service 2006; Holt et al. 2006; Bourlon et al. 2007; Besteman et al. 2003; Ghosh et al. 2003). Gong et al. (2007) proposed a design for a molecular water pump on the basis of MD simulations. Figure 4 shows the principle of the design. The model system is the uncapped armchair SWCNT. Three positive charges were arranged outside of the nanotube along the longitudinal direction, two of which were of symmetrical distribution versus the mid-point of the nanotube length. The third charge is positioned at the bottom of the nanotube. The water dipole orientation angle along the length of the channel follows the regular distribution: the water dipoles near the bottom end point down, whereas the water dipoles near the top end point up, showing a bipolar order. Such bipolar structure could produce biased electrostatic potential, which could push water molecules through the nanochannel. The authors compared the effectiveness of other systems (such as only one charge at the end of the nanotube) with the above system (called the main system) by analyzing the water dipole orientation angle and corresponding potential structure. For the main system, the change of the water dipoles orientation was slower than the other systems. Correspondingly, there was less deep potential well, and hence, the water permeation across the channel was easier. Meanwhile, it should be noted that the single-file structure of water inside the channel was crucial to the water dipole orientation distribution being shielded from thermal fluctuations. For wider channels, the thermal fluctuations will reduce the regular bipole distribution, which smoothes the asymmetries of the potentials and weakens the pumping abilities. 3.3 Particle separation or screening Selectivity in ion transport through nanochannels is an important branch of nanofluidics and has important applications. Banerjee et al. (2007) used MDS methods to demonstrate the utility of CNTs for water and ion separation 123 from a salt solution. Water was modeled using the simple point charge (SPC) model (Berendsen et al. 1981), ions using the so-called primitive intermolecular potential model (Chandrasekhar et al. 1984; Jorgensen et al. 1982), and CNTs using the parameters provided in Turner et al. (2002). The new contribution of this study is to describe the role of alternating spatial or temporal charge patterns in accomplishing the controlled intake of either ions or polar solvents such as water. The results showed that when spatially or temporally alternate charge patterns were applied to CNTs, the ion intake is minimal while allowing water passing through the CNTs. The implication is that CNTs could be used for the preferential intake of water molecules as waterpurifying devices (or conversely as ion concentrators for the solutions external to CNTs). The CNTs have been used for various sensing applications (Mashl et al. 2003; Li et al. 2004; Saito et al. 1992; Benedict et al. 1995; Chen et al. 2005; Li and Lin 2006), such as chemical sensors, for selective detection of nitrogen oxide and ammonia (Mahar et al. 2007). Xu and Aluru (2008) investigated the screening effects of semiconducting and metallic single-walled CNTs (SWCNTs) when the water molecules and various charged ions pass through the nanotubes. It is shown that metallic SWCNTs have much stronger screening abilities than semiconducting SWCNTs. The authors thought that relatively fewer efforts have been invested to study the polarizability and the screening effects of CNTs under electric fields generated by water dipoles and charged ions inside SWCNTs. In this article, the structure of water and ions inside the SWCNTs is obtained by performing MD simulations. The results from MD simulations are then used in tight-binding (TB) simulations to investigate the screening effects of metallic and semiconducting SWCNTs. The authors considered 5.3 nm long hydrogen-terminated metallic and semiconducting CNTs attached to reservoirs at both ends. Both the CNTs have similar diameters but different electric properties. This article demonstrated how to integrate the MD method with self-consistent TB (tight binding) calculations and ab initio methods. In the MD simulations, the initial systems were equilibrated for 100 ps under a constant pressure Microfluid Nanofluid (2010) 9:1011–1031 of 1 bar and then switched to NVT ensemble (the number of particles (N), the volume (V), of each system in the ensemble are the same, and the ensemble has a welldefined temperature (T), determined by the temperature of the heat bath.). Non-equilibrium MD simulations were performed for ion transport to obtain the coordinates of water molecules and the ions at intervals of 1 ps for TB simulation. Yang (2007) studied water structure in negatively charged model cylindrical nanopores, as well as the partitioning of positive ions of increasing size into the pore interior by using MD simulations. The authors used a piece of (5,5) armchair CNT pore comprising 100 carbons as a model cylindrical pore. Charges were distributed uniformly on the pore atoms. The aim of this article is not to model specifically a CNT but to use its cylindrical geometry as a model for a cylindrical pore. Simulations were performed using AMBER6.0, TIP3P model of water, and LJ description of carbon. Ions were represented as LJ spheres with ion charge placed at the center. Ion–water and ion– carbon LJ interactions were calculated using Lorentz– Berthelot mixing rules. The PME method was used to calculate the electrostatic interactions with grid spacing of 1 Å. Temperature and pressure were maintained at 300 K and 1 atm, respectively, using the Berendsen method. A time step of 1 fs was used in the simulations. This study found that the nanopores display selectivity toward partitioning of the larger cations over that for the smaller ions. The absence of the equilibrium kinetic simulations shows that the partitioning is significantly slower for Na ions compared to that for K and Cs ions, especially for lower charge densities on the pore. Thermodynamic and kinetic studies collectively suggest the presence of a barrier for partitioning of cations into the nanopore interior. Such a barrier is expected to slow down the flow rates of ions through nanopores. 4 Thermally driven flow When particles are added into a quiescent fluid that is subjected to a temperature gradient, they tend to migrate toward the cold end. The phenomenon is called thermophoresis. The origin of the thermophoretic force was first suggested by Maxwell (1879). The mechanism of the phenomenon is as follows. Since the molecules from the hotter region have a higher momentum tangential to the surface than those from the colder region, the total interaction may result in the thrust on the wall or the particle toward the cold region. In contrast to gases, the theory and the experiment of thermophoresis in liquids are far from well established even in the simplest case of a neutral particle in a non-conducting liquid. The difficulty for 1021 thermophoresis in liquids lies in the fact that there may be some mechanism responsible for the thermophoretic flow other than the kinetic interaction. Nevertheless, thermophoresis has practical implications in the fields such as the bore stability in the oil industry, the particle accumulation in a thermal boundary layer, and the transport of polymeric molecules in the microchannels. Han (2005) studied thermophoresis by MD. The authors believed that the kinetic theory on rarefied gases cannot predict a finite thermophoretic velocity observed in the liquids. Moreover, the intermolecular interaction potential plays a dominant role in the liquid state. Since the liquid atoms were under the influence of other particles’ potential and the fluid may become inhomogeneous at the solid– liquid interface, it is desirable to understand the dynamics in the solid–liquid interfacial region more clearly. A rigid and neutral particle with diameters of about 1 lm was immersed in quiescent non-conducting liquid. The temperature difference of 1 K is applied over a distance of 100 lm. The authors divided the computation domain into two sub domains: quiescent liquid domain far from the particle and creeping flow domain close to the particle. Thus, the computational cost would be reduced dramatically. In the quiescent liquid domain, analytical method was used. In the domain near the particle, MD method was used. By combining the two methods, it is shown that the temperature gradient applied to the non-conducting liquid induced the pressure gradient tangential to the neutral particle surface. The model didn’t consider any charge or charged materials, whose effects may be very significant in some cases. The dependence of the EDL as well as other properties on temperature must be accounted for in many cases. The phenomena may closely depend on the characteristics of the liquid structures developed in the interfacial region according to the given conditions. In the colder region, a higher bulk density and the stronger influence of the solid over the thermal mobility induce a higher rate of increase of density and a higher tangential pressure. As a result, a negative gradient of the tangential pressure with respect to temperature is developed and the particle moves to the cold end. Two years later, Han (2008) published another article. This article focused on the potential applications based on thermophoresis, such as pumping operation. Figure 5 shows the principle of the pump. Owing to the presence of solid, there always exists a pressure difference as shown by the dashed line in Fig. 5a. The pressure difference collectively results in an unbalanced force occurring only in the interfacial region Fig. 5b). The author studied four different types of systems using MD simulation. The four systems differed mainly in the channel configuration and the way that the gradient was applied. The author found some technical obstacles when designing the device and found 123 1022 Microfluid Nanofluid (2010) 9:1011–1031 Fig. 5 Illustration of the origin of thermophoretic force: a MD results of the tangential pressures with respect to the distance from the solid surface; b the integrated effect of the unbalanced pressure occurring near the surface (Han 2008) the solution for the difficulties. First, it was found that it is difficult in imposing a sufficiently large temperature gradient in the small scale. In this case, the features like thermal-contact-resistance at the interface need to be included in design considerations. The second was a limited flow-development under an increased viscous drag in the narrow channel. The first difficulty may be overcome by connecting a single pumping unit of a modest gradient with another in sequence. Assembling a large number of the units would provide a necessary pumping power with less limitation. The main idea of the novel design is as follows. The signs of the temperature gradient and the thermophoretic force alternate periodically. However, the author produced a significant amount of thermophoretic force using a weakly interaction wall or fluid in the region of temperature gradient of a certain sing (for example, minus in this case). Then, a larger amount of flow is developed by the forces in the region where the temperature gradient is of the opposite sign (plus in this case), and the net flow would be of a significant amount. The major part of the simulation system consists of liquid argon and platinum wall. The argon atoms interact with each others through the LJ potential. A harmonic potential was used for the atom–atom interactions in the Pt lattice. The thermally insulating slip wall was introduced, and their atoms were fixed in space in the same FCC lattice and interacted with fluid atoms through only the repulsive potential. When the temperature of the Pt walls needed to be controlled, the scheme based on the so-called ghost particle was used, which let the energy in and out to maintain a constant temperature. The ghost particles are located in boundaries of the simulation system and play the role of a thermal reservoir. By using the readjusting temperature method, only the velocity components in y and z directions, which are normal to the flow, are involved in the scaling to minimize the influence on particle drift. 123 5 Slip boundary condition The flow behavior at the submicron scale strongly depends on the boundary conditions at the solid–liquid interface. A number of experimental studies showed that the no-slip boundary condition is not valid in the submicron scale. The most popular Navier model relates the slip velocity and the shear rate with a proportionality coefficient, the slip length, which is determined by the linear extrapolation of the fluid velocity profile to zero. Many key parameters, such as wettability, surface roughness, coupled fluid structure, and shear rate, have effect on the magnitude of the slip length. There are extensive literature reporting research studies in this area (Churaev et al. 1984; Cottin-Bizonne et al. 2002; Joly et al. 2006; Pit et al. 2000; Schmatko et al. 2005; Bonaccurso et al. 2003; Sanchez-Reyes and Archer 2003; Schmatko et al. 2006, Vinogradova and Yakubov 2006; Migler et al. 1993; Horn et al. 2000; Craig et al. 2001; Choi et al. 2003; Koplik et al. 1989; Thompson and Robbins 1990; Thompson et al. 1995; Thompson and Troian 1997; Jabbarzadeh et al. 1999; Barrat and Bocquet 1999; Cieplak et al. 2001; Priezjev and Troian 2004; Galea and Attard 2004; Priezjev 2007; Khare et al. 1996; Koike and Yoneya 1998; Niavarani and Priezjev 2008a, b). Niavarani and Priezjev (2008a, b) investigated the slip phenomena in thin polymer films confined by flat or periodically corrugated surfaces by MD and continuum simulations. The LJ potential was employed for the wall–fluid interaction. In addition to the LJ potential, the neighboring monomers in the polymer chain interact with the finite extensible nonlinear elastic potential. In this study, to avoid a bias in the flow direction, the random force and the friction term were added to the equation of motion in the y direction. By using this method, the system will not be heated by the friction, which is not reasonable. The findings of this article are as follows. For atomically flat walls, the slip length Microfluid Nanofluid (2010) 9:1011–1031 passes through a shallow minimum at low shear rates and then increases rapidly at higher shear rates. In the case of periodic surface heterogeneities with the wavelength larger than the radius of gyration, the effective slip length decays monotonically with increase in the corrugation amplitude. At low Reynolds numbers, the inertial effects on slip boundary conditions are negligible. When the corrugation wavelengths are comparable to the radius of gyration, polymer chains stretch in the direction of shear flow above the crests of the surface corrugation, while the chains located in the grooves elongate perpendicular to the flow. As we all know, the substance consists of molecules which are arranged according to a certain kind of spatial structure. The arrangement of the molecules will form different crystal planes where the properties vary. From a microscopic point of view, the lattice arrangement of wall atoms at the liquid–solid interface will influence the liquid hydrodynamic characteristics. The surface orientation or lattice plane is a solid wall characterized by its Miller indices, which can be referred as a control parameter in nanochannel flows. Soong et al. (2007) studied Couette and Poideuille flows and explored the effects of wall lattice– fluid interactions on the hydrodynamic characteristics in nanochannels. The liquid confined between two parallel planes is argon and LJ potential was used to model the potential between fluid–fluid and fluid–wall atoms. The thermostat was applied only to the wall atoms and the fluid particles within the thin layer closest to the walls. There exists another method which applied the thermostat to the molecules in the entire fluid region. The authors employed a hypothetical solid wall of FCC lattice structure with a nonlinear spring model of interaction between wall atoms. Taking the FCC crystal as an example, as shown in Fig. 6, one can see that the molecules’ arrangement is different in 1023 different orientations. For FCC crystal lattice, there exist three types of surface orientation: h111i, h100i and h110i. The surface orientation will affect the flow behavior. Figure 7a and b shows the Couette and Poiseuille flows of the usual LJ fluids. For both flows, among the three lattice planes, the flow with a channel wall of the surface h100i is the most insensitive to the variation of the flow orientation. Khare et al. (1997) had compared the two approaches and recommended the former one. By using this method, the computational cost will be reduced. The MD results revealed that the flow orientation with respect to the wall lattice structure had very remarkable influence on the fluid slippage and in turn the velocity profile and the flow rate. With increasing wall–fluid interaction potential or wall density, the density layering could be enhanced. Among the three surface orientations, the slip length for the smoothest wall, FCC (111) is the largest, especially at a high wall density. By employing the wall surface FCC (110) as a model, it was found that flow along the molecular trenches of the surface led to the largest slip length. The MD simulations disclose that the influences of the flow orientation are quite distinct for various lattice planes. The fluid slippage and the flow rate will be different if the channel walls are of different lattice planes and/or the fluid flows over the wall at a different orientation angle. Yang (2006) carried out non-equilibrium MD simulations to investigate the effect of surface roughness and interface wettability on the nano-rheology and slip boundary condition in nanoscale. The roughness of the wall was performed by decorating the wall with periodic nanostrips of a certain height and width in x and y direction. The simulation results showed that the surface roughness suppresses the fluid slip for both of the hydrophilic and hydrophobic surfaces of nanochannels. For smooth and Fig. 6 Flow channel model, surface orientation, and flow orientation. a Nanochannel model, b face-centered cubic crystal lattice; c surface orientation or crystal plane; d flow orientation (Soong 2007) 123 1024 Microfluid Nanofluid (2010) 9:1011–1031 Fig. 8 Snapshots of liquid argon flowing through nanochannels. a Hydrophobic surface; b hydrophilic surface (Yang 2006) Fig. 7 Effects of flow orientation on a Couette and b Poiseuille flow velocity distributions (Soong 2007) hydrophobic surfaces, an air gap or nanobubble was found to appear at the fluid–solid interface, which resulted in the apparent slip velocity. The velocities of fluid molecules are enhanced with the driving force, while the nanostructures of fluid densities are independent of the driving force. The results showed that the slip length increased with the average fluid number density for hydrophobic nanochannels, while the effect of the average fluid number density on the slip length is not obvious for hydrophilic walls. The advantage of this model lies in that the wall molecules are not fixed, and the interactions among the wall molecules were determined by the LJ force and spring force. The surface wettability or wall–fluid interaction was adjusted by the length scale in LJ potential. In order to keep the temperature of the system constant, a Gaussian isokinetic thermostat is used to adjust the wall temperature to a constant value at each time step. For a nanochannel with a hydrophobic surface, as shown in Fig. 8a, an air gap or nanobubble exists near the fluid–solid interface. The experimental results of Tyrell and Attard (2001) were similar to the simulated results. For flow through a nanochannel with a hydrophilic surface, due to the strong solid– liquid interactions, Fig. 8b shows that the fluid molecules are orderly adsorbed in the vicinity of the solid interface. Priezjev and Troian (2006) investigated the influence of surface roughness on the slip behavior of a Newtonian liquid under a steady planar shear by employing three 123 different approaches, Stokes flow calculations, MD, and a statistical mechanical models. The LJ potential was adopted in this article to simulate the interaction between the particles. Each wall consists of two layers the atomic sites of which defined the (111) plane of a face-centered cubic (FCC) lattice. The equations of motion governing the coordinate and momentum values of the ith particle were integrated using a standard fifth-order Gear predictor–corrector scheme. The fluid was subjected to steady planar shear by translating the upper flat wall at a constant speed along the x-axis. The lower wavy wall remained stationary. The wavy wall was constructed by vertically displacing all particles within the (111) plane of the FCC lattice by a distance varying with a sine wave. The thermal reservoir and driven force were absent for the x- and z- components. By using this simplification, the computational cost can be reduced. However, this method cannot deal with the high shear rate condition. Therefore, in this article the shear rate never exceeded 0.033 s-1. The authors found a key parameter ka, where k is the wave-number, and a is the amplitude of a sinusoid wave. For Newtonian liquids subjecting to steady planar shear, the slip length decreased monotonically with increasing values of ka. A no-slip or even a negative-slip boundary condition could be obtained for sufficiently large values of ka. In small ka and large k, the authors found MD as a better method than Stokes method. The potential energy exerted by the wavy wall exhibits two competing wavelengths, one set by the lattice sites comprising a flat crystalline wall and the other set by the topological modulation describing a periodically roughened wall. In this study, the authors combined three methods spanning multiple length scales ranging from Microfluid Nanofluid (2010) 9:1011–1031 molecular to macroscopic dimensions and provided a detailed picture of the influence of periodic roughness on the slip length. Martini et al. (2008) tried to explore the molecular mechanisms of liquid slip. The MDS data and the results from an analytical model which support two mechanisms of flow slip were presented in this article. The authors distinguished two types of slip under the conditions of low levels of force and high levels of force. At low levels of force, the potential field generated by the solid creates a ground state which the liquid atoms preferentially occupy. Liquid atoms hop through this energy landscape from one equilibrium site to another. Visual evidence of the trajectories of individual atoms on the solid surface supports the view of localized hopping, independent of the dynamics outside a local neighborhood. The authors call this defect slip. At higher levels of force, the entire layer slips together, obviating the need for localized defects and resulting in the instantaneous motion of all atoms adjacent to the solid. The appearance of global slip leads to an increase in the number of slipping atoms and consequently an increase in the slip length. The wall atoms were tethered to its lattice site by a linear spring, which is more realistic. The model adopted a united atom model allowing both bond bending and torsion into consideration. A key factor that differentiates this study from others previously reported is that the fluid density was determined using a grand canonical Monte Carlo simulation. The advantage of this method lies in that as the chemical potential, rather than particle number, is fixed, the liquid is allowed to reach equilibrium if the channel were filled from a reservoir at the given temperature. 6 Other liquid–solid interactions In addition to the slip flow boundary conditions, a key topic in the studies of liquid flow in nanochannels by MDS is how the interaction between the fluid molecules and the solid walls will affect the liquid flow. Xu and Chen (2006) presented a MDS of molecules translocation through CNTs in vacuum and in aqueous environment. The authors compared the dynamic characteristics between the rigid molecules (C60 molecule and a finite segment of SWNT with a smaller diameter than the outer channel) and semi-flexible molecules (single-stranded and double-stranded DNA oligonucleotide molecules). The authors found that for rigid molecules the van der Waals interaction between the molecules and the CNT channel can induce rapid molecular oscillation inside the CNT channel in vacuum with little dissipation. The dynamics of spontaneous insertion into the water-solvated CNT channels was found to depend sensitively on the details of 1025 the interplay between van der Waals, hydrophobic and hydrogen bonding interactions in the hybrid nanotubewater molecule system. The authors used the AMBER99 force field to model atomic interactions in the buckyball molecule, the SWNT channel, and the ssDNA and dsDNA molecules. For studying nanotube molecule interaction, the carbon atoms of the nanotube and C60 molecule were modeled as uncharged LJ particles, corresponding to sp2 carbon–carbon bond length and bond angle in the SWNT and C60 molecules were maintained by applying harmonic potentials. For simulating the encapsulation process in aqueous environment, the CNT and the CNT- molecule complex were solvated in a water reservoir using the well known TIP3P water model. The study method was applied to evaluate electrostatic interactions, and sodium ions were added as counter ions to compensate negative net charges on the DNA molecules. Atomic interactions between the ssDNA and dsDNA molecules and the water solvent were modeled using the AMBER99 force field. The results showed that in the absence of water solvation, the van der Waals interaction between the molecule and the nanotube wall can induce a rapid spontaneous encapsulation of the molecule inside the nanotube channel. The encapsulation process is strongly impeded for nanotubes dissolved in water due to the competition between the van der Waals, hydrophobic and hydrogen bonding interactions in the nanotube/eater/molecule complex. Sofos et al. (2009) calculated the transport properties of liquid argon flowing through a nanochannel formed by krypton walls. In this study, wall atoms are constrained by LJ potential and an elastic spring potential, unlike other articles where the wall atoms are fixed. This approach can simulate more accurately the condition of wall atoms. An external driving force was applied along the x direction to every fluid particle. The magnitude of the external applied force was selected in order to avoid non-linear variations in fluid temperature induced by the flow. The wall atoms absorb the increase in kinetic energy of the fluid atoms and Nose–Hoover thermostats at the thermal walls were employed in order to keep the system’s temperature constant. The authors used two independent thermostats, one for the upper wall and another for the lower wall in order to achieve better thermostat of the wall atoms. The results show clearly the existence of a critical width of the nanochannel, below which the transport properties are affected in comparison with bulk properties. The authors explained the phenomenon as follows. For small width values, diffusion coefficient is highly anisotropic, the component normal to the wall being the smaller one. For the same width range, diffusivities are higher in the central layers than those close to the walls. Diffusivity is isotropic in planes parallel to the walls. The calculation of the diffusion 123 1026 coefficients in distinct layers parallel to walls revealed that diffusion is higher in layers close to the center of the channel and decreases for the layers adjacent to the walls. Shear viscosity values are affected by the channel width. The smaller the channel width, the larger its value. The existence of this critical width is attributed to the interaction of the fluid with the walls since at smaller channel widths the effect of the interactions of the walls on the fluid extends practically to all the fluid molecules in the channel, modifying considerably the fluid properties in comparison to the bulk. Sofos (2009) used NEMD to investigate the effect of periodic wall roughness on the flow of liquid argon through krypton nanochannels. The particle interactions were described by the LJ potential. The lower wall of the channel is smooth. The rough upper wall is constructed by ‘‘adding’’ extra wall atoms to form periodically spaced rectangular protrusions. The wall atoms are bound on FCC sites and remain in their original positions by an elastic spring force. The temperature is kept constant at with the application of Nose–Hoover thermostats. The authors presented a potential energy map that shows the existence of low and high potential energy regions inside the rectangular cavities. The potential energy field varies with the length of the grooves and affects the trapping of fluid atoms inside a roughness cavity. The number density contour plots show that fluid atoms tend to be localized inside rough wall cavities. As the cavities become narrower, this localization increases. The maximum value of streaming velocity in the center of the channel is not significantly affected by the presence of roughness. Velocity profiles are suppressed in the upper half of the channel where the rough wall is present. The result showed that as the wall cavities become narrower velocity values inside the cavities decrease and fluid atoms tend to be trapped inside them. As a result, slip on the boundary diminishes as fluid atoms are trapped inside the cavities. 7 Improvement in MD method The MDS method, as a distinguished flow physics simulation method from other simulation methods, can provide direct numerical experiments at atomistic level with various physical conditions. However, MDS scheme, as a numerical method, adopts some hypotheses which induce some deviations of the simulation results from the real physical process. Furthermore, limited by the computational ability, MDS method can only be applied for nanoscale systems. Therefore, how to introduce the rational assumptions which reflect the true physical process and decrease the computational consumption is the main research direction in order to improve MDS method. 123 Microfluid Nanofluid (2010) 9:1011–1031 7.1 Thermal wall models for viscous dissipation effect Owing to the large velocity gradient, heat generated in nanochannels by viscous dissipation cannot be neglected. In the real physical process heat generated by viscous dissipation would be conducted through the wall so that the system can keep a constant temperature. Therefore, in MD simulation, there are several special methods must be adopted in the MDS to keep the constant temperature (Andersen 1980; Berendsen et al. 1984; Nose 1984, 2002; Hoover 1985). A relaxing method (Buckingham 1938) was developed by coupling an external bath at a constant temperature to the system and modifying the velocities at each iteration through a differential feedback mechanism. This method was widely used. However, it does not produce canonical probability distribution in momentum space during simulations although it can maintain constant temperature. Nose–Hoover thermostat (Errington and Panagiotopoulos 1999; Zhang and Ely 2004; Galliero et al. 2006) is the most popular method, which can sustain the canonical ensemble distribution both in configuration and momentum space. This method generates canonical temperature fluctuations by affecting intermolecular forces from extended ensemble system. However, because the wall molecules are fixed, the results which derived from Nose–Hoover model show almost unphysical uniform temperature distribution. Kim et al. (2008) developed a novel thermal wall model by which the nature of thermal transport at the wall–fluid interface could be considered. The main idea of this method is as follows. The wall molecules were in balance places with thermal oscillations. Thermal oscillations of wall molecules influence the fluid by absorbing or supplying momentum and energy to fluid molecules via intermolecular interactions. The resulting force on a molecule of the wall consists of two parts. The first part is the interaction between fluid molecules and wall molecules, which can be described by LJ potential. The second part is the interaction with the wall molecules, which is the lattice bond springs attached to the wall molecules at their lattice positions. The forces between molecules depend on interparticles distances. Hence, thermal oscillations of wall molecules affect the forces acting on the fluid molecules through the variation of their position as well. In this way, the dissipation heat could be conducted to the walls. Walls were utilized as heat baths to maintain thermal equilibrium of the fluid. Then, for interactive thermal wall model, velocity scaling thermostat to each wall layer separately qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Tassignedwall with parameter glayer ¼ Tlayerwall at every time step was applied to keep the walls temperature constant. Because only the fluid flow and thermal properties were concerned, the thermostat applied on the walls does not influence the Microfluid Nanofluid (2010) 9:1011–1031 results. The temperature profiles derived from this method were compared with that derived from the Nose–Hoover thermostat method. The results of this method show distinct differences. When the interaction potential is weak or the wall has a large crystal bonding stiffness, temperature jumps at the wall–fluid interface were observed. However, because one cannot apply the infinitely small time interval in simulation, there exist possibilities that the fluid molecules would penetrate into the solid walls, which violate the true physics. This article did not say how to deal with this situation. If no special method is applied to ensure the fluid molecules would not entry the walls domain, then the law of conservation of mass will not be satisfied. Jakobtorweihen et al. (2006) described a novel algorithm that includes the effect of host lattice flexibility into MD simulations that traditionally use rigid lattices. This method uses a Lowe-Andersen thermostat for interface– fluid collisions to consider the most important aspects of flexibility. This method was applied for different guest molecules, temperatures, and types of CNTs. It has a crucial influence on the diffusive dynamics of the guest molecules. In this article, a CNT is modeled as a flexible framestudy as described by Walther et al. (2003). Bond vibrations were modeled with a Morse potential. Bond bending was controlled by a harmonic cosine potential. And bond torsions with a twofold torsion potential. However, in this study no intermolecular carbon–carbon interactions were considered. Intermolecular interactions were modeled with the shifted and truncated LJ 12–6 potential. The fluid molecules were modeled as united atoms. In order to avoid flexible host matrix simulations and ensure that transport occurs under a fluid–solid thermal equilibrium, the authors used the so-called fluid–solid thermal diffuse scattering (TDS) algorithm. Simulations with this new algorithm lead to the same self-diffusion coefficients and give the same heat transfer across the rigid interface as across a flexible interface. The results showed that the flexibility of a smooth pore cannot be neglected in the low loading regime. The small vibrations of the carbon atoms have a significant influence on the self-diffusion at low loading in CNTs. The loading dependence of self-diffusion for helium in a (20.0) CNT at 300 K was also investigated. Owing to an increased relative roughness (smaller molecule) of the CNT surface, the flexibility influence is less than that in the case of methane. 7.2 Potentials The LJ 12–6 potential is by far the most widely used potential for the intermolecular interactions in MD simulations. It allows a relatively quick computation of the interaction in MD. Furthermore, it has been shown to be able to reveal most of the features experimentally found in 1027 fluid states though it formulation is simple. However, the LJ 12–6 potential is a many-body empirical potential, other than a true representation of two-body interactions between argon atoms (Stone 1996). The LJ 12–6 potential represents the decay of repulsive interaction by an inverse 12 power dependence on intermolecular separation distance whereas intermolecular repulsion decays exponentially. In addition to LJ 12–6 potential, the general three-parameter LJ a–6 potentials and three-parameter exponential a–6 potentials (Buckingham 1938) which uses an exponential formulation of the repulsive part of the potential, have shown the improved efficiency over the standard LJ potential for thermodynamic (Errington and Panagiotopoulos 1999) and even for dynamic properties such as viscosity (Zhang and Ely 2004). Galliero et al. (2006) tried to find more rational potential formulation. h    iA general form of LJ potential is proposed as: a b e rr þ rr ; the exponential values of a, b can be varied, for example from 12 to 22, to obtain better potential formulations. They computed the viscosity and the pressure of spherical fluid particles interacting through LJ a–6 and exponential a–6 potentials by using MD simulation. They intended to test the ability of the various potentials to reproduce simultaneously the pressure and the viscosity of the same simple real fluids in various thermodynamic states. They found that LJ potentials and exponential potentials both can lead to accurate results by using the appropriate exponent a. However, the exponential potentials results are not superior to the LJ ones. Furthermore, the authors found that the best results were obtained by the LJ 12–6 and the exponential 14-6 potentials for each type of the real fluid that they simulated. 7.3 Reduction of the computation consumption It is well known that the applications of MDS were limited mainly because of the required large computation power. Considering the LJ 12–6 potential, the interaction between each two molecules decays vary fast to zero when their distance gets larger than several times of the equilibrium separation distance. Therefore, cut-off radium, which is usually about three times of the equilibrium distance, is defined so that the neighboring atoms outside this radius will be ignored. The identification of neighboring atoms and the calculation of inter-atomic distance and force consume most of the calculation time. Therefore, how to select the rational cut-off radius in order to reduce the computational consumption has attracted many attentions. Several algorithms about select rational cut-off radius, such as Verlet list, cell-linked list and their combinations (Frenkel and Smit 2002), have been proposed. Furthermore, several improved algorithms of identifying the neighboring atoms have been proposed (Mason 2005; Yao 123 1028 et al. 2004; Heinz and Hunenberger 2004; Sutmann and Stegailov 2006; Maximova and Keasar 2006). However, these algorithms are either valid only for specific systems or lack of validation. Wang et al. (2007) developed an estimation formula of the computation time for each MD algorithm calculation so as to find an optimized performance for each algorithm. With the formula proposed in this article, the best algorithm can be chosen based on different total number of atoms, system average density and system average temperature for the MD simulation. The new contribution of this study is to reduce the computation time by define appropriate Verlet radius, the cut-off radius and the average number density of atoms in a system, respectively. In order to establish the list in recording the neighboring atoms around each atom more efficiently during the MDS process, a commonly applied skill is to use a large radius called, the Verlet radius, and the sphere with the Verlet radius is called the interaction sphere. This list will be updated for every K time-steps. As shown in Fig. 9, if the list recording the neighboring atoms is established by radius R and will be repeatedly used in the following k time-steps, then time would be defined so that after (k - 1)t steps some neighboring atoms, which are initially located on the sphere with radius R and may penetrate the sphere, can be contained by such a list in the beginning. The Verlet-cell list algorithm combines the advantages of Verlet list and cell-linked list. In the beginning of each listupdating interval, it requires N operations to assign each atom to a certain cell. It also requires operations to calculate and judge the inter-atomic distance between the central atoms and neighboring atoms in surrounding cells and operations to establish the neighbor list with neighboring atoms inside the interaction sphere. Using this algorithm, the reduction factor is the same as Verlet list, while t is much smaller. Fig. 9 Illustration of the cut-off radius Rc and Verlet radius Rv as well as the corresponding number densities of atoms N00 and N0 (Wang et al. 2007) 123 Microfluid Nanofluid (2010) 9:1011–1031 8 Summary In this article, we reviewed physics of liquid flows in nanochannels simulated by MD method, including the fundamentals of MD simulation, flows under different driving forces, the boundary conditions, and the improved algorisms of MD methods. It is worthwhile to comment on the key issues which may be the promising research directions in the studies of the liquid flows in nanochannels by MD simulation. In the research of electro-osmotic flow in nanochannels, several conclusions are different from the classical PB theory, such as the laying density profile near the walls, the zeta potential varied with the walls roughness, and the diffusion decreased by the influence of roughness. Furthermore, the ions distribution predicted by MDS departs from the Boltzmann distribution, which is the base of the PB theory. The MD simulations reveal that fluid molecules near the wall are positioned with the preferred orientation, which cannot be predicted by the PB theory but can explain the various thermodynamic properties of the liquids in nanochannels. These conclusions may not reflect exactly the real physics, because some simplifying assumptions are used in MDS models. Nevertheless, the results tell us that the PB equations at least should be modified when used to solve the nanoscale fluid flow problems. The MDS method is a powerful tool to probe the distributions and structures of molecules and ions in nanoscale, which are the dominate factors controlling the thermodynamic parameters. Therefore, in the electroosmotic flow research, studying the structure and distribution of molecules and ions in nanochannels will be an important direction. Carbon nanotubes with charges can be used to develop the next generation devices for pumping, sensing and particle screening applications. The study of the potential CNT applications and the liquid molecular structure in CNTs are promising directions. The previous studies showed that the water molecular would have a dipole structure in the charged nanotubes. CNTs with charges in the proper positions can be used as the liquid pump and MDS has been applied to prove its feasibility. In addition, the structure of water and ions inside the CNTs has significant effects on the effectiveness of particle screening and partition. Thermophoresis has practical implications in the fields such as the bore stability in the oil industry, the particle accumulation in a thermal boundary layer, and the transport of polymeric molecules in the nanochannels. The studies of thermophoresis by using MD simulations are still in the early stage. For the application of using thermophoresis as a pumping method, studies showed that the obstacle is the small driven force. Although a cascaded method was proposed to overcome the above difficulty, Microfluid Nanofluid (2010) 9:1011–1031 however, thermophoresis is not a promising direction to drive the liquid in nanochannels compared to the electrokinetically driven flow. For the improved MD scheme, the focus is how to reflect the realistic flow physics and use less computer power. The key issues involve building better potential functions, the thermal wall model, and the cut-off radius. The rational potential function is the most important issue in the MD simulation. By far, the LJ 12–6 potential is the most popular potential function, which has been verified by lots of studies. However, as an empirical formula, it can only cover limited types of molecules interactions. Different types of interactions and potential expressions are urgently needed. Quantum method starts from the sub-atom level structure to predict the properties of atom level. Therefore, quantum theory may be used to develop the molecular interaction potential functions. The thermal wall model becomes more important because of the large thermal dissipation in nanoscale. The simulation will fail if the thermal wall model is not appropriate. By far, the most successful thermostatic model is the Nose model. It can keep the system temperature constant. However, recent studies show that if the wall molecules are fixed in Nose, then it may induce the velocity and temperature profile deviation from the realistic physics status. Oscillation wall molecules’model was proposed in the new studies. The dissipation heat could be dispersed from the liquid to the oscillated wall molecules. Then,thermostatic is applied to the wall molecules. The result of this model is more close to the realistic process. Acknowledgements The study is supported by the National Natural Science Foundation of China with the contract number 50906088, and the National Natural Science Foundation of China with the contract number 50825603. The authors wish to thank the support of the Natural Sciences and Engineering Research Council through a research grant to D. Li. References Alexander FJ, Garcia AL (1997) The direct simulation Monte-Carlo method. Comput Simul 11:588–593 Andersen HC (1980) Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys V72(4):2384–2393 Banerjee S, Murad S, Puri IK (2007) Preferential ion and water intake using charged carbon nanotubes. Chem Phys Lett 434:292–296 Barrat JL, Bocquet L (1999) Large slip effect at a nonwetting fluid– solid interface. Phys Rev Lett 82:4671–4674 Benedict LX, Louie SG, Cohen ML (1995) Static polarizabilities of single-wall carbon nanotubes. Phys Rev B 52:8541–8549 Berendsen HJC, Postma JPM, van Gunsteren WF, Hermans J (1981) Interaction models for water in relation to protein hydration. In: Intermolecular forces. Reidel, Dordrecht, The Netherlands, pp 331–342 Berendsen HJC, Postma JPM, van Gunsteren WF, Di Nola A, Haak JR (1984) Molecular dynamics with coupling to an external bath. J Chem Phys 81(8):3684–3690 1029 Berendsen HJC, Grigera JR, Straatsma TP (1987) The missing term in effective pair potentials. J Phys Chem 91:6269–6271 Besteman K, Lee JO, Wiertz FGM, Heering HA, Dekker C (2003) Enzyme-coated carbon nanotubes as single-molecule biosensors. Nano Letters 3:727–730 Bikerman JJ (1970) Physical surfaces. Academic Press, New York Bonaccurso E, Butt HJ, Craig VSJ (2003) Surface roughness and hydrodynamic boundary slip of a Newtonian fluid in a completely wetting system. Phys Rev Lett 90:144501 Bourlon B, Wong J, Miko C, Forro L, Bockrath M (2007) A nanoscale probe for fluidic and ionic transport. Nat Nanotechnol 2:104–107 Branam RD (2005) Molecular dynamics simulation of supercritical fluids. PhD thesis, Pennsylvania State University, Pennsylvania Buckingham RA (1938) The classical equation of state of gaseous helium, neon and argon. Proc R Soc Lond Ser A 168:264–283 Chandrasekhar J, Spellmeyer DC, Jorgensen WL (1984) Energy component analysis for dilute aqueous solutions of Li?, Na?, F-, and Cl- ions. J Am Chem Soc 106:903–910 Chaudhury MK, Whitesides GM (1992) How to make water run uphill. Science 256:1539–1541 Chen SC, Hseih WC, Lin MF (2005) Charge screening of singlewalled carbon nanotubes in a uniform transverse electric field. Phys Rev B 72:193412 Choi CH, Westin KJA, Breuer KS (2003) Apparent slip flows in hydrophilic and hydrophobic microchannels. Phys Fluids 15:2897–2902 Churaev NV, Sobolev VD, Somov AN (1984) Slippage of liquid over hydrophobic solid-surfaces. J Colloid Interface Sci 97:574– 581 Cieplak M, Koplik J, Banavar JR (2001) Boundary conditions at a fluid–solid interface. Phys Rev Lett 86:803–806 Cottin-Bizonne C, Jurine S, Baudry J, Crassous J, Restagno F, Charlaix E (2002) Nanorheology: an investigation of the boundary condition at hydrophobic and hydrophilic interfaces. Eur Phys J E 9:47–53 Craig VSJ, Neto C, Williams DRM (2001) Shear-dependent boundary slip in an aqueous Newtonian liquid. Phys Rev Lett 87:054504 Daiguji H, Yang PD, Majumdar A (2004) Ion transport in nanofluidic channels. Nano Letters 4:137–142 Darden T, York D, Pedersen L (1993) Particle mesh Ewald—an N. Log(N) method for Ewald sums in large systems. J Chem Phys 98:10089 de Gennes PG, Brochard-Wyart E, Quere DF (2003) Capillarity and wetting phenomena: drops, bubbles, pearls, waves. Springer, New York Dukhin SS, Derjaguin BV (1974) Electrokinetic phenomena, surface and colloid science, vol 7. Wiley, New York Ermakov SV, Jacobson SC, Ramsey JM (1998) Computer simulations of electrokinetic transport in microfabricated channel structures. Anal Chem 70:4494–4504 Errington JR, Panagiotopoulos AZ (1999) A new intermolecular potential model for the n-alkane homologous series. J Phys Chem B 103:6314–6322 Fan R, Wu YY, Li DY, Yue M, Majumdar A, Yang PD (2003) Fabrication of silica nanotube arrays from vertical silicon nanowire templates. J Am Chem Soc 125:5254–5255 Frenkel D, Smit B (2002) Understanding molecular simulation: from algorithms to applications, 2nd edn. Academic Press, San Diego Freund JB (2002) Electro-osmosis in a nanometer-scale channel studied by atomistic simulation. J Chem Phys 116:2194–2200 Galea TM, Attard P (2004) Molecular dynamics study of the effect of atomic roughness on the slip length at the fluid–solid boundary during shear flow. Langmuir 20:3477–3482 Galliero G, Boned C, Baylaucq A, Montel F (2006) Molecular dynamics comparative study of Lennard–Jones a-6 and 123 1030 exponential a-6 potentials: application to real simple fluids (viscosity and pressure). Phys Rev E 73:061201 Ghosh S, Sood AK, Kumar N (2003) Carbon nanotube flow sensors. Science 299:1042–1044 Gogotsi Y, Libera JA, Guvenc-Yazicioglu A, Megaridis CM (2001) In situ multiphase fluid experiments in hydrothermal carbon nanotubes. Appl Phys Lett 79:1021 Gong XJ, Li JY, Lu HJ, Wan RZ, Li JC, Hu J, Fang HP (2007) A charge-driven molecular water pump. Nat Lett 2:709–712 Haile JM (1992) Molecular dynamics simulation: elementary methods. Interscience, New York Han M (2005) Thermophoresis in liquids: a molecular dynamics simulation study. J Colloid Interface Sci 284:339–348 Han M (2008) Thermally-driven nanoscale pump by molecular dynamics simulation. J Mech Sci Technol 22:157–165 Heinz TN, Hunenberger PH (2004) A fast pairlist-construction algorithm for molecular simulations under periodic boundary conditions. J Comput Chem 25:1474–1486 Holt JK, Park HG, Wang YM (2006) Fast mass transport through sub2-nanometer carbon nanotubes. Science 312:1034 Hoover WG (1985) Canonical dynamics: equilibrium phase-space distributions. Phys Rev A 31:1695–1697 Horn RG, Vinogradova OI, Mackay ME, Phan-Thien N (2000) Hydrodynamic slippage inferred from thin film drainage measurements in a solution of nonadsorbing polymer. J Chem Phys 112:6424–6433 Hummer G, Rasalah JG, Noworyta JP (2001) Water conduction through the hydrophobic channel of a carbon nanotube. Nature 414:188–190 Jabbarzadeh A, Atkinson JD, Tanner RI (1999) Wall slip in the molecular dynamics simulation of thin films of hexadecane. J Chem Phys 110:2612–2620 Jakobtorweihen S, Lowe CP, Keil FJ, Smit B (2006) A novel algorithm to model the influence of host lattice flexibility in molecular dynamics simulations: loading dependence of selfdiffusion in carbon nanotubes. J Chem Phys 124:154706-1-13 Joly L, Ybert C, Trizac E, Bocquet L (2004) Hydrodynamics within the electric double layer on slipping surfaces. Phys Rev Lett 93:257805 Joly L, Ybert C, Bocquet L (2006) Probing the nanohydrodynamics at liquid–solid interfaces using thermal motion. Phys Rev Lett 96:046101 Jorgensen WL, Bigot B, Chandrasekhar J (1982) Quantum and statistical mechanical studies of liquids. 21. The nature of dilute solutions of sodium and methoxide ions in methanol. J Am Chem Soc 104:4584–4591 Khare R, dePablo JJ, Yethiraj A (1996) Rheology of confined polymer melts. Macromolecules 29:7910–7918 Khare R, de Pablo J, Yethiraj A (1997) Molecular simulation and continuum mechanics study of simple fluids in nonisothermal planar coquette flows. J Chem Phys 107:2589–2596 Kim CS (2007) Nonequilibrium molecular dynamics approach for nanoelectromechanical systems: nanofluidics and its applications. J Fluids Eng 129:1140–1146 Kim D, Darve E (2006) Molecular dynamics simulation of electroosmotic flows in rough wall nanochannels. Phys Rev E 73:0512031-12 Kim BH, Beskok A, Cagin T (2008) Thermal interactions in nanoscale fluid flow: molecular dynamics simulations with solid–liquid interfaces. Microfluid Nanofluid 5:551–559 Koike A, Yoneya M (1998) Chain length effects on frictional behavior of confined ultrathin films of linear alkanes under shear. J Phys Chem B 102:3669–3675 Koplik J, Banavar JR, Willemsen JF (1989) Molecular-dynamics of fluid-flow at solid-surfaces. Phys Fluids A 1:781–794 123 Microfluid Nanofluid (2010) 9:1011–1031 Kuo TC, Sloan LA, Sweedler JV, Bohn PW (2001) Manipulating molecular transport through nanoporous membranes by control of electrokinetic flow: effect of surface charge density and debye length. Langmuir 17:6298–6303 Li TS, Lin MF (2006) Electronic properties of carbon nanotubes under external fields. Phys Rev B 73:075432 Li Y, Lu D, Rotkin SV, Schulten K, Ravaiolli U (2004) Electronic structure and dielectric behavior of finite-length single-walled carbon nanotubes. In: 4th IEEE conference on nanotechnology, pp 273–275 Linke H, Aleman BJ, Melling LD, Taormina MJ, Francis MJ, DowHygelund CC, Narayanan V, Taylor RP, Stout A (2006) Selfpropelled Leidenfrost droplets. Phys Rev Lett 96:154502 Liu L, Zhao JB, Yin CY, Culligan PJ, Chen X (2009) Mechanisms of water infiltration into conical hydrophobic nanopores. Phys Chem Chem Phys 11:6520–6524 Ludwig KF (2004) Molecular dynamics simulations of sub- and supercritical injection. Master thesis, Pennsylvania State University, Pennsylvania Mahar B, Laslau C, Yip R, Sun Y (2007) Development of carbon nanotube-based sensors—a review. IEEE Sens J 7:266–284 Mahoney MW, Jorgensen WL (2000) A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J Chem Phys 112:8910–8922 Martini A, Roxin A, Snurr RQ, Wang Q, Lichter S (2008) Molecular mechanisms of liquid slip. J Fluid Mech 600:257–269 Mashl RJ, Joseph S, Aluru NR, Jakobsson E (2003) Anomalously immobilized water: a new water phase induced by confinement in nanotubes. Nano Letters 3:589 Mason DR (2005) Faster neighbour list generation using a novel lattice vector representation. Comput Phys Commun 170:31–41 Maximova T, Keasar C (2006) A novel algorithm for non-bonded-list updating in molecular simulations. J Comput Biol 13:1041–1048 Maxwell JC (1879) On stresses in rarified gases arising from inequalities of temperature. Philos Trans R Soc Lond 170:231–256 Micci MM, Oechsle SK, Mayer WOH (2000) Molecular dynamics simulations of primary atomization. In: 8th international conference on liquid atomization and spray systems, Pasadena, CA Migler KB, Hervet H, Leger L (1993) Slip transition of a polymer melt under shear-stress. Phys Rev Lett 70:287–290 Mitchell MJ, Qiao R, Aluru NR (2000) Meshless analysis of steadystate electro-osmotic transport. J MEMS 9:435–449 Miyamoto S, Kollman PA (1992) Settle—an analytical version of the shake and rattle algorithm for rigid water models. J Comput Chem 13:952–962 Moseler M, Landman U (2000) Formation, stability, and breakup of nanojets. Science 289:1165–1169 Niavarani A, Priezjev NV (2008a) Rheological study of polymer flow past rough surfaces with slip boundary conditions. J Chem Phys 129:144902 Niavarani A, Priezjev NV (2008b) Slip boundary conditions for shear flow of polymer melts past atomically flat surfaces. Phys Rev E 77:041606 Nose S (1984) A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys 81:511–519 Nose S (2002) A molecular dynamics methods for simulations in the canonical ensemble. Mol Phys 100:191–198 Pennathur S, Santiago JG (2005) Electrokinetic transport in nanochannels. 2. Experiments. Anal Chem 77:6782–6789 Pit R, Hervet H, Leger L (2000) Direct experimental evidence of slip in hexadecane: solid interfaces. Phys Rev Lett 85:980–983 Plimpton S (1995) Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 117:1–19 Priezjev NV (2007) Rate-dependent slip boundary conditions for simple fluids. Phys Rev E 75:051605 Microfluid Nanofluid (2010) 9:1011–1031 Priezjev NV, Troian SM (2004) Molecular origin and dynamic behavior of slip in sheared polymer films. Phys Rev Lett 92:018302 Priezjev N, Troian SM (2006) Influence of periodic wall roughness on the slip behaviour at liquid/solid interfaces: molecular-scale simulations versus continuum predictions. J Fluid Mech 554:25–46 Qiao R, Aluru NR (2003) Ion concentrations and velocity profiles in nanochannel electroosmotic flows. J Chem Phys 118:4692–4701 Qiao R, Aluru NR (2005) Atomistic simulation of KCI transport in charged silicon nanochannels: interfacial effects. Colloids Surf A 267:103–109 Qiao Y, Cao GX, Chen X (2007) Effects of gas molecules on nanofluidic behaviors. J Am Chem Soc 129:2355–2359 Qiao Y, Liu L, Chen X (2009) Pressurized liquid in nanopores: a modified Laplace-Yong equation. Nano Letters 9:984–988 Raviv U, Laurat P, Klein J (2001) Fluidity of water confined to subnanometre films. Nature 413:51–54 Regan BC, Aloni S, Ritchie RO, Dahmen U, Zettl A (2004) Carbon nanotubes as nanoscale mass conveyors. Nature 428:924–927 Saito R, Fujita M, Dresselhaus G, Dresselhaus MS (1992) Electronicstructure of chiral graphene tubules. Appl Phys Lett 60: 2204–2206 Sanchez-Reyes J, Archer LA (2003) Interfacial slip violations in polymer solutions: role of microscale surface roughness. Langmuir 19:3304–3312 Schmatko T, Hervet H, Leger L (2005) Friction and slip at simple fluid–solid interfaces: the roles of the molecular shape and the solid–liquid interaction. Phys Rev Lett 94:244501 Schmatko T, Hervet H, Leger L (2006) Effect of nanometric-scale roughness on slip at the wall of simple fluids. Langmuir 22:6843–6850 Service RF (2006) Desalination freshens up. Science 313:1088–1090 Shin H (2006) Non-equilibrium molecular dynamics simulations of nanojet injection. PhD thesis, Yonsei University, Seoul, Korea Sofos F, Karakasidis T, Liakopoulos A (2009) Transport properties of liquid argon in krypton nanochannels: anisotropy and nonhomogeneity introduced by the solid walls. Int J Heat Mass Transf 52:735–743 Soong CY, Yen TH, Tzeng PY (2007) Molecular dynamics simulation of nanochannel flows with effects of wall lattice–fluid interactions. Phys Rev E 76:036303-1-14 Spoel DVD, Lindahl E, Hess B, Buuren ARV, Apol E, Meulenhoff PJ, Tieleman DP, Sijbers ALTM, Feenstra KA, Drunen RV et al (2004) Gromacs user manual, version 3.2. University of Groningen, Groningen. www.gromacs.org Stein D, Kruithof M, Dekker C (2004) Surface-charge-governed ion transport in nanofluidic channels. Phys Rev Lett 93:035901 Stern MB, Geis MW, Curtin JE (1997) Nanochannel fabrication for chemical sensors. J Vac Sci Technol 15:2887–2891 Stone AJ (1996) The theory of intermolecular forces. Clarendon Press, Oxford Sun H, Ren P, Fried JR (1998) The COMPASS force field parameterization and validation for phosphazenes. Comput Theor Polym Sci 8:229–246 Sutmann G, Stegailov V (2006) Optimization of neighbor list techniques in liquid matter simulations. J Mol Liq 125:197–203 1031 Tas NR, Haneveld J, Jansen HV, Elwenspoek M, Berg AVD (2004) Capillary filling speed of water in nanochannels. Appl Phys Lett 85:3274–3726 Tenenbaum A, Ciccotti G, Gallico R (1982) Stationary nonequilibrium states by molecular dynamics—Fourier’s law. Phys Rev A 25:2778–2787 Thomas JA, McGaughey AJH (2009) Water flow in carbon nanotubes: transition to subcontinuum transport. Phys Rev Lett 102:184502 Thompson PA, Robbins MO (1990) Shear-flow near solids: epitaxial order and flow boundary-conditions. Phys Rev A 41:6830–6837 Thompson PA, Troian SM (1997) A general boundary condition for liquid flow at solid surfaces. Nature 389:360–362 Thompson PA, Robbins MO, Grest GS (1995) Structure and shear response in nanometer-thick films. Isr J Chem 35:93–106 Turner CH, Brennan JK, Pikunic J, Gubbins KE (2002) Simulation of chemical reaction equilibria and kinetics in heterogeneous carbon micropores. Appl Surf Sci 196:366–374 Tyrell J, Attard P (2001) Images of nanobubbles on hydrophobic surfaces and their interactions. Phys Rev Lett 87:176104 Vinogradova OI, Yakubov GE (2006) Surface roughness and hydrodynamic boundary conditions. Phys Rev E 73:045302 Wang DB, Hsiao FB, Chuang CH, Lee YC (2007) Algorithm optimization in molecular dynamics simulation. Comput Phys Comm 177:551–559 Werder T, Walther JH, Jaffe RL, Halicioglu T, Koumoutsakos P (2003) On the water–carbon interaction for use in molecular dynamics simulations of graphite and carbon nanotubes. J Phys Chem B 107:1345–1352 Wheeler TD, Stroock AD (2008) The transpiration of water at negative pressures in a synthetic tree. Nature 455:208–212 Whitby M, Quirk N (2007) Fluid flow in carbon nanotubes and nanopipes. Nature Nanotech 2:87–94 Whitesides GM (2006) The origins and the future of microfluidics. Nature 442:368–373 Xu Y, Aluru NR (2008) Carbon nanotube screening effects on the water–ion channels. Appl Phys Lett 93:043122-1-3 Xue YQ, Chen MD (2006) Dynamics of molecules translocating through carbon nanotubes as nanofluidic channels. Nanotechnology 17:5216–5223 Yang SC (2006) Effects of surface roughness and interface wettability on nanoscale flow in a nanochannel. Microfluid Nanofluid 2:501–511 Yang L (2007) Modeling the selective partitioning of cations into negatively charged nanopores in water. J Chem Phys 126: 084706-1-8 Yao ZH, Wang HS, Liu GR, Cheng M (2004) Improved neighbor list algorithm in molecular simulations using cell decomposition and data sorting method. Comput Phys Commun 161:27–35 Zhang HZ, Ely JF (2004) AUA model NEMD and EMD simulations of the shear viscosity of alkane and alcohol systems. Fluid Phase Equilib 217:111–118 Zhang WF, Li DQ (2007) Simulation of flow speed 3D nanochannel flow. Microfluid Nanofluid 3:417–425 123