林星宇， 李 凱， 孔 亮

(青島理工大學(xué)理學(xué)院， 青島 266033)

天然氣水合物是一種新型的清潔能源并廣泛分布于深海的沉積土體之中。中國(guó)已探明在南海海域存在著儲(chǔ)量巨大的天然氣水合物礦藏。開(kāi)采并利用天然氣水合物資源，對(duì)滿足中國(guó)日益增長(zhǎng)的能源消費(fèi)以及解決能源安全問(wèn)題具有重要的戰(zhàn)略意義[1]。天然氣水合物的賦存條件極為苛刻，其穩(wěn)定賦存于低溫高壓的環(huán)境下，外界輕微的擾動(dòng)便會(huì)打破其原有的平衡條件從而導(dǎo)致其急劇分解。若不加以控制，不但將破壞水合物礦藏，誘發(fā)海底滑坡等災(zāi)害，還會(huì)加劇溫室效應(yīng)。因此，安全有效地開(kāi)采天然氣水合物是一項(xiàng)極具挑戰(zhàn)的任務(wù)。

目前，學(xué)者們已經(jīng)提出多種開(kāi)采方案，主要開(kāi)采方法有降壓法、注熱法、化學(xué)試劑法和置換法[2-3]。其中，降壓開(kāi)采是一種相對(duì)經(jīng)濟(jì)、簡(jiǎn)單、高效的方式。天然氣水合物的開(kāi)采機(jī)理主要是人為打破水合物穩(wěn)定賦存的溫壓條件，將固態(tài)的水合物轉(zhuǎn)化為液態(tài)與氣態(tài)，再在地面上完成收集工作。

中外學(xué)者針對(duì)水合物的開(kāi)采模型開(kāi)展了大量研究。Kimoto等[4]對(duì)海底天然氣水合物的離解進(jìn)行了數(shù)值分析，考慮了固體到流體的相變，水和氣體的流動(dòng)、傳熱和地面變形；劉樂(lè)樂(lè)等[5]提出了一維水合物降壓開(kāi)采數(shù)學(xué)模型，該模型考慮了滲流場(chǎng)、變形場(chǎng)以及溫度場(chǎng)的多場(chǎng)耦合效應(yīng)，并對(duì)多孔介質(zhì)的變形進(jìn)行了敏感性分析；Sun等[6]提出了包含溫度場(chǎng)、滲流場(chǎng)、應(yīng)力場(chǎng)、化學(xué)場(chǎng)的耦合模型，且其建立的模型考慮了變形-滲流的雙向耦合過(guò)程；程遠(yuǎn)方等[7]建立了考慮氣-水兩相流動(dòng)的非等溫的開(kāi)采模型，通過(guò)對(duì)水合物礦藏的實(shí)例研究，分析了影響儲(chǔ)層應(yīng)力狀態(tài)的主要因素；胡曉慶等[8]以實(shí)際水合物儲(chǔ)藏為例，對(duì)水合物降壓開(kāi)采的產(chǎn)能進(jìn)行了模擬，分析了流固耦合作用對(duì)降壓開(kāi)采的動(dòng)態(tài)影響機(jī)制。盡管以上模型考慮了水合物降壓開(kāi)采時(shí)所涉及的主要因素，但對(duì)于水合物儲(chǔ)層的分解區(qū)域卻沒(méi)有明確給出。

ABAQUS是一款功能強(qiáng)大的有限元軟件，擅長(zhǎng)分析復(fù)雜的固體力學(xué)、結(jié)構(gòu)力學(xué)系統(tǒng)，特別是可以模擬高度非線性等問(wèn)題。同時(shí)，它提供了用戶子程序接口來(lái)自定義程序從而滿足各類用戶需求。因此，基于水合物降壓開(kāi)采基本方程，自主編寫(xiě)USDFLD用戶子程序，從而確定水合物儲(chǔ)層的分解區(qū)域，揭示水合物分解過(guò)程儲(chǔ)層物性參數(shù)變化規(guī)律，即在ABAQUS中實(shí)現(xiàn)降壓開(kāi)采條件下深海能源土斜坡的有效應(yīng)力以及變形場(chǎng)計(jì)算。

1 水合物降壓開(kāi)采基本方程

1.1 水合物分解特征方程

1.1.1 水合物相平衡曲線

天然氣水合物處于低溫高壓的環(huán)境當(dāng)中，溫度與壓力是影響水合物能否穩(wěn)定存在的關(guān)鍵。低溫高壓的賦存條件十分苛刻，輕微的干擾便可引起水合物分解。當(dāng)溫壓條件處于水合物相平衡曲線之下時(shí)，水合物便開(kāi)始分解。Ahmadi等[9]建立了溫度與壓力之間的相平衡方程，其表達(dá)式為

lgPeq=0.034 2(T-273.15)+0.000 5(T-273.15)2+6.480 4

(1)

式(1)中：Peq為溫度對(duì)應(yīng)的相平衡壓力，Pa；T為溫度，K。

1.1.2 水合物分解動(dòng)力學(xué)方程

水合物分解化學(xué)方程式為

(2)

式(2)中：n=5.75。

(3)

式(3)中：NH為單位體積水合物的物質(zhì)的量，mol；P為降壓開(kāi)采后的孔壓，Pa；N0為初始狀態(tài)下水合物的物質(zhì)的量，mol。

由水合物分解化學(xué)方程式可知：

(4)

(5)

而單位體積水合物分解消耗的水合物、以及產(chǎn)生的甲烷氣體與水的質(zhì)量速率分別為

(6)

(7)

(8)

式中：MH、MG、MW分別表示水合物、甲烷、水的摩爾質(zhì)量。

1.1.3 水合物分解超孔壓方程

天然氣水合物分解后，會(huì)產(chǎn)生大量的水和氣體。而當(dāng)水合物沉積層滲透率較低且排列緊密時(shí)，分解產(chǎn)生的壓力無(wú)法快速消散，從而導(dǎo)致超孔壓激增，進(jìn)而導(dǎo)致沉積土體的有效應(yīng)力降低，發(fā)生海底滑坡。因此需要對(duì)水合物分解之后所產(chǎn)生的超孔壓進(jìn)行計(jì)算。Nixon等[11]提出了超孔壓計(jì)算模型，其表達(dá)式為

(9)

式(9)中：Δσ′為水合物分解導(dǎo)致的有效應(yīng)力的改變量，Pa；Δu為孔壓的變化量，Pa；M為土體的側(cè)限壓縮模量，Pa；nsoil為土體的孔隙率；S為水合物的含水飽和度。

1.2 滲流場(chǎng)方程

1.2.1 質(zhì)量守恒方程

水相、氣相、水合物相和土骨架相的質(zhì)量守恒方程如式(10)～式(13)所示：

(10)

(11)

(12)

(13)

式中：下標(biāo)w、g、h、s分別表示水、氣、水合物、土骨架；ξw、ξg、ξh、ξs分別表示水、氣、水合物、土骨架的體積分?jǐn)?shù)，ρw、ρg、ρh、ρs分別表示水、氣、水合物、土骨架的密度，Vw、Vg、Vs分別代表水、氣、土骨架的達(dá)西速度；t為時(shí)間；xi為坐標(biāo)。

1.2.2 達(dá)西定律

達(dá)西定律又稱為線性滲流定律，其表達(dá)式為

(14)

式(14)中：Pex為水合物分解后產(chǎn)生的平均孔壓；μ為動(dòng)力黏滯系數(shù)；K為含水合物多孔介質(zhì)的滲透率；取值參考Masuda實(shí)驗(yàn)參數(shù)[12-13]。

1.3 變形場(chǎng)方程

飽和土是由固體土骨架與孔隙水組成的兩相體。當(dāng)受到外力后，總應(yīng)力可通過(guò)兩部分來(lái)承擔(dān)。其中由孔隙水來(lái)承擔(dān)，通過(guò)連通的孔隙水來(lái)傳遞的壓力稱為孔隙水壓力。而通過(guò)土骨架承擔(dān)，通過(guò)土顆粒之間的接觸來(lái)傳遞的壓力稱為有效應(yīng)力。其可表示為

σ=σ′+u

(15)

式(15)中：σ為土體的總應(yīng)力，Pa；σ′為土體的有效應(yīng)力，Pa；u為孔隙水壓力，Pa。

基于有效應(yīng)力原理，土骨架的平衡方程可以表示為

σij,j-(αδijPex),j+Fi=0

(16)

式(16)中：σij為總應(yīng)力張量；α為Biot系數(shù)；δij為克羅內(nèi)克張量；Fi為體力。(αδijPex),j反映滲流場(chǎng)對(duì)多孔介質(zhì)土骨架變形場(chǎng)的影響。

根據(jù)小變形理論，土的幾何方程寫(xiě)為

(17)

式(17)中：ui,j和uj,i為土骨架的位移；下標(biāo)i、j為約定符號(hào)。

將土體的彈性階段與塑性階段的本構(gòu)關(guān)系寫(xiě)成統(tǒng)一的形式，得

σ′ij=Dijklεkl

(18)

式(18)中：σ′ij為有效應(yīng)力張量；Dijkl為彈塑性矩陣；εkl為應(yīng)變張量；下標(biāo)i、j、k、l為約定符號(hào)。本構(gòu)模型選用修正劍橋模型。

1.4 有限元強(qiáng)度折減法

強(qiáng)度折減法提出了抗剪強(qiáng)度折減系數(shù)的概念，其原理可以解釋為在有限元計(jì)算中使邊坡土體抗剪強(qiáng)度逐漸降低至其破壞為止，此時(shí)程序根據(jù)彈塑性計(jì)算結(jié)果得到滑裂面，同時(shí)還可以得到邊坡的強(qiáng)度折減系數(shù)，此時(shí)的折減系數(shù)即為邊坡的安全系數(shù)[14]。

折減后的抗剪強(qiáng)度參數(shù)可表達(dá)為

c′=c/Fr

(19)

φ′=arctan(tanφ/Fr)

(20)

式中：c、φ分別為土體黏聚力和摩擦角；c′、φ′分別為土體折減后的黏聚力和摩擦角；Fr為折減系數(shù)，即邊坡的安全系數(shù)。

2 USDFLD用戶子程序的二次開(kāi)發(fā)與水合物儲(chǔ)層分解區(qū)域確定

對(duì)于降壓開(kāi)采水合物儲(chǔ)層分解區(qū)域的確定，需要以水合物降壓開(kāi)采基本方程為依據(jù)，通過(guò)將其編寫(xiě)到自定義的子程序中來(lái)完成。ABAQUS不僅能提供標(biāo)準(zhǔn)的有限元分析程序，同時(shí)還具有良好的開(kāi)放性，可利用它提供的用戶子程序接口來(lái)生成自定義的程序從而滿足用戶的需求，用戶可根據(jù)實(shí)際需求靈活定義荷載的分布模式、材料的本構(gòu)關(guān)系等。

2.1 USDFLD用戶子程序介紹

USDFLD是ABAQUS中常用來(lái)定義場(chǎng)變量的子程序，它允許在單元的積分點(diǎn)上定義場(chǎng)變量作為積分點(diǎn)變量的函數(shù)，也可以被定義為積分點(diǎn)上已有的計(jì)算結(jié)果數(shù)據(jù)的函數(shù)。同時(shí)也可以定義場(chǎng)變量相關(guān)的材料屬性。其編寫(xiě)規(guī)則簡(jiǎn)單、便捷，能夠滿足編寫(xiě)者對(duì)材料的各種自定義。USDFLD子程序接口如圖1所示。

圖1 USDFLD子程序接口Fig.1 The USDFLD subroutine interface

2.2 USDFLD用戶子程序場(chǎng)變量的定義

ABAQUS通常認(rèn)為場(chǎng)變量fi是節(jié)點(diǎn)數(shù)據(jù)，而當(dāng)ABAQUS開(kāi)始計(jì)算單元的應(yīng)力與剛度時(shí)，它會(huì)將節(jié)點(diǎn)數(shù)據(jù)fi插值到單元的積分點(diǎn)。而在利用USDFLD子程序計(jì)算單元的材料屬性前，在USDFLD子程序中被定義的這些數(shù)值會(huì)取代插值fi。因此，為了訪問(wèn)之前所定義的fi，需要在USDFLD中把它們存儲(chǔ)為相關(guān)解變量。

基于FORTRAN語(yǔ)言在子程序中定義場(chǎng)變量P、T分別為水合物儲(chǔ)層的孔壓與溫度。將水合物相平衡方程與水合物分解動(dòng)力學(xué)方程編入U(xiǎn)SDFLD子程序，通過(guò)設(shè)置條件語(yǔ)句對(duì)降壓后水合物儲(chǔ)層的孔壓進(jìn)行判斷。當(dāng)孔壓大小處于相平衡曲線之下時(shí)，水合物儲(chǔ)層處于不穩(wěn)定狀態(tài)便會(huì)開(kāi)始分解。

2.3 降壓開(kāi)采能源土分解程度的判斷

圖2為一維的水合物儲(chǔ)層，長(zhǎng)度為12 m，高度為1 m，儲(chǔ)層溫度為285 K，水合物飽和度為0.40。沿其左邊界設(shè)置降壓條件來(lái)模擬開(kāi)采井的降壓效果，降壓時(shí)間為48 h，降壓為13～7 MPa。在子程序接口中編寫(xiě)USDFLD用戶子程序，定義溫度常量與孔壓變量以及水合物相平衡方程。

圖2 水合物層降壓飽和度分布Fig.2 Distribution contour of saturation of hydrate layer under the condition of depressurization

當(dāng)溫壓條件處于相平衡曲線的不穩(wěn)定區(qū)域時(shí)，水合物開(kāi)始分解。水合物相平衡曲線如圖3所示，由子程序判斷得到的水合物儲(chǔ)層飽和度分布云圖如圖2所示。

圖3 水合物相平衡曲線Fig.3 Curve of hydrate phase equilibrium

由圖2可知，當(dāng)降壓開(kāi)采之后，靠近開(kāi)采井的水合物層分解完全。而遠(yuǎn)離開(kāi)采井的部分由于降壓效果的減弱使其孔壓再次處于相平衡曲線的穩(wěn)定區(qū)域，導(dǎo)致分解程度逐漸變小，水合物飽和度回歸到初始值。因此，可得到降壓開(kāi)采水合物儲(chǔ)層的分解區(qū)域。

2.4 相關(guān)物性參數(shù)的變化

水合物層相關(guān)的物性參數(shù)，如滲透率、彈性模量是水合物飽和度的函數(shù)。當(dāng)水合物層降壓分解之后，水合物飽和度發(fā)生了變化，繼而影響了滲透率以及彈性模量的大小。

采用含水合物多孔介質(zhì)滲透率模型[13]表示水合物的滲透率，其表達(dá)式為

K=K0[1-SH]N

(21)

式(21)中：K0為無(wú)水合物時(shí)多孔介質(zhì)的滲透率；SH為水合物飽和度；N為滲透率指數(shù)。

采用線性模型[5]反映水合物飽和度與彈性模量的關(guān)系，其表達(dá)式為

(22)

式(22)中：E為水合物沉積物的彈性模量；E0為無(wú)水合物時(shí)沉積物的彈性模量；E′為飽和度為SH0時(shí)水合物沉積物所對(duì)應(yīng)的彈性模量；SH0代表某一時(shí)刻的水合物飽和度。

在USDFLD用戶子程序中，設(shè)置相應(yīng)的狀態(tài)變量，編寫(xiě)滲透率以及彈性模量的函數(shù)方程。水合物層分解之后，其滲透率及彈性模量的分布分別如圖4、圖5所示。

由圖2、圖4、圖5可知，當(dāng)水合物分解后，產(chǎn)生甲烷氣體與水，水合物飽和度降低，土體的彈性模量也隨之降低，而能源土的滲透率增大。

圖4 水合物層滲透率分布Fig.4 Distribution contour of permeability of hydrate layer

圖5 水合物層彈性模量分布Fig.5 Distribution contour of elastic modulus of hydrate layer

3 實(shí)例分析

3.1 工況介紹

3.1.1 含水合物層斜坡的幾何模型及參數(shù)選取

選用南海北部的某個(gè)海底斜坡，幾何模型如圖6所示，為了模擬實(shí)際的海底斜坡工況，將斜坡劃分為上覆土層、水合物層、水合物周邊土層、下伏土層四部分。海底斜坡水平尺寸為4 500 m，高度為1 300 m，坡度為4°，上覆土層的厚度為60 m，水合物層長(zhǎng)度為2 506 m，厚度為30 m，埋深為150 m。通過(guò)USDFLD用戶子程序判斷得到水合物降壓開(kāi)采分解區(qū)長(zhǎng)度為152 m。下伏土層的厚度為860 m。溫度為278 K，初始孔隙率為0.35，水合物飽和度為0.25。各層土體的土性參數(shù)根據(jù)土工實(shí)驗(yàn)以及中外近年發(fā)表的典型數(shù)據(jù)[15-18]綜合確定，具體如表1所示。

圖6 海底斜坡幾何模型Fig.6 Geometric model of submarine slope

表1 土體力學(xué)參數(shù)

3.1.2 含水合物層海底斜坡的邊界條件

斜坡左右邊界為水平位移約束，底部為固定約束，斜坡頂面施加靜水壓力為13 MPa。

在水合物儲(chǔ)層的左邊界設(shè)置降壓井，施加降壓邊界條件，由初始的13 MPa降低至4 MPa。水合物完全分解所產(chǎn)生的超孔壓通過(guò)Grozic超孔壓模型計(jì)算為13 MPa。

3.2 結(jié)果分析

通過(guò)在ABAQUS中調(diào)用USDFLD子程序，判斷得到水合物儲(chǔ)層的分解區(qū)域，之后對(duì)海底斜坡施加應(yīng)力與位移邊界條件，得到含水合物層斜坡的有效應(yīng)力、塑性應(yīng)變分量、位移的分布云圖如圖7～圖9所示。

通過(guò)圖7～圖9可以看出，當(dāng)降壓開(kāi)采水合物時(shí)，水合物快速分解，在分解區(qū)會(huì)產(chǎn)生很大的超孔壓，由于地質(zhì)構(gòu)造，超孔壓一時(shí)難以消散。由有效應(yīng)力原理可知，土體的有效應(yīng)力迅速減小。在此基礎(chǔ)上，對(duì)海底斜坡進(jìn)行有限元強(qiáng)度折減法，以此來(lái)判斷海底斜坡的安全系數(shù)。

當(dāng)尚未進(jìn)行降壓開(kāi)采時(shí)，對(duì)含水合物層海底斜坡進(jìn)行強(qiáng)度折減，得知其安全系數(shù)為2.903，其塑性區(qū)貫通云圖如圖10(a)所示。

當(dāng)實(shí)施單井降壓開(kāi)采，即水合物完全分解之后，再對(duì)整個(gè)海底斜坡進(jìn)行強(qiáng)度折減，得到其安全系數(shù)為2.652，而塑性區(qū)貫通云圖如圖10(b)所示。

對(duì)以上結(jié)果(圖7～圖10)進(jìn)行分析可知，當(dāng)對(duì)含水合物層海底斜坡進(jìn)行單井降壓開(kāi)采時(shí)，水合物層的分解區(qū)相比水合物層分布長(zhǎng)度較小。水合物分解會(huì)導(dǎo)致水合物層分解區(qū)局部產(chǎn)生較大的超孔壓，從而對(duì)分解區(qū)局部的有效應(yīng)力影響較大。但從整個(gè)海底斜坡來(lái)看，單井降壓開(kāi)采對(duì)海底斜坡的穩(wěn)定性沒(méi)有顯著影響，其開(kāi)采前與開(kāi)采后的安全系數(shù)變化很小。

圖7 斜坡豎向有效應(yīng)力分布Fig.7 Distribution contour of vertical effective stress of slope

圖8 斜坡塑性應(yīng)變分量Fig.8 Contour of plastic strain component of slope

圖9 斜坡位移分布Fig.9 Distribution contour of displacement of slope

圖10 水合物斜坡塑性區(qū)貫通云圖Fig.10 Contour of plastic penetration area

4 結(jié)論

(1)歸納了水合物降壓開(kāi)采的基本方程，通過(guò)自主編寫(xiě)USDFLD子程序，對(duì)降壓開(kāi)采模型進(jìn)行了二次開(kāi)發(fā)，在ABAQUS中實(shí)現(xiàn)了降壓開(kāi)采條件下海底斜坡水合物儲(chǔ)層分解區(qū)域的判斷。

(2)得到了降壓開(kāi)采條件下天然氣水合物的擾動(dòng)范圍及飽和度的分布，滲透率與彈性模量的變化規(guī)律。在靠近開(kāi)采井的區(qū)域，水合物分解完全，其滲透率增大，彈性模量減小。

(3)在計(jì)算得到海底能源土斜坡水合物儲(chǔ)層分解區(qū)域基礎(chǔ)上，利用有限元強(qiáng)度折減法對(duì)斜坡的穩(wěn)定性進(jìn)行了數(shù)值計(jì)算，單井降壓開(kāi)采對(duì)斜坡的穩(wěn)定性影響不明顯。