多场作用下盐穴储气库腔体稳定性的数值模拟研究

宫丹妮 ，李景翠，万继方，庄晓谦，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月正式进入注气投产阶段，截至目前已经平稳运行十余年，其腔体形态为倒梨形。现根据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 m³，表面积为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 现场对比

对于盐穴储气库，溶腔可用体积以及腔壁位移量一直是现场人员以及研究学者重点关注的问题。溶腔可用体积是衡量储气库储气性能以及使用寿命的重要指标，而随多周期注采施工后的可用体积变化率则直接影响了储气库长时间运营的稳定性。现根据X井于2015年进行的第4 次声呐测腔结果，对比2010年声呐测腔结果(图10 左)。结果表明：X井溶腔可用体积由2010年的188 346.9 m³减至185 502.8 m³，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) 储气库运营初期腔内压最低时膨胀现象最明显，整体上看盐腔膨胀区有向外扩散的趋势，而后随着运营年限的增加而逐渐变小；第一次循环中腔壁突出区域出现明显的应力集中，此时受到的压应力最小，其受压状态有向受拉状态转变的趋势。结合以上分析得出腔体突出区域是储气库稳定性的薄弱区域，为储气库造腔工程形状优化提供理论依据。