天然氣水合物藏井眼蠕變縮徑研究

李 陽， 程遠方, 閆傳梁， 周曉輝

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

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

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

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

1 數學模型

1.1 流固耦合模型

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

σ′=σ-αpp

(1)

式(1)中：σ′ 為有效應力；σ為地層中任一點的總應力；α為有效應力系數；pp為孔隙壓力。

由于實際的多孔介質微觀組成復雜，很難用精確的數學關系將流體固體分開，因此所建立的模型為理想的連續系統，忽略巖石孔隙的微觀結構變化，但通過一些特殊變量的變化可以表征孔隙流體流動的過程，以此模擬流固耦合的過程。

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

(2)

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

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

(3)

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

1.2 水合物儲層蠕變模型

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

ε(t)=εe+εc

(4)

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

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

(5)

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

(6)

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

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

(7)

1.3 水合物儲層蠕變模型

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

(8)

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

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

2 數學模型

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

2.1 水合物儲層井眼模型

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

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

2.2 水合物儲層井眼模型

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

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

2.3 地應力平衡

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

圖3 地應力平衡Fig.3 Initial in-situ stress equilibrium

3 水合物層井眼縮徑規律

3.1 井眼鉆開時初始變形

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

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

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

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

3.2 水合物儲層蠕變特征

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

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

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

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

圖8 井壁應力分布圖Fig.8 Stress distribution on the wellbore wall

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3.2 地應力非均勻性對井眼縮徑的影響

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

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

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

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

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

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

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

圖17 水合物飽和度與彈性模量比值關系[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為某一飽和度下沉積物彈性模量。

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

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

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

4 結論

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

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

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

(4)蠕變發生一定時間后，井眼的應力分布將從不均勻分布向均勻分布發展，最終將消除地應力的不均勻對井眼應力分布的影響，從巖石受力角度井壁趨于穩定。

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