黑磷納米通道內壓力驅動流體流動特性*

張忠強 劉漢倫 范晉偉 丁建寧2) 程廣貴

1) (江蘇大學, 智能柔性機械電子研究院, 鎮江 212013)

2) (常州大學, 江蘇省光伏科學與工程協同創新中心, 常州 213164)

3) (大連理工大學, 工業裝備結構分析國家重點實驗室, 大連 116024)

1 引 言

近年來, 隨著微/納尺度技術的不斷發展, 越來越多的低維材料開始走出實驗室并應用于實踐中, 低維材料[1-3]在物理、化學、電學和光學等方面所表現出來的優良性能也推動著微/納流體器件的發展, 使得微/納流體器件的生產與應用成為可能, 例如: 微/納機電系統[4-6]、微/納生物傳感器[7]和微/納芯片[8]等.這些微納流體器件[9,10]可以進入宏觀流體器件無法到達的區域, 并具有低能耗、高靈敏度和高精度等特點, 但目前關于微/納流體器件的機理研究還處于探索階段, 需要進行相應的研究和補充, 而微/納流體器件中經常涉及的問題就是流體在納米通道內的流動特性[11,12]和流-固界面機理的研究[13,14].

黑磷[15-17]作為目前比較熱門的低維材料, 不僅具有特殊的光學、電學特性, 還具有高載流子遷移率及天然帶隙等特點[18].除此之外, 由于黑磷表面特殊的褶皺結構所導致的各向異性也是一個非常值得深入研究的特性.目前關于黑磷各向異性的研究主要集中于力學[19]、電學[20]和光學[21]方面.Chen等[19]利用原子力顯微鏡測量放置在凹槽上的不同手性的黑磷納米帶的力學特性, 結果發現沿鋸齒型方向和扶手椅型方向上測得的理想彈性模量分別為65和27 GPa, 得到黑磷的力學性能會受到各向異性影響的結論.Liu等[20]將黑磷應用于場效應管中, 通過將四對金屬接觸點按間距45°進行放置, 研究各向異性對黑磷場效應管的影響, 結果發現黑磷場效應管的漏極電流會隨著角度的變化而變化.Xia等[21]分別將厚度為2-5 μm的黑磷薄膜按30°的間距進行分塊, 并進行紅外偏振處理,結果發現不同角度黑磷的波長存在很大區別, 即各向異性可以對黑磷的光學特性產生影響.而目前關于各向異性對黑磷表面上流體流動特性影響的研究卻很少, 需要進行相應的研究, 故本文利用黑磷構建納米通道, 并研究各向異性等參數對黑磷納米通道內水分子流動特性的影響, 為黑磷的進一步發展以及在微/納流體器件中的應用提供了理論研究基礎.

2 模型和方法

本文選用的材料為單層黑磷, 黑磷是由磷原子通過sp3雜化形成褶皺狀的蜂窩結構, 其特殊的褶皺結構導致的各向異性是非常值得研究的內容, 單層黑磷模型的尺寸為6.56 nm × 4.00 nm, 磷原子數為720, 參考石墨烯和碳納米管手性[22,23]的概念, 將各向異性的概念應用于黑磷, 將黑磷溝槽走向與水分子流動方向之間的夾角定義為手性角θ,考慮構建周期性模型的需要, 邊界處的磷原子需要相互匹配, 目前, 本文只構造了 0°, 37.4°, 66.6°和90°四種使得邊界處黑磷原子能夠匹配的不同手性的單層黑磷模型, 圖1(a)為單層黑磷的模型圖.為了研究各向異性等參數對黑磷納米通道內水分子流動特性的影響, 將兩個平行的單層黑磷構造成黑磷納米通道, 并在黑磷納米通道內填充水分子組成本文使用的Poiseuille流系統模型圖, 如圖1(b)所示, 其中水分子的模型選用TIP4P模型, 不同于石墨烯的是, 單層黑磷的所有原子并不在同一平面上, 而是如圖1(b)所示的兩個平面, 這里使用藍色和綠色兩種顏色代表不同平面的黑磷原子.模擬系統在x和y兩個方向的尺寸分別為6.5和4.0 nm,兩個單層黑磷之間的垂直距離即為納米通道的寬度, 用H進行表示, 對黑磷納米通道內的水分子施加一個加速度gx, 形成Poiseuille流, 研究各向異性等參數對黑磷納米通道內水分子的Poiseuille流流動特性的影響.

本文采用分子動力學模擬的方法研究各向異性等參數對壓力驅動作用下黑磷納米通道內水分子的Poiseuille流流動特性的影響.整個模擬系統中,x和y方向設置為周期性邊界,z方向設置為固定邊界.首先需要對模擬系統進行弛豫作用, 使得整個系統達到平穩狀態, 系綜為正則系綜(NVT),時間步長為1 fs, 弛豫步數為1000000, 使得模擬系統的溫度穩定在300 K.由于采用的是假設恒溫條件, 所以x和y方向采用Langevin恒溫器保持298 K的恒溫, 并且不考慮黑磷-水界面間的熱蠕變和溫度躍變.勢函數選擇方面, 選用Stillinger-Weber (SW)[23]勢函數描述黑磷原子之間的相互作用, 水分子模型選取改進版的四位TIP4PEw模型[24], 通過 particle-particle particle-mesh(PPPM)方法計算水分子之間的長程庫侖勢作用,庫侖相互作用的截止距離為12 ?, 為了保持水分子的O-H距離和H-O-H鍵角的穩定性, 這里選用SHAKE算法[25].黑磷與水分子的相互作用采用Lennard-Jones (L-J)勢函數進行描述, 作用參數可以通過Lorentz-Berthelot混合算法計算,具體參數見表1[26,27].L-J作用的截止距離為10 ?.整個系統的模擬過程為5 ns (5 × 106步), 先用1 ns的時間進行弛豫, 使模擬系統達到平衡, 接著對黑磷納米通道內的水分子施加一個加速度, 使得水分子形成穩定的Poiseuille流, 此過程持續時間為2 ns,最后2 ns用來輸出和計算處理分子動力學模擬獲得的數據, 對這些數據進行處理得到需要研究的參數.

圖1 (a)單層黑磷模型圖, 其中手性角度θ指黑磷褶皺方向與水分子流動方向夾角; (b)黑磷納米通道內水分子流動的Poiseuille流模型圖Fig.1.(a) Monolayer black phosphorus models, chiral angle θ is the intersection angle between water flow direction adjacent the top plate and the ripple direction of BP monolayer; (b) poiseuille flow model of water molecules in black phosphorus nanochannels.

表1 L-J勢能函數的參數值Table 1.Parameter values of L-J potential function.

3 模擬結果討論與分析

3.1 驅動力對黑磷納米通道內水分子的邊界滑移的影響

采用如圖1(b)所示的Poiseuille流系統模型,選用手性角度為0°的單層黑磷模型, 將兩個相同單層黑磷平行放置, 間隙設置為6 nm, 即納米通道的寬度設置為6 nm, 對黑磷納米通道內水分子沿x方向施加不同加速度(gx=1.0 × 1012, 1.5 × 1012,2.0 × 1012m/s2), 探究驅動外力對黑磷納米通道內水分子Poiseuille流流體特性的影響.將黑磷納米通道內的水分子進行分層統計, 沿納米通道寬度方向按每層0.3 nm的高度對納米通道內的水分子進行分層, 通過分子動力學方法對每層內的水分子進行平均匯總, 得到如圖2所示的在不同驅動外力作用下的黑磷納米通道內水分子的速度分布圖, 從圖中可以看出, 在不同驅動外力下, 黑磷納米通道內的水分子沿納米通道寬度方向都近似為二次曲線分布, 滿足Poiseuille流的速度分布情況, 并且在邊界處出現輕微的滑移現象.已有研究表明[28,29],在恒外力條件下的牛頓平板Poiseuille流問題中,可以在宏觀N-S方程的基礎上考慮線性邊界滑移條件, 進行描述Poiseuille流中流體速度分布問題.確有文獻[30]對壁面速度運用高階多項式進行擬合,但考慮到本模擬中邊界處數據曲線相對于擬合曲線并未出現太大波動, 而且本文重點在于研究黑磷各向異性的影響, 故采用假定線性邊界條件進行處理.運用Origin軟件對圖中的曲線進行非線性擬合處理, 進而得到各曲線關于水分子流動速度與位置坐標之間擬合的二次曲線方程式, 通過將邊界處水分子(最外層水分子)的位置坐標代入方程式求得邊界滑移速度, 三種驅動外力對應計算得到的滑移速度Vs分別為6.933, 9.956和13.491 m/s, 顯而易見, 隨著驅動外力增加, 黑磷納米通道內邊界處水分子滑移速度也隨之增加, 從而對黑磷納米通道內水分子的流動特性產生影響.

圖2 模型手性角度為0°時水分子的速度分布Fig.2.The velocity distribution of water molecules when the chiral angle of the model is 0°.

數密度分布曲線可以清晰地顯示出納米通道內流體的分布情況, 對了解流體的流動特性具有很大幫助.同速度分布曲線的獲得方法相同, 數密度分布曲線同樣需要對黑磷納米通道內水分子進行分層劃分, 分別對四組模型內水分子按照層高為0.01 nm進行分層, 通過分子動力學模擬對每層的水分子數目進行統計, 輸出數據, 將每層獲得的數據除以所屬層的體積, 得到黑磷納米通道內水分子沿納米通道寬度方向的數密度分布圖, 如圖3所示, 圖中y軸數密度的單位代表著每平方納米內水分子的數目.可以看出, 曲線大致分為兩個部分,第一部分為兩邊靠近邊界處的曲線, 可以看出此時水分子的數目受到黑磷壁面的影響, 分布情況較為不均, 邊界處水分子發生密度堆積和分層現象, 出現“類固體”效應, 數密度分布曲線的第一個峰值即為邊界層, 而隨著距壁面距離的增加, 水分子的分布情況逐級平穩, 最后達到穩定狀態, 即圖中的第二部分--平穩曲線.從圖3的密度分布圖還可以看出, 壁面處的確出現了密度漲落, 從第一個峰值的放大圖可以看出, 水分子的數密度分布出現了輕微的差異, 但這種差異對結果影響很小, 而且水分子的速度曲線在邊界處的波動較小, 剪切應變率波動也就很小, 根據公式(τ為剪應力,μ為剪切黏性系數,為剪切應變率), 邊界處與通道中央處水分子的黏性差異較小.并且, 邊界層占通道寬度的比例較小, 因此本文仍采用通道寬度方向上的均勻化近似, 并用線性邊界滑移條件進行了處理.

圖3 黑磷納米通道內水分子沿通道寬度方向的數密度分布圖Fig.3.Number density distribution of water molecules along the channel width in the black phosphorus nanochannels.

3.2 各向異性對壓力驅動作用下黑磷納米通道內水分子流動特性的影響

在納米通道寬度不變的情況下, 研究各向異性對壓力驅動作用下黑磷納米通道內水分子Poiseuille流流動特性的影響, 分別使用手性角度分別為 0°, 37.4°, 66.6°和 90°的單層黑磷模型構造黑磷納米通道, 在通道內填充水分子組成四個不同的黑磷-水模擬系統, 分別對四個模擬系統中的水分子施加三種不同加速度(gx=1.0 × 1012, 1.5 × 1012,2.0 × 1012m/s2), 通過分子動力學方法模擬計算,研究各向異性對壓力驅動作用下黑磷納米通道內水分子流動特性的影響.

3.2.1 各向異性對壓力驅動作用下黑磷納米通道內水分子滑移特性的影響

如3.1節所述, 為了研究各向異性對邊界滑移特性的影響, 需要獲得納米通道內水分子沿納米通道方向的速度分布圖, 針對四種不同手性的模擬系統, 分別對黑磷納米通道內的水分子施加不同的加速度(gx=1.0 × 1012, 1.5 × 1012, 2.0 × 1012m/s2),得到如圖4所示的速度分布圖(圖4(a)-圖4(c)包含三種手性角, 圖2為一種), 綜合比較可以看出, 從四個不同手性模擬系統中計算得到的速度分布圖都近似為二次曲線分布, 符合納米通道內Poiseuille流的速度分布圖特征, 對四幅圖分別進行非線性擬合后得到在施加不同加速度下對應的四組非線性曲線, 使用與3.1節相同的方法, 分別計算得到不同手性模擬系統在施加不同加速度的條件下邊界處水分子的滑移速度VS(m/s), 如表2所列.從表中數據可以看出, 在不同手性模擬系統中, 隨著加速度值的增加, 邊界滑移速度都呈現增加的趨勢; 當施加的加速度值相同時, 隨著手性角度的增加, 黑磷納米通道內邊界處的水分子的滑移速度也呈現減小的趨勢.圖4(d)給出了流-固界面相互作用勢能云圖, 可以看出, 黑磷表面特殊的褶皺結構, 導致了其表面規則、粗糙的勢能分布, 從而影響了邊界處水分子的流動特性, 這是造成黑磷表面各向異性流動的內在原因[31].

表2 不同手性情況中, 不同加速度對應的水分子邊界滑移速度VS統計表Table 2.Statistical table of water molecule boundary slip velocity VS corresponding to different accelerations in different chiral conditions.

3.2.2 各向異性對壓力驅動作用下黑磷納米通道內水分子黏度系數的影響

納米通道內流體的黏度系數也是反映納米通道內Poiseuille流的一個重要參數, 尤其在許多的物理傳輸過程中扮演著非常重要的角色, 因此流體黏度系數的計算和研究也非常重要.不同于納米通道內Couette流中使用剪切應力比上剪切應變率的方法, 在納米通道內Poiseuille流的黏度系數通常通過對速度分布曲線擬合計算得到, 從圖4中可以發現黑磷納米通道內水分子的速度分布曲線近似為二次分布, 符合Poiseuille流特征.根據Stokes方程[32]可以得到關于水分子流動速度u(z)和位置坐標z之間的關系式

式中z1,z2為黑磷壁面坐標,μ為黏度系數.由于水分子密度參數ρ可以在模擬計算中得到, 加速度g為預先設定值, 水分子沿納米通道寬度方向速度分布關系式u(z)可以通過速度分布曲線擬合得到,故可以利用(1)式計算出水分子的黏度系數, 通過對公式兩邊求二次導數并變形可以得到

圖4 速度分布圖及勢能云圖 (a) 模型手性角度為37.4°時水分子的速度分布; (b) 模型手性角度為66.6°時水分子的速度分布;(c) 模型手性角度為90°時水分子的速度分布; (d)模型手性角度為90°時的勢能分布云圖Fig.4.Velocity distribution diagram and potential energy cloud diagram: (a) The velocity distribution of water molecules when the chiral angle of the model is 37.4°; (b) the velocity distribution of water molecules when the chiral angle of the model is 66.6°;(c) the velocity distribution of water molecules when the chiral angle of the model is 90°; (d) potential energy cloud diagram when the chiral angle of the model is 90°.

對速度分布曲線擬合得到的曲線ux(z)可以描述為

其中u0和a可以在擬合中直接得到,z0為通道中心坐標, 故d2ux/dz2=2a, (2)式可以最終轉變為

針對不同手性的模擬系統, 將不同的加速度值以及擬合速度分布曲得到的參數a代入(4)式中就可以計算得到相對應的黏度系數μ(mPa·s), 如表3所列.從表中可以發現, 不同手性的模擬系統在不同加速度條件下, 黑磷納米通道內水分子的黏度系數幾乎保持不變.為了更直觀地表現出這一結果, 將表中數據匯總繪制出黏度系數方差圖, 如圖5所示, 從圖中可以清晰地看出, 隨著手性角度的增加, 黑磷納米通道內水分子的黏度系數幾乎為一條直線, 不同手性的模擬系統中水分子的黏度系數值差別很小, 由此可以得出水分子的黏度系數值基本不隨手性角度變化而變化的結論.

表3 不同手性的模擬系統在不同加速度條件下的水分子黏度系數μ分布表Table 3.Distribution of water molecular viscosity coefficient μ of simulation systems with different chirality under different acceleration conditions.

3.3 通道寬度對壓力驅動作用下黑磷納米通道內水分子流動特性的影響

納米通道寬度也是影響納米通道內Poiseuille流流動特性的重要參數, 為了研究納米通道寬度的影響, 分別構建了通道寬度分別為H=3, 4, 5, 6 nm的四個模擬系統, 其中單層黑磷選用手性角度為90°的黑磷模型(即黑磷褶皺方向垂直于水分子的流動方向).四個模擬系統中的水分子被施加的加速度值gx都為2.0 × 1012m/s2, 通過分子動力學計算得到如圖6的速度分布圖.從圖中可以看出,不同孔寬模型內水分子的流動符合Poiseuille流動特征, 并且都在邊界處出現滑移現象, 通過對圖中曲線分別進行非線性擬合, 結合邊界處水分子坐標計算得到相應的滑移速度, 如表4所列, 從表中可以發現, 隨著通道寬度的增加, 通道內水分子的邊界滑移速度隨之增加.

圖5 不同手性的模擬系統在不同加速度條件下的水分子黏度系數方差分布圖Fig.5.Variance distribution of water molecular viscosity coefficient of a simulation system with different chirality under different acceleration conditions.

表4 不同納米通道寬度內水分子的邊界滑移表Table 4.Boundary slip of water molecules at different nanochannels widths.

圖6 不同納米通道寬度內水分子沿通道寬度方向速度分布圖Fig.6.Velocity distribution of water molecules along the width of different nanochannel widths.

3.4 黑磷層數對壓力驅動作用下黑磷納米通道內水分子流動特性的影響

黑磷層數也是研究黑磷納米通道內Poiseuille流中需要考慮的一個參數, 黑磷層與層之間的間隔為0.53 nm, 本分子動力學模擬中所選用的L-J勢函數的截斷半徑為1 nm, 故當黑磷層數超過兩層時, 最外層黑磷片對納米通道內的水分子的影響幾乎可以忽略不計, 故本節在單層黑磷納米通道的基礎上構造雙層黑磷片組成的納米通道, 在其中填入水分子構成Poiseuille流模擬系統, 并考慮其對納米通道內水分子流動特性的影響.

選取手性角度分別為0°和37.4°的黑磷片構建黑磷納米通道, 同樣對納米通道進行分層統計, 獲取水分子沿納米通道寬度方向的速度分布曲線, 并與單層黑磷納米通道模型系統進行對比, 得到如圖7所示的速度分布圖.可以看出, 當對納米通道內水分子施加相同加速度時, 黑磷層數對水分子的流動影響較小, 單層黑磷納米通道和雙層黑磷納米通道內的水分子的速度分布曲線存在輕微差異, 整體趨勢基本一致.比較兩種情況下納米通道內水分子的邊界滑移速度、黏度系數以及交互界面能等參數, 如表5所列, 結果表明雙層模型相對于單層模型邊界滑移速度出現輕微減小, 黏度系數基本保持不變, 通過比較兩種模型的交互界面能, 發現雙層模型與水分子的界面能較單層模型輕微增加, 因此邊界滑移速度的輕微減小是由于交互界面能的增加所造成的, 通過比較不同手性模型的模擬結果,不難發現黑磷的各向異性不受黑磷層數的影響.

圖7 不同層數模型對應的速度分布圖Fig.7.Velocity distributions corresponding to different layer models.

表 5 不同黑磷層數納米通道模型中流固界面參數對比Table 5.Comparison of the interfacial parameters for the models with different BP layers.

4 結 論

以黑磷材料構建納米通道, 基于Poiseuille模型研究了各向異性等參數對壓力驅動作用下黑磷納米通道內水分子的滑移及黏度特性的影響, 探究了各向異性對黑磷-水流固界面的影響作用, 揭示了驅動力、各向異性、納米通道寬度和黑磷層數等對納米通道內水分子Poiseuille流流動特性的影響規律.研究結果表明: 隨著驅動力的增加, 邊界滑移速度隨之增加; 各向異性也會對壓力驅動作用下納米通道內的水分子的流動特性產生影響, 具體表現為邊界滑移速度會隨著手性角度的增加而減小,而水分子黏度系數卻不受各向異性的影響; 最后,在加速度值保持不變的情況下, 研究納米通道寬度和黑磷層數對水分子流動特性的影響, 發現隨著納米通道寬度增加, 水分子滑移速度減小; 雙層模型中水分子的速度分布和沿通道寬度方向的數密度分布情況依舊與單層模型差異很小, 雙層黑磷與水分子之間交互界面能相對于單層黑磷也呈現增加的趨勢, 各向異性規律依然保持不變.