郭旭洋，金 衍，盧運虎，夏 陽，韋世明

（1.中國石油大學(xué)（北京）油氣資源與工程全國重點實驗室，北京 102249；2.中國石油大學(xué)（北京）石油工程學(xué)院，北京 102249； 3.中國石油大學(xué)（北京）理學(xué)院，北京 102249）

0 引言

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

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

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

1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

2 數(shù)值模型建立

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 降壓開采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)聚力劣化及儲層沉降程度隨降壓開采時間的關(guān)系。在該結(jié)果中，展示沉降程度的絕對值以表征降壓開采對地層沉降的影響大小，正值代表地層沉降的距離。在水平井筒及其以淺區(qū)域，地層沉降均隨著時間單調(diào)遞增，且水平井筒以淺的觀測點沉降值更高。內(nèi)聚力劣化結(jié)果顯示，井筒中心由于在0 時刻就施加壓降載荷，導(dǎo)致其初期就出現(xiàn)明顯的水合物分解、破壞及內(nèi)聚力劣化，隨后，由于水合物完全分解，并已進0 入塑性區(qū)，內(nèi)聚力變化不明顯。在井筒中心以前2 m 處，內(nèi)聚力在降壓開采5 d 后開始出現(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 所示為降壓開采過程中的壓力和水合物飽和度分布情況。結(jié)果顯示在20 d 內(nèi)，壓降有效波及至模型邊界，建立近井沉積物層內(nèi)壓差；水合物分解前緣擴展相對緩慢，20 d 內(nèi)擴展至短井筒以外5 m。

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

4 結(jié)論

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

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

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

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