多場作用下鹽穴儲氣庫腔體穩定性的數值模擬研究

宮丹妮 ，李景翠，萬繼方，莊曉謙，Beno?t Brouard，王琳琳

1 中國石油大學(北京)安全與海洋工程學院，北京 102249

2 中國石油集團工程技術研究院有限公司，北京 102206

3 Brouard工程公司，巴黎 75003

0 引言

金壇鹽巖儲氣庫是中國第一座正式投產使用的鹽穴儲氣庫，鹽穴儲氣庫工程技術經過十余年的不斷發展，目前已有20 余口井次進入投產運行階段[1-2]。然而鹽巖儲氣庫常年深埋地下，部分鹽腔形狀不規則，腔體內溫度壓力發生變化時其內壁會呈現應力集中，導致腔體局部位置破壞特征明顯[3]。儲氣庫長期力學穩定性研究中的難點在于鹽巖地下儲氣庫運營過程中會受到溫度場、應力場等多場共同作用，導致鹽腔變形過程變得復雜。因此本文采用一種新型三維有限元仿真模擬軟件LOCAS3D對金壇X生產井和現場實際注采壓力工況，并采用多場作用下的合理理論模型進行了數值模擬，將模擬結果與現場聲吶測腔結果對比，并預測儲氣庫運營期間腔體穩定性薄弱區域。

1 基本理論模型

1.1 鹽巖儲氣庫蠕變本構方程

實際上在長達數年蠕變的過程中，鹽巖瞬時彈性變形和初級蠕變變形與穩態蠕變相比較小。在工程中可以忽略用時較短的加速蠕變過程，因此鹽巖的蠕變行為可以通過穩態蠕變速率來衡量[5]。在國內外針對鹽巖蠕變性提出描述鹽巖蠕變的本構方程的基礎上[6-7]，金壇鹽巖蠕變行為考慮表征有效應力與溫度影響下的Norton-Hoff蠕變本構模型：

式中R=8.3142 J/(kg·K)，表示為氣體常數；T為絕對溫度，K；A為巖鹽材料常數，MPa-n/a；Q為鹽巖材料熱值常數，J/kg；n為應力指數，針對鹽巖一般取值為3～6。國內外學者通過對Avery Island鹽礦[8]取芯后的試樣進行的大量多階段蠕變測試得到的計算參數[9]：A=13000 MPa-n/a，n=3.14，Q/R=6495 K。

1.2 鹽巖受熱膨脹

鹽巖受熱膨脹的特征導致其在運行期間會產生一定量的膨脹，根據線性熱膨脹定律，考慮儲氣庫運營期間鹽腔總應變量可以表示為蠕變引起的應變與鹽巖熱膨脹后的應變之和，有：

式中εth表示為鹽巖熱應變張量的大小，無量綱；βT表示鹽巖的熱膨脹系數，其值取4×10-5/℃；εcr表示為因蠕變引起的應變大??；ε代表鹽腔總應變量的大小，無量綱。

1.3 儲氣庫失效準則

為了保證儲氣庫安全穩定的運營，運營期間要滿足無拉應力產生和無膨脹破壞區兩個原則：

(1)有效拉應力準則：由于鹽巖的抗拉強度Tsalt較低，其值為2 MPa左右甚至更小。通常規定有效拉應力為正值，有效壓應力為負值。Brouard等人認為當鹽腔內壓大于最小壓應力σmax時就會存在微裂縫的現象[10-12]：

(2)膨脹準則：通常為了更加直觀地反映出鹽體及其周圍產生的拉伸區，規定一個適用于多種標準的膨脹損傷安全系數FOS，當FOS ＜ 1 表明該區域有膨脹破壞的趨勢。由于膨脹損傷與加載路徑有關[13-15]，考慮洛德角的非線性損傷準則如下：

式中D1、D2、n為鹽巖材料參數，無量綱；標準應力σo=1 MPa；T0為待定鹽巖抗拉強度，MPa；ψ為洛德角，范圍在-30°～30°之間。

2 數值模擬

2.1 數值模擬方法

由于鹽巖儲氣庫長期埋于地下，運行期間會受到多種載荷的作用。本次模擬將采用LOCAS3D構建數值模型，通過有限元方法針對鹽巖這種流變性強的巖石進行數值計算。相比于FLAC3D與Abaqus這兩種常規模擬軟件，LOCAS在處理鹽腔有限元分析問題上具有一定的優勢，其能夠模擬出鹽腔的非線性和時變力學行為，進一步實現對鹽巖等一些巖土材料進行三維受力特性模擬和蠕變力學行為分析。

2.2 模型建立與基本參數

本次模擬基于金壇儲氣庫中代表井X井的現場實際數據，該井于2010年9月正式進入注氣投產階段，截至目前已經平穩運行十余年，其腔體形態為倒梨形?，F根據X井于2010年進入投產階段后的聲吶測腔結果建立真實三維幾何模型，由于鹽層滲透率極低，故不包含為本次模擬中主要的巖石物性參數。設定模型為一個長寬各為200 m，高為1200 m的長方體，包含4個巖層：0～425 m深的泥巖層I、425～590 m深的玄武巖層、590～988 m深的泥巖層Ⅱ和988～1200 m深的鹽巖層(圖1 左圖)，其數值模型基本力學參數與巖石蠕變計算參數的設置如表1、表2 所示。

表1 地層基本力學參數Table 1 Basic mechanical parameters of formation

表2 巖石蠕變計算參數Table 2 Calculation parameters of rock creep

圖1 金壇X井腔體幾何模型Fig. 1 Geometry model of Jintan X well cavity and loading pressure cycle

模型探明體積可達188 346.9 m3，表面積為16 637.81 m2。洞穴頂部埋深1015.67 m，總體高度為73.08 m。模型采用嵌套式雙重網格劃分，不同地層采用金字塔網格劃分，靠近鹽腔周圍網格元素采用六面體小網格劃分并嵌入幾何模型中，腔體壁面上平均單元大小為1.97 m。

2.3 初始條件

金壇鹽巖礦區的地質資料顯示，鹽巖礦區的巖層及其上下蓋層均為水平或近水平，地質構造應力不大。地面初始地應力取0 MPa，沿深度方向0～1200 m設置應力梯度，根據地層平均密度計算得到地層梯度，即埋深425 m處泥巖層I與玄武巖層界面處地應力為9.12 MPa，埋深590 m處玄武巖層與泥巖層Ⅱ界面處地應力為13.86 MPa，埋深988 m處泥巖層Ⅱ與鹽巖層界面處地應力為22.45 MPa，埋深1200 m處鹽巖層地應力為27.03 MPa。因金壇X溶腔中心處于地下1052 m深處，計算模型所處深度范圍在1015～1088.67 m，因此腔體頂部與腔體底部圍巖初始地應力設定為22.84 MPa與24.49 MPa。

模型假定地面溫度為定值，其值考慮為環境溫度T0≈10 ℃。地層溫度隨著深度的增加而增加，此模型域內溫度的分布規律可以表示為：

式中z為某一地層深度，m；TR(z)為該深度處地層的溫度。因此鹽腔中部埋深的地層溫度為41.56 ℃。假設同一深度地層溫度均勻不變，則初始條件如圖2 所示。

圖2 初始地層壓力與地熱梯度Fig. 2 Initial formation pressure and geothermal gradient

2.4 循環工況加載情況

模擬采用的15年現場注采工況如圖3 圖所示，一年內鹽腔在第91 天起經歷共12 天的最低壓力運營期(80.9 barg)，并在365 天后經歷30 天的最高壓力運營(149.2 barg)。由于腔體內壓達到最低后的幾個小時內發生膨脹破壞的可能性最大，因此選取6 個分析點，分別為運營期開端時刻(A)與結尾時刻(F)以及儲氣庫分別運營1年(B)、5年(C)、8年(D)和10年(E)最低壓力處對應時刻。

圖3 加載洞穴壓力循環和每年加載洞穴壓力循環Fig. 3 The cave pressure cycle and the cave pressure cycle in one year

計算中設定地面初始地應力取0 MPa，沿深度方向 0～1200 m 設置應力梯度，根據地層平均密度計算四個界面處地應力分別為9.12 MPa、13.86 MPa、22.45 MPa和27.03 MPa。金壇X溶腔中心處于地下1052 m深處，計算模型所處深度(1015～1088.67 m)較大，因此考慮地應力狀態為靜水應力狀態。

3 模擬結果分析

在周期運行壓力下的數值計算研究中，不同埋深處的腔壁變形量存在差異。從圖4 可以看出，越接近腔體中部埋深處的位移量越小，腔體上部變形量明顯增加并隨著埋深的減小而增加。同時隨著工作年限的不斷增加，腔壁位移變形量逐漸增加，而變形量增加的速率逐漸減小。為研究不同埋深處鹽腔腔壁的變形情況，考慮針對鹽腔上部、腔體中部、鹽腔下部3部分進行研究。對此X鹽腔采用以1036 m、1052 m、1061 m、1080 m這4 種不同界面深度為對象進行對比分析。

圖4 不同埋深處鹽腔腔壁變形情況Fig. 4 Deformation of salt cavity wall in different buried depths

3.1 鹽腔變形規律

鹽腔的不規則性導致腔壁處也呈現出不規則的變形，因此需要明確腔體產生大變形量的危險位置并分析該位置處的應力狀態。圖5 給出了儲氣庫運營10年后鹽腔腔壁位移情況。腔體中部位置處(埋深1052 m)的腔體直徑達到最大，在深度一定的情況下，選取了腔體邊緣處兩個坐標點(61，100)、(140，100)來分析腔壁位移情況。整體上看腔體變形量較小，儲氣庫具有較好的穩定性。腔體壁面總變形量較大區域出現兩處：一位于腔體中部上部埋深1036 m處，該區域在運行期間產生的最大位移平均為73.74 mm；二位于腔體中部下部埋深1080 m處，該區域產生的最大位移平均為55.52 mm。

圖5 橫剖面左、右腔壁上垂直向位移變化情況Fig. 5 Changes in vertical displacement on the left and right cavity walls of the transverse section

不同埋深處的橫剖面處的變形情況差別較大。如圖6 所示，大變形區域集中分布在鹽穴腔體上部(1036 m)和下部(1080 m)，腔體頂板的垂向位移方向向下、底板向上。同時相比于下部變形量，腔體頂板向下的垂直變形量更大，腔體中上部受到圍巖的壓力在上覆巖層的重力作用影響下變大。另外存在較大位移量的鹽穴腔體上部(1036 m)和下部(1080 m)屬于腔體內壁突出區域，而腔體中部平滑區域(1061 m)變形量最小。

3.2 鹽腔溫度演化規律

圖7 給出了儲氣庫投產后腔體內部流體溫度隨時間的變化曲線，整體上看腔體溫度呈現先增加后趨于穩定的狀態，造腔期至運營初期間溫度上升較快。水溶造腔采用了相對于鹽巖層較冷的地表水，當地表水注入鹽巖層后會出現溫度場擾動并重新分布的過程，此后腔內流體溫度逐年升高，地層與腔體流體之間的溫差隨時間推移而不斷變小。

圖7 腔內流體溫度隨時間變化曲線Fig. 7 The fluid temperature in the cavity changes with time

結合壓力加載情況(圖3)我們發現注氣加壓過程溫度出現上升趨勢，采氣降壓過程溫度下降趨勢。此時腔內流體壓力會隨著溫度場的擾動出現變化，即腔體內部的應力場重新分布。注氣過程中由于溫度的影響腔體受熱膨脹，圍巖間膨脹程度的差異導致腔體壁面產生壓應力；采氣過程中腔體冷卻收縮，受到圍巖的約束而產生拉應力。

3.3 鹽腔破壞規律

圖8 給出了儲氣庫運營第1年、第5年、第8年和第10年時壓力最小值處有效應力分布情況，在這4種情況下發現整個運營期間洞穴壁面上的受壓狀態有向受拉狀態發展的趨勢，腔體中上部1036 m處和中下部1080 m處應力向受拉狀態發展較快，腔體中部下部1059 m處有效壓應力最小且運營期間并無變化。

圖8 腔壁兩側(上)與整體(下)有效應力分布情況Fig. 8 Effective stress distribution on both sides of the cavity wall (left) and the whole (right)

對比圖8(下)上下兩組圖片發現，腔體中上部(1036 m)和中下部(1080 m)突出區域出現了明顯的應力集中，此時這兩個位置受到的壓應力最小，腔壁受力情況有向拉應力發展的趨勢。圖9 給出了儲氣庫分別運營第1年、第5年、第8年和第10年壓力最低處的安全系數橫縱剖面分布云圖，對比發現腔體突出區域安全系數小于1，其膨脹程度與平滑區域相比更加明顯，另外鹽腔整體膨脹程度出現逐年向外擴散的趨勢。在第一次循環壓力最低處出現大范圍膨脹區，并隨運營年限逐漸減小。說明運營初期1～2年間壓力最低處腔壁突出區域更易產生膨脹破壞，因此腔體局部突出區域是可能發生失效破壞的危險區域，該腔壁部位在重力作用下容易出現脫落現象。

圖9 第1、5、8、10年運營最低壓力處安全系數FOS分布云圖Fig. 9 The distribution cloud map of the safety factor FOS at the lowest operating pressure in the first 1, 5, 8 and 10 years

3.4 現場對比

對于鹽穴儲氣庫，溶腔可用體積以及腔壁位移量一直是現場人員以及研究學者重點關注的問題。溶腔可用體積是衡量儲氣庫儲氣性能以及使用壽命的重要指標，而隨多周期注采施工后的可用體積變化率則直接影響了儲氣庫長時間運營的穩定性?，F根據X井于2015年進行的第4 次聲吶測腔結果，對比2010年聲吶測腔結果(圖10 左)。結果表明：X井溶腔可用體積由2010年的188 346.9 m3減至185 502.8 m3，X井腔體運行第6年體積收縮率為1.51%。其中1080 m埋深處的發生較明顯的體積收縮，收縮量約78 mm。數值模擬金壇X井從2010年9月正式注氣投產后，在運營的10年期間儲氣庫體積收縮量隨時間的變化情況如圖9(右)所示。計算結果表明儲氣庫運營1年的體積收縮率約為0.23%，對比2015年腔體形態監測結果，發現該數值模擬結果與監測值之間誤差僅為0.021%，小于1%，說明該模擬方法具有較高的準確性和可靠性，為其他指標的分析提供了依據。

圖10 X腔體現場測腔(左)與模擬體積收縮情況(右)對比Fig. 10 Comparison of volume shrinkage between Sonar cavity test results (left) and Simulation (Right) of Cavity X

4 總結

本文針對儲氣庫在多場作用下長時間運營的穩定性問題，利用有限元方法實現了多場耦合基本模型的數值求解過程，并根據金壇鹽礦地層的結構特征、聲吶測腔數據以及現場加載壓力周期，采用嵌套式雙重網格劃分方式建立了X腔三維幾何模型。相比于傳統二維軸對稱模型，三維幾何模型更接近現場情況，數值計算后得出的結果更具有參考性。文中通過現場聲吶測腔結果驗證了前期模擬結果，得到了如下的結論：

(1)儲氣庫腔體內壁的位移變形量隨時間的增加而增加，而變形量增加的速率逐漸減小。由于上覆巖層重力作用該區域的位移變形更明顯，腔壁受到重力作用易發生脫落現象。建議儲氣庫建腔期要減少腔體突出面積的大小，運營期間應關注初期的腔壁突出區域；

(2)建腔初期和運營1～2年腔內流體溫度場重新分布，該階段溫度上升較快；采氣降壓過程中腔內溫度明顯下降，對應該過程中腔壁突出區域體積收縮情況明顯；

(3) 儲氣庫運營初期腔內壓最低時膨脹現象最明顯，整體上看鹽腔膨脹區有向外擴散的趨勢，而后隨著運營年限的增加而逐漸變??；第一次循環中腔壁突出區域出現明顯的應力集中，此時受到的壓應力最小，其受壓狀態有向受拉狀態轉變的趨勢。結合以上分析得出腔體突出區域是儲氣庫穩定性的薄弱區域，為儲氣庫造腔工程形狀優化提供理論依據。