Equilibrium Simulation of the Slip Coefficient in Nanoscale Pores
aa r X i v : . [ c ond - m a t . s o f t ] M a y Submitted to Physical Review Letters
Equilibrium Simulation of the Slip Coefficient in Nanoscale Pores
Vlad P. Sokhan ∗ National Physical Laboratory, Hampton Rd,Teddington TW11 0LW, United Kingdom
Nicholas Quirke † Department of Chemistry, Imperial College London, London NW7 2AY, United Kingdom (Dated: October 30, 2018)
Abstract
Accurate prediction of interfacial slip in nanoscale channels is required by many microfluidic applica-tions. Existing hydrodynamic solutions based on Maxwellian boundary conditions include an empiricalparameter that depends on material properties and pore dimensions. This paper presents a derivation of anew expression for the slip coefficient that is not based on the assumptions concerning the details of solid-fluid collisions and whose parameters are obtainable from equilibrium simulation. The results for the slipcoefficient and flow rates are in good agreement with non-equilibrium molecular dynamics simulation.
PACS numbers: PACS numbers: 47.15.gm, 47.11.Mn, 05.20.Jj, 47.90.+a, 66.20.+d, 68.08.–pc (cid:13)
Crown copyright, 2007 Kn = λ/L, is non-negligible, andis usually associated with low densities where λ is large. However, in narrow capillaries (where L is small), slip can be observed even at liquid densities. In the general case it is characterizedby a slip coefficient, l s [7, 8], which in the absence of temperature gradients relates the collectivemolecular velocity at the wall, the slip velocity, to the shear rate, u s = l s ∇ u [7]. Using kinetictheory Maxwell derived a microscopic expression for the slip coefficient [7], that can be written as l s = λ (cid:18) α − (cid:19) , (1)where λ = 2 η/ρ ¯ c is the mean free path, η is the shear viscosity, ρ is mass density of the fluid,and ¯ c is the mean speed of the molecules. Considering only two types of wall collision, specularand diffuse (Knudsen) reflection, he introduced a coefficient α that defines a fraction of specularlyreflected molecules. In more broad sense, α defines the fraction of the flux of tangential momen-tum transmitted in collisions and is often called the ‘accommodation coefficient’ (TMAC) [8]. Itsvalue is defined by the details of the solid-fluid interactions and (1) implies finite slip even forpurely diffusive reflections ( α = 1 ). For specularly reflecting surface the slip coefficient divergessince the fluid cannot grip onto the surface.Slip in nanoscale fluid flow depend on many parameters inculding surface roughness [6], elec-tric properties of the interface [9], wetting conditions [5], chemical patterning of the surface [10],and is a nonlinear function of the dynamic state [11]. Application of the generalized Navier-Stokes(NS) hydrodynamics to problems of fluid flow on the nanoscale is an attractive but highly non-trivial task even in simple cases such as plane Poiseuille flow. An accurate solution requires aknowledge of the material parameters of the fluid as a function of local density which deviatesfrom its bulk value in the proximity of the interface. As a result, the velocity profile in narrowpores also deviates from the macroscopic prediction [12]. The usual approach is to regard thefluid as incompressible and to replace the complex non-uniform flow problem plus simple no-slipboundary conditions with a simple flow problem but with boundary conditions that has been called‘exceedingly difficult’ for theoretical investigation [8], p. 96. The discontinuity in the flow fieldis introduced in this approach in the same way the surface excess was introduced by Gibbs inhis treatment of interface boundaries at equilibrium as illustrated in Fig. 1 for the case of plane2 IG. 1: Full line, velocity profile for the plane Poiseuille flow with surface slip; dashed line, extrapolatedvelocity profile. Vertical gray line marks the position of the solid. Dash-dotted line denotes the velocitygradient at the wall. Also shown are slip velocity, u s , and slip coefficient, l s . Poiseuille flow in z direction ( u = { , , u z } ) between walls at y = ± h , where the symmetry ofthe solution was taken into account and only one half was plotted. We emphasize that the finiteslip u s that appears in this approach is a purely artificial device introduced only to match the ap-proximate solution with the solution of full problem in the middle part of the channel. The fullsolution for the continuum velocity field does decay to zero at the wall, as required by continuityof stresses in classical hydrodynamics [13]. It should be noted however that the limiting value ofthe velocity field cannot be observed in any experiment or in a statistical particle-based simula-tion with a continuous solid-fluid potential since the velocity field, which is defined everywherethe fluid density is non-zero, cannot be measured from particle velocities in regions where par-ticles are not observed due to finite sampling of the statistical ensemble, i. e. where the particleBoltzmann factor is diminishingly small but non-zero.There were a number of attempts to estimate the TMAC from kinetic theory [14, 15] and using3olecular dynamics simulation of simple liquids [5, 16, 17]. However, the results obtained bydifferent groups and using different methods are inconcluisive [1], with large scatter in obtainedvalues of slip coefficient spanning two orders of magnitude. The agreement between recent ex-periments with various surfaces [18, 19, 20, 21], and and theory is also poor, with experimentalvalues for α generally higher than the corresponding theoretical expectations [15, 16, 18]. One ofthe limitations of the original Maxwell model is that it was developed to solve the half-space flow,or Kramers, problem [8] in the dilute gas regime and faces problems describing flow in pores offinite width. It is assumed that the TMAC is a local property independent of the pore width. Thefollowing argument however demonstrates why it could not be so. Viscous wall stress in planePoiseuille flow that is proportional to the velocity gradient at the wall scales as the pore widthfor the flow of gas of the same density. This stress is transmitted to the wall through collisionsthat are independent of the pore width and therefore the momentum transferred to the wall in eachcollision and consequently the TMAC depends on the pore width.In this paper we present a method of calculating the slip coefficient from equilibrium simulationthat does not require assumption about the type of wall collisions avoiding thus the necessity ofcalculating the TMAC. We illustrate the method on a simple case of gravity driven plane Poiseuilleflow in a pore of width H = 2 h and with acceleration due to an external field g = { , , g } assketched in Fig. 1. The method can be extended to a more general case of Poiseuille flow inducedby a pressure gradient using the equivalence between the pressure gradient and gravity-drivenforce in the direction of flow [13]. It is convenient to reformulate the slip velocity problem (theDirichlet boundary condition) in terms of interfacial viscosity (Neumann boundary condition), andthe lateral wall stress. This idea goes back to Navier [22] who obtained the boundary condition forthe velocity field on the basis of particle arguments (cf. last eqn. on p. 415 in Ref. [22]) as, η ∂u z ( y ) ∂y (cid:12)(cid:12)(cid:12)(cid:12) y w = βu z ( y w ) (2)where his parameter β is related to the interfacial shear viscosity η via l s = η/β. Note that for non-linear velocity profiles, as for Poiseuille flow, the slip coefficient, l s , (shown in figure) is differentfrom the slip length, defined as the distance from the wall where the extrapolated velocity profilevanishes. However, if the difference between them is small, simple geometric consideration allowsus to establish the relationship between the two types of boundary condition.Neglecting viscous dissipation, the wall shear stress in terms of the external driving force actingon fluid particles, which in this case is simply σ yz = ρgh [23], can be equated with the Stokesian4rag force per unit area exerted on the wall by the moving fluid [24], σ yz = − F S A = M u Aτ ≡ ρuhτ , (3)where M is the total mass of the fluid in the pore; u ≡ h − R h u z ( y ) dy is the mean fluidvelocity, and the relaxation time τ can be calculated in molecular dynamics simulation fromthe Langevin equation for the fluid subsystem considered as a single Brownian particle usingfluctuation-dissipation theorem. It has been shown recently [24] that the fluid velocity autocorre-lation function decays exponentially, C ( t ) = M − kT exp ( − t/τ ) . (4)This provides a simple way to determine τ in an equilibrium simulation by fitting a one-parameterexponential to velocity autocorrelation data. From the two expressions for the wall shear stress, asimple relationship between the fluid velocity and the acceleration due to the external force, u = τ g. (5)can be established. This surprisingly simple result shows that within the limits of linear regime[11] the rate of fluid flow in non-equilibrium steady state can be estimated from the characteristictime of the decay of fluctuations at equilibrium. Using it, we can also estimate both the slipvelocity and slip coefficient. Since the hydrodynamic solution in this case is given by a quadraticvelocity profile with slip, u z ( y ) = ρg η (cid:0) h − y (cid:1) + u s , (6)using (5) and the definition of u , one obtains for slip velocity u s = (cid:18) τ − ρh η (cid:19) g, (7)and by substituting it into the definition of the slip coefficient, l s = − ∂u/∂y | h u s , one obtainsfinally l s = τ ηρh − h . (8)This is the main result. It shows that the slip coefficient is independent of the external force (flux),but nonlinearly depends on the pore width, both directly and indirectly through the relaxation time τ ≡ τ ( h ) . The connection with the Maxwell’s result can be established by using the relationship5etween the TMAC and relaxation time, f , τ = ( f α ) − [24]. Taking the kinetic theory expres-sions for the wall collision frequency per particle, f = ¯ c/ h and viscosity, obtain the expressionfor the slip coefficient l s = λ α − h (9)that can be directly compared with Maxwell’s result (1). We note here that usage of kinetic theoryexpression for f introduces a significant error at liquid densities, and it was used in deriving (9)only to make the same level of approximation as was in deriving (1).In order to test the accuracy of (8), two series of equilibrium molecular dynamics runs were per-formed. There are two ways of controlling the Knudsen number in simulation: either by changingthe mean density of the system (the mean free path) or the pore width (the characteristic length),and we used both. In the first series, the relaxation time τ was calculated as a function of fluiddensity for a pore of fixed width. In the second, the pore width was varied while the normal wallpressure was kept constant. The slip coefficient estimated using (8) was compared with the valuesobtained directly in the parallel set of nonequilibrium molecular dynamics (NEMD) simulationsof the Poiseuille flow for the same systems. The system consists of supercritical Lennard-Jones(LJ) fluid, confined between two walls modeled by a rigid triangular lattice of atoms situated at y w = ± h, and periodically replicated along x and z axes to avoid edge effects. All interactionsin the system were of the LJ form, U ( r ) = 4 ε (cid:2) ( σ/r ) − ( σ/r ) (cid:3) , where ε and σ are the usualenergy and length parameters. Their values for the fluid-fluid interactions define the correspondingscales, and in the following reduced units [25] are used, denoted by the astersk. The solid-fluidinteraction parameters in these units were taken as ε sf = 0 . ε and σ sf = 0 . σ, which isappropriate for methane gas between carbon ‘rare-gas walls’. The surface number density of thesolid was n s σ = 1 . as in our earlier study [24] where more details about the system andthe numerical scheme can be found. All molecular dynamics calculations were performed usingthe classical molecular dynamics software package MDL [26]. The nonequilibrium steady-stateconditions were realized by placing the fluid in a uniform external field parallel to the walls andcoupling all fluid degrees of freedom to a Nos´e-Hoover thermostat at T = 2 . εk − , where k is the Boltzmann constant. The simulation cell was of dimensions . σ × H × L z , where thedimension in the flow direction, L z , was scaled with the density to keep the number of particles inthe system around 2500. The timestep was ∆ t ∗ = 7 . × − in LJ units and the integration timein each case was not less than . × (50M steps). The value of the acceleration due to theexternal force was varied between g ∗ = 4 · − and g ∗ = 0 . in order to keep the fluid velocity6elow u ∗ = 0 . . The steady states reported here are known to be well within the Newtonian regime[27] and the slip velocity in the linear regime [11]. To simplify the comparison we scaled all flowrelated properties to a common value g ∗ = 4 . · − .In the first series three pore widths were considered, H = 5 . σ, . σ, and . σ ; (inthe future denoted as σ , σ , and σ for brevity) and several number densities ranging from n ∗ = 0 . ( n ∗ = m − ρσ ) , corresponding to the rarefied fluid at Kn = 2 . for a narrow pore(where we used the kinetic theory expression for the mean free path and the hard sphere diameterfor LJ [28]), to a dense state of n ∗ = 0 . that for a wide pore gives Kn = 0 . , thus spanningmore than two orders of magnitude of Knudsen numbers. In all cases the Reynolds number definedas the ratio of inertial and viscous forces, Re = ρuH/η, was kept small, Re < .In order to calculate slip velocity using (8) we need the estimates of the shear viscosity in thechannel. It can be calculated in equilibrium moleuclar dynamics e.g. from the stress autocorrela-tion function using Green–Kubo relations [29]. For the bulk LJ fluid it has been recently accuratelyestimated in molecular dynamics simulation [29]. To simplify the calculations we used an empiri-cal equation of state for the viscosity [30] at a density equal to the mean density in the central partof the channel, where it is uniform to fit the simulated data of [29] at the required temperature.The accuracy of this procedure was estimated by comparing the obtained values with values cal-culated from the NS hydrodynamics using the quadratic fit to NEMD velocity profiles. The resultsfor three pore widths are presented in Fig. 2. The agreement for the two wider pores is excelent.Deviation from the bulk values for the third pore, H = 5 σ , is due to the overlap of adsorbed layersat two surfaces and as a consequence, to inaccuracy in determination of the corresponding bulkdensity. At low densities the viscosity markedly deviates from bulk values for all three pores when Kn ≥ .The ratio of the slip coefficients to the pore width, calculated in NEMD and estimated from therelaxation times and bulk viscosities at densities equal to that in the middle of pore using (8) arecompared on Fig. 3 for three pore widths. The relaxation time was estimated from the exponentialfit to the collective velocity autocorrelation function calculated in equilibrium molecular dynamicssimulaton using the procedure described in [24]. The statistical uncertainty in all cases is of theorder of symbol size and is slightly higher for wider pores since the longer relaxation times inthis case require more accurate estimation of the velocity autocorrelation time at long correlationtimes. For all pores the slip coeficients calculated using the two routes agree within statisticaluncertainties. For a given pore, both (1) and (9) predict linear dependence of the slip length on7 IG. 2: (color online). Comparison of the calculated shear viscosities for three pore widths as a function ofreduced number density in the centre of pore, symbols, with empirical EOS for bulk fluid, dashed line. the Knudsen number with zero and negative offset, correspondingly. For two wider pores thebehaviour is observed at high Knudsen numbers,
Kn > , but with positive offset that increaseswith pore width. Surprisingly, for the narrowest pore, the slip length appears to be independent ofthe Knudsen number (density) for Kn > . Since density enters only the first term in (8), this resltindicate that in the low density regime the relxation time in narrow pores increases linearly withdensity.The results also show that for the systems with the same Knudsen number the slip lengthincreases with the pore width. In order to establish whether there is a limiting value of the slipcoefficient a second series of calculations was performed at the normal pressure that correspondsto reduced density of n ∗ = 0 . up to the pore width H = 210 σ. Obtained values of the reducedslip coefficient are compared on Fig. 4 with the Maxwell’s results using (1). The results estimatedusing (9) agree with those obtained directly in NEMD within statistical uncertainty. They showthat at about H = 25 σ the slip coefficient reaches the limiting value l s = 4 . σ . Maxwell’s8 IG. 3: (color online). Slip coefficients calculated directly in NEMD (open symbols and full lines) andestimated from the relaxation time (eqn. (8)), filled symbols and dashed lines) as a function of Knudsennumber for three pore widths. Lines are drawn to guide the eye. theory, on the contrary, predicts linear scaling of the slip length with pore width and for the widestpore studied it gives the value l s = 77 . σ (note the difference in scales).This work was funded as part of the National Physical Laboratory’s Strategic Research Pro-gramme and EPSRC under grant GR/N64809/01. ∗ Electronic address: [email protected] † Also at Department of Physics, University College Dublin, Ireland[1] E. Lauga, M. P. Brenner, and H. A. Stone, in
Handbook of Experimental Fluid Dynamics , edited byC. Tropea, A. Yarin, and J. F. Foss (Springer, N.Y., 2007), Ch. 19.[2] C. Neto, D. R. Evans, E. Bonaccurso, H.-J. Butt, and V. S. J. Craig, Rep. Prog. Phys. , 2859 (2005). IG. 4: (color online). Slip coefficients (left scale) calculated directly in NEMD, triangles, and estimatedfrom the relaxation time (eqn. (8)), circles, as a function of pore width. For comparison, slip coefficientsestimated from the Maxwell’s theory, eqn. (1), diamonds, are also shown (right scale). Lines are drawn toguide the eye.[3] J.-L. Barrat and L. Bocquet, Faraday Discuss. , 119 (1999).[4] Y. Zhu and S. Granick, Phys. Rev. Lett. , 096105 (2001).[5] V. P. Sokhan, D. Nicholson, and N. Quirke, J. Chem. Phys. , 3878 (2001).[6] C. Cottin-Bizonne, J.-L. Barrat, L. Bocquet, and E. Charlaix, Nature Mater. , 237 (2003).[7] J. C. Maxwell, Philos. Trans. R. Soc. London , 231 (1879).[8] C. Cercignani, Mathematical Methods in Kinetic Theory (Plenum Press, NY / Macmillan, London,1969)[9] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, J. Chem. Phys. , 204716 (2006).[10] T. Qian, X.-P. Wang, and P. Sheng, Phys. Rev. E , 022501 (2005).[11] P. A. Thompson and S. M. Troian, Nature , 360 (1997).
12] K. P. Travis and K. E. Gubbins, J. Chem. Phys. , 1984 (2000).[13] G. K. Batchelor,
An Introduction to Fluid Dynamics (Cambridge University Press, Cambridge 1970).[14] D. Nicholson and S. K. Bhatia, J. Memb. Sci. , 244 (2006).[15] H. Ambaye and J. R. Manson, Phys. Rev. E , 031202 (2006).[16] V. P. Sokhan, D. Nicholson, and N. Quirke, J. Chem. Phys. , 8531 (2002).[17] G. Arya, H. C. Chang, and E. J. Maginn, Molec. Simul. , 697 (2003).[18] H. Legge, J. R. Manson, and J. P. Toennies, J. Chem. Phys. , 8767 (1999).[19] E. B. Arkilic, K. S. Breuer, and M. A. Schmidt, J. Fluid Mech. , 29 (2001).[20] O. V. Sazhin, S. F. Borisov, and F. Sharipov, J. Vac. Sci. Technol. A , 2499 (2001).[21] B.-Y. Cao, M. Chen, and Z.-Y. Guo, Appl. Phys. Lett. , 091905 (2005).[22] C. L. M. H. Navier, M´em. Acad. Sci. Inst. France , 389 (1827).[23] V. P. Sokhan, N. Quirke, and J. Greenwood, Molec. Simul. , 535 (2005).[24] V. P. Sokhan and N. Quirke, Molec. Simul. , 217 (2004).[25] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford 1987).[26] Molecular Dynamics Laboratory (MDL) is a classical molecular dynamics package for molecularsystems developed by V. Sokhan at NPL within the Strategic Research Fellowship Programme.[27] J.-L. Barrat and L. Bocquet, Phys. Rev. Lett. , 4671 (1999).[28] D. Heyes, J. Chem. Soc. Far. 2 , 705 (1988).[29] K. Meier, A. Laesecke, and S. Kabelac, J. Chem. Phys. , 3671 (2004).[30] L. V. Woodcock, AICHE J. , 438 (2005)., 438 (2005).