納米多孔石墨烯膜分離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的分離不具有選擇性。