水合物開發儲層溫度場擾動研究

董京楠，班凡生，車陽，曲帥，王子健，劉智強

1 中國石油集團工程技術研究院有限公司，北京 100081

2 中國石油大學(北京)石油工程學院，北京 102249

0 引言

天然氣水合物是以甲烷為主要成分，在高壓力、低溫度的條件下形成的一種籠型晶體，在全球范圍內有著巨大的儲量，也是公認的21世紀繼天然氣和石油能源之后最具開發潛力的清潔能源。在海底深水區和永久凍土深處廣泛存在，有著豐富的儲量[1-4]。在我國，除了東海、南海等海域均有大量的天然氣水合物礦藏，青藏髙原等凍土區域也已證實了天然氣水合物的存在，并成功獲取了水合物天然巖心樣品[5]。同時，據估計青藏高原凍土區天然氣水合物資源具有1.2×1012～2.4×1014m3的巨大的儲量[6]。

海洋天然氣水合物的開采不僅受溫度和壓力的影響，同時開采過程的地質構造作用、海平面的升降以及海底工程都會使得開發十分復雜?；谔烊粴馑衔镌诔爻合聵O其不穩定性的特性，目前所提出的天然氣水合物開釆方法基本都是利用降壓、升溫等手段打破沉積層水合物的相平衡條件，促使其分解，然后再收集溢出的天然氣。具體的開發手段主要包括降壓法、熱激發法、化學試劑法和CO2置換開采法等[7-13]。水合物開采過程，溫度壓力的控制決定著開采的效率，而開采對儲層的影響也不可忽視。由于海底表層沉積物大多屬于黏土礦物，其強度較低，儲層內的水合物晶體維持著基本骨架的穩定性。開采時溫度和壓力的變化會導致水合物分解成氣體和水，水合物對骨架的膠結作用降低，從而降低凍土骨架的強度，進而導致井壁失穩。而無論是海上還是陸地，水合物的大規模開釆必然會導致水合物沉積層的強度降低，繼而發生沉降，進而可能造成蓋層的失穩，甚至會引發地震等危害[14-17]。

對于水合物分解引起的一系列問題，其中包含相變的傳熱問題是其中首先要面對的。在經典的傳熱問題中，主要是考慮材料的本身熱導率、熱容對傳熱過程的影響。同時，大多數研究都是在常溫和高溫環境下的傳熱過程。而對于水合物的開發，第一其環境溫度較低，第二在開發過程中水合物會發生相變，同時吸收大量的熱，會嚴重影響儲層的溫度場。這些問題相對于傳統傳熱問題更加復雜，研究也相對較少，甚至是阻礙了水合物儲層開發的關鍵理論基礎之一。

本文主要研究在海洋區域天然氣水合物開采時，井筒周圍的地層溫度場變化。通過理論分析建立了低溫條件下含相變的巖石傳熱控制方程。在控制方程中，將比熱和潛熱合并成一個等效熱容參數。基于Gibbs-Thomson方程[18-19]，提出了一種利用壓汞數據得到孔隙內含冰量隨溫度變化的方法。采用所提出的模型，通過類比水合物開發過程水合物相變與冰水相變的本質，以水為介質來模擬水合物開采過程中井筒溫度升高對周圍地層的溫度擾動，為進一步研究開采過程中儲層的穩定性提供儲層內相變傳熱的理論基礎。

1 理論模型的建立

通常情況下巖石內的熱量的傳遞由熱導率決定，可以寫出一般形式的熱傳導控制方程(即Poisson方程[20])：

其 中T為 溫 度(℃)，t為 時 間(s)，ρ為 巖 石 密 度(kg·m-3)，c為 巖 石 比 熱(J·kg-1·K-1)，k是 巖 石 熱 導 率(W·m-1·K-1)，?是梯度算子。而當孔隙中發生相變時，熱量不僅會以比熱的形式傳遞(體系溫度發生變化，但組分相態不變)，還會以潛熱形式發生轉移(體系組分發生相變，但溫度不變)。因此，含相變的控制方程可以寫為：

其中L是水的潛熱(3.33×105J·kg-1)，ρi是冰的密度(kg·m-3)，θi孔隙中冰的體積分數。當假設凍融過程都是瞬時的，將dθi/dt= (dθi/dT)(dT/dt)帶入上式，則有

其中ρc-Lρidθi/dT這項定義為等效熱容，從等式可以看出，孔隙含冰量隨溫度的變化使得傳熱過程變得復雜。

為了求解這個非線性方程，首先需要確定含冰量與溫度間的關系(即dθi/dT)。冰在孔隙中發生相變，其相變溫度不是一個常數，而是與孔隙的尺寸有關。巖石孔隙中的相變可以由Gibbs-Thomson方程表示[18]

其中ri是被冰晶所侵入的最小入孔半徑，S是物質的摩爾熵(J·mol-1·K-1)，Vi表示冰的摩爾體積(m3·mol-1)，γiw是冰水間的界面應力(0.04 J·m-2)，α是冰水間的接觸角(在這里取α=0)。Gibbs-Thomson方程揭示出在一定溫度下，水只有在大于某一臨界尺寸的巖石孔隙中才能夠發生相變，孔隙越小，其間發生相變所需的溫度越低。

由Gibbs-Thomson方程可以得到冰晶侵入孔徑尺寸與溫度的關系。但是需要注意的是，由于分離壓的存在，在孔隙中，冰晶與壁面會存在一層較薄的液膜。這層液膜的厚度同樣與溫度有關[18]

其中ξ是表示分子間力的作用的特征長度(在這里取ξ～2.3?)，D是 擴 散 系 數，并 且D=γsi-γsw-γiw(取D～ 0.33 J·m-2)。因此，發生相變的真實孔隙半徑的大小是要略大于ri的(如圖1)，并且有以下關系

圖1 孔喉半徑與溫度壓力間關系。r(T)：考慮液膜厚度的Gibbs-Thomson方程表示的溫度和冰晶侵入孔徑關系；r(P)：壓汞法中Washburn方程表示的汞壓與孔徑關系Fig.1 Relation between the smallest pore access radius and temperature.r(T): relation between temperature and pore radius by Gibbs-Thomson equation considering the premelting film (G-T+e); r(P): relation between the pore radius and Hg pressure in mercury intrusion porosimetry (MIP).

所以，如果已知巖石的孔徑分布，含冰量與溫度間的關系就可以由方程(4)～(6)給出。在實驗方法中 ，孔徑分布數據一般是由壓汞法取得，壓汞法的原理就是Washburn方程

其中p是汞壓(MPa)，σHg是汞的表面張力(0.48 mN·m-1)，αHg是汞與固體壁面的接觸角(取cosαHg=0.765)。

結合方程(4)～(7)，可以得到一定孔徑下相變溫度與汞壓間的關系，

由上式可以看出，含冰量隨溫度的變化可以用壓汞法的數據來代替，即通過方程(8)將汞壓轉換成溫度。然而，由于液膜的存在，壓汞法求得的孔隙體積分數(θHg=VHg/V)要比同一孔隙中的含冰量(θi=Vi/V)要大。假設孔隙都是圓柱體孔隙，對應一定孔徑的含冰量與真實孔隙體積之間的關系為:

應該注意的是，以上關于冰含量與溫度間關系的推導是基于在每個孔隙中凍結和融化現象都獨立存在這樣一個假設。然而，由于結冰過程是在水中形成冰核，在此過程中產生冰水界面需要消耗表面能。因此，結冰過程往往是由外部孔隙向內逐步推進的，而不是在每個孔道中獨立發生。換句話說，結冰過程依賴于孔隙的連接方式，一個孔隙的凍結依賴于它的入口孔隙半徑大小。考慮一個入口孔徑較小的大孔，溫度降低時大孔不會先凍結，而需要等到入口孔喉結冰后，大孔才會結冰。相反，由于液膜的存在，融化過程不涉及新表面的生成，因此可以在每個孔隙中獨立發生。因此，上述模型僅適用于不受孔隙連通性影響的融化過程。

2 儲層的溫度場模擬

在水合物開采過程中，井筒內流體溫度相對較高，會對周圍儲層的溫度場產生擾動。針對這一問題，建立如圖2所示平面幾何模型。模型中心的井眼直徑設為200 mm，地層半徑為10 m且為均質地層。注意到水合物開采過程中，井筒內的流體通過熱對流方式將熱量傳遞給井筒壁面，其后再以熱傳導方式在儲層內傳遞，所以外邊界設為自由邊界，內邊界設為對流邊界。儲層初始溫度設為-4 °C，在使用隔熱管的情況下，井內流體大部分熱量被隔絕，那么設等效的流體溫度為5 °C。

圖2 有限元網格模型Fig.2 Finite element mesh model

利用實驗室得到的一組砂巖的壓汞實驗數據求得其孔徑分布曲線(如圖3a所示，這里用巖石A來表示，峰值孔徑約為10 μm)。同時，為了考慮其孔徑分布對傳熱過程的影響，將巖石A的孔徑分布曲線向右平移得到巖石B的孔徑分布圖(峰值孔徑約為50 μm)。再通過方程(8)和(9)得到巖樣A和B孔隙內含冰量與溫度間的關系(如圖3b所示)，巖石A和B的含冰量基本都在0 °C附近變化較快，推斷在此區間內的相變反應劇烈，而巖樣B含冰量大于巖石A則是由于巖樣B的孔隙度較大。結合表1給出的各組分物質的物性參數以及傳熱過程的相關系數，得到等效熱容值并代入，對地層的傳熱過程進行模擬。

表1 物性參數和傳熱系數[21]Table 1 Physical parameters and thermal coefficient in this work[21]

圖3 a.巖樣A和B的孔徑分布曲線；b.孔隙內含冰量與溫度間關系Fig.3 a.The pore distribution curve of sample A and B; b.ice content vs temperature in pores

在模擬結果中，基于模型軸對稱的特點，以軸心為原點，選取任一半徑為X軸。在1小時到30天之間的不同時刻，模擬得到地層溫度場的演化云圖(圖4a)。取其沿X軸方向的地層巖石溫度分布，得到圖4b所示結果。在不同的t時刻，溫度都呈梯狀分布，沿X軸方向逐漸降低。同時，由于同一熱源，初始相鄰區域溫差較大，溫度的傳播更快，其溫度梯度更陡；隨著熱量向外的擴散，傳播距離逐漸增大，同時溫差逐漸減小，地層溫度分布梯度逐漸變緩，可以預測一定地層的溫度分布曲線是斜率逐漸變緩的過程。

圖4 a.溫度場演化；b.不同時刻地層溫度分布Fig.4 a.the evolution of temperature field; b.temperature distribution of permafrost at different moments

對于巖樣A和B，取井筒處同一點的溫度變化曲線(如圖5所示)。地層溫度隨時間逐漸升高，當達到零度附近時A，B曲線均出現比較平緩的一段，隨后溫度繼續上升。B曲線在A曲線上方，表明其溫度上升速率比A快。同時，假設溫度隨時間變化的速率小于0.01以視為溫度曲線的平段區域(dT/dt＜0.01)，圖中曲線B的平段要比曲線A長。

圖5 巖石A和B的溫度隨時間的變化關系Fig.5 Temperature vs time for sample A and B

3 分析與討論

在模擬結果中，相變潛熱對地層的溫度變化有極大的影響。在地層的傳熱過程中，孔隙內冰不斷發生相變，平段區間的出現正是由于其孔隙中相變過程較為明顯，大量的熱量轉化為物質的潛熱，導致溫度的上升極其緩慢。從控制方程(3)中也可以看出，在其他參數的一定的條件下，當含冰量隨溫度劇烈變化時，即潛熱項的突然增大，等同于增大了整體的熱容值，系統吸收相同的熱量時其溫度變化遠沒有之前明顯。對于巖樣A和B，在相似的孔徑分布下，巖樣B的孔徑較大，其所對應的相變溫度就較大，即在相同的熱源下，巖樣A比B更早的發生相變，即B曲線前期的升溫速率要大于A。同時，在2種巖樣的含冰量隨溫度的變化圖中，零度附近B曲線的含冰量隨溫度的變化曲率大于A曲線，即在較小的溫度區間內，巖樣B孔隙中的含冰量變化等效增加的體系熱容較大，所以升溫更慢，溫度曲線的平段更長。

由上可知，曲線的平段在模型中是由含冰量隨溫度的變化決定的，而含冰量隨溫度變化的關系是由壓汞實驗的孔徑分布數據獲得的，那么，通過測量材料在低溫下的溫度變化曲線，根據溫度曲線平段的位置和長度就可以得到物質的孔徑分布。值得注意的是雖然壓汞法是實驗室中一種較為常見的測量孔隙分布的方法，它的測量范圍廣，操作簡單，但在理論模型中，利用壓汞法來反求出巖石孔隙中的含冰量與溫度間關系進行了兩次變換，有液膜和溫度壓力間這兩個過程轉換誤差的影響。因此，在條件允許的情況下，用更加精確的方法測量巖石孔徑分布或者直接測出巖石孔隙中含冰量隨溫度的變化關系都可以對模型進行優化。

4 結論

針對水合物開發過程由于相變現象吸收大量的熱，從而會造成儲層溫度場的擾動這一現象，類比冰水相變進行了儲層溫度場的擾動模型研究，并研究了不同孔徑分布對儲層溫度隨時間變化的規律的影響。主要得到了以下結論：

(1) 建立了含相變過程的傳熱控制方程。其中，將潛熱和比熱2種能量轉化形式合并為等效熱容，從而使控制方程簡化為類似于不含相變過程的經典Poisson方程。另外，通過壓汞數據可以預測含冰量與溫度的關系。

(2) 當儲層溫度升高時，其地層溫度變化曲線會在0 °C附近出現一段緩慢上升的階段。這主要是因為在這一階段相變劇烈，需要吸收大量熱量，從而使溫度變化并不明顯。

(3) 儲層的溫度演化受巖石本身的孔隙結構影響巨大。當孔隙度較大時，由于骨架的缺失，巖石的整體熱容會相應減少，從而使溫度變化速率升高。同時，孔徑分布也對溫度平段的位置和長度產生。反之，通過測量儲層的溫度變化曲線，也可以反演得到巖石的孔徑分布。