不同分子模型對(duì)甲烷水合物分解微觀特性表征

李佳，梁貞菊，王照亮，趙健，唐大偉

（1 中國(guó)石油大學(xué)（華東）新能源學(xué)院，山東青島266580； 2 大連理工大學(xué)能源與動(dòng)力學(xué)院，遼寧大連116024）

引 言

甲烷水合物是在低溫高壓下形成的一種類(lèi)似于冰的晶體包合物[1-2]，它以固態(tài)形式賦存于深海沉積物和永久凍土地區(qū)[3-4]。水分子通過(guò)氫鍵形成晶格主體，甲烷分子被包裹在水分子形成的籠形結(jié)構(gòu)中，甲烷分子與水分子之間通過(guò)范德華力作用結(jié)合[5-6]。據(jù)估計(jì)，全球儲(chǔ)存在海底和永久凍土水合物中的甲烷約為1×1015～1.2×1017m3，全球天然氣水合物的潛在產(chǎn)量非?？捎^。由于其特殊結(jié)構(gòu)，天然氣水合物對(duì)溫度壓力變化非常敏感，極易分解[7]。目前，國(guó)內(nèi)外尚未實(shí)現(xiàn)甲烷水合物的大規(guī)模開(kāi)發(fā)，主要原因是還未實(shí)現(xiàn)對(duì)甲烷水合物的可控分解。因此，研究水合物的分解動(dòng)力學(xué)特性具有重要意義，可為水合物的開(kāi)發(fā)利用提供重要依據(jù)。

近年來(lái)，國(guó)內(nèi)外對(duì)水合物的分解過(guò)程進(jìn)行了大量的實(shí)驗(yàn)和理論研究[8-13]。研究發(fā)現(xiàn)，水合物的分解反應(yīng)涉及氣、液、固三相變化，是一個(gè)復(fù)雜的吸熱過(guò)程。在實(shí)驗(yàn)上很難從宏觀角度測(cè)定甲烷水合物分解的反應(yīng)機(jī)理，而分子動(dòng)力學(xué)模擬能夠克服實(shí)驗(yàn)的局限性，可為水合物分解機(jī)理的微觀研究提供強(qiáng)有力的支持。國(guó)內(nèi)外學(xué)者已采用多種主、客體分子模型進(jìn)行了水合物分解的分子動(dòng)力學(xué)模擬研究[14-15]。常見(jiàn)的應(yīng)用于甲烷水合物分解的客體分子模型主要有OPLS-UA[16]聯(lián)合原子模型及TKM-AA[17-18]全原子模型。常見(jiàn)的水分子模型主要有TIP3P、SPC/E、TIP4P、TIP4P/Ice、TIP4P/2005、TIP4P-ew 和TIP5P 模型。由于TIP3P 模型很難正確反映水合物的結(jié)構(gòu)特性，而TIP5P 模型存在計(jì)算量過(guò)大等缺點(diǎn)，目前在水合物的模擬中很少采用TIP3P 和TIP5P 模型[17]，廣泛應(yīng)用的主體分子模型主要為SPC/E 模型及四點(diǎn)水分子模型。然而，由于力場(chǎng)參數(shù)不同，不同水模型預(yù)測(cè)的水的物理性質(zhì)也會(huì)存在很大差異。因此，了解哪些水模型適用于在給定的熱力學(xué)條件下測(cè)量某一特性非常重要。對(duì)于某一特定的水模型，其力場(chǎng)參數(shù)不會(huì)隨著熱力學(xué)條件的變化而變化。因而單一的水模型無(wú)法同時(shí)適用于描述水合物生成(低溫)以及分解(高溫)條件下的兩種不同的熱力學(xué)狀態(tài)[19]。因此，需要兩種不同的水模型分別模擬水合物的生成和分解兩個(gè)獨(dú)立的過(guò)程。根據(jù)文獻(xiàn)可知，TIP4P/Ice[20-21]、TIP4P/2005[21-23]和TIP4P-Ew[24]水模型是天然氣水合物生成研究的首選模型。而對(duì)于水合物的分解特性研究多采用TIP4P[19,25]和SPC/E[16,18,26-28]模型研究水合物的分解動(dòng)力學(xué)。即便如此，不同模型預(yù)測(cè)的水合物的屬性依然存在很大差異。因此，尋找適合于天然氣水合物分解研究的主客體分子模型，通過(guò)使用不同的水模型對(duì)水合物分解進(jìn)行比較研究，并確定水合物分解的最佳水模型，仍然是一個(gè)需要關(guān)注的模擬研究領(lǐng)域。

綜上所述，本文分別采用水分子-甲烷的TIP4P-UA、TIP4P-AA、SPC/E-AA 和SPC/E-UA 模型模擬甲烷水合物的分解動(dòng)力學(xué)特性。首先，模擬水合物分解過(guò)程中平衡溫度、序參數(shù)、不同區(qū)域勢(shì)能、質(zhì)量密度以及逸出水合物的甲烷分子數(shù)目的變化等分解動(dòng)力學(xué)參數(shù)，然后對(duì)水、界面、水合物區(qū)域進(jìn)行了劃分，計(jì)算出分解過(guò)程的活化能，進(jìn)而對(duì)比不同水分子、甲烷分子模型對(duì)分解的影響并研究不同模型下分解的微觀動(dòng)力學(xué)機(jī)制。

1 計(jì)算方法

1.1 模擬體系

首先，構(gòu)建SI 型甲烷水合物晶胞，其空間群為Pm3n，晶格長(zhǎng)度為11.877 ?(1 ?=0.1 nm)，根據(jù)X 射線衍射實(shí)驗(yàn)的結(jié)果確定C、O 原子的空間位置[29]，利用Bernal-Fowler 法則為C、O 原子添加H 原子以滿足整體凈偶極矩為零[30]。然后，將甲烷水合物晶胞沿X、Y、Z 方向擴(kuò)展為3×3×6 的SI 型甲烷水合物超晶胞結(jié)構(gòu)，沿Z 軸在甲烷水合物的兩側(cè)各放置1513個(gè)液態(tài)水分子。構(gòu)建的XZ 平面內(nèi)水-甲烷水合物-水初始模擬結(jié)構(gòu)如圖1 所示，X、Y、Z 方向均采用周期性邊界條件, ABC 區(qū)域?qū)⒃诤笪?.6 節(jié)詳細(xì)介紹。

圖1 初始模型Fig.1 Initial configuration of simulation system

1.2 模型選取

對(duì)水分子，選擇SPC/E、TIP4P 兩種水分子模型進(jìn)行計(jì)算。其中，SPC/E 是三點(diǎn)模型，三個(gè)作用位置分別對(duì)應(yīng)水分子的三個(gè)原子，三點(diǎn)模型具有計(jì)算簡(jiǎn)便高效的特點(diǎn)。TIP4P 為四點(diǎn)模型，該模型假定有一個(gè)虛擬原子位于H-O-H 角平分線上，距離氧原子0.125 ?。四點(diǎn)模型更為復(fù)雜，可以更合理地描述液態(tài)水的結(jié)構(gòu)和熱力學(xué)性質(zhì)[31]。由于SPC/E、TIP4P模型均為剛性模型，模擬中需固定水分子的鍵長(zhǎng)和鍵角。甲烷分子作為客體分子，最為簡(jiǎn)單的模型為單一作用點(diǎn)的OPLS-UA聯(lián)合原子模型。此外，對(duì)甲烷分子的模擬還選用了TKM-AA 全原子模型，該模型有五個(gè)相互作用點(diǎn)，分別對(duì)應(yīng)甲烷分子的五個(gè)原子。各種模型的參數(shù)如表1所示。

表1 水和甲烷模型參數(shù)Table 1 Parameters of water and methane models

分子間相互作用勢(shì)包括Lennard-Jones(L-J)勢(shì)及靜電勢(shì)，如式(1)所示

不同原子間的L-J 勢(shì)參數(shù)由Lorentz-Berthelot混合規(guī)則確定，如式(2)所示

1.3 模擬方法

采用LAMMPS 軟件進(jìn)行模擬。采用速度-Verlet 算法積分牛頓運(yùn)動(dòng)方程，設(shè)置時(shí)間步長(zhǎng)為1 fs，分子間相互作用的截?cái)喟霃皆O(shè)置為12 ?，采用PPPM(particle-particle particle-mesh)法處理長(zhǎng)程庫(kù)侖力，計(jì)算精度為10-6。系統(tǒng)的初始溫度設(shè)置為300 K，壓力設(shè)置為10 MPa，模擬步驟如下。

(1)對(duì)水合物做固定處理，將液態(tài)水置于NVT 系綜，采用Nose-Hoover 熱浴使之達(dá)到目標(biāo)溫度，溫度的阻尼系數(shù)設(shè)置為0.2 ps，運(yùn)行200 ps；

(2)繼續(xù)固定水合物，將液態(tài)水轉(zhuǎn)入NPT 系綜，設(shè)置目標(biāo)壓力為10 MPa，繼續(xù)弛豫使之達(dá)到目標(biāo)溫度，溫度的阻尼系數(shù)設(shè)置為1 ps，壓力的阻尼系數(shù)設(shè)置為2.5 ps，運(yùn)行200 ps；

(3)解除對(duì)水合物的固定，將整個(gè)系統(tǒng)置于NPT系綜進(jìn)行控溫控壓處理，運(yùn)行30 ps；

(4)將整個(gè)系統(tǒng)置于NVE系綜分解，運(yùn)行2 ns。

2 結(jié)果與討論

2.1 甲烷水合物的平衡溫度及分解熱

隨著分解的進(jìn)行，水合物不斷從液態(tài)水吸收熱量，液態(tài)水的動(dòng)能逐漸降低，溫度隨之下降。圖2為分解過(guò)程中溫度隨時(shí)間的變化。由圖2 可見(jiàn)，分解初始階段可用于提供熱量的動(dòng)能比較充足，液態(tài)水溫度下降較快，隨著可用動(dòng)能的消耗，液態(tài)水可提供的熱量有限，溫度緩慢下降。類(lèi)似的溫度變化趨勢(shì)在Alavi 等[18]的模擬中也得到了驗(yàn)證。此外，在恒定體積條件下，甲烷分子從籠形結(jié)構(gòu)中逸出導(dǎo)致體系壓力升高，這也可能導(dǎo)致水合物進(jìn)入相穩(wěn)定區(qū)，從而使得分解過(guò)程終止。達(dá)到穩(wěn)定狀態(tài)時(shí)，液態(tài)水的溫度稱(chēng)為平衡溫度。此處采用式(3)計(jì)算體系的平衡溫度[25]。

針對(duì)圖2 中的溫度分布(黑色曲線)，每20 個(gè)數(shù)據(jù)點(diǎn)取其平均值，作出溫度隨時(shí)間的波動(dòng)曲線(紅色曲線)，然后通過(guò)式（3）進(jìn)行擬合，可以得到體系的平衡溫度(藍(lán)色虛線)。由圖2 可見(jiàn)，在初始溫度為300 K、初 始 壓 力 為10 MPa 的 條 件 下，TIP4P-AA 和TIP4P-UA 的平衡溫度較高，分別為(292.63±3.71)K和(291.86±1.55)K。SPC/E-AA 和SPC/E-UA 預(yù)測(cè)的平衡溫度分別為(286.36±4.33)K 和(282.89±3.62)K。由于設(shè)置的模擬初始溫度較高，所以在相同的壓力下，平衡溫度越高，水合物越容易達(dá)到平衡狀態(tài)不再分解，因而平衡溫度越高，結(jié)構(gòu)越穩(wěn)定[32]。Goel[33]測(cè)得10 MPa 甲烷水合物的平衡溫度為286.11 K，可見(jiàn)當(dāng)水分子采用SPC/E 模型時(shí)得到的平衡溫度與實(shí)驗(yàn)值更加接近。相比SPC/E-UA，SPC/E-AA 模型得到的平衡溫度更精確，這是因?yàn)門(mén)KM-AA 全原子模型對(duì)CH4分子的各個(gè)原子均設(shè)置點(diǎn)電荷，以再現(xiàn)CH4氣體分子的電八極矩，同時(shí)將L-J 勢(shì)分配給CH4的中心碳原子，而對(duì)氫原子不設(shè)置L-J相互作用，與OPLS-UA 模型相比，更能精確地對(duì)甲烷分子進(jìn)行描述。

由于甲烷水合物的分解為吸熱反應(yīng)，其反應(yīng)方程式可以表示為

其中甲烷水合物的分解熱ΔHd可由Clausius-Clapeyron方程[34]求得，如下所示

式 中，Z 為 壓 縮 因 子，可 由Peng-Robinson 方程[35]求得，結(jié)合p-T 的對(duì)應(yīng)關(guān)系[36]，進(jìn)而可以求得采用TIP4P-AA、TIP4P-UA、SPC/E-AA 和SPC/E-UA模型時(shí)甲烷水合物的分解熱分別為53.85、53.27l、54.79 及56.68 kJ/mol。Rueff 等[37]運(yùn)用DSC 量熱法測(cè)得的甲烷水合物的分解熱為54.49 kJ/mol，孫志高等[38]由恒溫壓力搜索法測(cè)得的水合物的分解熱為56.22 kJ/mol，由此可見(jiàn)當(dāng)水分子采用SPC/E 模型計(jì)算甲烷水合物的分解熱時(shí)結(jié)果與實(shí)驗(yàn)值更加吻合。

2.2 分解過(guò)程F3序參數(shù)的變化

圖2 不同模型的平衡溫度Fig.2 Equilibrium temperature under different models

在模擬過(guò)程中，甲烷水合物逐層向內(nèi)分解，分解前沿區(qū)域由外向內(nèi)推移，水合物-液態(tài)水界面向水合物側(cè)移動(dòng)，可以通過(guò)計(jì)算氧原子的序參數(shù)[26]來(lái)區(qū)分類(lèi)水合物水和液態(tài)水，進(jìn)而確定水合物-液態(tài)水的界面位置。序參數(shù)F3的定義如下

由于水合物相水分子通過(guò)氫鍵形成四面體式排列，在理想水合物結(jié)構(gòu)中鄰近水分子的氧原子形成的夾角為104.25°。因此，可以通過(guò)任一時(shí)刻鄰近水分子中氧原子夾角的測(cè)量，來(lái)衡量水合物中氫鍵網(wǎng)絡(luò)偏離理想四面體式排列的程度。理想狀態(tài)下水合物相水分子的F3約為0，而液相中水分子的排列具有無(wú)序性，F(xiàn)3約為0.1[39]。圖3 為不同模型下體系中水分子的F3沿Z 軸的分布情況。可以判定F3≈0.1 的區(qū)域?yàn)橐簯B(tài)水區(qū)域，F(xiàn)3≈0 的區(qū)域?yàn)樗衔飬^(qū)域，在液態(tài)水區(qū)域至水合物區(qū)域的過(guò)渡段F3出現(xiàn)跳躍式下降，該過(guò)渡段被界定為界面區(qū)域。圖3 中用虛線對(duì)液態(tài)水、水合物及界面區(qū)域做了區(qū)分。可以看出，在初始時(shí)刻t=0 時(shí)水合物區(qū)域F3接近于0，水合物兩側(cè)界面區(qū)域分別位于35.3～41.5 ? 和到109.4～116.9 ?。分解至t=400 ps 時(shí)，采用TIP4PAA 模型組合時(shí)界面區(qū)域分別位于40.3～47.1 ? 和98.7～105.2 ?。采用TIP4P-UA 模型組合時(shí)界面區(qū)域分別位于41.83～46.9 ? 和101.6～106.7 ?，同時(shí)水合物區(qū)域的F3較初始時(shí)刻有所上升，這是由于隨著分解的進(jìn)行，液態(tài)水區(qū)域向水合物區(qū)域傳遞熱量，水合物內(nèi)部水分子活躍程度增加，籠形結(jié)構(gòu)出現(xiàn)扭曲變形造成的；而采用SPC/E-AA 模型組合時(shí)界面區(qū)域分別為48.7～55.1 ? 和88.7～94.8 ?，與TIP4P-AA 模型相比，該模型組合下，界面區(qū)域向內(nèi)推移，且水合物區(qū)域F3與TIP4P-AA 模型相比略微升高。這主要是由于TIP4P 水模型與SPC/E 水模型所攜帶的電荷量不同引起的。TIP4P 水中氧原子與氫原子的電荷量大于SPC/E 水中氧原子與氫原子的電荷量，氧原子與鄰近氫原子之間的庫(kù)侖力更大，形成的氫鍵更穩(wěn)固，因而籠形結(jié)構(gòu)不易遭到破壞；此外，采用SPC/E-UA 模型組合時(shí)界面區(qū)域分別為57.1～62.1 ? 和97.3～103.6 ?。Sloan 等[1]預(yù)測(cè)的水合物-水的界面厚度約為5 ?，而Myshakin 等[40]估算的界面厚度在5～7 ? 之間，以上通過(guò)F3預(yù)測(cè)的各水合物-水界面區(qū)域厚度均在5～7.5 ? 之間，這與Sloan 等[1]及Myshakin 等[40]的結(jié)果相吻合，從而也驗(yàn)證了此方法的正確性。

圖3 氧原子序參數(shù)沿Z軸變化Fig.3 F3 distribution along Z axis

圖4為不同模型下水合物中水的序參數(shù)F3在分解過(guò)程中的變化特征。由圖4 可見(jiàn)，四種不同模型下甲烷水合物在分解達(dá)到2000 ps 時(shí)整體的F3均未達(dá)到0.1，說(shuō)明至模擬結(jié)束時(shí)依然有部分水合物尚未分解。在0～2000 ps 內(nèi)，TIP4P 模型預(yù)測(cè)的F3變化比較平緩，分解結(jié)束時(shí)增至0.04 左右，說(shuō)明采用該模型時(shí)水合物中水分子的有序度更高，甲烷水合物分解相對(duì)較少。采用SPC/E模型時(shí)氧原子的F3曲線斜率大于采用TIP4P 模型時(shí)氧原子的F3曲線斜率，說(shuō)明SPC/E 模型下甲烷水合物的分解速率大于TIP4P 模型的預(yù)測(cè)結(jié)果?？蓮膭?shì)函數(shù)角度分析SPC/E 模型的快速分解機(jī)理。由表1 可知，SPC/E 模型的σ與TIP4P模型相差很小，然而與TIP4P相比，SPC/E模型的勢(shì)阱深度ε 更小，因而分子間范德華力作用較小，同時(shí)氧原子與氫原子攜帶電荷較少導(dǎo)致氫鍵作用小，因而SPC/E 模型中單個(gè)水分子受到的束縛更小，水分子更加靈活。保持水分子模型不變，對(duì)比甲烷分子的AA 及UA 模型，發(fā)現(xiàn)F3變化特性基本一致，說(shuō)明對(duì)甲烷分子的模擬選用TKM-AA 模型或OPLS-UA 模型對(duì)分解速率影響不大，這也表明在水合物的分解過(guò)程中水-水相互作用占主導(dǎo)地位。此外，圖4 還給出了水合物中水分子/甲烷分子的協(xié)同作用對(duì)分解的影響，由圖4可知，保持水分子模型及甲烷分子模型不變，當(dāng)水分子/甲烷分子協(xié)同作用程度增大為原來(lái)的1.5 倍時(shí)，F(xiàn)3曲線斜率明顯降低，因而水合物分解更慢。這主要是因?yàn)楫?dāng)主客體分子協(xié)同作用增強(qiáng)時(shí)，客體分子的局域化特征增強(qiáng)，主客體分子之間的耦合程度增強(qiáng)，能量傳遞阻力增大，因而籠形結(jié)構(gòu)穩(wěn)定性增強(qiáng)。

圖4 水合物中F3序參數(shù)隨時(shí)間的變化Fig.4 Time dependence of F3 in methane hdyrate

2.3 甲烷分子勢(shì)能的變化

由于每個(gè)甲烷分子受到的勢(shì)能源于模擬盒中所有其他分子對(duì)它的相互作用，根據(jù)圖3 中液態(tài)水區(qū)、水合物區(qū)、界面區(qū)的界定標(biāo)準(zhǔn)，分別對(duì)分解過(guò)程三個(gè)不同區(qū)域的單個(gè)甲烷分子的平均勢(shì)能進(jìn)行了計(jì)算。圖5 所示為單個(gè)甲烷分子勢(shì)能隨時(shí)間的變化。由圖5可見(jiàn)，隨著分解的進(jìn)行，甲烷分子的勢(shì)能逐漸增加，這是由于NVE 分解過(guò)程體系的動(dòng)能降低轉(zhuǎn)化為水合物分解所需熱量，從而導(dǎo)致了勢(shì)能的增加。對(duì)比不同區(qū)域的勢(shì)能可知液態(tài)水區(qū)單個(gè)甲烷分子勢(shì)能最大，界面區(qū)單個(gè)甲烷分子勢(shì)能次之，水合物區(qū)單個(gè)甲烷分子勢(shì)能最??；同時(shí)，不同區(qū)域甲烷分子勢(shì)能的增量也不一致，這也體現(xiàn)了溫度分布的不均勻性，即分解過(guò)程中不同區(qū)域存在溫度梯度和熱流梯度，從而驗(yàn)證了水合物分解過(guò)程傳熱的非平衡特性[18]。對(duì)比不同模型預(yù)測(cè)的勢(shì)能可知，SPC/E-AA 模型預(yù)測(cè)的單個(gè)甲烷分子勢(shì)能大于TIP4PAA 模型單個(gè)甲烷分子勢(shì)能。甲烷分子勢(shì)能越低，穩(wěn)定性越高，因而TIP4P分子模型穩(wěn)定性越高，這也導(dǎo)致了圖4中TIP4P模型預(yù)測(cè)的分解速度較慢。

2.4 分解過(guò)程密度分布的變化

圖6為初始時(shí)刻及分解結(jié)束時(shí)不同分子模型的構(gòu)象以及水合物中甲烷、水沿Z 軸的密度分布。在初始時(shí)刻，水合物中甲烷和水的密度分布比較規(guī)律，水的密度按照單雙峰呈現(xiàn)周期性變化，甲烷的密度按照高低峰呈現(xiàn)周期性變化。由圖6(b)ⅰ可知，對(duì)于完整的3×3×6型SI甲烷水合物而言，甲烷的密度分布共有12對(duì)高低峰組合。由于尚未分解，在左右兩側(cè)液態(tài)水區(qū)域甲烷和類(lèi)水合物水的密度均為0。隨著分解的進(jìn)行，籠形結(jié)構(gòu)不斷破裂，甲烷分子、水分子開(kāi)始從水合物中逸出，并擴(kuò)散到液態(tài)水中，甲烷和水分子密度分布的規(guī)律性逐漸減弱。由圖6(b)ⅱ～ⅴ可見(jiàn)，隨著分解的進(jìn)行，液態(tài)水區(qū)甲烷和水的密度開(kāi)始出現(xiàn)非零正值。分解結(jié)束時(shí)，對(duì)于TIP4P-AA 模型和TIP4P-UA 模型，由圖6（a）Ⅱ～Ⅴ可知，尚有四層水合物籠形結(jié)構(gòu)未分解，而對(duì)于SPC/E-AA 和SPC/E-UA 模型則只剩下約兩層沒(méi)有分解。因而，圖6(b)ⅱ～ⅴ中甲烷密度的高低峰殘存數(shù)分別為8 對(duì)、8 對(duì)、4 對(duì)、4 對(duì)。從圖6(a)Ⅳ所示的SPC/E-AA 模型預(yù)測(cè)的2 ns 時(shí)的構(gòu)象圖可知，甲烷分子發(fā)生聚集形成甲烷氣泡，這與圖6(b)ⅳ中甲烷的密度在34? 和108? 附近呈現(xiàn)峰狀結(jié)構(gòu)而水的密度在此處呈現(xiàn)谷狀結(jié)構(gòu)相對(duì)應(yīng)。這種現(xiàn)象的形成是由于隨著分解的進(jìn)行大量的甲烷分子快速釋放到液態(tài)水中，當(dāng)甲烷濃度超過(guò)液態(tài)水中甲烷的溶解度時(shí)，甲烷分子開(kāi)始聚集形成甲烷富集區(qū)。

圖5 單個(gè)甲烷分子勢(shì)能隨時(shí)間的變化Fig.5 Variation of potential energy of single methane molecule with time

2.5 分解過(guò)程甲烷分子數(shù)的變化特征

由于甲烷水合物中十四面體籠形結(jié)構(gòu)的直徑約為8 ?，結(jié)合Alavi 等[18]的研究方法，如果甲烷分子任一時(shí)刻的坐標(biāo)與初始時(shí)刻的坐標(biāo)距離超過(guò)8?，則認(rèn)為甲烷分子已經(jīng)從破碎的籠形結(jié)構(gòu)中逃逸出來(lái)。圖7為分解過(guò)程逸出的甲烷分子數(shù)的變化特征。由圖可見(jiàn)，SPC/E-AA 模型、SPC/E-UA 模型預(yù)測(cè)的逸出甲烷分子數(shù)大于TIP4P-AA 模型及TIP4P-UA 模型的預(yù)測(cè)結(jié)果，SPC/E-AA 模型與SPC/E-UA 模型預(yù)測(cè)的逸出甲烷分子數(shù)目基本保持一致，這意味著采用SPC/E 模型時(shí)水合物的分解速度快于TIP4P 模型，而甲烷分子模型對(duì)分解的影響甚小，這與根據(jù)圖4得出的結(jié)論一致。從另外三種模型下逸出甲烷分子數(shù)的變化曲線可以看出，分解初始階段甲烷的逸出速率較大，之后逐漸減緩。這是因?yàn)樵撃M在NVE 絕熱系綜進(jìn)行，水合物分解是一個(gè)吸熱過(guò)程[25,41]，隨著分解的進(jìn)行，NVE 系綜溫度逐漸下降，水合物分解所需驅(qū)動(dòng)力不足，傳熱傳質(zhì)阻力增大，因而分解速率不斷衰減。分解結(jié)束時(shí)，SPC/E-AA和SPC/E-UA 模型預(yù)測(cè)的逸出甲烷分子數(shù)約為290個(gè)，分解比例約為2/3，這與圖6 中分解結(jié)束時(shí)兩種模型的預(yù)測(cè)結(jié)果吻合較好，其中甲烷密度的殘存峰值對(duì)數(shù)為4對(duì)(殘存比例1/3)。對(duì)于TIP4P-AA 模型，紅色虛線標(biāo)記的位置處逸出的甲烷分子數(shù)出現(xiàn)下滑，這主要是由于甲烷水合物在分解過(guò)程中發(fā)生了二次生成效應(yīng)，當(dāng)甲烷分子濃度高且存在過(guò)冷驅(qū)動(dòng)力時(shí)逸出的甲烷分子會(huì)重新與水分子結(jié)合再次生成甲烷水合物。

根據(jù)Arrhenius 公式，甲烷水合物的分解速率可由甲烷分子數(shù)隨時(shí)間的變化來(lái)表示，如式(7)所示

由質(zhì)量守恒方程可得甲烷分子數(shù)與分解界面之間的關(guān)系為

將式(8)代入式(7)可得

根據(jù)式(7)～式(9)，計(jì)算可得甲烷水合物分解所需的活化能，其中采用SPC/E-AA 模型預(yù)測(cè)的活化能為83.8 kJ/mol，采用SPC/E-UA模型預(yù)測(cè)的活化能為84.1 kJ/mol，采用TIP4P-AA 模型預(yù)測(cè)的活化能為86.7 kJ/mol，采用TIP4P-UA 模型時(shí)活化能為86.2 kJ/mol。Clarke 等[42]測(cè)得甲烷水合物分解所需的活化能為81 kJ/mol，Myshakin 等[40]通過(guò)對(duì)水合物在不同溫度下反應(yīng)速度的擬合求得活化能為(82.2±2.1)kJ/mol，因而在NVE 系綜內(nèi)采用SPC/E 模型時(shí)所預(yù)測(cè)的活化能與上述測(cè)量值吻合較好。

2.6 擴(kuò)散系數(shù)

甲烷水合物在NVE 系綜分解時(shí)，原子從初始位置開(kāi)始不停移動(dòng)，ri(t)表示t 時(shí)刻時(shí)粒子的位置，粒子位移平方的平均值稱(chēng)為均方位移，即

當(dāng)t 值非常小時(shí)，R(t)呈指數(shù)增加；隨著t 值的增大R(t)近似直線分布。根據(jù)Einstein的擴(kuò)散定律

由式(11)可知，通過(guò)均方位移的斜率可求出擴(kuò)散系數(shù)。通過(guò)對(duì)甲烷水合物相平衡溫度及分解所需活化能的計(jì)算可知，當(dāng)水分子采用SPC/E 模型時(shí)的計(jì)算結(jié)果與實(shí)驗(yàn)值更加吻合，而TKM-AA 及OPLS-UA 模型對(duì)分解影響很小，因而此處采用SPC/E-AA 模型研究不同位置甲烷的擴(kuò)散規(guī)律。針對(duì)SPC/E-AA 模型預(yù)測(cè)的水合物的分解過(guò)程，將完整的3×3×6 水合物分為如圖1 所示的A、B、C 三個(gè)區(qū)域，A 區(qū)域?yàn)樽钔鈱铀衔?，C 區(qū)域?yàn)樽顑?nèi)層水合物，B 區(qū)域位于A 區(qū)域和C 區(qū)域之間。對(duì)這三個(gè)區(qū)域內(nèi)甲烷分子的均方位移進(jìn)行計(jì)算，可得到圖8 所示的不同區(qū)域甲烷的均方位移隨時(shí)間的變化關(guān)系。由圖8 可見(jiàn)，不同區(qū)域甲烷的均方位移存在很大差別，說(shuō)明分解過(guò)程中不同區(qū)域存在濃度梯度，體現(xiàn)了分解過(guò)程傳質(zhì)的非平衡特性[18]。與水接觸最近的A 區(qū)域內(nèi)甲烷分子的均方位移最大，對(duì)應(yīng)的甲烷分子擴(kuò)散系數(shù)為1.4244×10-5cm2/s。與A 區(qū)域相比，B區(qū)域甲烷分子的擴(kuò)散系數(shù)小一個(gè)量級(jí)，僅為6.6487×10-6cm2/s，因而B(niǎo) 層的傳質(zhì)阻力要遠(yuǎn)大于A 層。這主要是因?yàn)殡S著A 層的分解，在水合物和水的界面處形成一層準(zhǔn)液態(tài)膜，阻止B層內(nèi)甲烷分子的擴(kuò)散。此外，由圖6(a)Ⅳ體系構(gòu)象圖中可見(jiàn)，分解結(jié)束時(shí)C區(qū)域尚未分解，因此甲烷分子只是在平衡位置附近振動(dòng)及轉(zhuǎn)動(dòng)，其均方位移基本恒定在一個(gè)稍微大于零的值附近，均方位移隨時(shí)間的變化曲線斜率為0，這表明C 區(qū)域甲烷的擴(kuò)散系數(shù)為0。此外，由菲克定律可知：M = -D?ρ。其中等號(hào)左側(cè)項(xiàng)代表了質(zhì)量的傳遞率，右側(cè)第一項(xiàng)代表擴(kuò)散系數(shù)，右側(cè)第二項(xiàng)代表質(zhì)量濃度梯度，即傳遞過(guò)程的推動(dòng)力。從圖6 中甲烷的質(zhì)量密度的變化規(guī)律可以看出A 區(qū)域甲烷的質(zhì)量密度梯度最大，B 區(qū)域次之，因而甲烷水合物分解過(guò)程中B 區(qū)域分解時(shí)的推動(dòng)力小于A 區(qū)域的分解推動(dòng)力，同時(shí)根據(jù)圖8 的計(jì)算結(jié)果，A區(qū)域擴(kuò)散系數(shù)最大，B 區(qū)域次之，進(jìn)而導(dǎo)致A 區(qū)域分解時(shí)質(zhì)量的傳遞率最大，B 區(qū)域的質(zhì)量傳遞率小于A 區(qū)域。這也驗(yàn)證了圖7 逸出的甲烷分子隨時(shí)間的變化曲線中起始段斜率較大，之后逐漸減小的變化趨勢(shì)。

圖6 體系構(gòu)象和質(zhì)量密度Fig.6 System configuration and density distribution

圖7 逸出籠子的甲烷分子總數(shù)隨時(shí)間的變化Fig.7 Number of methane molecules fleeing out of methane hydrate as function of time

圖8 SPC/E-AA模型不同區(qū)域甲烷均方位移隨時(shí)間的變化Fig.8 Mean square displacement of methane at different zone under SPC/E-AA model

3 結(jié) 論

(1) 在初始溫度300 K、初始?jí)毫?0 MPa 條件下，采用SPC/E 模型模擬得到的NVE 系綜內(nèi)水合物分解的平衡溫度、活化能、分解熱與實(shí)驗(yàn)測(cè)量結(jié)果更加接近。

(2)根據(jù)不同相態(tài)水的F3序參數(shù)，給出了區(qū)分液態(tài)水、水合物水及界面的方法，判定的液態(tài)水-水合物界面區(qū)域厚度為5～7.5 ?。對(duì)液態(tài)水、界面、水合物區(qū)甲烷的勢(shì)能進(jìn)行計(jì)算，發(fā)現(xiàn)甲烷分子在水合物區(qū)的勢(shì)能最低，界面區(qū)勢(shì)能次之，液態(tài)水區(qū)勢(shì)能最高，體現(xiàn)了NVE系綜分解過(guò)程的非平衡傳熱特征。

(3)通過(guò)采用不同模型對(duì)分解過(guò)程中F3序參數(shù)、甲烷和水的密度分布及甲烷的逸出規(guī)律的模擬對(duì)比，說(shuō)明甲烷分子模型選用TKM-AA 或OPLS-UA模型對(duì)分解的影響甚微，而不同水分子模型對(duì)分解影響較為明顯。由于TIP4P 模型下分子勢(shì)能較低，結(jié)構(gòu)較穩(wěn)定，相同溫度、壓力條件下SPC/E 模型預(yù)測(cè)的分解速率大于TIP4P模型。

(4)計(jì)算得到不同區(qū)域甲烷分子的擴(kuò)散系數(shù)范圍為10-6～10-5cm2/s，不同區(qū)域甲烷的擴(kuò)散系數(shù)差距明顯，由外向內(nèi)傳質(zhì)阻力逐漸增加，體現(xiàn)了NVE 系綜下分解過(guò)程的非平衡傳質(zhì)特征。

符 號(hào) 說(shuō) 明

As——甲烷水合物與液態(tài)水接觸的截面積，m2

cv——水的比定容熱容，J/(kg·K)

D——擴(kuò)散系數(shù)，cm2/s

ΔHd——水合物的分解熱，kJ/mol

K0——反應(yīng)常數(shù)

I——分解界面的位置，?

M——質(zhì)量通量密度，kg/(m2·s)

MSD——均方位移，?2

mw0——液態(tài)水的初始質(zhì)量，kg

nCH4——甲烷分子數(shù)

p——平衡壓力，Pa

Q+,Q-——原子所帶的正負(fù)電荷量，e

qi,qj——原子所帶電荷，e

R——摩爾氣體常數(shù)，J/(mol·K)

ri(t)——t時(shí)刻時(shí)粒子的位置

T——平衡壓力對(duì)應(yīng)的平衡溫度，K

Teq——體系的平衡溫度，K

t——時(shí)間，ps

U(rij)——原子i與j之間的相互作用勢(shì)，4.18×103J/mol

Z——壓縮因子

ε——?jiǎng)葳迳疃龋?.18×103J/mol

ε0——真空介電常數(shù)，F(xiàn)/m

θjik——原子j，i，k形成的夾角

ρ——質(zhì)量密度，kg/m3

σ——原子間平衡間距，?

χ——標(biāo)度系數(shù)，用來(lái)調(diào)整范德華作用力強(qiáng)弱下角標(biāo)

i，j，k——氧原子編號(hào)