天然氣水合物降壓開采分解前緣移動數值研究

2020-12-31 03:12:36彭盈鈺蘇正劉麗華金光榮魏雪芹

海洋地質與第四紀地質 2020年6期

彭盈鈺，蘇正，劉麗華，金光榮，魏雪芹

1. 中國科學院天然氣水合物重點實驗室，中國科學院廣州能源研究所，廣州 510640

2. 中國科學院大學，北京 100049

天然氣水合物是由水分子與氣體分子組成的籠狀結構化合物[1]。在自然界中，水合物中的氣體分子主要是甲烷，因此是潛在的能源，同時也成為潛在的環(huán)境災害因素?？碧浇Y果表明，天然氣水合物主要存在于海洋沉積物以及陸地凍土地區(qū)，其全球儲量巨大，約為 2×1016m3[2]。

低溫高壓是天然氣水合物穩(wěn)定存在的必要條件，水合物的開采原理就是破壞水合物相平衡狀態(tài)，使其分解為水和氣體[3]。目前，常規(guī)的開采方法主要有3 種：降壓法[4]、注熱法[5]以及注入抑制劑[6]?，F場試采、實驗以及數值模擬結果表明，降壓法是最具潛力的開采方法，而其他開采方法適合作為輔助手段來提高水合物的開采效率[7]。

開采過程會破壞天然氣水合物的相平衡狀態(tài)，使儲層劃分為不同的區(qū)域：分解區(qū)域、未分解區(qū)域，在兩者之間存在一個特殊的過渡區(qū)域，即分解前緣。以水合物分解前緣為界的分解區(qū)域與未分解區(qū)域之間的流體性質、儲層地質條件以及地層溫度壓力規(guī)律等都不相同[8]。在開采過程中水合物分解前緣移動規(guī)律與水合物開采動態(tài)有著緊密關聯，分解前緣能直接反映了水合物開采特征。另外，由于水合物在儲層中起到了一定的膠結支撐作用，隨著分解前緣的移動，地層穩(wěn)定性降低，可能引發(fā)地質災害[9]。因此，在水合物開采過程中，研究水合物分解前緣移動規(guī)律具有重要指示意義。

目前，對于天然氣水合物分解前緣移動問題的理論研究多基于Stefan 邊界理論，即將水合物分解過程類比冰消融過程，是一個伴隨相變的傳熱過程。Makogon 首次借鑒Stefan 問題，計算得到了降壓分解水合物過程中壓力分布的自相似解[10]。Verigin利用Stefan 移動邊界問題，建立了一維半無限大水合物藏降壓開采模型，模型考慮分解前緣兩側的氣相流動以及分解前緣處氣體質量守恒，根據模型與Stefan 問題的相似性對模型進行線性化自相似求解[11]。Ji 等在Verigin 模型基礎上，建立水合物藏降壓開采數學模型，假設水相靜止，考慮了氣相流動以及溫度的變化，模型認為對流傳熱的作用比熱傳導強，分解區(qū)和水合物區(qū)的能量守恒方程中考慮了熱對流以及節(jié)流和氣體絕熱效應，分解前緣處沒有考慮水合物分解吸熱的作用，將模型方程線性化處理后，自相似求解了分解前緣隨時間的移動[12]。喻西崇等借鑒Ji 等提出的數學模型，利用自相似原理推導出分解前緣移動表達式與溫度、壓力分布表達式[13]。唐良廣，李剛等將水合物分解過程看作移動邊界問題，建立了水合物層溫度分布的一維傳熱模型，模型考慮了分解區(qū)和水合物區(qū)的熱傳導以及分解前緣處的能量守恒，根據自相似求解得到不同時刻水合物藏溫度分布以及分解前緣的位置[14]。張旭輝等建立二維熱傳導模型，研究水合物儲層有熱水管垂直穿過時水合物最大分解范圍，采用分離變量法對模型進行求解，結果表明分解前緣的最大移動距離隨溫度的增大而增大[15]。劉樂樂建立水平一維降壓-加熱數學模型，將模型有限差分離散后數值計算得到分解相變陣面的位置與時間的平方根呈正比[16]。李明川等建立注熱移動界面的三相一維傳質模型，數值差分計算得到分解前緣移動速度前期較高，后逐漸降低[17]。另外，隨著數值模擬器的成熟，Long 和 Tjok 利用 HydrateResSim 模擬水合物藏降壓開采分解前緣移動，結果表明分解前緣的平均速度隨絕對滲透率的增大而增大[18]。鄭如意等數值模擬研究模型邊界條件、滲透率、初始水合物飽和度、總熱導率、井筒加熱溫度和井底壓力等對水合物分解前緣移動的影響，發(fā)現分解前緣移動速度與滲透率、總熱導率、井筒加熱溫度和邊界供熱成正比，相反，增大水合物初始飽和度和井底壓力會降低分解前緣移動速度。此外，分解前緣移動規(guī)律也會隨著參數的變化而變化[19]。

總的來說，前人采用數學建模以及數值模擬等多種方法對水合物分解前緣的移動規(guī)律進行了大量研究。與模型數值解相比，解析計算模型求解方便，計算過程也有助于深入了解某些物理變化的重要性。但大部分解析計算依據自相似原理，以分解前緣的移動與時間的平方根呈線性關系為前提假設，由此得到的分解前緣移動只在平均速度上存在差異，而分解前緣位置隨時間變化規(guī)律已經確定，這將影響分解前緣移動規(guī)律探究。此外，利用現有數值模擬器探究水合物分解前緣移動規(guī)律時，要綜合網格壓力、相平衡壓力以及水合物飽和度來判斷分解前緣移動位置，這個過程相對復雜繁瑣，并且單一儲層模型的探究不具有普適性。

據此，本文建立了水合物降壓分解一維三相數學模型，區(qū)別于傳統(tǒng)自相似求解假設前提（即分解前緣的移動與時間的平方根呈線性關系），利用量級分析與偏微分方程無維化轉化方法，解析計算探究了分解前緣隨時間的移動規(guī)律，并由分解前緣移動計算相關產氣量，對水合物開采動態(tài)進行了簡單快速評估。

1 概念模型

模型為Class3 水合物藏。儲層頂底板滲透率低，壓力傳遞較慢，將頂部層和底部層認為是定溫定壓邊界，允許發(fā)生熱量和流水交換。

水合物分解是一個吸熱的相變過程，會導致地層溫度降低，此時，周圍環(huán)境就會向水合物層傳遞熱量。這種熱補償除了維持水合物分解所需熱量外，還可以彌補吸熱反應導致的地層溫度下降。對于具有一定厚度的水合物層，可忽略頂底層圍巖傳入熱量的影響，水合物分解相變前緣所吸收的熱量主要來自于單位體積儲層內能和未分解區(qū)熱傳導[20]。

隨著時間的推移，水合物分解由近井區(qū)域向外擴散。假設水合物分解不是發(fā)生在整個儲層內，而是發(fā)生在一定的狹窄區(qū)域內，可以將該區(qū)域視為一個表面，即所謂的水合物分解前緣[12]。它將儲層分為兩個部分，即水合物分解區(qū)與水合物未分解區(qū)。水合物分解區(qū)域自由氣體在壓力梯度的驅動下向井內流動，而分解前緣則向相反的方向移動。

根據前述假設，水合物分解過程可簡化為分解前緣隨時間向外移動的過程，而分解前緣就是天然氣水合物發(fā)生分解的臨界面，其厚度忽略不計，分解前緣處的地層壓力即為該地層溫度下的水合物相平衡壓力[21]。在水合物層發(fā)生降壓分解時，儲層尺度內可采用水合物平衡分解模型，分解前緣界面（S（t）位置）把水合物層分為兩個區(qū)域（圖 1）：已分解氣水區(qū)（r0＜x＜S（t））與未分解水合物區(qū)（S（t）＜x＜∞）。

根據概念模型的簡化，本文進一步假設為：

（1）模型考慮三相（水合物、甲烷、水）兩組分（甲烷、水），不考慮甲烷和水合物的溶解；

（2）由于氣體與水之間的壓力差不大，忽略毛細管壓力的影響；

（3）在水合物開采中，擴散作用貢獻小于對流作用，忽略氣水擴散對水合物分解的影響；

（4）水合物分解所需能量主要包括分解前緣所在區(qū)域的單位體積儲層內能以及未分解區(qū)域傳導熱。

圖1 天然氣水合物開采分解前緣遷移概念模型Fig.1 Conceptual model of migration of decomposition front of gas hydrate

2 數學模型

在數學模型中，將水合物儲層看作x 方向的一維流體場。描述多孔介質中水合物分解過程的主要方程包括水合物分解區(qū)域和分解前緣處的質量守恒方程、能量守恒方程。其中，分解前緣移動相關方程參考Stefan 模型對冰水自由邊界的描述[22]，其他方程則類似于Tsypkin[23-24]、Ahmadi[25]等提出的模型方程。

2.1 模型主要方程

水合物分解區(qū)域質量守恒方程：

當分解區(qū)與未分解區(qū)的壓力梯度都比較小時，水合物分解區(qū)域與未分解區(qū)域的能量守恒方程簡化為（體現水合物分解區(qū)域與未分解區(qū)域溫度變化關系）：

式中：下標1，2 表示分解區(qū)與未分解區(qū)，λ 是導熱系數，c是比熱容，為拉普拉斯算子。

水合物分解前緣處質量守恒方程為：

水合物分解前緣處能量守恒方程為（體現水合物分解前緣上的溫度變化）：

2.2 模型輔助方程

分解前緣壓力為水合物分解為氣水的相平衡壓力[26]：

計算熱傳導系數方程：

2.3 模型初邊值條件

模型初始條件與邊界條件如表1 所示。

2.4 求解結果

我們將分解區(qū)氣相質量守恒方程（1）與分解前緣氣相質量守恒方程（3）作量級分析（具體步驟見附錄A）。量級分析的思想是，如果方程是基于無量綱和歸一化變量的表現形式，方程不同項的系數能用來度量這些項的重要程度[27]。無量綱化后的方程（1）中第二項系數遠大于第一項系數1，可忽略方程（1）中第一項對方程的影響。另外，由于氣水相運動黏度相差100 倍，在相同壓力梯度下，水合物分解后，達到傳輸平衡，水相飽和度認為是不隨時間變化。因此，在這個模型中水合物分解區(qū)流體流動簡化為擬定常流動，得到分解區(qū)壓力傳導關系：

表1 初始條件與邊界條件Table 1 Initial conditions and boundary conditions

根據分解區(qū)壓力傳導方程（7）與分解前緣質量守恒方程（3）得到分解前緣隨時間移動規(guī)律：

由式（8）可知分解前緣的移動與流體相滲透率，水合物分解相平衡壓力與井底壓力之間的差值有關，與時間的平方根呈線性關系。

將分解區(qū)與未分解區(qū)能量守恒方程作無量綱轉換后代入積分得到（具體步驟見附錄A）：

將式（9），（8）代入式（4）并與式（5）聯立得到一個用于求解分解前緣相平衡壓的超越方程組：

單位橫截面積生產井產氣速率為：

由式（11）可以看到，產氣速率與分解前移動距離成反比，說明隨著天然氣水合物的分解，產氣量逐漸減少。

當分解前緣移動到S（t）時，單位橫截面積生產井總產氣體積為：

實際生產井總產氣體積為：

式中：h為井射孔有效長度為井孔半徑。

3 討論

3.1 模型驗證

Yousif 通過水合物砂巖樣品降壓分解實驗，探究了分解前緣移動現象[28]。以Yousif 實驗為依據，將實驗中所設置的相關參數代入模型中，并通過對超越方程（10）以及式（5），（8）聯立求解，得到分解前緣移動位置，實驗結果與模型計算結果進行對比，如圖2 所示，實驗數據與模型計算結果近似，從而驗證了模型結果的可靠性。

3.2 示例分析

以南海神狐海域天然氣水合物藏[29-31]實際參數為例，模型計算所需的基本參數如表2 所示。不考慮水合物分解過程中冰的生成對模擬結果的影響，井底壓力為 3 MPa。

圖2 實驗與模型下的分解前緣移動Fig.2 Model match of experimental results for location of hydrate dissociation front

表2 相關物性參數Table 2 The correlated parameters used for calculation

通過表2 所示模型參數，進行水合物儲層降壓開采模型計算，得到200 d 內分解前緣隨時間變化規(guī)律見圖3。在天然氣水合物降壓開采過程中，當地層壓力低于水合物相平衡壓力時，水合物開始分解，并出現水合物分解前緣。在開采60、120、200 d后，模型計算分解前緣隨時間移動距離分別約為33.35、47.17、60.90 m。從圖 3a 中，可以看到隨著開采時間的推移，水合物分解前緣移動曲線斜率變小，說明分解前緣移動速度降低。

根據水合物分解前緣移動速率變化（圖3b），看出水合物分解前緣移動速率在生產前期達到最大值，隨開采時間推移，移動速率變慢最后將趨于平穩(wěn)。出現這一現象的原因在于，當井底降壓開始時，水合物平衡狀態(tài)被打破，此時壓降開始傳遞，井壓與地層壓力差作為水合物分解的主要驅動力。隨著在地層中壓力傳遞以及水合物分解過程能量消耗，導致地層溫度下降，水合物相平衡壓力與井壓之間的壓差減小，分解過程隨之變慢；在開采后期，儲層能量不足，水合物分解主要依靠儲層熱量的傳導，分解前緣移動速率保持較低的平穩(wěn)狀態(tài)，在這種情況下分解主要受儲層熱物理性質的影響。此時，我們應該考慮儲層注熱等技術進一步促進水合物分解。

圖3 水合物開采特征a. 分解前緣隨時間移動規(guī)律，b. 分解前緣移動速率隨時間變化，c. 產氣速率隨時間變化，d. 總產氣量隨時間變化Fig.3 Characteristics of hydrate productiona. The moving of decomposition front with time，b. the velocity of decomposition front varies with time，c. the rate of gas production varies with time， d. the volume of gas production varies with time

圖3c 給出了甲烷產氣速率（單位高度和單位寬度的產量）隨時間變化，我們可以看到，生產井（單位高度和單位寬度）在200 d 內,最大產氣速率約為250 m3/d，后期產氣速率接近 20 m3/d。水合物分解氣體體積變化與分解前緣移動有關，即井口產氣速率受到水合物分解前緣移動的影響，因此，產氣速率隨時間變化趨勢與分解前緣移動速率一致。同時，如圖3d 所示，在200 d 內甲烷總產氣量（單位高度和單位寬度的產量）為 18 000 m3。

圖4 模型計算總產氣量與試采結果對比Fig.4 The model calculated results vs the test result

為了進一步分析模型計算結果，將模型結果與2017 年南海神狐海域第一次水合物試采情況進行比對。根據中國地質調查局的報告，連續(xù)產氣60 d后試采結束，累產氣量3 ×105m3，平均日產氣 5 000 m3[32]。另外，天然氣水合物儲層厚度為40 m，假設井射孔貫穿整個水合物層，井孔半徑為0.15 m。我們將60 d試采數據進行曲線擬合，得到圖4 所示模型計算總產氣量與試采總產氣量變化。從圖4 可以看到，在開采10 d 內的兩者產氣量比較吻合，隨著水合物開采過程的進一步推進，模型計算總產氣值高于與實際開采結果，兩者差異增大，但整體相對誤差在可接受范圍內。模型計算結果和實際試采結果之間存在較大差異的原因主要在于3 個方面，首先本文模型中認為井射孔貫穿整個水合物層，水合物降壓分解出的甲烷氣體能快速有效的向井口方向移動，減小了孔隙中氣體積聚所導致的地層壓力增大的影響，從而有利于分解甲烷氣的產出、后續(xù)壓降傳遞和水合物的進一步分解，同時，模型忽略了實際開采中氣體的溶解與擴散等作用的消耗；其次本模型沒有考慮水合物二次生成、冰的形成等對水合物分解過程以及甲烷氣流動的影響，在實際開采中，當儲層局部溫度低于冰點以下時，會有冰的生成以及水合物的再次形成，固相物質的出現會降低地層孔隙度和地層絕對滲透率，從而影響壓降的傳遞，水合物分解減慢，并阻礙流體向井口方向流動，導致產氣量減小，最后，實際的天然氣水合物藏是一個三維空間，水合物分解過程是發(fā)生在三維空間內的物理化學變化，本文僅是一維簡單模型，不可避免與實際情況存在差距。因此，本文模型計算所得天然氣水合物降壓開采總產氣量高于實際開采值，是對總產氣量變化的樂觀預測，能簡單快速地為實際開采提供大方向的參考。

3.3 敏感性分析

為了分析水合物儲層相關參數對水合物分解前緣移動距離的影響，本文采用單次單因子敏感性分析方法，除了變量參數外，其他參數均保持在前述的參考數值[33]。

3.3.1 儲層初始溫度

地層初始溫度是影響天然氣水合物開采的重要儲層參數[34]。探究不同地層溫度下單井降壓開采天然氣水合物分解前緣移動變化規(guī)律。為了不改變儲層初始壓力，將溫度變化控制在滿足水合物相平衡條件內。

在地層初始壓力與生產井壓不變的情況下，不同地層溫度下（286.15～288.15 K）水合物分解前緣移動變化如圖5 所示。可以看到，儲層溫度的改變顯著影響水合物分解前緣的移動，當儲層溫度增大時，分解前緣移動距離增大，儲層溫度變化1 K 時，200 d 分解前緣移動距離在參考情況（287.15 K）基礎上變化33%。出現這一現象的原因是，天然氣水合物分解是一個吸熱的過程，儲層初始溫度越高，儲層所能提供水合物分解的能量越多，有利于水合物分解；此外，隨著儲層溫度的升高，相應的水合物相平衡壓力增大，水合物分解速度更快。因此，當水合物儲層溫度很低時，多考慮提高儲層溫度作為輔助手段來促進水合物分解，提高產氣量。

圖5 不同儲層初始溫度下水合物分解前緣移動距離Fig.5 The location of hydrate dissociation front at different initial reservoir temperatures

3.3.2 儲層絕對滲透率

儲層滲透率是反映流體運移能力的重要水力學參數[35]。最新的調查數據表明，南海儲層類型主要為黏土質粉砂—低滲透粉砂，地層絕對滲透率平均為2～5 md[36]。據此，敏感性分析的地層滲透率范圍為 2～5 md。

對于不同地層滲透率，圖6 顯示了分解前緣移動隨時間變化。正如預期的那樣，隨著地層滲透率的減小，分解前緣移動速率變慢，移動距離減小。地層滲透率變化對水合物分解有明顯的影響，較高的地層滲透率可以提高流體運移速率，有利于地層整體壓降傳遞，為天然氣水合物分解提供更大的驅動力，水合物分解前緣移動距離隨之增大，產氣量隨之增加。當儲層滲透率提高 1 md（以 2.9 md 情況為基準），分解前緣移動距離較參考情況增大17%。

3.3.3 儲層孔隙度

另一個重要的地層參數是孔隙度，而滲透率又是與孔隙度有關函數，在研究儲層孔隙度變化對水合物分解前緣移動影響時，根據Kozeny-Carman 方程[37]（孔隙度的三次方與滲透率之間成正比），在改變地層孔隙度的同時，地層滲透率也隨之改變。

隨著地層孔隙度的增大，分解前緣移動速率減小，移動距離也隨之減?。▓D7）。在孔隙度為0.38時，此時水合物相平衡壓力為3.05 MPa，井口和井口前端之間的壓差很小，這意味著巖石的熱量不足以用于水合物分解，此時只能通過水合物未分解區(qū)域的熱量流入來獲得額外的能量。在這種情況下，分解前緣移動由儲層熱物理參數決定。這是因為地層孔隙度很大時，含天然氣水合物的沉積物單位體積潛熱也較大，但由于天然氣水合物含量高，比熱低，熱導率低，使得地層的熱導率變小，因此,地層溫度下降更快，從而不利于水合物分解。

圖6 不同儲層絕對滲透率下水合物分解前緣移動距離Fig.6 The location of hydrate dissociation front at different absolute permeability

圖7 不同儲層孔隙度下水合物分解前緣移動距離Fig.7 The location of hydrate dissociation front at different reservoir porosity

4 結論

（1）通過參數量級分析，將天然氣水合物分解滲流場簡化為擬定常場，積分求解得到分解前緣隨時間移動規(guī)律，同時將求解溫度場變化方程進行無維化轉換，推導出溫度分布方程，并根據Yousif 對含甲烷水合物砂巖樣品降壓分解實驗進一步驗證了模型結果的可靠性。

（2）通過南海神狐海域相關參數的實例分析，發(fā)現多孔介質中水合物分解前緣移動速率隨時間減小，移動距離與時間的開平方呈線性關系；井口氣體產量隨時間減小，最后趨于穩(wěn)定。另外，對比了南海水合物試采結果與模型計算的總產氣量，發(fā)現模型計算值高于實際開采結果，且兩者相對誤差在25%范圍內。

分析誤差出現的原因，主要是建立的模型忽略了氣體的溶解與擴散等作用，且沒有考慮水合物二次生成、冰的形成等對水合物分解過程以及甲烷氣流動的影響。因此，本文通過研究水合物分解前緣移動規(guī)律對產氣速率、產氣量進行了較為樂觀的預測，為水合物開采潛力評價提供了一種新的簡單快速的計算方法，同時對開采動態(tài)評估給出大方向的參考。

（3）通過對地層初始溫度、絕對滲透率以及孔隙度敏感性分析發(fā)現，地層初始溫度和滲透率與水合物分解前緣移動距離之間成正相關關系，初始地層溫度對水合物分解過程影響顯著，另外，地層孔隙度越大，分解前緣移動速率反而降低，移動距離減小，井口與分解前緣壓差減小，此時分解前緣移動速率由儲層熱物理參數決定。

附錄 A 模型求解

將分解區(qū)氣相質量守恒方程（1）中基本物理量表示為特征值與無量綱數的形式：