安國印 王斌 王超 陳春花 路崢 史志峰 王同濤

1. 華北石油管理局有限公司江蘇儲氣庫分公司；2. 中石化川氣東送天然氣管道有限公司；3. 中國科學(xué)院武漢巖土力學(xué)研究所

金壇鹽穴儲氣庫作為我國第一座鹽穴儲氣庫，從選址到建成運行已經(jīng)經(jīng)歷了20 余年。目前已經(jīng)累計建成鹽穴單腔40 余個，形成工作氣量約12 億m3。鹽穴建造原理：通過井眼向鹽巖地層中注入淡水溶解鹽巖地層然后再將鹵水排出形成地下空間［1］。以金壇為例，鹽穴儲氣庫的體積一般可以達到十幾萬到幾十萬方［2］，其高度可以達到150 m、直徑達到80 m 左右。因此，確保這么大的地下空間結(jié)構(gòu)的穩(wěn)定性是鹽穴儲氣庫建設(shè)和運行過程中需要解決的關(guān)鍵問題。2015 年聲吶測腔結(jié)果表明金壇B 鹽穴頂部發(fā)生了較大范圍的頂部失穩(wěn)破壞，失穩(wěn)破壞區(qū)域長度約為50 m、厚度2～4 m［3］。該次失穩(wěn)破壞嚴重威脅到了套管鞋位置處的密封安全，如果不采取措施存在再次發(fā)生失穩(wěn)破壞的可能性。由于鹽穴儲氣庫失穩(wěn)破壞既表現(xiàn)出很強的塑性破壞特征［4-6］(腔體體積收縮、底部隆起和腔壁大變形)，又表現(xiàn)出很強的脆性破壞特征(頂部脫落、片幫和開裂)，如何準確、定量地對鹽穴儲氣庫的穩(wěn)定性進行評價面臨諸多挑戰(zhàn)。同時，鹽穴儲氣庫穩(wěn)定性受到的影響因素較多［7-10］，包括：腔體幾何參數(shù)、地層結(jié)構(gòu)參數(shù)和運行參數(shù)等，且這些參數(shù)存在一定的耦合關(guān)系。

為了確保鹽穴儲氣庫的運行安全，一般會利用帶壓聲吶測腔技術(shù)對鹽穴儲氣庫的腔體形狀進行定時測量，并根據(jù)測量結(jié)果對腔體的穩(wěn)定性進行評價，進而對腔體運行參數(shù)進行優(yōu)化。2019 年，在對金壇A 鹽穴儲氣庫進行聲吶測腔時，發(fā)現(xiàn)與2017 年聲吶測腔結(jié)果相比，該腔體體積收縮率達到了3.86%，遠大于金壇鹽穴儲氣庫其他腔體體積收縮率。為了對該腔體體積收縮率偏大的原因進行分析，進而對該鹽穴儲氣庫運行參數(shù)進行優(yōu)化，擬開展以下工作：建立該腔體的三維地質(zhì)力學(xué)模型并對其2017—2019年的運行工況進行模擬；對該腔體的穩(wěn)定性進行評價，找出導(dǎo)致腔體體積收縮率過大的原因。

1 A 鹽穴儲氣庫背景介紹

金壇A 鹽穴儲氣庫于2005 年9 月開始鉆井，同年10 月完鉆，累計井深1 198 m，其中?244.5 mm完井套管柱下深為1 012 m。2007 年10 月開始造腔作業(yè)，2015 年12 月完成造腔。2015 年12 月下入聲吶對該溶腔進行聲吶測腔，發(fā)現(xiàn)該腔體在上部明顯出現(xiàn)縮徑并且形成了大量懸挑結(jié)構(gòu)，威脅到后期注氣排鹵和運行安全。通過利用天然氣阻溶技術(shù)對該懸挑部分進行修復(fù)，并于2017 年5 月對修復(fù)完成的A 溶腔進行注氣排鹵前的最后一次聲吶測腔。該次聲吶測腔結(jié)果表明：該腔最大半徑為57.5 m，位于深度1 072.0 m、方位160～340°處；最大直徑為73.8 m，位于深度1 124.0 m、方位175～355°處。腔體總體積約為19.96×104m3，腔頂?shù)茁裆罘謩e為1 026.8 m 和1 130.3 m。2017 年7 月，該腔開始注氣排鹵作業(yè)并于2017 年11 月完成注氣排鹵，正式投產(chǎn)運行。為了對該腔的運行安全進行評價，2019 年4 月進行了帶壓測腔作業(yè)，測腔結(jié)果表明腔體總體積約為19.19×104m3，可知體積收縮率約為3.86%，明顯偏大。根據(jù)金壇鹽穴儲氣庫前期聲吶測腔結(jié)果可知，對于類似腔體運行1.6 a 左右體積收縮率一般都小于1%［11］。圖1 給出了A 溶腔2017 年和2019 年聲吶測腔結(jié)果。通過對比圖1 中2017 年和2019 年A 腔聲吶測腔數(shù)據(jù)可以發(fā)現(xiàn)由于A 腔體形狀較為復(fù)雜，初步認為2019 年聲吶測腔時這些位于腔體上部位置和腔底的空間未能被探測到是導(dǎo)致體積收縮率測量值偏大的主要因素。

根據(jù)前期鉆井資料可知A 腔所在含鹽系地層頂部埋深約為996 m、底部埋深約為1 195 m，整個含鹽系地層厚度為199 m，其中厚度大于2 m、泥質(zhì)含量大于99%的夾層有4 個，它們的頂?shù)茁裆罘謩e為：1 077.6～1 079.6 m、1 099.8～1 102.2 m、1 160.6～1 163.2 m 和1 186.0～1 188.4 m。從圖1 中可以看出A 腔體形狀較為不規(guī)則，腔壁上形成了很多懸挑塊體結(jié)構(gòu)。A 腔在不同方位上的半徑分布也表明其在不同方位上腔體半徑相差較為懸殊。同時，圖1也表明A 腔頂部較為平直且尺寸較大，最大位置處直徑可以達到60 m 左右，存在頂部失穩(wěn)掉塊的可能性。

圖1 A 腔2017 年和2019 年聲吶測腔結(jié)果Fig. 1 Sonar cavity test results of Cavity A in 2017 and 2019

2 三維地質(zhì)力學(xué)模型及邊界條件

根據(jù)圖1 中的聲吶測腔結(jié)果可知A 腔體形狀較為復(fù)雜，如果在建立三維地質(zhì)力學(xué)模型時完全按照聲吶測腔實際數(shù)據(jù)進行建模將會導(dǎo)致模型網(wǎng)格數(shù)過多、計算效率低、甚至計算結(jié)果不收斂。因此，本文在三維地質(zhì)力學(xué)建模過程中采用平均半徑尺寸進行建模，這樣既可以保證腔體體積基本與原測量值一致(誤差小于2%)、腔體形狀基本與原測量形狀相似，又可以消除腔壁不規(guī)則對計算結(jié)果的影響。圖2給出了A 腔三維地質(zhì)力學(xué)模型及其邊界條件。該模型為一個長方體，其長、寬、高分別為800 m、400 m和700 m，A 溶腔位于模型中心，腔體高度約為100 m、最大直徑約為60 m(y方向的半徑為30 m)?？梢杂嬎愕玫紸 腔體邊界與模型邊界之間的最小距離約為5～8 倍A 腔體在相應(yīng)方向上的尺寸，能夠?qū)⑦吔缧?yīng)對計算結(jié)果的影響降到最低。考慮到模型的對稱性，本文建立了1/2 三維地質(zhì)力學(xué)模型以提高計算效率。該地質(zhì)力學(xué)模型包括泥巖層、鹽巖層和夾層。其巖石力學(xué)參數(shù)取值如表1 所示。

金壇鹽巖蠕變特征符合Norton-Hoff 穩(wěn)態(tài)蠕變模型，其數(shù)學(xué)表達式為

圖2 A 腔三維地質(zhì)力學(xué)模型及其邊界條件Fig. 2 3D geomechanics model of Cavity A and its boundary conditions

表1 金壇鹽巖巖石力學(xué)性能參數(shù)Table 1 Rock mechanical parameters of Jintan salt rock

模型中一共包括3 個夾層，它們頂部埋深分別為1 077.6 m、1 099.8 m 和1 160.6 m，厚度分別為2 m、2.4 m 和2.6 m，其他夾層由于距離腔體距離較遠在本模型中未考慮。模型底部作用有固定約束條件，即計算過程中模型底部不發(fā)生任何方向上的位移。模型頂部作用有上覆巖層壓力，本次模擬計算中將埋深750 m 以上的地層簡化為相應(yīng)的上覆巖層壓力(σv)作用在模型頂部。金壇地區(qū)平均上覆巖層壓力梯度約為23 kPa/m，可以計算得到上覆巖層壓力σv為17.25 MPa(圖2)。模型的4 個垂直面上作用有水平方向上的位移約束，限制計算過程中模型在水平方向上發(fā)生位移，而垂直方向上可以自由移動。根據(jù)金壇地區(qū)前期地應(yīng)力測試結(jié)果可知：鹽巖地層中同一深度各個方向上的初始地應(yīng)力和地應(yīng)力梯度相等。按照這一原則對模型施加初始地應(yīng)力條件。腔體內(nèi)施加有內(nèi)壓，由于主要對A 溶腔投產(chǎn)后的穩(wěn)定性進行評價，因此內(nèi)壓按照實際內(nèi)壓監(jiān)測值選取。圖3 給出了2017 年投產(chǎn)至2019 年聲吶測腔時間范圍內(nèi)內(nèi)壓監(jiān)測值和數(shù)值模擬擬合值變化曲線。該內(nèi)壓是套管鞋深度位置處的壓力，是根據(jù)井口壓力監(jiān)測值計算得到。在數(shù)值模擬計算過程中采用分段函數(shù)對套管鞋位置處壓力監(jiān)測值進行擬合，然后施加到腔體內(nèi)壁上。模擬計算時間從2017 年11 月至2019 年4 月，約1.6 a。

圖3 A 腔套管鞋位置處壓力隨時間變化Fig. 3 Variation of the pressure of Cavity A at the casing shoe over the time

為了消除邊界效應(yīng)對于計算結(jié)果的影響，整個計算模型尺寸比較大，而研究重點是腔體周圍的巖體，因此在網(wǎng)格劃分時采用放射性網(wǎng)格，即：距離腔體越近網(wǎng)格尺寸越小，距離越遠尺寸越大。這樣既能夠保證研究區(qū)域具有較高的計算精度，又能夠提高計算效率。在模型網(wǎng)格劃分時使用到了四面體網(wǎng)格、六面體網(wǎng)格和金字塔型網(wǎng)格。為了避免網(wǎng)格質(zhì)量對計算結(jié)果的影響，在計算之前對所有的網(wǎng)格質(zhì)量進行檢測，消除低質(zhì)量的網(wǎng)格對計算結(jié)果的影響。該模型一共包括379 671 個節(jié)點、65 705 個單元。考慮到FLAC3D在計算巖土結(jié)構(gòu)大變形等問題中具有較為顯著的優(yōu)勢，將計算模型導(dǎo)入到FLAC3D進行計算，利用Tecplot 對計算結(jié)果進行后處理。

3 穩(wěn)定性評價結(jié)果分析及運行參數(shù)優(yōu)化

3.1 體積收縮率

圖4 給出了A 腔從2017 年11 月開始注氣投產(chǎn)到2019 年4 月帶壓測腔這一時間段體積收縮率計算值隨時間變化關(guān)系曲線。體積收縮率定義為鹽穴體積減少量與2017 年聲吶測腔體積之比，是衡量鹽穴儲氣庫穩(wěn)定性和可用性的重要指標之一。體積收縮率隨著運行時間逐漸增加，體積收縮率增加速率隨著運行壓力增加而降低。計算結(jié)果表明：A 腔運行1.6 a 后體積收縮率約為0.54%，是一個較為合理的數(shù)值。而2019 年聲吶測腔數(shù)據(jù)表明：與2017 年的相比，體積收縮率約為3.86%。通過對比圖1 中的聲吶測腔結(jié)果初步推斷造成A 腔體積收縮率較大的原因為測量誤差，而不是腔體真實收縮造成的。2019 年聲吶測腔數(shù)據(jù)與2017 年相比，腔體上部和底部區(qū)域出現(xiàn)較為明顯變形，即測量不到2017 年測到的那些空間。而其他位置(埋深1 040～1 120 m)處腔體基本上沒有變形，因此排除了A 溶腔只在腔體上部和下部變形的可能性。

圖4 體積收縮率隨時間的變化Fig. 4 Variation of volume shrinkage over the time

3.2 變形量

圖5 給出了運行1.6 a 后A 鹽穴圍巖中變形量分布云圖。計算結(jié)果表明在當前運行壓力條件下運行1.6 a 后腔壁整體變形量較小，一般在10 cm 以內(nèi)。這與2019 年聲吶測腔數(shù)據(jù)較為一致(不考慮腔體上部和底部未檢測到的部位)，說明本文建立的模型具有較高的計算精度。從變形量角度，A 鹽穴儲氣庫具有較好的穩(wěn)定性。從圖5a 中可以看出總變形量(DISP)較大區(qū)域出現(xiàn)在A 鹽穴腔體中下部，主要因為該部分直徑較大，受到的載荷也較大。圖5b中計算結(jié)果表明垂直變形量(ZDISP)在腔頂部向下、底部向上，這與實際情況相符。同時，可以看出溶腔頂部向下垂直變形量最大，因此腔頂套管鞋位置處套管安全是受到垂向變形控制的。

圖5 運行1.6 a 后A 腔圍巖中變形量分布云圖Fig. 5 Distribution contour of deformation in the surrounding rock of Cavity A after operation of 1.6 a

3.3 塑性區(qū)

Mohr—Coulomb 是巖土工程中最為常用的準則之一，具有物理意義明確、表達式簡單和參數(shù)容易利用室內(nèi)實驗獲取等優(yōu)點，在巖土工程中得到了較為廣泛的應(yīng)用。圖6 給出了利用Mohr-Coulomb 準則獲得的A 鹽穴圍巖中塑性區(qū)分布云圖?？芍苄云茐膮^(qū)域主要集中在腔頂、腔體上部的懸挑部位和腔底。頂部較為平直，跨度較大(約30 m，建模時按照評價值選取)，在地應(yīng)力外擠載荷作用下容易發(fā)生局部破壞。當局部破壞發(fā)生后，在重力作用下局部破壞區(qū)域?qū)那豁斆撀湎聛?，嚴重時可能會威脅到套管鞋的密封安全，后期需要重點監(jiān)控。腔體上部的懸挑部位體積較大，也存在失效破壞的可能性。由于該懸挑部分體積較大，發(fā)生破壞時可能會造成腔體上部應(yīng)力場重分布，造成腔頂大面積失穩(wěn)。因此，在后期運行過程中建議加強監(jiān)控，預(yù)防鹽穴頂部局部破壞與上部懸挑結(jié)構(gòu)失穩(wěn)破壞之間連帶破壞，以保證整個腔體的安全。而鹽穴底部塑性區(qū)發(fā)生破壞的可能性較小，因為在自重作用下以及底部結(jié)構(gòu)對承載較為有利，塑性區(qū)不會發(fā)生脫落和片幫等破壞。

圖6 運行1.6 a 后利用Mohr-Coulomb 準則計算得到A 腔圍巖中塑性區(qū)分布云圖Fig. 6 Distribution contour of plastic zone in the surrounding rock of Cavity A after operation of 1.6 a calculated according to Mohr-Coulomb criteria

3.4 等效應(yīng)變

等效應(yīng)變作為衡量鹽巖蠕變變形破壞的評價指標之一，其定義、計算公式和安全閾值見文獻［8,11］。圖7 給出了運行1.6 a 后A 腔圍巖中等效應(yīng)變(ES)分布云圖。A 腔圍巖中等效應(yīng)變較小，遠小于3%的臨界值。從等效應(yīng)變角度，A 腔具有較好的穩(wěn)定性。計算結(jié)果也表明：2019 年聲吶測腔過程中未能夠精確獲得腔體形狀是導(dǎo)致體積收縮率較大的主要因素，而不是實際腔體體積收縮導(dǎo)致的。

3.5 剪脹安全系數(shù)

圖7 運行1.6 a 后A 腔圍巖中等效應(yīng)變(ES)分布云圖Fig. 7 Distribution contour of equivalent strain (ES) in the surrounding rock of Cavity A after operation of 1.6 a

鹽巖在破壞過程中表現(xiàn)出較為典型的剪脹破壞特征，利用剪脹安全系數(shù)對鹽穴儲氣庫穩(wěn)定性進行評價已經(jīng)在工程領(lǐng)域得到了較為廣泛的推廣和應(yīng)用。剪脹安全系數(shù)定義和計算方法見文獻［8,11］。圖8 給出了運行1.6 a 后A 腔圍巖中剪脹安全系數(shù)(SF)分布云圖。從圖8 中可以看出A 腔圍巖中整體安全系數(shù)較高，絕大部分區(qū)域的安全系數(shù)大于1.0，腔體上部懸挑部位存在剪脹安全系數(shù)小于1.0的區(qū)域，是可能發(fā)生失效破壞的危險區(qū)域。同時，計算結(jié)果還表明夾層對局部區(qū)域的剪脹安全系數(shù)影響較大，但是整體安全系數(shù)較高。從剪脹安全系數(shù)角度，在目前運行條件下A 腔具有較好的穩(wěn)定性。

圖8 運行1.6 a 后A 腔圍巖中剪脹安全系數(shù)(SF)分布云圖Fig. 8 Distribution contour of dilation safety factor (SF) in the surrounding rock of Cavity A after operation of 1.6a

4 結(jié)論

(1)根據(jù)金壇鹽巖地層結(jié)構(gòu)特征、聲吶測腔數(shù)據(jù)和注采氣運行參數(shù)建立了A 鹽穴儲氣庫穩(wěn)定性評價的三維地質(zhì)力學(xué)模型，對該腔體穩(wěn)定性進行了分析，確定出該腔體體積收縮率過大的原因。

(2)計算結(jié)果表明：A 腔體整體穩(wěn)定性較好，腔體體積收縮率約為0.54%，不存在發(fā)生較大面積腔壁失穩(wěn)的風險；由于A 腔體頂部結(jié)構(gòu)跨度大且較為平直，存在頂部局部失穩(wěn)掉塊的風險，建議后期加強監(jiān)測，確保頂部失穩(wěn)在可控范圍內(nèi)。

(3)通過對比2017 年和2019 年聲吶測腔結(jié)果認為導(dǎo)致A 腔出現(xiàn)體積收縮率過大的主要原因是聲吶測量誤差而不是真實的體積收縮。A 腔體形狀較為復(fù)雜，腔壁上存在大量不規(guī)則空間，這些空間在聲吶測腔時不能夠被精確測量到，從而導(dǎo)致出現(xiàn)A 腔體體積收縮率過大。