恒泰艾普软件培训系列教材-应力场模拟和裂缝预测

恒泰艾普软件培训系列教材Stress Field Modeling应力场模拟和裂缝预测恒泰艾普石油勘探开发技术有限公司目录一、 模块功能二、 原理和方法三、 参数和使用说明四、 应用关键和应用技巧一、模块功能利用地质、钻井和测井资料，计算拉梅常数和剪切模量等参数，建立地质模型、力学模型及数学模型，运用三维有限差分数值模拟方法对应力场进行模拟，研究构造、断层、地层厚度、区域应力场等地质因素与裂缝分布的关系，预测与构造有关的裂缝分布及发育程度。

二、原理和方法（一）应力场的概念地壳中或地球体内，应力状态随空间点的变化，称为应力场，或构造应力场。

应力场一般随时间变化，但在一定地质阶段相对比较稳定。

研究应力场，就是研究应力分布的规律性，确定地壳上某一点或某一地区，在特定地质时代和条件下，受力作用所引起的应力方向、性质、大小以及发展演化等特征。

随着地质演化，一个地区常常经受多次不同方式的地壳运动，导致同一地区内，呈现出受不同时期不同形式地应力场作用所形成的各种构造极其叠加或改造的复杂景观。

因此，只有最近一期地质构造事件，未经破坏或改造，才能确切地反映这个时期的应力场。

应力场可按空间区分为全球、区域和局部地应力场；按时间区分为古地应力场和今地应力场；按主应力作用方式区分为挤压、拉张和剪切地应力场。

（二）地质模型和应力场关系地支模型的建立是做好应力场模拟的先决条件，首先将储层的目的层连同上下盖层和覆盖层作为一个岩石块体的隔离体来计算，然后从地质的角度提出构造成因，构造裂缝的特征，构造应力场的宏观特征及断层发育史。

我们现在研究的构造应力场主要在早白垩世构造伸展期与晚白垩世构造反转期形成，因此研究的地质体应为相应时期的古构造图。

对于挤压构造，应取受挤压之前的古构造作为地质体；而对于伸展构造，考虑到伸展作用的长期性及伸展对构造缝所形成的控制作用，应取伸展之后的古构造作为地质体。

在此基础上，恢复古构造剖面图，推断地质隔离体的受力方向及大小，设定边界条件并提出反演应力场及裂缝的地质标准。

但是，由于地质体是一个十分复杂的地下岩石块体，其地壳中各种地质构造形态，类型，成因是在漫长历史时期的地质演化过程中形成的，这种复杂的地质演化过程不可能恢复，只能用相对静止的观点和相对简化的方法去处理构造与古应力场的问题。

考虑到储油构造中主要的一种类型是背斜构造，它是由地壳受到挤压发生弯曲，或由于基底隆起使沉积地层上拱而形成。

在从历史时期的连续的接受应力而形成现代的构造体系，我们可以近似的以现在的构造应力场来作为古构造应力场的一个发展模型，来模拟在这样的应力场的条件下会导致怎样的构造裂缝体系（三）基本原理本软件以弯曲薄板作为油层构造模拟的力学模型，用二维的方法来处理油层构造，这种方法计算方便，人工干预少，对于应力场模拟可以将边界作为自由边界处理，不需要另外考虑边界条件。

该方法主要以背斜构造作为模型进行分析，得出可以用构造面上一点的最大曲率值作为该点裂缝发育程度的判据，而以最小主曲率方向指示可能出现的张裂缝走向，这样就将构造裂缝的分布问题化为构造面的主曲率计算问题。

设以薄板中面为z=0的坐标面，规定按右手规则，以平行于大地坐标为X，Y坐标，以向上为正。

沿X ，Y 正方向的位移分别为u x ，u y ，沿Z 方向的位移为扰度w(x,y)。

在直角坐标系中，如图：X,u xY , u y图 1，薄板模型示意图1. 基本方程直角坐标系中的变形几何方程：wu zu y ux u z u z u x u y u y u x u z y z yz z x xz z z y x xy y y x x =∂∂+∂∂=∂∂+∂∂=∂∂=∂∂+∂∂=∂∂=∂∂=)(),(),(,,γγεγεε (1) 由薄板理论可知，有：yw z u x w zu y x ∂∂=∂∂=, (2) 且有：2222,,2x y xy w wz z z 2w x y x εεγ∂∂===∂∂∂y∂∂ (3)定义曲率变形分量：y x wyw x w xy y x ∂∂∂−=∂∂−=∂∂−=22222,,κκκ (4) 因此，应变分量可写为：xy xy y y x x z z z κγκεκε2,,−=−=−= (5)2. 物理本构关系(广义虎克定律)yzyzx y z z xz xz z x y y xyxy z y x x EEE E E E τνγσσνσετνγσσνσετνγσσνσε)1(2)],([1)1(2)],([1)1(2)],([1+=+−=+=+−=+=+−=(6)其逆关系为：kkxzxz z z yz yz y y xy xy x x G G G G G G εθγτλθεσγτλθεσγτλθεσ==+==+==+=,2,2,2 (7)λ，G 为 拉梅(Lame)常数，G 也就是剪切模量(Shear modulus)，E 杨氏模量(Yong modulus)。

θ为体积应变。

将前面的式代入，得到：xy xy x y y y x x G v E v E γτεενσεενσ1),(1),(122=+−=+−=(8a) 因而有：xy xy xy x y y y x x EzG v E v Ez κνκτκκνσκκνσ)1(2),(1),(122+−=−=+−−=+−−= (8b)将地层厚度t=2z 代入上式，得到由曲率分量表示的地层面上的应力分量：xy xy x y y y x x Etv Et v Et κντκκνσκκνσ)1(2),()1(2),()1(222+−=+−−=+−−= (9) 由上式可知，当地层面向上凸时，曲率大于零，正好对应上凸地层面受拉张应力，张应力为正。

为了与地质力学符号相符，这里采用压应力为正，张应力为负的符号约定。

曲率小于零，表示地层上凸。

求出该点的沿坐标的应力后，就可求出其主应力及其方向：22min 22max )2(2,)2(2xyyx yx xyyx yx τσσσσστσσσσσ+−−+=+−++=(10)max σ与X 轴的夹角α，min σ与X 轴的夹角β：xy xg t τσσα−=max )(, yxy g t σστβ−=min )( (11)因此，若能得到地层面的扰度方程或其面上点的曲率，就可以估算其上的应力场，进而计算由此应力产生的裂缝。

3. 地层曲率的计算 地层趋势面的拟合由前面理论可知，若能求出地层面的曲率分量，就可以求出其上的应力场。

采用趋势面拟合方法拟合地层面的趋势函数，进而计算其上点的曲率分量。

采用最小二乘法拟合趋势面。

设趋势面的待定系数的函数为：25423210),(y a xy a x a y a x a a y x w +++++= (12)由层面散点处的坐标值(x,y,z),，建立最小二乘方程，对一个散点：),(i i i i i y x w z −=ε (13a)5,4,3,2,1,0;02==∂∂j a ji ε (13b)当用n 个散点拟合一个趋势面时，可得到拟合方程组，解此方程组，： 其中求和号表示，即对1..n 点求和。

解此线性方程组，就可得到趋势面函数。

∑ni趋势面的曲率计算式：42522322,2,2a y x wa yw a x w xy y x −=∂∂∂−=−=∂∂−=−=∂∂−=κκκ (14)4. 裂缝参数计算1) 曲率参数由上述解方程组可得到拟合曲面系数a 3，a 4，a 5，由(14)式可得到该点处的曲率。

2) 应变参数由(5)式可得到应变值。

其中，z=t/2。

再由(9)，(10)，(11)可分别计算出相应的应力，主应力和主应力方向。

3) 应力参数由(9)，(10)，(11)可得到主应力及其方向。

（四） 应力场应用流程FRS 系统提供的应力场数值模拟，考虑了储层岩石的厚度、岩性，并且考虑了储层受构造控制的裂缝分布等因素，基于弹性薄板理论，在计算构造的曲率分布基础之上，进一步计算构造的应力、应变分布，然后根据构造的应力应变场，就可以对储层裂缝的发育程度及展布关系进行分析。

构造曲率：表示构造面梯度变化的快慢 最大主应变：表示形变的大小 张应变（+值）：与裂缝密度有关 压应变（-值）：表示底层压实变形 最大主应力：压应力（+值）：平行于裂缝方向 张应力（-值）：垂直于裂缝方向三、参数和使用说明模块参数面板和功能描述如下：图1应力场数值模拟功能：利用已知的构造数据、纵波速度数据、厚度数据、密度数据等来创建应力场模型 (图1)。

使用提示：以上提到的数据均为网格数据，是在FRanaly 模块中产生的。

地震层位的网格数据是必须的，其它数据可以使用常量来代替。

1) Horizon surface 选择地震层位的网格数据。

2) Velocity of P-wave 选择与地震层位对应的网格化的纵波速度数据3) Thickness 加载层厚度数据，厚度数据是指两层之间的厚度或指定时间窗口范围内的厚度。

每个网格中的厚度值在空间上是变化的。

可以用厚度的平均值来代替厚度网格，相当于常量网格。

4) Constant Thickness 使用常量的地层厚度。

网格中的厚度值相同，空间上无变化。

5) Vp/Vs 纵波横波速度比数据。

可以使用常量代替。

6) Constant Vp/Vs 使用常量的纵横波速度比。

网格中的比值相同，空间上无变化。

7) Rock Density 岩石密度数据。

可以使用常量代替。

8) Constant Density 使用常量的地层密度。

网格中的密度值相同，空间上无变化。

9) Moving Window Length 滑动时窗的长度。

只有时窗内的采样点才用于计算。

计算的结果作为时窗内中间点的输出值。

时窗长度越大，分辨率越低，计算时间越长。

时窗越短，分辨率越高。

建议时窗长度为5－10个采样点。

10) Fault 加载所选层位上的断层数据。

断层属性影响计算结果。

11) Boundary 加载所选层位的边界数据。

12) Principal Stress 输出主应力网格数据文件名，网格指示最小主压应力或最大主张应力。

13) Principal strain (x 100.0) 输出主应变网格数据文件名，指示最大主张应变（正值）或最小主压应力。

14)Principal curvature (1.0/km) 输出曲率网格数据文件名，指示曲率的分布。

15) Principal stress orientation 输出主应力方向网格文件名，指示最大主应力的方向。

零度表示X轴正方向。

应力场方向和强度统计功能：统计最大主应力的方向。

根据当前的矢量图，统计指定的范围内的背景变量输出标量的应力方向的玫瑰图（图2）。

图 2使用提示：1)Angle Interval (degree)指定统计的最大主应力方向的角度间隔（单位：度）。