黃翔 ，劉芳 *，巨鑫

1 同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092

2 同濟(jì)大學(xué)巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092

0 引言

天然氣水合物是由天然氣(主要是甲烷)分子與水分子在特定溫度和壓力條件下形成的一種籠形冰狀晶體。據(jù)估計(jì)，全球天然氣水合物儲(chǔ)量是煤炭、石油、常規(guī)天然氣資源總碳量的2 倍[1]，被認(rèn)為是未來綠色替代能源。中國(guó)天然氣水合物總資源量可達(dá)8.4×1013m3，主要分布在南海海域、青藏高原凍土區(qū)、東海海域及東北凍土區(qū)域，其總量約占我國(guó)天然氣水合物總資源量的77.7%[2]。因此，天然氣水合物的資源利用對(duì)我國(guó)具有非常重大的戰(zhàn)略意義。

天然氣水合物的開采利用帶來了一系列工程問題和挑戰(zhàn)。天然氣水合物廣泛賦存于水深為300～4000 m的海底大陸架[1]，垂直分布深度從海床表面到海床以下幾百米。在日益增長(zhǎng)的能源需求和科技進(jìn)步的驅(qū)動(dòng)下，海洋資源的開發(fā)和利用逐漸進(jìn)入深水，浮式深水基礎(chǔ)設(shè)施(如浮船式生產(chǎn)儲(chǔ)油卸油系統(tǒng)、半潛式生產(chǎn)或系泊系統(tǒng)、張力腿海洋平臺(tái)和立柱式海洋平臺(tái)等)[4]將廣泛分布，其基礎(chǔ)結(jié)構(gòu)可能穿越水合物富集區(qū)。水合物開采將導(dǎo)致其賦存土體力學(xué)性質(zhì)顯著劣化[5]，儲(chǔ)層應(yīng)力狀態(tài)改變，進(jìn)而影響深水基礎(chǔ)結(jié)構(gòu)的承載性能，給深水基礎(chǔ)設(shè)施的安全帶來隱患。但是這方面的研究目前還鮮有報(bào)道。

本文以浮式海洋平臺(tái)常用的錨板基礎(chǔ)為例，采用數(shù)值模擬方法研究水合物開采過程對(duì)臨近深水基礎(chǔ)抗拔承載性能的影響。含水合物沉積物是多組分多相多孔介質(zhì)復(fù)雜系統(tǒng)，水合物的開采過程涉及物質(zhì)變化、熱量傳遞、質(zhì)量交換和介質(zhì)變形，伴隨熱場(chǎng)、流場(chǎng)、力場(chǎng)和化學(xué)場(chǎng)的耦合變化。針對(duì)水合物開采的熱—流—力—化學(xué)(Thermo-Hydro-Mechanic-Chemical,THMC)耦合數(shù)值模擬研究技術(shù)仍然存在諸多挑戰(zhàn)。全耦合方法在降低計(jì)算成本、解決收斂困難方面尚待突破；將熱—流—化學(xué)(THC)分析與力學(xué)分析進(jìn)行解耦的時(shí)序耦合模擬仍然是主流方法，主要采用非等溫多相流分析程序與已有力學(xué)分析的開源或商業(yè)程序進(jìn)行耦合，用于分析水合物開采的工程力學(xué)效應(yīng)，如水合物開采過程的海床變形[6]和井壁穩(wěn)定性[7]等。本文借鑒了上述研究思路，基于美國(guó)勞倫斯—伯克利國(guó)家實(shí)驗(yàn)室開發(fā)的TOUGH + HYDRATE多相流分析程序[8]和大型有限元分析商業(yè)程序ABAQUS[11]，實(shí)現(xiàn)了水合物開采過程中水合物儲(chǔ)層—深水錨板基礎(chǔ)相互作用的THMC時(shí)序耦合數(shù)值模擬。

1 數(shù)值模型

1.1 單向時(shí)序耦合的模擬框架

圖1 所示為本文THMC時(shí)序耦合的基本框架。首先進(jìn)行時(shí)間離散，在每個(gè)時(shí)步內(nèi)THC耦合分析和力學(xué)分析分步依次進(jìn)行，分別由TOUGH+HYDRATE和ABAQUS完成，前者忽略土骨架變形的影響，因此屬于單向時(shí)序耦合方法。TOUGH+HYDRATE是針對(duì)水合物開采所開發(fā)的熱流分析程序，基于質(zhì)量與能量守恒方程求解水合物分解過程中的相態(tài)演變、質(zhì)量變化、熱量交換和流場(chǎng)變化，可以考慮水合物生成和分解的平衡和動(dòng)態(tài)模型[8]，該程序已被廣泛用于天然氣水合物的開采模擬[9-10]。

在THC耦合分析中，含水合物沉積物系統(tǒng)的多物理場(chǎng)耦合基本控制方程組(即質(zhì)量守恒方程和能量守恒方程)可寫成統(tǒng)一積分,形式如下 ：

當(dāng)上式為質(zhì)量守恒方程時(shí)，κ表示物質(zhì)，即甲烷(m)或水(w)；當(dāng)上式為能量守恒方程時(shí)，κ表示熱量。Mκ，F(xiàn)κ和qκ分別表示質(zhì)量(或熱量)增量、液體通量(或熱通量)和質(zhì)量(或熱量)的源或匯。Ω為單位方向向量為n的閉合面積Г所圍成的任意體積。求解式(1)可以得到孔壓場(chǎng)演變過程以及水合物分解鋒面擴(kuò)展。

如圖1 所示，每個(gè)時(shí)步內(nèi)，THC分析所得到的孔壓和水合物飽和度的空間分布傳遞給力學(xué)分析模塊，根據(jù)當(dāng)前水合物飽和度值更新水合物賦存土體的力學(xué)參數(shù)(如彈性模量E、黏聚強(qiáng)度c和內(nèi)摩擦角φ)。由于溫度變化范圍小，相較于滲流場(chǎng)和化學(xué)場(chǎng)而言，溫度場(chǎng)對(duì)力場(chǎng)的影響較小，在此未予考慮。

圖1 THMC單向時(shí)序耦合示意圖Fig. 1 Illustration of THMC one-way sequenced coupled process

圖2 模型示意圖和儲(chǔ)層參數(shù)Fig. 2 Schematic illustration of the model and reservoir parameters

1.2 數(shù)值建模

1.2.1 THC分析模型

本文模擬中國(guó)南海東沙某區(qū)域的豎井降壓開采水合物的過程及其對(duì)臨近錨板基礎(chǔ)承載性能的影響，場(chǎng)地條件參考2013 年GMGS2 調(diào)查航次8 號(hào)鉆井?dāng)?shù)據(jù)，開采條件和模型參數(shù)參考文獻(xiàn)[13]。圖2 為計(jì)算模型的示意圖。該場(chǎng)地的水深為800 m，海床溫度為4 ℃，地溫梯度為0.045 ℃/m，含兩個(gè)水合物儲(chǔ)層，分別位于海床以下9 m和63 m，厚度分別為14 m和35 m，淺部?jī)?chǔ)層的水合物飽和度為12%，深部?jī)?chǔ)層含厚度為3 m的高飽和度夾層(水合物飽和度55%)，該儲(chǔ)層除夾層外的水合物飽和度為40%。該場(chǎng)地各土層的參數(shù)見表1。

如圖3 所示，為了提高計(jì)算效率，THC分析和力學(xué)分析的計(jì)算域有所不同。其中，THC分析部分采用了以開采豎井為軸的二維軸對(duì)稱模型，計(jì)算域的半徑為250 m，高120 m，頂部和底部各設(shè)置厚0.1 m的隔水隔熱層作為邊界。開采井半徑為0.1 m，深98 m，采用降壓開采，壓強(qiáng)由靜水壓力(在8.00～8.98 MPa范圍)降為4.5 MPa。

假設(shè)開采井內(nèi)的流體運(yùn)動(dòng)也滿足達(dá)西流，開采井單元的孔隙度設(shè)為1，穿越儲(chǔ)層部分的完井段不設(shè)套管，滲透系數(shù)取 5.0×10-6m2，其余部分井筒設(shè)套管，滲透系數(shù)取1.0×10-18m2[14]。THC的計(jì)算網(wǎng)格劃分見圖3，對(duì)水合物層及靠近開采井區(qū)域進(jìn)行了加密。

表1 THC分析的模型參數(shù)表Table 1 Parameters used in the THC coupled analysis

圖3 THC模型網(wǎng)格劃分Fig. 3 Domain discretization of THC model

1.2.2 力學(xué)分析模型

圖2 虛線所示為力學(xué)分析的三維計(jì)算域，長(zhǎng)100 m、寬30 m、高98 m，底部邊界固定，側(cè)面邊界只允許豎向位移。計(jì)算網(wǎng)格見圖4，采用四結(jié)點(diǎn)線性四面體單元，錨板基礎(chǔ)附近進(jìn)行網(wǎng)格局部加密以提高計(jì)算精度。

方形錨板基礎(chǔ)寬3 m、厚0.3 m、埋深20 m，利用對(duì)稱性，取1/2 錨板，對(duì)稱面置于計(jì)算域側(cè)面，考慮錨板與開采井距離分別為20 m和40 m兩種工況。錨板采用彈性模型，楊氏模量200 GPa，泊松比0.3。為了避免非線性接觸帶來的收斂困難，將錨板基礎(chǔ)與土體的接觸模式設(shè)置為綁定約束，即假設(shè)兩者不發(fā)生相對(duì)運(yùn)動(dòng)。

地層的力學(xué)參數(shù)參考文獻(xiàn)[7]，土體采用摩爾—庫(kù)倫彈塑性模型[18]，按照水合物飽和度值分為4 類土體，主要參數(shù)匯總見表2。根據(jù)室內(nèi)水合物土體三軸壓縮試驗(yàn)結(jié)果[19-21]，隨著水合物飽和度的變化，內(nèi)摩擦角變化不明顯。鑒于錨板承載力的大小主要取決于土體的強(qiáng)度參數(shù)(即內(nèi)摩擦角和黏聚強(qiáng)度)，且剪脹角和泊松比變化規(guī)律的試驗(yàn)數(shù)據(jù)欠缺，本文采用了相同的泊松比和剪脹角，并假設(shè)土骨架大體相似，選用相同內(nèi)摩擦角，但考慮了水合物分解所引起的黏聚強(qiáng)度和彈性模量的降低。

在水合物開采過程中，分解鋒面持續(xù)由井周向外運(yùn)移，分解區(qū)域逐漸擴(kuò)展，在分解區(qū)域內(nèi)的土體強(qiáng)度和變形參數(shù)發(fā)生變化。由于分解鋒面處的水合物飽和度變化梯度較大，分解鋒面范圍較窄，為簡(jiǎn)化計(jì)算，忽略分解鋒面處飽和度變化影響，預(yù)設(shè)水合物飽和度閾值(6%)，將水合物儲(chǔ)層分為水合物分解區(qū)和未分解區(qū)，并分別采用相應(yīng)的材料參數(shù)。

力學(xué)分析包括9 個(gè)分析步，第1 個(gè)分析步為地應(yīng)力平衡，第2 個(gè)和第9 個(gè)分析步分別為開采前和開采30 年后停采靜水壓條件下的錨板拉拔試驗(yàn)?zāi)M，第3～8 個(gè)分析步分別經(jīng)歷不同開采時(shí)間(1 年、2 年、5年、10 年、20 年、30 年)下的錨板拉拔試驗(yàn)?zāi)M。錨板拉拔試驗(yàn)?zāi)M采用位移控制加載，模擬錨板在豎向拉拔情況下的地層變形及破壞情況，得到錨板的力與位移曲線，從而確定錨板在水合物開采過程中的拉拔承載性能的變化規(guī)律。

圖4 力學(xué)分析模型Fig. 4 The model for mechnical analysis

表2 力學(xué)參數(shù)表Table 2 Mechanical parameters used in the study

1.3 分析模塊的交互實(shí)現(xiàn)

THC和力學(xué)分析模塊的數(shù)據(jù)交互程序采用Python語言編寫，將THC分析所得到的孔壓場(chǎng)和水合物飽和度計(jì)算結(jié)果傳遞給力學(xué)分析模型，分析流程見圖5。

THC和力學(xué)分析的計(jì)算網(wǎng)格不同，需將前者結(jié)果插值到后者所采用的計(jì)算網(wǎng)格中，完成數(shù)據(jù)交互。即將力學(xué)分析模型中的任意節(jié)點(diǎn)以開采井為軸旋轉(zhuǎn)并映射到THC分析所采用的二維軸對(duì)稱平面上，然后在該平面上采用反距離權(quán)重插值法[12]，用通過二維軸對(duì)稱網(wǎng)格單元中心所得到的變量值進(jìn)行插值，獲得待求節(jié)點(diǎn)的相應(yīng)變量值，插值公式如下：

其中，S0為待求點(diǎn)節(jié)點(diǎn)坐標(biāo)，Si為已知點(diǎn)坐標(biāo)，λi為已知點(diǎn)Si的插值權(quán)重，d0i為待求點(diǎn)與已知點(diǎn)之間的距離。

圖 5 THC和力學(xué)分析模塊的交互示意圖Fig. 5 Communication between the THC and mechanical analysis

圖6 THC分析結(jié)果：(a) 產(chǎn)氣率；(b) 累積產(chǎn)氣量Fig. 6 Results of the THC analysis: (a) the gas production rate ; (b) the cumulative volume of produced gas

2 結(jié)果分析

2.1 THC模擬結(jié)果

2.1.1 產(chǎn)氣速率及累積產(chǎn)氣量

圖6 為產(chǎn)氣率和累積產(chǎn)氣量與時(shí)間的關(guān)系，與文獻(xiàn)[13]的分析結(jié)果大致吻合。如圖6 所示，水合物開采初期的產(chǎn)氣率較高，開采2 個(gè)月左右達(dá)到最高產(chǎn)氣率，約7400 STm3/d，之后產(chǎn)氣率迅速降低，開采7 年后降至4500 STm3/d左右，隨后產(chǎn)氣率趨于平穩(wěn)，開采15 年后產(chǎn)氣率逐漸降低，開采30 年后產(chǎn)氣率降至3800 STm3/d左右。而產(chǎn)氣量持續(xù)上升，30 年的總產(chǎn)氣量約為5.5×107STm3。

2.1.2 水合物分解鋒面擴(kuò)展

圖7 為水合物開采中孔隙壓力分布情況，圖中虛線框?yàn)榱W(xué)分析域范圍。最大降壓區(qū)位于開采井附近，隨著時(shí)間推移，由于水合物分解，儲(chǔ)層表觀滲透率有所增加，下儲(chǔ)層井周降壓梯度逐漸增加。

圖8 為水合物開采中水合物飽和度分布的變化。隨著開采進(jìn)行，水合物分解鋒面逐漸向遠(yuǎn)處延伸。由于溫度較高的海水沿壓強(qiáng)梯度從下方匯入，加速了下儲(chǔ)層水合物的分解，下儲(chǔ)層分解鋒面的擴(kuò)展速度明顯快于上儲(chǔ)層。

2.2 力學(xué)模擬結(jié)果

以開采10 年后的錨板基礎(chǔ)拉拔數(shù)值模擬實(shí)驗(yàn)為例，當(dāng)錨板豎向位移分別為0.3 m、0.6 m、0.9 m時(shí)的豎向有效應(yīng)力分布情況見圖9。錨板基礎(chǔ)承受拉拔荷載時(shí)，錨板上方的土體受壓，有效應(yīng)力增大，而錨板下方土體受拉，有效應(yīng)力較小。錨板對(duì)周圍土體的影響范圍有限，土體的受壓區(qū)隨著錨板基礎(chǔ)豎向位移的增大而增大。

圖10 給出了不同井距和開采時(shí)間下錨板基礎(chǔ)拉拔荷載—位移曲線，荷載均隨著位移增加而增加。考慮錨板基礎(chǔ)抗拔承載性能的充分發(fā)揮，選取位移為30%錨板寬度的拉拔荷載作為錨板基礎(chǔ)的極限抗拔承載力(如圖10 中的虛線所示)。不同井距情況下極限抗拔承載力隨著開采過程的變化見圖11。

如圖11 所示，在水合物開采初期，當(dāng)井距為20 m時(shí)，錨板的抗拔承載力明顯提高，因?yàn)榻祲洪_采導(dǎo)致井周土體有效圍壓增加，土體摩擦強(qiáng)度有所提高，進(jìn)而提高了錨板基礎(chǔ)的抗拔承載力。當(dāng)井距為40 m時(shí)，由于錨板基礎(chǔ)距離開采井較遠(yuǎn)，降壓區(qū)影響范圍有限，錨板周圍土體強(qiáng)度尚未受到降壓開采的影響，因此錨板抗拔承載力基本不變。

開采10 年后，不同井距下的錨板抗拔承載力均有所降低，錨板周圍的儲(chǔ)層由于水合物逐漸分解，導(dǎo)致儲(chǔ)層土體黏聚強(qiáng)度持續(xù)降低，引起承載力下降。開采20 年后，距井20 m處的錨板周圍儲(chǔ)層的水合物已全部分解，且土體應(yīng)力基本穩(wěn)定，因此錨板基礎(chǔ)的抗拔承載力基本不變；但距井40 m處儲(chǔ)層水合物還在持續(xù)分解，因此該處錨板基礎(chǔ)的承載力繼續(xù)下降。

圖7 水合物開采過程中孔隙壓力分布演化關(guān)系圖Fig. 7 Evolution of spatial distribution of pore pressure during gas production

圖8 水合物開采過程中水合物飽和度分布演化關(guān)系圖Fig. 8 Temporospatial distribution of hydrate saturation during gas production

圖9 開采10 年后的土體豎向有效應(yīng)力分布圖Fig. 9 Spatial distributions of the vertical effective stress at 10 years

停止開采后，水合物開采所產(chǎn)生的超靜孔壓逐漸消散，土體中的孔壓恢復(fù)至初始靜水壓狀態(tài)，由于孔壓回升導(dǎo)致土體中有效圍壓降低，土體摩擦強(qiáng)度降低，因此不同井距的錨板基礎(chǔ)承載力均明顯降低。經(jīng)過長(zhǎng)期開采后停采，不同井距處儲(chǔ)層的水合物均完全分解，儲(chǔ)層黏聚強(qiáng)度和應(yīng)力水平都大致相當(dāng)，因此停采后不同井距錨板基礎(chǔ)承載力降至相同水平，且低于水合物開采前的承載力值。拔承載力變化同時(shí)受控于兩種機(jī)制：(1)降壓引起土體有效圍壓增加導(dǎo)致承載力提高；(2)水合物分解引起土體黏聚強(qiáng)度降低導(dǎo)致承載力下降。當(dāng)錨板距離開采井較近時(shí)，第一種機(jī)制在開采初期占主導(dǎo)作用，因此錨板承載力有所提高；隨著水合物分解鋒面隨時(shí)間擴(kuò)展，第二種機(jī)制的影響逐漸顯現(xiàn)，錨板承載力降低。當(dāng)錨板距離開采井較遠(yuǎn)時(shí)，第二種機(jī)制占主導(dǎo)。當(dāng)水合物停止開采后，孔壓對(duì)承載力的影響機(jī)制逐漸消失，由于不可逆的儲(chǔ)層黏聚強(qiáng)度喪失，錨板基礎(chǔ)的承載力明顯低于開采前初始承載力。

圖10 不同開采時(shí)間下錨板基礎(chǔ)的荷載-位移曲線：(a) 距井20 m；(b) 距井40 mFig. 10 Load-displacement curves of the plate anchor at different times when the anchor is at (a) 20 m and (b) 40 m away from the well

圖11 不同井距下極限抗拔承載力隨開采時(shí)間的變化Fig. 11 The evolving ultimate uplift capacity of the anchor under different well-anchor distances

3 結(jié)論

本文開發(fā)了含水合物多孔介質(zhì)的多相流分析程序TOUGH+HYDRATE和有限元商業(yè)軟件ABAQUS的交互模塊，實(shí)現(xiàn)了水合物開采過程中儲(chǔ)層海床與錨板基礎(chǔ)相互作用的THMC單向時(shí)序耦合分析，根據(jù)南海水合物富集區(qū)的典型地層條件建立了儲(chǔ)層—錨板相互作用數(shù)值分析模型，研究了降壓豎井開采水合物過程中錨板基礎(chǔ)抗拔承載性能的變化規(guī)律，為深海水合物富集區(qū)的油氣開采和基礎(chǔ)設(shè)施平臺(tái)的深水基礎(chǔ)設(shè)計(jì)提供了理論參考。本文的主要結(jié)論如下：

(1)降壓法開采水合物過程中，錨板基礎(chǔ)的抗拔承載力的變化受控于兩種不同機(jī)制，一是孔壓下降引起儲(chǔ)層摩擦強(qiáng)度增加從而提高錨板抗拔承載力，二是水合物分解引起儲(chǔ)層黏聚強(qiáng)度減小從而降低錨板抗拔承載力。

(2)當(dāng)錨板基礎(chǔ)距離開采井較近時(shí)，上述兩種機(jī)制交替起主導(dǎo)作用，因此抗拔承載力隨著開采時(shí)間的推移先明顯提高隨后逐漸降低。

(3)當(dāng)錨板基礎(chǔ)距離開采井較遠(yuǎn)時(shí)，降壓的影響機(jī)制并不明顯，因此錨板抗拔承載力在開采初期大致保持穩(wěn)定，長(zhǎng)期開采之后明顯降低。

(4)當(dāng)水合物長(zhǎng)期開采后停采，降壓對(duì)錨板抗拔性能的影響機(jī)制逐漸消失，但由于水合物分解導(dǎo)致錨板基礎(chǔ)周圍儲(chǔ)層土體的黏聚強(qiáng)度已喪失，與開采前初始抗拔承載性能相比，停采后的錨板基礎(chǔ)抗拔承載力降低至20%左右。