纳米多孔石墨烯膜分离N2/CH4的密度泛函计算与分子动力学模拟

贾晓霞， 杨江峰， 王 勇， 李晋平

(1.太原理工大学 精细化工研究所， 山西 太原 030024； 2.气体能源高效清洁利用山西省重点实验室， 山西 太原 030024)

天然气，与其他化石燃料(如煤和石油等)相比，每能量单位所产生的温室气CO2较少，被认为是一种清洁能源[1-2]。而我国能源以煤炭为主体，如今面临着降低碳排放的巨大压力，倡导使用天然气这一低碳的气体能源将成为大势所趋，但我国的常规天然气储量较为匮乏，近年来的进口数量达到1011m3/a，因此进行非常规天然气的开发与利用是对我国能源短缺的有力补充[3]。在我国非常规天然气中，煤层气储量丰富，达到了36×1012m3 [4]，但我国煤层气的低渗透性使得地面抽采极其困难，而通过井下抽采的煤层气中会混入大量空气，其中N2体积分数达到40%～60%，N2的存在降低了煤层气的燃烧热值[5-6]，因此开发高效的CH4/N2分离技术是发展煤层气资源的关键步骤。CH4和N2是比较难分离的一对气体混合物：它们分子动力学相近(CH4和N2直径分别为 0.38nm和0.36nm)，N2略小，扩散能力高于CH4；而 CH4本身极化率高于N2，吸附剂对其的吸附势能和吸附热普遍高于N2，传统的分离技术对CH4/N2的分离效果普遍较差[7]，采用吸附分离虽有较好的分离效果，但却受吸附剂稳定性及高成本问题的局限不利于其工业化发展。

在气体分离领域，相对传统分离技术，膜分离由于成本低、易操作、占地面积小以及易维护等优点，是未来气体分离技术的发展方向之一[8-9]。其中，石墨烯(Graphene)是由sp2杂化的碳原子构成的单层蜂窝状二维网格结构物质，具有良好的机械性能和化学稳定性、原料来源丰富等优点。由于石墨烯仅具有单原子层的厚度，可克服膜分离材料在厚度方面对气体通量的限制，是用于气体分离的理想膜材料[10-13]。然而无缺陷石墨烯由于其致密的结构，对于气体分子具有不可渗透性，但可通过去除石墨烯上的部分碳原子形成纳米孔，基于分子筛分效应实现分离气体的目的[14-15]。目前已有大量实验和理论研究表明，纳米多孔石墨烯具体很好的气体渗透性和选择性[16-19]。Li等[20]制备了具有结构缺陷的氧化石墨烯膜，通过实验测得其对H2/CO2和H2/N2混合气体的选择性系数分别达到3400和900；Zhou等[21]通过密度泛函(DFT)计算表明H修饰的多孔石墨烯能够很好地从CO2、CO 和CH4中分离H2；Liu等[22]通过分子动力学(MD)模拟了N和H修饰的多孔石墨烯对CO2/N2的分离，结果表明，CO2穿过膜的速率远高于N2。

在多孔石墨烯中，由于部分碳原子的缺失导致悬空键的出现，其结构具有很强的化学活性，而H饱和可以起到稳定结构的作用。因此本研究选用了不同孔径大小的H饱和多孔石墨烯，通过密度泛函理论计算和分子动力学模拟相结合的方法，研究了CH4/N2混合气体的分离过程，以期为材料的实际应用以及新材料的设计提供参考。

1 计算方法与模型

图1为H饱和纳米孔石墨烯模型。如图1所示，通过去除部分C原子挖孔的方式建立了纳米多孔石墨烯模型，并用H原子对纳米孔进行饱和。考虑了3种不同孔径的纳米孔，依据去除C原子数量分别命名为H-pore-10、 H-pore-13和H-pore-16。

事实上，实验制备这一孔型的纳米多孔石墨烯存在一定的难度，但是基于目前的实验进展我们认为还是有望实现的。例如，Koenig等[23]利用紫外光刻蚀的方法在机械剥离的石墨烯上制备了直径为0.34 nm 的纳米孔。

由于石墨烯上部分C原子的去除，石墨烯的结构发生变化，因此在MD模拟之前，首先通过DFT计算对所设计的多孔石墨烯模型进行了结构优化。该结构优化计算是通过Materials Studio软件中的DMol3模块实现的。模型的设计如下：对石墨烯进行6×6超胞以模拟无限平面，并在Z轴方向加2 nm的真空层，以避免周期性镜像结构的影响。通过广义梯度近似(GGA)来描述体系的交换关联能，并利用PBE关联梯度修正泛函对模型进行充分优化[24-25]。此外，采用双数值极化基组(DNP)，并将原子轨道的截断半径设置为0.6 nm。能量收敛判据为10-5Ha，原子受力的收敛判据为0.02 Ha/nm，电子自洽计算的能量收敛判据为10-6Ha。在描述气体分子与石墨烯膜间的相互作用时，笔者使用了色散力校正DFT(DFT-D2)方法。

通过MD模拟，笔者研究了纳米多孔石墨烯的气体分离选择性和气体通量。MD模拟是通过Materials Studio软件中Forcite模块实现的。石墨烯膜中的原子以及各气体分子的力场参数从Materials Studio 中COMPASS力场选择，COMPASS力场是第一个基于从头计算、能够准确预测孤立态和凝聚态分子的分子结构、构象、振动、热力学性质的分子力场[26]。模拟模型的构建参考文献[8]中的方法，2个优化好的纳米多孔石墨烯片(3 nm×3 nm)将整个模拟盒子的高(12 nm)分割成3等份，构成3个相等的模拟盒子，石墨烯上的所有C原子保持固定，而孔边缘的H原子允许运动。在中间的盒子放入54个CH4和54个N2分子，将此模拟盒子命名为气室，两边的盒子命名为侧室(扩散区域)。在模拟过程中，如果气体分子出现在侧室，则此气体分子被认为是穿过气体。在3个维度上使用周期性边界，模拟采用NVT系综，通过Andersen恒温方法将模拟温度控制在298 K，模拟时间为20 ns，步长为1.0 fs，5 ps输出1帧。气体与气体以及气体与石墨烯之间的相互作用由Lennard-Jones (L-J) 势能模型和库仑作用力相结合的方式描述。所有 L-J 作用点的交互作用参数均使用Lorentz-Berthelot 混合规则计算，库伦相互作用采用Ewald算法；L-J 势能截断半径取1.2 nm。

2 分子动力学模拟及气体穿过石墨烯膜的能垒分析

首先，对多孔石墨烯进行了结构优化，并通过计算石墨烯电子密度等势面的方式，考察了石墨烯纳米孔的尺寸。图2为纳米多孔石墨烯膜的电子密度等势面图(20 e/nm3)。孔的尺寸可定义为电子等势面孔内边缘的最短距离(s)和最长距离(l)的平均值，如表1所示。

图2 纳米多孔石墨烯电子密度等势面Fig.2 The resulting pore electron density isosurfaces of the nanoporous graphene sheets Isovalue of 20 e/nm3(a) H-pore-10; (b) H-pore-13; (c) H-pore-16

Graphene membranes/nml/nmAverage pore diameter/nmH-pore-100.2560.3700.313H-pore-130.4060.4060.406H-pore-160.3860.5520.469

sandlrepresent the short and long lengths of the pores in the graphene sheets shown in Fig.2

通过MD模拟，笔者考察了H饱和纳米多孔石墨烯对CH4/N2的分离性能。初始状态下，气体分子聚集在模拟盒子的中间气室内，随着模拟时间的变化，气体可以扩散到两边的侧室，扩散到侧室的气体即为穿透分子。通过MD模拟计算得到的石墨烯膜H-pore-10、 H-pore-13和H-pore-16中穿透分子数目随模拟时间变化曲线如图3所示。图4为经过20 ns模拟后气体分子分布状况的快照。可以发现，在H-pore-10中，经过20 ns的MD模拟后，仍无任何气体穿过H-pore-10进入到侧室，说明H-pore-10对于CH4和N2不具有穿透性。这是由于CH4和N2的分子动力学直径(分别为0.38 nm和0.36 nm)均大于H-pore-10的孔径(0.313 nm)造成的。在H-pore-13中，可以观察到，随着模拟时间增加，N2分子不断穿过石墨烯膜进入到侧室，而未观察到任何CH4气体分子穿过，说明H-pore-13可以高选择性地分离N2/CH。决定分离膜性能的2个重要参数分别是分离选择性系数和气体通量，因此笔者对H-pore-13的N2通量也进行了考察。基于图3中N2随时间变化的穿透数曲线，根据文献[27]中的线性回归方法，可得到H-pore-13的N2摩尔流量为8.58×10-5mol/(s·m2·Pa)，即2.56×105GPU(1 GPU=3.35×10-10mol/(s·m2·Pa)，典型聚合物分离膜仅为～100 GPU)。综合上述分析，H-pore-13在天然气纯化方面具有广阔的应用前景。对于H-pore-16，由于其孔径尺寸(0.469 nm)远大于CH4和N2的分子动力学直径，CH4和N2均能较为容易地穿过该石墨烯膜，其对N2/CH4无分离性能。

图3 气体分子穿透数随时间变化的曲线Fig.3 The number profiles of permeated molecules through the graphene membranes as a function of simulation time (a) H-pore-10; (b) H-pore-13; (c) H-pore-16

图4 MD 模拟20ns后N2/CH4混合气体穿过石墨烯膜的快照Fig.4 Snapshots of the N2/CH4 mixture permeating through graphene after a MD simulations time of 20 ns C—Gray; N—Blue; H—White

为进一步探究多孔石墨烯分离气体混合物的机理，笔者通过DFT计算了气体与石墨烯间的相互作用能，得到相互作用能随气体分子在石墨烯上吸附高度变化的曲线，从而推算出气体穿过多孔石墨烯膜所需的能垒。相互作用能的计算方法：以石墨烯孔中心为原点，分子质心到原点的距离为X轴，每隔0.02 nm固定一次气体分子，并计算气体分子与石墨烯间的相互作用能Eint，Eint通过如下公式进行计算：

Eint=Esum-(Egraphene+Egas)

(1)

其中，Esum为石墨烯与气体分子的总能量，eV；Egraphene和Egas分别为孤立石墨烯和气体分子的能量，eV。气体在石墨烯孔中心的构型方式对能垒影响较大，笔者尝试了让气体以不同构型穿过石墨烯膜，发现当CH4分子中的3个H原子指向石墨烯平面，而另外1个H原子与C原子的连线垂直于石墨烯平面时，N2分子轴线垂直于石墨烯平面，气体分子穿过石墨烯膜所产生的能垒最小。图5为CH4和N2分子穿过石墨烯膜孔时气体分子与石墨烯间相互作用能曲线随气体分子到石墨烯吸附高度的变化。可以定义气体从最稳吸附位置到达石墨烯孔中心位置时最大作用能与最小作用能间的差值为能垒，得到的气体扩散能垒数值如表2所示。

图5 CH4和N2分子穿过石墨烯膜孔时气体与石墨烯膜间的相互作用能曲线随气体分子吸附高度的变化Fig.5 Interaction energy profiles for CH4 and N2 passing through the pores of graphene as a function of adsorption height(a) H-pore-10; (b) H-pore-13; (c) H-pore-16

表2 CH4和N2分子穿过石墨烯膜的能垒(Ebarrier)Table 2 Potential energy barriers (Ebarrier) for CH4 and N2passing through the graphene

相比于H-pore-13和H-pore-16，H-pore-10中CH4和N2穿过需要克服较大的能垒，分别为1.28 eV和0.56 eV，因此在MD模拟中CH4和N2均不易穿过H-pore-10。H-pore-13和H-pore-16中，由于孔径足够大，CH4和N2分子与石墨烯上孔边缘的原子的距离较远，其相互作用为吸引力。N2穿过H-pore-13所产生的能垒(0.05 eV)要远低于CH4穿过该膜所产生的能垒(0.13eV)，因此在MD模拟中H-pore-13对N2/CH4分离具有较高的选择性。对于H-pore-16，CH4和N2分子穿过膜所产生的能垒相等且较小，在MD模拟中CH4和N2均较易穿过该石墨烯膜，H-pore-16对N2/CH4不具有分离选择性。

为进一步理解气体穿过能垒的差异原因，笔者对石墨烯膜与气体分子间的相互作用力进行了分析。图6为石墨烯膜与气体分子相互作用力在极值点时的电子密度等势面。石墨烯膜与气体分子间的相互斥力与电子密度等势面的交叠有关[28]。可以发现，在H-pore-10中，由于其较小的孔径，气体分子与石墨烯的电子密度等势面存在很大交叠，其相互作用力表现为斥力，气体穿过需要克服很大的能垒。在H-pore-13中，结合图5的相互作用能曲线，可以发现，由于CH4与石墨烯膜间存在电子密度等势面交叠，而N2与该石墨烯膜无交叠，石墨烯膜对CH4的斥力要大于对N2的斥力，因此导致CH4穿过该石墨烯膜需要克服较大的能垒。在H-pore-16中，该石墨烯膜与气体分子间无电子密度等势面的交叠，其对气体分子的斥力较弱，这是CH4和N2穿过能垒几乎无差别的原因之一。石墨烯膜与气体分子间的非键相互作用力还包括静电力。为此，笔者分析了H-Pore-13中静电力对气体穿过能垒的影响。调整H-pore-13上的电荷，使其各个原子的电荷均为零，这样可排除静电力对膜分离效果的影响。对这一调整至电中性的石墨烯膜进行了气体分离MD模拟，得到随时间变化的穿透分子数曲线，如图7 所示。可以看出，N2分子的穿透数依然随时间的增加而不断增加，而仅有少许数目CH4分子穿过石墨烯膜。这说明静电力不是影响H-pore-13气体分离性能的主要因素。H-pore-13之所以能够选择性地分离N2/CH4，主要是由于N2/CH4与H-pore-13之间不同电子密度等势面交叠程度导致的气体穿过能垒差异。

图6 气体穿过石墨烯膜纳米孔时的电子密度等势面Fig.6 Electron-density isosurface for gases passing through the nanopores of the graphene membranesIsovalue of 50 e/nm3; (a) N2/H-pore-10; (b) CH4/H-pore-10; (c) N2/H-pore-13; (d) CH4/H-pore-13; (e) N2/H-pore-16; (f) CH4/H-pore-16

图7 在调整电荷至电中性的H-pore-13中气体混合物穿透数随时间的变化曲线Fig.7 The number profiles of permeated molecules for gas mixture through the neutralized H-pore-13 graphene membranes as a function of simulation time

3 结 论

通过DFT计算和MD模拟相结合的方法，考察了H饱和的多孔石墨烯对CH4/N2的分离性能。

(1) 由于H-pore-10孔径过小，气体分子与石墨烯的电子密度等势面存在很大交叠，CH4和N2穿过H-pore-10所产生的能垒较大，均不易穿过。

(2) H-pore-13中，由于CH4和N2不同的构型，CH4与石墨烯膜间存在电子密度等势面交叠，使其与石墨烯膜间的斥力较大，导致CH4和N2穿过H-pore-13所产生的能垒具有差异，进而使H-pore-13能够很好地分离N2/CH4。

(3) 由于H-pore-16孔径过大，CH4和N2与石墨烯间的相互作用力较弱，均能较易穿过，其对CH4/N2的分离不具有选择性。