正渗透过程中的分子动力学模拟研究

黄政宇，周宁玉,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。