含水层储气库注入阶段盖层力学完整性数值模拟分析

杨军伟， 仲家锐， 贾善坡， 滕振超， 杨典森

(1. 东北石油大学 土木建筑工程学院，黑龙江 大庆 163318； 2. 哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001； 3. 东北石油大学 非常规油气研究院，黑龙江 大庆 163318； 4. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071 )

0 引言

天然气以清洁、高效等特点成为现代化城市燃气发展的主要方向[1-3]。地下储气库的加速建设可弥补中国天然气储存能力不足，保障国家能源安全并完善天然气产业链。根据储层、地质等特点，地下储气库主要分为油气藏型、含水层型、盐穴型和废弃矿坑型。油气藏型储气库改建是最常用、最经济的一种储气方法，但在大型工业城市中心和大城市附近，并非都有适合建设地下储气库的枯竭油气田，且长距离管线输送对安全供气的要求很高。建造含水层型地下储气库所需地层结构分布较广泛，一般在大城市周围可以找到合适的建库区域[4-6]。大量天然气注入深层地质构造产生巨大的超压和应力状态变化，可能危及盖层的力学完整性，导致天然气泄漏。因此，有必要开展对含水层储气库注入阶段盖层力学完整性的研究，以确保注入的天然气能在预期寿命内安全储存在地下含水层储气库中。

人们研究含水层储气库多集中于CO2封存时盖层的稳定性。对于阿尔及利亚的Salah天然气项目，RUTQVIST J等[7]建立弹性三维模型， 探讨CO2注入后盖层隆起的原因和机理，但模型存在很多不确定因素，不适合实际应用，同时忽略天然气在注入过程中产生的应力变化对盖层剪切及拉伸破坏的影响；VILARRASA V等[8]通过实验测试具有代表性的原位条件下CO2储存的岩石类型，将测得的参数作为数值模拟的输入数据，模拟CO2注入对盖层剪切破坏的影响，但忽略注入过程对盖层变形的影响；在恒应力条件下，FAVERO V等[9]研究CO2注入对盖层流体力学性能的影响，但忽略注入过程引起的应力变化对盖层地质力学完整性的影响；HADIAN P等[10]应用4种实验技术，评估超临界CO2注入对盖层岩石物性的影响，但忽略注入过程对盖层力学性能的影响；在弹性状态下，VILARRASA V等[11]研究CO2注入位置和断层渗透率对断层稳定性等地质力学的影响，以及注入井周围盖层稳定性的影响，但未考虑盖层参数对盖层力学完整性的影响；KIM K等[12]应用数值模拟，研究CO2注入对盖层岩石地质力学性质的影响，但未考虑注入速率对盖层地质力学性质的影响；LI Chao等[13]选择热弹性力学本构模型，采用热—流—固耦合模拟方法，研究深层含水层注入大量CO2后盖层的变形和地质力学的不稳定，但未考虑地层参数对盖层变形和地质力学不稳定影响的敏感性；贺凯[14]应用实验和数值模拟相结合的方法，在CO2地质封存过程中，对储盖层体系的地质力学作用影响进行研究，但未考虑对盖层力学作用影响的敏感性。对于天然气注入过程中盖层的地质力学变化的研究相对较少，而应用数值模拟技术模拟含水层储气库天然气注入对盖层力学完整性的研究更少。

笔者利用ABAQUS软件，建立含水层储气库天然气注入阶段的数值模型，探讨含水层储气库在天然气注入过程中盖层的力学行为特征，分析盖层参数、天然气注入速率和地应力因数对盖层力学完整性的影响，为我国含水层地下储气库盖层稳定性的研究提供技术支持和借鉴，也为工程实际提供参考。

1 地质力学模型

1.1 计算条件

以某含水层目标为研究对象，地层为相对完整型背斜构造，储层为砂岩，盖层为泥岩，前期筛选结果表明在该区块建设含水层储气库成功率较高。取99块储层岩心样品进行物性测试，孔隙度最小为2.5%，最大为15.4%，平均为10.0%；渗透系数最小为4.00×10-10m/s，最大为1.29×10-6m/s，平均为1.00×10-7m/s。取16块盖层岩心样品进行物性测试，孔隙度最小为6.4%，最大为26.3%，平均为10.0%；渗透系数最小为1.22×10-13m/s，最大为3.75×10-9m/s，平均为1.00×10-12m/s；突破压力大于10.00 MPa，满足静态密封性要求。

对储层和盖层岩心样品进行三轴压缩实验，不同围压下储层和盖层岩石的应力—应变曲线分别见图1和图2。由图1-2可知，储层岩石的弹性模量为4.69～24.29 GPa，平均为11.31 GPa，泊松比为0.21～0.39，平均为0.29；盖层岩石的弹性模量为5.59～9.13 GPa，平均为7.18 GPa，泊松比为0.10～0.14，平均为0.11。对岩石应力—应变数据进行处理，得到储层和盖层岩石对应的强度参数：内摩擦角分别为28.2°和25.0°；黏聚力分别为17.23和20.53 MPa。

图1 不同围压下储层岩石应力—应变曲线(80 ℃)Fig.1 Stress-strain curves of reservoir rocks under different confining pressures(80 ℃)

图2 不同围压下盖层岩石应力—应变曲线(常温)Fig.2 Stress-strain curves of caprock rocks under different confining pressures(room temperature)

对储层和盖层部分岩心样品进行间接拉伸实验。储层和盖层试样的间接拉伸破坏曲线见图3。由图3可知，储层岩石抗拉强度为7.44～13.13 MPa，平均为10.31 MPa；盖层岩石抗拉强度为4.15～7.44 MPa，平均为5.52 MPa。

由于地层岩心测试资料较少，参照储层、盖层岩石力学实验数据测试结果及相关文献资料[15-16]确定模型计算参数(见表1)。

图3 储层和盖层试样间接拉伸破坏曲线Fig.3 Indirect tensile failure curves of reservoir and caprock samples

表1 模型计算参数

1.2 模型建立

研究区构造相对平缓，属于微幅构造且储层厚度较大，为水平构造，数值模拟采用二维轴对称模型[17]。根据地质资料建立二维模型(见图4)，模型为矩形，长度为5 000 m，由盖层、储层和底层组成，盖层、储层、底层厚度分别为63、276、193 m(见图4(a))。盖层上部有深度为3 248.0 m的地质层，以等效覆盖层表示。天然气沿储层中一点以0.035 m3/s的恒定注入速率向储层注入6个月，注入点深度为3 487.0 m。

由于研究目标是盖层地质力学变化而不是孔隙压力本身，模拟采用单相流体流动，流体相为非必要限制条件[18]。

以ABAQUS软件中的流—固耦合计算模块为求解器，求解流—固耦合问题，通过黏土Biot一维固结过程的数值模拟与理论解对比，对数值模型进行验证[19]。模型网格划分采用CAX4RP单元(四结点轴对称四边形单元)、双线性位移、双线性孔隙压力。在储层和盖层中，越靠近注入点网格划分越细，越远离注入点网格划分越粗，二维模型有限元网格划分见图4(b)。

假设3种相互正交的主应力中有一种具有垂直方向，垂向应力的大小由上覆岩体和流体的质量决定[20]。不同深度的垂向应力可表示为

σv=ρgz,

(1)

式中：σv为垂向应力；ρ为覆盖层平均密度；g为重力加速度；z为深度。考虑覆盖层岩石的自重，由式(1)可得盖层顶部的垂向应力为83.1 MPa(ρ=2 608 kg/m3，g=9.81 m/s2，z=3 248 m)。初始水平应力为各向异性，模型的地应力因数为0.76。初始孔隙压力梯度为10.50 MPa/km，在天然气注入开始时，注入点的初始孔隙压力为36.61 MPa。

数值模拟采用总孔隙压力的方法，力学边界：上边界为上覆压力；下边界为法向位移约束；左边界为轴对称边界；右边界为法向约束。渗透边界：上边界为渗透边界；下边界不排水；左边界不排水；右边界为封闭边界。

图4 二维模型的物理模型及网格划分Fig.4 Physical model and mesh generation of two dimensional model

1.3 盖层破坏准则

岩石破坏准则为给定岩石破坏时应力状态需要满足的条件，在条件外岩石不能承受施加的载荷[21]。为确定盖层在天然气注入过程中是否稳定，需要定义盖层的破坏准则。如果应力状态低于破坏准则，则盖层是稳定的。采用莫尔—库伦破坏准则[22]模拟岩石破坏的表达式为

τ=c+(σn-pf)tanφ,

(2)

式中：τ为剪切面上的剪应力；c为黏聚力；σn为剪切面上的正应力；pf为孔隙压力；φ为内摩擦角。岩石的黏聚力和内摩擦角见表1。岩石剪切面上的剪应力和正应力分别为

(3)

(4)

式(3-4)中：σ1为最大主应力；σ3为最小主应力；θ为岩石破坏面与σ3方向的夹角[23-24]。当τ≥c+(σn-pf)tanφ时，盖层发生剪切破坏。

岩石拉伸破坏发生在拉伸应力超过岩石抗拉强度时。材料的抗拉强度是材料破裂之前能承受的最大拉应力。根据材料第一强度理论，拉伸破坏准则[25-26]为

(5)

2 数值模拟结果

为使初始应力场和模型边界条件之间达到力学平衡，应进行地应力平衡，其计算结果的准确性决定模拟结果的可靠性。采用ODB文件导入法对模型进行地应力平衡，得到初始地应力条件下垂向位移数量级为10-7m；垂向应力为压应力，自上而下逐渐增大；孔隙压力自上而下随深度线性增加，初始地应力平衡效果较好。

2.1 孔隙压力

离注入点水平距离10 m处取一条垂向测线CE，天然气注入后，不同注入时间的测线CE孔隙压力变化曲线见图5。由图5可知，盖层的孔隙压力随深度的增加而增加，最大孔隙压力位于注入6个月时的盖层底部(深度为3 311.0 m)，为40.92 MPa，比初始孔隙压力增加6.15 MPa。模型最大孔隙压力位于天然气注入6个月时的注入点(深度为3 487.0 m)，为59.86 MPa，是初始孔隙压力的1.635倍。

当天然气注入6个月时，离注入点不同水平距离取4条垂向测线，得到测线AB、CD、FG、HI的孔隙压力变化曲线(见图6)。由图6可知，离注入点越近的垂向测线孔隙压力变化越大，离注入点较远的垂向测线孔隙压力变化较小。

图5 不同注入时间的测线CE孔隙压力变化曲线Fig.5 Curves of pore pressure variation of the CE under different injection time

图6 不同水平距离的垂向测线孔隙压力变化曲线Fig.6 Curves of pore pressure variation of the vertical measuring lines under different horizontal distances

不同注入时间的储、盖层界面孔隙压力增量变化曲线见图7。由图7可知，最大孔隙压力增量发生在天然气注入6个月时的注入点附近，为6.118 MPa，随离注入点水平距离增加，孔隙压力增量变小。

盖层底部测点J孔隙压力随注入时间变化曲线见图8。由图8可知，随注入时间增加，测点J孔隙压力增大，且孔隙压力在天然气注入1个月内迅速上升，变化速率较大；随注入时间继续增加，孔隙压力的增加速率逐渐变小。

图7 不同注入时间的储、盖层界面孔隙压力增量变化曲线Fig.7 Curves of pore pressure variation at the reservoir-caprock interface under different injection time

图8 盖层底部测点J孔隙压力随注入时间变化曲线Fig.8 Curve of pore pressure variation of the measuring point J at the bottom of the caprock with the injection time

2.2 盖层位移

孔隙压力增大引起的应力变化导致储层孔隙膨胀，储、盖层界面上移。天然气注入1个月和6个月时，储、盖层界面位移变化曲线见图9。由图9(a)可知，天然气注入6个月时，注入点附近出现最大垂向位移，为0.073 m，且随离注入点距离增加，垂向位移逐渐变小；由图9(b)可知，天然气注入6个月时，在离注入点水平距离1.0 km处产生最大水平位移，为0.032 m，且水平位移在离注入点一定水平距离内随距离增加而变大，超过一定水平距离后随距离增加而变小。

图9 不同注入时间的储、盖层界面位移变化曲线Fig.9 Curves of displacement at the reservoir-caprock interface under the different injection time

储、盖层界面最大位移随注入时间变化曲线见图10。由图10可知，随注入时间增加，储、盖层界面最大垂向和最大水平位移逐渐增大，且位移增加的速率在注入1个月内变化显著，之后位移增加的速率逐渐减小。

图10 储、盖层界面最大位移随注入时间变化曲线Fig.10 Curves of maximum displacement at the reservoir-caprock interface with injection time

2.3 测线CE应力

不同注入时间的测线CE应力变化曲线见图11。由图11可知，盖层在同一注入时间内垂向应力和水平应力随深度增加先减小后增大；相同深度时，盖层的垂向应力和水平应力随注入时间的增加而增大，且垂向应力和水平应力未出现正值，即该工况下盖层未发生拉伸破坏。

2.4 盖层剪切破坏

采用莫尔—库仑失效准则[22]，得到测点J在天然气注入过程中应力状态变化(见图12)。由图12可知，天然气注入过程中，莫尔圆没有达到盖层强度包络线。盖层在天然气注入6个月时没有发生破坏，并保持稳定。盖层强度包络线的截距较大且斜率较陡，预计莫尔圆不会触及盖层破坏线。如果采用更高的注入速率，则最小主应力可能变为负值。因此，需要进行参数敏感性分析，确定天然气注入速率对盖层剪切破坏的影响。

图11 不同注入时间的测线CE应力变化曲线Fig.11 Curves of stress at measuring line CE under different injection time

图12 测点J应力莫尔圆状态变化Fig.12 State change disgram of stress Mohr circle at a measuring point J

3 参数敏感性

对盖层参数、天然气注入速率和地应力因数[28-29]进行敏感性分析，评价参数对盖层变形及应力变化的敏感性。参数敏感性分析输入数据见表2，其中A0、B0、C0、D0和E0为原工况数据。取盖层底部为水平测线，研究参数变化对测线最大位移的影响，分析各参数对盖层变形的敏感性；在盖层底部取测点J，研究参数变化对测点J应力的影响，讨论各参数对盖层拉伸和剪切破坏的敏感性。

表2 盖层参数、天然气注入速率及地应力因数敏感性分析输入数据

3.1 盖层弹性模量

天然气注入6个月时，盖层底部最大位移随盖层弹性模量变化曲线见图13。盖层底部测点J的正应力和剪应力随弹性模量变化曲线见图14。由图13-14可知，盖层弹性模量从3.59 GPa增加到10.77 GPa，盖层底部最大水平位移减小3.760×10-3m，最大垂向位移减小2.280×10-3m；测点J的正应力增加0.41 MPa，剪应力降低0.27 MPa。因此，盖层弹性模量对盖层的最大水平位移和最大垂向位移影响较小，且对最大水平位移的影响比最大垂向位移的大；对盖层拉伸和剪切破坏趋势影响较小。

图13 盖层底部最大位移随盖层弹性模量变化曲线Fig.13 Curves of the maximum displacement at the bottom of the caprock with Youngs modulus of caprock

图14 测点J应力随盖层弹性模量变化曲线Fig.14 Curves of stress at measuring point J varying with Youngs modulus of caprock

3.2 盖层泊松比

天然气注入6个月时，盖层底部最大位移随盖层泊松比变化曲线见图15。盖层底部测点J正应力和剪应力随泊松比变化曲线见图16。由图15-16可知，盖层泊松比从0.10增加到0.35，盖层底部最大水平位移增大0.640×10-3m，最大垂向位移增大0.660×10-3m；测点J正应力降低0.08 MPa，剪应力增加0.04 MPa。因此，盖层泊松比对盖层的最大水平位移和最大垂向位移影响较小，可以忽略不计；对盖层拉伸破坏和剪切破坏趋势几乎没有影响。

3.3 盖层渗透系数

天然气注入6个月时，盖层底部最大位移随盖层渗透系数变化曲线见图17。盖层底部测点J正应力和剪应力随渗透系数变化曲线见图18。由图17-18可知，盖层渗透系数从10-14m/s增加到10-8m/s，盖层底部的最大水平位移增大7.080×10-3m，最大垂向位移减小6.940×10-3m；测点J的正应力增加0.01 MPa，剪应力增加0.59 MPa。因此，盖层渗透系数对盖层的最大水平位移和最大垂向位移影响较大；对盖层拉伸破坏趋势几乎没有影响，对剪切破坏趋势影响较小。

图15 盖层底部最大位移随盖层泊松比变化曲线Fig.15 Curves of the maximum displacement at the bottom of the caprock with Poisson ratio of caprock

图16 测点J应力随盖层泊松比变化曲线Fig.16 Curves of stress at measuring point J varying with Poisson ratio of caprock

图17 盖层底部最大位移随盖层渗透系数变化曲线Fig.17 Curves of the maximum displacement at the bottom of the caprock with permeability coefficient of caprock

3.4 天然气注入速率

天然气注入6个月时，盖层底部最大位移随注入速率变化曲线见图19。盖层底部测点J正应力和剪应力随注入速率变化曲线见图20。由图19-20可知，天然气注入速率从0.025 m3/s增加到0.055 m3/s，盖层底部最大水平位移增大27.010×10-3m，最大垂向位移增大62.490×10-3m；测点J的正应力降低4.30 MPa，剪应力降低0.66 MPa。因此，天然气注入速率对盖层最大位移、拉伸破坏和剪切破坏趋势影响较大，且对拉伸破坏的影响比剪切破坏的大。

图18 测点J应力随盖层渗透系数变化曲线Fig.18 Curves of stress at measuring point J varying with permeability coefficient of caprock

图19 盖层底部最大位移随天然气注入速率变化曲线Fig.19 Curves of the maximum displacement at the bottom of the caprock with injection rate of gas

图20 测点J应力随天然气注入速率变化曲线Fig.20 Curves of stress at measuring point J varying with injection rate of gas

3.5 地应力因数

天然气注入6个月时，盖层底部最大位移随地应力因数变化曲线见图21。盖层底部测点J正应力和剪应力随地应力因数变化曲线见图22。由图21-22可知，地应力因数由0.66增加到1.26，盖层底部的最大水平位移和最大垂向位移几乎没有改变；测点J正应力增加15.02 MPa，剪应力先降低7.51 MPa，后增加7.01 MPa。因此，地应力因数对盖层最大位移几乎没有影响，对盖层拉伸破坏和剪切破坏趋势影响很大。

图21 盖层底部最大位移随地应力因数变化曲线Fig.21 Curves of the maximum displacement at the bottom of the caprock with coefficient of stress

图22 测点J应力随地应力因数变化曲线Fig.22 Curves of stress at measuring point J varying with coeffcient of stress

盖层渗透系数和天然气注入速率对盖层变形的影响显著，天然气注入速率和地应力因数对盖层应力的影响显著。因此，在实际工程实践中，规划深层地质构造天然气注入方案时，应首先考虑盖层渗透系数、天然气注入速率和地应力因数等因素。

4 结论

(1)以拟建含水层储气库为依托工程，建立含水层储气库地质力学模型，采用ABAQUS软件对天然气注入阶段的盖层力学完整性、盖层变形及应力敏感性进行模拟分析。

(2)注入早期的注入点附近引起的孔隙压力变化和地质力学变化最大，是盖层中岩石破坏风险最大的关键部位。

(3)对盖层变形最为敏感的参数为天然气注入速率，其次为盖层的渗透系数；对盖层应力变化最为敏感的参数为地应力因数，其次为天然气注入速率；盖层泊松比和模型右侧边界类型对盖层变形和应力变化的影响不显著。