張 華，梁海峰，邵偉強(qiáng)，張錫彥

（太原理工大學(xué) 化學(xué)工程與技術(shù)學(xué)院，山西 太原 030024）

水合物是由水分子組成的籠狀晶體，包裹一個(gè)或多個(gè)客體分子[1]。天然氣水合物也稱甲烷水合物，儲(chǔ)量巨大、燃燒清潔，對(duì)于解決全球能源和環(huán)境問題具有戰(zhàn)略意義，對(duì)其開采已成為國內(nèi)外研究熱點(diǎn)。水合物開采方法包括降壓法、注熱法、注化學(xué)劑法和二氧化碳置換法等[2-4]，其中降壓法、注熱法是相對(duì)實(shí)用的開采方式。

甲烷水合物儲(chǔ)存地質(zhì)條件復(fù)雜，目前主要是通過實(shí)驗(yàn)和數(shù)值模擬等方式開展研究。YOUSIF等[5]通過實(shí)驗(yàn)研究了多孔介質(zhì)中水合物的分解特性，分析了滲透率和孔隙率對(duì)水合物分解過程的影響。宋永臣等[6]模擬了天然氣水合物在降壓和邊界傳熱條件下的分解過程，對(duì)不同位置處壓力、溫度和產(chǎn)氣率等的變化進(jìn)行了研究。DAI等[7]建立了水合物沉積物在不同孔隙結(jié)構(gòu)和賦存狀態(tài)下的滲透率變化模型。WANG等[8-10]建立了減壓、注熱和減壓聯(lián)合注熱的解析模型，發(fā)現(xiàn)降壓與注熱聯(lián)合開采能夠有效地提高水合物的分解速率。WANG等[11]利用微流體裝置研究了孔隙中甲烷水合物的形成和分解，并且在形成過程中觀察到穩(wěn)定和不穩(wěn)定的水合物晶體。

上述研究均采用宏觀數(shù)值模擬方法，不能體現(xiàn)多孔介質(zhì)孔隙結(jié)構(gòu)及水合物賦存方式的影響。格子Boltzmann（LBM）方法從介觀角度出發(fā)，可以直接簡單地處理流體的流動(dòng)過程及流體與周圍環(huán)境之間的相互作用。LBM方法首先運(yùn)用于水合物滲透率研究的模擬中。YU等[12]采用LBM方法研究了水合物存在對(duì)多孔介質(zhì)滲透率的影響。喻西崇等[13]和JI等[14-15]采用LBM方法研究了不同飽和度水合物單相流和多相流的滲透率。HOU等[16]基于LBM方法研究了礦物顆粒和水合物孔隙尺度分布等因素對(duì)滲透率的影響，并提出了“控制滲流通道”的概念。DONG等[17]結(jié)合有限元和LBM方法模擬了水合物在巖石中的不同膠結(jié)方式，發(fā)現(xiàn)膠結(jié)方式對(duì)巖石的電阻率、彈性模量和滲透率有一定影響。

近年來LBM方法逐步被運(yùn)用于水合物分解過程的模擬。ZHANG等[18]基于LBM方法描述了甲烷水合物分解的相變和不同孔隙下的孔隙結(jié)構(gòu)演化，并對(duì)裂解鋒的推進(jìn)和溫度分布進(jìn)行了一些研究。吉夢霞等[19]和苗旭[20]基于LBM方法研究水合物分解特性，并將分解過程看作恒溫。YANG等[21-22]基于LBM方法建立了水合物分解孔隙尺度的數(shù)學(xué)模型，研究了多孔介質(zhì)中甲烷水合物解離過程的影響因素，分析了傳質(zhì)傳熱和水飽和度對(duì)甲烷水合物解離過程的影響。綜上，LBM方法運(yùn)用于水合物分解過程的模擬較少，且大多未考慮分解過程中的傳熱，以及水合物分解的固相消融影響。

本文基于LBM方法，建立描述甲烷水合物分解的動(dòng)力學(xué)模型、流體流動(dòng)模型以及傳質(zhì)和傳熱模型，運(yùn)用固相更新（VOP）法追蹤水合物固相變化，研究了不同邊界條件（加熱和絕熱）、不同影響因素（邊界溫度、孔隙率和賦存狀態(tài)）下，水合物的傳熱機(jī)理和產(chǎn)氣特性。

1 模型與模擬

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

Kim-Bishnoi水合物分解動(dòng)力學(xué)模型被廣泛應(yīng)用于多孔介質(zhì)水合沉積物的模擬[23]，該模型水合物分解的驅(qū)動(dòng)力為三相反應(yīng)的平衡逸度fed與局部氣體的逸度f之差。水合物產(chǎn)氣率表示如下：

式中，nt為甲烷物質(zhì)的量，mol；t為反應(yīng)時(shí)間，s；As為水合物反應(yīng)表面積，m2；kd為水合物反應(yīng)動(dòng)力學(xué)速率，mol/(m2·Pa·s)。

式中，k0為固有反應(yīng)動(dòng)力學(xué)常數(shù)，1.25×105mol/(m2·Pa·s)；ΔEa/R為9752.73 K；ΔEa為反應(yīng)活化能，J/mol；R為摩爾氣體常數(shù)，8.314 J/(mol·K)；T為反應(yīng)溫度，K。

一般情況下，式(1)中平衡逸度fed和局部氣體的逸度f用相應(yīng)溫度T的平衡壓力pe和甲烷氣體壓力p代替（單位為Pa），表示如下：

1.2 基于LBM的甲烷水合物分解模型

1.2.1 流體流動(dòng)模型

在本研究中，因?yàn)樗衔锓纸膺^程會(huì)產(chǎn)生水和甲烷氣體，故采用ZHAO等[24]在Shan-Chen模型基礎(chǔ)上改進(jìn)的氣液兩相流流動(dòng)大密度比偽勢模型，引入新的勢函數(shù)來表示氣液兩相間的相互作用力。各組分的密度分布函數(shù)表示如下：

式中，fki(x,t)為某組分k在位置x和時(shí)間t下第i個(gè)速度方向上的密度場分布函數(shù)；t為反應(yīng)時(shí)間，s；δt為時(shí)間步長；為平衡分布函數(shù)；ei為粒子的離散速度；τkf是粒子的無量綱弛豫時(shí)間。

在本研究中流體流動(dòng)采用二維的D2Q9，空間上9個(gè)速度方向的離散速度ei表示如下：

式中，c為晶格速度，c=δx/δt，δx為晶格間距。

D2Q9平衡分布函數(shù)表示如下：

式中，ρ為組分密度，kg/m3；權(quán)系數(shù)ω0=4/9，ω1～4=1/9和ω5～8=1/36；格子聲速。流體密度ρ和流速u的關(guān)系如下：

為了解決偽勢模型中氣液兩相大密度比帶來的不穩(wěn)定問題，兩組件的相互作用勢可表示如下：

式中，φ1為氣體的相互作用勢；φ2為液體的相互作用勢；0＜a0＜1，本研究中a0=0.005；ρ10=-0.0008/log(a0)；ρ20=0.0003。

新的總相互作用勢如下：

式中，G為組分間的相互作用力，N；ψ為勢函數(shù)；a和b為參數(shù)，分別取值1和2。

1.2.2 傳質(zhì)模型

傳質(zhì)方程中濃度分布函數(shù)表示如下：

式中，gki為組分k的濃度分布函數(shù)；為組分k的平衡濃度分布函數(shù)；τkg為濃度場組分k的無量綱弛豫時(shí)間；Dk為擴(kuò)散系數(shù)；Ck為組分濃度，mol/L。

1.2.3 傳熱模型

甲烷水合物分解的同時(shí)會(huì)不斷吸收熱量。本研究采用KARANI等[25]提出的附加源項(xiàng)方法，處理具有多個(gè)演化和移動(dòng)界面的甲烷水合物分解過程中的耦合傳熱。

式中，hi為某時(shí)間段的溫度分布函數(shù)；為平衡溫度函數(shù)；τh為溫度場的無量綱弛豫時(shí)間；ST為源相，K。

傳熱的D2Q5模型能夠準(zhǔn)確地恢復(fù)到宏觀方程且能提高計(jì)算效率，簡化的D2Q5晶格模型是通過放棄原始D2Q9模型對(duì)角方向（i=5～8）的離散速度而得到的，平衡的溫度函數(shù)表示如下：

溫度場的熱擴(kuò)散系數(shù)和宏觀溫度由以下確定：

在流體與水合物固體界面處，水合物逐漸被分解，周圍的熱量被吸收，而在其他的地方?jīng)]有分解，所以源相ST=0，ST的表達(dá)式如下：

計(jì)算水合物分解吸收的熱量采用MASUDA等[26]所給出的表達(dá)式：

式中，mH為水合物產(chǎn)氣率，g/s；MH為水合物的摩爾質(zhì)量，g/mol；參數(shù)r=56599 J/mol；參數(shù)d=16.744 J/(mol·K)。

1.2.4 邊界條件

在水合物液固界面處發(fā)生非均相反應(yīng)，把水合物分解看成一級(jí)反應(yīng)，參照文獻(xiàn)[27]中可以得出分解速率為：

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

ZHANG等[28]提出了一種LBM濃度邊界的一般回彈格式，濃度邊界條件和濃度梯度計(jì)算如下：

式中，CR為邊界點(diǎn)處的濃度，mol/L；CF為相鄰點(diǎn)流體點(diǎn)的濃度，mol/L；|Δx|為兩點(diǎn)間距離，m。

分解界面處的濃度計(jì)算如下：

確定邊界濃度后，根據(jù)ZHANG等[28]的反彈公式，可以確定邊界條件的分布函數(shù)為：

1.2.5 水合物分解固相更新法

模擬甲烷水合物分解的過程中，KANG等[29]采用開發(fā)的VOP法來不斷地更新甲烷固體水合物的體積變化。VOP法是一種前沿跟蹤方法，可用于跟蹤不斷演化的固液界面，廣泛應(yīng)用于反應(yīng)輸運(yùn)問題。本研究中，只有流體與水合物分解界面處的固體水合物需要進(jìn)行處理。將每個(gè)節(jié)點(diǎn)表示為一個(gè)控制體積，水合物節(jié)點(diǎn)的體積通過以下公式計(jì)算：

式中，As為表面反應(yīng)面積，m2；為摩爾體積，m3/mol ；V為水合物體積，m3。

水合物體積根據(jù)每一個(gè)時(shí)間步長進(jìn)行更新：

1.2.6 模擬實(shí)現(xiàn)過程

模擬甲烷水合物的分解過程，主要步驟如下：

（1）將單位無量綱化，初始化計(jì)算區(qū)域；

（2）基于Kim-Bishnoi模型，構(gòu)建水合物分解的LBM流體流動(dòng)、傳質(zhì)和傳熱模型；

（3）更新水合物結(jié)構(gòu)并計(jì)算相關(guān)信息；

（4）不斷重復(fù)步驟（2）、（3）直至水合物完全分解。

基于LBM法進(jìn)行水合物分解過程，分解、流動(dòng)、傳質(zhì)和傳熱機(jī)理的研究。建立多孔介質(zhì)模型，該模型計(jì)算的物理尺寸為0.22 m×0.10 m，采用220×100的格子單位，每一個(gè)圓形顆粒半徑為8（格子單位），如圖1所示。

圖1 多孔介質(zhì)的孔隙結(jié)構(gòu)Fig.1 Pore structure of porous media

孔隙率是模擬多孔介質(zhì)沉積物的重要因素之一，對(duì)于圖1中排列的多孔介質(zhì)模型，孔隙率計(jì)算式如下：

式中，ε為圓的直徑DP與圓心距D的比值。

2 模擬結(jié)果與討論

2.1 不同邊界傳熱方式對(duì)甲烷水合物分解的影響

圖2和圖3分別為懸浮狀態(tài)下邊界溫度為280.3 K和邊界絕熱下的溫度云圖。由圖2可知，當(dāng)邊界加熱時(shí)，邊界處和開口處的水合物由于受到邊界溫度和開口壓力的影響，率先分解完成，然后向中心推進(jìn)。邊界的熱量也隨著分解的進(jìn)行逐漸向中心傳遞，當(dāng)水合物完全分解，邊界的溫度并沒有均勻分布，邊界的傳熱范圍十分有限。由圖3可知，當(dāng)邊界絕熱時(shí)，溫度一直降低，直到水合物停止分解；受甲烷和水向出口處流動(dòng)的影響，開口處的溫度較高。

圖2 邊界溫度280.3 K下的溫度云圖Fig.2 Temperature cloud diagram with boundary temperature 280.3 K

圖3 邊界絕熱下的溫度云圖Fig.3 Temperature cloud diagram with boundary adiabatic

2.2 不同邊界溫度對(duì)甲烷水合物分解的影響

甲烷水合物分解是吸熱反應(yīng)，因此邊界傳熱對(duì)其有一定影響。圖4給出了懸浮狀態(tài)下不同邊界溫度對(duì)水合物分解的影響。由圖4(a)可知，不同邊界溫度下，水合物分解前期的累計(jì)產(chǎn)氣量和產(chǎn)氣率變化趨勢接近，說明前期主要由壓力控制，邊界溫度的影響很?。浑S著溫度的傳遞，邊界處溫度高分解完成的時(shí)間較短，且完全分解的累計(jì)產(chǎn)氣量為4925 cm3；但隨著邊界溫度增高，其對(duì)水合物分解的影響不斷減小。而邊界絕熱時(shí)，累計(jì)產(chǎn)氣量為1832.8 cm3，僅為完全分解時(shí)的3/8，因此在水合物降壓分解過程中需要通入充足熱量，以便水合物可以完全分解。由圖4(b)可知，當(dāng)邊界加熱時(shí)，溫度短暫上升，由于水合物分解吸收熱量大于邊界提供的熱量，溫度開始下降，隨著邊界加熱的不斷進(jìn)行，溫度又上升到邊界溫度。而邊界絕熱時(shí)，外界不能提供熱量，平均溫度下降到273.2 K左右后不再下降，因此在水合物開采過程中需要考慮水結(jié)冰和水合物再生成的情況。

圖4 不同邊界溫度對(duì)甲烷水合物分解的影響Fig.4 Effect of different boundary temperatures on decomposition of methane hydrate

2.3 不同孔隙率對(duì)甲烷水合物分解的影響

孔隙率是模擬多孔介質(zhì)沉積物的重要因素之一，反映了甲烷水合物飽和度和多孔介質(zhì)形態(tài)。模擬了邊界溫度278.3 K和邊界絕熱條件下，0.37、0.45和0.53這3種孔隙率對(duì)水合物分解的影響。由圖5(a)可知，當(dāng)邊界傳熱時(shí)，孔隙率較大產(chǎn)氣率下降速度較快，水合物分解完成時(shí)間較短，說明孔隙率較大有利于流體流動(dòng)和熱量傳遞，促進(jìn)了水合物分解。反應(yīng)在時(shí)間步長11000步之前，孔隙率較大的產(chǎn)氣率也大，然而隨著反應(yīng)的進(jìn)行，孔隙率較大的產(chǎn)氣率反而變小，這是由于孔隙率較大的產(chǎn)氣量較大，反應(yīng)表面積減小。而邊界絕熱時(shí)，由于熱量的影響，水合物較早停止分解，產(chǎn)氣率和累計(jì)產(chǎn)氣量變化趨勢幾乎一致，可以忽略孔隙率對(duì)水合物分解的影響。由圖5(b)可知，隨著孔隙率的增大，初始溫度下降的越快，最低點(diǎn)溫度也越低，達(dá)到最低點(diǎn)溫度的時(shí)間也不斷地提前，之后快速上升到邊界溫度。

圖5 不同孔隙率對(duì)甲烷水合物分解的影響Fig.5 Effect of different porosity on decomposition of methane hydrate

2.4 不同賦存方式對(duì)甲烷水合物分解的影響

甲烷水合物的賦存方式反映了水合物在地層中與周圍多孔介質(zhì)的相互關(guān)系，不同的賦存方式影響著水合物的基本特性。水合物在多孔介質(zhì)中的主要賦存形式可以分為兩種：懸浮和包裹（圖1）。本節(jié)模擬了懸浮和包裹兩種水合物賦存狀態(tài)下，邊界溫度278.3 K和絕熱邊界條件下的情況。由圖6(a) 可知，絕熱邊界條件下，包裹狀態(tài)累計(jì)產(chǎn)氣量為4126.9 cm3，懸浮狀態(tài)累計(jì)產(chǎn)氣量為1832.8 cm3，分別是完全分解總量的5/6和3/8，可見包裹狀態(tài)有利于提高水合物的產(chǎn)氣量。由圖6(b)可知，懸浮和包裹狀態(tài)下降的溫度接近，表明包裹狀態(tài)水合物更多是由壓力控制分解的。邊界加熱時(shí)，包裹狀態(tài)下的分解速率較快，但溫度的下降和上升也十分迅速，而且最低溫度下降到273.2 K左右。故在包裹狀態(tài)下開采時(shí)，應(yīng)注意溫度影響。在邊界加熱條件下，由于表面積的影響，包裹狀態(tài)水合物在時(shí)間步長3060步內(nèi)分解完全；而懸浮狀態(tài)水合物在時(shí)間步長30240步內(nèi)分解完全，遠(yuǎn)遠(yuǎn)大于包裹狀態(tài)。

圖6 不同賦存方式對(duì)甲烷水合物分解的影響Fig.6 Effect of different occurrence modes on decomposition of methane hydrate

3 結(jié)論

基于LBM方法，建立了水合物分解的Kim-Bishnoi動(dòng)力學(xué)模型，流體流動(dòng)、傳熱和傳質(zhì)模型。從邊界絕熱、邊界加熱兩個(gè)方面研究了不同影響因素下甲烷水合物的傳熱機(jī)理和產(chǎn)氣特性，得到以下結(jié)論。

（1）邊界加熱和邊界絕熱下水合物分解的累計(jì)產(chǎn)氣量分別為4925.0 cm3和1832.8 cm3，邊界加熱有利于提高水合物的分解速率和產(chǎn)氣量。隨著邊界溫度的增加，水合物完成分解所用時(shí)間變化并不明顯。

（2）孔隙率較大有助于流體流動(dòng)，從而加快邊界溫度的傳遞。而絕熱條件下，受溫度的影響，孔隙率對(duì)水合物影響較小。開采過程中，孔隙率較大時(shí)，需要同時(shí)考慮降壓和加熱開采。

（3）絕熱條件下，包裹狀態(tài)水合物累計(jì)產(chǎn)氣量為4126.9 cm3，懸浮狀態(tài)累計(jì)產(chǎn)氣量為1832.8 cm3，分別是完全分解總量的5/6和3/8。包裹狀態(tài)不僅可以縮短開采時(shí)間，還可以提高產(chǎn)氣量。