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

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

引 言

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

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

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

1 計算方法

1.1 模擬體系

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

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

1.2 模型選取

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

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

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

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

1.3 模擬方法

采用LAMMPS 軟件進行模擬。采用速度-Verlet 算法積分牛頓運動方程，設(shè)置時間步長為1 fs，分子間相互作用的截斷半徑設(shè)置為12 ?，采用PPPM(particle-particle particle-mesh)法處理長程庫侖力，計算精度為10-6。系統(tǒng)的初始溫度設(shè)置為300 K，壓力設(shè)置為10 MPa，模擬步驟如下。

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

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

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

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

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 甲烷分子勢能的變化

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

2.4 分解過程密度分布的變化

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

圖5 單個甲烷分子勢能隨時間的變化Fig.5 Variation of potential energy of single methane molecule with time

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

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

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

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

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

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

2.6 擴散系數(shù)

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

符 號 說 明

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

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

D——擴散系數(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-——原子所帶的正負電荷量，e

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

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

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

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

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

t——時間，ps

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

Z——壓縮因子

ε——勢阱深度，4.18×103J/mol

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

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

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

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

χ——標度系數(shù)，用來調(diào)整范德華作用力強弱下角標

i，j，k——氧原子編號