邊界傳熱條件下甲烷水合物降壓分解的LBM模擬

張 華，梁海峰，邵偉強，張錫彥

（太原理工大學 化學工程與技術學院，山西 太原 030024）

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

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

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

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

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

1 模型與模擬

1.1 甲烷水合物分解的動力學模型

Kim-Bishnoi水合物分解動力學模型被廣泛應用于多孔介質水合沉積物的模擬[23]，該模型水合物分解的驅動力為三相反應的平衡逸度fed與局部氣體的逸度f之差。水合物產氣率表示如下：

式中，nt為甲烷物質的量，mol；t為反應時間，s；As為水合物反應表面積，m2；kd為水合物反應動力學速率，mol/(m2·Pa·s)。

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

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

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

1.2.1 流體流動模型

在本研究中，因為水合物分解過程會產生水和甲烷氣體，故采用ZHAO等[24]在Shan-Chen模型基礎上改進的氣液兩相流流動大密度比偽勢模型，引入新的勢函數來表示氣液兩相間的相互作用力。各組分的密度分布函數表示如下：

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

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

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

D2Q9平衡分布函數表示如下：

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

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

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

新的總相互作用勢如下：

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

1.2.2 傳質模型

傳質方程中濃度分布函數表示如下：

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

1.2.3 傳熱模型

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

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

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

溫度場的熱擴散系數和宏觀溫度由以下確定：

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

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

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

1.2.4 邊界條件

在水合物液固界面處發生非均相反應，把水合物分解看成一級反應，參照文獻[27]中可以得出分解速率為：

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

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

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

分解界面處的濃度計算如下：

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

1.2.5 水合物分解固相更新法

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

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

水合物體積根據每一個時間步長進行更新：

1.2.6 模擬實現過程

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

（1）將單位無量綱化，初始化計算區域；

（2）基于Kim-Bishnoi模型，構建水合物分解的LBM流體流動、傳質和傳熱模型；

（3）更新水合物結構并計算相關信息；

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

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

圖1 多孔介質的孔隙結構Fig.1 Pore structure of porous media

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

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

2 模擬結果與討論

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

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

圖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 不同邊界溫度對甲烷水合物分解的影響

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

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

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

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

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

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

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

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

3 結論

基于LBM方法，建立了水合物分解的Kim-Bishnoi動力學模型，流體流動、傳熱和傳質模型。從邊界絕熱、邊界加熱兩個方面研究了不同影響因素下甲烷水合物的傳熱機理和產氣特性，得到以下結論。

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

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

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