石墨烯及二氧化钛孔道内小分子扩散的反应分子动力学模拟

，,2， ，

(1.浙江工业大学 化学工程学院，浙江 杭州 310014；2.浙江工业大学 理学院，浙江 杭州 310023)

纳米孔道广泛存在于生物、化学和化工材料等领域，例如离子通道、催化材料以及多孔材料吸附剂[1-2].流体分子在这些孔道中的行为由于限域作用表现出与宏观性质完全不同的特征，对物质的输运和结构变化等动力学性质有着重要影响[3-7]，因此分子在纳米孔道的扩散传递一直为人们所关注，也是现代化学工程的学科前沿.纳米孔道多种多样，但一维纳米管和二维纳米孔道是孔道组成的基本单元.目前，该领域主要采用实验观测和计算模拟来研究受限客体分子在孔道内的扩散传递行为，为薄膜分离、生物体蛋白通道和能源存储等过程提供机理解释[8-11].如通过实验研究小分子在氧化石墨烯表面扩散传递过程中发现该材料具有高效的筛分能力[12]，能快速分离尺寸不同的气体分子[13].然而在很多情况下，采用原位实验技术研究物系的传递过程相对较为困难，计算机模拟作为一种有效的研究手段在研究纳微尺度的反应、扩散和传递等领域展现了其强大的优势[14].如采用密度泛函计算研究氢分子在石墨烯表面的吸附现象[15]，揭示了氢分子在该体系中的稳定构型；采用粗粒化分子动力学方法研究水在碳纳米管非均相聚集团簇中的扩散过程[16]，为复杂三维结构体系的模拟提供思路.采用量化计算方法研究反应扩散过程对计算能力要求较高，基于经典力场的分子动力学模拟并不能将有机小分子在扩散过程中与孔道材料上可能的化学反应进行关联，而孔道上的化学反应会对传递过程产生明显影响.在Vanduin教授等提出ReaxFF反应力场后[17-18]，纳微尺度上的反应扩散模拟成为计算模拟领域的一个新亮点，如Kim等发展了适用于研究水在金红石(110)表面反应扩散的反应力场[19].目前，采用基于反应力场分子动力学方法研究纳米孔道中的有机小分子反应扩散传递性能的工作还很少，而这对于进一步研究工业应用中常见有机小分子在纳孔中的传递十分重要.

从目前研究情况来看，扩散传递模拟主要针对水、氧气和二氧化碳等无机小分子[6,20-21]，而针对工业常用的有机小分子(苯、苯酚、甲烷)的研究还较少.苯、苯酚和甲烷都是重要的化工原材料，苯通过Fe/MIF的微孔沸石能够一步氧化得到苯酚[22-23]，符合当今原子经济和绿色化工的理念；甲烷在纳米孔隙下的微观机理研究能够为开发页岩气提供重要的理论指导.为了解释传输孔道材料的表面性质对有机小分子的孔道扩散行为的影响，使用石墨烯和锐钛矿(101)孔道材料作为研究对象并将这2种体系进行对比，石墨烯是一种重要的二维碳材料，拥有优异的电学、光学和力学性质，但难以在常温下与小分子进行反应，而纳米二氧化钛有较强的催化活性，促使小分子在其表面进行化学反应.因此，笔者选用石墨烯或二氧化钛(101)受限空间为纳米孔道模型，利用分子动力学模拟对苯、苯酚和甲烷分子在其中的分布和反应扩散性质进行研究.

1 模拟方法

1.1 ReaxFF反应力场

ReaxFF反应力场以原子间的键级(Bond order, BO)关系为基础，通过建立体系能量与键级的关系而确定.该力场的系统能量被划分成多个部分[17]，包含了键能、配位校正能、三体键角应变能、四体扭转角应变能、库伦作用能和范德华力作用能等，具体表达形式为

Esystem=Ebond+Eover+Eunder+Eval+Epen+Etors+Econj+EvdWaals+Ecoulomb

(1)

键级BOij定义为原子间距离rij的函数，并且使用化学键中单键键级、双键键级和三键键级的定义将键级分为有化学意义的3个部分[17]，即

(2)

其中：r0为平衡键长；Pbo1，Pbo2，Pbo3，Pbo4分别为ReaxFF力场回归计算出的键级经验参数.

基于反应力场的反应动力学模拟方法是连接量子化学和分子动力学方法的桥梁，为大规模研究凝聚态性质及处理其中可能存在的化学反应过程提供了可能，并且基于反应力场的分子动力学方法已经得到了广泛应用.

1.2 模拟细节

分子动力学模拟过程是通过开源的LAMMPS软件实现，模拟中石墨烯的尺寸设为44.3 Å×42.6 Å，在周期性边界的立方模拟盒子中，2层石墨烯的间距设置分别为10，12，14，16，18，20 Å.锐钛矿(101)的尺寸为43.5 Å×45.3 Å，采用周期性边界的单斜模拟盒子(α=90°，β=90°，γ=110.30°)，并且固定所有钛原子，层间距的设置与石墨烯体系一致.层间分子(苯、苯酚、甲烷)的初始密度为0.2 g/cm3，模拟时间步长为0.25 fs，总时长500 ps，模拟中采用正则系综(NVT)系综，并使用Nose-Hoover控温算法恒温在300 K.其他相关参数都按照文献[24].

2 结果与讨论

2.1 苯分子在纳米孔道内的吸附与扩散

2.1.1 苯在受限空间内的分布情况

苯分子在石墨烯和二氧化钛纳米孔道中模拟500 ps后的分子分布情况如图1所示，在两种孔道材料体系中，苯分子均趋向于聚集在孔道材料附近.为了进一步考察分子分布，图2显示了沿z轴方向的原子数密度分布情况.以图2(a)体系为例，苯分子在孔道中呈对称分布，且苯分子在石墨烯近表面的浓度远高于远离石墨烯的区域，该高密度层被认为是吸附层[6].实际上在苯分子与石墨烯表层还存在一个过渡损耗层，其原因是苯分子作为非极性分子遇到疏水表面，石墨烯中的碳原子和苯分子距离过近时产生强烈的排斥作用，最终苯分子处在近石墨烯表面的平衡位置，运动能力明显减弱[6，25-26].在吸附层外，范德华力明显减弱，分子运动能力较强，分布比较均匀.进一步分析不同孔隙宽度对分布的影响，发现随着孔隙宽度增加，分子越趋向石墨烯表面聚集，而在远离石墨烯表面的区域分子分布较少.对比发现，苯分子在孔道宽度小于16 Å的石墨烯孔隙中，两侧的吸附层几乎可以充满孔道区域，而当宽度较大时，2个吸附层中间出现明显的低密度区域.在二氧化钛体系图2(b)中，二氧化钛为亲水表面，但由于苯分子比较稳定，不在锐钛矿二氧化钛(101)表面发生反应，仅构成物理吸附，并限制了苯分子在孔道的运动.此外，苯分子在二氧化钛表面不存在像石墨烯体系中的低密度过渡层，随着狭缝宽度增大，近表面区域的分子数密度增长速度放缓，虽然远离表面区域的分子数密度随着狭缝宽度的增加而减小，但该区域一直有苯分子存在.

图1 苯分子在受限空间内的分布示意图Fig.1 Snapshots of benzene molecules confined in graphene or TiO2 (101) nano-slits

图2 受限空间内垂直孔道方向苯分子的数密度分布Fig.2 Number densities distributions along z-direction of benzene molecules confined in nano-slits

2.1.2 动力学性质

分子在受限空间里的扩散传递过程与自由空间里明显不同，我们的关注点在于受限空间尺寸对扩散系数的影响.均方位移(Mean-square displacement，MSD)为系统内粒子移动距离平方的统计平均，通过Einstein关系可以计算出分子的扩散系数，这也是衡量分子扩散能力的重要依据.其表达式为

(3)

其中：t为选取计算扩散系数的时间段长度；N为体系内扩散分子的总原子数；rj(t)为原子j在t时刻的原子坐标；D为扩散系数；MSD为均方位移.

图3(a)为苯分子在不同厚度的石墨烯孔道中MSD随时间的演变过程，选取200～300 ps的部分估算斜率(与扩散系数成正比).结果显示：当石墨烯孔道厚度在10～14 Å时，苯分子的扩散系数随着孔道宽度的增大而增加；进一步扩大宽度时，扩散系数逐渐减小.孔道宽度在10 Å时，扩散系数最小，这是因为孔道宽度小，苯分子更频繁地与孔道表面碰撞，表现出扩散受阻.当孔道宽度达到14 Å时，扩散系数明显增大，模拟结果与文献[27]吻合.对于二氧化钛体系，从图3(b)看到在200～300 ps的变化不大，间距相对较小的体系扩散系数较小.这是受到苯分子自身作用半径的影响，被限制在平衡位置.

图3 不同受限尺寸下苯分子的均方位移Fig.3 MSD of benzene molecules confined in nano-slits

2.2 苯酚分子在纳米孔道内的吸附与扩散

2.2.1 苯酚在受限空间内的分布情况

同样沿z方向统计500 ps后体系中原子数密度的分布情况(图4).在石墨烯体系中，分子分布呈对称分布，且同样存在中间低密度区域，近石墨烯表面的苯酚分子运动能力较弱，在远离石墨烯表面的苯酚分子受到的范德华力减弱，分子运动能力较强，分布比较均匀.分析不同孔道宽度对分布的影响，发现随着孔隙宽度增加，分子越来越向近石墨烯表面聚集.在远离石墨烯表面的区域，分子数量越来越少，出现明显的两级分化.苯酚分子在孔道宽度达到14 Å的石墨烯夹缝中，只出现了两个波峰，而波峰之间的区域几乎没有苯酚分子.说明在狭缝宽度大于14 Å时，分子几乎全部聚集到了两个吸附层.

从图4(b)发现：每条曲线都有2个紧靠表面的波峰，但是不呈对称分布，波峰宽度也各不相同，这是由于苯酚分子较大，羟基H与桥位O结合使得分子构型不稳定，甚至与界面呈90°垂直，且有新产物生成，影响原子数密度的测量.2个波峰之间的位置始终有分子存在，且会随着层间距的增加数密度逐渐减小.比较图5(a，b)，发现苯酚分子在石墨烯体系的分布随缝隙宽度的增加，波峰越来越高，有更多分子聚集在吸附层，但在二氧化钛体系的分布与缝隙宽度关系比较模糊，主要是因为发生化学反应生成了新物质.

2.2.2 扩散行为分析

图5(a)为苯酚分子在不同厚度的石墨烯孔道中MSD随时间变化的情况，苯酚分子在孔道厚度10～14 Å时，扩散系数随着孔道宽度增加而增加，这是因为孔道宽度小，苯酚分子更频繁地与孔道表面碰撞，表现出扩散受阻.当孔道厚度为14 Å时扩散速率达到最大.而对于二氧化钛体系，苯酚分子几乎都吸附在孔道材料上，难以扩散，仅当孔道尺寸增加到20 Å时扩散系数才有略微增大的趋势.且在160 ps之前增长迅速，之后趋于0，可能的原因是孔道宽度的增大，远离表面的苯酚分子慢慢吸附到孔道表面附近.

从图4(b)的z轴分布上可以看到：在二氧化钛层附近有明显的苯酚分子聚集，此外根据MSD曲线，扩散相较于同尺寸石墨烯受限空间体系呈现明显的下降趋势，据此推测苯酚分子与二氧化钛层之间有较强相互作用，且在模拟条件下，苯酚分子的羟基H可以与二氧化钛的桥位O形成化学键，因此苯酚分子在二氧化钛体系中可能是由于化学反应的发生约束了分子在孔道中的扩散.

图5 不同受限尺寸下苯酚分子的均方位移Fig.5 MSD of phenol molecules confined in nano-slits

2.3 甲烷分子在纳米孔道内的吸附与扩散

2.3.1 甲烷在受限空间内的分布情况

从图6(a)的数密度分布可以看出：2个波峰成对称分布，波峰位置都是靠近孔道两侧，说明甲烷分子全部都聚集在2个石墨烯孔道表面附近，而z=0 Å附近的位置上几乎没有甲烷分子分布.波峰高度随着层间距的增加而逐渐增强，而波峰宽度基本相同，说明甲烷分子几乎全部聚集在2个吸附层.图6(b)为甲烷分子在二氧化钛孔道中的分布情况，看到每一条曲线都有2个波峰，波峰位置同样靠近孔道两侧，说明甲烷分子与表面有相对较强作用力.在远离孔道表面的区域，存在着少量甲烷分子；波峰强度随着孔道宽度的增大而增强；波峰宽度较大，说明甲烷分子与孔道材料相互作用力较弱.

图6 受限空间内垂直孔道方向甲烷分子的数密度分布Fig.6 Number densities distributions along z-direction of methane molecules confined in nano-slits

2.3.2 扩散行为分析

图7(a)为甲烷分子在不同间距的石墨烯孔道中的均方位移，可以发现均方位移与时间近似成正比，扩散行为由于受限作用，产生了尺寸现象，并在10 Å左右出现较大的扩散系数[28]，之后扩散随着孔道宽度的增大而逐渐减小.石墨烯和二氧化钛表面是2种表面性质截然不同的表面，比较甲烷分子石墨烯孔道和二氧化钛孔道图7(b)的模拟结果，发现甲烷分子在二氧化钛孔道的运动更受限，而在二氧化钛孔道内的扩散系数规律与分子在石墨烯孔道相反，随孔道的增大，扩散逐渐增强，在18 Å的时候达到最大.

图7 不同受限尺寸下甲烷分子的均方位移Fig.7 MSD of methane molecules confined in nano-slits

2.4 苯、苯酚和甲烷分子在纳米孔道中反应、扩散传递过程的比较

由于分子的物理结构及化学反应活性不同，苯、苯酚和甲烷在石墨烯孔道中的扩散行为也完全不同.比较图3(a)和图5(a)可以发现：苯分子在16 Å的石墨烯孔道中有相对最大的扩散系数；苯酚分子则在14 Å的石墨烯孔道中有最大的运动活性；而甲烷分子由于体积相对最小，且不会发生化学反应，在越小的石墨烯孔道中，有更大的扩散性质，在10 Å的孔道中扩散系数达到最大值4.383×10-8m2/s.

对于分子在锐钛矿二氧化钛孔道间的情况，由于苯酚分子的羟基H会与二氧化钛的表面桥位O成键，限制分子的运动，所以苯酚分子在孔道内的扩散系数远远小于苯分子和甲烷分子在孔道内的扩散系数.同样由于甲烷分子与孔道材料不发生化学反应，因此其扩散系数最大.

3 结 论

采用基于反应力场的分子动力学模拟方法，考察了3种有机小分子(苯、苯酚和甲烷)在间距10～20 Å的受限空间内的反应、扩散传递行为，并得到了不同孔道材料及不同受限空间尺寸下小分子流体的均方位移及密度分布.通过对模拟结果的分析，发现在相对惰性的石墨烯孔道内，苯和苯酚分子的扩散系数在孔道宽度范围内会有一个极大值，且都呈现出先增后减的趋势，说明存在最优扩散孔道尺寸；而甲烷分子的扩散系数随着孔道宽度的增加单调减小.在相对活性较大的锐钛矿二氧化钛孔道内，苯酚分子由于发生化学反应的原因扩散明显受限，而甲烷分子在孔道内的扩散则十分活跃.此外，受限空间内的有机小分子均趋向于在孔道材料两侧积聚.

参考文献：

[1] 朱英红,李艳芳,朱颖,等.MnO2/CNTs电催化剂及其对C—H键的选择性催化活化[J].浙江工业大学学报,2015,44(2):128-132.

[2] 黄海凤,戎文娟,陈虹宇,等.酮类有机废气在Y分子筛上吸附性能的研究[J].浙江工业大学学报,2014,42(5):513-518.

[3] HUMMER G, RASAIAH J C, NOWORYTA J P. Water conduction through the hydrophobic channel of a carbon nanotube[J]. Nature,2001,414(6860):188-190.

[4] SANSOM M S, BIGGIN P C. Biophysics: water at the nanoscale[J]. Nature,2001,414(6860):156-159.

[5] MASHL R J, JOSEPH S, ALURU N, et al. Anomalously immobilized water: a new water phase induced by confinement in nanotubes[J]. Nano letters,2003,3(5):589-592.

[6] CICERO G, GROSSMAN J C, SCHWEGLER E, et al. Water confined in nanotubes and between graphene sheets: a first principle study[J]. Journal of the american chemical society,2008,130(6):1871.

[7] WANG L, ZHAO J, FANG H. Water clusters confined in nonpolar cavities by ab initio calculations[J]. The journal of physical chemistry C,2008,112(31):11779-11785.

[8] YANG H Y, HAN Z J, YU S F, et al. Carbon nanotube membranes with ultrahigh specific adsorption capacity for water desalination and purification[J]. Nature communications,2013,4(2220):1-8.

[9] KIM J, LU W, QIU W, et al. Ultrafast hydration dynamics in the lipidic cubic phase: discrete water structures in nanochannels[J]. The journal of physical chemistry B,2006,110(43):21994-22000.

[10] ZANGI R. Water confined to a slab geometry: a review of recent computer simulation studies[J]. Journal of physics: condensed matter,2004,16(45):S5371.

[11] ZHENG J, LI L, TSAO H K, et al. Strong repulsive forces between protein and oligo (ethylene glycol) self-assembled monolayers: a molecular simulation study[J]. Biophysical journal,2005,89(1):158-166.

[12] JOSHI R, CARBONE P, WANG F C, et al. Precise and ultrafast molecular sieving through graphene oxide membranes[J]. Science,2014,343(6172):752-754.

[13] LI H, SONG Z, ZHANG X, et al. Ultrathin, molecular-sieving graphene oxide membranes for selective hydrogen separation[J]. Science,2013,342(6154):95-98.

[14] 何若愚,裘敏杰,陈骏,等.光谱和分子模拟法研究洛伐他汀与牛血清蛋白相互作用[J].浙江工业大学学报,2014,42(3):349-354.

[15] ARELLANO J, MOLINA L, RUBIO A, et al. Density functional study of adsorption of molecular hydrogen on graphene layers[J]. The journal of chemical physics,2000,112(18):8114-8119.

[17] VANDUIN A C, DASGUPTA S, LORANT F, et al. ReaxFF: a reactive force field for hydrocarbons[J]. The journal of physical chemistry A,2001,105(41):9396-9409.

[18] CHENOWETH K, VAN DUIN A C, GODDARD W A. ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation[J]. The journal of physical chemistry A,2008,112(5):1040-1053.

[19] KIM S Y, KUMAR N, PERSSON P, et al. Development of a ReaxFF reactive force field for titanium dioxide/water systems[J]. Langmuir,2013,29(25):7838-7846.

[20] LIU B, WU R, BAIMOVA J A, et al. Molecular dynamics study of pressure-driven water transport through graphene bilayers[J]. Physical chemistry chemical physics,2016,18(3):1886-1896.

[21] WEI M J, ZHOU J, LU X, et al. Diffusion of water molecules confined in slits of rutile TiO2(110) and graphite (0001)[J]. Fluid phase equilibria,2011,302(1):316-320.

[22] JIA J, PILLAI K S, SACHTLER W M. One-step oxidation of benzene to phenol with nitrous oxide over Fe/MFI catalysts[J]. Journal of catalysis,2004,221(1):119-126.

[23] HENSEN E, ZHU Q, JANSSEN R, et al. Selective oxidation of benzene to phenol with nitrous oxide over MFI zeolites: 1. on the role of iron and aluminum[J]. Journal of catalysis,2005,233(1):123-135.

[24] MARTI J, SALA J, GUARDIA E. Molecular dynamics simulations of water confined in graphene nanochannels: from ambient to supercritical environments[J]. Journal of molecular liquids,2010,153(1):72-78.

[25] ARGYRIS D, TUMMALA N R, STRIOLO A, et al. Molecular structure and dynamics in thin water films at the silica and graphite surfaces[J]. The journal of physical chemistry C,2008,112(35):13587-13599.

[26] MAMATKULOV S I, KHABIBULLAEV P K, NETZ R R. Water at hydrophobic substrates: curvature, pressure, and temperature effects[J]. Langmuir,2004,20(11):4756-4763.

[27] 叶宏飞.水在碳基纳米尺度通道内的粘性,扩散和剪切流动行为的模拟研究[D].大连：大连理工大学,2012.

[28] JIAO S P, XU Z P. Selective gas diffusion in graphene oxides membranes: a molecular dynamics simulations study[J]. ACS applied materials & interfaces,2015,7(17):9052-9059.