《1. 引言》

1. 引言

随着世界人口的持续增长和工业的不断发展，淡水资源短缺已成为人类面临的一个重大问题[1]。与多效蒸馏、多级闪蒸和蒸汽压缩蒸馏等技术相比，反渗透（RO）技术被认为是缓解这一问题的有效工具之一，这主要因为其能耗较低[2]。然而，RO技术有其自身的瓶颈，包括水通量和离子截留之间此消彼长的关系。因此，研究人员努力设计出新的RO膜，以使其在不牺牲离子截留的情况下最大限度地提高水通量[3]。自20世纪70 年代研制出第一批商用复合RO膜以来[4]，研究人员通过优化分离层的结构和化学性质等方式，在提高其性能方面取得了重大进展[5]。然而，仅仅使用在过去几十年中获得的经验知识，似乎很难在不牺牲离子截留的情况下显著提高水通量[6]。为了解决这个问题，有必要在分子水平上加深对离子截留和水通过RO膜的传输机制的了解和认识。与实验方法相比，分子动力学（MD）模拟提供了一种观察纳米限制条件下粒子运动的强大方法[7]。

平板RO膜通常为薄膜复合结构，该机构由三层组成：超薄分离层、多孔夹层和无纺基质。超薄分离层决定了RO膜的渗透性和选择性，因此，大多数分子动力学模拟都集中于这一层。例如，通过MD模拟首次对RO 膜的物理特性进行研究[8–11]。研究人员还非常关注分离层的结构对水通量的影响，以及水穿过该层的分子机制。Ding等[12]对水通过RO膜的过程进行模拟，发现被限制在膜中的水分子形成了由环状和线性聚集体构成的相互连接的氢键网络。Gao等[13]证实RO膜的非平衡条件可以通过非平衡分子动力学（NEMD）模拟法进行模拟。Shen等[14]进一步观察到，随着连通渗透自由体积分数的增加，水的输送得到增强。Luo等[11]使用平衡MD模拟法预测来自自由能和离子路径的离子截留。在这些模拟中，通常使用带有外部水或溶液储层的水合膜。这种模型一般能够成功地模拟整个传递过程。此外，这些模拟的结果为改变单体的化学结构和聚合过程以获得更好的性能提供了有意义的指导。然而，这些模拟研究未区分RO过程总阻力中的界面阻力和内部阻力，从而导致实验水通量的估算和预测不可靠。

当水分子通过RO膜时，传递阻力由两部分产生，即水分子进出膜（界面阻力）和水在膜内传递（内部阻力） [15]。因此，这两部分的阻力应被分别对待，因为只有内部阻力与膜厚度有关，而且在实验条件（详见第2.1节）下估算水通量时，应只乘以这一部分阻力。为了分别获得界面阻力和内部阻力，需要建立两个模型，即有储水层的模型和无储水层的模型。无储水层的模型可模拟膜内水的流动行为，从而直接获得内部阻力。总阻力可以通过有储水层的模型获得，因为该模型可对整个过程进行模拟，包括水分子进入和离开RO膜以及水分子在RO膜内的传递。通过模拟分别获得总阻力和内部阻力之后，可以通过从总阻力中减去内部阻力来得到界面阻力。然后，用膜厚度系数乘以内部阻力（而不是总阻力），再加上界面阻力得到实验条件下膜的总阻力，即可估算出实验条件下膜厚度大于100 nm的纯水通量。

本文运用稳态非平衡分子动力学（SS-NEMD）模拟，通过直接在水分子的氧原子上施加外力（EF）来研究不同亲水性的膜的阻力分布。进行SS-NEMD模拟之后，计算出阻力，然后估算出实验条件下的水通量。通过将内部阻力和界面阻力分离，本文研究了亲水性和阻力之间的关系。此外，为了确认本文提出的方法的可靠性，将具有不同成分的聚酰胺（PA）RO膜的模拟水通量和离子截留与早期已发表研究[16]中的实验结果进行了比较。

《2. 方法》

2. 方法

《2.1. 聚酰胺 RO 膜的阻力组成》

2.1. 聚酰胺 RO 膜的阻力组成

通常将模拟结果与实验结果进行直接比较。然而，由于计算能力的限制，与实验中使用的实际厚度（大于 100 nm）相比，模拟采用的膜厚度相当小（5~10 nm）。因此，模拟得到的阻力通常是通过考虑厚度系数和水通量来估算的，如式（1）所示。

式中，J 为纯水通量；ΔP 为压降；µ 为水的动态黏度；R 为传递阻力；d 为膜厚度。

如果只采用一个模型（有储水层），则不可能从膜的总阻力（Rinterfacial + Rinterior）中区分出界面阻力和内部阻力。所以，在使用式（1b）进行估算的过程中会错误地放大界面阻力。因此，有必要区分界面阻力和内部阻力，具体如式（2）所示。

式中，Rinterfacial为界面阻力；Rinterior 为内部阻力。很明显，只有 Rinterior 与膜厚度成正比。

《2.2. RO 膜的构造》

2.2. RO 膜的构造

聚酰胺RO膜的组成对阻力的影响也值得关注，因为该影响规律将有助于设计出高性能聚酰胺RO膜。我们之前的研究[17]证明，聚酰胺分离层中残留的羧基对水分子有很强的亲和力，导致水的传递阻力增加。由于用典型的均苯三甲酰氯（TMC）和间苯二胺（MPD）界面聚合得到的聚酰胺层的化学性质不能被灵活的调节，因此，本研究选择了两种类型的酰卤单体来构建聚酰胺RO膜，即有四个酰卤基团的联苯四酰氯（BTEC）和有两个酰卤基团的间苯二甲酰氯（IPC）。通过调节 BTEC和IPC的比例，可以轻松地调节残余羧基的含量。由此可揭示膜的传递阻力，特别是分开内部阻力和界面阻力对聚酰胺RO膜组成的依赖性。如上所述，在两个模型（即一个有储水层的模型和另一个无储水层的模型）的模拟中分别采用具有一定亲水性的膜。

本研究采用大尺寸原子或分子大规模并行模拟器（LAMMPS）程序[18]进行SS-NEMD模拟。用适用于聚合物和无机材料的聚合物一致性力场（PCFF）模拟原子的相互作用[19]。系统的总势能（Etotal）可按式（3）进行定义。

式中，Evalence包括键伸展、角度、扭转和平面外能量；Ecross-term 包括键长度和角度变化；非键能（Enon-bond）分为L-J (9~6)范德华作用能和静电力作用能[20,21]。水模型也用PCFF进行表征。

继我们之前的研究[17]，我们期望在本研究中比较水分子进入（和离开）膜的阻力和水在膜内流动的阻力。因此，我们构建了有储水层和无储水层的膜，以便于分别计算总阻力和膜内的阻力。我们在模拟中使用的单体 3,3′,5,5′-BTEC、IPC和MPD与Zhao等[16]的实验研究使用的单体相同，如图1所示。首先用500个单位的MPD、 n 个单位的BTEC和m 个单位的IPC组成聚合物，并将其随机放置在5.5 nm×5.5 nm×5 nm的矩形盒中，以达到合适的密度[22]。聚合物的重复单元通过–COCl和–NH2 基团之间的相互作用形成的酰胺键而连接在一起。因此，我们基于启发式距离标准[23]在聚合物一端桥接– COCL，在另一端桥接–NH2。在系统中插入1600个水分子，以保持23%的含水量[10]。之后，进行能量最小化，以使分子均匀分布。此外，将未反应的–COCl基团人工转化为–COOH基团，从而形成无储水层的最终膜。

《图1》

图1. BTEC、IPC与MPD聚合生成聚酰胺的化学方程式。

使用上述相同的建模方法制造具有不同BTEC与 IPC比率的四种膜。为了确保各种膜具有相同的厚度和密度，我们采用尺寸相同的模拟盒。这四种膜的组成见表1，其中通过将IPC的质量除以酰卤的总质量计算出IPC的质量分数。为了得到具有储水层的膜，在水合膜的两侧再放置2000个水分子。在左侧，40个钠（Na）离子和40个氯（Cl）离子被水分子溶解，以模拟海水[24]。在所有三个方向上均采用周期性边界。在能量最小化之后，该系统在300 K和0.1 MPa下，在z方向上按等温-等压系综（NPT）运行1 ns，获得了正确的水密度。最终构型如图2所示。然后进行SS-NEMD模拟。

《表1》

表1 不同IPC含量的聚酰胺膜的组成

《图2》

图2. 水合聚酰胺（有储水层）最终构型的快照。用红色标记聚酰胺原子，用蓝色标记聚酰胺离子。

应该注意的是，这种方法只有在膜的结构是均匀的情况下才有效。在本模拟中，我们在所有三个方向上均采用周期性边界条件，并且我们建立的模型实际上是缩小版的实验性聚酰胺膜。因此，在预测膜的实验厚度时，模拟模型可以被视为均匀结构。

《2.3. 非平衡分子动力学模拟的细节》

2.3. 非平衡分子动力学模拟的细节

NEMD模拟是一种探索流体通过多孔介质传输的技术[25,26]。有两种方法可以实现压力梯度，一种是在颗粒上施加EF，另一种是使用活动壁（MW）作为活塞[27]。利用EF法，压降可以与作用于每个水分子上的恒定作用力关联起来。这种方法的一个优点是，研究人员可发现无储液层的膜内部的动力学特性，以防止界面效应对动力学行为产生影响[17,28]。Ritoes等[29]和Borg等 [30]提出了一种改进的EF法，他们通过仅对横截面均匀的区域施加高斯分布力，在两个储层之间施加压力差。这种方法非常适合复杂几何形状的模拟，特别是对于产生恒定压降的横截面积不均匀的情况。在MW法中，通过在上游壁上施加较大的力而在下游壁上施加较小的力来产生压降。事实上，这种方法也可以保持恒定的压力，但它打破了流动方向的周期性。

我们用分子模拟来计算总传递阻力，并将其分成两部分：界面阻力和内部阻力。因此，有必要建立一个无外部储层的分子模型。对于无储水层的模型，对所有水分子施加作用力是产生压降的唯一方法。此外，聚酰胺RO膜的水力渗透率极低，因此很难通过仅在聚合物膜外施加力来进行NEMD模拟。此外，Ding等[31]对该方法进行了验证，验证结果表明模拟值与实验值吻合较好。在本研究中，我们对流体原子（包括水分子和离子）施加EF，以模拟压降。为了防止出现虚假的旋转动力学行为，对水氧原子而不是水分子本身施加了EF [31]。在测量传递特性时，我们跟踪了模拟时间内通过膜的流体原子的位置。由于在模拟中采用的是三维周期条件，所以进水侧和渗透侧没有明确的边界。因此，离子截留是通过计算通过聚酰胺膜区域的离子数来确定的。此处模拟的压降比实验压降高两个数量级。这是因为较小的力可以较低信噪比。因此，在测量流体运动时，需要进行长时间的模拟来减少不确定性，这就使得计算成本增加。

利用正则系综（NVT）在300 K下调节聚酰胺膜的温度。如上所述，我们通过在水分子上施加力来获得恒定的能量输入。为了防止水温意外升高，根据不包括质心速度的温度，使水温在300 K下保持恒定，这在 NEMD模拟[32]中被广泛应用。Appendix A中的图S1提供了温度结果随带有储水层的IPC0模拟时间变化的示例。温度稳定的原因是通过这种方法可以很好地消除额外的能量。

在无储水层的情况下，我们对水分子施加0.1722、 0.1865、0.2000、0.2152和0.2296 kcal·mol–1 ·Å–1 （注： 1 kcal = 4184 J）的力（Fz）；在有储水层的情况下，我们对水分子施加0.0492、0.0533、0.0574、0.0615和 0.0656 kcal·mol–1 ·Å–1 的力，以沿z 方向产生不同的压降（600、650、770、800 MPa）。为了限制水流在流动水分子施加的力的作用下沿z 方向流动，并尽可能地保持聚酰胺的柔性，我们在模拟过程中固定了几个均匀分布在系统中的聚酰胺原子[14]。模拟进行了7 ns，其中前2 ns 的数据被视为达到稳态的数据，这些数据不会被用于进一步分析。按照每10 000步保存一次数据且时间步长为1 fs的方式采集轨迹。在整个模拟过程完成之后，我们使用这些轨迹分析了RO膜中水分子的结构和动力学特性。

《3. 结果和讨论》

3. 结果和讨论

《3.1. 膜的结构特性》

3.1. 膜的结构特性

由于结构性质必然影响RO膜中水分子的动态性质，因此有必要弄清楚RO膜的结构性质是否取决于单体组成。首先研究储水层和水合膜的密度，如图3所示。

《图3》

图3.对于具有不同IPC含量的膜，聚酰胺和水分子的密度随膜位置在z方向（膜厚度方向）上的变化而变化。（a）IPC0；（b）IPC20；（c）IPC40；（d） IPC60。“Dry polyamide”代表聚酰胺原子的密度（不包括其周围的水分子）。

显而易见，在所有情况下，水合膜都有三个区域：主体水区域、聚合物/水界面区域和内部区域。主体水区域的水密度保持在1 g·cm–3 左右，聚合物/水界面区域的水密度急剧下降，而内部区域的水密度最终达到约0.24 g·cm–3 。这些发现与已发表的模拟研究[17,33]的结果一致。聚合物/水界面的形成有助于聚酰胺的溶胀。在内部区域，密度在合理范围内的振荡受到膜内空隙的影响。从图3可以明显看出，在不同情况下，聚酰胺的密度分布并不完全相同。这种现象归因于聚酰胺膜中水簇的不同分布。尽管密度分布略有不同，但水合聚酰胺的平均密度几乎相同（约1.24 g·cm–3 ）。

应该注意的是，图3（a）的界面宽度明显大于其他三种膜的界面宽度，表明IPC0膜具有更高的溶胀度。这可能源于IPC0更亲水的特性。同时，较高的溶胀度可能会由于含水量较高而降低界面阻力。为了验证IPC0膜有较高的溶胀度，我们还考察了10 nm厚的IPC0膜。密度分布见Appendix A中的图S2，图中代表的溶胀度相似。

为了进一步证实这四种膜具有相似的内部孔径，我们计算了它们的孔径分布（PSD）。用Poreblazer v3.0测量孔径分布，该软件集成了其他地方描述的蒙特卡罗积分程序[34]。每种膜的孔径分布如图4所示。所有膜的孔径几乎都在4~10 Å范围内。由于孔隙是由聚酰胺形成，并且充满了水分子，因此最小孔径必须大于单个水分子的直径（约3.2 Å）。最大孔径小于10 Å，属于聚酰胺RO膜的平均孔径范围[35]。

每种膜的平均孔径的统计数据如图4所示。在图4 （a）和（c）中，孔径分布均匀，因为各孔径的比例基本相同。在图4（b）中，0.65 nm和0.95 nm的孔径占据了峰值截面。在图4（d）中，孔径主要分布在0.6~0.7 nm。尽管孔径的分布彼此不同，但按孔径分布的积分计算，每种膜的平均孔径约为0.72 nm。平均孔径与RO膜的孔径非常吻合——约为0.70 nm，这是其他研究[36]得到的估算值。由于水通量几乎与孔隙半径的平方呈线性关系 [37]，所以，每种膜的相似孔径确保了孔径因素对于具有各种IPC含量的膜几乎有相同的传递阻力。此外， 0.72 nm的平均孔径允许水分子渗入膜并阻碍离子通过 [38]。这四种膜的透水性和离子截留特性将在下文进行讨论。

《图4》

图4. 具有不同IPC含量的膜的孔径分布。（a）IPC0；（b）IPC20；（c）IPC40；（d）IPC60。

《3.2. 纯水通量》

3.2. 纯水通量

通过对膜施加五个压降（600、650、700、750和 800 MPa）以确定通过膜的水分子的动态特性。很明显，尽管由于随机的热运动而产生小的波动，水通量几乎与压降成正比。因此，较低压降下的水通量，如5.5 MPa 实验条件下的水通量可以被直接估算得到[14]。

对于单体组成的每种情况，我们采用了两种模型，即有储水层的模型和无储水层的模型。通过式（2a），由有储水层的模型可计算总阻力。内部阻力也可以由无储水层的模型计算。接下来，可以通过从总阻力中减去内部阻力来计算界面阻力。从图5（a）中可以看出，对于无储水层和有储水层的情况，斜率分别为370和140，然后可计算内部阻力和界面阻力，如图6所示。

《图5》

图5. 具有不同IPC含量的膜的水通量。（a）IPC0；（b）IPC20；（c）IPC40；（d）IPC60。

《图6》

图6. 不同IPC含量的膜（厚度为5 nm）的模拟阻力。

当IPC含量增加时，无储水层情况下的斜率上升到 491，而有储水层情况下的斜率下降到123。对于无储水层的情况，上升的斜率表明渗透率提高以及内部阻力降低。对于有储水层的情况，下降的斜率表明渗透率降低以及总阻力增加。

《3.3. 传递阻力分析》

3.3. 传递阻力分析

根据式（2）可知，利用两种不同的模型得到水流速度和压降之间的关系后，可以直接求得内部阻力和总阻力。界面阻力可通过从总阻力中减去内部阻力求得，因为在本模拟中，有储水层的膜和无储水层的膜具有相同的膜厚度。之后，模拟厚度条件下的内部阻力和界面阻力可以被求得。实验厚度条件下的内部阻力可以通过将厚度系数（即实际实验中使用的膜和模拟中使用的膜之间的厚度比）乘以模拟内部阻力求得。最后，通过确定实验厚度条件下界面阻力和等效内部阻力之和，得到实验厚度条件下的总阻力。值得一提的是，该过程是可靠的，因为假设膜厚度对界面阻力的影响可以忽略不计 [39]。图6给出了膜厚度为5 nm时的模拟阻力。

由于每个BTEC单元与MPD聚合后会留下一些残余羧基，而IPC单元则不会出现这种情况，所以IPC含量的增加将导致所产生的聚酰胺分离层中残余羧基的数量减少。我们之前的研究[17]证明，聚酰胺分离层中的羧基对水分子有很强的亲和力。因此，更少的羧基会吸引更少的水分子，因此水分子的内部阻力会减小。本研究得到的计算阻力再次证实了这一结论。如图6所示，随着IPC含量的增加（这意味着BTEC含量降低），内部阻力持续下降。

界面阻力的变化与内部阻力的变化相反。随着IPC 含量的增加，界面阻力明显增加，这与大多数分子动力学模拟研究结论[24]相一致。这种现象源于羧基和水分子之间的强大亲和力。残留在聚酰胺层表面的羧基越少，吸引的水分子数就越少，导致水分子对膜表面的润湿性越差，从而界面阻力也越大。

此外，图6揭示了对于不同的IPC含量，界面阻力达到总阻力的62%以上，表明决定内部阻力的膜厚度对 RO膜的总阻力的影响较小。这些模拟都是基于厚度为 5 nm的聚酰胺层，然而，真实RO膜的聚酰胺层的厚度通常约为200 nm。为了估算真实RO膜的阻力，一种典型的方法是根据不同厚度，将5 nm的聚酰胺层的阻力乘以系数40。就内部阻力而言，这种方法是合理的，因为内部阻力仅取决于膜内部的水传递的距离，即聚酰胺层的厚度。

界面阻力由表面性质决定，不随聚酰胺层的厚度而变化。因此，这种方法将大大高估界面阻力，从而高估总阻力。基于这种不可靠的方法可知，界面阻力在 RO膜的总阻力中也起着主要作用，因为在这种简单的估算中，内部阻力和界面阻力的比率并不随聚酰胺层的厚度而变化。这一结论显然与实验观察结果不一致，因为在实验观察中，聚酰胺层越薄造成的水渗透率越高 [40–42]。

通过构建带有储水层、厚度为10 nm和IPC0含量的膜进一步验证RO膜实验水通量的估算方法。同样在膜上施加五个压降（600、650、700、750和800 MPa）。水通量和压降之间的关系如图S3所示。斜率约为98，由此可计算出5.5 MPa时的水通量约为539 L·m–2 ·h–1 。由式（2）可知，基于5 nm情形的结果，可以求得厚度为10 nm时的水通量为559 L·m–2 ·h–1 。也就是说，误差仅为3.7%，这验证了第2.1节中提出的估算过程的正确性。

具有200 nm聚酰胺层的真实RO膜的阻力还可以通过分别考虑界面阻力和内部阻力的贡献率来进行估算。与5 nm聚酰胺层的界面阻力相比，其界面阻力保持不变，而将5 nm聚酰胺层的内部阻力乘以系数40，可以获得200 nm聚酰胺层的内部阻力。如图7所示，内部阻力在总阻力中起主要作用，并且对于所有四种膜，界面阻力的贡献率下降到10%以下。

《图7》

图7. 不同IPC含量的膜（厚度为200 nm）的外推阻力。

根据这一结论，RO膜的总阻力可以通过降低聚酰胺层的内部阻力来有效降低。为此，一种常用的方法是制备更薄的聚酰胺层——主要通过改变界面聚合的参数——这与实验观察结果一致。然而，由于界面聚合过程速度极快，聚酰胺层的厚度只能在相对有限的范围内减小。除了调整聚酰胺层的厚度之外，还可以通过改变聚酰胺层的化学组成来有效降低内部阻力。具体方法是，对由不同比例的混合单体产生的具有200 nm恒定厚度的聚酰胺层进行外推。由图7同样可知，总阻力随着 IPC含量的增加而下降，这表明对于由具有更高IPC含量的配方制备的膜，可以预期其水通量更大。令人印象深刻的是，这样的认识与参考文献[16]的实验结果非常吻合。

最近，Zhao等[16]进行了一项有趣的研究，通过在聚酰胺的界面聚合中使用两种酰卤单体混合物（BTEC 和IPC）与MPD反应，研究了聚酰胺层的化学组成对 RO性能的影响。他们发现，在不牺牲离子截留的情况下，水通量随着IPC含量的增加而稳定增加。例如，当 IPC含量从0%增加到60%时，水通量增加了34%。为了直接与这些实验结果进行比较，表2列出了5.5 MPa压降下200 nm厚的聚酰胺分离层的膜的估算水通量。很明显，水通量随着IPC含量的增加而增加。同时，没有观察到离子穿过聚酰胺分离层的膜，甚至没有离子能够渗透到聚酰胺区域。因此，所有情况下的离子截留都达到 100%。此外，当IPC含量从0%增加到60%时，水通量从48 L·m–2 ·h–1 增加至63 L·m–2 ·h–1 。也就是说，水通量增加了31.25%，这与实验结果[16]高度吻合——不仅在趋势方面，在纯水通量的增加方面也是如此。

《表2》

表2 在实验条件下（5.5 MPa及200 nm厚）估算的纯水通量，以及不同IPC含量的膜的离子截留

《4. 结论》

4. 结论

本文通过分别考虑界面阻力和内部阻力研究了水在聚酰胺RO膜中的传递阻力。按有储水层和无储水层两种条件，建立了模拟这两种阻力对传递阻力贡献率的模型。内部阻力取决于聚酰胺层的厚度和化学性质。当厚度很小时，如模拟时通常采用的5 nm，内部阻力起次要作用，而界面阻力占总阻力的62%以上。重要的是，我们发现对于较厚的聚酰胺层，这两个贡献率应该不会成比例地增加，因为界面阻力与聚酰胺层的厚度无关。由于真实的RO膜的聚酰胺层厚度通常为200 nm，因此界面阻力的贡献率显著降低到10%以下。为了研究聚酰胺层的化学性质对这两种阻力的影响，我们研究了由两种酰卤混合物与MPD通过界面聚合制备的RO膜。随着聚酰胺层中残留羧基数量的减少，内部阻力降低，而界面阻力增加。然而，由于内部阻力起主要作用，总阻力随着PA层中残留羧基数量的减少而降低。这些结论在趋势和通量增强方面均与实验结果高度吻合。这些模拟结果加深了我们对RO膜中水传递的认识，一方面可以可靠地预测水通量，另一方面可以为合理设计高通量RO膜提供指导。

《致谢》

致谢

感谢国家重点研究发展计划（2017YFC0403902）、国家基础研究计划（2015CB655301）、国家自然科学基金（21825803）、江苏省自然科学基金（BK20190085、 BK20150063）、江苏省高校优秀创新团队计划、江苏省高校学术优先发展计划（PAPD）的资助。同时还要感谢南京工业大学高性能计算中心和无锡国家超级计算中心为我们提供的计算资源支持。

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Yang Song, Mingjie Wei, Fang Xu, and Yong Wang declare that they have no conflict of interest or financial conflicts to disclose.

《Appendix A. Supplementary data》

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2020.03.008.