吉夢霞，梁海峰，張 強(qiáng)，楊 琨

（太原理工大學(xué) 化學(xué)化工學(xué)院，山西 太原 030024）

天然氣水合物又稱甲烷水合物，其主要由天然氣和水在低溫高壓下形成的籠型化合物，主要存在于海底沉積層和永久凍土區(qū)。 由于其儲(chǔ)量巨大、能量密度高且清潔無污染，對其開采已成為解決當(dāng)前能源問題的重要選擇。 天然氣水合物開采方法主要有降壓、注熱和注入化學(xué)抑制劑法，其中降壓法是通過泵吸作用降低環(huán)境壓力使天然氣水合物分解，也是最直接經(jīng)濟(jì)的開采方法。 中、日、美等多個(gè)國家已經(jīng)進(jìn)行了天然氣水合物試開采，并取得了一定的進(jìn)展[1-3]。 但商業(yè)開采還面臨一系列問題，目前主要通過實(shí)驗(yàn)和數(shù)值計(jì)算來對水合物分解特性、影響因素及分解對巖層力學(xué)特性影響進(jìn)行研究。 多孔介質(zhì)骨架結(jié)構(gòu)和水合物的賦存狀態(tài)是重要影響因素。Zhang等[4]、Zhan等[5]構(gòu)建了不同粒徑和孔徑的多孔介質(zhì)，分析了多孔介質(zhì)對水合物分解特性、甲烷回收率以及穩(wěn)定壓力的影響；Wen等[6]對多孔介質(zhì)中甲烷水合物分解進(jìn)行了可視化微觀觀測，探究了顆粒尺寸對甲烷水合物完全分解時(shí)間的影響。 Wang等[7]通過實(shí)驗(yàn)研究粘結(jié)狀態(tài)與孔隙填充狀態(tài)下的水合物分解對沉積物變形的影響，分析水合物形態(tài)和分布對產(chǎn)氣和平均傳熱速率的影；Avinash等[8]提出了預(yù)測儲(chǔ)層中天然氣水合物形成和分解的數(shù)學(xué)模型，模型中考慮了多孔介質(zhì)和海水對天然氣水合物動(dòng)力學(xué)的綜合影響，可以反映多孔砂粒及其內(nèi)部孔的尺寸和形狀的不規(guī)則性與水合物分解的關(guān)系；Xue等[9]通過建立孔隙網(wǎng)絡(luò)模型研究孔隙率與孔道微觀結(jié)構(gòu)對流體在孔隙中流動(dòng)行為的影響，為天然氣水合物的分解研究提供新思路。 Chen等[10]研究了不同飽和度水合物的分解速率， 并分別探究懸浮和包裹狀態(tài)下的水合物分別對應(yīng)的飽和度-滲透性關(guān)系。 Hou等[11]進(jìn)行了礦物顆粒和水合物的孔隙分布對滲透率變化影響的介觀模擬；喻西崇等[12]用格子Boltzmann模型對多孔介質(zhì)中水合物生成與分解過程中飽和度的變化對多孔介質(zhì)滲透率的影響進(jìn)行了模擬。

目前研究多孔介質(zhì)形態(tài)對水合物分解產(chǎn)生的影響多為實(shí)驗(yàn)研究，對于數(shù)值模擬而言，微觀和宏觀模擬都存在各自的局限性。 宏觀模擬關(guān)注水合物以及多孔介質(zhì)的物性， 忽略了多孔介質(zhì)的結(jié)構(gòu)；而微觀方法側(cè)重分解過程水合物的穩(wěn)定性，忽略了孔隙結(jié)構(gòu)與水合物的生成方式；而介觀模擬為滲透率研究，尚未與水合物的分解結(jié)合分析各影響因素。本文采用格子Boltzmann方法（Lattice Boltzmann Method,LBM）模擬多孔介質(zhì)中甲烷水合物分解，研究多孔介質(zhì)的孔隙率和骨架結(jié)構(gòu)以及水合物懸浮和包裹兩種狀態(tài)下對分解的影響。

1 多孔介質(zhì)模型的構(gòu)建

本研究通過改變多孔介質(zhì)顆粒的半徑，構(gòu)造孔隙率不同的多孔介質(zhì)，模擬孔隙率對多孔介質(zhì)中水合物分解的影響；通過改變多孔介質(zhì)顆粒的骨架結(jié)構(gòu)， 模擬孔隙結(jié)構(gòu)對多孔介質(zhì)中水合分解的影響；通過模擬水合物懸浮和包裹兩種賦存狀態(tài)，分析不同賦存狀態(tài)對水合物分解的影響。

1.1 多孔介質(zhì)孔隙率及骨架結(jié)構(gòu)

孔隙率和骨架結(jié)構(gòu)是多孔介質(zhì)的重要參數(shù)，對于正方形規(guī)則排列的多孔介質(zhì)（如圖1所示）整體與局部孔隙率的大小一致，令ε為：

孔隙率φ可以表達(dá)為[13]：

式中，Dp為圓的直徑，m；D為圓心距，m。

通過改變特征單元的形狀得到不同的骨架結(jié)構(gòu)。 在相同的孔隙率下特征單元呈正方形比菱形具有更大的滲流尺寸、更小的迂曲度。

圖1 多孔介質(zhì)顆粒正方形排列

1.2 水合物賦存狀態(tài)

水合物在多孔介質(zhì)中的賦存狀態(tài)可以分為兩類：（1）水合物作為骨架結(jié)構(gòu)的一部分；（2）水合物作為孔隙充填物。 分別選取其中的懸浮和包裹為例，如圖2(a)、(b)所示。

圖2 水合物賦存狀態(tài)

2 甲烷水合物分解動(dòng)力學(xué)模型

對甲烷水合物分解過程做以下簡化處理：在整個(gè)分解過程中不考慮熱量傳遞；不考慮水合物的二次生成；并將多孔介質(zhì)與水合物均簡化為圓形。

2.1 Kim-Bishnoi動(dòng)力學(xué)模型

通過降低環(huán)境壓力促使甲烷水合物分解，采用Kim-Bishnoi動(dòng)力學(xué)模型對甲烷水合物分解進(jìn)行模擬[14]。 水合物分解的動(dòng)力為三相平衡狀態(tài)的逸度fed與甲烷逸度f的差值， 水合物分解速率可以表示如下：

式中，nt為水合物中CH4的量，mol；t為時(shí)間，s；kd為水合物動(dòng)力學(xué)系數(shù)，mol/（m2·Pa·s）；As為水合物總體表面積，m2；fed為三相平衡逸度，Pa；f為甲烷逸度，Pa。

kd和平衡壓力pe可以表達(dá)為：

式中，k0為水合物分解的本征動(dòng)力學(xué)常數(shù)，mol/（m2·Pa·s）；R為摩爾氣體常數(shù)，J/(mol·K）；ΔEa為反應(yīng)的活化能，kJ/mol；Te為平衡溫度，K。

2.2 反應(yīng)邊界條件

甲烷水合物的分解反應(yīng)過程可以看作水合物表面消融的過程，水合物分解速率表達(dá)式為[15]：

式中，G為甲烷氣體密度，kg/m3；M為水合物分解產(chǎn)生氣體的摩爾質(zhì)量，g/mol；為氣體黏度，m2/s。

Zhang等[16]將濃度梯度表示為：

式中，D為擴(kuò)散系數(shù)，m2/s；CR為邊界點(diǎn)處的濃度，mol/L；CF為相鄰點(diǎn)處的濃度，mol/L；|Δx|為兩點(diǎn)之間的距離，m。

邊界點(diǎn)處濃度表達(dá)為：

式中，cs為格子聲速，m/s。

結(jié)合邊界濃度得到未知點(diǎn)處的分布函數(shù)：

2.3 VOP表面更新

對于水合物在分解過程中邊界的處理，本文采用固相更新的VOP方法[17]，即將水合物區(qū)域初始體積設(shè)定為1，其余部分初始體積設(shè)定為0，隨著分解反應(yīng)的進(jìn)行，水合物區(qū)域體積不斷地減小，當(dāng)減小為0時(shí)分解反應(yīng)結(jié)束。 水合物區(qū)域體積隨反應(yīng)時(shí)間變化的關(guān)系式如下：

水合物在分解過程中通過下式更新：

3 大密度比偽勢模型

偽勢模型也稱作Shan-Chen模型，該模型通過偽勢函數(shù)描述兩相間的相互作用力，格子Boltzmann方程描述了具有離散速度的流體粒子運(yùn)動(dòng)的過程，以D2Q9模型為例， 假設(shè)流體中分別存在不同組分，各組分的例子分布函數(shù)為：

式中，x為格子上的一個(gè)格點(diǎn)；ci為流體粒子的速度集合；Δt為離散時(shí)間步長（格子單位）；t為當(dāng)前時(shí)間步長（格子單位）；fki為速度分布函數(shù)；fkieq為某一組分的平衡態(tài)分布函數(shù)；為粒子的弛豫時(shí)間（格子單位），與黏度相關(guān)。

粒子分布函數(shù)與混合流體的密度關(guān)系如下：

式中，ρk為粒子的密度，kg/m3；ρ為混合流體的密度，kg/m3。

根據(jù)Shan和Chen[18]有效密度函數(shù)可以用指數(shù)形式表達(dá)為：

通過勢函數(shù)得到狀態(tài)方程：

式中，a、b是參數(shù)， 其值分別取1、2；R為氣體常數(shù)，J/(mol·K）。

通過式(16)，可以得到勢函數(shù)表達(dá)式如下：

為了可以形成穩(wěn)定的兩相，Zhao[19]引入了新的勢函數(shù)來描述兩相間的相互作用力：

式中，ρ1為液相密度，kg/m3；ρ2為氣相密度，kg/m3；0＜a0＜1，ρ10、ρ20取決于模擬的條件,在本文中，a0=0.005，ρ10=-0.0008/log(a0)，ρ20=0.0003。

得到新的作用力表達(dá)式如下：

式中，參數(shù)G11是液相間的相互組用力，N；G12是液相與氣相間的相互作用力，N；G22是氣相間相互作用力，N；w為權(quán)系數(shù)；e為格子步長。

通過以上描述，得到有效模擬具有高密度比的兩相流動(dòng)模型。

Shan和Chen[20]對基于偽勢模型改進(jìn)，平衡速度被重新定義為：

其中為了保持無作用力碰撞時(shí)的動(dòng)量守恒，即：

可以得到u′的表達(dá)式為：

在Shan-Chen模型中速度u的表達(dá)式為：

為了減少誤差， 在Shan-Doolen模型中u取碰撞前后的平均值，即：

4 模擬方法及結(jié)果

在本文模擬中， 計(jì)算區(qū)域物理尺寸為：0.36 m×0.18 m，格子單位為360×180。 綜合上面所提到的多孔介質(zhì)構(gòu)造以及格子Boltzmann方法模擬多孔介質(zhì)中水合物的分解，主要步驟如下：（1）構(gòu)建大密度比偽勢LBM模型；（2）結(jié)合水合物分解動(dòng)力學(xué)模型，即Kim-Bishnoi動(dòng)力學(xué)模型；（3）采用LBM處理邊界，通過VOP更新相界面，直至水合物完全分解。

分析多孔介質(zhì)孔隙率以及骨架結(jié)構(gòu)和水合物賦存狀態(tài)對水合物產(chǎn)氣率和累計(jì)產(chǎn)氣量的影響，通過反應(yīng)過程中達(dá)姆科勒數(shù)（Da）的變化表征水合物的分解速率。達(dá)姆科勒數(shù)（Da）是表征化學(xué)反應(yīng)速率與擴(kuò)散速率比值的準(zhǔn)數(shù)，即達(dá)姆科勒數(shù)（Da）值越大反應(yīng)速率越大[21]。

式中，a為反應(yīng)的分子計(jì)算系數(shù)；γ為反應(yīng)速率，mol/(L·s)；L為模型長度，m；E為混合流體速度，m/s；C為反應(yīng)物濃度，mol/L。

4.1 分解過程

分解過程通過甲烷氣相密度變化反映，以多孔介質(zhì)半徑為15、水合物半徑為10（格子單位）為例。圖3、圖4分別為水合物在t=100、1000、2000（格子單位） 時(shí)分別處于懸浮和包裹狀態(tài)的分解過程示意圖，可以看到處于包裹狀態(tài)的水合物在2000（格子單位）時(shí)間步長時(shí)完全分解，而處于懸浮狀態(tài)下還有大量水合物未完成分解。

圖3 懸浮狀態(tài)下水合物分解氣相密度（格子單位）

圖4 包裹狀態(tài)下水合物分解氣相密度（格子單位）

分解過程中Da準(zhǔn)數(shù)較小，由此可以得出擴(kuò)散速率較快，反應(yīng)速率決定整個(gè)分解過程的快慢。 同時(shí)整個(gè)區(qū)域甲烷密度較為均勻，說明水合物在多孔介質(zhì)中反應(yīng)速率基本相同，水合物分解反應(yīng)均勻發(fā)生。

4.2 影響因素分析

4.2.1 多孔介質(zhì)孔隙率的影響

孔隙率是多孔介質(zhì)的重要參數(shù)，與多孔介質(zhì)顆粒的大小、形狀、排列等有關(guān)。 保持多孔介質(zhì)顆粒的圓心距不變，通過改變顆粒半徑改變孔隙率。 設(shè)定多孔介質(zhì)顆粒半徑R為10、13、15、18、20（格子單位，下同），對應(yīng)孔隙率分別為84.49%、73.79%、65.11%、49.76%和37.98%。產(chǎn)氣率和累計(jì)產(chǎn)氣量隨時(shí)間的變化規(guī)律如圖5所示。 在分解初期， 水合物總量一定時(shí)， 孔隙率較大的多孔介質(zhì)在分解前期產(chǎn)氣率較高，隨著分解的進(jìn)行孔隙率較大的多孔介質(zhì)產(chǎn)氣率反而低于孔隙率較小的多孔介質(zhì)。 考慮隨著反應(yīng)的進(jìn)行影響產(chǎn)氣率的主導(dǎo)因素由孔隙率變?yōu)樗衔锏目偙砻娣e，并且前期產(chǎn)生較多的甲烷氣體也會(huì)一定程度抑制甲烷水合物的后續(xù)分解。 綜合整個(gè)分解趨勢，較大的孔隙率更有利于水合物的分解。 由于初始水合物總量相同，所以分解完成后的累計(jì)產(chǎn)氣量一致。

孔隙率較大的多孔介質(zhì)中水合物較先完成分解，孔隙率大于50%時(shí)影響還不是非常明顯，孔隙率小于50%時(shí)，完全分解所需要的時(shí)間大大增加，半徑R=20的多孔介質(zhì)中水合物完全分解的時(shí)間是R=10的多孔介質(zhì)的3倍多。 分解過程中的Da準(zhǔn)數(shù)在分解初期達(dá)到最大值并迅速下降， 中期Da值較為穩(wěn)定，表明水合物的分解前期是擴(kuò)散控制，中期擴(kuò)散與反應(yīng)相當(dāng)，后期為化學(xué)反應(yīng)控制的過程。 在水合物分解中期孔隙率較大的多孔介質(zhì)Da值較高，表明孔隙率的影響在分解中期表現(xiàn)明顯，持續(xù)較高的分解速率使得孔隙率較大的多孔介質(zhì)中的水合物更快地完全分解。

圖5 懸浮狀態(tài)下孔隙率對水合物分解的影響

4.2.2 多孔介質(zhì)骨架結(jié)構(gòu)的影響

多孔介質(zhì)骨架結(jié)構(gòu)反映了孔隙的迂曲度，在模擬骨架結(jié)構(gòu)對水合物分解的影響時(shí)，保持多孔介質(zhì)孔隙率一定，僅改變顆粒的排列，在本文中以正方形和菱形排列為例。 分別模擬粒徑為10和15時(shí)兩種排列方式下水合物的產(chǎn)氣特性，如圖6所示。

圖6 懸浮狀態(tài)下排列方式對水合物分解的影響

從圖6可看出， 多孔介質(zhì)排列對產(chǎn)氣特性的影響不明顯，在分解前期，正方形排列在產(chǎn)氣率和一段時(shí)間內(nèi)累計(jì)產(chǎn)氣量增量上略大于菱形排列，此時(shí)控制反應(yīng)進(jìn)行的主要因素是多孔介質(zhì)骨架結(jié)構(gòu)，菱形排列的迂曲度更大，會(huì)在一定程度上阻礙甲烷氣體的流動(dòng)，減緩整個(gè)反應(yīng)進(jìn)程，但影響不是十分明顯。 隨著分解反應(yīng)的進(jìn)行迂曲度的影響逐漸減弱，水合物表面積的影響增強(qiáng)。

半徑R=10時(shí)菱形排列下水合物完全分解的時(shí)間是正方形排列的1.06倍，R=15時(shí)約為1.12倍，多孔介質(zhì)半徑的增加對產(chǎn)氣特性的影響更顯著，同時(shí)迂曲度也會(huì)增大，對水合物分解的抑制作用就更加明顯。 分解初期達(dá)姆科勒數(shù)達(dá)到最大值并迅速下降，隨著分解的進(jìn)行Da值趨于穩(wěn)定。 但中期菱形排列較正方形排列Da值更大，孔隙迂曲度的抑制作用主要表現(xiàn)在水合物分解的初期。

4.2.3 水合物賦存方式的影響

水合物的賦存狀態(tài)大致可以分為兩類：水合物作為骨架結(jié)構(gòu)的一部分與孔隙填充的一部分。 在水合物作為骨架結(jié)構(gòu)時(shí)，可以是膠結(jié)、包裹、支撐三種狀態(tài)。 本文模擬水合物孔隙填充與包裹兩種賦存狀態(tài)對水合物產(chǎn)氣特性的影響， 如圖7所示。 分解初期，包裹狀態(tài)下的水合物分解速率約是懸浮狀態(tài)下水合物的4倍。 包裹狀態(tài)下水合物的產(chǎn)氣率迅速下降直至水合物完成分解，累計(jì)產(chǎn)氣量也迅速達(dá)到最大值。 相比而言，懸浮狀態(tài)下的水合物無論是產(chǎn)氣率還是累計(jì)產(chǎn)氣量變化都較為平緩。 包裹狀態(tài)的水合物比表面積更大，更有利于分解反應(yīng)。

多孔介質(zhì)顆粒半徑R=15的水合物懸浮狀態(tài)完全分解時(shí)間是包裹狀態(tài)的10倍多，水合物的賦存狀態(tài)對產(chǎn)氣特性影響非常明顯。 通過分析水合物分解過程中達(dá)姆科勒數(shù)的變化，包裹狀態(tài)下的水合物Da準(zhǔn)數(shù)在短時(shí)間變化明顯，沒有明顯的穩(wěn)定階段。

圖7 賦存狀態(tài)對水合物分解的影響

5 結(jié)語

（1）甲烷水合物的分解過程與多孔介質(zhì)骨架結(jié)構(gòu)及水合物的賦存狀態(tài)等因素有關(guān)。 較大的孔隙率有利于水合物的分解，孔隙率小于50%時(shí)，分解時(shí)間顯著增加，Da準(zhǔn)數(shù)也明顯減小。 因此孔隙率小于50%時(shí)降壓分解不具有顯著優(yōu)勢， 可以考慮降壓注熱聯(lián)合開采縮短開采時(shí)間。

（2）較為狹窄的滲流通道會(huì)一定程度上阻礙水合物的分解， 阻礙作用在水合物分解初期明顯，實(shí)際開采中可以考慮前期注熱降壓聯(lián)合開采提高產(chǎn)氣率，中后期僅采用降壓開采降低成本。

（3）包裹狀態(tài)較懸浮狀態(tài)更有利于水合物分解，等量水合物懸浮狀態(tài)分解速率遠(yuǎn)遠(yuǎn)大于包裹狀態(tài)，在實(shí)際開采中可以優(yōu)先選擇水合物處于包裹狀態(tài)的區(qū)域開采；對不同的水合物賦存方式對應(yīng)選擇更加適合的開采方法同時(shí)兼顧效率與成本。