董京楠，班凡生，車陽(yáng)，曲帥，王子健，劉智強(qiáng)

1 中國(guó)石油集團(tuán)工程技術(shù)研究院有限公司，北京 100081

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

0 引言

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

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

對(duì)于水合物分解引起的一系列問(wèn)題，其中包含相變的傳熱問(wèn)題是其中首先要面對(duì)的。在經(jīng)典的傳熱問(wèn)題中，主要是考慮材料的本身熱導(dǎo)率、熱容對(duì)傳熱過(guò)程的影響。同時(shí)，大多數(shù)研究都是在常溫和高溫環(huán)境下的傳熱過(guò)程。而對(duì)于水合物的開(kāi)發(fā)，第一其環(huán)境溫度較低，第二在開(kāi)發(fā)過(guò)程中水合物會(huì)發(fā)生相變，同時(shí)吸收大量的熱，會(huì)嚴(yán)重影響儲(chǔ)層的溫度場(chǎng)。這些問(wèn)題相對(duì)于傳統(tǒng)傳熱問(wèn)題更加復(fù)雜，研究也相對(duì)較少，甚至是阻礙了水合物儲(chǔ)層開(kāi)發(fā)的關(guān)鍵理論基礎(chǔ)之一。

本文主要研究在海洋區(qū)域天然氣水合物開(kāi)采時(shí)，井筒周圍的地層溫度場(chǎng)變化。通過(guò)理論分析建立了低溫條件下含相變的巖石傳熱控制方程。在控制方程中，將比熱和潛熱合并成一個(gè)等效熱容參數(shù)?；贕ibbs-Thomson方程[18-19]，提出了一種利用壓汞數(shù)據(jù)得到孔隙內(nèi)含冰量隨溫度變化的方法。采用所提出的模型，通過(guò)類比水合物開(kāi)發(fā)過(guò)程水合物相變與冰水相變的本質(zhì)，以水為介質(zhì)來(lái)模擬水合物開(kāi)采過(guò)程中井筒溫度升高對(duì)周圍地層的溫度擾動(dòng)，為進(jìn)一步研究開(kāi)采過(guò)程中儲(chǔ)層的穩(wěn)定性提供儲(chǔ)層內(nèi)相變傳熱的理論基礎(chǔ)。

1 理論模型的建立

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

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

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

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

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

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

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

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

圖1 孔喉半徑與溫度壓力間關(guān)系。r(T)：考慮液膜厚度的Gibbs-Thomson方程表示的溫度和冰晶侵入孔徑關(guān)系；r(P)：壓汞法中Washburn方程表示的汞壓與孔徑關(guān)系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).

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

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

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

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

應(yīng)該注意的是，以上關(guān)于冰含量與溫度間關(guān)系的推導(dǎo)是基于在每個(gè)孔隙中凍結(jié)和融化現(xiàn)象都獨(dú)立存在這樣一個(gè)假設(shè)。然而，由于結(jié)冰過(guò)程是在水中形成冰核，在此過(guò)程中產(chǎn)生冰水界面需要消耗表面能。因此，結(jié)冰過(guò)程往往是由外部孔隙向內(nèi)逐步推進(jìn)的，而不是在每個(gè)孔道中獨(dú)立發(fā)生。換句話說(shuō)，結(jié)冰過(guò)程依賴于孔隙的連接方式，一個(gè)孔隙的凍結(jié)依賴于它的入口孔隙半徑大小?？紤]一個(gè)入口孔徑較小的大孔，溫度降低時(shí)大孔不會(huì)先凍結(jié)，而需要等到入口孔喉結(jié)冰后，大孔才會(huì)結(jié)冰。相反，由于液膜的存在，融化過(guò)程不涉及新表面的生成，因此可以在每個(gè)孔隙中獨(dú)立發(fā)生。因此，上述模型僅適用于不受孔隙連通性影響的融化過(guò)程。

2 儲(chǔ)層的溫度場(chǎng)模擬

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

圖2 有限元網(wǎng)格模型Fig.2 Finite element mesh model

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

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

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

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

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

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

圖5 巖石A和B的溫度隨時(shí)間的變化關(guān)系Fig.5 Temperature vs time for sample A and B

3 分析與討論

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

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

4 結(jié)論

針對(duì)水合物開(kāi)發(fā)過(guò)程由于相變現(xiàn)象吸收大量的熱，從而會(huì)造成儲(chǔ)層溫度場(chǎng)的擾動(dòng)這一現(xiàn)象，類比冰水相變進(jìn)行了儲(chǔ)層溫度場(chǎng)的擾動(dòng)模型研究，并研究了不同孔徑分布對(duì)儲(chǔ)層溫度隨時(shí)間變化的規(guī)律的影響。主要得到了以下結(jié)論：

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

(2) 當(dāng)儲(chǔ)層溫度升高時(shí)，其地層溫度變化曲線會(huì)在0 °C附近出現(xiàn)一段緩慢上升的階段。這主要是因?yàn)樵谶@一階段相變劇烈，需要吸收大量熱量，從而使溫度變化并不明顯。

(3) 儲(chǔ)層的溫度演化受巖石本身的孔隙結(jié)構(gòu)影響巨大。當(dāng)孔隙度較大時(shí)，由于骨架的缺失，巖石的整體熱容會(huì)相應(yīng)減少，從而使溫度變化速率升高。同時(shí)，孔徑分布也對(duì)溫度平段的位置和長(zhǎng)度產(chǎn)生。反之，通過(guò)測(cè)量?jī)?chǔ)層的溫度變化曲線，也可以反演得到巖石的孔徑分布。