郭旭洋，金 衍，盧運(yùn)虎，夏 陽(yáng)，韋世明

（1.中國(guó)石油大學(xué)（北京）油氣資源與工程全國(guó)重點(diǎn)實(shí)驗(yàn)室，北京 102249；2.中國(guó)石油大學(xué)（北京）石油工程學(xué)院，北京 102249； 3.中國(guó)石油大學(xué)（北京）理學(xué)院，北京 102249）

0 引言

天然氣水合物資源量大，常存在于陸域凍土和海底沉積層內(nèi)［1］。海域天然氣水合物主要通過(guò)降壓法、熱激法和注抑制劑法進(jìn)行開(kāi)采［2-5］。天然氣水合物試采現(xiàn)場(chǎng)常通過(guò)降壓法降低天然氣水合物儲(chǔ)層壓力，引起儲(chǔ)層內(nèi)的天然氣水合物相變［6-7］。然而，降壓法開(kāi)采誘發(fā)天然氣水合物分解的同時(shí)也會(huì)導(dǎo)致儲(chǔ)層力學(xué)性質(zhì)發(fā)生變化。由于海域水合物儲(chǔ)層多為弱膠結(jié)沉積物，導(dǎo)致儲(chǔ)層骨架出現(xiàn)變形、破壞、地層沉降等現(xiàn)象，影響儲(chǔ)層的穩(wěn)定性，嚴(yán)重時(shí)可能進(jìn)一步誘發(fā)地質(zhì)災(zāi)害，嚴(yán)重影響天然氣水合物降壓法開(kāi)采的安全性和效率，需要開(kāi)展針對(duì)性研究［8-10］。

多場(chǎng)耦合理論被廣泛運(yùn)用于海域天然氣水合物鉆采誘發(fā)的儲(chǔ)層巖土力學(xué)與巖石力學(xué)響應(yīng)研究中。Rutqvist 和Moridis［11］提出了一種聯(lián)合TOUGH+HYDRATE 和FLAC3D 軟件的模擬方法并結(jié)合Mohr-Coulomb 判據(jù)表征了水合物開(kāi)采誘發(fā)的儲(chǔ)層穩(wěn)定性問(wèn)題。這種方法能夠分別發(fā)揮TOUGHHYDRATE 在水合物分解、氣水流動(dòng)、熱傳遞方面的優(yōu)勢(shì)和FLAC3D 在計(jì)算固體骨架變形方面的優(yōu)勢(shì)。該方法未將巖石固體骨架變形的求解與水合物相變分解放在同一數(shù)值系統(tǒng)求解。在此基礎(chǔ)上，J.Kim 等［12］提出了一種雙向流固耦合算法，提高了孔隙壓力和骨架變形的求解精度。Sanchez 等［13］提出了一種全耦合方法，聯(lián)合室內(nèi)物理實(shí)驗(yàn)和天然氣水合物開(kāi)采現(xiàn)場(chǎng)尺度模擬，表征了降壓法開(kāi)采天然氣水合物過(guò)程中的流固熱化耦合規(guī)律。通過(guò)流固熱化多場(chǎng)耦合模擬，能夠預(yù)測(cè)降壓過(guò)程中的孔隙壓力和溫度傳遞時(shí)空演化規(guī)律，表征水合物分解過(guò)程中儲(chǔ)層巖石骨架的應(yīng)力釋放和體積變化。改變降壓方式和工程參數(shù)會(huì)對(duì)水合物儲(chǔ)層的變形產(chǎn)生明顯的影響［14］。在天然氣水合物儲(chǔ)層流固熱化耦合模型的控制方程推導(dǎo)和模型構(gòu)建的過(guò)程中，儲(chǔ)層骨架力學(xué)性質(zhì)對(duì)模型穩(wěn)定性和準(zhǔn)確性的影響較大，內(nèi)聚力和剛度等典型的強(qiáng)度和變形參數(shù)的取值需要進(jìn)行準(zhǔn)確的判斷。在多場(chǎng)耦合條件下，這些力學(xué)參數(shù)會(huì)影響儲(chǔ)層骨架的應(yīng)力應(yīng)變特征、剪切破壞規(guī)律和水合物分解及二次生產(chǎn)［15］。當(dāng)天然氣水合物儲(chǔ)層的壓力和有效應(yīng)力發(fā)生顯著變化后，塑性變形和塑性破壞等現(xiàn)象出現(xiàn)，并可能誘發(fā)儲(chǔ)層出砂，影響天然氣水合物儲(chǔ)層的流動(dòng)特性和井筒的穩(wěn)定性［16-17］。水合物儲(chǔ)層出砂過(guò)程可以通過(guò)計(jì)算流體力學(xué)方法和離散元法建立數(shù)值模型表征，水合物顆粒和砂顆粒的運(yùn)移過(guò)程中，顆粒受力情況復(fù)雜，受到流動(dòng)方向的影響十分顯著［18］。水平井降壓開(kāi)采誘發(fā)的儲(chǔ)層流固熱化響應(yīng)與直井開(kāi)采具有差異，水平井的布井方式和布井?dāng)?shù)量均會(huì)對(duì)多場(chǎng)耦合響應(yīng)造成影響，導(dǎo)致井周儲(chǔ)層內(nèi)的巖石骨架變形破壞規(guī)律出現(xiàn)不同［19］。李陽(yáng)等［20］提出了一種適用于南海神狐海域水合物儲(chǔ)層的多場(chǎng)耦合井壁穩(wěn)定性分析方法，為鉆井過(guò)程中的鉆井液密度、鉆井液溫度和鉆井時(shí)間的優(yōu)化提供了理論參考。孫金等［21］利用Duncan-Chang 模型描述水合物儲(chǔ)層的變形破壞特征，并結(jié)合滲流場(chǎng)、溫度場(chǎng)、變形場(chǎng)進(jìn)行了耦合建模，發(fā)現(xiàn)降壓開(kāi)采的井底壓力、儲(chǔ)層厚度、彈性模量、儲(chǔ)層滲透率等因素均控制降壓開(kāi)采誘發(fā)海床沉降的程度。由于水合物開(kāi)采過(guò)程中儲(chǔ)層沉降可能導(dǎo)致井筒出現(xiàn)位移和變形，井筒-土壤耦合模型能夠同時(shí)表征水合物儲(chǔ)層變形和井筒的應(yīng)力應(yīng)變規(guī)律，為安全試采提供工程參數(shù)的參考［22］。天然氣水合物儲(chǔ)層的流固熱化建模過(guò)程中，空間建模常采用二維或三維模型。二維模型能夠表征特定平面內(nèi)的儲(chǔ)層響應(yīng)，而且數(shù)值模型求解效率高。三維建模則能夠表征近井地帶多個(gè)方向的壓力、應(yīng)力應(yīng)變和溫度的變化，能夠比二維模型考慮更多的滲流流型，但也會(huì)顯著提高算力需求，導(dǎo)致數(shù)值模型的求解速度降低［23］。

已有研究說(shuō)明，多場(chǎng)耦合理論能夠?yàn)樘烊粴馑衔飪?chǔ)層降壓法開(kāi)采誘發(fā)的儲(chǔ)層壓力、變形、溫度、破壞等響應(yīng)提供數(shù)學(xué)建模依據(jù)，為海域天然氣水合物安全、高效鉆采提供工程參數(shù)優(yōu)化的參考。本文提出了一種海域天然氣水合物儲(chǔ)層全耦合流固熱化數(shù)值模型，模型考慮了水合物分解、熱傳導(dǎo)和儲(chǔ)層骨架彈塑性變形與破壞本構(gòu)。基于全耦合模型，提高了各物理場(chǎng)數(shù)值求解時(shí)的穩(wěn)定性和準(zhǔn)確性，分析典型天然氣水合物儲(chǔ)層和降壓參數(shù)條件下的近井儲(chǔ)層力學(xué)響應(yīng)規(guī)律，量化相關(guān)因素對(duì)海域天然氣水合物儲(chǔ)層降壓開(kāi)采誘發(fā)的多場(chǎng)耦合響應(yīng)及儲(chǔ)層沉降的影響特征。

1 控制方程

本研究通過(guò)連續(xù)介質(zhì)表征天然氣水合物儲(chǔ)層，考慮水合物分解過(guò)程儲(chǔ)層內(nèi)的質(zhì)量守恒、能量守恒、靜力平衡和水合物分解動(dòng)力學(xué)［19，23-24］。根據(jù)質(zhì)量守恒和滲流理論，儲(chǔ)層中氣和水的滲流表示為：

式中：φ——儲(chǔ)層孔隙度；Sw、Sg——各相的飽和度；ρw、ρg——各相的密度；vw、vg——各相的速度；sw、sg——各相的流入或流出質(zhì)量速度。

降壓法開(kāi)采過(guò)程中，水合物分解速度通過(guò)動(dòng)力學(xué)模型表示：

式中：RMH——水合物分解速率；kd——反應(yīng)速率；MMH——摩爾質(zhì)量；pe——相平衡壓力；pg——?dú)庀鄩毫?；As——反應(yīng)比表面積。

在多場(chǎng)耦合過(guò)程中，儲(chǔ)層內(nèi)溫度和壓力的變化會(huì)導(dǎo)致儲(chǔ)層固體骨架和孔隙流體發(fā)生體積變化，該過(guò)程表示為：

式中：x——與壓力變化相關(guān)的壓縮系數(shù)；β——與溫度變化相關(guān)的膨脹系數(shù)；p——壓力；T——溫度。

在表征儲(chǔ)層變形過(guò)程中，通過(guò)應(yīng)力張量和靜力學(xué)平衡的控制方程進(jìn)行計(jì)算：

式中：σ——應(yīng)力張量。

方程未考慮重力對(duì)應(yīng)力變化的作用。

由于本研究關(guān)注近井地帶的沉積物損傷與力學(xué)性質(zhì)劣化，對(duì)于儲(chǔ)層固體骨架采用彈塑性本構(gòu)：

式中：C——彈性張量；ε——應(yīng)變；εp——塑性應(yīng)變；α——有效應(yīng)力系數(shù)；I——二階單位張量。

近井地帶的儲(chǔ)層固體骨架破壞判據(jù)通過(guò)Mohr-Coulomb 模型和Drucker-Prager 模型表示，判斷降壓法開(kāi)采不同階段、不同位置的彈、塑性狀態(tài)及塑性破壞程度。

通過(guò)有限元方法對(duì)上述方程進(jìn)行空間離散，通過(guò)后向歐拉法進(jìn)行時(shí)間離散，構(gòu)建數(shù)值系統(tǒng)進(jìn)行求解。求解時(shí)采用全耦合法，使用稀疏矩陣直接求解器進(jìn)行矩陣求解。

2 數(shù)值模型建立

基于上述控制方程，建立海域天然氣水合物儲(chǔ)層降壓法開(kāi)采二維模型，針對(duì)近井短水平井降壓開(kāi)采的問(wèn)題進(jìn)行模擬和分析。圖1 所示為近井地帶40 m 長(zhǎng)度和40 m 厚度區(qū)域，水平井筒長(zhǎng)度2 m 位于模型中心。滲流場(chǎng)中，頂部、底部、左側(cè)、右側(cè)邊界均為封閉邊界，無(wú)流體流入和流出；溫度場(chǎng)中，4 個(gè)邊界均為隔熱邊界，無(wú)熱量流入和流出；應(yīng)力場(chǎng)中，左側(cè)、右側(cè)、底部邊界均為輥支撐邊界，頂部為自由邊界以表征沉降作用。由于全耦合模型對(duì)數(shù)值穩(wěn)定性要求高，如圖1 所示采用加密網(wǎng)格處理。儲(chǔ)層固體骨架楊氏模量為40 GPa，泊松比為0.2，骨架密度為2600 kg/m3，滲透率為2 mD，孔隙度為0.15，初始天然氣水合物飽和度為40%，初始孔隙壓力為14 MPa，井底生產(chǎn)壓力為3 MPa。儲(chǔ)層強(qiáng)度參數(shù)中，初始內(nèi)聚力為0.861 MPa，內(nèi)摩擦角為12°，降壓過(guò)程中，內(nèi)聚力變化可以表示為：

式中：c0——初始內(nèi)聚力；α、β——系數(shù)。

3 結(jié)果與討論

根據(jù)建模參數(shù)，模擬降壓法20 d 內(nèi)誘發(fā)的儲(chǔ)層壓力、溫度、有效應(yīng)力、水合物飽和度、沉降等參數(shù)變化規(guī)律，為海域天然氣水合物降壓開(kāi)采過(guò)程中井周地層的多場(chǎng)耦合響應(yīng)規(guī)律提供理論參考。

3.1 壓力、溫度和水合物分解規(guī)律

圖2 所示為降壓開(kāi)采0.5、5、10、20 d 時(shí)儲(chǔ)層內(nèi)的壓力分布情況。結(jié)果顯示水平井筒降壓能夠較快波及到近井地帶，降壓開(kāi)采初期（例如0.5 和5 d），由于水平井筒幾何形狀的影響，滲流形態(tài)為橢圓形流動(dòng)，導(dǎo)致壓降區(qū)域的擴(kuò)展不是嚴(yán)格的徑向流。當(dāng)降壓開(kāi)采時(shí)間進(jìn)一步增加，流動(dòng)形態(tài)逐漸偏向于徑向流，壓降均勻向儲(chǔ)層邊界波及。這說(shuō)明采用水平井筒進(jìn)行天然氣水合物降壓開(kāi)采時(shí)，儲(chǔ)層壓降規(guī)律不同于點(diǎn)狀壓降誘發(fā)的徑向流，而是前期為橢圓流動(dòng)，后期為近似徑向流的狀態(tài)。5 d 的壓力分布結(jié)果顯示，水平井筒降壓開(kāi)采能夠較快降壓前緣波及至近井地帶，至20 d 時(shí)儲(chǔ)層壓力已經(jīng)顯著下降。這說(shuō)明水平井筒降壓能夠建立近井地帶壓差，促進(jìn)水合物分解和向井內(nèi)流動(dòng)。在模擬時(shí)間內(nèi)，壓降能夠有效波及至井周20 m 范圍內(nèi)。

圖3所示為溫度隨降壓開(kāi)采時(shí)間變化。與壓力時(shí)空演化規(guī)律類似，溫度變化在初期也呈現(xiàn)近似橢圓形擴(kuò)展，后期近似徑向擴(kuò)展。這是由于氣水滲流和溫度在儲(chǔ)層的傳遞耦合程度更高導(dǎo)致的。對(duì)比圖2 的壓力分布和圖3 的溫度分布發(fā)現(xiàn)，壓降波及程度遠(yuǎn)大于降溫范圍，這是由于降壓誘發(fā)水合物分解，水合物分解吸熱、儲(chǔ)層降溫、溫度梯度導(dǎo)致熱傳遞等過(guò)程隨后發(fā)生，因此儲(chǔ)層壓降比溫度降低波及更遠(yuǎn)。

圖3 不同降壓時(shí)間的溫度分布Fig.3 Temperature in the reservoir at different depressurization time steps

圖4 所示為天然氣水合物飽和度隨著降壓開(kāi)采的進(jìn)行的變化規(guī)律，水合物分解區(qū)域隨著降壓時(shí)間的增加而增加，且分解區(qū)形狀近似橢圓形，橢圓形長(zhǎng)軸方向與水平井筒走向一致。這一結(jié)果與圖2 的壓降規(guī)律存在關(guān)系：壓降導(dǎo)致相平衡破壞，水合物開(kāi)始分解，由于降壓初期孔隙壓力降低近似橢圓形流動(dòng)，在長(zhǎng)軸方向更加明顯，導(dǎo)致水合物分解在橢圓長(zhǎng)軸（水平井筒方向）更加明顯。隨著降壓開(kāi)采時(shí)間增加，水合物分解區(qū)逐漸擴(kuò)大，降壓20 d 后，水合物分解前緣最遠(yuǎn)擴(kuò)展至井筒外3 m 左右。通過(guò)分析色標(biāo)梯度發(fā)現(xiàn)，降壓0.5、5 和10 d 后，水合物分解前緣明顯，由完全分解至未分解過(guò)渡很快；降壓20 d 后，水合物分解前緣相對(duì)平緩，水合物飽和度由0 過(guò)渡至40%的距離略微變長(zhǎng)。

圖4 不同降壓時(shí)間的天然氣水合物飽和度分布Fig.4 Natural gas hydrate saturation in the reservoir at different depressurization time steps

通過(guò)圖2～4 所示不同時(shí)刻的儲(chǔ)層壓力、溫度、水合物飽和度分布發(fā)現(xiàn)，在降壓開(kāi)采后的特定時(shí)間內(nèi)，壓力場(chǎng)、溫度場(chǎng)和水合物分解場(chǎng)呈現(xiàn)出相關(guān)性，水平井筒降壓開(kāi)采導(dǎo)致的多場(chǎng)耦合響應(yīng)應(yīng)當(dāng)區(qū)別于其他井型開(kāi)采誘發(fā)的儲(chǔ)層響應(yīng)。在溫度、壓力、水合物飽和度變化規(guī)律中，壓降波及最遠(yuǎn)，溫度降低波及的程度次之，水合物分解前緣推進(jìn)的速度最慢。這一現(xiàn)象是由于各物理場(chǎng)不同的控制機(jī)理導(dǎo)致的。

3.2 井周儲(chǔ)層力學(xué)響應(yīng)規(guī)律

降壓開(kāi)采誘發(fā)近井地帶儲(chǔ)層的多場(chǎng)耦合響應(yīng)，在已有溫度、壓力、水合物分解等特征的基礎(chǔ)上，進(jìn)一步分析近井儲(chǔ)層力學(xué)響應(yīng)規(guī)律。

圖5 所示為降壓開(kāi)采0.5 d 和20 d 的儲(chǔ)層σx′和σz′分布情況，分別表示沿水平井筒方向和垂向的有效正應(yīng)力。此應(yīng)力由降壓開(kāi)采誘發(fā)，能直接反映各方向儲(chǔ)層骨架受到的壓縮狀態(tài)應(yīng)力分布。結(jié)果顯示，由于水平井筒直接施加壓降，水平井筒區(qū)域有效應(yīng)力更高，且最大值隨著壓降時(shí)間的增加而上升。σz′在水平井筒兩端更易出現(xiàn)應(yīng)力集中現(xiàn)象。結(jié)果顯示，在近井地帶較小范圍內(nèi)建模的條件下，模型邊界均會(huì)由于壓力變化而出現(xiàn)有效應(yīng)力上升，且有效應(yīng)力上升的特征在水平方向和垂直方向均有不同。整體上，由于海域天然氣水合物常位于深水淺層，上覆地層應(yīng)力相對(duì)于水平方向地應(yīng)力較高，也會(huì)導(dǎo)致σx′和σz′分布出現(xiàn)差異。

圖5 不同降壓時(shí)間的水平有效正應(yīng)力和垂向有效正應(yīng)力分布Fig.5 Effective stresses in the horizontal and the vertical directions in the reservoir at different depressurization time steps

圖6所示為近井地帶儲(chǔ)層破壞規(guī)律云圖，通過(guò)塑性應(yīng)變（體積應(yīng)變）表示。由于水平井筒兩端存在應(yīng)力集中現(xiàn)象，井筒兩端出現(xiàn)較明顯的塑性應(yīng)變。由于降壓開(kāi)采引入較大的生產(chǎn)壓差，0.5 d 就誘發(fā)了明顯的井周塑性變化，塑性區(qū)主要出現(xiàn)在井筒和井周10 m 范圍內(nèi)，模型邊界塑性變形不明顯。隨著降壓開(kāi)采時(shí)間的增加（20 d），井周塑性區(qū)及塑性變形進(jìn)一步增加，但波及范圍并未擴(kuò)展至模型邊界，塑性破壞仍以近井區(qū)域?yàn)橹鳌?/p>

圖6 不同降壓時(shí)間的近井地帶塑性破壞區(qū)分布Fig.6 Plastic regions in the near-well area at different depressurization time steps

圖7 所示為降壓開(kāi)采20 d 后近井地帶內(nèi)聚力和沉降分布情況。由于近井地帶水合物分解及塑性演化，儲(chǔ)層內(nèi)聚力出現(xiàn)明顯降低，強(qiáng)度發(fā)生明顯劣化，內(nèi)聚力劣化的區(qū)域與水合物分解區(qū)域高度重合，體現(xiàn)了天然氣水合物對(duì)儲(chǔ)層力學(xué)強(qiáng)度的主控作用。水平井筒對(duì)應(yīng)區(qū)域的上覆地層出現(xiàn)垂向沉降，水平井筒正下方的下伏地層沉降程度略低于水平井筒以外區(qū)域，這一現(xiàn)象受到降壓誘發(fā)的多場(chǎng)耦合作用控制，在水平井筒上方儲(chǔ)層骨架壓縮明顯促進(jìn)沉降，在水平井筒下方降壓誘發(fā)的儲(chǔ)層骨架壓縮減緩地層沉降的程度。

圖7 降壓開(kāi)采20 d 后的近井區(qū)域內(nèi)聚力分布及沉降特征Fig.7 Plastic regions in the near-well area at different depressurization time steps

圖8 所示為水平井筒中心（0 m，0 m）及其以淺2 m 處（0 m，2 m）的內(nèi)聚力劣化及儲(chǔ)層沉降程度隨降壓開(kāi)采時(shí)間的關(guān)系。在該結(jié)果中，展示沉降程度的絕對(duì)值以表征降壓開(kāi)采對(duì)地層沉降的影響大小，正值代表地層沉降的距離。在水平井筒及其以淺區(qū)域，地層沉降均隨著時(shí)間單調(diào)遞增，且水平井筒以淺的觀測(cè)點(diǎn)沉降值更高。內(nèi)聚力劣化結(jié)果顯示，井筒中心由于在0 時(shí)刻就施加壓降載荷，導(dǎo)致其初期就出現(xiàn)明顯的水合物分解、破壞及內(nèi)聚力劣化，隨后，由于水合物完全分解，并已進(jìn)0 入塑性區(qū)，內(nèi)聚力變化不明顯。在井筒中心以前2 m 處，內(nèi)聚力在降壓開(kāi)采5 d 后開(kāi)始出現(xiàn)下降，表明水合物分解和塑性區(qū)演化影響到該處，隨后內(nèi)聚力在2 d 內(nèi)持續(xù)下降，直到水合物完全分解，其數(shù)值趨于穩(wěn)定。

圖8 水平井筒中心及其以淺2 m 處的內(nèi)聚力劣化及沉降演化規(guī)律Fig.8 The evolution patterns of cohesion and subsidence at the center of the horizontal wellbore and at 2 m above it

圖9 所示為降壓開(kāi)采過(guò)程中的壓力和水合物飽和度分布情況。結(jié)果顯示在20 d 內(nèi)，壓降有效波及至模型邊界，建立近井沉積物層內(nèi)壓差；水合物分解前緣擴(kuò)展相對(duì)緩慢，20 d 內(nèi)擴(kuò)展至短井筒以外5 m。

圖9 不同時(shí)刻近井區(qū)域壓力和水合物飽和度分布Fig.9 The distribution of pressure and hydrate saturation near the wellbore at different time steps

4 結(jié)論

本研究通過(guò)海域天然氣水合物儲(chǔ)層流固熱化耦合建模，形成了數(shù)值穩(wěn)定性較好的全耦合模型，分析了水平井筒降壓開(kāi)采誘發(fā)的近井地帶壓力、溫度、變形、破壞區(qū)變化規(guī)律，考慮了內(nèi)聚力劣化及塑性破壞對(duì)全耦合過(guò)程的影響。得到以下結(jié)論：

（1）降壓法開(kāi)采過(guò)程中，近井地帶流動(dòng)形態(tài)受到水平井筒的幾何尺寸影響，前期呈現(xiàn)近似橢圓形，后期呈現(xiàn)徑向特征；壓降波及程度較溫度降更廣，而水合物分解區(qū)擴(kuò)展速度最慢。

（2）降壓誘發(fā)的有效正應(yīng)力在水平井筒方向和垂向呈現(xiàn)不同的響應(yīng)特征，水平井筒兩端的位置出現(xiàn)應(yīng)力集中現(xiàn)象。這種響應(yīng)特征的差異化也與深水淺層的上覆地層應(yīng)力和水平地應(yīng)力場(chǎng)有關(guān)。

（3）降壓法誘發(fā)近井地帶的水合物分解和塑性區(qū)域演化對(duì)地層的強(qiáng)度劣化及沉降均有直接影響。內(nèi)聚力的劣化與水合物分解和塑性區(qū)演化關(guān)聯(lián)程度高。在水平井筒以淺處，地層沉降效果被水平井筒降壓增強(qiáng)；在水平井筒以深處，地層沉降一定程度上被水平井筒降壓所緩解，沉降程度可以降低約5 mm。