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

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

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

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

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

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

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

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

相較于金屬材料,巖石材料內部結構復雜,并且由于不同巖石材料之間存在性質差異,因此巖石材料的力學行為也各不相同。前人研究巖石性質對工程產(chǎn)生的影響時,常將巖石應力-應變曲線簡化為理想彈塑性、理想彈脆性和理想應變軟化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段,巖石殘余應力階段。因此,為了更為恰當?shù)乇磉_巖石的力學行為,應當準確描述巖石應力出現(xiàn)拐點前后即應變軟化階段的巖石應力-應變曲線情況。

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

1.2 水泥環(huán)應變軟化模型建立

假設水泥環(huán)在初始屈服后應力狀態(tài)均符合莫爾-庫侖強度準則,即:

(1)

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

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

(2)

其中:

(3)

式中:ε為水泥環(huán)軸向應變。

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

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

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

(4)

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

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

(5)

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

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

(6)

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

(7)

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

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

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

(8)

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

聯(lián)立上式,可得峰前應變硬化階段水泥環(huán)黏聚力的演化規(guī)律為

(9)

式中:η為軟化參數(shù)。

同理,可得峰后應變軟化階段水泥環(huán)黏聚力的演化規(guī)律為

(10)

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

2 模型驗證

2.1 模型計算參數(shù)確定

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

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

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

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

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

(11)

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

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

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

表2 不同圍壓和軟化狀態(tài)下水泥環(huán)的黏聚力值

2.2 水泥環(huán)試樣三軸壓縮數(shù)值模擬

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

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

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

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

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

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

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

3 地下儲氣庫井筒水泥環(huán)完整性分析

3.1 數(shù)值模型建立

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

表3 模型材料參數(shù)

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

(12)

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

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

表4 水泥環(huán)接觸面參數(shù)

模型邊界條件設置如下:

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

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

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

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

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

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

(13)

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

(14)

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

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

3.2 模擬計算結果與分析

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

3.2.1 水泥環(huán)最大剪應變增量分析

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

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

圖11 不同井筒注采氣次數(shù)下水泥環(huán)最大剪應變增量 云圖(應變軟化本構模型)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) 隨著注采次數(shù)的增加,水泥環(huán)內側發(fā)生的最大剪應變增量呈增加趨勢。

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

3.2.2 水泥環(huán)最大主應力分析

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

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

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

由圖12和圖13可知:

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

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

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

3.2.3 水泥環(huán)徑向水平位移分析

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

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

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

4 結 論

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

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

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