ABAQUS中不同断层带的煤层地应力场反演

孔 鹏，任广聪，薄海江，马新仿，刘羽欣，高海滨

(1.中联煤层气有限责任公司 煤层气研发中心，山西 太原 030000； 2.中国石油大学(北京) 石油工程学院，北京 102249； 3.中海油田服务股份有限公司 油田生产事业部，天津 300000； 4.中联煤层气有限责任公司 晋城分公司，山西 晋城 048000)

引 言

近年来，我国煤层气开发取得了一定的成果，柿庄区块是我国煤层气勘探程度较高的煤层气区块之一，由于煤层自身性质的限制，水力压裂改造技术是煤层气高效开发的必须手段[1-2]。水力压裂形成的裂缝沿最大水平主应力扩展[3-4]，因此，全面认识煤层初始地应力场分布对于预测水力裂缝走向、压裂设计等具有重要意义[5]。通过Kaiser声发射实验[6]和古地磁实验[7-8]可以准确地获取地应力的大小和方向，但是室内实验只能测试有限位置的数据，可以根据这些有限数据点进行全区地应力反演，得到整个区块的地应力大小和方向分布。地层中普遍有一定规模的断层发育，断层的存在会影响地应力的连续分布状态[9-10]，因此，模型中对于断层处理方式的不同会对计算结果产生较大影响。

郭怀志[11]最早提出利用有限元数学模型回归分析方法求解初始地应力场，奠定了多元回归计算地应力场的基础，后来发展出偏最小二乘法、位移反分析方法、神经网络法和灰色关联法等。近年来，许多学者使用不同分析软件进行地应力场研究，包括ABAQUS、ANSYS、FLAC3D等，回归方法从多元线性回归发展到支持向量回归、多约束方法、p值法[12-14]。但是这些方法注重对算法的改进，忽略了储层自身的地质状况和岩石力学性质等的影响。

使用有限元方法建立的地质模型中，一般以断层带的形式处理断层，即在断层附近加密网格，赋予不同属性，目前使用的断层带处理方法主要包括单一性质和过渡断层带[12]。断层本身是一种较大规模的天然裂隙[15-16]，因此，可以考虑通过损伤力学模型来反映断层，利用ABAQUS中的粘结单元(cohesive)来模拟断层，考虑岩石的力学强度等性质，使结果更真实可靠。

本文使用大型有限元分析软件ABAQUS，建立地质模型，使用3种方法处理断层，分别为单一性质断层带、过渡断层带和粘结单元断层带，通过对3种方法结果的对比，分析各种方法的优劣性。

1 地应力反演求解原理

利用边界载荷反分析方法求解地应力场，如图1所示，将工区两边铰接固定，另外两边共分为m份，每次以相同大小的载荷进行加载，则一共加载m个载荷，在目标井位置，如果所有载荷的叠加效果与实际测得地应力大小和方向相同，则达到反演目的，从而获取整个工区的应力分布场。

每个目标井的应力约束条件包括最大水平主应力、最小水平主应力和最大水平主应力方向，但是在ABAQUS中无法直接将地应力方向作为约束条件输入，为方便计算，需要将主应力大小和方向转化为正应力S11、S22、S12，且拉力为负，压力为正[17-18]，具体关系如下：

(1)

(2)

(3)

其中：S11为X轴的正应力，MPa；S22为Y轴的正应力，MPa；S12为XY平面上的切应力，MPa；σH为最大水平主应力，MPa；σh为最小水平主应力，MPa；θ为X轴与最大水平主应力的夹角，(°)。

因此，假设共有n口目标井，则加载m个载荷时，每个载荷都会对目标井产生影响，m个载荷对1#井的应力影响为

(4)

同理，对第n口目标井有，

(5)

根据上述过程，m个载荷对n口目标井加载后，产生系数矩阵

(6)

可以看出，当m=3n时，该系数矩阵有唯一解，可以直接求解系数；当m<3n时，该系数矩阵有多解，可以通过多元线性回归分析，将常数项设为0，求取系数的大小。

2 模型建立

柿庄南区块位于沁水盆地，地层整体为单斜构造，局部发育有褶曲、断层和陷落柱，其中，断层以正断层为主。稳定发育的主要煤层为二叠系下统山西组的3号和石炭系上统太原组的15号2套煤层，研究工区为3号煤层。3号煤层埋深在600～750 m，呈南浅北深，中部埋深大于两侧的趋势，局部大于900 m，平均厚度6 m。具有储层压力低、渗透率低、饱和度低、非均质性强的特点。平均孔隙度为4.35%，渗透率介于(0.04～1.10)×10-3μm2，储层压力2.4～3.01 MPa，储层温度21.7～27.5 ℃。煤层(原煤)含气量一般在13～20 m3/t，含气量相对较高。

基于柿庄南区块3号煤层物性、岩石力学数据及煤层气压裂井分布，在ABAQUS中建立工区，大小为4 km×4 km，工区外围增加1 km宽度，消除边界效应的影响。选择well 1、well 2、well 3、well 4、well 5、well 6共6口井作为约束井，井眼位置如图2所示。工区中存在2个断层，分别为断层Ⅰ和断层Ⅱ。其中，well 1与断层Ⅰ相邻，well 2与断层Ⅱ距离较近。

约束井地应力大小和方向的准确性对反演结果至关重要。地应力与泊松比、杨氏模量等有关，又因本区块受构造运动的影响，水平主应力很大部分来源于地质构造运动产生的构造应力，因此，采用组合弹簧模型计算地应力剖面较为合适[19]。水力裂缝沿最大水平主应力方向扩展，因此，地应力方向根据邻井水力压裂微地震监测确定。约束井主应力大小、方向以及正应力、切应力如下表所示。

组合弹簧模型计算方法如下：

(7)

(8)

(9)

式中，σH为最大水平主应力，MPa；σh为最小水平主应力，MPa；β1为最大构造应力系数；β2为最小构造应力系数；μs为测井纵波时差，μs/m；pp为孔隙压力，MPa；η为有效应力系数；σv为垂向应力，MPa；ρb为岩石密度，g/cm3；ρ0为处理研究井段顶界至井口的地层密度平均值，g/cm3；H0为处理井段顶界深度，m；H为处理井段底界深度，m。

约束井主应力大小、方向及式(1)—式(3)计算的正应力、切应力如表1所示。

图3(a)为单一性质断层带，即在断层位置加密网格，并且赋予相同属性，整个模型中存在断层带和基质岩石2种材料；图3(b)为过渡断层带，即将断层区域分为过渡区和中心区，过渡区和中心区赋予不同属性，整个模型中存在过渡带、断层带和基质岩石3种材料；图3(c)为粘结单元断层带，即将断层区域分为过渡区和cohesive粘结单元区，利用粘结单元模拟断层行为，中间线条表示cohesive粘结单元。图3(a)和(b)是目前研究中使用的断层带处理方法，这两种方法虽然将断层和基质的性质区分开来，但是并没有考虑断层带的力学行为，因此，建立图3(c)的cohesive粘结单元断层带。

如图4所示，将模型划分网格，岩石基质和断层过渡区域网格类型均为渗流/应力单元(CPE4P)，cohesive粘结单元网格类型为孔隙应力粘结单元(COH2D4P)。使用过渡网格，工区采用小尺寸网格，断层区域局部加密，提高计算精度；边界区域采用大尺寸网格，减少计算时间。由于断层带处理方式不同，整体网格数量有微弱差异。

使用渗流/应力单元可以考虑地层孔隙和地层流体的影响，粘结单元从力学损伤角度考虑更加符合断层自身性质。为了控制岩石力学之外的变量对结果的影响，孔隙度取5%，渗透率取0.1×10-3μm2。其中，岩石基质杨氏模量、泊松比数据来源于室内三轴压缩实验结果，过渡区域取值实验结果的0.8倍，中心区域和粘结单元取值实验结果的0.6倍。根据室内三轴实验结果，柿庄南煤岩杨氏模量为4.13～8.03 GPa，泊松比0.35～0.44，因此，模型基质取值如表2所示。

表2 材料属性Tab.2 Properties of rock

3 模型求解

(1)计算过程

以单一性质断层带为例展示求解过程，如图5所示，将工区两边铰接固定，另外两边共划分为8等份，每次施加20 MPa的压力，并分别记录8次施加应力后目标井位置的S11、S22、S12值(表3)。

表3 单次载荷计算结果Tab.3 Calculation results of loads

根据式6，使用多元线性回归方法计算每个载荷对应的系数，计算结果见表4。将校正后的8个载荷同时加载，计算结果即为初始地应力场分布。

表4 线性回归结果Tab.4 Linear regression results

同理，可以求得过渡断层带、粘结单元断层带2种处理方式下的载荷系数。

(2)计算结果

图6分别为单一性质断层带、过渡断层带和粘结单元断层带3种模型的最大水平主应力计算结果。可以看出，最大水平主应力大小在断层带区域及附近存在差异，尤其是断层端部，应力大小与附近区域存在突变，表现出非连续性的分布。断层位置应力值小于周围基质区域，断层Ⅰ处最大水平主应力一般小于10 MPa，而周围基质区域大于10 MPa。另外，过渡断层带模型中，虽然应力分布在断层带附近出现突变，但是不会出现错位，即应力大小连续分布，在断层处出现“间断”；单一性质断层带和粘结单元断层带则会出现分布错位。

图7分别为单一性质断层带、过渡断层带和粘结单元断层带3种模型的最小水平主应力计算结果。可以看出，3种模型最小水平主应力大小总体一致，在断层Ⅱ下端，应力分布出现部分差异。断层位置最小水平主应力偏小，一般小于6 MPa。另外，过渡断层带模型中，应力分布在断层带附近只出现突变，但是不会出现错位。

图8和图9分别为单一性质断层带、过渡断层带和粘结单元断层带3种模型的最大水平主应力、最小水平主应力方向的计算结果。

3种模型计算的应力方向一致，最大水平主应力方向与最小水平主应力垂直，断层处应力方向发生突变，最小水平主应力方向与断层走向一致。

(3)结果分析

为了确定不同断层带对计算结果的影响，通过后处理分别读取3种情况下断层Ⅰ和断层Ⅱ附近的well 1井和well 2井水平应力大小和方向，并与实际测试值对比，结果如表5和表6所示。a表示单一断层带，b表示过渡断层带，c表示粘结单元断层带。

表5 水平应力计算结果与实际测试值对比Tab.5 Comparison between calculating and testing horizontal stress

表6 水平应力计算值与实际测试值之差Tab.6 Difference between calculating and testing horizontal stress

可以看出，3种断层带处理方法计算的最大水平主应力和最小水平主应力绝对值比较接近，与实际值均存在误差。

其中，断层Ⅰ附近的well 1井水平主应力绝对值偏小，水平应力差值也偏小。过渡断层带处理方式的水平应力差计算结果最接近实际值，偏小0.12 MPa；单一断层带计算结果与实际值相差最大，偏小0.34 MPa；粘结单元断层带计算结果次之，偏小0.25 MPa。

断层Ⅱ附近的well 2井水平主应力绝对值偏大，水平应力差值也偏大。其中，单一断层带的水平应力差计算结果最接近实际值，偏大0.35 MPa；过渡断层带与实际值偏差最大，偏大0.49 MPa；粘结单元断层带计算结果次之，偏大0.38 MPa。

从3种断层带处理方法计算出的地应力方向可以看出，粘结单元断层带结果最接近实际方向。Well 1井计算结果为粘结单元断层带最准确，单一断层带次之，过渡断层带方向偏差最大；well 2井计算结果为粘结单元最准确，单一断层带次之，过渡断层带方向偏差最大。因此，可以判断，粘结单元计算的方向最准确。

4 结 论

(1)在断层带尤其是断层端部，应力大小与附近区域存在突变，表现非连续性分布，地应力方向在断层带发生突变。

(2)计算地应力大小时，粘结单元断层带计算结果的误差均为中等，单一断层带和过渡断层带计算结果不稳定；计算地应力方向时，粘结单元计算的水平应力方向最准确，单一断层带次之，过渡断层带计算的水平应力方向偏差最大。

(3)同时考虑地应力大小和方向误差，粘结单元断层带处理方法计算方向最准确，其原因为粘结单元将断层视为天然裂隙处理，考虑损伤演化。模拟计算煤层地应力场分布时，应关注断层等地质构造的岩石力学性质，并将之与岩石基质性质区分开来，使结果更符合实际情况。