实习指导书：石油构造分析有限元模拟

实习六　构造应力场有限元数值模拟

一、实习目的和意义

构造应力场的数学模拟是伴随20世纪60年代后期高速、大内存电子计算机的出现而出现的，特别是把计算数学中的有限单元法用于构造应力场分析后，数学模拟得到迅速发展和广泛应用。到今天，有限单元法已成为构造应力场定量研究的一个重要手段。

数学模拟主要是采用数学力学方法，对构造模型的应力场、位移场、应变场、位移速度场和应变速率场等进行定量分析。其主要目的是为了揭示构造应力场分布趋势及展布规律，为构造成因和演化的研究提供定量依据。已有的大量研究表明：构造应力场的有限元数值模拟方法可以解决众多的石油地质问题，主要集中在3个方面（唐永等，2012）：①构造应力场与油气运移之间的关系，通过应力场来分析流体势和油气运移的方向；②构造应力场与岩石破裂准则相结合预测断裂的发育；③构造应力场与流体渗流场耦合分析（流固耦合），讨论构造应力场对流体运移的贡献。

二、实习内容

（1）了解构造应力场有限元数值模拟的基本原理、方法和流程。

（2）学习并掌握有限元数值模拟软件的基本操作和处理过程，能进行简单的构造平面和剖面的有限元数值模拟分析。

（3）根据所给的地质模型以及地质背景资料，构建几何模型、物理模型，进行网格单元的划分、给定边界条件和约束条件，并进行构造应力场的数值模拟分析。

（4）对数值模拟结果进行分析，完成课程实习报告。

三、实习所用资料

（1）研究区构造平面和剖面形态、断层和褶皱的分布、地层格架等。

（2）岩石力学参数、受力作用方式及边界条件。

四、实习所用软件

Dassault Systems Simulia Corp．公司的商业软件Abaqus V6．7版，无任何解题规模及运算限制；Abaqus软件学生教学版，有节点限制。

五、实习步骤

（1）根据地质背景资料和构造形态，构建几何模型。

构造模拟目标体的建立是构造模拟难易程度（计算时间长短、边界条件定义、接触之间的准确约束），以及后期是否建模成功的关键。应以地质背景、构造形态和地质观测分析的数据为主要依据，对目标体进行抽稀。抽稀实际上就是研究层次的聚焦、研究范围的甄选。通常的情况下在研究目标区域构造应力场的分布状况时，我们只需要抓住大区域的主体断层和褶皱即可满足研究的需求，其他小断层和褶皱以及对区域构造应力场分析影响不大的小构造可以在建模时忽略。但在研究目标区域中某条主体断层发育及其对周围地层的影响程度时，必须把所检测到的断层信息以及围岩中发育的断层派生构造尽可能在构造模型中体现出来，即三维模型的建立过程中，在考虑抽稀聚焦的情况下，要尽可能保证后期计算分析时的信息量，这样计算所得到的结果才能合理化、准确化。

地质体演化的真实状况是难以恢复的。应力场分析计算所得到的结果仅仅是与构造变形特征较相符的情况之一，所以我们计算所得到的结果是一种认识上最优、解释较为合理的地质过程，是对该现象的一种有效阐述。这种阐述的合理性、可信性在很大程度上依赖前期抽稀所获得的目标体。一般来说，在目标体的建立上遵循以下几个标准：①研究实体离散块体尽可能少，否则会增加Assembly的难度和工作量；②目标体网格划分及插值方法选择合理，在积分计算所耗时间少；③边界条件易于添加；④较少的接触设置；⑤在后面的计算中能将重点关注区域的应力精确地表达出来。

几何模型的构建主要包括3个方面：①目标选择，确定模型的大小和模拟形态；②对目的层段进行粗化和归并：由于目的层段原始的沉积环境及在成岩过程中所经历的物理化学作用较为复杂，致使在垂向上发育的岩石类型变化较为频繁，为了提高计算速度和模拟结果的精度，可按照岩性特征将整个目的层段进行粗化合并以简化；③优选主要的褶皱和断裂。

（2）物理模型的建立。

几何模型仅仅展示了模拟区构造和地层格架的空间形态，真正能够参与计算的是物理结构模型。建立物理模型一般包括4个方面：①材料属性的确定，即岩石力学参数，包括杨氏模量、泊松比、赫尔抗压强度等；②网格单元划分，为了突显重点区的应力特征以及减少人为误差，采用人工控制与自动相结合的方法，兼顾计算机容量和复杂程度，从地质模型中抽象出力学模型，进行了格化；③边界条件的确定，根据构造作用方式及演化，给出边界约束条件；④确定作用力的大小，从宏观和微观两个方面对构造应力大小进行解析，宏观法利用断层或节理共轭角通过数学解析来获取古应力值（陈庆宣，1996），微观法利用断裂中方解石脉的机械双晶求取应力值（万天丰，1988）。

（3）利用有限元模拟软件Abaqus对构造应力场进行分析计算，处理相关图件，并对应力场模拟结果进行分析。

六、实习指导

1．构造应力场数值模拟方法及原理

一定范围内的岩石在地壳运动的作用下，其岩石内部会出现不同性质、不同方向的应力以抵抗外力作用，这些在空间上彼此关联、时间上持续性的应力所作用的一定范围称为构造应力场。构造应力场中的应力分布和变化是连续而有规律的（陈庆宣，1996）。我们研究构造应力场的目的就是为了揭示研究范围内应力展布规律及其对区域中各种构造发育的控制效应，进而更深入叙述区域构造应力场的性质、大小，然后结合岩石破裂准则、应变能大小来推测可能出现的构造和裂缝，分析构造的成因和演化。

目前结构力学上的有限元方法已经是计算形变和应力的很有效的数值方法。地质构造应力场的数值模拟主要就是应用二维和三维的有限元方法，量化分析构造应力场的大小及展布。但是由于研究的目标体是一个十分复杂的地下岩石块体，它的复杂性和数值表达的困难主要体现在两方面：第一，地壳中种类繁杂的地质构造形态、构造类型、多因素的构造成因是在漫长历史时期地质演化过程中形成的，这种复杂的地质演化过程基本不可能用单一的方法和手段恢复。所以我们只能用相对静止的观点和方法去简化处理岩石体的问题，因此，模拟的目标既是实实在在的又难以用确切的数值表现出来，即计算过程中的条件和参数很难获得准确的数据。第二，研究目标覆盖范围很大，一般都有几十平方千米，并且相对埋深较大，况且影响控制岩石物理性质和地质构造特征的因素是多种多样的，不同地区有不同的特点，即便是同一地区，甚至是同一层位其特征也可能有差异。因此，地质体呈复杂的非均质，而对过于复杂的非均质体，我们在模拟计算中由于较多因素的限制，不可能将目标体考虑的非常详细，只能用相对均匀的岩体区近似地表示实际的地质体。基于以上的种种局限，数值模拟中的模型误差就难以避免，但只要我们将这些误差局限在一定范围内，仍能很好地表现出我们所要描述的地质体的各种属性值。

有限元方法是将弹性连续体离散成为有限个单元，用这些有限个单元近似的表示目标层，其目的是：有限个单元更容易计算，更便于计算机进行数值计算，有利于较好地量化各种属性。针对每个单元通过固体力学变分原理推导各个单元刚度矩阵，将节点位移同节点处所施加的应力的对应关系建立起来，再将各个单元的刚度矩阵集合起来，得到各节点位移的线性方程组，通过计算求解便得到岩石储集体各个单元节点的变形位移情况。有限元法作为一种数值方法，简单来说就是利用插值方法建立复杂问题的近似公式。应用样条插值函数（也称为形函数），将目标体内部各种几何特征的关系转化为离散的数值问题，每个单元上的节点位移值是通过位移形函数与该单元节点的变形位移来计算的。

当然每种数值模拟都是建立在一定的假设的基础上的，有限元实际运用计算也不例外。为了突出问题的本质，并使问题得以简单化和抽象化，一般只考虑弹性变形，在弹性力学中，提出4个基本假设，线性有限元模拟就是在这些假设的基础上而得以顺利实施的，并广泛地运用到岩体力学、土木工程以及结构工程分析上面。①岩层内的物质连续性（Continuity），认为物质中没有空隙，因此可以采用连续函数来描述对象，即在前面所讲到的重叠与开裂的问题。因此在研究过程中，将断层及其附近的区域看成断裂带，与基质形成一个连续的整体，方便后面的网格化及施加边界条件。②岩层内部的物质均匀性（Homogeneity），这里所阐述的均质性，仅仅是一个相对的概念。由于沉积、成岩、构造、时间、生物等的作用，研究区域的储集岩层不可能是均一的，有些研究区的岩石属性非均质相当强烈。但这种非均质性是有一定尺度的，只要研究的标度选择适当，一定范围内的岩层属性变化可以忽略，对整个研究影响不是很大。为此可以将该范围内的岩体识为均质体，可以用一种材料属性来表示区域内的岩石特征。③受力变形主要是线弹性变形（Linear elasticity），岩体变形与外力作用的关系呈线性，外力除去后，岩石储集体可以恢复原状，因此在数值模拟过程中描述岩石力学性质的方程式也是线性的。④小变形（Small deformation）精度的需要，为了使应力场模拟获得较好的结果，在研究过程中会将整个目标体划分成一定规格的细网格。较细的网格决定了研究目标不能发生较大的变形，否则将与前面的连续性假设相矛盾。细网格说明相邻之间的节点距离较短，若发生大的变形，就极易发生节点之间的重叠，这和晶格位错是相同的道理。小变形问题，在数学公式上表现为高阶小量（一般就是在二阶以上的）在计算时可将其忽略。这样可以大大减少计算时间，节约计算资源。

基于以上的假设，有限元的计算主要集中在以下3个方程：静力学平衡方程、几何学的几何方程、材料力学的本构方程。这3个方程将岩体的位移、应变、应力三者较好地串联在一起。其中计算节点位移的方程是由总体刚度矩阵所表示的线性代数方程组，然后由位移场和几何方程计算获得形变场，最后由形变场和本构方程计算应力场。

1）平衡方程（外力、应力）

在外力作用下处于平衡状态时，其内部各单元节点处的应力状态是各不相同的。在研究区域内任一点均可找到一个微六面体，将其应力分解为σx、σy、σz、τxy、τyz、τzx。同样外力F也可以根据坐标分解为F x、F y、F z［式（6－1）］。根据理论力学，在外力作用下平衡状态的物体，其表面各点处的应力分量应当与作用在该处的外力分量相平衡，平衡方程较好地将边界静力条件与施加在目标体单元上的应力联系建立起来。

式中：σx、σy、σz——x、y、z方向上的主应力，Pa；

　　　τxy、τyz、τzx——平面xy、yz、zx上的剪切应力，Pa；(www.chuimin.cn)

　　　F x、F y、F z——作用于某一点沿x、y、z的外力分量，N。

　　　2）几何方程（位移、应变）

当目标岩层受到外界作用力的影响，构造岩层内部各单元点的位置会发现变化，这种变化包括两种形式：①岩层内各点的绝对位置虽然均有变化，但任意两点之间的相对距离却始终未变，目标体仅仅是整体位置发生了变动，而岩层本身的形状、体积都没有改变，即在变化过程中无形变产生，这对构造应力场的反演一般无较大的借鉴价值；②岩体内部任意两点之间的相对距离发生了改变，从而使其形状和尺寸发生了变化，岩层整体产生了变形，这种变形变位对恢复构造应力场价值较大。因此研究岩层在外力作用下的变形规律，只需要研究物体内各点的相对位置变动情况。当然我们所选择的目标体必须是连续体，也就是在受力变形过程中，自始至终保持连续性而无任何重叠和开裂现象产生。

由于刚体位移不产生变形，下面将仅仅讨论变形位移。若球形岩体在应力作用下变形成为椭球体，其内部一点M变形后到达M′位置，则其变形位移MM′可以转化为x、y、z方向的3个分量u、v、w［式（6－2）］。岩体内各点的位移不同，因此位移分量应是点的位置坐标函数：

式中：u、v、w——x、y、z三轴方向的3个位移分量，量纲依据所研究的尺度来确定。

在三维空间内（图6－1），受力岩体的某一点沿三轴x、y、z方向上的线应变εx、εy、εz，以及过M点参照x、y、z轴所取3个相互垂直棱边所夹直角的改变量，也即剪应变γxy、γyz、γzx分别为：

式中：εx、εy、εz——x、y、z方向的线应变，无量纲；

　　　γxy、γyz、γzx——xy、yz、zx面上的剪应变，无量纲。

式（6－2）表明岩层各处的变形位移，但还不能直观地说明目标体变形的强烈程度。式（6－3）则较好地建立了位移分量和应变分量之间的关系。

图6－1　变形位移示意图

（据唐永等，2011）

3）本构方程（应变、应力）

应力与应变是相辅相成的，对某种具体的岩石储集体，在一定围压作用下，应力和应变之间有一定的关系，这种关系是可以量化确定的。根据Ladanyi D（1970）致密灰岩破裂实验，在弹性变形过程中，岩石将能量以变形能的形式保留，一旦外力卸载，岩体会将应变能释放出来，变形也将完全恢复；应力和应变成线性关系变化，并且一一对应。由于在实验室条件下，不能完全模拟岩层所经历的地质条件的变化，所以岩层的塑性特征很难获取，一般在研究岩层应变能时只考虑弹性形变，本次研究过程中将目标岩层在一定范围内视为各向同性体。在此前提下，应力应变关系为：

式中：E——杨氏弹性模量，Pa；

　　　μ——泊松比，无量纲。

2．构造应力场数值模拟的基本流程

构造应力场建模可以归结为3个方面（图6－2）：①目标体的建立（前处理来完成）；②数值计算（求解器中数值分析来完成）；③模拟结果输出及分析（后处理来完成）。其主体思路：建立目标岩体的几何模型将目标岩体介质划分成不同的有限个单元（二维问题可以划分三角形、矩形或任意多边形等，三维模型可以划分多面体），并对单元进行分析，然后将每个单元结合起来进行整体研究，最后将计算结果可视化（图6－3）。

图6－2　Abaqus软件有限元数值模拟处理过程

（据沈传波等，2004）

图6－3　构造应力场数值模拟流程图

（据唐永等，2011）

3．构造应力场数值模拟成果图示例（图6－4～图6－9）

图6－4　宣汉－达县地区（a）有限元分析物理模型和（b）施加的边界条件

（据唐永等，2012）

图6－5　宣汉－达县地区飞仙关组燕山晚期最大主应力分布图

（据唐永等，2012）

图6－6　广西灵山地区最大主应力分布图

（据唐永等，2014）

图6－7　广西灵山地区最大主应力迹线图

（据唐永等，2014）

图6－8　某盆地地震测线有限元数值分析网格化模型

图6－9　某盆地地震测线有限元数值模拟最大主应力分布