热力耦合作用下地下盐岩储气库注采运行稳定性研究

彭井宏，周军，胡承强，黄薪宇，梁光川

（1.西南石油大学石油与天然气工程学院，四川 成都 610500；2.中国石油西南油气田分公司安全环保与技术监督研究院，四川 成都 610500）

地下盐岩储气库具有安全性好、运行费用低和注采效率高等特点［1-3］，在天然气调峰方面发挥着重要作用。在地下盐岩储气库连续循环注采运行过程中，腔体内的气体不断被压缩和膨胀，导致腔体内压力和温度不断发生变化［4］，直接引起腔体围岩温度、应力分布的变化。温度变化时，由于腔体围岩内部不同的矿物颗粒热膨胀系数不同，因而发生内部变形约束，产生附加的热应力可能会导致腔体围岩的损伤破坏［5-6］。 因此，开展热力耦合下的盐岩储气库运行稳定性研究具有十分重要的意义。

目前，针对地下储气库稳定性的研究已经取得了一定进展。大部分研究都是先进行室内实验或者理论分析，而后将所得的实验结果或理论结果进行数值模拟，通过数值模拟结果分析地下盐岩储气库的稳定性。实验研究方面：李银平等［7］开展了盐穴地下储气库水溶造腔管柱的动力稳定性实验；邓朝福等［8］对平顶山和淮安2 处矿井的纯盐岩进行了单抽加卸载实验，以研究纯盐岩变形破坏过程的基本特征；井岗等［9］为确定盐穴储气库上限压力，在金坛储气库进行了小型水力压裂地应力测试。数值模拟方面：刘建平等［10］利用数值模拟手段研究了水平盐穴储气库在不同埋深下腔体长期运行的体积收缩情况；李景翠等［11］以金坛区块典型井为例，运用Locas 软件从腔体稳定性和体积收缩率两方面进行了数值模拟；王粟等［12］利用FLAC 3D 软件模拟了注采运行过程中盐岩性质、矿柱稳定性、注采气对腔体稳定性的影响。

在进行地下盐岩储气库稳定性研究时，国内大部分学者都注重压力对腔体的影响，往往忽略了温度也是影响地下盐岩储气库稳定性的重要因素［13］。 尽管已经有部分学者发现温度场对于地下盐岩储气库稳定性的影响不可忽略，但开展的研究工作主要聚焦于腔体温度分布和热力传导等理论研究［14-17］，因此，有必要针对注采过程中产生的热应力对盐岩储气库稳定性影响进行进一步研究。

本文首先分析了盐岩储气库运行过程中的稳定性判别指标，随后根据变质量基本方程建立了注采过程中腔内温度和压力的计算模型，在此基础上，对腔体进行了不同注采时间和注采速率下的数值模拟，判断腔体围岩损伤情况，从而判断腔体的稳定性。

1 稳定性判别标准

腔体围岩的完整性与整个地下盐岩储气库稳定性的关系非常紧密，可以说，研究地下盐岩储气库稳定性主要就是研究腔体围岩的损伤情况［18］。 本文在地下盐岩储气库运行稳定性的研究中，采用拉张损伤判据、剪切损伤判据和膨胀损伤判据3 个判据作为判定稳定性的标准［19］。

1）拉张损伤判据。 拉张损伤是指盐岩在受到拉张应力作用时，在非常低的强度下会产生破裂。根据第一强度理论，围岩不发生拉张损伤的条件为

式中：σ1为腔体围岩所受最大主应力，MPa；σt为腔体围岩抗拉强度，MPa（为避免围岩脆性破坏，围岩不允许出现拉张应力，故σt取值为0）；σx，σy分别为x，y方向上的正应力分量，MPa；τxy为切应力分量，MPa。

2）剪切损伤判据。 剪切损伤是指当盐岩所受的剪应力超过盐岩的剪切强度，岩体产生破裂。在实际工程中常使用Mohr-Coulomb 准则描述：

式中：φ 为内摩擦角，（°）；σ3为腔体围岩所受最小主应力，MPa；C为黏聚力，MPa。

3）膨胀损伤判据。 膨胀损伤是指在单轴和三轴压缩下，岩体的应变随着应力的提高，先增加到某一程度，然后又突然减小，使得盐岩产生微裂纹的现象。该现象从宏观上会导致盐岩体积增加。 线性膨胀损伤判据η＞1 时，盐岩发生膨胀损伤。 η 的表达式为

式中：J2为偏应力张量的第二不变量，MPa；a，b为膨胀系数（由岩体性质决定，本文取a=0.274，b=-10.8）；I1为应力张量的第一不变量，MPa。

2 温度和压力计算模型

为了模拟地下盐岩储气库注采过程中的围岩损伤情况，需要得到注采运行过程中腔体内的温度、压力变化规律作为数值模拟的边界条件。地下盐岩储气库大多采用“单腔单井”的建设模式，将地下腔体看作一个系统，所有质量和能量的进出均是发生在井筒与腔体的连接处（见图1）。

图1 地下盐岩储气库示意Fig.1 Schematic diagram of underground salt rock gas storage

地下盐岩储气库在运行过程中涉及的参数众多，若全面考虑建立的模型会非常复杂，故对地下盐岩储气库的部分参数和运行状态进行简化。假设：1）地下盐岩储气库不是新建投产的，不考虑水溶造腔过程中围岩形成的“冷带”；2）腔内流体全为天然气，忽略腔体底部的盐水［20］；3）腔内天然气温度、压力分布均匀；4）气体流入和流出盐腔的过程视为稳定流动；5）腔内天然气为理想气体，满足理想气体状态方程［21］。

地下盐岩储气库注采过程中，腔体部分的变质量系统基本方程由气体状态方程（式（4））、系统热力状态方程（式（5）、式（6））、系统质量守恒方程（式（7））和系统能量守恒方程（式（8））构成［22］：

式中：p为气体压力，Pa；V为气体体积，m3；m为气体质量，kg；T为气体温度，K；U为系统总内能，J；t为时间，s；cV，cp分别为气体质量定容热容、质量定压热容，J/（kg·K）；H为系统总焓，J；Δmi，Δme分别为流入、 流出系统的气体质量变化量，kg；ΔQ为系统与外界交换的热量，J；h为比焓，J/kg；v为流入、流出系统的控制体速度，m/s；Δm为控制体质量变化量，kg。

对式（1）—式（5）进行联立推导，可以得到注采过程中地下盐岩储气库腔体内压力和温度的微分方程，再对微分方程进行积分推导， 便可获得压力和温度的计算方程［23］。

采气过程中，腔内温度、压力方程：

式中：m01为采气开始时腔内气体质量，kg；ρ 为气体密度，kg/m3；q1为 采 气 速 率，m3/s；Y为 气 体 对 流 换 热 系数，W/（m2·K）；Y0为t=0 时的气体对流换热系数，取值为0.01 W/（m2·K）；A为腔壁面积，m2；R为气体常数，J/（kg·K）；p01为采气开始时腔内的压力，MPa；T01为采气开始时腔内的温度，K；d为井筒直径，mm。

注气过程中，腔内温度、压力方程：

式中：m02为注气开始时腔内气体质量，kg；q2为注气速率，m3/s；Tin为注入气体的温度，K；T02为注气开始时腔内的温度，K；p02为注气开始时腔内的压力，MPa。

3 数值模拟模型

本文数值模拟计算采用FLAC 3D 有限差分软件。地下盐岩储气库埋深为1 000～1 200 m，整个计算区域大小为500 m×500 m×700 m，其中盐岩层厚度为200 m，泥岩层总厚度为500 m，平均分布在盐岩层上下两端，盐岩层中含有2 层厚度为2 m 的泥岩夹层，模型采用鸡蛋形，即溶腔上半部分为椭球体，下半部分为球体（见图2）。

图2 地下盐岩储气库物理模型Fig.2 The physical model of underground salt rock gas storage

由于建立的地下盐岩储气库为轴对称模型，因此选用整个模型的一半进行后续的模拟计算，节省数值模拟的时间。 数值模型使用FLAC 3D 软件的命令流建立，共计生成24 056 个节点和21 888 个单元，单元类型均为六面体（见图3）［4］，选用岩体力学领域广泛应用的Mohr-Coulomb 准则作为模型计算准则［24］。数值模型静力学参数（见表1）均取自我国第一座已投入运行的地下盐岩储气库——金坛储气库。

表1 地下盐岩储气库静力学参数Table 1 Statics parameters of underground salt rock gas storage

图3 地下盐岩储气库数值模型及网格划分Fig.3 Numerical model and grid division of underground salt rock gas storage

在对地下盐岩储气库进行数值模拟计算前，需要先对所建模型的初始地应力场进行求解，通常称为静力平衡求解。地下盐岩储气库属于深埋工程，初始地应力场通常为构造应力场和自重应力场的叠加。 将盐岩的侧压系数设为1，泥岩的侧压系数设为0.85，结合已经设定的边界条件，采用快速应力边界法进行初始地应力场的求解。参考文献［25］中的盐岩地层温度实验，将地下盐岩储气库的地温梯度设为0.025 5 ℃/m，地表温度设为25 ℃。地下盐岩储气库的热力参数如表2 所示。在得到满意的地应力场后，即计算模型达到了静力平衡，开始进行初始温度场的设置。在计算模型顶端设置固定温度为42.85 ℃（316.00 K），底端设置固定温度为60.70 ℃（333.85 K），选用热传导的温度模型，求解得到初始温度场。

表2 地下盐岩储气库热力参数Table 2 Thermal parameters of underground salt rock gas storage

在进行数值模拟分析前，需要使用之前推导的腔内温度、压力模型作为数值模拟计算的边界条件［26-27］。地下盐岩储气库的运行参数和天然气的物性参数如表3 所示［28-29］。 根据文献［30］分析金坛储气库的注采数据， 可以得到单个盐腔在日常调峰时的平均注采量（40×104m3/d）和应急平均注采量（80×104m3/d），因此，在研究地下盐岩储气库运行过程稳定性时， 选择注采速率为40×104，60×104，80×104，100×104，120×104m3/d 的工况进行研究。 根据注采过程中注采速率、注采时间的不同和是否考虑温度影响的情况（不考虑温度影响对应压力作用， 考虑温度影响对应热力耦合作用），本文共设计了120 种工况。

表3 地下盐岩储气库运行参数和气体物性参数Table 3 Operating parameters and gas physical property of underground salt rock gas storage

4 注气过程模拟

选择注气速率分别为40×104，60×104，80×104，100×104，120×104m3/d 的工况，将计算得到的不同注气速率、 注气时间下的腔内温度和压力作为边界条件，在FLAC 3D 软件中进行数值模拟计算， 分析腔体围岩是否发生拉张、剪切、膨胀损伤，研究注气速率和注气时间对地下盐岩储气库注气运行稳定性的影响。

4.1 温度和压力

注气前，腔内气体已与围岩进行了充分的换热，此时，腔内气体与围岩表面的温度均取储气库中部的地层温度323.65 K。注气前，储气库会经历一个低压储气的阶段， 设定该阶段的压力为9 MPa。 运用式（11）、（12），可以计算出不同注气时间、注气速率时的腔内温度和压力。

不同注气条件下的腔内压力、温度变化趋势如图4 所示。 从图4a 可以看出：同一注气速率下，随着注气时间的增加，腔内压力呈线性增加，当注气速率为40×104m3/d 时，腔内压力增加到12.6 MPa，平均每天增加0.303 MPa；同一注气时间下，腔内压力随着注气速率的增大而增大；注气时间越长，注气速率对腔体压力的影响越明显。从图4b 可以看出：相同注气速率下，腔内温度随注气时间的增加呈线性增加，当注气速率为60×104m3/d 时，腔内温度增加到327.96 K，平均每天增加0.36 K；相同注气时间下，腔内温度随注气速率的增大而增大；注气时间越长，注气速率对腔体温度增加的影响越明显。

图4 不同注气条件下腔内压力和温度的变化Fig.4 Changes of pressure and temperature in cavity under different gas injection conditions

4.2 拉张损伤

拉张损伤通过围岩所受最大主应力进行判断。 对于地下盐腔而言，所受最大主应力小于0 则说明没有发生拉张损伤［18］。 注气速率40×104m3/d、注气时间12 d时，腔体围岩分别在压力作用、热力耦合作用下的拉张损伤情况如图5 所示。

图5 注气过程中的腔体围岩拉张损伤云图Fig.5 Tension damage nephogram of surrounding rock in cavity during gas injection

不同注气条件下的地下盐岩储气库腔体围岩所受最大主应力如图6 所示。 在FLAC 3D 软件中，最大主应力为正值时表示拉应力，为负值时表示压应力。

图6 不同注气条件下腔体围岩所受最大主应力Fig.6 Maximum principal stress of surrounding rock in cavity under different gas injection conditions

从图6 可知：注气过程中，腔体围岩均受压应力作用，无拉应力存在，因此可以判定，在设定的注气工况下，腔体围岩不发生拉张损伤；相同注气速率下，最大主应力随着注气时间的增加而减小，说明围岩所受压应力越来越大；相同注气时间下，注气速率越大，最大主应力越小，说明压应力越来越大；相同注气速率和时间下，热力耦合作用下的腔体围岩最大主应力始终小于压力作用下的数值，该差值随着注气时间的增加而增加。

4.3 剪切损伤

在FLAC 3D 软件中，可以通过查看模型的塑性区状态判断出现剪切损伤的区域。 塑性区是指该区域岩体承载的压力超过了其极限承载力，导致岩体出现不可恢复的区域。压力作用、热力耦合作用下的腔体围岩塑性区见图7。 塑性状态符号“none”“shear”“tension”分别表示正常无损伤、出现剪切损伤、出现拉张损伤，时间符号“n”“p”分别表示当前循环、以前循环。

图7 注气过程中的腔体围岩塑性区Fig.7 Plastic zone of surrounding rock in cavity during gas injection

4.4 膨胀损伤

注气速率120×104m3/d、注气时间12 d 时，腔体围岩分别在压力作用和热力耦合作用下的膨胀损伤情况见图8。不同注气条件下的腔体围岩膨胀损伤判据值见图9（空心图式为压力作用，实心图式为热力耦合作用，下同）。 由图可知，注气过程中，腔体围岩并未发生膨胀损伤。压力作用时，随着注气时间的增大，膨胀损伤判据值呈线性下降；相同注气时间下，随着注气速率的增大，膨胀损伤判据值也在下降。热力耦合作用时，膨胀损伤判据值总体也呈下降趋势，但下降范围很小。由此可见，热力耦合作用下的腔体围岩膨胀损伤判据值总体要高于压力作用下。

图8 注气过程中的腔体围岩膨胀损伤云图Fig.8 Expansion damage nephogram of surrounding rock in cavity during gas injection

图9 注气过程中的腔体围岩膨胀损伤判据Fig.9 Expansion damage criterion value of surrounding rock in cavity during gas injection

综上所述，在注气速率40×104～120×104m3/d、注气时间0～12 d 条件下，腔体围岩不发生拉张损伤。 注气导致腔体内的温度和压力增加，腔体围岩受热膨胀被更深处围岩约束，受到向内的压应力，随着注气速率和注气时间的增加，出现膨胀损伤的可能性反而变小。但是在注气时，无论是否添加温度场，腔体围岩均出现了剪切损伤，并且热力耦合作用下的腔体围岩塑性区远远大于压力作用下的腔体围岩塑性区。因此，在对地下盐岩储气库进行运行稳定性分析时， 需要充分考虑温度的影响，并且着重关注岩体的剪切损伤情况。

5 采气过程模拟

对地下盐岩储气库采气过程中的运行稳定性进行研究时， 选择采气速率分别为40×104，60×104，80×104，100×104，120×104m3/d 的工况，将计算所得的腔体温度和压力代入FLAC 3D 软件中作为数值模拟计算的边界条件，得到不同采气时间和采气速率下的腔体围岩损伤情况，并对这些损伤情况进行分析，判断储气库的稳定性。

5.1 温度和压力

采气前，腔内气体已与围岩进行了充分换热，此时腔内气体与围岩表面的温度均取储气库中部的地温温度323.65 K。 在高压储气阶段，腔体的压力设定为16 MPa，该压力为采气过程中腔内初始压力。 不同采气时间、 采气速率下的腔内压力和温度变化情况见图10。从图可知：相同采气速率下，随着采气时间的增加，腔内压力和温度均呈线性减小；相同采气时间下，采气速率越大，腔内压力和温度下降越大。

5.2 拉张损伤

采气速率120×104m3/d、采气时间12 d 时，腔体围岩分别在压力作用、热力耦合作用下的拉张损伤云图如图11 所示。不同采气条件下腔体围岩所受最大主应力如图12 所示。

图11 采气过程中腔体围岩拉张损伤云图Fig.11 Tension damage nephogram of surrounding rock in cavity during gas production

图12 不同采气条件下腔体围岩所受最大主应力Fig.12 Maximum principal stress of surrounding rock in cavity under different gas production conditions

由图12 可知：相同采气速率下，腔体围岩所受压应力随着采气时间的增加而减小； 采气速率120×104m3/d 时，腔体围岩在只有压力作用的情况下并未出现拉张损伤，而在热力耦合作用下，腔体围岩所受应力在12 d 时从压应力转变成拉应力， 大小为0.228 9 MPa，此时围岩已经发生了拉张损伤；相同采气时间下，采气速率越大，腔体围岩所受压力越小；此外，热力耦合作用下腔体围岩所受最大主应力的变化速率要高于仅有压力作用下的变化速率。

值得注意的是，压力作用下采气速率为120×104m3/d 时，最大主应力（压应力）的变化速率要略高于其他采气速率下的情况（见图12a）。 这是由于在不考虑热力耦合作用时，腔体应力主要受腔内气体压力影响。当采气速率较高时，腔内压力的下降速率加快，从而导致腔体压应力也逐渐增加。此外，在热力耦合作用下采气第2 天时，最大主应力曲线出现了较显著的拐点（见图12b）。 这是由于采气前，腔内气体处于相对稳定的状态，当以某个速率开始采气时，腔内气体的平衡被打破，气体压力和温度急剧变化，但由于数值模拟过程考虑的是以恒定速率进行采气，所以在随后的采气过程中，最大主应力值又趋于平稳。 而在压力作用下，初始采气时刻也存在轻微的拐点（见图12a），但由于仅考虑压力的作用，其波动幅度较热力耦合作用低很多。

5.3 剪切损伤

在压力作用下，腔体围岩只有在采气速率120×104m3/d、采气时间12 d 时，才有塑性区的出现（见图13）。该区域主要集中在腔体与泥岩夹层的交界处。热力耦合作用下时，腔体围岩塑性区出现在整个围岩，并且在腔体顶部、底部和泥岩夹层交界处均出现“shear-n”状态，表示在采气过程中，腔体在这3 个部位容易出现剪切损伤。和注气过程一样，热力耦合作用下的塑性区远远大于压力作用下的塑性区。

图13 采气过程中的腔体围岩塑性区Fig.13 Plastic zone of surrounding rock in cavity during gas production

5.4 膨胀损伤

采气速率120×104m3/d、采气时间12 d 时，腔体围岩分别在压力作用、热力耦合作用下的膨胀损伤情况如图14 所示。 由图可知，热力耦合作用下、采气速率120×104m3/d 时，腔体围岩的膨胀损伤判据出现了大于1 的区域。膨胀损伤主要集中在腔体侧部，出现在腔体围岩和泥岩夹层的交界处。

图14 采气过程中的腔体围岩膨胀损伤云图Fig.14 Expansion damage nephogram of surrounding rock in cavity during gas production

不同采气条件下的腔体围岩膨胀损伤判据值如图15 所示。 由图可知：压力作用下，随着采气时间和采气速率的增大，膨胀损伤判据值不断增加，但均小于1，因此不会出现膨胀损伤；热力耦合作用下，采气速率在40×104～80×104m3/d 时，随着采气速率的增加，腔体损伤判据值增加较小，但采气速率大于100×104m3/d 时，腔体损伤判据值在采气第8 天时开始快速增大，特别是在采气速率为120×104m3/d 时，膨胀损伤判据值的增大趋势非常显著，并最终超过上限值1，达到1.146。这是由于在膨胀损伤判据的表达式中，该值与腔体所受应力成指数关系。当采气速率较大时，腔体所受应力的变化速率也相应增大，最终导致了膨胀损伤判据值急剧增加。

图15 采气过程中的腔体围岩膨胀损伤判据Fig.15 Expansion damage criterion value of surrounding rock in cavity during gas production

综上所述，在采气速率40×104～120×104m3/d、采气时间0～12 d 的条件下，腔体围岩发生了拉张损伤、剪切损伤和膨胀损伤。 采气导致腔体内的温度和压力减小， 腔壁附近的围岩温度较更深部的围岩温度更低，腔壁围岩收缩，对腔体围岩产生向外的拉应力，当腔体围岩的内压和热应力大于地层压力时，腔体就会发生拉张损伤。 腔体围岩的剪切损伤大都发生在腔体和泥岩夹层的交界处， 这是由于泥岩夹层和盐岩的导热系数不同，温度变化引起腔体围岩变形不均，因此产生剪切损伤。 采气过程中，腔内温度减小，围岩应力张量第一不变量增加，主应力之间的差值减小，偏应力张量的第二不变量减小，从而产生膨胀损伤。

6 结论

1）基于变质量热力学，结合气体状态方程、系统热力状态方程、系统质量守恒方程和系统能量守恒方程，建立了地下盐岩储气库腔内气体温度和压力的计算模型。 考虑到运行过程中的稳定性判据主要由腔体是否发生损伤决定， 最终选定注采过程中腔体稳定性的判别标准为拉张损伤、剪切损伤和膨胀损伤。

2）在注气过程中，腔体围岩不发生拉张损伤和膨胀损伤，但是会发生剪切损伤，并且热力耦合作用下的腔体围岩剪切损伤区域远远大于压力作用下的剪切损伤区域。因此，在对地下盐岩储气库进行注气稳定性分析时，需要充分考虑温度的影响，并且着重关注岩体的剪切损伤情况。

3）在采气过程中，腔体围岩发生拉张损伤、剪切损伤和膨胀损伤。在采气速率120×104m3/d、采气时间12 d 时，热力耦合作用下腔体围岩所受最大主应力为拉应力（0.228 9 MPa），此时膨胀损伤判据值为1.146，大于上限值1。 与注气过程相比，采气过程中出现的腔体损伤更多，应该对采气过程中地下盐岩储气库腔体的稳定性重点关注。