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

?

正滲透過程中的分子動(dòng)力學(xué)模擬研究

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

（1. 后勤工程學(xué)院 國(guó)防建筑規(guī)劃與環(huán)境工程系，重慶 401311； 2. 重慶市水工業(yè)與環(huán)境工程技術(shù)研究中心，重慶 401311）

采用分子動(dòng)力學(xué)模擬方法，利用Materials Studio 6.0軟件模擬計(jì)算了H2O、Na+、Cl-在正滲透膜內(nèi)的擴(kuò)散系數(shù)。同時(shí)模擬研究了三種不同環(huán)境溫度、溶液濃度下，三種粒子在正滲透膜內(nèi)擴(kuò)散系數(shù)的變化規(guī)律。實(shí)驗(yàn)結(jié)論與模擬結(jié)果一致，即溶液濃度不變的情況下，擴(kuò)散系數(shù)隨著溫度的升高而增大；環(huán)境溫度不變的情況下，擴(kuò)散系數(shù)隨著溶液濃度的增大而減小。

正滲透膜；擴(kuò)散系數(shù)；分子動(dòng)力學(xué)；MS模擬

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

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

1 模擬方法

1.1 晶格的構(gòu)建

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

1.2 分子力場(chǎng)

在整個(gè)模擬過程中，模擬皆是基于分子力學(xué)和分子力場(chǎng)的DISCOVER模塊實(shí)現(xiàn)的。所有的計(jì)算都是在COMPASS力場(chǎng)下進(jìn)行的。該力場(chǎng)把有機(jī)分子體系的力場(chǎng)與無機(jī)分子體系的力場(chǎng)相統(tǒng)一，能夠應(yīng)用于有機(jī)和無機(jī)小分子、高分子，一些金屬離子、金屬氧化物和金屬的模擬研究[8]。

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

靜電相互作用表示分子中各原子靜電荷的庫(kù)倫相互作用對(duì)勢(shì)能的貢獻(xiàn)，被描述成：

1.3 模擬設(shè)置

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

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

在模擬不同NaCl溶液濃度對(duì)水分子、鈉離子、氯離子在CTA膜內(nèi)的擴(kuò)散系數(shù)的影響時(shí)。濃度1代表晶格當(dāng)中包含285個(gè)水分子和1個(gè)NaCl分子，濃度2代表晶格當(dāng)中包含285個(gè)水分子和2個(gè)NaCl分子，濃度3代表晶格當(dāng)中包含285個(gè)水分子和3個(gè)NaCl分子

在濃度不變的情況下，改變?nèi)芤旱臏囟龋?98，308，318 K），研究不同NaCl溶液溫度對(duì)水分子、鈉離子、氯離子在CTA膜內(nèi)的擴(kuò)散系數(shù)的影響。

晶格構(gòu)建成功后，利用“Forcite”模塊，對(duì)于獲得的擴(kuò)散模擬模型，首先使用“Geometry Optimization”進(jìn)行結(jié)構(gòu)優(yōu)化，接著選擇“Anneal”進(jìn)行模擬退火，得到穩(wěn)定的晶格模型。之后，在“Dynamica”內(nèi)分別進(jìn)行20 ps的NVT及20 ps的NPT動(dòng)力學(xué)模擬，達(dá)到弛豫結(jié)構(gòu)的效果。在模擬過程的最后進(jìn)行40ps的NVT動(dòng)力學(xué)模擬，輸出內(nèi)容選擇“ALL”，不僅得到平衡后的晶格模型結(jié)構(gòu)，同時(shí)得到分子的運(yùn)動(dòng)軌跡，為接下來的分析計(jì)算做好鋪墊。

2 分析與討論

2.1 分析與討論

在單位時(shí)間內(nèi)通過垂直于擴(kuò)散方向的單位截面積的擴(kuò)散物質(zhì)流量，即擴(kuò)散通量J（Diffusion flux）與該截面處的濃度梯度成正比，其數(shù)學(xué)表達(dá)式[10]為：

式中，稱為擴(kuò)散系數(shù)，為擴(kuò)散物質(zhì)的體積濃度，d/d為濃度梯度，“—”號(hào)表示擴(kuò)散方向與濃度梯度方向相反。擴(kuò)散系數(shù)相當(dāng)于濃度梯度為1的擴(kuò)散通量，由此可見擴(kuò)散運(yùn)動(dòng)會(huì)隨著值的增大而變得快速。擴(kuò)散系數(shù)的計(jì)算如下：

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

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

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

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

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

當(dāng)發(fā)生愛因斯坦擴(kuò)散時(shí)，我們可以說得到的擴(kuò)散系數(shù)是有效可靠的。圖3是上述均方位移曲線圖的對(duì)數(shù)曲線圖。圖中，水分子的值一開始小于1，接著無限接近于1?？芍肿拥臄U(kuò)散發(fā)生了愛因斯坦擴(kuò)散。同理，鈉離子和氯離子的都發(fā)生了愛因斯坦擴(kuò)散，由此接下來計(jì)算得到的擴(kuò)散系數(shù)是可靠地。在本文的后續(xù)的模擬研究中，水分子、鈉離子、氯離子的擴(kuò)散都發(fā)生了愛因斯坦擴(kuò)散，計(jì)算出的擴(kuò)散系數(shù)都有效。

對(duì)得到的均方位移曲線，利用前面介紹的方法進(jìn)行計(jì)算，可以得到CTA膜內(nèi)鈉離子的擴(kuò)散系數(shù)為0.033 8 A2/ps，氯離子的擴(kuò)散系數(shù)為0.056 5 A2/ps，水分子的擴(kuò)散系數(shù)為0.188 8 A2/ps。

圖3 H2O、Na+、Cl-的MSD-t對(duì)數(shù)曲線

2.2 環(huán)境溫度、NaCl溶液濃度對(duì)擴(kuò)散系數(shù)的影響

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

（a）H2O

（b）Na+

（c）Cl-

圖4（a）為三種不同環(huán)境溫度、NaCl溶液濃度下，水分子在CTA膜內(nèi)的擴(kuò)散系數(shù)。溶液濃度為濃度1的情況下，溫度為298 K時(shí)水分子的擴(kuò)散系數(shù)為0.188 8 A2/ps，溫度由298 K升至308 K時(shí)，擴(kuò)散系數(shù)由0.188 8 A2/ps變?yōu)?.2054 A2/ps，溫度升至318 K時(shí)，擴(kuò)散系數(shù)擴(kuò)大至0.248 2 A2/ps。在溶液濃度為濃度2的情況下，對(duì)應(yīng)三種溫度（298、308、318 K）下，水分子在CTA膜內(nèi)的擴(kuò)散系數(shù)分別為0.176 1、0.186 2、0.239 4 A2/ps。濃度3時(shí)，對(duì)應(yīng)擴(kuò)散系數(shù)也有著同樣增大的趨勢(shì)。由此可見，在溶液濃度一定的情況下，水分子在正滲透膜（CTA膜）內(nèi)的擴(kuò)散系數(shù)隨溫度的升高而增大。在溫度為298 K的情況下，溶液濃度從濃度1升至濃度2、濃度3時(shí)，對(duì)應(yīng)的擴(kuò)散系數(shù)從0.1888 A2/ps減小至0.1761、0.1691 A2/ps。溫度為308 K時(shí)，三種不同濃度對(duì)應(yīng)的擴(kuò)散系數(shù)分別為0.205 4、0.186 2、0.185 7 A2/ps，逐漸減小。溫度為318 K時(shí)，擴(kuò)散系數(shù)的變化有著同樣的趨勢(shì)。由此可見，在環(huán)境溫度條件一定時(shí)，溶液的濃度越高，水分子在正滲透膜內(nèi)的擴(kuò)散系數(shù)越小。首先，這是與理論相符合的。其次，通過相關(guān)實(shí)驗(yàn)驗(yàn)證，當(dāng)?shù)蜐舛鹊柠}溶液作為原料液時(shí)，在宏觀方面也表現(xiàn)為，隨著溫度的增大或濃度的降低，水通量也相應(yīng)隨之增大[12]。

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

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

3 結(jié)論

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

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

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

［3］金昊白. 接枝親和配基的聚乙烯醇凝膠與組織型纖溶酶原激活劑的相互作用[D].浙江大學(xué),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］劉清芝, 楊登峰,胡仰棟. 水和鹽分子在反滲透膜內(nèi)擴(kuò)散過程的分子模擬[J]. 高等學(xué)?；瘜W(xué)學(xué)報(bào), 2009, 30(3): 568-572.

［7］陶長(zhǎng)貴,等.氧氣在聚丙烯內(nèi)吸附和擴(kuò)散的分子模擬[J]. 物理化學(xué)學(xué)報(bào), 2009,25(7): 1373-1378.

［8］ 梁太寧,楊小震. 分子力場(chǎng)——COMPASS簡(jiǎn)介[J]. 化學(xué)通報(bào):網(wǎng)絡(luò)版, 2001(1).

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

［10］景蓓蓓. 化學(xué)發(fā)光用于界面?zhèn)髻|(zhì)研究初探[D].陜西師范大學(xué),2011.

［11］趙青. 摻雜對(duì)鈦酸鋇基陶瓷材料的結(jié)構(gòu)、性能及氧擴(kuò)散的影響, 遼寧科技大學(xué),2012.

［12］陳璐斌,等. 正向滲透除硼試驗(yàn)研究[J]. 后勤工程學(xué)院學(xué)報(bào), 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-），男，江蘇丹陽人，碩士研究生，研究方向：水處理技術(shù)與裝備。E-mail：838967960@qq.com。

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