李 陽， 程遠(yuǎn)方, 閆傳梁， 周曉輝

(中國石油大學(xué)(華東)石油工程學(xué)院， 青島 266580)

天然氣水合物具有資源量大、能量密度高和低污染等優(yōu)點(diǎn)[1]，被公認(rèn)為是具有良好前景的重要后續(xù)清潔能源[2]，全球天然氣水合物的碳儲量相當(dāng)于現(xiàn)已探明的石油、煤和天然氣等常規(guī)化石能源總含碳量的2倍[3-5]。受水合物相平衡條件和地溫梯度的限制，天然氣水合物主要分布在深海和永久凍土區(qū)，埋存深度淺，沉積物膠結(jié)程度差[6-8]。各國學(xué)者對于水合物的力學(xué)性質(zhì)展開了大量研究[9-11]，發(fā)現(xiàn)水合物儲層具有極強(qiáng)的蠕變特性[12]。由于水合物在儲層中起到了很強(qiáng)的膠結(jié)作用，因此水合物飽和度的不同會對其蠕變特性產(chǎn)生非常大的影響。對于不同的水合物藏，其水合物飽和度必然不同，同時鉆井和開采等工程作業(yè)也會導(dǎo)致水合物分解，影響其飽和度和蠕變特性[13]。

由于鉆井工程的擾動，井眼附近會產(chǎn)生應(yīng)力集中，對于蠕變性地層，應(yīng)力的增大會使其產(chǎn)生蠕變變形。蠕變對鉆井最大的影響是會造成井眼縮徑，影響起下鉆作業(yè)。目前對于蠕變地層井壁穩(wěn)定的研究主要在鹽巖地層和軟泥巖地層，近年也開始關(guān)注脆性巖石蠕變。殷有泉等[14]采用理論計算和數(shù)值模擬的方法，分析了蠕變地層中套管載荷，認(rèn)為套管載荷的黏彈性解最終趨近于彈性解。Willson等[15]、Poiate等[16]建立了二維有限元或有限差分模型，分析了水平地應(yīng)力作用下鹽膏層中套管載荷。Machay等[17]、Wang等[18]等建立了三維有限元或有限差分模型，分析了三維地應(yīng)力作用下鹽膏層中套管力學(xué)狀態(tài)。Machay等[19]又提出了分步模擬方法，對鉆井和固井等過程井眼的力學(xué)狀態(tài)進(jìn)行了模擬分析。Lao等[20]建立了鹽膏層中斜井井筒的有限元模型，分析了井斜角和方位角對井眼縮徑和套管力學(xué)狀態(tài)的影響規(guī)律。Li等[21]建立了熱流固三場耦合的凍土地層井眼收縮力學(xué)模型，并分析了鉆井液溫度、地層溫度及井眼打開時間對井眼狀態(tài)的影響。

這些不同地層井眼縮徑的研究對于水合物地層鉆井作業(yè)過程中井眼的變形分析有很大的借鑒意義，但水合物地層的蠕變特性及其影響因素都與其他地層有很大的不同，目前關(guān)于水合物儲層蠕變特性對井眼縮徑規(guī)律的影響還鮮有研究。因此，現(xiàn)從應(yīng)力的角度對鉆穿水合物儲層時井眼的蠕變縮徑規(guī)律進(jìn)行研究，并分析水合物儲層井眼蠕變縮徑的影響因素及其作用規(guī)律。

1 數(shù)學(xué)模型

1.1 流固耦合模型

地層的應(yīng)力場與滲流場是相互作用、互相影響的。應(yīng)力場作用使得巖石骨架發(fā)生變形，改變地層的滲流特征，同時滲流場的變化又會改變地層的孔隙壓力分布從而影響巖石骨架承受的有效應(yīng)力，孔隙壓力對地層有效應(yīng)力影響的計算式為

σ′=σ-αpp

(1)

式(1)中：σ′ 為有效應(yīng)力；σ為地層中任一點(diǎn)的總應(yīng)力；α為有效應(yīng)力系數(shù)；pp為孔隙壓力。

由于實(shí)際的多孔介質(zhì)微觀組成復(fù)雜，很難用精確的數(shù)學(xué)關(guān)系將流體固體分開，因此所建立的模型為理想的連續(xù)系統(tǒng)，忽略巖石孔隙的微觀結(jié)構(gòu)變化，但通過一些特殊變量的變化可以表征孔隙流體流動的過程，以此模擬流固耦合的過程。

對于固相材料的應(yīng)力變化采用虛功原理表征[22]，即

(2)

式(2)中：V為計算區(qū)域體積；Sa為牽引力作用下的表面積；σ為總應(yīng)力矩陣；δε為虛應(yīng)變率矩陣；ts為表面牽引力矢量；f為體力矢量；δv為虛速度矢量；：為雙點(diǎn)積。

通過拉格朗日法將固相的虛功方程離散到網(wǎng)格之中，是流體流經(jīng)這些網(wǎng)格，流體滿足連續(xù)性方程，如式(3)所示，流體的滲流過程遵循Darcy定律。

(3)

式(3)中：pf為孔隙流體密度；φ為多孔介質(zhì)孔隙度；vfp為滲流速率矩陣；n為外表面外法線方向向量。

1.2 水合物儲層蠕變模型

蠕變分析是用于表征模型受力發(fā)生彈性變形后基于其黏塑性屬性而發(fā)生的一種后效反應(yīng)。單軸蠕變實(shí)驗中，當(dāng)材料處于彈性階段進(jìn)行蠕變實(shí)驗時，其總應(yīng)變可分解為彈性應(yīng)變和蠕變應(yīng)變之和，即

ε(t)=εe+εc

(4)

式(4)中：ε(t)為總應(yīng)變；εe為彈性應(yīng)變；εc為蠕變應(yīng)變，是溫度T、時間t以及應(yīng)力σ的函數(shù)，即

εc=f(σ,t,T)=f1(σ)f2(t)f3(T)

(5)

采用時間硬化法則時，σ為有效應(yīng)力，f1(σ)采用Norton冪律法則，則其初始蠕變和第二蠕變階段可表達(dá)[23]為

(6)

式(6)中：A、m、n為蠕變參數(shù)，由蠕變實(shí)驗獲得。

將式(5)對時間微分并化簡即可獲得時間硬化法則，即

(7)

1.3 水合物儲層蠕變模型

實(shí)驗研究表明，低溫水合物具有較強(qiáng)的蠕變特性，其蠕變曲線與富冰凍土在溫度接近0 ℃時蠕變曲線接近[24-25]，其蠕變本構(gòu)模型為

(8)

式(8)中：ε為應(yīng)變；δ為軸向差應(yīng)力；θ為負(fù)溫；t為時間；A、B、C、D為實(shí)驗擬合參數(shù)。擬合結(jié)果如圖1所示。

圖1 凍土水合物層蠕變模型擬合結(jié)果[24]Fig.1 Creep curves of the sediments with gas hydrates[24]

2 數(shù)學(xué)模型

采用ABAQUS有限元軟件對所建立的數(shù)學(xué)模型進(jìn)行計算。ABAQUS軟件是一個功能強(qiáng)大的非線性有限元軟件，廣泛應(yīng)用于巖土力學(xué)的相關(guān)計算[26-27]。同時ABAQUS擁有豐富材料數(shù)據(jù)庫可滿足大多數(shù)工程需要，UMAT子程序自主定義任意材料屬性。

2.1 水合物儲層井眼模型

水合物儲層開采流固耦合力學(xué)模型如圖2所示，儲層半徑為2.5 m，AB邊為井眼，井眼半徑為0.25 m，分別在AB邊、CD邊和DE邊設(shè)置孔壓邊界，CD、DE同時設(shè)置位移邊界；BC邊和AE邊為軸對稱邊界。模型單元類型為CPE4P，即4節(jié)點(diǎn)平面應(yīng)變雙線性位移、雙線性孔壓單元，模型總網(wǎng)格數(shù)為3 200。

圖2 流固耦合力學(xué)模型示意圖Fig.2 Schematic diagram of fluid-solid coupling mechanics model

2.2 水合物儲層井眼模型

假設(shè)該水合物層所在深度為500 m。模擬所需數(shù)據(jù)如表1[28-30]所示。

表1 地層基本參數(shù)[28-30]Table 1 Basic parameters of the model[28-30]

2.3 地應(yīng)力平衡

地應(yīng)力平衡是用于模擬地層在地應(yīng)力作用下的狀態(tài)，在上覆巖層壓力、最大水平地應(yīng)力、最小水平地應(yīng)力作用下地層處于應(yīng)力平衡狀態(tài)，當(dāng)使用ABAQUS模擬地應(yīng)力平衡不會出現(xiàn)明顯的位移和應(yīng)變時，說明所確定的地應(yīng)力參數(shù)是合理的。地應(yīng)力平衡通常是使用ABAQUS模擬巖土問題時的第一步。圖3顯示了模型模擬的結(jié)果，土體內(nèi)的最大位移為10-18m的量級，這表明在Geostatic步驟中給出的所有初始設(shè)置在隨后的分析中是有效的。

圖3 地應(yīng)力平衡Fig.3 Initial in-situ stress equilibrium

3 水合物層井眼縮徑規(guī)律

3.1 井眼鉆開時初始變形

使用鉆頭鉆開地層的實(shí)質(zhì)是地應(yīng)力的釋放過程，在ABAQUS模擬巖土開挖過程最常用的方法是“生死單元法”，通過賦予土體單元“死”的特性來表示外部作用使這一部分土體離開地層。同時對井眼處施加鉆井液液柱壓力，在地應(yīng)力、鉆井液液柱壓力、孔隙壓力作用下，地層會出現(xiàn)一個初始的彈性變形，其應(yīng)力分布如圖4所示，位移分布如圖5所示。同時由于鉆井液液柱壓力大于地層壓力，鉆井液會向地層發(fā)生滲流，使井周孔隙壓力發(fā)生改變。

圖4 井眼鉆開后初始米塞斯應(yīng)力分布Fig.4 Initial Mises stress distribution after drilling

圖5 井眼鉆開后地層初始位移Fig.5 Initial displacement of the formation after drilling

由地層初始應(yīng)力分布和位移分布可以看出，地層的井口會出現(xiàn)一個明顯的應(yīng)力集中，且在最小水平地應(yīng)力方向上應(yīng)力最大，在最大水平地應(yīng)力方向上位移最大，在最小水平地應(yīng)力方向上位移最小，此時井眼將呈現(xiàn)為橢圓形。

3.2 水合物儲層蠕變特征

對于非蠕變地層來說，井壁穩(wěn)定的判斷是一個瞬時或者短時效過程。但對于蠕變地層來說，井眼的變形是一個持續(xù)的過程，隨著蠕變的進(jìn)行井眼將發(fā)生縮徑，同時井眼的形狀也將不斷改變。使用ABAQUS模擬蠕變過程時使用VSICO分析步進(jìn)行分析，井眼打開240 h后其應(yīng)力分布與位移分布分別如圖6、圖7所示。

圖6 蠕變后地層米塞斯應(yīng)力分布Fig.6 Mises stress distribution after creep

圖7 蠕變后地層位移Fig.7 Displacement after creep

由圖6可以看出，經(jīng)240 h蠕變后，地層近井地帶應(yīng)力分布與井眼剛打開時的應(yīng)力分布出現(xiàn)很大變化，選擇井壁處1/4圓周路徑上的點(diǎn)的應(yīng)力繪制在同一坐標(biāo)軸，如圖8所示，起點(diǎn)為最小水平地應(yīng)力方向上的點(diǎn)，既0°方向為最小水平地應(yīng)力方向，90°為最大水平地應(yīng)力方向。

圖8 井壁應(yīng)力分布圖Fig.8 Stress distribution on the wellbore wall

由圖8可以看出，井眼打開初期只發(fā)生彈性變形階段，井壁處出現(xiàn)應(yīng)力集中。井周一周在0°方向應(yīng)力最大，在90°方向應(yīng)力最小，這也符合非均勻地應(yīng)力下直井井眼應(yīng)力應(yīng)力分布情況。此時在0°方向應(yīng)力大小為11.7 MPa，90°方向應(yīng)力大小為10.6 MPa，最大應(yīng)力與最小應(yīng)力比值約為1.1，這與初始設(shè)置的地應(yīng)力比相同。當(dāng)井眼鉆開240 h以后應(yīng)力分布發(fā)生了大的改變，總體上與初始時刻相比井壁處應(yīng)力有所減小，這是由于井眼縮徑所導(dǎo)致，井眼縮徑將減輕應(yīng)力集中程度。在最0°方向應(yīng)力減小值最大，減小為為初始時刻的87%。同時井周應(yīng)力差減小，最大應(yīng)力處任為0°處但最小應(yīng)力分布處為45°處，二者應(yīng)力比值約為1.03小于初始時刻1.1，90°與0°方向應(yīng)力比為1.005 5，井周應(yīng)力分布已近似相等。

對于蠕變地層來說應(yīng)力分布的變化必將導(dǎo)致蠕變速率的變化，從而影響井眼的變形。由圖5和圖7可以看出，蠕變發(fā)生前井眼受非均勻地應(yīng)力作用在最大水平地應(yīng)力出發(fā)生最大變形，在最小水平地應(yīng)力處變形最小，選擇蠕變發(fā)生前后1/4圓周地層變形如圖9所示。

由圖9可以看出，蠕變前后井壁位移發(fā)生變化，但總體趨勢相同。0°方向(最小水平地應(yīng)力方向)最小，90°(最大水平地應(yīng)力方向)最大，但最小位移與最大位移比值卻發(fā)生很大變化。蠕變發(fā)生前，井壁最大位移為1.43 cm，最小位移為1.25 cm，最大位移與最小位移的比值約為1.14；蠕變發(fā)生后井壁最大位移為6.79 cm,最小位移為6.7 cm，最大位移與最小位移的比值為1.01。這表明蠕變時井眼不同方位會以不同的蠕變速率進(jìn)行縮徑，圖10所示為最大水平地應(yīng)力方向和最小水平地應(yīng)力方向井眼縮徑隨時間變化曲線。

圖9 蠕變前后井周位移分布Fig.9 Wellbore displacement distribution before and after creep

圖10 井眼縮徑曲線Fig.10 Wellbore shrinkage curve

圖11為最小水平地應(yīng)力方向與最大水平地應(yīng)力方向蠕變速率曲線，從圖11中可以看出，不同方位井壁上的點(diǎn)初期縮徑速率很大并隨著時間增加而下降，這是由于水合物衰減蠕變的特征決定。圖12所示為不同時刻最小水平地應(yīng)力方向井眼縮徑速率與最大水平地應(yīng)力方向縮徑速率比值，從圖12中可以看出，初始時刻最小地應(yīng)力方向縮徑速率與最大主應(yīng)力方向縮徑速率的比值約為1.17，略大于初始最大水平地應(yīng)力與最小水平地應(yīng)力比值。

圖11 蠕變速率曲線Fig.11 Wellbore shrinkage rate

圖12 0°方向與90°方向縮徑速率比值Fig.12 The ratio of shrinkage rate at the direction of 0° to 90°

這主要是受井眼的應(yīng)力分布及井眼的初始位移所決定。從圖8看出，彈性階段結(jié)束后井壁出現(xiàn)應(yīng)力集中，在最小水平地應(yīng)力方向最大，最大水平地應(yīng)力方向最小，二者比值約為1.1。水合物沉積物的蠕變速率受應(yīng)力影響，在井眼不同方位地層在初始時刻二者所受應(yīng)力比值為1.1，故二者的應(yīng)變速率比值也約為1.1。同時由圖9可知，初始時刻最大水平地應(yīng)力方向井壁由于彈性變形產(chǎn)生初始縮徑量大于最小水平地應(yīng)力方向，因此二者的縮徑速率比值將略大于1.1。

由圖12可以看出，隨著蠕變的進(jìn)行，0°與90°方向的縮徑速率比值不斷減小逐漸趨近于1。這是由于初始時刻0°縮徑量小于90°方向，但其縮徑速率大于90°方向，這將導(dǎo)致0°與90°方向井眼縮徑差距將不斷減小，井眼形狀的改變將影響經(jīng)井眼周圍應(yīng)力分布，減小兩個方向的應(yīng)力差，使得縮徑速率的比值隨時間減小并趨近于1。這也解釋了井眼打開240 h后，井眼周圍最大應(yīng)力與最小應(yīng)力的比值遠(yuǎn)小于蠕變發(fā)生前并趨近于1。由此對于水合物地層的井眼形狀變化如圖13所示。

井眼打開初期在地應(yīng)力作用下井眼由于彈性變形將呈現(xiàn)為橢圓形，長軸為最小水平地應(yīng)力方向，短軸為最大水平地應(yīng)力方向。隨著蠕變的進(jìn)行，井筒趨近于圓。井周應(yīng)力分布初始時刻在最小水平地應(yīng)力方向與最大水平地應(yīng)力方向應(yīng)力集中程度差異越大，隨著蠕變的進(jìn)行，應(yīng)力分布將趨于均勻。

3.3 水合物儲層井眼縮徑影響因素分析

由于實(shí)際地層中水合物層的深度、地應(yīng)力的比值水合物飽和度都影響水合物層的井周應(yīng)力分布，因此建立一系列模型模擬不同因素對于水合物層井眼縮徑規(guī)律的影響。

3.3.1 水合物層深度對井眼縮徑的影響

選擇地應(yīng)力比為1∶1.1、鉆井液密度為1.05 g/cm3、井眼直徑為0.25 m模擬不同深度下水合物層井眼縮徑。圖14所示為不同深度下蠕變240 h后的井眼最大縮徑量。

由圖14可以看出，井眼的縮徑量隨水合物層所在深度的增加而增大，既水合物層所在深度越深，最終的井眼尺寸越小。這是由于地層深度越大，地應(yīng)力的值越大，雖然地應(yīng)力的比值相同，但應(yīng)力差值增大；同時地應(yīng)力與井眼液柱壓力的差值越大從而導(dǎo)致井眼的縮徑量隨深度的增加而增加。

圖14 不同深度井眼縮徑量Fig.14 Wellbore shrinkage of the hydrate reservoir in different depth

3.3.2 地應(yīng)力非均勻性對井眼縮徑的影響

選擇地層深度為500 m、鉆井液密度為1.05 g/cm3、井眼直徑為0.25 m模擬不同地應(yīng)力比條件下水合物層的井眼縮徑。圖15(a)和圖15(b) 分別為不同地應(yīng)力比值下0°和90°蠕變240 h前后的井眼變形量。

圖15 蠕變前后井眼縮徑量Fig.15 Wellbore shrinkage before and after creep

由圖15可以看出，地應(yīng)力比值的增大會增大井眼的縮徑程度。對井眼形狀來說，蠕變發(fā)生前井筒微小變形是由地應(yīng)力的擠壓作用而造成的彈性變形，地應(yīng)力比值越大，造成井口的不同方向變形不同程度越大，井眼向橢圓變化。但地層蠕變的發(fā)生后井眼的變形將不再是持續(xù)原本的彈性變形方向，而是從橢圓向圓發(fā)展，地應(yīng)力差值越大蠕變同一段時間后地層的橢圓程度越大。

圖16所示為初始時刻0°方向井眼縮徑速率與90°方向比值。從圖16可以看出，地應(yīng)力比值越大，井眼不同方向的縮徑速率比值越大，既蠕變初期的井眼收縮的不均勻性越大。一方面是由于地應(yīng)力比值越大，初始時刻最小主應(yīng)力方向應(yīng)力集中程度較最大主應(yīng)力方向越大；同時初始位移差也越大使得蠕變速率的不一致性隨地應(yīng)力的增大而增大。

圖16 不同地應(yīng)力比下的縮徑速率比值Fig.16 Ratio of wellbore shrinkage rate at different in-situ stress ratios

3.3.3 水合物飽和度對眼縮徑的影響

由于水合物飽和度與地層蠕變目前研究較少，但水合物飽和度對地層的其他力學(xué)性質(zhì)的影響已有一個較為清晰的研究。因此，從彈性模量的角度研究水合物飽和度對地層井眼縮徑的影響。實(shí)驗得出水合物沉積物彈性模量和飽和度為0的沉積物彈性模量比值與其飽和度曲線，如圖17所示[31]。

圖17 水合物飽和度與彈性模量比值關(guān)系[31]Fig.17 Relationship between hydrate saturation and elastic modulus ratio[31]

擬合曲線可用三次多項式表示為

E/E0=8.85Sh3+1.54Sh2+1.249 4Sh+1

(9)

式(9)中：E0為水合物飽和度為0時沉積物彈性模量；E為某一飽和度下沉積物彈性模量。

建立水合物飽和度與水合物沉積物彈性模量的關(guān)系，由此模擬水合物層井眼縮徑規(guī)律。圖18所示為水合物飽和度與井眼縮徑量關(guān)系。

圖18 水合物飽和度與井眼縮徑Fig.18 Hydrate saturation and wellbore shrinkage

隨著水合物飽和度的降低井眼的縮徑量逐漸增加，當(dāng)水合物飽和度由30%變?yōu)?0%時，蠕變縮徑量增加約14%。這是由于飽和度的降低水合物層的彈性模量降低，初始彈性變形階段彈性位移增大；其次彈性模量減小會增大井眼的變形能力，使井眼的縮徑速率更快，最終使井眼的縮徑量增加。

4 結(jié)論

(1)井眼打開時，水合物蠕變地層會在非均勻地應(yīng)力作用下發(fā)生彈性變形，在最大水平地應(yīng)力方向位移最大，最小水平地應(yīng)力方向最小。同時井眼應(yīng)力集中出現(xiàn)不均勻現(xiàn)象，在最小水平地應(yīng)力方向應(yīng)力最大，最大水平地應(yīng)力方向最小，二者比值約為地應(yīng)力比值。

(2)井眼發(fā)生蠕變縮徑時，不同方位的縮徑速率不同，最小水平地應(yīng)力方向縮徑速率最大，最大水平地應(yīng)力方向縮徑速率最小。

(3)隨著蠕變的發(fā)展，井眼的縮徑速率隨時間發(fā)生指數(shù)衰減，當(dāng)井眼鉆開時間較長時，最大縮徑速率與最小縮徑速率比直接近1，最大縮徑量與最小縮徑量也接近1，井眼逐漸呈現(xiàn)為圓形等速縮徑。

(4)蠕變發(fā)生一定時間后，井眼的應(yīng)力分布將從不均勻分布向均勻分布發(fā)展，最終將消除地應(yīng)力的不均勻?qū)蹜?yīng)力分布的影響，從巖石受力角度井壁趨于穩(wěn)定。

(5)隨著水合物層所在深度的增加，井眼的蠕變量將增大；井眼的縮徑量將隨著地應(yīng)力比值的增大而增大，同一時刻井眼的縮徑速率、不均勻性也隨地應(yīng)力比值的增大而增大；隨著水合物飽和度的降低，井眼的縮徑量增加。