壓力驅(qū)動流體納米管膜傳遞的分子動力學(xué)模擬

馮 劍

1 前言

人工固態(tài)納米孔膜以及生物管道中流體的流動揭示了很多新的現(xiàn)象，如碳納米管中流體的快速流動，非線性傳遞以及或氮化硼管道的大滲透性等[1]。近十幾年來，該領(lǐng)域備受研究者關(guān)注，有了眾多實驗[2-3]、理論和計算機(jī)模擬方面[4-6]的研究工作，并得出了很多研究結(jié)果，尤其是水在碳納米管或其他納米管中出現(xiàn)的流動增強(qiáng)現(xiàn)象。人們更感興趣是這些現(xiàn)象的起因以及如何控制，這方面依然存在較大的爭議[7]。目前人們認(rèn)為流體通過納米管膜體系傳遞主要存在兩種阻力，一種是流體力學(xué)阻力，另一種是熱力學(xué)阻力[8]。流體力學(xué)阻力由流體在體相到納米管口流線發(fā)生變形而引起。由于空間位阻，流體在納米管口的流體性質(zhì)會發(fā)生一定變化。熱力學(xué)阻力是因為流體由體相到納米管發(fā)生熵減小，從而造成該過程是自由能增大的非自發(fā)過程。如何區(qū)分這兩類阻力以及識別它們在納米管膜傳遞中發(fā)揮的作用，目前依然是待解難題。

流體在納米管膜體系的傳遞是一個非平衡態(tài)問題，很多性質(zhì)與模型選擇和初始條件設(shè)定關(guān)系密切，從而造成很多理論和模擬結(jié)果間存在一定的差異性。如驅(qū)動流體在納米管膜體系傳遞就有多種模型，一種是為每個流動粒子施加固定力的“重力模型”，一種是僅對儲液池某區(qū)域粒子施加作用力[7]，還有是在納米管膜體系某處設(shè)置具有選擇性的反射壁等[9]。

當(dāng)前的計算機(jī)模擬研究中，也有針對流體從儲液池進(jìn)入納米管這一特定過程的傳遞模擬[10]，但依然很少關(guān)注在該過程中驅(qū)動壓力和體系結(jié)構(gòu)的變化。本文使用分子動力學(xué)研究流體在活塞驅(qū)動下以較低的流速通過納米管膜體系時體系發(fā)生的現(xiàn)象。推動活塞的壓力驅(qū)動模型是真實的物理過程，從而避免了對體系粒子設(shè)置不適宜的驅(qū)動力而影響模擬結(jié)果。雖然活塞驅(qū)動的連續(xù)傳遞模型較難構(gòu)建，較難適用于長時間尺度的連續(xù)模擬。但對于短時間的流體模擬是有效的，它能夠用于流體由體相區(qū)到納米管等介觀相區(qū)的阻力識別。

2 模型與方法

本文使用的是分子動力學(xué)模擬是基于Langevin的恒溫方法，方程形式如下：

其中Ui是粒子i與其它所有粒子相互作用能，相互作用形式采用Lennard-Jones勢能來表示。粒子的作用參數(shù)采用粗?；腗ARTINI力場中的相關(guān)參數(shù)[11]，各物理量單位也與MARTINI力場選擇一致。是摩擦系數(shù)，Wi是隨機(jī)力。

圖1 納米管膜體系傳遞示意圖

為了研究流體在圓柱形納米管中的流動現(xiàn)象，本文采用圓柱形納米管膜系統(tǒng)。模擬盒子xyz三邊長度分別是9.4nm，9.4nm和32.9nm。在納米管中心z軸方向放置一個長度為10的納米管，其中=0.47nm。納米管由MARTINI力場中的Nda類型粒子組成。納米管直徑分別為4、5、6、7，直徑是納米管粒子的內(nèi)緣距離。在納米管兩個邊緣，即5處有一個虛擬分隔膜，其虛擬粒子勢能參數(shù)也采用Nda類型。由于這個分隔膜在本文模擬中主要其分隔作用，它與溶劑粒子的作用勢能形式依然采用Lennard-Jones勢能形式，而不是其積分勢。本文采取的算法是，當(dāng)溶劑粒子與虛擬壁面作用時，在其與虛擬壁面垂直方向，處于虛擬壁面內(nèi)部距離壁面0.5有一個虛擬粒子。溶劑粒子參數(shù)采用MARTINI力場中的P4類型，它是該力場下的水分子模型。根據(jù)MARTINI力場，溶劑粒子間的相互作用參數(shù)為5.0 kJ/mol，溶劑與納米管及虛擬粒子間的相互作用參數(shù)為4.0 kJ/mol。由于納米管粒子固定，納米管粒子間無相互作用。溶劑粒子間的相互作用參數(shù)小于其與納米管和虛擬粒子的相互作用，這表明納米管及虛擬壁面是疏水的。

整個體系被納米管及膜分為三個區(qū)域，納米管一側(cè)充滿溶劑粒子的稱之為上游儲液池，另一側(cè)真空區(qū)為下游儲液池。10200個溶劑粒子隨機(jī)插入上游儲液池。模擬盒子z方向邊緣也使用兩個虛擬壁面，其中上游的虛擬壁面可以移動，它是驅(qū)動流體流動的活塞，對壁面附近的粒子施加作用力。下游的虛擬壁面始終固定。這兩個虛擬壁面的虛擬粒子類型依然是Nda型，兩個虛擬壁面間沒有相互作用。

3 結(jié)果與討論

圖2 四個體系活塞驅(qū)動壓力與時間的關(guān)系

四種不同直徑的納米管膜體系，通過勻速向下游移動虛擬活塞，統(tǒng)計出不同時刻溶劑粒子施加在活塞上的力，其反作用力就是活塞施加的力，進(jìn)而根據(jù)活塞面積計算出活塞的驅(qū)動壓力。圖2是四條曲線分別對應(yīng)四種直徑的納米管膜體系中活塞的驅(qū)動壓力與時間的變化關(guān)系。四條曲線自上至下分別對應(yīng)直徑為4、5、6和7的納米管。該四條曲線均表明，在活塞持續(xù)驅(qū)動過程中，驅(qū)動壓力越來越大，達(dá)到一定時間后，驅(qū)動壓力會在一定值附近漲落，基本保持穩(wěn)定。納米管徑越大，這個穩(wěn)定的壓力值越小。

圖2中的四條曲線均能看出，在穩(wěn)定的壓力前，驅(qū)動壓力隨時間的變化關(guān)系并非是單調(diào)增加的，均多次出現(xiàn)壓力達(dá)到某值后又急劇降落的情形，并在低到一定值后又緩慢升高，整條曲線呈現(xiàn)鋸齒狀。要解釋這種情況本文統(tǒng)計了不同時刻進(jìn)出納米管的粒子數(shù)目。

圖3 驅(qū)動壓力與納米管內(nèi)粒子數(shù)與時間的關(guān)系

四種直徑納米管內(nèi)粒子數(shù)與時間的變化關(guān)系相似。圖3是直徑為4的納米管中的粒子數(shù)(上面的曲線)以及驅(qū)動壓力(下面的曲線)隨模擬時間的變化關(guān)系。從該圖可以看出，納米管中的粒子數(shù)隨著時間增大，逐漸增大。但粒子數(shù)出現(xiàn)明顯的波動現(xiàn)象，即不斷發(fā)生粒子進(jìn)入納米管和返回上游儲液池的現(xiàn)象。并且，納米管中的粒子數(shù)與模擬時間的曲線出現(xiàn)多個平臺現(xiàn)象，即在一定時間內(nèi)納米管中的粒子數(shù)保持相對穩(wěn)定，這個平臺對應(yīng)不同的體系狀態(tài)。比較圖中兩條曲線發(fā)現(xiàn)，當(dāng)納米管內(nèi)粒子數(shù)波動較大時，驅(qū)動壓力的波動也較大。

進(jìn)一步通過繪制體系的粒子結(jié)構(gòu)圖發(fā)現(xiàn)，在模擬初期就有少量溶劑粒子進(jìn)入納米管，但由于體相流體的吸引作用，少量粒子并不會脫離液面。隨著驅(qū)動流體的活塞不斷移動，進(jìn)入納米管流體的粒子數(shù)增加，進(jìn)而會發(fā)生蒸發(fā)和凝聚現(xiàn)象，如在t=280ps時觀測到多個粒子脫離液面，在360ps時這些粒子又返回液面，見圖4所示。流體從液面蒸發(fā)再凝結(jié)到液面的現(xiàn)象反復(fù)出現(xiàn)。這兩個時間的壓力分別是443.0bar和316.4bar，即溶劑分子從表面蒸發(fā)會造成壓力降低，而溶劑分子向表面凝結(jié)壓力會上升，這是因為粒子在液體表面凝結(jié)加劇了納米管道中流體的擁擠程度，使得上游溶劑分子進(jìn)入納米管更加困難。

圖4 t=280ps(a)和t=360ps時體系的結(jié)構(gòu)

另外，在t=12120ps時壓力有一個峰值(5525.5bar)，而在12240ps時壓力降到一個較低值(3789.8bar)。在這個時間段內(nèi)納米管內(nèi)粒子數(shù)最大值為99個，而最小值有93個，也就是說最多有6個粒子回流到上游儲液池。當(dāng)粒子回流到上游儲液池造成上游密度增大，使得壓力增大。正是由于納米管內(nèi)粒子的回流和納米管道內(nèi)蒸發(fā)的氣體分子在液體表面凝聚造成了驅(qū)動壓力的增大。

進(jìn)入納米管的粒子之所以能回流到上游儲液池熵因素發(fā)揮了很大的作用。粒子由儲液池進(jìn)入納米管是熵值減小的過程，該過程是非自發(fā)的，需要在外力做功的情況下才能實現(xiàn)。也就是說，溶劑粒子由上游儲液池進(jìn)入納米管遭受更大的阻力，這種因熱力學(xué)因素產(chǎn)生的阻力稱作熱力學(xué)阻力。圖5給出了上游儲液池(上面曲線)和直徑為4的納米管中(下面曲線)溶劑濃度隨時間的變化關(guān)系。由該圖可以看出在傳遞過程中，上游的儲液池中溶劑的濃度始終高于納米管內(nèi)的濃度。由此可以看出，熱力學(xué)阻力并非因濃度差而產(chǎn)生。僅從濃度考慮，似是納米管內(nèi)溶劑的化學(xué)勢小于上游儲液池的化學(xué)勢能。事實上，納米管內(nèi)的溶劑系統(tǒng)是一個非均相系統(tǒng)，溶劑化學(xué)勢還要考慮納米管壁面這個外場的作用。

圖5 上游儲液池與納米管內(nèi)溶劑濃度與時間的關(guān)系

圖6 上游儲液池與納米管內(nèi)粒子數(shù)與時間的關(guān)系

隨著流體傳遞過程的持續(xù)進(jìn)行，溶劑粒子通過納米管并在下游管口側(cè)出現(xiàn)溶劑堆積，溶劑粒子從納米管進(jìn)入下游儲液池將受到更大的阻力，從而在納米管內(nèi)出現(xiàn)更多的堆積粒子，納米管內(nèi)的溶劑的密度也會增大。圖2的壓力以及圖5上游儲液池溶劑濃度在模擬后期基本維持不變，顯示了流體通過納米管基本達(dá)到穩(wěn)流狀態(tài)。這也可以從圖6的上游儲液池溶劑粒子數(shù)在模擬后期對時間的線性衰減看出。通過對圖6后期的數(shù)據(jù)擬合發(fā)現(xiàn)，這個穩(wěn)流的流量約為1.8105mol m-2s-1。圖2模擬后期基本穩(wěn)定的壓力正是在相應(yīng)納米管膜體系中維持穩(wěn)流所需要的壓力。在相同推動速度下，納米管徑越大達(dá)到穩(wěn)定流動需要的推動力越小。

圖7給出了四個體系后期的穩(wěn)定壓力與納米管徑三次方倒數(shù)的變化關(guān)系。從該圖可以發(fā)現(xiàn)，穩(wěn)定壓力與納米管直徑立方的倒數(shù)有較好的線性關(guān)系?？紤]到，本文模擬體系的下游是真空狀態(tài)，這個穩(wěn)定壓力就是流體上下游的壓力差，這個線性關(guān)系與文獻(xiàn)結(jié)果吻合[10]。

圖7 納米管體系穩(wěn)定壓力與管徑的關(guān)系

4 結(jié)論

本文使用分子動力學(xué)研究基于壓力驅(qū)動的流體在不同直徑納米管膜中的傳遞過程。在壓力驅(qū)動過程中，壓力的漲落主要來自于流體粒子的回流。在流體流經(jīng)納米管之初，處于表面分子的蒸發(fā)對壓力漲落也有貢獻(xiàn)。直徑越大的納米管驅(qū)動壓力越小，驅(qū)動壓力與納米管直徑三次冪指數(shù)的倒數(shù)成正比。處于納米管內(nèi)流體的密度要小于上游儲液池的密度，造成流體粒子向上游回流可以認(rèn)為是熵在起作用。