正滲透過程中的分子動力學模擬研究

黃政宇，周寧玉,2，謝朝新,2，郭 豪

?

正滲透過程中的分子動力學模擬研究

黃政宇1，周寧玉1,2，謝朝新1,2，郭 豪1

（1. 后勤工程學院 國防建筑規劃與環境工程系，重慶 401311； 2. 重慶市水工業與環境工程技術研究中心，重慶 401311）

采用分子動力學模擬方法，利用Materials Studio 6.0軟件模擬計算了H2O、Na+、Cl-在正滲透膜內的擴散系數。同時模擬研究了三種不同環境溫度、溶液濃度下，三種粒子在正滲透膜內擴散系數的變化規律。實驗結論與模擬結果一致，即溶液濃度不變的情況下，擴散系數隨著溫度的升高而增大；環境溫度不變的情況下，擴散系數隨著溶液濃度的增大而減小。

正滲透膜；擴散系數；分子動力學；MS模擬

正滲透過程因具有能耗低、膜污染輕、設備簡單以及對許多污染物截留效率高等優點，在海水淡化、飲用水處理和廢水處理中具有良好的應用前景[1]。水通量是衡量正向滲透過程的重要指標之一，而水通量的大小取決于水分子在膜中的擴散系數，為此研究人員對擴散系數進行了廣泛的研究。目前，擴散系數主要通過實驗方法測定，但實驗條件要求嚴格，操作時間較長，對設備的要求很高，測定難度較大。隨著當前蒸蒸日上的計算機技術，利用分子模擬技術探究微觀世界漸漸成為一種熱潮。計算機模擬不僅能提供定性的描述，而且能模擬出高分子材料結構與性能的定量的結果。[2]采用分子動力學模擬正向滲透過程中粒子的擴散系數可有效克服擴散系數測試方法中存在的難題。

分子動力學（Molecular Dynamics）[3]是常用的分子模擬方法之一，用于模擬分子體系內各分子的運動。通過求解牛頓力學方程，得到體系中各粒子微觀狀態隨時間的變化，接著講得到的微觀狀態量對時間做平均量計算，即可得到體系的一系列宏觀性質及各粒子的空間分布情況。目前，研究人員利用分子動力學對擴散過程進行了研究，如Kotelyanskii和Wagner等[4,5]研究了在加NaCl和不加NaCl兩種情況下，水在水合的無定形聚酰胺膜中的擴散機理,并計算了水在膜中的擴散系數;劉清芝[6]采用分子動力學模擬方法研究了水和鹽分子在反滲透復合膜內的擴散過程，得到了結構單體與膜分離性能的關系;陶長貴[7]探究了不同聚合度的聚丙烯內氧氣的擴散過程，同時分析了其擴散的機理。但利用分子動力學模擬正向滲透過程中粒子的擴散還鮮有報道。本文通過分子動力學方法，模擬研究了在不同溫度、濃度下，水分子、鈉離子和氯離子在正滲透膜中的擴散，為研究正向滲透擴散系數提供了有效的方法，探索水溶液特性對正向滲透過程的影響規律。

1 模擬方法

1.1 晶格的構建

本文采用Accelry公司的MS（Materials Studio 6.0）軟件，研究水分子、鈉離子、氯離子在正滲透膜三醋酸纖維膜（CTA）內的擴散。首先，在MS Visualizer環境中初步構建小分子（水分子、鈉離子和氯離子）及正滲透膜（CTA），其中膜包含一條高分子鏈（含10個單體），然后采用Discover模塊內的“Smart Minimizer”，對小分子和該高分子鏈進行能量最小化以獲得最穩定的分子構型，進行能量最小化時，收斂水平一般使用“Fine”。得到每個分子的優化結構后，利用“Amorphous Cell”將高分子鏈和小分子構建擴散體系模型。每個體系模型內僅放置一個CTA高分子鏈以減少鏈段效應帶來的相互作用力影響。邊界條件采用周期性邊界條件。在構建晶格時還需設置溶液的溫度和模型的密度以及幾何構型數量。圖1為CTA膜在溫度為298 K，濃度為1%的NaCl溶液中的模型結構，即模擬模型內包含1個NaCl分子和285個水分子，以及一個膜的高分子鏈，幾何構型數量為10，最終密度為1.0 g/mL。用同樣的方法構建不同溫度、濃度下的模型結構，用于后續研究。

1.2 分子力場

在整個模擬過程中，模擬皆是基于分子力學和分子力場的DISCOVER模塊實現的。所有的計算都是在COMPASS力場下進行的。該力場把有機分子體系的力場與無機分子體系的力場相統一，能夠應用于有機和無機小分子、高分子，一些金屬離子、金屬氧化物和金屬的模擬研究[8]。

該力場對于振動頻率的計算誤差比其他模塊要小至少50%。力場中的非鍵合勢是由范德華勢（vdW）和靜電相互作用勢兩種勢能組合而成。其中范德華勢的描述為[9]：

靜電相互作用表示分子中各原子靜電荷的庫倫相互作用對勢能的貢獻，被描述成：

1.3 模擬設置

在模擬過程中，選擇“atom-based”計算范德華相互作用力，同時，計算靜電相互作用力時選擇“Ewald”法。本文研究的體系均為小分子體系，因而計算兩種作用力時以原子基團為基礎，即選用“Group-based”模塊。這樣不僅對結果精確度的影響可以忽略不計，而且大大減少了模擬機時。

圖1 298 K，溶液濃度1%的晶格模型

在模擬不同NaCl溶液濃度對水分子、鈉離子、氯離子在CTA膜內的擴散系數的影響時。濃度1代表晶格當中包含285個水分子和1個NaCl分子，濃度2代表晶格當中包含285個水分子和2個NaCl分子，濃度3代表晶格當中包含285個水分子和3個NaCl分子

在濃度不變的情況下，改變溶液的溫度（298，308，318 K），研究不同NaCl溶液溫度對水分子、鈉離子、氯離子在CTA膜內的擴散系數的影響。

晶格構建成功后，利用“Forcite”模塊，對于獲得的擴散模擬模型，首先使用“Geometry Optimization”進行結構優化，接著選擇“Anneal”進行模擬退火，得到穩定的晶格模型。之后，在“Dynamica”內分別進行20 ps的NVT及20 ps的NPT動力學模擬，達到弛豫結構的效果。在模擬過程的最后進行40ps的NVT動力學模擬，輸出內容選擇“ALL”，不僅得到平衡后的晶格模型結構，同時得到分子的運動軌跡，為接下來的分析計算做好鋪墊。

2 分析與討論

2.1 分析與討論

在單位時間內通過垂直于擴散方向的單位截面積的擴散物質流量，即擴散通量J（Diffusion flux）與該截面處的濃度梯度成正比，其數學表達式[10]為：

式中，稱為擴散系數，為擴散物質的體積濃度，d/d為濃度梯度，“—”號表示擴散方向與濃度梯度方向相反。擴散系數相當于濃度梯度為1的擴散通量，由此可見擴散運動會隨著值的增大而變得快速。擴散系數的計算如下：

稠密流體中單一分子的運動軌跡不是一個簡單單一的，而是近似于數學上的隨機漫步。這種路徑由愛因斯坦在對布朗運動的研究中進行了著名的分析，發現做隨機漫步的粒子的移動距離的平方的平均數與時間成正比。[11]模擬研究中，利用分子動力學對擴散系數進行模擬計算。動力學模擬過程記錄粒子的運動軌跡，通過“Forcite”模塊中的“Analysis”對上述得到的運動軌跡進行利用，在“Dynamic”分類里可以選擇分析“mean squared displacement”，得出粒子的均方位移曲線MSD。

本文計算擴散系數的方法就是利用MSD曲線進行計算。利用軟件生成的MSD對時間的曲線，應用線性擬合的方法，將曲線擬合成一條直線，記為=+，其中斜率為。斜率表示MSD對時間微分的比率，可以近似于式（4）中的微分。由于MSD的值在定義上已經對擴散原子做了平均，因此公式（4）可以簡化為：=/6。

圖2為298 K，濃度1的水分子、鈉離子與氯離子在CTA膜內的均方位移曲線：

圖2 H2O、Na+、Cl-均方位移曲線

在計算出擴散系數之前，要先確認即將得到的擴散系數是不是可靠有效的。確認的方法可以采用判斷函數lg（MSD）與lg的斜率值的大小。<1時，屬于異常擴散；=1時屬于愛因斯坦擴散[9]。

當發生愛因斯坦擴散時，我們可以說得到的擴散系數是有效可靠的。圖3是上述均方位移曲線圖的對數曲線圖。圖中，水分子的值一開始小于1，接著無限接近于1。可知水分子的擴散發生了愛因斯坦擴散。同理，鈉離子和氯離子的都發生了愛因斯坦擴散，由此接下來計算得到的擴散系數是可靠地。在本文的后續的模擬研究中，水分子、鈉離子、氯離子的擴散都發生了愛因斯坦擴散，計算出的擴散系數都有效。

對得到的均方位移曲線，利用前面介紹的方法進行計算，可以得到CTA膜內鈉離子的擴散系數為0.033 8 A2/ps，氯離子的擴散系數為0.056 5 A2/ps，水分子的擴散系數為0.188 8 A2/ps。

圖3 H2O、Na+、Cl-的MSD-t對數曲線

2.2 環境溫度、NaCl溶液濃度對擴散系數的影響

圖4為在三種溫度（298、308、318 K）和三種NaCl溶液濃度（濃度1、濃度2、濃度3）下，水分子、鈉離子、氯離子在CTA膜中的擴散系數。

（a）H2O

（b）Na+

（c）Cl-

圖4（a）為三種不同環境溫度、NaCl溶液濃度下，水分子在CTA膜內的擴散系數。溶液濃度為濃度1的情況下，溫度為298 K時水分子的擴散系數為0.188 8 A2/ps，溫度由298 K升至308 K時，擴散系數由0.188 8 A2/ps變為0.2054 A2/ps，溫度升至318 K時，擴散系數擴大至0.248 2 A2/ps。在溶液濃度為濃度2的情況下，對應三種溫度（298、308、318 K）下，水分子在CTA膜內的擴散系數分別為0.176 1、0.186 2、0.239 4 A2/ps。濃度3時，對應擴散系數也有著同樣增大的趨勢。由此可見，在溶液濃度一定的情況下，水分子在正滲透膜（CTA膜）內的擴散系數隨溫度的升高而增大。在溫度為298 K的情況下，溶液濃度從濃度1升至濃度2、濃度3時，對應的擴散系數從0.1888 A2/ps減小至0.1761、0.1691 A2/ps。溫度為308 K時，三種不同濃度對應的擴散系數分別為0.205 4、0.186 2、0.185 7 A2/ps，逐漸減小。溫度為318 K時，擴散系數的變化有著同樣的趨勢。由此可見，在環境溫度條件一定時，溶液的濃度越高，水分子在正滲透膜內的擴散系數越小。首先，這是與理論相符合的。其次，通過相關實驗驗證，當低濃度的鹽溶液作為原料液時，在宏觀方面也表現為，隨著溫度的增大或濃度的降低，水通量也相應隨之增大[12]。

圖4（b）為三種不同環境溫度、NaCl溶液濃度下，Na+在CTA膜內的擴散系數。溶液濃度為濃度1的情況下，溫度為298 K時Na+的擴散系數為0.033 8 A2/ps，溫度上升至308 K和318 K時對應的擴散系數為0.038 1和0.046 1 A2/ps。在溶液濃度為濃度2的情況下，對應三種溫度（298、308、318 K）下，Na+在CTA膜內的擴散系數分別為0.028 1、0.035 3、0.039 5 A2/ps，逐漸增大。濃度3時，對應的三種溫度下Na+的擴散系數有著同樣的變化趨勢。由此可見，在溶液濃度一定的情況下，Na+在正滲透膜（CTA膜）內的擴散系數隨溫度增大而增大。在溫度為298 K的情況下，三種溶液濃度（濃度1、濃度2、濃度3）對應的Na+擴散系數由0.033 8 A2/ps減小至0.028 1、0.020 1 A2/ps。溫度為308 K時，三種不同濃度對應的擴散系數分別為0.038 1、0.035 3、0.029 6 A2/ps，不斷減小。溫度為318 K時，三種不同濃度對應的擴散系數有著同樣的趨勢。由此可見，在環境溫度條件一定時，溶液的濃度越高，Na+在正滲透膜內的擴散系數越小。

圖4（c）為三種不同環境溫度、NaCl溶液濃度下，Cl-在CTA膜內的擴散系數。溶液濃度為濃度1的情況下，溫度由298 K增大到308和318 K，Cl-的擴散系數由0.056 5 A2/ps增大到0.061 5 s、0.066 7 A2/ps。在溶液濃度為濃度2的情況下，對應三種溫度（298、308、318 K）下，Cl-在CTA膜內的擴散系數分別為0.027 4、0.036 4、0.045 0 A2/ps，不斷增大。濃度3時，Cl-的擴散系數在不同溫度下的變化趨勢與前兩種濃度一致。由此可見，在溶液濃度一定的情況下，溫度越高，Cl-在正滲透膜（CTA膜）內的擴散系數越大。在溫度為298 K的情況下，溶液濃度從濃度1增大為濃度2、濃度3時，擴散系數由0.056 5 A2/ps減小為0.027 4、0.009 1 A2/ps。溫度為308 K時，三種不同濃度對應的擴散系數分別為0.061 5、0.036 4 s、0.019 9 A2/ps。溫度為318 K時，三種不同濃度對應的擴散系數有著與其他溫度下同樣的變化趨勢。由此可見，在環境溫度條件一定時，溶液的濃度越高，Cl-在正滲透膜內的擴散系數越小。

3 結論

本文利用Materials Studio軟件，采用分子動力學的模擬方法，計算出了三種不同溫度和三種不同溶液濃度下水分子、鈉離子、氯離子在三醋酸纖維膜（CTA）內的擴散系數。結果表明：溶液濃度不變的情況下，隨著溫度的升高擴散系數不斷增大；溫度不變的情況下，隨著溶液濃度的增大擴散系數不斷減小。同時，此次有效的動力學模擬對后續的正滲透過程的模擬研究起到了指導作用，驗證了模擬研究的可行性與可信性。

［1］ 王毅,等. 三醋酸纖維(CTA)正滲透膜活性層朝向對正滲透過程的影響[J]. 四川環境, 2014, 33(2): 17-24.

［2］ 李春艷，劉華,劉波濤.分子模擬的方法及其應用. 當代化工，2011,40(5):494-497.

［3］金昊白. 接枝親和配基的聚乙烯醇凝膠與組織型纖溶酶原激活劑的相互作用[D].浙江大學,2011.

［4］Kotelyanskii, M.J., N.J. Wagner , M.E. Paulaitis, Atomistic simulation of water and salt transport in the reverse osmosis membrane FT-30[J]. Journal of Membrane Science, 1998,139(1): 1-16.

［5］Kotelyanskii, M.J., N.J. Wagner , M.E. Paulaitis, Molecular dynamics simulation study of the mechanisms of water diffusion in a hydrated, amorphous polyamide[J]. Computational & Theoretical Polymer Science, 1999. 9(3-4): 301-306.

［6］劉清芝, 楊登峰,胡仰棟. 水和鹽分子在反滲透膜內擴散過程的分子模擬[J]. 高等學校化學學報, 2009, 30(3): 568-572.

［7］陶長貴,等.氧氣在聚丙烯內吸附和擴散的分子模擬[J]. 物理化學學報, 2009,25(7): 1373-1378.

［8］ 梁太寧,楊小震. 分子力場——COMPASS簡介[J]. 化學通報:網絡版, 2001(1).

［9］ 岳雅娟,等. 有機分子在聚乙烯膜中擴散過程的分子動力學模擬[J]. 化工學報, 2012. 63(1): 109-113.

［10］景蓓蓓. 化學發光用于界面傳質研究初探[D].陜西師范大學,2011.

［11］趙青. 摻雜對鈦酸鋇基陶瓷材料的結構、性能及氧擴散的影響, 遼寧科技大學,2012.

［12］陳璐斌,等. 正向滲透除硼試驗研究[J]. 后勤工程學院學報, 2016,32(3): 62-67.

Molecular Dynamics Simulation in the Process of Forward Osmosis

1,1, 2,1, 2,1

（1. Department of Architecture Planning and Environment Engineering,Logistical Engineering University, Chongqing 401311, China； 2. Water Industry and Environment Engineering Technology Research Center, Chongqing 401311, China）

The diffusion coefficients of H2O,Na+,Cl-in the forward osmosis membrane were calculated by molecular dynamics simulation using Materials Studio 6.0.At the same time, the variation law of diffusion coefficients of three particles in the forward osmosis membrane under three different ambient temperature and solution concentration conditions was simulated. The simulation results were in agreement with the results of experiments: the diffusion coefficient increased with the increase of temperature when the concentration of the solution was constant; the diffusion coefficient decreased with the increase of the concentration of the solution when the ambient temperature was constant.

Forward osmosis membrane; Diffusion coefficients; Molecular dynamics; MS simulation

TQ 018

A

1671-0460（2017）10-2022-04

2017-07-12

黃政宇（1993-），男，江蘇丹陽人，碩士研究生，研究方向：水處理技術與裝備。E-mail：838967960@qq.com。

周寧玉，副教授。E-mail：415525869@qq.com。