SiO2對甲烷水合物生成初期熱導率的影響

摘 要：針對目前多孔介質在甲烷水合物生成過程中微觀導熱機理研究的不足，基于SiO2組成的多孔介質縫隙模型，提出采用分子動力學方法進行模擬研究。根據模擬體系的微觀結構，采用平衡態分子動力學（equlibrium molecular dynamics，EMD）方法和非平衡態分子動力學（nonequilibrium molecular gynamics，NEMD）方法研究均相溶液在SiO2縫隙內的熱導率變化過程，分析SiO2完整縫隙和缺陷縫隙分子構象圖和密度分布云圖，獲得SiO2對甲烷水合物生成初期熱導率的影響規律，并進行微觀機理分析。結果表明，采用EMD方法得到含SiO2縫隙的均相溶液熱導率均值為0.53 W/（m·K），均相溶液的熱導率為0.74 W/（m·K）；采用NEMD方法得到均相溶液、完整SiO2縫隙、缺陷SiO2縫隙的熱導率分別為0.986、0.581和0.439 W/（m·K）。采用EMD方法和NEMD方法均驗證了SiO2縫隙的存在有利于甲烷水合物的生成和儲存，熱導率模擬值均接近實驗測量值，含有縫隙的SiO2可明顯降低模擬體系的熱導率，有利于甲烷水合物的生成和儲存。

關鍵詞：化學熱力學；甲烷水合物；密度分布；EMD；NEMD；SiO2縫隙；熱導率

中圖分類號：O642

文獻標識碼：A"" DOI：10.7535/hbkd.2025yx01007

收稿日期：2024-07-07；修回日期：2024-09-28；責任編輯：胡姝洋

基金項目：

國家自然科學基金（51906018）

第一作者簡介：

劉遠超（1977—），男，黑龍江五常人，副教授，博士，主要從事微納尺度傳熱、能源高效利用等方面的研究。

E-mail：liuyuanchao@bipt.edu.cn

Effect of SiO2 on thermal conductivity of methane

hydrate formation in early stage

LIU Yuanchao" ZHAO Tingting" LIU Xinhao" WANG Yuanning2

（1.School of Mechanical Engineering， Beijing Institute of Petrochemical Technology， Beijing 102617， China;

2.China National Petroleum Corporation Ningxia Petrochemical Branch， Yinchuan， Ningxia 750026， China）

Abstract：In view of the lack of research on the microscopic thermal conductivity mechanism in the formation of methane hydrate in porous media， a molecular dynamics method was proposed to simulate the porous media gap model composed of SiO2. According to the microstructure of the simulated system， EMD method and NEMD method were proposed to study the thermal conductivity change process of the homogeneous solution in the SiO2 gap， and the molecular conformation diagram and density distribution cloud diagram of SiO2 intact gap and defect gap were analyzed. The influence law of SiO2 on the thermal conductivity in the initial stage of methane hydrate formation was obtained and the microscopic mechanism was analyzed. The results show that the mean thermal conductivity of the homogeneous solution with SiO2 gap obtained by EMD method is 0.53 W/（m·K）， and that of the homogeneous solution is 0.74 W/（m·K）. The thermal conductivity of homogeneous solution， intact SiO2 gap and defective SiO2 gap obtained by NEMD method is 0.986， 0.581 and 0.439 W/（m·K）， respectively. The EMD method and NEMD method are used to verify that the existence of SiO2 gaps is conducive to the formation and storage of methane hydrate， and the simulated thermal conductivity values are close to the measured values. SiO2 containing gaps can significantly reduce the thermal conductivity of the simulated system， which is conducive to the formation and storage of methane hydrate.

Keywords：chemical thermodynamics; methane hydrate; density distribution;EMD; NEMD; SiO2 gap; thermal conductivity

甲烷水合物是一種具有特殊物理化學性質的清潔能源，研究甲烷水合物生成過程的熱力學特性具有重要意義，而探究海底沉積物對甲烷水合物生成初期熱導率的影響有助于揭示甲烷水合物生成過程的微觀機理。甲烷水合物中的甲烷濃度遠遠超過其在水中的溶解度，且在海底沉積物中能夠穩定存在。因此，以沉積物為介質，采用水合物形式儲存甲烷氣體前景可觀^[1]。

目前，多孔介質中甲烷水合物的生成過程在國內外已有研究。王山榕等分析了多孔介質的粒徑對甲烷水合物生成過程的影響，結果表明，多孔介質的粒徑越小，越有利于水合物生成。EM等^[3]研究了高嶺土納米管中水合物的生成，認為原始黏土納米管的孔隙尺寸和表面化學性質均符合氣體儲存材料的要求，可以作為吸附氣體的載體。方翔宇等^[4]發現粒徑為30 nm和50 nm的SiO2均對水合物生成及生長起促進作用，在粒徑相同時，SiO2質量分數越大，促進效果越明顯。萬麗華等^[5]研究發現，當溫度高于Debye溫度時，熱導率隨溫度升高變化較小。王思敏^[6]研究發現，水合物熱導率隨水合物籠子占有率下降而降低。ROSENBAUM等^[7]采用單側瞬態平面儀技術測定出壓縮甲烷水合物樣品在261.5～277.4 K時的熱導率為0.68 W/（m·K）。陳強^[8]采用粒徑為0.18～0.35 nm的天然海砂作為多孔介質，使用交叉熱線法測得含飽和水的沉積物體系熱導率為0.77 W/（m·K）。

盡管多孔介質中水合物的生成過程已有研究，然而分子層面的機理分析亟待補充。為了研究沉積物中水合物的生成過程及沉積物對熱導率的影響，本文擬采用分子動力學模擬方法，以含縫隙SiO2表示沉積物，分別采用平衡態分子動力學（equlibrium molecular dynamics，EMD）和

非平衡態分子動力學（nonequilibrium molecular gynamics，NEMD）方法模擬純水體系、完整SiO2縫隙和缺陷SiO2縫隙中不同溫度條件下均相溶液中水合物的生成情況，并對模擬體系的導熱過程進行微觀結構及機理分析。

1 模型建立

采用Materials Studio軟件搭建模擬體系，如圖1所示。模擬體系上下兩側是SiO2組成的多孔介質縫隙，中間部分是水分子和甲烷分子組成的均相溶液。

根據X射線晶體衍射結果，SI甲烷水合物晶胞包含 46個水分子、8個甲烷分子^[9]。本文在此晶胞分子數基礎上擴大16倍，構建由水分子和甲烷分子組成的均相溶液，溶液中包含736個水分子、128個甲烷分子。從Materials Studio晶體庫中導入SiO2晶體結構后，使用Buildgt; symmertygt;supercell功能進行擴胞，搭建SiO2納米簇，形成規則的SiO2縫隙。圖1 a）中SiO2層高為0.7 nm，圖1 b）是有缺陷的SiO2縫隙和均相溶液組成的模擬體系，與圖1 a）所示的體系作對比，研究不同結構的SiO2縫隙對甲烷水合物生成過程的熱力學特性影響。圖1 b）的SiO2層高為1.4 nm，缺陷SiO2縫隙缺口大小為0.267 2 nm×3 nm×1 nm。

采用Lammps軟件進行并行模擬計算，水分子采用TIP4P/Ew^[10]水模型，用SHAKE算法約束水分子的鍵長和鍵角，甲烷分子采用OPLS-AA^[11]力場，利用Lennard-Jones勢能函數描述SiO2分子間及原子間的相互作用。使用EMD方法計算溫度分別為250、255、265、270、275、285 K時，均相溶液和含SiO2縫隙均相溶液的熱導率值；采用NEMD方法計算冷源溫度為250 K、熱源溫度為275 K時均相溶液熱導率隨時間的變化規律，在30 ps的NVT弛豫之后進行20 ns的NPT模擬，時間步長為0.01 fs。

2 結果分析

2.1 密度云圖分布

SiO2納米簇可吸附體系中游離的甲烷分子，其表面有利于甲烷分子的聚集，為甲烷水合物的生成和氣體的儲存提供了有利條件^[12]。為了觀察SiO2縫隙對甲烷分子的吸附作用，圖2顯示了10 ns時SiO2縫隙內均相溶液的構象，SiO2縫隙表面處聚集了大量的甲烷分子，圖3為其密度分布云圖。

圖2中紅色代表氧原子，黃色代表硅原子，藍色代表碳原子，藍色虛線表示水分子之間的氫鍵。在模擬時間段內，體系中還未大量生成規則的甲烷水合物，而是在水分子較為密集的區域生成不穩定水合物團簇。圖2 a）中大量甲烷氣體聚集在SiO2縫隙的一側，水分子則涌向另一側，生成不規則的水合物團簇。圖2 b）中，甲烷分子分布較為均勻，不規則的SiO2表面增加了甲烷的吸附面積，為生成甲烷水合物提供更多的成核點^[13]。

結合分子構象圖和密度分布云圖可以發現，圖3 a）中間部分的分子密度小，甲烷分子和水分子分別分布在SiO2縫隙兩側，圖3 b）中水分子和甲烷分子分布較為均勻。這是由于SiO2縫隙存在缺陷，與均相溶液接觸面積增大，因此水分子和甲烷分子的運動范圍更廣，生成甲烷水合物的區域較為分散，有利于甲烷氣體儲存。綜合圖2 a）和圖3 a）可以觀察到水合物不規則團簇在固液界面處大量生成，甲烷氣體吸附在SiO2縫隙表面，運動幅度變小，與游離的水分子充分接觸，有利于甲烷水合物生成，滿足水合物界面成核假說^[14]。

2.2 熱導率計算

分子動力學計算熱導率的常用方法有EMD和NEMD。EMD方法的理論依據是線性響應理論，通過模擬系統的平衡態進而獲得熱導率，主要通過Green-Kubo方法實現。NEMD方法則是以施加溫差或施加熱流的方法對系統施加擾動，建立非穩態導熱過程，利用Fourier導熱定律計算熱導率^[15]。相比之下，EMD方法理論上采用周期性邊界條件可以實現計算無限大的體系，但是該方法得到收斂的熱流量需要較長的弛豫時間，計算量較大。NEMD則具有較好的收斂性，但早期的NEMD方法在能量守恒方面存在一些問題，經過修正得到了改善^[16-17]。

2.2.1 采用EMD方法計算熱導率

為驗證模型的可靠性，首先采用EMD方法，通過熱導率自相關函數是否收斂來判斷模型的合理性，計算常壓下均相溶液在溫度分別為250、255、265、270、275、285 K時的熱導率值。根據上節分析，溫度為275 K時，模擬體系主要為不規則的水合物，狀態為固態。固體中的能量傳輸依賴于晶格振動形式的原子運動。在非導體中，能量通過晶格波（聲子）傳遞；在導體中，能量通過晶格波和自由電子2種方式傳遞。

聲子是“晶格振動的簡正模能量量子”，熱傳導的過程是聲子間的能量交換過程。溫度高的地方，質點晶格振動較強，鄰近的質點溫度低，晶格振動弱，通過質點間的相互作用，振動較弱的質點在振動較強的質點影響下，振動加劇，這個過程熱量發生了轉移，從高溫傳到低溫，產生熱傳導。

EMD模擬是利用Green-Kubo方法^[18-19]將熱流自關聯函數對一定關聯時間進行積分以得到熱導率。在兩體勢作用下，

熱通量J和熱導率k計算見式（1）和（2）。

J=1V∑i

eivi+12∑ilt;j

［Fij（vi+vj）］rij ，（1）

kuv（t）=VkBT2∫t0〈Ju（0）Jv（t）〉dtn！r！（n-r）！ ，（2）

式中：ei為每個原子的能量；V為模擬體系的體積；vi、vj為速度矢量；Fij為粒子i受到粒子j的作用力；rij為粒子i、j之間的距離；kB為玻爾茲曼常數；T為溫度；〈Ju

（0）Jv（t）〉是熱流自相關函數（HCACF）在不同時間起點的平均。

圖4中均相溶液和完整SiO2縫隙內均相溶液的熱流自相關函數均收斂，證明本文模型和模擬過程的可靠性。均相溶液在16 ns時收斂，而SiO2縫隙內均相溶液收斂時間更短，在14 ns左右，說明SiO2縫隙的存在可以降低均相溶液的熱導率，有利于水合物的儲存。

采用EMD方法模擬了不同溫度下SiO2縫隙對均相溶液熱導率的影響，如圖5所示。從圖5中的數據可以得到，SiO2縫隙內均相溶液熱導率值為0.53 W/（m·K），均相溶液的熱導率為0.74 W/（m·K），模擬結果處于0.45～0.79 W/（m·K）之間，與ROSENBAUM等^［7］的實驗測量結果較為接近，且在相同溫度下，SiO2縫隙內的均相溶液熱導率值低于均相溶液熱導率值。與實驗測量的SI甲烷水合物熱導率值相比^［7］，模擬的SiO2縫隙內熱導率值明顯偏低，原因是SiO2縫隙增加了均相溶液的聲子散射，導致熱導率明顯降低。

SiO2體系熱導率降低，表明富含SiO2的沉積物均可作為良好的甲烷水合物儲存介質。圖5中均相溶液與SiO2縫隙內均相溶液熱導率的變化趨勢與溫度的相關性不明顯，這與實驗得到的結論一致。此結論表明2種體系在模擬溫度下主要以不穩定水合物團簇的形式存在，因此熱導率變化對溫度的依賴性較小。

2.2.2 采用NEMD方法計算熱導率

EMD方法適用于給定初始溫度、不受外界其他變量干擾的體系。為研究給定溫差對均相溶液熱導率的影響，本文采用NEMD方法，通過在模擬體系兩側（Z方向）外部施加溫差向體系施加擾動，達到穩態后統計結構中的熱流分布，進而計算出中間部分水合物的熱導率。非平衡態分子動力學熱導率的計算公式如式（3）所示，其中A為傳熱面積，ΔT為溫度梯度。

κ=-JA×ΔT 。（3）

熱傳導模型如圖6所示，冷、熱源設置在Z方向兩側，溫度分別為250和275 K，SiO2層為固定端，傳熱區為均相溶液。

為了計算模擬時間段內體系的溫度梯度，分別將純水體系、完整SiO2體系和缺陷SiO2體系沿Z方向分為30、52和66份，保持傳熱區長度接近，每100 000步更新1次，一共更新100次。

圖7為3種體系（黑色表示純水體系，藍色代表完整SiO2體系，紅色代表缺陷SiO2體系）在升溫過程的能量-時間變化的曲線，能夠反映出體系內水合物生成是否達到平衡。能量變化曲線波動小，幾乎接近直線，說明3種體系在模擬10 ns后接近穩定狀態。

能量曲線的斜率代表熱流，斜率越大，導熱系數越高。從圖7可以看出，體系兩側存在溫差時，單位時間體系內以導熱的方式通過中間部分傳遞的熱量從大到小依次排序為均相溶液、缺陷SiO2縫隙和完整SiO2縫隙。通過含縫隙SiO2的熱量較低，有利于甲烷水合物儲存。有缺陷的SiO2縫隙對熱量傳遞也存在一定的削弱作用。均相溶液不受SiO2縫隙的限制，熱量傳遞不受影響，與有SiO2縫隙的體系相比，水合物生成過程不穩定。

圖8為10 ns時完整SiO2縫隙傳熱區的溫度變化規律曲線，可以看出傳熱區域內的溫度隨位移變化曲線可近似擬合成直線，其斜率為該時刻的溫度梯度。

在求出各時刻的溫度梯度后，即可得到體系的熱導率隨時間變化的曲線，如圖9所示。圖9中，藍色代表均相溶液，黑色代表缺陷SiO2體系，紅色代表完整SiO2體系。從圖9可以看出，均相溶液中的傳熱區熱導率波動較大，熱導率平均值為0.986 W/（m·K）。完整SiO2體系傳熱區熱導率波動較為平緩，熱導率平均值為0.408 96 W/（m·K），與文獻[11]的實驗測量值接近。缺陷SiO2體系的傳熱區熱導率變化幅度處于兩者之間，熱導率平均值為0.896 W/（m·K）。

在波動趨于平穩后，3種體系最終的熱導率分別為0.986 W/（m·K）、0.581 W/（m·K）和0.439 W/（m·K）。

導致熱導率隨時間波動較大的主要原因包括：熱流從熱源處向冷源處傳遞，導致原子速度變大，從而引起體系內部溫度劇烈波動，體系處于非穩態，粒子運動較強；模擬過程中產生了聲子倒逆現象，SiO2縫隙的存在減緩了聲子倒逆現象，因此波動小于均相溶液。均相溶液處于非穩態時間較長，不利于水合物的生成和儲存。

與EMD方法相比，NEMD方法收斂速度更快。兩者在計算含SiO2縫隙均相溶液的熱導率時得到的結果相近，均接近實驗測量的海底沉積物縫隙中水合物的熱導率值，但能量波動較大，導致NEMD模擬測量的均相溶液熱導率高于實驗測量的SI型水合物的熱導率值（0.49～0.77 W/（m·K）），說明NEMD方法在計算均相溶液熱導率能量守恒方面仍存在問題。

3 結 語

采用分子動力學方法研究了海底沉積物對甲烷水合物生成初期熱導率的影響，模擬了完整SiO2縫隙和缺陷SiO2縫隙內均相溶液熱導率隨溫度和時間的變化，主要結論如下。

1）通過分析10 ns時SiO2完整縫隙和缺陷縫隙分子構象圖和密度分布云圖，發現缺陷SiO2表面增加了甲烷的吸附面積，為生成甲烷水合物提供更多的成核點，同時缺陷SiO2縫隙內甲烷水合物生成更均勻，儲氣率更高。

2）采用EMD方法得到了均相溶液和SiO2縫隙內均相溶液的熱導率。均相溶液和SiO2縫隙內均相溶液熱導率的變化趨勢與溫度的相關性不明顯。SiO2縫隙內均相溶液熱導率均值為0.53 W/（m·K），均相溶液的熱導率為0.74 W/（m·K），熱導率模擬值接近實驗測量值。

3）采用NEMD方法得到了均相溶液、完整SiO2縫隙和缺陷SiO2縫隙內均相溶液熱導率在固定溫差下隨時間的變化情況。穩定后，3種體系的熱導率分別為0.986、0.581和0.439 W/（m·K），SiO2縫隙內的均相溶液熱導率接近實驗測量值，而均相溶液處于非穩態時間較長，不利于甲烷水合物的生成和儲存。

4）采用EMD方法和NEMD方法均驗證了SiO2縫隙的存在有利于甲烷水合物的生成和儲存，缺陷SiO2縫隙內更利于甲烷水合物儲存。

本文模擬了甲烷水合物生成初期的密度分布和熱導率變化，模擬時間較短，體系內未生成大量的規則籠狀甲烷水合物，同時本文搭建的SiO2縫隙結構簡單，而實際多孔介質結構更為復雜。因此，在后續研究中，擬搭建更加完善的多孔介質模型，并適當增加模擬時間，研究不同類型規則的甲烷水合物的生成情況。

參考文獻/References：

[1] NAIR V C，GUPTA P，SANGWAIJ S.Natural gas production from a Marine clayey hydrate reservoir formed in seawater using depressurization at constant pressure，depressurization by constant rate gas release，thermal stimulation，and their implications for real field applications[J].Energy amp; Fuels，2019，33（4）：3108-3122.

王山榕，劉衛國，楊明軍，等.多孔介質內甲烷水合物生成動力學研究[J].熱科學與技術，2019，18（3）：173-178.

WANG Shanrong，LIU Weiguo，YANG Mingjun，et al.Kinetics of methane hydrate formation in porous media[J].Journal of Thermal Science and Technology，2019，18（3）：173-178.

[3] EM Y，STOPOREV A，SENENOV A，et al.Methane hydrate formation in halloysite clay nanotubes[J].ACS Sustainable Chemistry amp; Engineering，2020，8（21）：7860-7868.

[4] 方翔宇，寧伏龍，歐文佳，等.井筒條件下疏水納米SiO2對水合物形成的影響[J].中國科技論文，202 16（4）：370-376.

FANG Xiangyu，NING Fulong，OU Wenjia，et al.Effect of hydrophobic nano-SiO2 on hydrate formation in wellbore[J].China Sciencepaper，202 16（4）：370-376.

[5] 萬麗華，梁德青，吳能友，等.甲烷水合物導熱機理的分子動力學模擬[J].中國科學（化學），201 42（1）：52-59.

WAN Lihua，LIANG Deqing，WU Nengyou，et al.Molecular dynamics simulation on the mechanisms of thermal conduction in methane hydrates[J].Scientia Sinica Chimica，201 42（1）：52-59.

[6] 王思敏.基于分子動力學的天然氣水合物生長過程模擬及熱導率計算[D].哈爾濱：哈爾濱工業大學，2019.

WANG Simin.Simulation of the Growth Process and Calculation of Thermal Conductivity of Natural Gas Hydrate Based on Molecular Dynamics[D].Harbin：Harbin Institute of Technology，2019.

[7] ROSENBAUM E J，ENGLISH N J，JOHNSON J K，et al.Thermal conductivity of methane hydrate from experiment and molecular simulation[J].Journal of Physical Chemistry B，2007，111（46）：13194-13205.

[8] 陳強.多孔介質中天然氣水合物動態聚散過程的熱物性響應及熱分析應用研究[D].青島：中國海洋大學，2014.

CHEN Qiang.Thermal Properties Response and Thermal Analysis Technique Application of Natural Gas Hydrates During the Dynamic Formation and Dissociation Process in Porous Media[D].Qingdao：Ocean University of China，2014.

[9] 劉瞻.水合物成核過程及動力學抑制劑對其影響的模擬研究[D].北京：中國石油大學，2019.

LIU Zhan.Molecular Simulation of Hydrate Nucleation and the Influence of Kinetic Hydrate Inhibitors on It[D].Beijing：China University of Petroleum，2019.

[10]SIANI P，FRIGERIO G，DONADONI E，et al.Modeling zeta potential for nanoparticles in solution：Water flexibility matters[J].The Journal of Physical Chemistry C，Nanomaterials and Interfaces，202 127（19）：9236-9247.

[11]ROBERTSON M J，YUE Qian，ROBINSON M C，et al.Development and testing of the OPLS-AA/M force field for RNA[J].Journal of Chemical Theory and Computation，2019，15（4）：2734-2742.

[12]NGUYEN N N，NGUYENA V，STEELK M，et al.Interfacial gas enrichment at hydrophobic surfaces and the origin of promotion of gas hydrate formation by hydrophobic solid particles[J].Journal of Physical Chemistry C，2017，121（7）：3830-3840.

[13]李香璇，崔衛，馬挺，等.多孔通道內CO2水合物生成特性實驗研究[J].工程熱物理學報，202 45（6）：1773-1779.

LI Xiangxuan，CUI Wei，MA Ting，et al.Experimental study on the formation characteristics of CO2 hydrate in porous media[J].Journal of Engineering Thermophysics，202 45（6）：1773-1779.

[14]HAWTIN R W，QUIGLEY D，RODGER P M.Gas hydrate nucleation and cage formation at a water/methane interface[J].Physical Chemistry Chemical Physics，2008，10（32）：4853-4864.

[15]曹炳陽.一種模擬熱導率的非平衡分子動力學方法[J].計算物理，2007，24（4）：463-466.

CAO Bingyang.A nonequilibrium molecular dynamics method for thermal conductivity simulation[J].Chinese Journal of Computational Physics，2007，24（4）：463-466.

[16]FAYAZ-TORSHIZI M，XU Weilun，VELLA J R，et al.Use of boundary driven nonequilibrium molecular dynamics for determining transport diffusivities of multicomponent mixtures in nanoporous materials[J].The Journal of Physical Chemistry B，202 126（5）：1085-1100.

[17]IIDA S，TOMOSHI K.Free energy and kinetic rate calculation via non-equilibrium molecular simulation：Application to biomolecules[J].Biophysical Reviews，202 14（6）：1303-1314.

[18]CALDER N L，KARASIEV V V，TRICKEY S B.Kubo-greenwood electrical conductivity formulation and implementation for projector augmented wave datasets[J].Computer Physics Communications，2017，221：118-142.

[19]PANOUKIDOU M，WAND C R，CARBONE P.Comparison of equilibrium techniques for the viscosity calculation from DPD simulations[J].Soft Matter，202 17（36）：8343-8353.