Slip length of confined liquid with small roughness of solid-liquid interfaces Li Wan∗ and Yunmi Huang Department of Physics, Wenzhou University, Wenzhou 325035, P. R. China (Dated: January 4, 2017) We havestudied theslip length of confined liquid with small roughness of solid-liquid interfaces. DyadicGreenfunctionandperturbationexpansionhavebeenappliedtogetthesliplengthquantita- tively. Inthesliplength,botheffectsoftheroughnessoftheinterfacesandthechemicalinteraction between the liquid and the solid surface are involved. For the numerical calculation, Monte Carlo method has been used to simulate the rough interfaces and the physical quantities are obtained 7 statistically over the interfaces. Results show that the total slip length of the system is linearly 1 proportional to the slip length contributed from the chemical interaction. And the roughness of 0 the interfaces plays its role as the proportionality factor. For the roughness, the variance of the 2 roughness decreases the total slip length while thecorrelation length of theroughness can enhance theslip length dramatically to a saturation value. n a PACS numbers: 47.10.ad, 47.15.-x, 47.61.-k, 87.85.Rs,47.63.mf J Keywords: Slip length, confined liquid, rough interface, correlation length 3 ] I. INTRODUCTION And for the interfaces with small roughness, the rough n interfaces can be treated as flat ones with the roughness y d The dynamics of liquid confined in nano structures involved in the slip length of the flat interfaces. So, for - have attracted interests due to their potential appli- the latter, a question arises. What is the slip length u cations in nanofluidics [1–6]. It has been considered b of the interfaces with the roughness involved? In l f that the Navier-Stokes(NS) equation is still valid for this study, we will focus on the interfaces with small . s the description of the fluid hydrodynamics even if the roughness and address the problem to get the b relating c i liquid is scaleddownto a few molecularlayers[7–14]. In to the roughness in a quantitative way. s order to solve the NS equation for the confined liquid, y h a proper boundary conditions(BC) should be applied The slip length b is defined as the ratio of the shear p at the solid-liquid interfaces. A commonly used BC viscosity η of liquid to the friction coefficient κ of the [ is the no-slip BC, which requires that the fluid has solid surface by b = η/κ. It is known that the η is zero velocity at the interfaces. Such BC increases the position dependent in the liquid [7, 11, 16, 37, 38]. Es- 1 v hydrodynamic resistance to the fluid when the size of pecially, the viscosity of liquid in the region close to the 9 the confined liquid is decreased, which brings difficulties interfaces may be enhanced compared to the viscosity 6 for the applications of the nanofluidics [9]. In order in the region away from the interfaces. The discrepancy 5 to overcome the difficulties, the slippage of liquid at of the viscosities in the two regions leads to the shift of 0 the solid surfaces has been introduced to decrease the the exactlocationofthe BC appliedawayfromthe solid 0 hydrodynamic resistance [2, 9]. The BC with the surfaces by a few molecular layers [16, 39]. However, . 1 slippage of the liquid at the interfaces is known as the in this study, we consider that the fluid is Newtonian 0 Navier boundary condition (NBC). The NBC requires and the η is constant in the whole computational 7 that the velocity components of the liquid does not domain for simplicity, since we only focus on the ef- 1 vanishattheinterfacesinthe movingdirection. Inorder fectoftheroughnessoftheinterfacesonthesliplengthb. : v to quantify the slippage of liquid, a slip length b has i been introduced in the NBC [15–24]. Generally, the X The κ is used to show the resistance of the solid no-slip BC can be recovered from the NBC by setting r surfaces to the liquid. Generally, two mechanisms a b = 0. And, larger the slip length b, more slippery the are contributed to the κ. One mechanism is of the liquid at the interfaces. Since the slip length b used in chemical interaction between the liquid molecules and the NBC determines the solutions to the NS equation the solid molecules at the solid-liquid interfaces, like for the confined liquid, many authors have proposed that the attractive or repulsive force to aqueous liquids various structures for the interfaces to enhance the slip depends on whether the solid surfaces are hydrophobic length b [25–36]. But for the various structures, the or hydrophilic respectively. The other mechanism is fundamental structure is the roughness, which is the due to the fine structures of the interfaces, such as the nature of the interfaces. For the interfaces with large roughness of the interfaces. The collusions between roughness, one has to solve the NS equation by using the liquid molecules and the rough interfaces bring mathematical tools such as the finite element method. the friction force to the liquid. We denote the friction coefficient contributed from the chemical interaction by κ and the contribution of roughness of the interfaces w ∗ [email protected] by κ in the following. For an ideally flat interface, the r 2 chemical interaction contributes totally to the κ with A. Fluid Velocity κ = κ and κ = 0. However, for a real interface, the w r roughness is inevitable and κr is nonzero. In such a ThecompleteNavier-Stokes(NS)equationforthefluid real case, the friction forces respecting the two mecha- in the channel has the form of nisms then are obtained by f = κ v and f = κ v w w s r r s d~v 4η with v the fluid velocity at the interfaces in the η∇×∇×~v+ρ =−∇p+ρ~g−ρ~v·∇~v+ ∇(∇·~v). (1) s dt 3 moving direction. Due to the additivity of the friction forces,wegetthetotalfrictioncoefficientasκ=κ +κ . Here,~v isthe fluidvelocityandρthe liquiddensity. The w r ρandthe η havebeenassumedtobe uniformforthe liq- uidin the channel, since we focus onthe effects fromthe The determination of friction force of the solid sur- solid surfaces only. In the right hand side of the eq.(1), faces, hence the b, for the confined liquid has been a the first two terms represent the force of the gradient long-standing question and has been heavily studied by of pressure and the gravitational force respectively. The manyauthors[7,11,18,21,25–36,40–43]. Experimental third term is for the convectand the final term is due to and theoretical results, as well as the Molecular Dynam- the compressibilityofthe fluid. Inthe eq.(1), the equal- ics Simulations, have been reported. In those studies, ity of ∇2~v =∇(∇·~v)−∇×∇×~v has been applied. We the roughness of the solid surfaces causing the no-slip think the roughness of the interfaces is not too large to BC has been pointed out and later confirmed by the bringtheconvecteffectinthechannelandthefluidisin- Molecular Dynamics Simulations [18, 35, 36]. However, compressible. Thus,wedropoffthelasttwotermsonthe theroleoftheinterfaceroughnessintheslipBCactually righthandsideofeq.(1). Andforsimplicity,weusef~to is still lack. A very interesting result basing on the represent the two force terms by f~= −∇p+ρ~g. Based linear response theory has been obtained, in which the on the above considerations, the complete NS equation structure factor of the liquid in the region close to the is reduced to be interfaces is involved in the κ [9, 16]. However, the ρd~v 1 structure factor does not have an explicit relation to ∇×∇×~v+ = f~. (2) η dt η the roughness of the interfaces in that study. In order to complete the understanding of the flow boundary Now we apply the time Fourier transformations of conditions for confined liquid, it is necessary to bridge U~(ω) = +∞~veiωtdt and F~(ω) = +∞f~eiωtdt on the −∞ −∞ the slip length b to the roughness of the solid surfaces, both sideRs of the eq.(2) and transfoRrm the equation in which is our goal in this study. the frequency domain to be ρ 1 ∇×∇×U~ −i ωU~ = F~. (3) η η ↔ We introduce a dyadic Green function G(~r,~r′,ω), satis- fying the following equation ↔ ρ ↔ ↔ ∇×∇×G−i ωG= Iδ(~r−~r′), (4) η II. THEORY ↔ with I the dyadic unit and ~r(~r′) the position vector in the channel. We consider the equations of one fixed fre- We consider a model that simple liquid is confined by quency and drop off the notation of ω for convenience. two parallel solid plates. We set the coordinate with the ↔ Taking the dot product of G on both sides of the eq.(3), z axis normal to the plates and the x axis parallel to and taking the dot product of U~ on both sides of eq.(4), the moving direction of the liquid. The two plates both finally subtracting the two equations after the dot prod- are infinitely long to get rid off the ending effects of the ucts, we obtain structure. The two solid-liquid interfaces are located at z = +h and z = −h for the upper and lower plates [∇×∇×U~]·↔G−U~·[∇×∇×↔G]= 1F~·↔G−U~·↔Iδ(~r−~r′). (5) respectively. The friction force applied on the fluid is η considered to be opposite to the moving direction of the Integrating both sides of the eq.(5) over the computa- xaxis,whilethecomponentofthefrictionforcealongthe tional domain, we get the solution of U~ for the liquid ydirectionisaveragedtobezero. Forsimplicity,wethink with the form of the interfacesareflatalongthe y direction,andthenthe ↔ roughness of the interfaces is the function of x only. We U~(r′)=− [(nˆ×U~(r))·(∇×G(r,r′)) use the function ζ (x) to describe the roughness of the ZZ u upper interface anduse ζd(x) forthe lowerone. The two +(nˆ×∇×U~(r))·↔G(r,r′)]ds roughness functions both are averaged to be zero over 1 ↔ the interfaces at their locations of z = +h and z = −h + F~(r)·G(r,r′)dV. (6) respectively. ZZZ η 3 On the right hand side of the eq.(6), the first integral of non-ideal interfaces into the friction coefficient κ of the ds is the integral over the total surface of the com- ideally flat interfaces. The total friction coefficient κ RpRutational domain, which is obtained by the applying of the flat interfaces after the average then is the sum of the Green theorem on the volume integral of the left of κ and κ as mentioned in the introduction. The w r hand side of the eq.(5). The Green theorem has been total slip length b of the flat interfaces after the average shown in the Appendix A. The nˆ is the unit vector can be written as b = b b /(b +b ) with b = η/κ w r w r w w outward normal to the interfaces. Note that the direc- and b = η/κ by using the equalities of κ = κ +κ r r w r tionofthenˆvariesontheinterfacesduetotheroughness. and b = η/κ. It can be seen that for flat interfaces with b → ∞, the contribution of chemical interaction r Forsimplicity,weassumetheinletofthechannelalong dominates the friction force with b = b . On the other w thexdirectionisperiodictotheoutlet. Thustheintegral side, when b →∞, the friction force then is dominated w over the inlet and outlet cross-section surfaces would be by the roughness of the interfaces with b = b . So, after r cancelled. Similar cancellation is hold for the integral getting the total slip length b of the confined liquid, over the cross-section surfaces normal to the y direction the rough solid-liquid interfaces can be replaced by flat due to the uniformity along the direction. In this way, ones in an effective way, which simplifies the problems the term of surface integral in the eq.(6) remains as the of the fluid dynamics in real systems with the roughness integral over the solid-liquid interfaces only. Without having been involved in the b. In the following study, losing the generality, we will focus on only the upper we will set the b as a priori since the b depends on w w interface to figure out the contribution of the roughness the natural interaction between the solid surfaces and to the slip length, since the similar treatment can be fluid molecules, and focus on the informations of the applied to the lower interface. For convenience, we use roughness. K to represent the integrand of the surface integral by We start from the upper rough interface. For conve- ↔ ↔ K =−(nˆ×∇×U~)·G−(nˆ×U~)·(∇×G). (7) nience,wedenoteτˆtheunitvectorontheinterfacealong the tangential direction in the x,z plane and perpendic- ular to the nˆ. The τˆ can be expressed by τˆ = nˆ×U~×nˆ. |U~×nˆ| B. Boundary Condition Thelocalforcebalanceontheroughinterfaces,hencethe NBC, can be found generally to be The essential of the NBC is that the friction force ap- plied on the fluid per unit length should be balanced by U~ ·nˆ =0, (8) the shear stressof the fluid at the interface. The friction ηnˆ·(∇U~)·τˆ=−κ U~ ·τˆ, (9) force and the shear stress have the same magnitude but w the opposite directions. We denote the velocity compo- nˆ·(∇U~)·(τˆ×nˆ)=0. (10) nent of the fluid parallel to the interfaces by v . Then || In the general NBC, eq.(8) means the component of in the case that the solid-liquid interfaces are ideally flat, the shear stress is along the x direction and has U~ normal to the interfaces vanishes, showing that the been well defined as the product of the liquid viscosity liquid can not penetrate the solid surface. In the eq.(9), η and the normal gradient of the v to the interface, the term on the left hand side is the shear stress which || while the friction force is along the −x direction and is is along the τˆ direction and is on the plan normal to the product of the friction coefficient κ and the v at nˆ. The term on the right hand side of the eq.(9) is || the interface. Thus, for the ideal interfaces, the force the friction force. The final condition eq.(10) means no balance can be expressed by η∂v|| = −κ v valid at stress along the y direction according to our model. ∂z w || the interfaces, with the minus meaning the friction force It is emphasized here that the slip length used in the resists the fluid flow. This is exactly the NBC. Here, eq.(9) is the b = η/κ and no κ is involved. This is only the κ appears in the NBC, since no any effect of w w r w because the force balance considered now is hold locally. roughness of the interfaces contributes to the friction Theinvolvementofκ inthesliplengthbwillbeobtained coefficient in such ideal case. The slip length then is r aftertheaveragingoftheeq.(6)overtheroughinterfaces defined as b =η/κ . w w by using the general NBC. For the application of the generalNBC,weneedtomodifytheeq.(8,9,10)tohave In the non-ideal case that the solid-liquid interfaces the following forms are not ideally flat but with roughness, the NBC can be expressed locally since the force balance is still hold lo- U~ ·nˆ =0, (11) callyonthe interfaces. We willaveragethe fluidvelocity 1 eq.(6) over the rough interfaces to get an effective fluid U~ ·[(nˆ·∇)U~]=− U~ ·U~, (12) velocity. From an alternative view, the effective fluid bw velocitycanbesolvedfromtheNSequationbyimposing [(nˆ·∇)U~]·[nˆ×U~]=0. (13) an effective NBC(ENBC) on ideally flat interfaces. The ENBC has absorbed the effect of roughness of the TheabovemodificationcanbefoundintheAppendixB. 4 C. Resolution of K D. Perturbation Expansion WeneedtoresolutetheK definedineq.(7)intoseveral Wehaveusedthefunctionζ (x)todescribetherough- u terms to figure out the effect of the roughness. Before nessoftheupperinterface. Forconvenience,weomitthe the resolution, we decompose the vector ∇×U~ into two subscripttodenotetheroughnessbyζ(x). Thelocalunit vectors with one vector along the y axis and the other vectors of the interface can be expressed in the term of vector in the x,z plan, having ζ(x), reading [(∇×U~)·(nˆ×U~)](nˆ×U~) 1 ∂ζ ∇×U~ = nˆ = (− , 0, 1), (nˆ×U~)·(nˆ×U~) N ∂x 1 ∂ζ (nˆ×U~)×(∇×U~)×(nˆ×U~) τˆ= (1, 0, ), (21) + . (14) N ∂x (nˆ×U~)·(nˆ×U~) with By using the eq.(11, 12, 13,14), the vector n×∇×U in ∂ζ the K can be split into three terms, reading N = ( )2+1. (22) r ∂x nˆ×∇×U~ =T~1+T~2+T~3, (15) Now we need to average the fluid velocity U~ over the with rough interface to get an effective velocity W~ =< U~ >. −U~ ·U~(nˆ×nˆ×U~) Alternatively, the fluid vector W~ can be solved from the T~1 = , NS equation by imposing the ENBC on a flat interface. b (nˆ×U~)·(nˆ×U~) w The ENBC for such flat interface should have the form −nˆ·((U~ ·∇)U~)(nˆ×nˆ×U~) of T~2 = , (nˆ×U~)·(nˆ×U~) ∂W W x x W =W =0, =− . (23) (nˆ×U~)[nˆ·((nˆ×U~)·∇)U~] z y ∂z b T~3 = . (nˆ×U~)·(nˆ×U~) Here, the total slip length b has been used to comprise the two effects of the interface, the chemical interaction The details for the above derivation can be found in the and the roughness. The dyadic Green function for the AppendixC.WewritetheU~ intheformofU~ =Mτˆwith ↔ effective system with the flat interface is denoted by g. M the magnitude of U~. Then, we can further simplify ↔ The boundary conditions for the g are the same as the the terms of T~1,T~2,T~3 to be eq.(23)onlybyreplacingtheW withtheg whereiand i ij Mτˆ j represent the x, y and z coordinates. Note that g is T~1 = , (16) ↔ ij b the component of g, meaning that the influence of the w T~2 =Mτˆ(nˆ·(τˆ·∇)τˆ), (17) source of i component is↔on the object of j component. We expand the U~ and G of the rough interface in the T~3 =M(nˆ×τˆ)(nˆ·((nˆ×τˆ)·∇)τˆ). (18) termsofW~ and↔g oftheflatinterfacerespectively. Since In the above simplification, the equality of nˆ·τˆ= 0 has the friction force consideredis parallelto the x direction been applied. The similar treatment can be applied to and the y component of the friction force is negligible, the vector nˆ×U~ in the K, leading to weonlyconsiderthexcomponentofU~ fortheexpansion andneglecttheyandzcomponents. Thelattertwocom- nˆ×U~ =T~4 =Mnˆ×τˆ. (19) ponentsareaveragedto be zeroovertheroughinterface. ↔ Thus, the integrandK defined in the eq.(7) can be writ- In this way, for the G, we only consider the components ten as ofGxi. Wehavenotedthattheroughnessoftheinterface should not be too large, or the roughness will bring the ↔ ↔ K =−(T~1+T~2+T~3)·G−T~4·(∇×G). (20) convection and turbulence effects to the system. Under the condition of |ζ/b| << 1, the perturbation expansion It can be understood that for an ideally flat interface of the U~ ·xˆ and the G read where the τˆ is nota variablebut a constantvectoralong xi the x direction on the interface, the two terms of T~2 and ∂Wx ζ Mτˆ·xˆ=W +ζ =W (1− ), T~3 bothdisappearbecausetheoperator∇appliedonthe x ∂z x b constantvectorτˆcauseszerointheeq.(17,18). Andthen ∂g ζ xi only the terms of T~1·↔G and T~4·(∇×↔G) remain in the Gxi =gxi+ζ ∂z =gxi(1− b). (24) K for the ideally flat interface. Therefor, the terms of In above derivation, the ENBC of eq.(23) have been ↔ (T~2+T~3)·G in the K are originatedfrom the roughness applied. of the interface. 5 In the following, we expand the K in the terms of W the correlation function is Gaussian, which reads and g. For the first term T~1·↔G of the K, we get < ζ(x′)ζ(x′ +x) >= σ2e−x2/l2. Here, σ is the variance of the roughness. The x is the distance between two T~1·Gxixˆˆi= Wx(1− ζ)2gxiˆi. (25) positions on the interface for the correlation. And the l bw b is the correlationlength, meaning that the two positions with their distance larger than l lose their correlation. By substituting the eq.(21,24) into the second term of Theinterfacewithalongerl hasasmoothertopography. the K, we get In the expression eq.(29) of b, the information of l has T~2·Gxixˆˆi=Wx(1− ζ)2(nˆ·(τˆ·∇)τˆ)gxiˆi been reflected by the factor of ∂∂x2ζ2 in the A. We have b used Mont Carlo method to simulate various rough in- ζ ∂2ζ terfacesbyvaryingσ andl. The latticeparameterofthe =W (1− )2 N−3g ˆi (26) x b ∂x2 xi solid surface is used for the dimensionless. The physical quantities in this work are evaluated statistically over We observe that the vectors T3 and T4 both are along the rough interfaces. Fig. 1 shows that the < N > is they direction,whichhavenegligiblecontributiontothe functional of the correlation length l of the interface effect of the roughness as mentioned in this work. Thus, by varying the roughness variance σ. In the figure, we we will figure out the roughness information only from the two terms of eq.(25) and eq.(26). E. Slip Length 100 Thesurfaceintegraloftheeq.(25)andeq.(26)obtained =4 from the eq.(6) can be effectively equivalent to =32 W >10 =64 < [T~1·Gxixˆˆi+T~2·Gxixˆˆi]ds>= xgxiˆidxdy. N Z Z Z Z b < (27) The right hand side of the above equation is the sur- face integral over the flat interface after the average and 1 canbe understoodfromthe eq.(16)with the eq.(17) and eq.(18) vanishing at the flat interface. Here, b is used to 0 60 120 180 240 300 comprisethe botheffects ofthe chemicalinteractionand l theroughnessfortheflatinterface. Therefore,wegetthe FIG. 1: Statistical average of N over interfaces, noted equation for the total slip length b, reading by <N >, as a function of the correlation length l with N ζ ∂2ζ ζ 1 various roughness variances σ. The N has been defined < (1− )2+ N−2(1− )2 >= (28) b b ∂x2 b b in the text by eq.(22). And the <N > is plotted in the w logarithmic scale. by using the ds = Ndxdy. The factor of W g ˆidxdy onRbRoth sides RofRthe equation are canceled, x xi find that the < N > deviates from the number of one meaning that the b is a character of the confined liquid, in the range of l close to zero. Larger the σ, stronger andisindependent onthe velocityoffluidandtime. For simplicity, we denote A = 1 ∂2ζ + N . After averaging the deviation of the < N >. With the increasing of N2∂x2 bw l, the < N > decreases and approaches to the unit, both sides of the eq.(28) statistically over the interface, indicating that the interface is becoming smooth with we solve the slip length b as ∂ζ going to zero and τˆ deviating from the unit vector ∂x 2<ζA>+1+ (2<ζA>+1)2−4<A><ζ2A> xˆ slightly. A rougher interface with a larger σ need b= . a larger l to smooth the interface and to recover the p 2<A> < N > to the unit. In the fig.1, the deviation of the (29) < N > actually means that the l and the σ both play This is the main result of this work. their roles in determining the slip length. And the l gradually weakens its effect in the determination with the increasing of l. Finally, only the σ remains to de- III. RESULTS AND DISCUSSIONS terminethesliplengthifthelgoestothelimitofinfinity. The topography of the rough interface is described Fig. 2 shows the slip lengths as functional of l with by the ζ(x) together with the correlation func- various σ. It should be noted that in the eq.(24) the tion of the roughness. For simplicity, we set that |ζ/b| << 1 is used for the perturbation expansion 6 instead of |ζ| <<1. Thus, we use σ/b to represent the w variance of the roughness in the plot instead of σ itself, 1.2 (a) considering that the b is equivalent to b in the extreme w 1.0 case of σ =0. We apply σ/b ≤0.4 in the plot to verify w the perturbation expansion. Fig.2(a) and (b) are for 0.8 the total slip length b basing on the eq.(29), while the w b b / 0.6 bw=10 =0.1 fig.2(c) and (d) are for the slip length br contributed =0.2 from the roughness only. The b is obtained from 0.4 = /bw =0.3 1 = 1 − 1 . The slip lengths are rrenormalized by the =0.4 br b bw 0.2 b in the figures for clarity. The b has been indicated w w in each figure. In the fig.2(a) and (b), it can be found 0.0 0 50 100 150 200 250 300 that the b increases dramatically with the increasing of l thel andthengoestoasaturationvalue. Thesaturation value of the b is independent on the l, but decreases 1.2 with the increasing of σ as expected. The increasing of (b) b with the l emphasizes that in order to increase the 1.0 b it is essential not only by choosing proper material 0.8 for the channel fabrication, but also by smoothing the w solid surfaces in the fabrication. The b in the fig.2(c) b / b 0.6 and (d) are calculated from the data ofrfig.2(a) and (b) =0.1 0.4 respectively, showing the similar behaviors of b in the bw=80 =0.2 =0.3 fig.2(a) and (b). In the fig.2(d), we have appended the 0.2 = /bw =0.4 scale of l to 1000 because the b /b does not reach the r w 0.0 saturation at the l = 300. The fig. 2(c) and (d) show 0 50 100 150 200 250 300 l that the slip length contributed from the roughness is a function of l and σ in the range of l close to zero, and finally is a function of σ only when the l goes to the 90 (c) limit of infinity. We need to note that in the fig.2 the scaling behaviors of the slip lengths are dependent on the ratio of σ to b rather than the σ itself. And for the 60 w w =0.1 interfaces with large roughness where the perturbation b / br bw==1/b0w ==00..23 expansion is invalid, say σ/bw > 0.4, rigorous numerical =0.4 calculations such as the finite element method should 30 be applied to solve the NS equation, which is out of the scope of this work. 0 0 50 100 150 200 250 300 Before we figure out the relation between the the l b(b ) and the b , we firstly make a simple analysis. We r w approximate the following equalities of < A >= 1/b , w 90 < ζA >=< ζ ∂2ζ > and < ζ2A >=< ζ2 > /b and (d) N2∂x2 w find that b can be expressed by b ∼ b < f(ζ, ∂2ζ) >. w ∂x2 That means b is linearly proportional to b and the 60 w w roughness of the interface plays its role as the propor- b / br bw=80 ==00..12 tionality factor. Now we check the proportion in fig.3 =0.3 with various σ. Here, the correlation length is fixed 30 = /bw =0.4 at l = 200 in the plot and σ is used instead of the σ/b . Fig.3(a) shows that the ratio of the b to the b w w goes up to a constant with the increasing of b . The 0 w 0 200 400 600 800 1000 constant ratio confirms the linear proportion of the b to l the b and the roughness plays as the proportionality w FIG. 2: Slip lengths as functions of the correlation factor. Larger the σ, smaller the ratio, meaning that length l with various roughness variances σ. (a) and (b) the roughness weakens the effect of bw. Fig.3(b) is for are for the total slip length b of the system while the (c) the br calculated from the data of fig.3(a), also showing and (d) are for the slip length b contributed from the the linear behavior as the b. The relationship between r roughness of interfaces. In the plot, the slip lengths and the br and the bw indicates that the roughness of the σ are all renormalizedby the slip length b contributed interface can not play the role itself in determining the w from the chemical interaction of the interfaces. b, but need to be with the bw together. That means if b is infinite, b must be infinite too and then is the b. w r 7 It has been revealed that some biological nanopores havepeculiar features in the transportof water. Such as IV. CONCLUSION intheaquaporinchannelsasthekeycomponentofmany biological processes, the permeability of water is larger In this work, we have studied the slip length of the than what is expected in classicalnanofluidic framework by three orders of magnitude [43]. It has been suggested confinedliquidwithsmallroughnessofsolid-liquidinter- faces. The quantitative expression of the slip length has that the hydrophobic surface of the aquaporin channels is the potential mechanism to the large permeability of been obtained, in which the roughness of the interfaces and the chemical interaction between the liquid and the water in the channels. To supplement the suggestion, we think that the large correlation length of the chan- solid are involved. In the study, perturbation expansion of the interface roughness is required. The results we nels may contribute to the water transport as have been studied in this work, even though how the positions are have obtained are threefold. First, the slip length b of the confined liquid is linearly proportional to the slip correlated in the aquaporin channels is not clear. What length b contributed from the chemical interaction. is more, the biological nanopores are active matters, in w Second, the roughness of the interface plays its role which the correlation between the two positions on the as a proportionality factor for the b and the b . And nanopore interfaces may behave in a complicated way w only with the chemical interaction together does the ratherthanwhatdoesin the artificialnanofluidic frame- roughnesshaveits effect onthe b. Finally,largevariance works. Theconsiderationofthecorrelationlengthinthe of the roughness decreases the b and the increasing of biological nanopores is the extension of this work and the correlation length of the roughness increases the under the research. slip length in a dramatical way to a saturation value. Those results propose that we need to improve our techniques to smooth the solid surfaces by enlarging the correlation length of the roughness as well as Appendix A decreasing the variance of the roughness for the channel fabricationbesidestofindpropermaterialswithlargeb . w The second kind Green theorem is ↔ ↔ [(∇×∇×U~)·G−U~ ·(∇×∇×G)]dV ZZZ 1.0 (a) = [(nˆ×U~)·(∇×↔G)+(nˆ×∇×U~)·↔G]ds. (A1) ZZ 0.9 w0.8 b b / l=200 0.7 Appendix B =10 =20 0.6 =40 We need to modify the NBC eq.(8, 9, 10)for the appli- 0.5 cationasthefirststepinsplittingthe K. The eq.(8)can 0 200 400 600 800 1000 bw be written in the form of U n =0, (B1) 80 a a (b) 70 60 l=200 with the Einstein’s notation. Here, the subscript a represents the coordinates x,y,z. w 50 =10 b b / r 40 ==2400 Then, we substitute the τ = nˆ×U~×nˆ into the eq.(9), 30 |U~×nˆ| 20 and get rid of the common denominators on both sides 10 of the equation. In this way, we get the modified eq.(9) reading 0 0 200 400 600 800 1000 bw 1 nˆ·(∇U~)·[nˆ×U~ ×nˆ]=− U~ ·[nˆ×U~ ×nˆ]. (B2) FIG. 3: (a)Proportion of b to bw. (b)Proportion of br to bw b . In the plot, the b and b have been renormalized by w r the b for clarity. w The left hand side of the eq.(B2) can be simplified by 8 using the Einstein’s notation Thus, for the eq.(C1), we have nˆ·(∇U~)·[nˆ×U~ ×nˆ] ∂U =(na~ia)·(~ib∂xc~ic)·(ǫefg~ienfǫghkUhnk) nˆ×[(∇×U~)·(nˆ×U~)](nˆ×U~) b ∂U 1 =ǫ ǫ n en U n =[− U~ ·U~ −nˆ·[(U~ ·∇)U~]](nˆ×nˆ×U~). (C2) efg ghk b∂x f h k b b w ∂U And e =(δ δ −δ δ )n n U n eh fk ek fh b f h k ∂x b ∂U =n U e =U~ ·[(n·∇)U~]. (B3) b e ∂x b (nˆ×U~)×(∇×U~)×(nˆ×U~) Here, the eq.(B1) has been applied. Similarly, the right ∂U ∂U c b hand side of the eq.(B2) can be simplified as − 1 U~ ·U~. =ǫbefianeUf(ncUa∂x +naUc∂x bw b c So the eq.(B2) can be transformed as ∂U ∂U c b −n U −n U ) a c c a ∂x ∂x 1 b c U~ ·[(n·∇)U~]=− U~ ·U~. (B4) =U~nˆ·[(nˆ×U~)·∇]U~ +nˆ(nˆ×U~)·[(U~ ·∇)U~] b w −nˆU~ ·[(nˆ×U~)·∇]U~ −U~[(nˆ·∇)U~ ·(nˆ×U~)]. This equation will be used for the splitting of the K. (C3) The eq.(10) then can be transformed to be ∂U nˆ·(∇U~)·(U~ ×nˆ)=ǫbefnaneUf b =0 The last term of the above equation is dropped off by ∂x a using the eq.(B8). For the eq.(C3), we have ⇒[(n·∇)U~]·[n×U~]=0. (B5) For clarity, we collect the modified NBC here U~ ·nˆ =0, (B6) nˆ×(nˆ×U~)×(∇×U~)×(nˆ×U~) 1 =(nˆ×U~)[nˆ·[(nˆ×U~)·∇]U~], (C4) U~ ·[(nˆ·∇)U~]=− U~ ·U~, (B7) b w [(nˆ·∇)U~]·[nˆ×U~]=0. (B8) By using the equality of nˆ × nˆ = 0. In this way, the Appendix C nˆ×∇×U~ can be split into three terms by Before we resolve the K, we need to decompose the ∇×U by two vectors. One vector is along the nˆ ×τˆ. The other vector is in the y,z plan. The decomposition −U~ ·U~(nˆ×nˆ×U~) T1 = , has been shown in the eq.(14). Now we simplify the nu- b (nˆ×U~)·(nˆ×U~) w merators of the equation. The factors of the terms are −nˆ·((U~ ·∇)U~)(nˆ×nˆ×U~) simplified as the following. T2 = , (nˆ×U~)·(nˆ×U~) ∂U ∂U (∇×U~)·(nˆ×U~)=n U k −n U k (nˆ×U~)[nˆ·((nˆ×U~)·∇)U~] j k∂xj k j ∂xj T3 = (nˆ×U~)·(nˆ×U~) . =U~ ·[(nˆ·∇)U~]−nˆ·[(U~ ·∇)U~] 1 =− U~ ·U~ −nˆ·[(U~ ·∇)U~]. (C1) b w [1] J.C.T.EijkelandA.vandenBerg,Nanofluidics: whatis [2] Lyderic Bocquet and Elisabeth Charlaix, Nanofluidics, itandwhatcanweexpectfromit? MicrofluidNanofluid., from bulk to interfaces, Chem. Soc. Rev.39,1073(2010) 1,249(2005) [3] Z. Siwy and A. Fulinski, Fabrication of a Synthetic Nanopore Ion Pump, Phys. Rev.Lett., 89,198103(2002) 9 [4] R. Karnik, R. Fan, M. Yue, D. Li, P. Yang and A. Ma- [24] Y. Zhu and S. Granick, Rate-Dependent Slip of New- jumdar, Electrostatic Control of Ions and Molecules in tonian Liquid at Smooth Surfaces, Phys. Rev. Lett., Nanofluidic Transistors, NanoLett., 5,943(2005) 87,096105(2001) [5] R. B. M. Schasfoort, S. Schlautmann, J. Hendrikse and [25] Tingyi Liu, Chang-Jin Kim, Turning a surface super- A.vandenBerg,Field-EffectFlowControl forMicrofab- repellent even to completely wetting liquids, Science , ricated Fluidic Networks,Science, 286,942(1999) 346,1096(2014) [6] R.Karnik,C.Duan,K.Castelino,H.DaigujiandA.Ma- [26] J.Ou,B.PerotandJ.P.Rothstein,Laminardragreduc- jumdarn,Rectification of Ionic Current in a Nanofluidic tion in microchannels using ultrahydrophobic surfaces, Diode, NanoLett., 7,547(2007) Phys. Fluids, 16,4635(2004) [7] D. Y. Chan and R. G. Horn, The drainage of thin [27] J. Ou and J. P. Rothstein, Direct velocity measurements liquid films between solid surfaces, J. Chem. Phys., of the flow past drag-reducing ultrahydrophobic surfaces, 83,5311(1985) Phys. Fluids, 17,103603(2005) [8] U.RavivandJ.Klein,FluidityofBoundHydration Lay- [28] C. H. Choi, U. Ulmanella, J. Kim, C. M. Ho and ers, Science, 297,1540(2002) C. J. Kim, Effective slip and friction reduction in [9] L. Bocquet and J. -L. Barrar, Flow boundary conditions nanograted superhydrophobic microchannels, Phys. Flu- from nano- to micro-scales, Soft matter, 3,685(2007) ids, 18,087105(2006) [10] A. Maali, T. Cohen-Bouhacina, G. Couturier and J. - [29] C. H. Choi and C. J. Kim, Large Slip of Aqueous Liquid P. Aime, Oscillatory Dissipation of a Simple Confined Flow over a Nanoengineered Superhydrophobic Surface, Liquid,Phys. Rev.Lett., 96,086105(2006) Phys. Rev.Lett., 96,066001(2006) [11] J. -M. Georges, S. Millot, J. -L. Loubet and A. Tonck, [30] M.CalliesandD.Quere,On water repellency, SoftMat- Drainage of thin liquid films between relatively smooth ter, 1, 55 (2005) surfaces, J. Chem. Phys., 98,7345(1993) [31] R.Truesdell,A.Mammoli, P.Vorobieff,F.vanSwoland [12] T. Becker and F. Mugele, Nanofluidics: Viscous Dis- C. J. Brinker, Drag reduction on a patterned superhy- sipation in Layered Liquid Films, Phys. Rev. Lett., drophobic surface, Phys.Rev.Lett., 97, 044504(2006) 91,166104(2003) [32] L.Bocquet,P.TabelingandS.Manneville,Commenton [13] Y. Leng and P. T. Cummings, Fluidity of Hydration ”Large slip of aqueous liquid flow over a nanoengineered Layers Nanoconfined between Mica Surfaces, Phys. Rev. superhydrophobic surface”, Phys. Rev.Lett.,97, 109601( Lett., 94,026101(2005) 2006) [14] T. -D. Li, J. Gao, R. Szoszkiewicz, U. Landman and [33] P. Joseph, C. Cottin-Bizonne, C. Y. J.-M. Benoit, C. E. Riedo, Structured and viscous water in subnanometer Journet, P. Tabeling and L. Bocquet,Slippage of wa- gaps, Phys.Rev.B., 75,115415(2007) ter past superhydrophobic carbon nanotube forests in mi- [15] G. K. Batchelor, an Introduction to Fluid Dynamics crochannels, Phys.Rev.Lett., 97, 156104(2006) (Cambridge Univeristy Press, Cambridge, 1967) [34] E. Lauga and H. A. Stone, Effective slip in pressure- [16] L. Bocquet and J. -L. Barrat, Hydrodynamic boundary driven Stokes flow ,J. Fluid Mech., 489, 55(2003) conditions, correlation functions, and Kubo relations for [35] S. Richardson, On the no-slip boundary condition, J. confined fluids, Phys.Rev. E., 49,3079(1994) Fluid Mech., 59, 707(1973) [17] R.Pit,H.HervetandL.Leger,DirectExperimental Evi- [36] C.Cottin-Bizonne,J.-L.Barrat,L.BocquetandE.Char- denceofSlipinHexadecane: SolidInterfaces,Phys.Rev. laix, Low-friction flows of liquid at nanopatterned inter- Lett., 85,980(2000) faces, Nat. Mater., 2, 237(2003) [18] C.Cottin-Bizonne,C.Barentin,E.Charlaix,L.Bocquet [37] H. W. Hu, G. A. Carson, and S. Granick, Relaxation and J. -L. Barrat, Dynamics of simple liquids at hetero- time of confined liquids under shear, Phys. Rev. Lett., geneous surfaces: Molecular-dynamics simulations and 66,2758(1991) hydrodynamic description,Eur.Phys.J.E,15,427(2004) [38] J.N.Israelachvili,P.M.McGuiggan,andA.M.Homola, [19] N. V. Priezjev and S. M. Troian, Molecular Origin and Dynamic Properties of Molecularly Thin Liquid Films, Dynamic Behavior of Slip in Sheared Polymer Films, Science, 240,189(1988) Phys.Rev.Lett., 92,018302(2004) [39] ShuyuChen,HanWang,TiezhengQian,andPingSheng, [20] A.Jabbarzadeh,P.HarrowellandR.I.Tanner,Lowfric- Determining hydrodynamic boundary conditions from tionlubrication between amorphous walls: unravelingthe equilibrium fluctuations, Phys. Rev.E, 92,043007(2015) contributions of surface roughness and in-planedisorder, [40] P. A. Thompson and S. M. Troian, A general bound- J. Phys. Chem., 125,034703(2006) ary condition for liquid flow at solid surfaces, Nature, [21] P. A. Thompson and M. O. Robbins, Shear flow near 389,360(1997) solids: Epitaxial order and flow boundary conditions, [41] P.A.Thompson, G.S.Grest and M.O.Robbins,Phase Phys.Rev.A, 41,6830(1990) transitions and universal dynamics in confined films, [22] C.Cottin-Bizonne,B.Cross,A.SteinbergerandE.Char- Phys. Rev.Lett., 68,3448(1992) laix,BoundarySliponSmoothHydrophobicSurfaces: In- [42] E.D.Smith,M.O.RobbinsandM.Cieplak,Frictionon trinsic Effects and Possible Artifacts, Phys. Rev. Lett., adsorbed monolayers, Phys. Rev.B, 54,8252(1996) 94,056102(2005) [43] T.Walz,B.L.Smith,M.L.Zeidel,A.EngelandP.Agre, [23] M. Cieplak, J. Koplik and J. R. Banavar, Boundary Biologicallyactive two-dimensional crystals of aquaporin Conditions at a Fluid-Solid Interface, Phys. Rev. Lett., CHIP, J. Biol. Chem., 269,1583(1994) 86,803(2001)