趙 莉,方 彬,劉 冰,鄭海霞,展凱云
(1. 中國石油大學(華東) 理學院,山東 青島 266580;2. 中國地質大學(武漢) 工程學院,湖北 武漢 430074)
開設綜合性研究型實驗課程是一種新穎的教學模式,有助于激發學生的求知欲,培養創新能力、創新思維,擴寬知識面,提升綜合運用各種交叉學科知識的能力[1-3]。分子動力學模擬是近年來飛速發展的一種微觀尺度仿真方法,是研究凝聚態物理的有力工具,其在對水合物微觀乃至宏觀性質的研究中取得了重大突破,如對熱力學和力學性質、形成與分解等的研究方面。這是一門結合物理、數學和化學的綜合技術,它以經典力學、量子力學和統計力學為基礎,利用數值求解經典力學運動方程得到體系的運動軌跡,并給出相應的結構信息與性質。本文所述的實驗依托縱向科研項目,結合分子動力學的科研熱點和專業發展前沿,以水合物生長動力學為研究載體,探索了水合物晶體不同晶面對水合物生長動力學的影響,生動展示了固體力學的基礎理論,拓展了學生的科研視野。
水合物為客體小分子與水分子在低溫高壓條件下形成的籠型化合物,客體分子位于由水分子構成的氫鍵網絡(籠狀結構)中,最常見的客體分子為甲烷,即形成甲烷水合物[4]。自然界中永久凍土和深海沉積物中蘊藏大量甲烷水合物資源,其中海洋沉積物中占據90%,保守估計為傳統化石能源(煤、石油、天然氣)的2倍[5-6]。由此可見,水合物資源具有巨大的能源潛力,是未來可替代傳統能源的一種清潔能源。2017年,我國在南海海域水合物試采中取得了重大突破。目前,水合物已被我國列為第173個礦種,它的產業化是我國未來能源戰略的重要組成部分。除此之外,水合物特殊的結構特性使其在二氧化碳封存[7]、氣體混合物分離[8-9]、儲氫[10-11]等領域具有廣泛的應用。
現階段,由于開采技術的限制以及復雜的沉積環境等不利因素,人們對水合物的研究依舊停留在實驗室以及試開采階段。作為潛在的可替代能源,天然氣水合物資源的調查勘探以及未來的開發利用,都離不開對其形成條件和規律的研究。然而,水合物資源的開采研究是一項復雜的過程,涉及物理化學、工程地質學、環境學、機械制造等多種學科。本文主要針對的是水合物的生長研究,該過程對了解水合物在地層中的賦存模式具有至關重要的作用,是能源工業中天然氣成藏研究的重要課題。Larsen等人[12]通過單晶生長實驗觀察到,生長速度慢的晶面發育較好,而生長速度快的晶面通常在生長過程中會消失,不同晶面生長的差異性直接影響到水合物的分布模式。但是,到目前為止該過程的微觀機理尚不明確。
計算機、Material Studio軟件(MS)、Gromacs開源軟件包、VMD可視化軟件,以及Origin繪圖軟件。
(1)實驗開始之前,教師將固體力學相關書籍和學習資料、Gromacs和MS操作手冊、VMD可視化軟件以及Origin繪圖軟件等的使用說明發送給學生,引導學生自主學習。使學生充分了解固體力學的相關知識理論,掌握晶體類型以及晶面定義,掌握分子動力學模擬的理論基礎和Gromacs輸入文件的構建,能夠對可視化軟件及數據處理軟件進行靈活運用。
(2)教師提前準備關于水合物生長研究的相關文獻發給學生,讓學生提前閱讀,并在此基礎上,引導學生自主查閱資料,理清在文獻閱讀中所遇到的問題,找出文獻中的研究空白。
(3)學生通過閱讀以上資料,以及自主查閱文獻了解整個模擬流程和理論基礎,并嘗試進行實驗設計。
3.2.1 模型構建
模型構建是計算機模擬的關鍵過程,在本實驗設計中,準確構建水合物單晶胞是研究的基礎與關鍵。本文選取的是I型甲烷水合物,其模型搭建過程分兩步:第一步從X射線晶體學研究結果中獲取水合物晶體主體水分子中氧原子的初始位置;第二步在已經確定好位置的氧原子上插入相應的氫原子,并調整氫原子方向使其滿足 Bernal-Fowler法則,同時保證單晶胞偶極矩和勢能最小[13]。該晶胞包含8個甲烷分子與46個水分子(比例為1∶5.75)。在此單晶胞的基礎上,通過 MS可視化軟件切出我們所關注的三個晶體面(100)、(110)、(111),如圖1所示,紅色表示水分子,可通過氫鍵作用組成籠型結構,客體分子甲烷(灰色)則占據晶籠內部空間。相應的晶胞參數如表1所示。另外,對應單晶胞的三個晶面分別構建超晶胞:(100)晶面構建 3×3×2 超晶胞;(110)晶面構建 3×2×3超晶胞;(111)晶面構建1×2×4超晶胞。隨后在Z方向將超晶胞體系延伸10 nm,留下10 nm空白空間,并在該區域隨機放入甲烷與水分子,甲烷與水分子比例為1∶5.75,滿足I型甲烷水合物比例(1∶5.75),三個晶面體系初始構型的參數與甲烷水分子數如表2所示。

表1 不同晶面所對應的晶體參數以及切除時的表面向量

表2 三個晶面體系初始構型的參數與甲烷水分子數
3.2.2 模型優化及計算模擬
在利用分子動力學模擬之前,需要對初始構型進行優化,使其能量最小,保證其初始構型穩定,隨后進行100 ns的分子模擬用于數據采集與分析。采用帶隨機項的速度重新標度實現溫度耦合,并利用Parrinello-Rahman擴展-集成壓力控制對模擬盒子向量進行重新標度[14-15]。恒溫以及恒壓常數分別設置為0.2 ps和2.0 ps,并利用Particle-Mesh Ewald (PME)求和方法處理長程靜電相互作用[16]。原子間的范德瓦爾斯相互作用用 Lennard-Jones勢來描述,截斷半徑為12 ?,該半徑應小于仿真體系向量的一半。采用跳蛙算法用于牛頓運動方程的積分,時間步長設置為2 fs[17]。在X、Y、Z三個方向均考慮周期性邊界條件,水分子用 Tip4p/2005模型描述[18],甲烷分子采用OPLS-UA勢函數[19],模擬采用NPT系綜和周期性邊界條件。溫度壓力分別維持在260 K和50 MPa,因為在該溫壓條件下水合物生長驅動力較大,有利于生長并縮短模擬時間。
為了便于結果討論,我們將初始構型分為兩個部分,如圖2所示,左邊為水合物晶核,右邊為水氣混合液相。同時,為了更好地表征水合物的生長過程,右邊水氣混合相被均勻分割為10等份,并分別標記為L1—L10。

圖2 水合物生長動力學初始模型(以100表面體系為例)
為了確定在生長過程中哪些水分子和甲烷分子屬于水合物相,必須要有判斷標準。一般情況下,水氧原子可以在冰或水合物相中與相鄰水的4個氧原子形成四面體單元的頂點,F3參數基于104.25° O—O—O角排列,這與液相水氧原子中是非常不同的[20]。因此,在水合物生長模擬過程中,可以用Baez和Clancy提出的F3階參數表征水分子的局域狀態[21]。該算法可以提供固態水分子氫鍵網絡與標準四面體結構的偏離。

式中,θjik是3個水氧原子之間的角度,第i個水氧原子為中心,j和k水氧原子在i原子兩邊且位于i水氧原子第一水殼層以內(O—O原子徑向分布函數的第一極小值,半徑3.5 ?)。無序結構的F3值大于四面體結構的F3值,如籠形水合物。因此在我們的模擬中定義,當F3小于0.05時,水分子屬于水合物相。此外,當晶籠中多個宿主水分子脫離分子間氫鍵約束時,甲烷分子很容易從不完全水合物氫鍵籠中逃逸。因此,假設甲烷分子(5.5 ?)在第一水合物殼層內的水合物相(F3<0.05,定義如上)周圍的水分子數目大于15,則認為甲烷分子仍處于水合物的多面體晶籠中。
水合物生長過程是一個相變過程,水氣混合液相逐漸變為水氣包裹且具有固定對稱性的晶體結構,水合物生長過程是以水合物晶面為基底,以層方式逐步向液相生長的(見圖3)。水合物/液相界面先由水分子形成半晶籠,然后吸附客體甲烷分子用以完成生長。可以看出,在開始階段,由于甲烷分子隨機分布于水氣液相中,水合物生長較快,隨后由于甲烷分子在水中溶解度較低,而此時甲烷分子又過飽和,于是甲烷分子快速聚集成納米氣泡,受周期性邊界條件的影響,液相中甲烷氣呈現圓柱體狀,此時,水合物晶體生長驅動力較大,部分甲烷分子克服界面張力擴散至水合物晶體表面新形成的空腔中,該晶籠完成客體分子生長,但此過程需要時間較長,生長過程也較為緩慢。

圖4 模擬時間結束水合物生長體系快照(Y-Z方向)
從圖4可以看出,不同晶面水合物體系生長快慢不同。同為100 ns生長時間,100表面體系液相并未完全生長,而110及111表面體系生長均已完成。此外,100表面體系在模擬結束時,由于生成水合物量較少,大量甲烷氣體匯聚并單獨成為一相,整個體系有水合物相和液相,呈氣相分離狀態。110及 111兩個晶面體系生長完成,所有液相均已轉化為固相,但同時也有殘余的少量氣體分子聚集,說明生長過程中并不是所有的晶籠都被客體甲烷分子占據,存在部分未占據的水合物晶體。相比 110體系,111體系殘余氣體數更少。為更好地描述水合物生長的快慢,我們統計不同液相層L1—L5水合物相水氣分子數,如圖5所示。
圖5所示水合物生長是從L1—L5依次生長的。相比之下,100表面體系只完成到L2層的生長,110、111表面體系已經完成了L1—L5層的生長。另外,110與111的晶面體系雖然都完成了生長,但其完成時間卻有很大差異。110表面體系大概在 90 ns左右完成L5層的生長,而111表面體系在60 ns左右就完成L5層水合物生長。相比之下,111表面生長速率最大,110表面次之,而100表面最慢,這與實驗觀察所得結果基本一致,生長速率越慢,實際存在越多。

圖5 L1—L5層中水合物相水分子以及甲烷分子數隨時間的變化
為進一步闡明不同水合物晶面生長快慢的內在機制,統計了三種晶面在Z軸方向上的密度變化,如圖6所示。而表3中統計了不同晶面上水分子以及甲烷分子的數量。可以看出,不同晶面的密度以及分子數是不同的。就密度來說,111面密度最大,110面次之,而 100面最小;從分子數上說,111面所包含的甲烷數最多。因此,晶體生長快慢與表面密度以及甲烷分子數有很大關系,表面甲烷數越多,生長越快。晶面結構與主客分子數共同決定了生長速率。

圖6 不同晶面體系在Z軸方向上的密度變化

表3 不同晶面水合物晶體單位體積粒子數
本實驗屬于理論探究型實驗,以實際問題為背景進行理論研究,并構建了水合物晶面生長的動力學模型。本實驗采用的是分子尺度常用的仿真方法——分子動力學模擬。實驗中,通過問題分析構建不同分子尺度模型,結合統計力學及漲落耗散等理論綜合分析仿真結果,從納米級尺度上模擬水合物晶體的生長過程,并解釋了不同水合物晶面生長快慢的內在機制。學生在教師的指導下完成整個模擬實驗過程,不僅了解了科研選題,進一步掌握了物理化學、化學工程、能源技術和環境技術等理論知識及研究的基本方法,同時還可激勵學生的求知欲望和探索精神。此外,考慮課時量及學生情況的不同,還對實驗內容進行了深層次的拓展,使模擬內容更加接近實際應用。擴展內容如下:
(1)由于 I型水合物晶體對稱性為Pm3n,除文中提到的幾種常見晶面外,還有多種晶面,學生可自行學習并使用軟件構建初始構型,進行動力學模擬。
(2)由于真實環境下水合物地層溫度壓力會有所不同,學生可以研究不同溫壓條件下水合物的生長規律及不同晶面在該生長過程中的作用。
(3)在真實開采或流動保障過程中,會在流體中酌情添加促進劑或抑制劑,一般為有機或無機小分子(特殊情況下也有聚合物等),學生可自行查閱文獻,自主設計實驗,研究水合物在有添加劑情況下的生長情況及對晶面的影響。
(4)除傳質過程外,水合物生長過程為放熱過程,在考慮動力學過程的同時,學生也可通過更換系綜來研究傳熱快慢對生長的影響。此外,水合物在相變過程中液固變化導致的導熱系數改變,也可對水合物后續生長產生影響。
實驗課后,要求學生對實驗原理進行概括總結,對實驗數據進行整理分析,并按照科研論文的格式撰寫實驗報告。此外,學生還需闡述對該課題的理解,并提出實驗過程需改進的地方。對于優秀學生開展的拓展實驗內容部分,由教師協助對報告進行修改,并以科技論文形式向有關專業期刊投稿。
本綜合性創新實驗通過分子動力學模擬,研究了I型水合物晶體不同晶面對生長動力學的影響。研究發現,111面最快,110面次之,100面的生長最慢。該結論與實驗結果一致,且對研究甲烷水合物的勘探開發利用具有重要意義。
該實驗包含的知識體系復雜,涉及對物理化學、經典力學、分子化學、凝聚態物理、數學等多個交叉學科的綜合運用,有助于培養學生的科研創新能力、問題解決能力、團結協作能力,以及知識體系綜合運用能力。到目前為止,魏琦等6位本科生在專業期刊發表了相關的科研論文;劉曉旭、楊光瑞等4位學生通過自主查閱文獻擬定了相關的大學生創新科研項目并獲批。實踐表明,通過設計性實驗的綜合訓練,可有效鍛煉學生的科研探索精神,并大大縮短其融入科研團隊的時間。