南怡伶，孔憲，李繼鵬，盧滇楠

?

納米狹縫中水流動非平衡分子動力學模擬

南怡伶，孔憲，李繼鵬，盧滇楠

（清華大學化學工程系，化學工程聯合國家重點實驗室，北京 100084）

采用非平衡分子動力學模擬(non-equilibrium molecular dynamics simulation)方法研究了不同間距納米狹縫之間水的流動行為。研究了納米狹縫間距、壁面性質和外部壓力對水流動速度徑向分布、有效黏度、壁面速度和滑移長度的影響，討論了Navier-Stoke (N-S)方程的適用性。研究結果表明，N-S方程僅適用于3 nm以上的孔道；狹縫尺寸的增加和施加壓力的增加均會使得管內流速增加，而造成表觀黏度降低以及滑移長度增加。壁面親水性的增加僅使得滑移長度降低，表觀黏度并沒有發(fā)生較大變化。

納米狹縫；分子動力學模擬；流動；有效黏度；滑移長度；非平衡分子動力學；表面

引 言

納米狹縫中水的流動特性是生物化學分析[1]、水鹽分離[2]和能量轉化[3]等技術的共性科學基礎，也是納米流體流動領域研究的前沿和基礎。納米狹縫是一種較常見的納米結構之一，是一種狹縫間距小于100 nm，且狹縫長度遠大于其間距的幾何結構[4]。隨著石墨烯等材料在膜科學中越來越多的應用，狹縫孔道的表面性質，狹縫尺寸等結構因素和外場壓力等操作因素對納米狹縫中水流動行為的影響對于膜材料設計和膜分離優(yōu)化至關重要。

區(qū)別于宏觀的水，在納米狹縫中水的平衡結構和傳遞行為均具有一些獨特的性質。在較小的狹縫孔中，水呈現層狀分布，并且隨著狹縫的增加，層狀分布逐漸變弱[5-7]。Algara-Siller等[8]和Nair等[9]發(fā)現了在石墨烯狹縫孔中，水可以形成“方形冰”結構。實驗[9-11]和模擬[6, 12-13]均發(fā)現，石墨烯或者氧化石墨烯狹縫孔中的水要比主體相中的水傳遞速度快，并且受管壁和管的尺寸影響較大。這是因為在宏觀流動下，水的流動主要受水本身性質的影響；而在納米狹縫中，隨著狹縫比表面積的增加，水的流動顯著受到狹縫壁面表面力的影響。因此水的流動同時受到壁面表面力和水自身性質的雙重影響，從而導致水的相圖、結構和動力學等方面產生許多獨特的性質[14]。

一般認為，在納米狹縫孔中，傳統(tǒng)的流體動力學理論如Navier-Stoke (N-S)方程不再適用。因此以前的納米狹縫孔的模擬工作，很少利用N-S方程進行關聯[6, 12-13]。然而，根據目前的研究結果，在大于流體分子尺寸10倍的孔道內，N-S方程具有一定的適用性。Travis 等[15]模擬了在兩種不同間距納米狹縫中，Lenard-Jones（LJ）流體的泊肅葉流動。研究表明，當納米孔道直徑小于10倍LJ流體分子直徑時，水的速度分布不能用N-S方程來描述。Bo?an等[16]模擬了黏土狹縫中（2～9 nm）水的流動，也得到類似的結論。Sendner等[17]通過改變液體/固體的相互作用力來改變接觸角，得到了納米孔道流動模型中滑移長度與水潤濕性的關系式。其研究發(fā)現黏性邊界條件只能用于弱液體/壁面相互作用。Huang等[18]研究了Couette流動中表面電荷密度對滑移長度的影響，結果表明隨著表面電荷密度的增加，滑移長度降低。

本文利用非平衡分子動力學模擬（NEMD）技術研究了納米狹縫中水的流動行為。系統(tǒng)研究了狹縫大小、驅動壓力和表面親水性3個因素對納米狹縫中的水流動行為的影響，進一步驗證了N-S方程的適用性。對于適合用N-S方程進行描述的體系，分析了上述因素對于流體的表觀黏度和滑移長度這兩個表征流體動力學行為參數的影響。上述研究對于基于流動納米器件的分子設計可提供一定的理論指導。

1 模型和模擬方法

1.1 模型

石墨烯狹縫孔不僅實驗中可以制備，其理想的結構也使其成為理論研究中的普遍模型[19-20]。因此采用石墨烯作為納米狹縫的壁面。水在石墨烯納米狹縫中的流動模型如圖1所示。納米狹縫壁面由結構相同的兩片armchair構型的石墨烯平板組成，石墨烯平板的長度和寬度分別為 6.140 nm和3.403 nm。模擬中采用了周期性邊界條件，在該尺寸下，模擬盒子中的石墨烯與其周期性鏡像的石墨烯正好可以組成一張無限大的完美石墨烯片層。SPC/E[21]水分子填充于兩個石墨烯平板之間的狹縫中。為了避免方向上的鏡像對于體系的影響，在石墨烯平板兩側各留出了1 nm的真空層。為了研究狹縫間距（）的影響，構造了石墨烯間距分別為1、2、3、4、5、6和10 nm的7個水流動模型。不同間距模型所添加的水分子個數如表1所示，充填后的水密度均為0.997 g·cm-3。為了研究壁面親疏水性對于水流動的影響，改變了C原子的LJ參數[22]，建立了C1、C2兩種新的原子。其中C1的LJ參數是SPC/E水模型中O的參數，C2的相互作用勢參數是C1的2倍。

表1 不同模擬條件下納米狹縫中水分子個數

1.2 模擬方法

分子動力學模擬所用的軟件為NAMD[23]。模擬過程中，碳原子位置固定。采用Particle Mesh Eward (PME)方法[24]計算電荷相互作用。LJ和靜電相互作用的截斷距離均為1.2 nm。模擬的步長為1.0 fs，每隔0.5 ps保存1幀軌跡。

在非平衡模擬之前，首先進行時長為6 ns的NVT系綜的平衡模擬。平衡模擬中，利用Langevin方法將系統(tǒng)的溫度維持在300 K。系統(tǒng)穩(wěn)定后，便作為初始構象進行NVE系綜下的非平衡動力學模擬，即不對水分子進行控溫。該方法可以使得水分子的運動符合牛頓經典力學，是模擬受限空間中流體力學的合適方法[25-26]。在非平衡動力學模擬中，每個水分子的氧原子在方向上均施加一個固定大小的力，從而在方向上形成水的定向流動。由于體系的均一性及周期性，這對應于水在無限長納米狹縫孔的流動。非平衡模擬時長為10 ns，非平衡模擬中所有模型均在2.5 ns之后達到穩(wěn)定。

表2 石墨烯與水分子間相互作用參數

1.3 分析方法

納米狹縫中水分子的密度分布()計算公式如下

式中，()為=[,+D)范圍內水分子的平均密度（根據氧原子的坐標進行判斷）；m,H2O為第個水分子的質量；A為模擬盒子底面積，其值為20.56 nm2；D為計算水分子密度沿軸分布時所取的步長，這里D=0.1 nm。

狹縫中不同處水分子在方向上的速度分布u()以氧原子在=[,+D)區(qū)間平均速度來計算，即

式中，u,i為第個水分子中O原子在方向上的速度；(z-)和(+D2-z)均為Heaviside階躍函數。本文中，D2=0.2 nm。

根據N-S方程對狹縫內水分子的速度分布進行擬合。本文模型中，僅方向存在穩(wěn)態(tài)流動，因此N-S方程可以化簡為

假定剪切黏度和流動方向壓力梯度與無關，則其解的形式為

式中，1和2是與無關的常數，在如圖1所示坐標定義條件下，1=0。

在光滑的疏水壁面，水的邊界速度往往不為0，可以用滑移長度（）來表征流體水在壁面處的運動[17]。值越大，則表明壁面對水流動的阻力越小。的定義為

= u(0)/′(0) (5)

式中，u(0)為壁面處流體水的速度，而′(0)=du(0)/d為流體水在壁面速度u(0)對的一階導數。

2 模擬結果與討論

2.1 驅動壓力對于狹縫內水分子密度分布的影響

首先比較了驅動壓力對狹縫內水分子密度分布的影響，結果如圖2所示。圖中狹縫寬度分別為3.0、5.0和10.0 nm。

圖2結果表明，在沒有外界驅動壓力條件下，不同納米狹縫內水分子均呈現典型的近壁有序、遠壁無序的特點，這與傳統(tǒng)平衡動力學模擬是一致的[6, 27-28]。當狹縫寬度大于3.0 nm之后，狹縫中心水密度已經與體相水沒有區(qū)別。在流動狀態(tài)下，狹縫內水分子的密度分布趨勢與非流動狀態(tài)相同。值得注意的是，流動對于貼近壁面的第1個水分子層產生顯著影響，隨著驅動壓力的增加，第1個水分子層峰高降低且趨向壁面。這說明在流動狀態(tài)下，壁面處水的結構有序性會降低。

2.2 流動特性

2.2.1 狹縫間距對流動的影響 研究了狹縫間距對流體水平均速度和速度分布的影響，結果如圖3和圖4所示。在模擬過程中，施加于水的驅動壓力為1035 Pa。

圖3結果表明，當流體處于穩(wěn)定流動狀態(tài)下，整體上流體的平均流速隨著狹縫寬度的增加而增加。但在2 nm以下的孔道內，流體的平均速度會有一些偏移。這是由于2 nm以下的孔道內，水的流動不再是連續(xù)介質流動[29]。

圖4給出了不同尺寸的孔道中，水在方向的速度分布。結果表明，當板間距大于等于3 nm時流形呈拋物線分布，符合宏觀流體中的N-S方程。這與Bo?an等[16]研究水在不同尺寸的黏土孔道中的結果一致。而低于3 nm時，流形沒有呈明顯的拋物線形，說明N-S方程不再適用。在低于2 nm的納米狹縫中，水分子之間致密的氫鍵網絡會強化水分子間相互作用，從而導致類似平推流的均勻速度分布，而在3 nm及其以上的納米狹縫中，在遠離壁面的水分子中心氫鍵網絡不再致密和規(guī)整，這種差異性導致N-S方程對水分子和其他LJ流體適用性的差異。

進一步利用式(4)對狹縫間距分別為3、4和5 nm的速度分布進行擬合，得到納米狹縫中水的表觀黏度和滑移長度，結果如表3所示。

表3 納米狹縫中水的表觀黏度

表3結果表明，流體的表觀黏度均比純水體系的黏度(～100×10-5Pa·s)高。隨著狹縫間距的增加，狹縫內水的表觀黏度降低。這可以解釋圖3的結果，即狹縫間距增加導致的表觀黏度的降低有助于流體流動速度的增加。從表3還可看出，滑移長度與狹縫間距呈正相關。這說明隨著狹縫間距的增加，狹縫壁面對于水流動的阻礙作用會逐漸減弱。

2.2.2 驅動壓力對流動的影響 進一步研究了驅動壓力對狹縫內水流動平均速度和速度分布的影響，結果如圖5和圖6所示。狹縫寬度為5.0 nm，根據前面的研究，狹縫內水的流動分布可以用N-S方程來擬合。

圖5結果表明，狹縫內水的平均流速隨著驅動壓力的增加而增加，且呈線性關系，這符合N-S方程的描述，并且與之前的模擬結果一致[12, 29]。

狹縫內水的速度分布如圖6所示，速度分布均呈拋物線分布，可以用式（4）進行擬合，得到狹縫內水的剪切黏度和滑移長度，結果如表4所示。

表4表明，狹縫內水的表觀黏度隨著驅動壓力的增加而降低。而滑移長度隨著驅動壓力的增加而增加。這說明隨著驅動壓力的增加，流體的運動速度增加，壁面對流體的阻礙作用逐漸削弱。

表4 驅動壓力對狹縫內水的表觀黏度和滑移長度的影響

2.2.3 壁面親水性對流動的影響 進而通過改變納米狹縫中C原子的LJ勢能參數，研究了壁面親疏水性對于水的流動行為的影響，結果如圖7所示。針對狹縫寬度5 nm的模型，驅動壓力為1035 Pa。需要指出的是，這里親水性壁面并不包括長程靜電相互作用。因為前期多壁碳納米管的研究中，對多壁碳管外壁進行電荷修飾亦會對碳管內水分子流動行為產生顯著影響[30]。

圖7結果表明，隨著親水性的增加，管內流速降低，但仍然呈拋物線分布。高親水性壁面阻礙了水的流動，從而造成流速的降低。利用式（4）進行擬合，得到了狹縫內水的剪切黏度和滑移長度，結果如表5所示。

表5 壁面親疏水性對水的表觀黏度和滑移長度的影響

表5結果表明，隨著壁面親水性的增強，滑移長度降低，表明壁面對于水流動的阻力增加。然而表觀黏度并不會隨著壁面親疏水性的變化有明顯的改變。這說明對于表觀黏度屬于流體性質，壁面對其影響不大。

3 結 論

本文采用非平衡分子動力學模擬的方法研究了納米狹縫中水分子的流動特性。在納米狹縫中，壁面附近水分子密度高，而距離壁面約1.0 nm處，水分子密度逐漸趨于宏觀水的密度。水的流動使得壁面處水的結構有序性降低，具體體現在施加壓力形成水的定向流動后，壁面第1個水分子的峰高降低，并且向壁面更加靠近。當狹縫寬度小于3.0 nm時，速度分布呈平推流的模型。只有當狹縫寬度大于等于3.0 nm時，水的速度分布符合二次型分布，水的流動行為才可用N-S方程進行描述。不過，N-S方程的參數（黏度及滑移長度）不再是固定值。狹縫寬度和施加的壓力會影響水的表觀黏度和滑移長度。表觀黏度與狹縫寬度或施加壓力呈負相關，滑移長度與狹縫寬度或施加壓力呈正相關。這表明隨著狹縫寬度的增加或者施加壓力的增加，壁面對于水流動的影響逐漸減弱。壁面親水性的增加會使得滑移長度降低，但是表觀黏度不會有明顯變化。在未來的研究中，如何從理論上描述滑移長度和表觀黏度的變化，對于構建納米孔道流動理論是至關重要的。

符 號 說 明

Axy——模擬盒子底面積，nm2 d——石墨烯平板間距，nm mi,H2O——第i個水分子的質量 Nwater——狹縫中水分子個數 ——流動方向的壓力梯度 Dp——驅動壓力，Pa q——帶電量，e ux(z0)——壁面處流體水的速度，m·s-1 ux,i——第i個水分子中O原子在x方向上的速度，m·s-1 Dz——沿z軸所取的步長，nm e——LJ勢能參數，kJ·mol-1 l——滑移長度，nm m——狹縫中水的剪切黏度，Pa·s Q(zi-z)——Heaviside階躍函數 r(z)——納米狹縫中水分子的密度分布

References

[1] GUAN W, FAN R, REED M A. Field-effect reconfigurable nanofluidic ionic diodes[J]. Nature Communications, 2011, 2(2): 506.

[2] JAIN T, GUERRERO R J S, AGUILAR C A,Integration of solid-state nanopores in microfluidic networkstransfer printing of suspended membranes[J]. Analytical Chemistry, 2013, 85(8): 3871-3878.

[3] DAVIDSON C, XUAN X. Electrokinetic energy conversion in slip nanochannels[J]. Journal of Power Sources, 2008, 179(1): 297-300.

[4] EIJKEL J C T, BERG A V D. Nanofluidics: what is it and what can we expect from it?[J]. Microfluidics and Nanofluidics, 2005, 1(3): 249-267.

[5] ZHU Y, ZHANG Y, SHI Y,Lubrication behavior of water molecules confined in TiO2nanoslits: a molecular dynamics study[J]. Journal of Chemical & Engineering Data, 2016, 61(12): 4023-4030.

[6] 趙夢堯, 楊雪平, 楊曉寧. 石墨烯狹縫受限孔道中水分子的分子動力學模擬[J]. 物理化學學報, 2015, 31(8): 1489-1498. ZHAO M Y, YANG X P, YANG X N. Molecular dynamics simulation of water molecule in confined graphene silt[J]. Acta Physico-Chimica Sinica, 2015, 31(8): 1489-1498.

[7] WEI M J, ZHANG L, LU L,Molecular behavior of water in TiO2nano-slits with varying coverages of carbon: a molecular dynamics simulation study[J]. Physical Chemistry Chemical Physics, 2012, 14(48): 16536-16543.

[8] ALGARA-SILLER G, LEHTINEN O, WANG F C,Square ice in graphene nanocapillaries[J]. Nature, 2015, 519(7544): 443-445.

[9] NAIR R R, WU H A, JAYARAM P N,Unimpeded permeation of water through helium-leak-tight graphene-based membranes[J]. Science, 2012, 335(6067): 442-444.

[10] JOSHI R, CARBONE P, WANG F C,Precise and ultrafast molecular sieving through graphene oxide membranes[J]. Science, 2014, 343(6172): 752-754.

[11] HUANG H, SONG Z, WEI N,Ultrafast viscous water flow through nanostrand-channelled graphene oxide membranes[J]. Nature Communications, 2013, 4(4): 2979.

[12] YANG X, YANG X, LIU S. Molecular dynamics simulation of water transport through graphene-based nanopores: flow behavior and structure characteristics[J]. Chinese Journal of Chemical Engineering, 2015, 23(10): 1587-1592.

[13] YANG X, DAI H, XU Z. Water permeation and ion rejection in layer-by-layer stacked graphene oxide nanochannels: a molecule dynamic simulation[J]. The Journal of Physical Chemistry C, 2016, 120(39): 22585-22596.

[14] KARNIADAKIS G E, BESKOK A, ALURU N. Microflows and Nanoflows: Fundamentals and Simulation[M]. Berlin: Springer Science & Business Media, 2006.

[15] TRAVIS K P, TODD B, EVANS D J. Departure from Navier-Stokes hydrodynamics in confined liquids[J]. Physical Review E, 1997, 55(4): 4288.

[16] BOT?AN A, ROTENBERG B, MARRY V,Hydrodynamics in clay nanopores[J]. The Journal of Physical Chemistry C, 2011, 115(32): 16109-16115.

[17] SENDNER C, HORINEK D, BOCQUET L,Interfacial water at hydrophobic and hydrophilic surfaces: slip, viscosity, and diffusion[J]. Langmuir, 2009, 25(18): 10768-10781.

[18] HUANG D M, COTTIN-BIZONNE C, YBERT C,Aqueous electrolytes near hydrophobic surfaces: dynamic effects of ion specificity and hydrodynamic slip[J]. Langmuir, 2008, 24(4): 1442-1450.

[19] LIU B, WU R, BAIMOVA J A,Molecular dynamics study of pressure-driven water transport through graphene bilayers[J]. Physical Chemistry Chemical Physics, 2016, 18(3): 1886-1896.

[20] HO T A, STRIOLO A. Molecular dynamics simulation of the graphene-water interface: comparing water models[J]. Molecular Simulation, 2014, 40(14): 1190-1200.

[21] BERENDSEN H J C, GRIGERA J R, STRAATSMA T P. The missing term in effective pair potentials[J]. The Journal of Physical Chemistry, 1987, 91(24): 6269-6271.

[22] HUMMER G, RASAIAH J C, NOWORYTA J P. Water conduction through the hydrophobic channel of a carbon nanotube[J]. Nature, 2001, 414(6860): 188-190.

[23] KALé L, SKEEL R, BHANDARKAR M,NAMD2: greater scalability for parallel molecular dynamics[J]. Journal of Computational Physics, 1999, 151(1): 283-312.

[24] DARDEN T, YORK D, PEDERSEN L. Particle mesh Ewald: an-log() method for Ewald sums in large systems[J]. The Journal of Chemical Physics, 1993, 98(12): 10089-10092.

[25] MA M, GREY F, SHEN L,Water transport inside carbon nanotubes mediated by phonon-induced oscillating friction[J]. Nat. Nano, 2015, 10(8): 692-695.

[26] BERNARDI S, TODD B D, SEARLES D J. Thermostating highly confined fluids[J]. The Journal of Chemical Physics, 2010, 132(24): 244706.

[27] DASILVA L B. Structural and dynamical properties of water confined in carbon nanotubes[J]. Journal of Nanostructure in Chemistry, 2014, 4(2): 1-5.

[28] NAKAMURA Y, OHNO T. Structure of water confined inside carbon nanotubes and water models[J]. Materials Chemistry and Physics, 2012, 132(2/3): 682-687.

[29] THOMAS J A, MCGAUGHEY A J. Water flow in carbon nanotubes: transition to subcontinuum transport[J]. Physical Review Letters, 2009, 102(18): 184502.

[30] 陳其樂, 孔憲, 盧滇楠, 等. 外壁荷電性質對雙壁碳納米管中水分子運動行為的影響[J]. 化工學報, 2014, 65(1): 319-27. CHEN Q L, KONG X, LU D N,Effect of charge properties of outer wall on the motion of water molecules in double-walled carbon nanotubes[J]. CIESC Journal, 2014, 65(1): 319-327.

Non-equilibrium molecular dynamics simulation of water flow inside nano-slit

NAN Yiling, KONG Xian, LI Jipeng, LU Diannan

(State Key Laboratory of Chemical Engineering,Department of Chemical Engineering, Tsinghua University,Beijing100084,China)

Flow behavior of water in nano confinement is essential to various application fields, including water purification, desalination, energy conversion, DNA sequencing,. It has been recognized that traditional hydrodynamics theory like Navier-Stokes (N-S) is no longer applicable in systems with lower dimension. To assess the limits of N-S equation, the molecular dynamics simulation is used to study water flow behavior in nano-slit. The nano-slit is formed by two parallel graphene sheets separated by a certain distance. Flow rate profiles of water in nano-slit with different distance between two graphene sheets show that when the distance between two graphene sheets is less than 3 nm, N-S equation cannot describe the flow behavior correctly. This means, N-S equation is applicable for channels with size larger than 3 nm, which is about ten times the diameter of water molecule. For pores in which N-S equation is applicable, effective viscosity and slip length were obtained by fitting the flow rate profiles with N-S equation. The influences of pore size, driving force, and wall hydrophobicity on the flow behavior were also investigated with emphasis on the effective viscosity and slip length. With the increases of slit pore size or driving force, water average velocity increases, accompanied by an increase in effective viscosity and a decrease in slip length. The increase of the wall hydrophilicity of the nano-slit results in a decrease of slip length while imposes no obvious effect on the effective viscosity.

nano-slit; molecular dynamics simulation; flow; effective viscosity; slip length; non-equilibrium molecular dynamics; surface

10.11949/j.issn.0438-1157.20161527

TQ 028.8

A

0438—1157（2017）05—1786—08

盧滇楠。

南怡伶（1994—），女，碩士研究生。

國家自然科學基金項目（21476125）。

2016-10-31收到初稿，2017-01-17收到修改稿。

2016-10-31.

LU Diannan, ludiannan@ tsinghua. edu.cn

supported by the National Natural Science Foundation of China (21476125).