双层氧化石墨烯纳米体系中受限水的介电常数

胡 波，朱昊辰

（污染控制与资源再利用国家重点实验室,长江水环境教育部重点实验室,同济大学环境科学与工程学院,上海 200092）

近年来，随着世界人口的快速增长以及城市化、工业化的迅速进展，水资源所受到的影响十分巨大，污染问题层出不穷，人均淡水资源也因此短缺［1］.为了改善水资源循环和匮乏的现状，自20世纪下半叶，海水淡化已被广泛认为是实现可持续、安全的生活淡水供给的有效办法［2~4］.其中，反渗透技术已成为解决水资源短缺的主要方法之一.针对目前市场上的主流产品，聚合物膜存在对长期暴露于强氧化剂，特别是对氯的抵抗性较差的问题［5］，石墨烯基材料除了拥有前者高离子选择性的优点外［6］，还具有快速的水输送通道［7］、可控的层间间距［8］、高水渗透性［9，10］、良好的比表面积和较高的化学及机械稳定性等优势［11］，因此越来越多地被用于纳米多孔膜的制备和改性［12，13］.石墨烯的性能可以通过表面和边缘功能化以及制造缺陷、选择性允许水分子通过和排斥溶解盐或颗粒来调整碳主链的结构.氧化石墨烯（Graphene oxide，GO）在膜的基底和边缘存在含氧官能团，表现出许多类石墨烯的性质，这些含氧官能团的存在使氧化石墨烯膜具有更好的物理化学性质.富含氧的官能团为其后续的功能改性或与其它材料混合提供了便利.此外，氧化石墨烯上大量的孔洞和带电基团也为分离膜材料的制备提供了基础.但是，根据目前国内外的研究进展，多层石墨烯和氧化石墨烯层的实际运用相较于其潜在的应用范围而言并不理想.主要原因之一是还未完全了解在纳米尺度下涉及分离过程的物理行为，通过应用分子动力学模拟来揭示宏观实验无法解释的物理现象，为其更好地运用于分离领域提供了可行性.在所有的物理性质中，介电常数作为控制流体在受限系统中传质特性的参数起着重要作用.如，介电行为的变化直接影响着海水淡化应用中纳滤、反渗透和纳米流体的过程［14~16］.在纳米尺度下，特别是在氧化石墨烯层这样的受限体系中，水溶液的结构和动力学特性完全不同于处于开放状态下的水溶液，系统的规模也类似于分子间相互作用的尺度大小.

近年来，通过分子动力学（Molecular dynamics，MD）模拟研究受限水在不同基质表面附近的介电行为已有报道.如Renou等［17，18］通过一系列分子动力学模拟，研究了柱状二氧化硅纳米孔表面电荷密度对介电常数的影响，他们观察到，介电常数随表面电荷密度的增加而降低，这是由于强烈的电磁场限制了水分子的自由旋转.实际上，介电常数可以通过纳米通道内水分子偶极矩的波动来确定.受限条件下介电常数的降低，是由于水-基底界面附近水偶极子转动自由度的降低，削弱了水分子对外加电场的响应能力.计算介电常数时需要考虑一些因素.在线性、均匀、各向同性介质中，其是一个常数；但在像氧化石墨烯纳米通道的线性各向异性介质中，其是一个二阶张量；而在非均匀介质中，其又是介质内部位置的函数.为此，已经建立了多种理论构型来计算圆柱或平板限制几何中液体的局部介电常数.如Ghoufi等［16，19，20］已经证明了限制在二氧化硅纳米管中的水具有显著的介电各向异性.我们［21］对双层石墨烯膜的研究已证实这些发现，研究表明，受限水的轴向介电常数大于主体相中的介电常数，而径向介电常数却要小得多.也与Fumagalli等［22］通过原子力显微镜测量了两层石墨烯薄片与六方氮化硼之间纳米通道中纯水的介电常数的结论一致.

如前所述，水的介电性质在很大程度上取决于受限相中的结构和骨架的大小.对于GO双层膜，膜表面的氧化程度对受限水的介电性能有重要影响.本文通过分子动力学模拟研究了不同孔径和氧化程度的双层石墨烯中纯水的介电特性，并进一步探讨了石墨烯氧化程度和孔径大小对介电常数变化的影响.研究结果对于理解分子尺度上异常的基础物理行为以及对新型膜过滤和脱盐材料纳米孔的渗透改性具有重要意义和价值.

1 分子动力学模拟模型的构建

分子动力学模拟使用LAMMPS软件［23］完成.实际的石墨烯纳米孔道形状不规则，很难根据实际情况真实还原，可以将其等效成两层单层石墨烯膜平行叠加，层间距离为d，分别为0.6，0.9和1.2 nm；每个孔道两侧为两个水盒子，在x，y和z方向上的尺寸分别为3.6，3.6和3.6 nm（图1），每个水盒子由1560个水分子组成，密度为1000 kg/m3.不同孔径的石墨烯纳米通道有固定的长度6 nm，在两个水盒子和纳米通道之间的界面处放置了4片垂直于xy平面的石墨烯壁以限制水分子在纳米通道边界上的运动.

Fig.1 3D image of a nanochannel formed by two layers of graphene flanked by two storage pools

在氧化石墨烯纳米片上构建的通道是通过在石墨烯双分子层表面随机添加不同浓度的羟基实现的，羟基的浓度（Cr）即氧化程度的表达式为Cr=nOH/nC=0，10%，20%，30%，35%，35%，40%［其中，0代表纯石墨烯（PG），nOH和nC分别表示羟基和碳原子的数目］.利用基于点电荷的刚性非极化模型TIP4P/2005对水分子进行建模，其中H和O所带的电荷分别为0.5564 e和0 e［24］；利用包括范德华（vdW）和库伦静电相互作用的全原子势CHARMM27［25］力场描述原子之间的相互作用；利用Paticle-particle-particlemesh（PPPM）法［26］计算截断半径为1.2 nm的长程静电相互作用，均方根精度为10−4；利用12-6 Lennard-Jones（LJ）势函数模拟截断半径为1 nm的vdW相互作用力（V），定义如下：

式中：rij（nm）为两个原子之间的距离；εij（kJ/mol）表示相互作用强度的势阱的深度；σij（nm）为i和j之间的核心距离.根据文献［25，27，28］，将力场参数列于表1.

为了保证系统内水的密度为1000 kg/m3，首先在两侧的水盒子采用常温常压（NPT）系综分子动力学方法模拟300 K，0.1 MPa条件下水的密度（等温等压模拟采用Berendsen法［29］），以1 fs为时间步长使水分子自由运动600 ps.然后整个分子动力学模拟选取等体积等温（NVT）系综，在300 K温度下，在时间步长为1 fs的情况下，先进行了一次300 ps的短期运行，以便快速使系统接近平衡.当系统接近平衡时，为了确保平衡，同样以1 fs为时间步长进行5 ns的分子动力学平衡模拟.最后，再以1 fs为时间步长进行最后10 ns的模拟进程，并在此阶段采集数据.

Table 1 LJ potential and charge of the atoms used in the simulation*

所有分子动力学模拟采用周期性边界条件，并通过平衡分子动力学（EMD）进行，使用SHAKE算法［30］和Nosé-hoover恒温器［31，32］.

2 介电常数的计算

介电常数是反映极性溶剂微观形态的宏观概念，通常由长程关联所产生的电场和由氢键关联的极化之间的关系决定.根据Bonthuis等［33］提出的线性响应理论，在目前的分子动力学（MD）模拟中，水的介电常数是从总偶极矩（M）的波动中计算得出的.极化密度的变化［ΔP（r）］是由以下给出的波动-耗散介电函数的线性化形式得到：

式中：β为逆热能；F为一个外部均相电场；P（r）为瞬时极化密度；表示在不同构型上的平均值.

在准二维受限系统中，ε是一个张量，并且对受限体的形状和大小有着很强的依赖性.因此，水的介电特性在本质上是各向异性的，并表现出各种异常.在平板几何体系中，受限介质仅在z方向上不均匀（即无对称性），因此局部对角张量可以通过由与壁面平行（ε||）和与壁面正交（ε⊥）来定义.

根据Bonthuis的理论，在平板的体系中平行方向上介电常数的表达式如下：

式中：ε0为真空介电常数；ΔE||为电场的变化，与位移场的局部变化（D）相关.根据麦克斯韦方程ΔE||=E||，这意味着平均切向电场是恒定的.式（2）中的平行方向处的F||指的是E||.因此式（4）可以通过考虑式（2）来重新排列：

对于垂直方向的介电常数，其平均位移场为常数，ΔD⊥（z）=D⊥.因此，垂直处的介电常数可表示如下：

均相场F⊥对应于正交情况下的D⊥/ε0.因此，结合式（2）和式（6）得到：

值得注意的是，式（5）和式（7）是计算非均相系统在同一方向上介电张量的波动公式，且只适用于受限体系中介电常数的计算［33，34］.

因此，在给定的坐标系中，计算了两种不同的极化涨落所导致的不同方向处的介电常数.对于整体的平均静态介电常数（εavg）的计算如下：

3 结果与讨论

3.1 层间距及氧化程度对径向介电常数的影响

图2（A）~（C）比较了不同层间距的氧化石墨烯通道中，膜表面氧化程度对于受限水溶液径向介电常数［（ε（r）］的影响.需要注意的是，由于TIP4P/2005水模型的介电常数（εbw=60）比实验值（εbw=78）要小，故以两者比值的方式来呈现介电变化的规律.该方法可以更准确地比较受限系统中介电的行为.在受限状态下，特别是在最窄的孔道（d=0.6 nm）中，水的介电常数比块体（bulk）水中的介电常数低［图2（A）］.在0.9和1.2 nm孔道中，由于受限条件下水分子的分布导致靠近界面的局部介电常数增大［图2（B）和（C）］.此外，在纯石墨烯膜上添加羟基对水在受限体系中的介电性能有显著影响.为了更好地分析其中的原因，图2（D）展示了在不同纳米通道中水的径向平均介电常数与块体水的介电常数（εbw）的比值和氧化程度的函数关系.很明显，受限水的介电常数随纳米通道空间的减小而减小，因为在较窄的通道环境中，水分子会表现出更强的有序取向，减缓水分子的偶极波动.随着膜表面羟基浓度的增加，观察到较窄的层间距对介电常数的影响比较宽的层间距对介电常数的影响更大.当氧化程度达到40%时，0.6 nm孔道的介电常数比纯石墨烯的介电常数降低了77%，而1.2 nm孔道的介电常数降低了43%.随着孔道宽度的减小，氧化程度的增加使介电常数进一步降低.结果表明，当氧化程度达到最大值时，0.6 nm通道中水的介电常数相对于块体水的介电常数下降到10%.

Fig.2 Distribution of the ratio of the radial permittivity[ε(r)]of confined aqueous solution to that of water in the bulk system with the degree of oxidation in the double⁃layer graphene channels with d=0.6 nm(A),0.9 nm(B)and 1.2 nm(C),the ratio of their mean radial permittivity(εavg)to εbw with the degree of oxidation(D)

另一个有趣的现象是，对于最宽的孔道（d=1.2 nm），受限水的介电常数随膜表面氧化程度的增加而减小，而对于相对较窄的纳米孔道（即d=0.6和0.9 nm），介电行为会呈现出非单调的趋势.在0.6 nm孔道内，当Cr=20%时，介电常数出现了“跳跃点”；而在0.9 nm纳米孔道内，该“跳跃点”出现在了Cr=35%.为了进一步探究“跳跃点”可能会出现的孔径范围，考虑到水分子的尺寸，选取孔径为0.8 nm的孔道进行了相同的模拟实验，实验结果如图2（D）所示.在Cr=30%的情况下，该现象也出现在0.8 nm通道中.这表明不同氧化程度对介电行为的影响受纳米孔道间距的影响很大.纳米孔道间距越窄，介电常数对氧化程度的变化越敏感.同时，也可推断，在一定的孔径范围内，出现跳跃点所对应的氧化程度会随着孔道宽度的增加而增加.在0.6 nm通道中，随着羟基浓度的升高，介电常数开始下降，这是因为羟基的加入使通道内的空间迅速压缩，导致水分子的旋转受到强烈的限制.随着氧化程度的进一步增加，孔道壁面的亲水性增强，水分子的分布和取向变得更加不规则，从而导致水偶极子的波动增大.对于相对较宽的0.9 nm孔道，只有当孔道中膜表面氧化程度较高时“跳跃点”才会出现.这表明在较宽的PG孔道中，石墨烯表面OH基团浓度的增加（Cr≤30%时）所引起的孔道空间压缩对介电常数的降低起着主导作用.虽然随着氧化程度的增加，膜表面的亲水性增强，但在相对较低的氧化程度下，空间压缩效应是导致介电常数降低的主要原因.然而，当Cr=35%时情况出现了反转，与0.6 nm通道Cr=20%时发现的现象十分类似.当孔道宽度进一步加宽到1.2 nm时，随着氧化程度的增加，石墨烯表面并没有出现“跳跃点”.

以上结果表明，纯石墨烯通道中的介电行为不仅与氧化程度有关，而且与通道本身的宽度有关.当孔道较窄时，即使氧化程度低也会让介电性能呈现出非单调的变化趋势；相反，当孔道较宽时，这种变化也可以发生于氧化程度足够高的孔道中.为了验证这一猜想，在1.2 nm孔道中，将石墨烯表面的氧化程度一直提升到最高（80%）（即增加了Cr=50%，Cr=60%，Cr=65%，Cr=70%，Cr=80%）.结果发现，在Cr=70%时，介电常数发生了不均匀性，并同样出现了与0.6和0.9 nm的类似的“跳越点”.

3.2 受限水的局部密度与介电常数的相关性

图3（A）~（C）分别显示了在不同氧化程度下石墨烯通道内水密度的径向分布.可见，水分子的径向分布随通道间距的增大和石墨烯表面氧化程度的降低而变宽.此外，随着双层石墨烯膜壁上氧化程度的增加，在所有剖面上观察到的峰都变得更平缓，这可能是由于纳米通道内有效空间的减少.对于0.6 nm的通道［图3（A）］，观察到在如此狭窄的通道中只能容纳一层水分子，这主要是因为水分子的尺寸约为0.3 nm.这一发现与Algara-Siller等［35］的实验结果一致.图3（B）和（C）分别清楚地显示了0.9和1.2 nm通道中存在双层和三层水分子，与之前的同行实验和理论研究一致［27，35，36］.0.9和1.2 nm纳米通道之间的区别在于通道中心水分子的浓度（ρ）.前者表明双水层的水分子很少，而后者占据了石墨烯双层之间较多的水分子，模糊了两层之间的界面，意味着水分子的分布就如同在块体水中一样.

Fig.3 Radial profile of the water density as a function of various oxidation concentration in the double⁃layer graphene channels with d=0.6 nm(A),0.9 nm(B)and 1.2 nm(C)

通常，水分子的分布峰随氧化程度的增大而减小，而与纳米通道的大小无关，同时水分子浓度整体的分布形状随氧化程度的变化而保持不变.然而，进一步观察发现，对于0.6和0.9 nm的通道，水分子分布峰的强度并不随氧化浓度的增加而单调减弱.这与图2（D）中介电常数变化的情况完全一致（对于0.6和0.9 nm通道，峰值密度分别为ρ0.60%>ρ0.620%>ρ0.610%>ρ0.630%>ρ0.635%>ρ0.640%和ρ0.90%>ρ0.910%>ρ0.920%>ρ0.935%>ρ0.930%>ρ0.940%）.相反，在最宽的1.2 nm孔道中无此现象，这也与介电常数的情况相对应.综合考虑图2（D）和图3（A）~（C），可见，在0.6和0.9 nm等相对较窄的纳米通道中，通道尺寸引起的限制效应与膜表面亲水性引起的水分子随机分布是相互竞争的关系.通道越窄，因羟基加入而产生的亲水效应和受限效应之间的竞争就越早出现.这就是0.6 nm通道的介电常数在Cr=20%时发生异常变化，而0.9 nm通道内直到Cr=35%时才出现“跳跃点”的原因.然而，对于最宽的d=1.2 nm通道，由于空间足够大，即使在通道壁表面存在高浓度的羟基（即Cr=40%），也不会出现这种竞争效应.此时介电常数的不断降低主要是由于通道中水分子的旋转空间受到不断增多的羟基限制.

3.3 水分子与膜壁面夹角的概率分布

如前所述，介电常数的降低是由偶极迁移率下降时水分子的优先取向所引起的.为了解释这种物理现象，图4（A）~（D）给出了水分子偶极矢量（）与石墨烯膜表面所在的平面法线（）之间的夹角θ（图5）和其出现概率［P（θ）］的变化关系.为了更有效地比较孔道尺寸对介电常数的影响，图4（A）分别绘制了在不同孔径的纯石墨烯纳米孔道中水分子夹角θ的分布.结果表明，水分子在0.6 nm通道内具有高度的优先取向，峰的位置位于θ=90°处，说明在狭窄的通道中，水分子的构象倾向于与通道壁面垂直.随着孔道间距的增大，0.6 nm孔道处的尖峰消失，曲线变宽.水分子在0.9 nm通道中的分布峰位于θ=100°附近，水分子的平面也被认为与0.6 nm通道一样垂直于石墨烯壁面，略有约10°的波动.对于间隔更大的通道，如1.2 nm通道，与0.6 nm通道相比，水分子的偶极矩几乎失去了它的优势.θ角的最大概率分布在50°～78°之间，这意味着在较宽的通道中水分子有更多的构象选择.其实，纳米通道中自由空间的减少会导致水分子转动自由度的下降，进而抑制偶极波动，从而导致介电常数的降低.这些结果都完全验证了如图2（D）所示的不同纯石墨烯孔道介电常数的变化趋势，并与Giri等［37］的研究结论相符.

Fig.4 Probability of angular distribution of water molecules for PG nanochannels(A)and GO nanochannels with a size of 0.6 nm(B),0.9 nm(C)and 1.2 nm(D)

Fig.5 Schematic diagram of angle(θ)between the dipole vector of water molecules()and the normal of the graphene surface plane()

接下来，集中讨论氧化程度对不同大小通道中θ角分布的影响.在较窄的通道内（d=0.6，0.9 nm），水分子分布的峰值随氧化程度的增加而保持一致，而在大角度处曲线陡度变化较大［图4（B）和（C）］.并且对于这两个纳米通道，变化也并不是单调的.在0.6 nm通道内，当氧化程度增加到20%时，其曲线陡度与Cr=10%的曲线陡度发生了反转的情况；同样的现象发生在0.9 nm通道内，但是出现曲线陡度反转现象时膜表面的氧化程度分别为30%和35%.

这些发现与所研究的纳米通道内的介电行为和水分子的局部密度密切相关.对于间距较大（d=1.2 nm）的纳米通道，随着氧化程度的增加，曲线变得更加陡峭，表明水分子的取向受通道壁面上羟基的影响很大［图4（D）］.另一方面，注意到随着氧化程度的增加，角分布θ的峰值从50°～78°向100°～135°方向移动.但在较窄的通道中没有观察到这种峰值移动的现象，因为此时通道中已经没有足够的可用空间来允许水分子的其它构象存在.即在d=1.2 nm的纯石墨烯和氧化石墨烯纳米孔道中，尽管羟基的加入降低了介电常数，但其限制效应并未对介电行为产生严重的影响.这也是在这个相对较大的系统内，水分子取向的角分布是石墨烯氧化程度单调函数的原因.

3.4 氢键动力学分析

相比于开放环境，受限体系中介电常数的变化行为依赖于石墨烯通道的空间尺寸及氧化程度，表明在受限体系中，水分子的结构组织形式会有不同的变化.通常，介电性质与水分子的偶极波动有关，而水分子的偶极波动又受氢键网络的组成和强度的影响.因此，计算了PG/GO基质与水之间每个水分子的氢键数目以及它们相应的偶极矩关联时间（表2）.为此，应用了几何判据［38］，在r（O…H）<2.5 nm和150°<∠O…H—O<180°的条件下，考虑了两个水分子之间的氢键.通过计算氢键的平均寿命（τHB），即偶极弛豫时间，研究了氢键的动力学.τHB可以通过假设偶极矩的自相关函数的德拜指数行为来计算：C（t）=exp（−t/τHB），其中C（t）描述了偶极矩的自相关函数，并派生为C（t）=〈h（t0）h′（t+t0）〉/〈（h（t0））2〉，其中h为当两个被标记分子从t0时刻到t时刻形成氢键时的单位，否则为零.C（t）描述了氢键构型的弛豫时间，没有考虑氢键断裂时的动力学.由表2可见，氢键的数量和相应的偶极弛豫时间随孔道间距的增加而衰减，与氧化程度无关.对于最大的孔道（d=1.2 nm），PG纳米孔道的弛豫时间仅为0.72 ps，而对于相对较窄的孔道（d=0.6 nm），τHB=3.23 ps.表明水在更宽通道内比在小通道内通过得更快，说明纳米通道的渗透率很高.由表2也可见，水分子的偶极波动受到石墨烯表面OH基团的强烈限制，导致介电常数的降低.此外，随着nHB和τHB的减少，0.6 nm孔道中的氧化程度达到了20%.但是随着石墨烯表面OH基团浓度的进一步增加，此现象也不再出现.事实上，刚开始添加到PG表面的OH基团会导致水分子与膜表面迅速形成氢键，导致偶极迁移率和介电常数降低.在这个狭窄的纳米通道中，当氧化程度达到20%时，相邻的OH基团会在PG表面自动排列，形成氢键网络，从而削弱了其与通道中水分子结合形成水-基质氢键的能力.当Cr=40%时，在可用空间极大减少的情况下形成高数目的氢键，并在纳米孔道内形成刚性结构，造成了体积排斥效应［20］的加剧和介电常数的下降.当氧化程度增加到35%时，在更宽的0.9 nm通道中也出现了同样的现象.当氧化程度更高时，预计会出现刚性结构，水的双层结构也会开始消失［图3（B）］，诱导水分子向纳米通道的中心聚集.在1.2 nm处氢键数目最低，动态稳定性最好（对应的衰减速率最快），表明水分子在大的纳米通道中更不稳定，同时也比0.6和0.9 nm纳米通道中的有序性更差.这些结果强调了调节层间距对于理解水在双层GO膜所构成的纳米通道内的渗透及其传输机制的重要意义.

Table 2 Number of hydrogen bond per water molecule between GO substrate and water and their corresponding dipolar relaxation time(τHB)for both PG and GO nanochannels with various oxidation degree

4 结 论

对PG和GO纳米孔道中的水溶液进行了分子动力学模拟，研究了不同间隔的双层石墨烯膜的氧化程度对其介电性能的影响.受限水的介电常数随纳米通道空间的减小而减小，因为在较窄的通道环境中，水分子会表现出更强的有序取向.随着膜表面氧化程度的增加，较宽的层间距对介电常数的影响比较窄的纳米孔道的影响更大.当氧化程度达到40%时，0.6 nm孔道中受限水的介电常数比纯石墨烯中的介电常数降低了77%，而1.2 nm孔道中受限水的介电常数降低了43%.随着孔道宽度的变窄，GO膜表面氧化程度的增加使介电常数进一步降低.当氧化程度达到最大值时，0.6 nm孔道中水的介电常数相对于块体水中的介电常数下降到了10%.对受限水的局部密度和水分子取向的分析都解释了上述结果的原因和机理，并与所研究的纳米通道内水的介电行为密切相关.研究结果再次印证了改变层间距对调节水的渗透性及其传输机理的重要意义，并为以GO为基底的先进材料的制备及其在水处理中的应用提供了理论基础.