基于應變軟化模型的地下儲氣庫井筒水泥環完整性研究

馮臏輝,姜婷婷,尚修瑞,張建華

(1.武漢理工大學安全科學與應急管理學院,湖北 武漢 430070;2.武漢理工大學資源與環境工程學院,湖北 武漢 430070)

天然氣作為一種清潔能源,在工業生產、交通運輸和生活起居等領域均有廣泛的應用。根據國家發改委和國家能源局印發的《“十四五”現代能源體系規劃》,到2025年,全國集約化布局的儲存天然氣能力達到550億m3至600億m3,對于我國天然氣行業而言,“十四五”的主要任務之一是加快地下儲氣庫建設[1]。地下儲氣庫具有儲氣規模大、成本低的優點,是我國儲備天然氣首選方式[2]。地下儲氣庫可靠運行的重要前提之一是井筒密封特性良好。因此,保證井筒在不同工況下的密封性,尤其是保證井筒水泥環的完整性對于保障地下儲氣庫平穩運行具有重要意義[3]。

對于地下儲氣庫井筒水泥環完整性失效問題,國內外學者開展了一系列研究。如:Zinkhan等[4]提出地下儲氣庫井筒由套管-水泥環-地層構成;Vrlstad等[5]、舒剛等[6]基于儲氣庫注采井水泥環試驗設備,對儲氣庫注采井水泥環破壞機理與影響因素進行了研究;殷有泉等[7]、趙新波等[8]基于彈性力學,將井筒問題簡化為平面問題給出了非均勻地應力和熱固耦合下套管-水泥環-井筒組合體應力場的解析解;趙效鋒等[9]基于彈塑性力學,將井筒問題簡化為多層組合厚壁筒問題,得出了固井界面微環隙的計算方法。數值模擬方法由于其試驗條件可控且具有便于進行重復試驗的優勢,在地下儲氣庫井筒穩定性研究中得到了廣泛應用,如:Valov等[10]采用完全耦合的線性熱-孔隙-彈性模型,描述了流體壓力和非均質地應力對套管的壓縮作用,并基于有限元方法研究了導致地下儲氣庫井筒水泥環失效的因素;Yuan等[11]通過建立三維地層-水泥環-套管彈性組合模型,研究了不同作業參數對兩膠結面的范式等效應力和組合體總位移的影響;Bu等[12]針對套管-水泥-地層界面易產生微環隙問題,基于ABAQUS有限元軟件研究了套管內壓降低過程中水泥與套管界面的拉應力,以及水泥硬化收縮過程中水泥與地層界面的拉應力。

以上關于地下儲氣庫井筒水泥環完整性的研究中,水泥環本構模型為線彈性或理想彈塑性,然而水泥環在圍壓作用下發生破壞存在塑性形變[13],且有應變軟化現象[14],故在地下儲氣庫井筒完整性研究中運用能反映水泥環破壞全過程的水泥環本構模型對于井筒水泥環密封失效機理研究具有重要的意義。本文基于水泥環力學行為特性,考慮水泥環在圍壓條件受載時的應變軟化特征,建立了水泥環應變軟化本構模型,并通過建立套管-水泥環-地層組合體數值模型,根據地下儲氣庫作業特點,對地下儲氣庫井筒水泥環密封失效原因進行了分析,驗證了所建立的水泥環本構模型的可靠性。

1 考慮峰前硬化階段的巖石應變軟化模型建立

1.1 典型巖石應力-應變曲線

相較于金屬材料,巖石材料內部結構復雜,并且由于不同巖石材料之間存在性質差異,因此巖石材料的力學行為也各不相同。前人研究巖石性質對工程產生的影響時,常將巖石應力-應變曲線簡化為理想彈塑性、理想彈脆性和理想應變軟化3種形式,如圖1所示。圖1(a)中,對于理想彈塑性材料,在彈性階段,材料應力與應變呈線性變化,在塑性階段,材料應力保持不變,而材料應變隨時間推移逐漸增大;圖1(b)中,對于理想彈脆性材料,材料應力達峰值之后迅速跌落,隨后材料進入塑性流動階段;圖1(c)中,對于理想應變軟化材料,與彈脆性材料相比,材料在應力達峰值之后有一個應變軟化階段,隨后再進入殘余階段。

圖1 典型的巖石應力-應變曲線[15]Fig.1 Typical simplified rock stress-strain curves[15]

以上簡化的巖石應力-應變曲線盡管應用廣泛,但也存在不合理之處,即未考慮材料的峰前硬化階段的力學行為。如圖2所示,以典型的巖石全應力-應變曲線為例[16],該典型的巖石全應力-應變曲線一般劃分為5個階段:①OA段,巖石內部微裂隙初始受壓閉合階段;②AB段,巖石線彈性變形階段;③BC段,巖石應變軟化階段;④CD段,巖石應力迅速跌落階段;⑤DE段,巖石殘余應力階段。因此,為了更為恰當地表達巖石的力學行為,應當準確描述巖石應力出現拐點前后即應變軟化階段的巖石應力-應變曲線情況。

圖2 典型的巖石全應力-應變曲線Fig.2 Typical complete stress-strain curve of rocks

1.2 水泥環應變軟化模型建立

假設水泥環在初始屈服后應力狀態均符合莫爾-庫侖強度準則,即:

(1)

式中:σ1和σ3分別為水泥環最大和最小主應力(MPa);c和φ分別為水泥環的黏聚力(MPa)和內摩擦角(°)。

用軸向塑性應變εp作為軟化參數,在水泥環出現屈服后,將其應力看作軸向應變ε的函數,水泥環軟化階段應力-應變關系函數可表示如下[17]:

(2)

其中:

(3)

式中:ε為水泥環軸向應變。

如圖3所示,將水泥環黏聚力與軟化參數之間的函數關系根據巖石應力-應變曲線劃分為多段并進行分析,以描述水泥環黏聚力的演化規律。

對于水泥環初始屈服至強度峰值軟化階段任意微元段,有:

E0(dε-dεp)=Δσ=σ′(ε)dε

(4)

式中:E0為水泥環初始屈服前的彈性模量(MPa);σ′(ε)為在軸向應變ε處水泥環的應力-應變曲線的斜率。

水泥環初始屈服前的彈性模量E0可用初始屈服點與原點連線的斜率來表示:

(5)

式中:σ0為水泥環初始屈服強度(MPa);ε0為水泥環初始屈服強度處對應的軸向應變。

因此,對于峰前應變硬化階段,三軸壓縮過程中水泥環軸向應變ε1與軸向塑性應變εp之間的關系可表示為

(6)

同理,對于峰后應變軟化階段,有:

(7)

將巖石應力-應變全過程曲線劃分為4階段,如圖4所示。

圖4 4階段巖石應變軟化模型[18]Fig.4 4-stage strain softening model of rock[18]

對于圖4中表達的巖石應變軟化模型,峰前應變硬化階段巖石應力-應變曲線關系式為

(8)

式中:σM、σ0分別為水泥環峰值強度和初始屈服強度(MPa);εM、ε0分別為水泥環峰值強度和初始屈服強度處對應的軸向應變。

聯立上式,可得峰前應變硬化階段水泥環黏聚力的演化規律為

(9)

式中:η為軟化參數。

同理,可得峰后應變軟化階段水泥環黏聚力的演化規律為

(10)

由上兩式可知:水泥環應力-應變曲線呈線性階段時,水泥環黏聚力便與軟化參數呈線性變化,而水泥環應力-應變曲線呈非線性階段時,可將其處理為若干直線段聯合以代替該曲線。

2 模型驗證

2.1 模型計算參數確定

為了證明本文所建立的考慮峰前硬化階段巖石應變軟化模型的合理性,需要基于巖石強度參數的演化規律,對巖石應力-應變全過程曲線通過數值模擬試驗進行驗證。文獻[19]給出了水泥環三軸壓縮試驗資料,包括水泥環的彈性模量、泊松比以及不同圍壓下三軸壓縮應力-應變曲線,并對三軸壓縮試驗曲線進行簡化,見圖5。

圖5 不同圍壓下水泥環三軸壓縮試驗應力-應變曲線Fig.5 Stress-strain curves of the triaxial compression test of the cement ring under different confining pressures

設簡化試驗曲線每一構成部分的兩端點坐標分別為(ε′,σ′)、(ε″,σ″),且初始屈服點坐標為(ε0,σ0)。由式(5)計算得到不同圍壓下水泥環初始屈服前的彈性模量E0,見表1。

表1 不同圍壓下水泥環初始屈服前的彈性模量E0

對于水泥環初始屈服至峰前硬化階段,由式(6)計算得到水泥環軸向塑性應變εp,并由下式計算軟化參數η:

(11)

式中:ε″和ε′分別為簡化的各段水泥環應力-應變曲線結束時和開始時的軸向應變;σ″和σ′分別為簡化的各段水泥環應力-應變曲線結束時和開始時的軸向應力(MPa);η為軟化參數。

由式(2)計算得到峰前硬化階段水泥環的黏聚力值。

同理,可計算出峰后軟化階段水泥環的黏聚力。將上述計算結果匯總于表2。

表2 不同圍壓和軟化狀態下水泥環的黏聚力值

2.2 水泥環試樣三軸壓縮數值模擬

本文利用FLAC3D軟件建立了標準圓柱體水泥環試樣數值模型(直徑為50 mm,高為100 mm),模型材料密度為1 900 kg/m3,彈性模量為5.99 GPa,泊松比為0.20,內聚力為5.30 MPa,內摩擦角為24.5°,抗拉強度為3.61 MPa,對圓柱形網格模型頂、底面用位移控制軸向加載,側面施加垂直壓應力以模擬圍壓,所建立的數值模型如圖6所示。結合表2中的數據,利用建立的水泥環應變軟化模型,對圍壓分別為0、5、10、15和20 MPa條件下的水泥環試樣進行了三軸壓縮數值模擬,并將水泥環三軸壓縮試驗應力-應變曲線與數值模擬曲線進行了對比,其結果見圖7。

圖6 標準圓柱體水泥環試樣FLAC3D數值模型示意圖Fig.6 Schematic diagram of FLAC3D numerical model for standard cylindrical cement ring sample

圖7 不同圍壓下水泥環三軸壓縮試驗應力-應變曲線 與數值模擬曲線的對比Fig.7 Comparison between stress-strain curves of cement ring triaxial compression test and numerical simula- tion curves under different confining pressures

由圖7可知:采用水泥環應變軟化模型計算得到的水泥環應力-應變模擬曲線經歷了彈性變形階段、峰前硬化階段、峰后軟化階段和殘余階段,隨著圍壓的增大,水泥環峰值強度不斷增加,模擬曲線變化趨勢與試驗曲線基本一致;模擬曲線彈性階段在初始屈服點結束,初始屈服點與原點連線的斜率對應于水泥環發生屈服前的彈性模量,這與本文假設相符;此外,對于低圍壓下水泥環應力-應變曲線,至峰值強度后,水泥環應力跌落,以殘余強度維持一定形態,而在高圍壓下水泥環應力-應變曲線經歷了峰值強度,也不會出現急劇的應力跌落現象,而是表現為延性,其強度緩慢降低。

文獻[20]指出邊坡巖土體剪應變增量最大處相對于其他位置最容易發生破壞變形。因此,可以通過巖石最大剪應變增量的變化情況來判斷模型對應力作用的敏感程度。為了確定圍壓對水泥環試樣破壞程度的影響,采用莫爾-庫侖本構模型和應變軟化本構模型,經相同的位移加載計算得到的試樣最大剪應變增量變化曲線,如圖8所示。

圖8 不同圍壓下水泥環試樣內部各點最大剪應變增量變化曲線Fig.8 Curves of maximum shear strain increment at each point inside the cement ring sample under different confining pressures

由圖8可知:采用應變軟化模型計算得到的水泥環試樣最大剪應變增量均大于采用莫爾-庫侖模型計算值,這是因為采用應變軟化模型時,水泥環試樣對應力作用更敏感,更易發生破壞;通過對比兩組曲線的峰值發現,在同等的加載條件下,采用兩種本構模型計算得到的水泥環試樣最大剪應變增量的差值隨圍壓的增加而減少,這主要是由于采用應變軟化模型時水泥環試樣產生的最大剪應變增量減小,并在圍壓的作用下,水泥環試樣的強度弱化速度減慢,抵抗變形的能力增強。由此可見,對地下儲氣庫井筒水泥環完整性進行評價采用應變軟化模型更為合適。

3 地下儲氣庫井筒水泥環完整性分析

3.1 數值模型建立

本研究先使用ANSYS軟件建立地下儲氣庫井筒模型并進行網格劃分,模型長、寬和高分別為1 m、1 m和0.05 m,模型單元為60 420個、節點為73 338個;然后將所建立的地下儲氣庫井筒模型通過軟件接口導入FLAC3D軟件中,生成FLAC3D井筒網格模型,如圖9(a)所示,其中內部第一個圓環表示套管,最外層表示地層,套管與地層之間的圓環部分為水泥環。對構成井筒的套管賦予線彈性模型,地層設為泥巖且采用莫爾-庫侖模型[21],并分別采用莫爾-庫侖模型和應變軟化模型對水泥環的最大剪應變增量、最大主應力和內部位移進行了模擬計算與分析,模型材料參數見表3。水泥環與套管和井筒之間的接觸面采用FLAC3D軟件提供的分界面單元模擬[22],分界面剛度10倍于周邊單元體的最大等效剛度,其計算公式如下:

表3 模型材料參數

圖9 地下儲氣庫井筒網格模型示意圖Fig.9 Wellbore grid model of underground gas storage

(12)

式中:Kn為分界面法向剛度(GPa);Ks為分界面切向剛度(GPa);K為體積模量(GPa);G為剪切模量(GPa);ΔZmin為法向上與分界面單元體相接觸的網格最小寬度(m)。

套管-水泥環和水泥環-地層的接觸面參數見表4,所建立的井筒接觸面模型見圖9(b)。

表4 水泥環接觸面參數

模型邊界條件設置如下:

1) 套管、水泥環、地層的下表面及地層四周固定三向位移。

2) 根據上覆巖層的平均容重和深度,對前后左右四個面分別施加20 MPa的垂直應力以代表水平地應力,對頂面施加21 MPa的垂直應力代表上覆巖層重力[19]。

3) 對套管內表面施加一個通過FISH函數添加的使用正弦函數作為返回值的呈周期變化的法向應力來作為套管循環載荷。

在模擬計算時,首先考慮到巖層自重應力的影響,進行初始地應力的平衡;然后用空模型開挖出井筒;最后計算在套管內壓循環荷載作用下,水泥環位移、應力場的變化。

在地下儲氣庫注采循環作業中,井筒內部會受到交替的高壓和低壓作用。在注氣階段,水泥環受到的應力增大并且產生變形;而在采氣階段,水泥環受到的應力減小但僅有彈性變形恢復,即水泥環經加載后再卸載有塑性變形存在[23],并且在持續交替壓力下,水泥環累積塑性變形不斷增加。因此,水泥環最大剪應變增量、最大主應力和徑向位移會隨地下鹽穴儲氣庫注采作業次數發生變化,其中:

1) 最大剪應變增量集中處反映水泥環發生剪切破壞的可能位置,FLAC3D軟件中定義應變增量張量是與節點位移相關的物理量,在一個無限小的時間dt內,單元應變增量計算公式為[20]

(13)

2) FLAC3D軟件中正值表示拉應力,最大主應力反映水泥環受拉情況,單元最大主應力計算公式為[22]

(14)

式中:J2為偏應力張量第二不變量;θσ為應力洛德角(°);σm為平均應力(MPa)。

3) 水泥環徑向位移反映水泥環內部發生形變的程度,為水泥環徑向各點在受力產生形變后所在位置與原始位置坐標點的差值。

3.2 模擬計算結果與分析

地下鹽穴儲氣庫壽命一般為30～50年[24],為了衡量地下儲氣庫整個壽命周期內井筒水泥環的完整性,將井筒注采氣次數設為15、30、45和60次。

3.2.1 水泥環最大剪應變增量分析

本研究分別采用莫爾-庫侖本構模型和應變軟化本構模型模擬計算得到不同井筒注采氣次數下水泥環的最大剪應變增量云圖,見圖10和圖11所示。

圖10 不同井筒注采氣次數下水泥環最大剪應變增量 云圖(莫爾-庫侖本構模型)Fig.10 Contour of maximum shear strain increment of cement sheath under different gas injection and production times (Mohr-Coulomb constitutive model)

圖11 不同井筒注采氣次數下水泥環最大剪應變增量 云圖(應變軟化本構模型)Fig.11 Contour of maximum shear strain increment of cement sheath under different gas injection and production times (strain softening constitutive model)

由圖10和圖11可知:

1) 隨著注采次數的增加,水泥環內側發生的最大剪應變增量呈增加趨勢。

2) 井筒注采次數相同時,采用莫爾-庫侖本構模型模擬計算得到的水泥環最大剪應變增量小于采用應變軟化本構模型的計算值;在經歷15次注采循環后,采用上述兩種本構模型計算得到的水泥環最大剪應變增量分別為0.146和0.194;在經歷60次注采循環后,采用上述兩種本構模型計算得到的水泥環最大剪應變增量分別為0.217和0.263。

3.2.2 水泥環最大主應力分析

本研究采用莫爾-庫侖本構模型和應變軟化本構模型計算得到的不同注采氣次數下水泥環模擬最大主應力云圖,見圖12和圖13。

圖12 不同井筒注采氣次數下水泥環最大主應力云圖 (莫爾-庫侖本構模型)Fig.12 Contour of maximum principal stress of cement sheath under different gas injection and production times (Mohr-Coulomb constitutive model)

圖13 不同井筒注采氣次數下水泥環最大主應力云圖 (應變軟化本構模型)Fig.13 Contour of maximum principal stress of cement sheath under different gas injection and production times (strain softening constitutive model)

由圖12和圖13可知:

1) 對于兩種本構模型的計算結果,最大主應力的最大值位置均在水泥環內側,且均為壓應力。

2) 水泥環內側有拉應力集中現象,拉應力集中區域主要呈兩點分布且隨著循環加卸載次數的增加而擴大,但均小于水泥環的抗拉強度(圖12),即此時水泥環不會出現徑向裂隙。

3) 最大主應力集中區域呈多點狀均勻分布在水泥環內側,隨著循環加卸載次數的增加,該區域逐漸擴大并有連結成環的趨勢,且有最大主應力集中區域的應力值大于水泥環的抗拉強度(圖13),對于井筒則會出現套管-水泥環分離現象,進而使水泥環密封失效。

3.2.3 水泥環徑向水平位移分析

在完成30次注采氣循環后,水泥環徑向不同位置處的水平位移圖,見圖14。

圖14 水泥環徑向不同位置處的水平位移Fig.14 Radial horizontal displacement of cement sheath

由圖14可以看出:在同等條件下,采用應變軟化本構模型計算得到的水泥環徑向水平位移普遍大于采用莫爾-庫侖本構模型的計算值,且采用應變軟化本構模型計算得到的水泥環徑向水平位移在增加至峰值后維持在一定水平,其變化細微,而采用莫爾-庫侖本構模型計算得到的水泥環徑向水平位移增加至峰值后逐漸減小;此外,首個水泥環徑向監測點位于套管與水泥環交界面,該監測點與其鄰近監測點水泥環徑向水平位移差值明顯大于其他相鄰監測點差值。

4 結 論

1) 本文基于莫爾-庫侖強度準則,考慮峰前硬化階段和峰后軟化階段對水泥環應變軟化的影響,選取軸向塑性應變為軟化參數,獲取了水泥環黏聚力隨軸向塑性應變的變化規律,以此建立了水泥環應變軟化本構模型,模擬結果表明:運用該模型模擬得到的水泥環應力-應變曲線與試驗曲線一致程度較高,能準確描述水泥環峰前硬化階段和峰后軟化階段的力學行為。

2) 對水泥環試樣發生最大剪應力增量的分析結果表明:水泥環的抗壓強度受圍壓的影響,相比于莫爾-庫侖本構模型,應變軟化本構模型更能表達出水泥環的強度參數隨應力增加而發生軟化的特性。

3) 建立了套管-水泥環-地層組合體數值模型,對水泥環賦予莫爾-庫侖本構模型和應變軟化本構模型,模擬分析了地下儲氣庫井筒水泥環的完整性,結果表明:相比于采用莫爾-庫侖本構模型的計算結果,采用應變軟化本構模型時,水泥環內部產生的徑向水平位移較大,水泥環所受的應力也更大,此時水泥環更易發生塑性變形,產生微環隙。