基于cohesive單元海域天然氣水合物儲層水力壓裂模擬

蔣亞峰，田英英，李小洋，韓澤龍，張欣，趙明，劉聃，尤志偉，梁金強(qiáng)，黃偉

（1.中國地質(zhì)科學(xué)院勘探技術(shù)研究所，河北 廊坊 065000；2.清華大學(xué)，北京 100083；3.中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局，廣東 廣州 511466）

0 引言

天然氣水合物是由水和甲烷氣體在低溫高壓條件下組成的結(jié)晶化合物［1］，普遍分布在深海海底沉積物和陸地永久凍土中［2］，據(jù)初步預(yù)測，賦存于水合物中的天然氣資源量達(dá)到3×1015m3，其中可開采量達(dá)到3×1014m3，被視為是21世紀(jì)一種重要的潛在的清潔替代能源［3-5］。由于海域天然氣水合物富集程度較低，分布范圍廣，因此，與其他常規(guī)化石能源相比，從天然氣水合物儲層中開采天然氣難度極大。目前，世界上僅有中國和日本2個(gè)國家進(jìn)行了海域天然氣水合物試開采。2013年，日本首次在Nankai海槽開展了海域天然氣水合物試采作業(yè)，經(jīng)過6 d作業(yè)，產(chǎn)氣量達(dá)11.95×104m3，平均日產(chǎn)氣量2×104m3［6-8］。2017年在相同海域，根據(jù)第一次試采的相關(guān)經(jīng)驗(yàn)，再次開展了第二次試開采。其中，第一口井由于井底嚴(yán)重出砂只開采了12 d，產(chǎn)氣4.1×104m3；第二口試采井持續(xù)開采了24 d，累計(jì)產(chǎn)氣22.25×104m3［9］。2017年，我國在南海神狐海域成功完成了天然氣水合物試采工程，實(shí)現(xiàn)了60 d內(nèi)持續(xù)產(chǎn)氣30.9×104m3［10］。上述3次海域水合物試采工程均采用單垂直井降壓法開采。2020年，我國首次采用水平井降壓法對南海神狐海域天然氣水合物進(jìn)行第二次開采，累計(jì)產(chǎn)氣達(dá)149.86×104m3，日均產(chǎn)氣3.57×104m3［11］。雖然通過多次嘗試在總產(chǎn)氣量和產(chǎn)氣速率上取得了極大的進(jìn)步，并驗(yàn)證了相關(guān)鉆井技術(shù)、開采工藝和相關(guān)裝備的可行性，但是與商業(yè)化開采標(biāo)準(zhǔn)還有一段距離。

與日本水合物試采區(qū)的砂質(zhì)地層對比，我國南海海域水合物主要賦存于泥質(zhì)粉砂地層中，其賦存環(huán)境更加復(fù)雜，泥質(zhì)含量高導(dǎo)致滲透率低，嚴(yán)重制約了降壓開采的效果，從而使得常規(guī)的降壓開采難以獲得理想的產(chǎn)氣量［12］。其實(shí)，很早之前便已有研究人員提出將水力壓裂應(yīng)用于海域水合物儲層改造中，利用水力壓裂形成的高導(dǎo)流裂縫增加了降壓效果和水合物分解面積，從而提高產(chǎn)氣量［13］。張永勤等［14］提出利用旋轉(zhuǎn)噴射鉆井將改造劑注入到水合物儲層中形成蜂窩狀支柱，從而起到提供流動通道和支撐地層的作用。彭賽宇等［15］提出通過對儲層灌注特殊的泡沫水玻璃漿液改造劑，形成方狀骨架結(jié)構(gòu)，從而達(dá)到提高儲層滲透率和增大產(chǎn)氣量的效果。KONNO等［16］針對含水合物砂質(zhì)沉積物，開展了不同側(cè)重的水力壓裂試驗(yàn)。含甲烷水合物砂質(zhì)沉積物水力裂縫表現(xiàn)出拉伸破壞特征，擴(kuò)展方向垂直于最小主應(yīng)力方向，與完整固結(jié)巖石壓裂行為相似。雖然相關(guān)理論提出了很多，但是，針對水合物儲層的水力壓裂裂縫的起裂機(jī)理、擴(kuò)展規(guī)律方面的研究較少。因此，研究水合物儲層水力壓力過程中裂縫的起裂、擴(kuò)展和形態(tài)，對后續(xù)分析儲層水合物分解速率以及推動海域天然氣水合物的安全高效開采具有重要意義?；诖耍狙芯吭谇叭死碚摰幕A(chǔ)上，利用abaqus軟件建立基于cohesive單元的二維水力壓裂數(shù)值模型，研究壓裂過程中的儲層基質(zhì)變形、裂縫起裂和擴(kuò)展、裂縫形態(tài)等參數(shù)，以驗(yàn)證采用cohesive單元法分析水力壓裂及水合物儲層采用水力壓裂進(jìn)行改造的可行性。

1 水力壓裂基本方程

1.1 cohesive單元簡介

在Abaqus軟件中，通過在原始單元內(nèi)插入cohesive單元，將其分為上下兩面，其中上表面和下表面的節(jié)點(diǎn)距離表示裂縫寬度，中間面的節(jié)點(diǎn)表示壓裂切向流動的狀態(tài)，如圖1所示。該方法可以模擬壓裂過程中的儲層基質(zhì)變形、裂縫起裂和擴(kuò)張、裂縫形態(tài)等參數(shù)，并且在頁巖氣和煤層氣等非常規(guī)能源儲層水力壓裂模型分析中得到了成功運(yùn)用［17-19］，為后續(xù)研究其他特殊地層水力壓裂裂縫擴(kuò)展規(guī)律提供了基礎(chǔ)［20］。

圖1 cohesive單元示意Fig.1 Schematic diagram of a cohesive unit

1.2 水力壓裂裂縫擴(kuò)展方程

CHEN［20］和ZOU等［21］提出，采用牽引擴(kuò)展定律描述水合物儲層發(fā)生破壞前的線性彈性行為、裂縫起裂和擴(kuò)展過程，如圖2所示。隨著流體的注入，儲層介質(zhì)牽引力不斷增大，壓裂儲層發(fā)生線性彈性行為，裂縫尖端開始不斷延伸，當(dāng)注入流體達(dá)到內(nèi)部彈性模量Tmax時(shí)，擴(kuò)展達(dá)到臨界值δ0。當(dāng)其大于δ0時(shí)，牽引力隨著裂縫的擴(kuò)展而減小即發(fā)生損傷演化行為，直到分離達(dá)到臨界值δf（壓裂的裂紋尖端），此時(shí)牽引力或注入流體強(qiáng)度消失，裂紋尖端完全停止延伸［22］。

圖2 牽引擴(kuò)展定律Fig.2 Law of traction propagation

在單一類型荷載下，GONG等［23］也指出牽引擴(kuò)展定律在分層過程中具有線性荷載和線性退化兩大特征：

式中：k0——cohesive單元表面剛度；δ0、δf——分別為單元界面的初始形變和完全破壞的相對位移；D——截面破壞損傷因子，其取值在0和1之間，當(dāng)D為0時(shí)，表示材料完好無損，當(dāng)D為1時(shí)，表示材料已完全破壞。

應(yīng)用二次應(yīng)力準(zhǔn)則則可以判斷水力壓裂裂縫是否開裂，即當(dāng)cohesive網(wǎng)格單元任一方向的三向應(yīng)力與其對應(yīng)的強(qiáng)度極限之比等于1時(shí)，水合物儲層開始破裂：

水合物儲層在完全起裂后，牽引擴(kuò)展定律中的損傷演化模型如下：

同時(shí)通過線性位移擴(kuò)展準(zhǔn)則，可以推出損傷因子D的公式：

1.3 流體流動方程

圖3為水合物儲層壓裂流體的切向流動和法向流動示意圖。由圖3可知，根據(jù)水合物儲層壓裂的流固耦合關(guān)系模型，可以看出壓裂液主要由切向流動和法向流動2部分組成，分別導(dǎo)致裂縫的擴(kuò)展和壓裂液濾失。

圖3 流體的切向流動和法向流動示意Fig.3 Tangential and normal flow patterns of fluid

流體切向流動方程：

流體法向流動方程：

式中：q——切向流量；Δp——單元長度方向壓力梯度；ω——裂縫寬度；μ——壓裂液粘度；qt、qb——分別為單元上下表面的壓裂液流速；ct、cb——分別為上下表面的滲濾系數(shù)；pt、pf、pb——分別為上、中、下表面的流體壓力。

1.4 孔隙介質(zhì)固體顆粒的控制方程

水力壓裂過程會改變地層中巖石的受力場，對于巖石這類的多孔介質(zhì)，其應(yīng)力平衡方程可以應(yīng)用虛功原理描述：固定時(shí)間內(nèi)巖石的虛功，等于作用在整個(gè)巖石上的體力與面力產(chǎn)生的虛功的疊加［24］。

巖體介質(zhì)的平衡方程為［25］：

式中：δε——虛應(yīng)變率矩陣；t——面力矩陣；δv——虛速度矩陣；f——體力矩陣；dV——單元體。

2 水力壓裂有限元模型

本文的研究對象為水合物儲層，與常規(guī)油氣儲層相比巖性較軟，其特殊的賦存形態(tài)形成了復(fù)雜的地層結(jié)構(gòu)，加之不穩(wěn)定性、沉積物和應(yīng)力狀態(tài)等多種因素影響，水力壓裂過程極其復(fù)雜，目前相關(guān)的研究均處于起步階段。因此，本文借助泥巖、頁巖水力壓裂相關(guān)模型基礎(chǔ)，并只選取雙向單縫為研究目標(biāo)，簡化模型。

2.1 數(shù)值模擬參數(shù)選取

2.1.1 儲層彈性模量和泊松比的選取

Lijith等［26］等借助大量的擬合數(shù)據(jù)結(jié)果，提出了水合物沉積物的彈性模量計(jì)算公式：

式中：E——水合物儲層彈性模量；——有效圍壓；SH——水合物飽和度；a——代表原始含天然氣水合物沉積物的初始剛度的常數(shù)，大小取決于沉積物的組成及地層壓力；b——代表受水合物形態(tài)、水合物形成條件及飽和度影響的常數(shù)。

WILLIAM等［27］通過單軸壓裂試驗(yàn)測得水合物儲層等溫楊氏模量7800 MPa，泊松比0.317；DVORKIN等［28］測得天然氣水合物儲層等溫楊氏模量為7900 MPa。通過對比，本次選取彈性模量7800 MPa，泊松比0.317。

2.1.2 滲透率參數(shù)的選取

DAI等［29］基于Kozeny-Carman模型提出了水合物儲層的滲透率模型，該模型與Kumar等［30］、Johnson等［31］和Kleinberg等［32］的實(shí)驗(yàn)結(jié)果取得較好的吻合，其公式為：

式中：K′—水合物儲層的滲透率；K0——水合物中沉積物的滲透率；SH——水合物飽和度。

2.1.3 孔隙比的計(jì)算

水合物儲層的有效孔隙度受水合物飽和度的影響，具體公式為：

式中：φ′——水合物儲層的有效孔隙度；φ——無水合物狀態(tài)下的孔隙度。

孔隙比是指水合物儲層中有效孔隙占水合物與沉積物體積之和的百分比，其計(jì)算公式為：

式中：e——水合物儲層孔隙比。

2.2 二維水力壓裂有限元模型

為了更好地分析儲層在水力壓裂過程中的裂縫擴(kuò)展及應(yīng)力、應(yīng)變變化規(guī)律，對二維水力壓裂數(shù)值模擬做出如下假設(shè)：

（1）水合物儲層為均質(zhì)地層；

（2）孔隙壓力為3.3 MPa；

（3）儲層上下、左右邊界固定，不產(chǎn)生位移。

同時(shí)，分別選用20 m×20 m、100 m×100 m兩種尺度的模型進(jìn)行模擬，其中20 m×20 m仿真模型主要用于研究裂縫及周邊的應(yīng)力、應(yīng)變，100 m×100 m模型主要用于研究裂縫擴(kuò)展方向、延伸機(jī)理方面的研究。此外還設(shè)置了水平方向和豎直方向的最大地應(yīng)力σ=5 MPa與最小地應(yīng)力σz=4 MPa，模型中心位置設(shè)置為壓裂液的注入點(diǎn)，cohesive預(yù)設(shè)單元為裂縫發(fā)展方向，對稱分布在井筒兩端，如圖4所示。2個(gè)模型均以（x，z）為坐標(biāo)的格式離散劃分成100×100=10000個(gè)cohesive網(wǎng)格塊。模型選取的主要參數(shù)，如表1所示。

表1 水合物儲層水力壓裂模型主要參數(shù)Table 1 Main parameters of the hydrate reservoir hydraulic fracturing model

3 結(jié)果與分析

本文分別對20 m×20 m、100 m×100 m兩個(gè)模型尺寸下的裂縫擴(kuò)展進(jìn)行了數(shù)值模擬，得出了水合物儲層水力壓力過程中裂縫擴(kuò)展的最大半長及最大寬度，通過2種尺寸下應(yīng)力、應(yīng)變、孔壓的分布狀態(tài)分析最優(yōu)模擬結(jié)果，最后對比不同時(shí)間段裂縫擴(kuò)展的速率，得出最優(yōu)壓裂時(shí)長，為后續(xù)的水合物水力壓裂模擬提供指導(dǎo)。

3.1 裂縫擴(kuò)展尺寸

模擬結(jié)果，如表2所示，2種模型尺寸下裂縫的半長均為6 m，大尺寸模型裂縫最大寬度為5.8 mm，小尺寸模型裂縫最大寬度5.5 mm。2組模擬結(jié)果基本吻合，驗(yàn)證了cohesive單元法應(yīng)用于水合物水力壓裂研究的可行性，為研究水合物儲層裂縫變化、擴(kuò)展機(jī)理提供了良好的基礎(chǔ)。

表2 模擬結(jié)果Table 2 Simulation results

3.2 仿真結(jié)果分析

模型建立的過程中，我們分別設(shè)立了（x，z）2個(gè)方向的邊界值，即水平方向和豎直方向的最大地應(yīng)力σx=5 MPa與最小地應(yīng)力σz=4 MPa。在模擬的過程中，隨著壓裂液的不斷注入，裂縫不斷擴(kuò)張的同時(shí)，其內(nèi)部的壓裂液沿法向切力的方向作用于裂縫內(nèi)表面，并不斷向外輻射。如果應(yīng)力位移傳遞量設(shè)置過小，很容易導(dǎo)致應(yīng)力不能有效傳遞，在邊界位置產(chǎn)生變形，影響模擬結(jié)果的準(zhǔn)確性。

圖5為2種模型尺寸裂縫擴(kuò)展模擬結(jié)果，從圖5可以明顯看到20 m×20 m裂縫擴(kuò)展應(yīng)變圖中應(yīng)力在傳播的過程中在邊界位置立即出現(xiàn)變形，而100 m×100 m模型中應(yīng)力呈扇形均勻的擴(kuò)散，模擬效果更加真實(shí)；同時(shí)，可以明顯看出20 m×20 m裂縫擴(kuò)展尖端應(yīng)力和孔壓圖中受位移傳遞量的限制存在發(fā)散現(xiàn)象，而100 m×100 m模型中裂縫延展尖端的應(yīng)力和孔壓變化不受邊緣的限制，模擬的區(qū)域圖形呈現(xiàn)規(guī)律性和對稱性。所以對比模擬結(jié)果后裂縫最大寬度在100 m×100 m模型下的5.8 mm更加準(zhǔn)確。

圖5 裂縫擴(kuò)展模擬結(jié)果Fig.5 Simulation results of crack propagation

3.3 裂縫寬度隨時(shí)間的變化規(guī)律

圖6為注入點(diǎn)裂縫寬度隨注入時(shí)間與壓力的變化特征，從圖中可以看出裂縫的發(fā)生形態(tài)與注入壓力的關(guān)系，隨著注入壓力的急劇增加，儲層開始起裂，裂縫瞬間張開，此時(shí)注入壓力最大，達(dá)到了25 MPa，然后再經(jīng)歷一段波動后趨于穩(wěn)定，注入壓力約為18 MPa。裂縫寬度在起裂后，以小幅度波動的方式緩慢變大。原因在于，當(dāng)開始向儲層注入壓裂液時(shí)，在注入點(diǎn)附近很容易形成憋壓導(dǎo)致裂縫瞬間增大。起裂發(fā)生后，裂縫瞬間沿預(yù)設(shè)單元擴(kuò)展方向裂開。注入壓力出現(xiàn)微降，在水平地應(yīng)力的作用下，裂縫寬度出現(xiàn)短暫的收縮，后續(xù)隨著注入液的補(bǔ)充，開始逐步趨于緩慢的增長狀態(tài)，形成這種“階梯式”的裂縫擴(kuò)展形態(tài)。

圖6 裂縫寬度隨時(shí)間與壓力的變化規(guī)律Fig.6 Variation law of crack width vs time and injection pressure

4 結(jié)論

本文通過利用abaqus軟件建立基于cohesive單元的二維水力壓裂數(shù)值模型，研究壓裂過程中裂縫起裂和擴(kuò)展、裂縫形態(tài)等參數(shù)，得出了水合物儲層水力壓裂裂縫擴(kuò)展的最大半長及最大寬度。具體結(jié)論如下：

（1）采用100 m×100 m和20 m×20 m大、小兩尺寸模型分別建立相同參數(shù)下的水合物二維水力壓裂裂縫擴(kuò)展模型，最終得出2組基本相同的數(shù)據(jù)結(jié)果：裂縫長度均為6 m，最大寬度分別為5.8、5.5 mm，論證了cohesive單元法分析水合物水力壓裂方案的可行性。

（2）在二維水合物儲層水力壓裂模擬中，通過對2種尺寸模型數(shù)值模擬結(jié)果的分析，在利用cohesive單元法分析水合物水力壓裂模型中應(yīng)盡量構(gòu)建尺寸較大的模型，這樣更利于應(yīng)力的傳遞及裂縫擴(kuò)展的模擬，得出的實(shí)驗(yàn)結(jié)果更加準(zhǔn)確。

（3）裂縫在擴(kuò)展的過程中會伴隨著階梯性延展，一方面，由于剛開始注壓時(shí)，較高的排量很容易在注入點(diǎn)附近形成憋壓導(dǎo)致裂縫瞬間增大，并出現(xiàn)了短暫的波動。同時(shí)，隨著預(yù)設(shè)單元擴(kuò)展方向迅速裂開，此時(shí)裂縫中的壓裂液內(nèi)壓會出現(xiàn)微降，在水平地應(yīng)力的作用下，裂縫寬度出現(xiàn)短暫的收縮的現(xiàn)象，后續(xù)隨著注入液的補(bǔ)充，開始逐步趨于緩慢的增長。這種急速擴(kuò)張—收縮的現(xiàn)象形成了水合物儲層“階梯式”的水力壓裂擴(kuò)展規(guī)律。該研究能夠?yàn)楹Ｓ蛱烊粴馑衔飪铀毫烟峁┮欢ǖ睦碚撝笇?dǎo)，提高水合物儲層的滲流能力。