我国天然气储备量存在较大的地域差异性，使用量存在典型的季节差异性，北方储量远大于南方储量，冬季用量远高于夏季用量，地下储气库(Underground Gas Storage，UGS)作为天然气峰调的重要手段，对我国天然气战略具有重大意义。枯竭油气藏储气库和盐岩储气库均对地域要求较高，依赖于大型油气田或者地下盐矿，而含水层储气库选址范围较大，可以作为一种重要的选择类型。含水层储气库依据天然断层的空间分布，构造多孔介质储层来储蓄天然气，通过固井技术实现周期性采注，运行过程中多孔介质的流体压力增大，固体压力减小，储层膨胀会使断层部位(包括不活动断层)的有效应力发生改变，容易开启或破裂，导致天然气泄漏。因此，断层安全性是含水层储气库研究的重要课题。

断层安全性评价经历了只针对泥岩对接和泥岩涂抹层的定性评价(泥岩涂抹因子法SSF^[1]、涂抹泥灰比法SGR^[2-3])、半定量的Knipe图解法^[4]和利用断层—储层排替压力差法^[5-7]的精细化定量评价三个阶段。常用的断层评价方法均没有从整体考虑注采荷载对断层安全性的影响，并且需要大量实验数据支持，因此，以数值模型为基础，模拟储气库断层应力状态的研究逐渐增多。例如：赵斌等^[8]采用气水二相渗流理论及数值模拟方法，结合储气库构造参数建立了含水层储气库采注效应模型；危齐等^[9]构建呼图壁地下储气库的ANSYS有限元模型，ZHANG等^[10]建立格洛斯特盆地简化三维储层模型，研究注采压力对储气库的影响；HUANG等^[11]提出三维离散断裂网络—断层(3D DFN-fault)模型，考虑了断层滑动引起的剪切位移和正常位移。王则等^[12]在对断层破碎带进行识别和发育模式划分的基础上，利用测井资料、三轴应力试验获得岩石力学参数，基于有限元法模拟了3种断层破碎带发育模式。谭丽娟等^[13]以断层体为研究对象，提出了断层封闭综合指数的概念，用断面正压力、断层带泥质含量、断层带碎屑物成岩胶结作用三大影响因素量化乘积表示。赵乐强等^[14]研究了准噶尔盆地西北地区断层内流体活动过程及对断层启闭性的影响。景紫岩等^[15]基于相似试验提出了基于断裂带泥岩连续涂抹断层封闭性评价参数新算法。尽管国内外对断层模型的研究不断深入，但是对含水层储气库断层的研究尚浅，并未形成一套完整的断层安全性评价体系。

基于此，笔者根据含水层储气库的特性，借助ANSYS有限元模拟软件提出含水层储气库断层安全性评价方法，以3D实体模型来模拟含水层储气库在不同运行压力条件下断层处的应力状态，构建合理的断层安全性判据，依据静力学分析方法评价断层的安全性。另外，合理预测含水层储气库的最大运行压力可以有效防止不活动断层发生失效破坏^[16]，以苏北盆地东台坳陷大丰—兴化探区白驹储气地质体为例，构建3D实体模型，在合理评价目标断层安全性的同时，基于断层安全性判据预测其最大运行压力，为储气库的安全运行提供参考依据。

1.   ANSYS断层安全性评价方法

1.1   断层安全性评价思路

ANSYS软件是美国ANSYS公司研制的一款大型通用有限元分析(Finite Element Analysis，FEA)软件，包括前处理模块、计算模块和后处理模块三部分，通过构建研究对象的3D实体模型，既可以进行瞬态力学分析，也能够模拟实体模型在不同荷载条件下的静态力学状态。由于含水层储气库注采过程较短，为了研究断层运营期的长期稳定性，需借助ANSYS软件的前处理模块构建3D实体模型，计算模块模拟若干个不同内压下断层应力状态，后处理模块实现断层的静态力学分析，探究运行压力对断层主应力的影响(图 1)。

图  1  断层安全性评价思路

Figure  1.  Diagram for fault safety evaluation

下载: 全尺寸图片 幻灯片

借助ANSYS软件构建出有效的3D实体模型是目标断层安全性评价的基础，因此有必要对3D实体模型初始平衡状态下地层、储层和断层的应力状态进行有效性验证，避免模型失真。此外，如何对已建模型施加实际运行荷载，从而模拟不同储层内压；如何根据断层处应力状态判断断层安全性，从而预测运行压力的极限值，均是实际应用中需要解决的重要问题。笔者将以白驹储气地质体实际地层和断层的探测数据为基础构建ANSYS有限元模型，具体从模型构建、模型有效性验证和断层应力模拟及分析等方面展开研究。

1.2   断层安全性评价标准

含水层储气库运行过程中，断层处应力状态是影响断层安全性的主要因素，也是断层安全性评价的重要依据。国内外不少专家和学者针对断层应力状态的研究提出大量基于岩石力学和流体力学的计算方法，并分析应力与断层封闭性定量关系。王珂等^[17]推导出了三向地应力状态下断面正应力和剪应力的计算公式，并分析了地应力大小和方向对断层封闭性的影响。CUSS等^[18]对饱和高岭石和球黏土断层泥断层复活电位进行实验研究，发现一旦凿槽内的孔隙压力足以克服作用于断层上的正应力，就会发生滑动。YIN等^[19]对高发育走滑断裂带地应力非均质性进行三维有限元模拟，验证了断层规模、断层密度、断层走向和断层带效应对沿走滑断裂带的应力非均质性有重要影响。

断层应力状态影响断层的紧闭系数，紧闭系数越大，断层封闭性越强，而断层应力对封闭性影响主要体现在正应力性质、应力非均质性和采注应力变化上。对于多孔介质储层而言，注采之前，断层主要承受的应力为地层施加的压应力；注气之后，断层内会出现孔隙压力，对断层的作用表现为拉应力。不断向储层注气，内部孔隙压力逐渐增大，孔隙压力对断层作用的拉应力与地层压力对断层作用的压应力相互抵消，导致断层所受的有效应力降低。随着储层内压不断增大，最终断层所受的压应力会被运行所产生的拉应力完全抵消，断层将达到开启或者破裂的临界状态。因此，将断层的这种临界状态作为断层安全性评价的标准，即“断层出现拉张应力”，通过分析不同运行压力条件下断层应力云图，找出断层处X/Y/Z方向应力中最先出现拉应力时所对应的运行压力，作为最大运行压力的预测值，超过这个值断层可能发生拉张破坏，造成断层开启或者破裂。

2.   ANSYS模型实例

白驹储气地质体位于苏北盆地东台坳陷东北部的大丰—兴化探区，自晚白垩世以来构造演化经历了晚白垩世区域坳陷成盆期、古近纪幕式断陷成盆期和新近纪以来的整体沉降坳陷成盆期，其古近纪拉张断陷成盆期分为初始断陷—断陷—坳陷、强烈断陷、收敛以及抬升剥蚀等4个幕式阶段。依据岩性、地震反射特征，古近系可划分为T₂⁰、T₂³、T₃³、T₃⁴、T₃⁵、T₄⁰六个关键地层界面，其中最优储气地层位于T₃⁵面和T₄⁰面之间，为含水构造多孔介质储层。储层断层发育较多，且断层面起伏较大，根据实际地质资料对基础数据进行筛选，不断优化模型，构造出有效的3D模型。

2.1   大丰—兴化探区地质概况

大丰—兴化探区北抵建湖隆起，西临高邮凹陷，南靠溱潼凹陷和海安凹陷，东边与小海凸起和裕华凸起接壤。该地区的勘探主体为白驹凹陷，受几条近东西向大断层的控制，形成凹凸相间的构造格局，自北向南依次为大丰次凹、施家舍断裂带、洋心次凹、草埝断阶带、吴堡低凸起和冯家墩次凹(图 2)。大丰—兴化探区地层非常破碎，构造复杂，火山岩较发育，岩性横向变化非常快。主要沉积了一套远端辫状河三角洲环境的地层，岩性为灰色厚层粉砂岩、含砾砂岩、砂砾岩与深灰色、灰黑色泥岩互层，而在凹陷深凹地区推测为半深湖—深湖沉积环境。

图  2  苏北盆地东台坳陷大丰—兴化探区构造

Figure  2.  Structural map of Dafeng-Xinghua exploration area, Dongtai Depression, Subei Basin

下载: 全尺寸图片 幻灯片

白驹凹陷在充填过程中由于受应力场转换以及沉积作用间断和湖平面突变等事件的影响，发育了各级不同的并在凹陷内可追踪对比的等时地层界面。笔者依据岩性和地震反射特征，运用高分辨率层序地层学方法，对大丰—兴化探区进行了层序界面划分，为后续建模提供实践地质参数，以建立更加准确的三维地质模型，避免模型失真。通过仔细追踪和对比，在白驹凹陷古近系共识别出6个具对比意义的等时地层界面(图 3)。

图  3  苏北盆地东台坳陷大丰—兴化探区层序地层综合划分柱状图

据文献[20]修改。

Figure  3.  Comprehensive division of sequence stratigraphy in Dafeng-Xinghua exploration area, Dongtai Depression, Subei Basin

下载: 全尺寸图片 幻灯片

2.2   数据处理及模型创建

通过Surfer15.0软件地质建模中的插值函数对坐标点进行筛选，以地层T₂⁰为例，其数据筛选步骤如下所示：

(1) 将源文件导入到Surfer15.0软件，选定插值函数及X、Y坐标间隔；

(2) 对这些坐标点自动筛选、排序，并且可以直接绘制出地质等高线图，生成三维地质曲面(图 4)；

图  4  地层T₂⁰数据处理效果

Figure  4.  Processing of data of stratum T₂⁰

下载: 全尺寸图片 幻灯片

(3) 输出地质曲面的数据点坐标，达到数据筛选的目的；

(4) 通过Surfer to ANSYS程序将坐标点编译成ANSYS软件可识别格式。

Surfer15.0软件中的插值函数筛选得到的地质曲面较大程度吻合实际地层面和断层面，对实体模型的精度影响不大。采用该方法，筛选出T₂⁰、T₂³、T₃³、T₃⁴、T₃⁵和T₄⁰共6个地层界面的有效坐标点，并选出9条需重点关注的目标断层，按照“点点成线、线线成面、面面成体”的基本顺序在ANSYS软件中创建出初步地层—断层实体模型(图 5)。

图  5  初步地层—断层实体模型

Figure  5.  Preliminary stratum-fault entity model

下载: 全尺寸图片 幻灯片

2.3   实体模型调整及优化

初步模型的地层多是由小三棱柱体组成，虽然能够较好地反映地层起伏趋势，但过于细化，不利于网格划分，需对地层进行适当调整，将有断层贯穿的地层设计成一个规则体，并保留地层的实际起伏状态。另外，为了更真实地模拟储气地质体地质条件，将T₂⁰往上延伸得到上覆层，T₄⁰底往下延伸得到下覆层。

初步实体模型的断层多为半贯穿状态，后续平面切割和网格划分存在较大难度，需在不改变断层起伏状态的前提下适当调整半贯穿断层，即调整断层不明显倾角和起伏较大的边际点(图 6)。形成图 7所示的地层—断层实体模型，该模型东西方向长14 979.79 m，南北方向长10 769.28 m，占地面积为161.32 km²，埋深范围为650~2 480.13 m，厚度为1 830.13 m。研究区域的断层数为9条，地层包括上覆层、T₂⁰—T₂³层、T₂³—T₃³层、T₃³—T₃⁴层、T₃⁴—T₃⁵层，T₃⁵—T₄⁰层(储层)和下伏层。

图  6  优化后断层模型

Figure  6.  Optimized fault model

下载: 全尺寸图片 幻灯片

图  7  优化后地层—断层模型

Figure  7.  Optimized stratum-fault model

下载: 全尺寸图片 幻灯片

2.4   参数设置及网格划分

模拟运行荷载前，需对已建地层—断层实体模型进行属性和材料的定义，确定各项参数的赋值，比如弹性模量、泊松比、密度等。结合相关资料和地震剖面的岩性特征，对地层及断层的物理力学参数进行赋值，参数值来源于岩石力学实验及相关参考文献^{[9, 21]}，见表 1。

表  1  模型参数

Table  1.  Model parameter

目标层		单元属性	弹性模量/GPa	泊松比	密度/(kg·m-3)
上覆层(盖层)		SOLID185实体单元	47.82	0.32	2 070
地层	T20—T23		46.62	0.30	2 060
	T23—T33
	T33—T34
	T34—T35
储层	T35—T40		44.90	0.28	2 050
下伏层			45.92	0.33	2 090
断层	F008/F022/F023/ F026/F063/F072/ F080/F083/F115		47.62	0.33	2 090

下载: 导出CSV 

| 显示表格

网格划分是模型求解的关键，是要将已建好的地质模型划分为一个个可以定义属性、材料、参数等的微小单元，后续对这些单元求解计算来反映储气库的应力状态。已建模型网格划分具有以下几个特点：(1)地层和断层统一定义为实体属性，均采用SOLID185实体单元，改变实体单元材料的属性即可区分地层、储层和断层；(2)为提高网格划分的精度，提前将模型相关线段进行划分，划分大小为100 m；(3)中间区域存在许多贯穿断层，采用Free规则自由划分，其他区域均采用Sweep规则体进行划分(图 8)。

图  8  模型网格划分效果

Figure  8.  Division effect of model grid

下载: 全尺寸图片 幻灯片

3.   模型有效性验证

3.1   初始应力平衡分析

通过模型重力加载来模拟实际地质条件的初始应力状态，即约束模型边界，Z方向施加模型自重求解初始平衡后的应力场。一方面可以验证已建模型初始平衡后应力数值是否符合实际地质条件；另一方面在模型初始应力平衡后的基础上再次对模型进行不同内压的加载计算，能够有效模拟运行之后储气库的应力情况。地层—断层初始平衡后X/Y/Z方向的应力状态如图 9所示，应力场呈层状分布，随着埋深的增加而不断增大，其中储层对应的X、Y方向的应力值在20 MPa左右，Z方向的应力值在40 MPa左右，与储气库运行之前储层的实际内压相符，表明已建模型可以模拟实际运行之前的真实情况。另外，储层部分对应的断层X方向应力值(第一主应力)在17 MPa左右，Y方向应力值(第二主应力)在17.3 MPa左右，Z方向应力值(第三主应力)在40 MPa左右，均处于正常范围。

图  9  模型初始平衡后地层和断层X/Y/Z轴方向应力云图

Figure  9.  Stress nephogram along X/Y/Z axis of strata and fault after initial equilibrium of the model

下载: 全尺寸图片 幻灯片

3.2   岩石轴压实验验证

为进一步验证模型初始平衡模拟结果的准确性，抽取丰探1号井埋深分别为1 426.50和1 630.35 m的岩石样品做单轴压缩试验，测定岩石的单轴抗压强度(RC)分别为35和40 MPa，而初始平衡对应埋深的应力为31.1和35.4 MPa，小于岩石的单轴抗压强度，即岩石不会发生破裂。另外，实际储气库地层岩石会受到不同围压的作用，为了验证围压对岩石抗压强度的影响，抽取丰探9号井埋深2 776.30 m的岩石样本，制作A、B试样，增加不同的围压(2.0，4.0 MPa)进行三轴加载试验(图 10)，试验表明：相同区域和埋深的岩石随着增大围压，抗压强度也随着增加，并且该试样的抗压强度均超过90 MPa，而模型初始平衡后对应埋深的应力只有50 MPa左右，远小于岩石的抗压强度，表明已建模型能够较准确地模拟实际地质情况。

图  10  岩石轴压试验应力曲线

Figure  10.  Stress curves of rock axial compression test

下载: 全尺寸图片 幻灯片

4.   运行过程断层安全性评价

4.1   运行荷载加载原理

含水构造的储层部分不是常见的枯竭油气藏或者大型的地下盐穴，而是特殊的多孔介质储层。由于储层不是空腔，ANSYS软件无法直接施加面荷载来模拟储层内压，经过综合考虑，决定用“位移约束”来模拟“储层内压”。

对于多孔介质储层而言，向孔隙中注入气体后产生膨胀，对固体介质产生了压力。注气之前，储层内的地层压力(P)只由固体压力(P_s)组成；注气之后，储层内的地层压力变为由气体产生的流体压力(P_f)和固体压力共同组成，并且二者之间存在着固定的转化方式。由于地层压力总体保持不变，当不断向储层注气时，储层的流体压力增大，固体压力减小，储层内的固体介质就会向外膨胀，产生相应的位移量。

初始平衡：

注气之后：

不断注气：

式中：P为地层压力(保持不变)，MPa；P_s0为初始平衡储层内压(固体压力)，MPa；P_s为固体压力，MPa；P_f为流体压力(即注气压力)，MPa；n为孔隙率，依据实际取值20%；ΔP_s为固体压力减小量，MPa；E为储层的弹性模量，MPa；μ为泊松比；ε为产生位移的应变。

取特定的注气压力，通过公式(1)、(2)计算相应固体压力减小量，而弹性模量和泊松比是已知，由应力—应变公式(3)可以求出产生位移的应变值，再由应变乘以储层的厚度就可以得到储层产生的竖直位移量，由应变乘以储层水平方向的长度就可以得到储层产生的水平位移量，由应变乘以储层的宽度就可以得到储层产生的横向位移量。基于ANSYS软件计算模块功能，通过“位移约束法”约束多孔介质储层上下表面各节点相应的位移量，求解计算，即可模拟特定注气压力下模型的应力状态。

4.2   断层应力状态分析

经过综合计算，取特定注气压力梯度(3.15 MPa)对已建模型进行加载模拟。即运行压力分别取3.15，6.30，9.45，12.60，15.75，18.90，22.05，25.20，28.35和29.90 MPa，对模型储层上下表面约束相对于位移量，求解计算得出断层在X/Y/Z轴方向的应力状态，如图 11所示，为运行压力在28.32和29.90 MPa条件下，目标断层的应力云图。

图  11  运行压力28.32 MPa和29.90 MPa条件下断层X/Y/Z轴方向应力云图

Figure  11.  Stress nephogram along X/Y/Z axis direction of fault at running pressure of 28.32 MPa and 29.90 MPa

下载: 全尺寸图片 幻灯片

从第一主应力图来看，可以发现与储层接触的上下段断层的第一主应力呈现一定的减小规律，从数值上进行分析得出每增加3.15 MPa左右的运行压力，与储层接触的上下段断层的第一主应力会相应减小0.46 MPa左右，整体上是呈线性减小的规律。

从第二主应力图来看，可以发现与储层接触的上下段断层的第二主应力没有明显的变化规律，从数值上进行分析得出随着运行压力的增大，与储层接触的上下段断层的第二主应力值会出现小幅度的增大。

从第三主应力图来看，可以发现与储层接触的上下段断层的第三主应力基本保持着稳定的状态，从数值上进行分析得出随着运行压力的增大，与储层接触的上下段断层的第三主应力值基本维持在-12 MPa左右，没有出现明显的增减。

模型加载运行荷载时，储层部分需要进行开挖处理之后，再对储层上下表面的节点约束相应位移量模拟储层的膨胀，会导致处于储层部分的断层两端向外拉，逐步抵消来自模型自重产生的压应力，与储层部分接触的上下段断层会更早地出现拉应力。因此，记录不同运行压力条件下与储层部分接触的上下段断层最小压应力值，并绘制变化图(图 12)，即可判断目标断层的安全性。

图  12  不同运行压力下断层最小应力变化

Figure  12.  Minimum stress variation of fault under different operating pressure

下载: 全尺寸图片 幻灯片

与初始平衡状态相比，断层的各个主应力明显减小。断层X轴方向最小压应力呈线性规律减小，断层Y轴方向最小压应力先是缓慢增长，增长幅度并不大，当运行压力达到15.75 MPa后开始出现小幅度减小，28.35 MPa后逐渐趋近于平缓；断层Z轴方向最小压应力一直都保持着一个稳定的值。含水层储气库运行压力的增大，对断层水平方向应力的影响最大，对断层竖直方向应力的影响次之，对断层垂直方向应力的影响最小，并且当断层水平方向应力出现拉应力时，断层竖直和垂直方向远远没有出现拉应力。

4.3   最大运行压力预测

考虑压力对断层水平方向应力场的影响，以此来判定运行压力的极限值，当运行压力达到28.35 MPa时，断层F008与储层接触的上下段的最小压应力值达到0.17 MPa(图 13)。再增加一个压力梯度，即取运行压力29.9 MPa，断层F008与储层接触的下端出现拉应力，此时拉应力值为0.07 MPa(图 14)，即断层F008可能出现拉张破坏，断层有开启或破裂的风险。

图  13  运行压力28.35 MPa条件下断层F008的水平方向应力云图

Figure  13.  Horizontal stress nephogram of fault F008 with operating pressure of 28.35 MPa

下载: 全尺寸图片 幻灯片

图  14  运行压力29.90 MPa条件下断层F008的水平方向应力云图

Figure  14.  Horizontal stress nephogram of fault F008 with operating pressure of 29.90 MPa

下载: 全尺寸图片 幻灯片

因此，判断运行压力的极限值取值范围在28.35~29.90 MPa，经过区间等比例置换计算出运行压力的极限值为29.50 MPa，为保证运营期断层的安全性，建议其运行压力不能超过29.50 MPa。

5.   讨论与展望

断层性质及其封闭性受控于断裂面充填物性质等地质条件，要想采用本方法进行断层的封闭性分析，首先要对断层、目的地层的本构模型获得清晰的认识，进而分析断裂面充填物性质的影响。从安全角度考虑，假如断层未受到储气库建设扰动前处于封闭状态，那么从断层是否出现拉应力判断其封闭性是否受到破坏是个比较好的选择。注气的工作制度、生产参数，如注气压力、流量等对断层封闭性也有影响，这些因素通常只能转换成等效应力才能考虑按照本方法进行分析。此外，本方法主要是从静力的角度进行分析，由于注气压力等变化通常很缓慢，可以认为是一个准静态变化过程，采用本方法分析依旧具有一定的可行性。

6.   结论

本文以苏北盆地东台坳陷大丰—兴化探区白驹储气地质体为研究对象，旨在探究运营过程中断层的安全性，以3D实体模型的形式模拟不同运行压力下断层的应力状态，通过分析得出如下结论：

(1) 基于ANSYS软件的静态力学分析构建了一套完整的断层应力模拟、分析和评价体系，并以实例论证了其可行性。首先，运用ANSYS软件进行数据测试、参数赋值、网格划分等工作，构建完整的3D实体模型，并通过初始应力平衡计算和断层岩石力学实验的结果验证了模型的有效性。其次，根据实际情况确立位移约束的荷载施加方案，并以3.15 MPa为梯度模拟出不同运行压力条件下断层处的应力状态。最后，考察与储层部分接触的上下段断层的最小压应力值，以断层出现“拉张应力”为安全性判别标准进行断层安全性评价。

(2) 通过考察不同运行压力条件下断层在水平、竖直和垂直方向的最小压应力的变化情况，可以得出结论：运行压力对断层水平方向应力影响程度较大，对断层竖直和垂直方向影响程度较小；随着运行压力的不断增大，断层水平方向有开启或者破裂的风险，当运行压力超过29.50 MPa时，断层水平方向将会出现拉应力，可能发生拉张破坏，影响断层的安全性。因此，为保证运营期断层的安全性，建议实际运行压力不超过29.50 MPa。