CN113820750B - 基于弹塑性力学定量预测泥岩构造裂缝的方法 - Google Patents

Abstract

本申请提供一种基于弹塑性力学定量预测泥岩构造裂缝的方法，涉及石油勘探开发技术领域。其中，该定量预测泥岩构造裂缝的方法包括：单井厚层泥岩段构造裂缝的分析和统计；进行岩石力学实验及力学参数测井解释；计算脆性指数及恢复古力学参数；泥岩屈服‑破坏复合强度准则建立；泥岩构造裂缝力学参数模型推导；地质模型构建及古应力状态恢复；三维古应力场模拟及裂缝定量预测。本申请技术方案解决了脆性指数约束下的泥岩屈服‑破裂复合强度准则建立的难题，适合于任何含泥岩隔层型碎屑岩储层的裂缝定量预测工作。

Description

基于弹塑性力学定量预测泥岩构造裂缝的方法

技术领域

本申请涉及石油勘探开发技术领域，尤其涉及一种基于弹塑性力学定量预测泥岩构造裂缝的方法。

背景技术

随着油气勘探开发的不断进行与开采技术的持续进步，复杂断块油气藏、岩性油气藏、非常规致密砂岩油气藏以及裂缝性油气藏受到人们越来越多的关注；由于我国含油气盆地以陆相沉积地层为主，频繁的湖平面变化，造成了垂向上砂岩、泥岩的频繁交互变化，在多期构造运动作用下，不同厚度的泥岩段内也发育了一系列不同力学性质的裂缝，统计发现主要以低角度剪切缝为主，且裂缝内极少充填，两壁常出现光滑的摩擦镜面。

目前，针对裂缝的大量研究主要集中储层段内，逐渐形成了较为成熟的地质、地球物理和力学识别和预测方法，但对泥岩和页岩裂缝的研究仍停留在单井描述和平面预测阶段，相对于砂岩储层，泥岩由于掺杂了灰质、硅质、石膏等矿物成分，在受力变形过程中表现出不同程度的复杂塑性力学行为，给井距大、取芯少、泥岩厚地区的裂缝预测工作带来了很大难题。

发明内容

本申请实施例的目的是提供一种基于弹塑性力学定量预测泥岩构造裂缝的方法，以实现对泥岩和页岩裂缝的预测研究。

为解决上述技术问题，本申请实施例提供如下技术方案：

本申请提供一种基于弹塑性力学定量预测泥岩构造裂缝的方法，该基于弹塑性力学定量预测泥岩构造裂缝的方法包括：单井厚层泥岩段构造裂缝的分析和统计；进行岩石力学实验及力学参数测井解释；计算脆性指数及恢复古力学参数；泥岩屈服-破坏复合强度准则建立；泥岩构造裂缝力学参数模型推导；地质模型构建及古应力状态恢复；三维古应力场模拟及裂缝定量预测。

在本申请的一些变更实施方式中，所述单井厚层泥岩段构造裂缝的分析和统计，具体为：

收集现场成像测井资料和钻井取芯资料，统计泥岩段地层产状和构造裂缝参数；

绘制单井裂缝发育柱状图，明确厚层泥岩段裂缝分布的差异性，针对不同裂缝发育程度的深度点采样进行镜下薄片观察，识别矿物类型及含量，分析其对裂缝密度的影响；

其中，所述地层产状包括层理倾角和裂缝产状，所述裂缝参数包括裂缝线密度参数、裂缝体密度参数、裂缝开度参数及裂缝长度参数。

在本申请的一些变更实施方式中，所述进行岩石力学实验，具体为：

在泥岩厚层段，钻取柱塞岩样，且定义岩心钻取方向与轴向垂直，进行巴西劈裂试验、三轴压缩试验和直剪试验，测试并计算泥岩的Ⅰ型和Ⅱ型断裂韧性，同时记录各类试验的应力-应变曲线；

其中，所述三轴压缩试验用于获取岩石力学参数，所述岩石力学参数包括密度、弹性模量、泊松比、抗压强度、屈服强度、内摩擦系数及内摩擦角参数，所述巴西劈裂试验用于获得岩石的抗压强度，所述直剪试验用于获取岩石的抗剪强度、内摩擦系数以及内摩擦角；

在具有层理面的泥岩段，钻取柱塞岩样，且定义岩心钻取方向与轴向平行，进行巴西劈裂试验和直剪试验，同时记录各类试验的应力-应变曲线；

其中，所述巴西劈裂试验用于获得所述层理面的抗拉强度，所述直剪试验用于获得所述层理面的抗剪强度、内摩擦系数和内摩擦角。

在本申请的一些变更实施方式中，所述力学参数测井解释，具体为：

通过全波列测井和偶极横波测井计算得到横波时差曲线，并与常规测井资料进行相关性分析，拟合泥岩段的横波时差预测公式(1)，并通过公式(2)、(3)、(4)、(5)、(6)计算单井动态岩石力学参数；

Δt_s＝0.948Δt_p+57.42 (1)

式中：Δt_s为横波时差(μs/m)，Δt_p为纵波时差(μs/m)；

S_c＝E_d0.008V_sh+0.0045(1-V_sh) (4)

式中：E_d为动态杨氏模量(GPa)，ρ_b为体积密度(g/cm³)，μ_d为动态泊松比，S_c为单轴抗压强度(MPa)，C为岩石内聚力(MPa)，为层理面的内摩擦角(°)，V_sh为泥质的体积含量；

基于所述单井动态岩石力学参数结果进行断裂韧性计算，得出针对泥岩的Ⅰ型和Ⅱ型断裂韧性计算模型(7)(8)；

K_ⅡC＝0.0466P+0.1674S_t-0.1851 (8)

式中：P为围压(MPa)，S_t为互层岩石的抗拉强度(MPa)。

在本申请的一些变更实施方式中，所述计算脆性指数及恢复古力学参数，具体为：

通过公式(9)(10)对所述单井动态岩石力学模型进行连续矫正，制作单井静态力学参数图；

E_S＝2.902E_d-14.92 (9)

μ_S＝4.59μ_d-1.07 (10)

式中：E_S为静态杨氏模量(GPa)，μ_S为静态泊松比；

结合密度测井资料、阵列声波或偶极声波所测量的纵、横波数据和杨氏模量、泊松比参数，优选泊-杨法根据公式(11)计算泥岩段的脆性指数；

式中：BI_e为弹性参数脆性指数，E_max和E_min分别为某一层段杨氏模量的最大和最小值，μ_max和μ_min分别为某一层段泊松比的最大和最小值，E为杨氏模量GPa；

收集地震叠前反演得到的研究区的P波阻抗、S波阻抗和密度数据体，并以单井测井解释结果为约束，计算现今应力条件下的弹性模量、泊松比，并运算得到泥岩段的脆性指数体；

收集所述研究区的深度域地震解释剖面，所述剖面垂直于主体构造走向，且所述剖面包含目的层和断裂解释结果，采用构造平衡剖面法恢复断距和二维褶皱形态，确定构造强烈活动期或关键造缝期；

建立所述研究区的泥岩层段构造地质模型，在划分断层相对级别的前提下，以一级、二级断层为控制，采用断层相关褶皱理论进行关键造缝期的古构造形态恢复；

通过三轴力学实验拟合建立不同围压条件下的力学参数变化曲线，将围压等效为埋藏深度下的上覆地层压力，根据所述古构造形态恢复的结果得到关键造缝期的泥岩层段古埋深，结合力学参数变化曲线计算该时期的古力学参数和脆性指数三维体，并按照脆性指数的分布趋势和区间特征将其进行等级划分。

在本申请的一些变更实施方式中，所述泥岩屈服-破坏复合强度准则建立，具体为：

判断三向应力状态下层理面的滑移和破坏情况；

判断三向压应力状态下层理面的滑移和破坏情况；

判断基岩的屈服和破裂状态，在三向挤压应力状态下，采用弹塑性屈服准则和塑形硬化准则判断岩石的屈服临界状态和破裂状态；

判断基岩中的裂缝倾角。

在本申请的一些变更实施方式中，所述判断三向应力状态下层理面的滑移和破坏情况，具体为：

当最小主应力σ₃<0或σ₃为张应力时且与地层的夹角为45°-90°时，采用公式(12)判断层理面是否发生滑移和破坏情况；

式中：F^t为岩石发生张破裂的应力条件，F^S为岩石发生剪破裂的应力条件，σ_t为层理面的抗张强度(MPa)，μ为岩石的泊松比，P为围压(MPa)；

所述判断三向压应力状态下层理面的滑移和破坏情况，具体为：

当水平挤压力σ₁与地层的夹角为30°-60°时，采用公式(13)判断层理面是否发生滑移和破坏；

式中：c_j为层理面的内聚力(MPa)，为层理面的内摩擦角(°)，θ_h为水平挤压力σ₁与地层的夹角(°)；

所述采用弹塑性屈服准则和塑形硬化准则判断岩石的屈服临界状态和破裂状态，具体为：

采用公式(14)判断岩石的屈服临界状态和破裂状态；

式中：I₁为应力第一不变量，J₂为应力第二不变量，a和K为与岩石内摩擦角和内聚力相关的实验常数，为岩石的内摩擦角(°)。

在本申请的一些变更实施方式中，所述泥岩构造裂缝力学参数模型推导，具体为：

选取岩石内表征单元体进行分析并简化，得到公式(15)；

式中：w_总为当前应力状态下所述单元体受力吸收的总应变能密度(J/m³)，w_d为岩石屈服后发生塑性变形所消耗的塑性应变能密度(J/m³)，w为存储的可释放弹性应变能密度(J/m³)，w_e为想产生裂缝必须克服的弹性应变能密度(J/m³)，V表征所述单元体的体积(m³)，S_f为新增裂缝的表面积(m²)，J为产生单位面积裂缝所需要的能量(J/m²)；

根据能量守恒原理和最大应变能密度理论，得到公式(16)；

式中：D_vf为所述单元体内缝隙的体密度(m²/m³)，E为性模量(GPa)，μ为泊松比；

推导建立三向挤压应力状态下泥岩和层理面裂缝密度、开度的力学模型，采用拉-剪复合破裂、莫尔－库仑破裂准则和德鲁克-普拉格准则进行岩石破裂判别，若满足破裂条件，通过公式(17)进行裂缝参数计算；

式中：J₀为零围压时裂缝表面能(J/m²)，ε₀为最大弹性张应变，σ₃为最小张应变，K₀为比例系数，w为应变能密度(J/m³)，w_f为裂缝应变能密度(J/m³)，w_e为新增单位裂缝表面积必须克服的应变能(J/m³)，σ_p为岩石破裂应力(MPa)，D_lf为所述单元体内缝隙的线密度(l/m)，b为裂缝开度m，E为杨氏模量GPa，为岩石的内摩擦角(°)，θ为破裂角(°)，C₀为岩石的内聚力MPa，L₁、L₃为最大主应力、最小主应力方向的单元网络长度m，b为裂缝开度m；

其中，所述最大弹性张应变为对应于0.85σ_c时的应变，所述最小张应变为对应0.45σ_c时的张应变。

在本申请的一些变更实施方式中，所述地质模型构建及古应力状态恢复，具体为：

收集地震解释、地质分层和测井解释结果，构建所述研究区的断层模型以及储层地质模型，以进行面模型重构，并通过面模型重构进行三维模型构建，完成三维地质力学模型的建立；

对所述研究区的钻井岩心进行声发射试验取样，每组样品至少包括4个试样；

将加工好的所述试样进行重复加载实验，加载速率为0.1MPa/s，测出所述试样在受载过程中的声发射信号随载荷的变化曲线，在第二次加载声发射载荷变化曲线上，找出声发射抹录不尽点，参照所述抹录不尽点的载荷值，在第一次加载声发射载荷曲线中，确定出凯塞尔点，取所述凯塞尔点和抹录不尽点载荷的平均值，即为岩心在地下所受的最大正应力，所述凯塞尔点的数量代表构造运动的期次；

根据构造形迹力学方法，进行裂缝的分期配套，进行裂缝走向统计分析，找出共轭裂缝的优势组系，以其共轭夹角平分线作为关键造缝期古应力场最大主应力方向；

其中，在所述声发射试验取样中，当试样的数量为4个时，其中一个所述试样取自品行井眼轴线的方向，其余所述试样取自水平面内互成45°角的三个方向。

在本申请的一些变更实施方式中，所述三维古应力场模拟及裂缝定量预测，具体为：

基于所述三维地质力学模型和关键造缝期古应力测试结果，设置模型力学边界，设定加载速率为0.1MPa/s，进行弹塑性有限元三维应力场模拟，并反复尝试，直到古应力场的最大主应力分量与次级挤压断层分布趋势、最小主应力分量与次级伸展断层分布一致、剪应力分量与走滑断层分布趋势一致为止；

基于得到的破裂准则、屈服准则、裂缝密度、开度和产状力学模型通过软件获取每个节点的最大主应力、最下主应力、中间主应力、剪应力、最大主应变、最下主应变、中间主应变、屈服强度参数，并经过计算得到每个节点裂缝的体密度、线密度和开度值；

以基岩力学层为单位进行单井预测裂缝线密度、开度与实测值对比，当平均吻合率低于90％时，修正裂缝力学参数、调整边界条件重新模拟三维应力场，直至平均吻合率达到90％以上为止，完成厚层泥岩段构造裂缝的三维定量预测。

相较于现有技术，本申请提供的基于弹塑性力学定量预测泥岩构造裂缝的方法，该方法解决了脆性指数约束下的泥岩屈服-破裂复合强度准则建立的难题，适合于任何含泥岩隔层型碎屑岩储层的裂缝定量预测工作；将任何软、硬岩石中的裂缝产生都归结为可释放弹性应变能的释放，推导建立了泥岩裂缝参数的力学模型，依托有限元模拟平台，可有效预测碎屑岩储层中构造裂缝的空间分布特征，包括裂缝密度、开度和产状，为泥岩裂缝性油气藏的勘探开发提供了地质依据，为避免钻井过程中井漏和井壁失稳事故的发生提供了保障，减少了勘探开发的风险和成本。

附图说明

通过参考附图阅读下文的详细描述，本申请示例性实施方式的上述以及其他目的、特征和优点将变得易于理解。在附图中，以示例性而非限制性的方式示出了本申请的若干实施方式，相同或对应的标号表示相同或对应的部分，其中：

图1示意性地示出了本发明实施例提供的种基于弹塑性力学定量预测泥岩构造裂缝的方法的流程示意图；

图2示意性地示出了本发明实施例提供的种基于弹塑性力学定量预测泥岩构造裂缝的方法的泥岩p-q塑性屈服曲线中参数对初始屈服面的影响示意图；

图3示意性地示出了本发明实施例提供的种基于弹塑性力学定量预测泥岩构造裂缝的方法的泥岩内变量等效塑性剪应变和硬化函数的关系示意图；

图4示意性地示出了本发明实施例提供的种基于弹塑性力学定量预测泥岩构造裂缝的方法的单位体积泥岩样品中能量的分配关系示意图。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施方式。虽然附图中显示了本公开的示例性实施方式，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

需要注意的是，除非另有说明，本申请使用的技术术语或者科学术语应当为本申请所属领域技术人员所理解的通常意义。

实施例一

参考附图1-附图4所示，本发明的实施例一提出一种基于弹塑性力学定量预测泥岩构造裂缝的方法，该基于弹塑性力学定量预测泥岩构造裂缝的方法包括：

步骤1：单井厚层泥岩段构造裂缝的分析和统计；

步骤2：进行岩石力学实验及力学参数测井解释；

步骤3：计算脆性指数及恢复古力学参数；

步骤4：泥岩屈服-破坏复合强度准则建立；

步骤5：泥岩构造裂缝力学参数模型推导；

步骤6：地质模型构建及古应力状态恢复；

步骤7：三维古应力场模拟及裂缝定量预测。

具体的，本方案针对塑性变形复杂的厚层泥岩段，采用岩心观察和测井解释进行单井构造裂缝的统计分析，采用岩石力学实验和弹塑性力学理论推导建立泥岩屈服破裂准则以及裂缝参数表征模型，采用地震解释和地层对比构建三维地质模型，在关键造缝期应力状态恢复的前提下，定量预测三维古应力场分布及裂缝空间分布。

进一步的，步骤1：单井厚层泥岩段构造裂缝的分析和统计，具体包括：

步骤1.1：由于成岩缝多数被矿物充填，收集现场成像测井(FMI)资料和钻井取芯资料，统计泥岩段地层产状和构造裂缝参数，包括层理倾角和裂缝产状、线密度、体密度、开度、长度参数；

步骤1.2：基于观察、统计结果，绘制单井裂缝发育柱状图，明确厚层泥岩段裂缝分布的差异性，针对不同裂缝发育程度的深度点采样进行镜下薄片观察，识别矿物类型及含量，分析其对裂缝密度的影响。

进一步的，步骤2：进行岩石力学实验及力学参数测井解释，具体包括：

步骤2.1：针对泥岩厚层段，钻取柱塞样，定义岩芯钻取方向与轴向垂直，参照岩石力学试验规范，岩样加工成尺寸设定为5cm×2.5cm×2.5cm的圆柱形岩样，并开展巴西劈裂试验、单轴/三轴压缩实验及直剪试验；其中，三轴压缩实验所用的岩样采用相近深度、相近岩性的三块为一组，并设置不同的围压值(0MPa、5MPa、10MPa……60MPa)，获取岩石力学参数，包括密度、弹性模量、泊松比、抗压强度、屈服强度、内摩擦系数和内摩擦角；通过巴西劈裂试验测试获得岩石的抗拉强度，同时在泥岩内设置初始裂缝，长度为2a，裂缝与加载方向夹角为θ，测试并计算出泥岩的I型和II型断裂韧性，记录各类岩石力学实验的应力-应变曲线；

步骤2.2：针对具有层理面的泥岩段，定义岩芯钻取方向与岩芯轴向平行，开展巴西劈裂试验和直剪试验；其中，通过巴西劈裂试验测试获得层理面的抗拉强度，通过变围压下的直剪试验获得层理面的抗剪强度、内摩擦系数和内摩擦角，同时记录各类岩石力学实验应力-应变曲线；

步骤2.3：针对常规测井资料缺少横波时差曲线的特点，用全波列测井和偶极横波测井计算得到横波时差曲线，并与常规测井资料进行相关性分析，拟合泥岩段的横波时差预测公式，即公式(1)：

Δt_s＝0.948Δt_p+57.42 (1)

式中：Δt_s为横波时差(μs/m)，Δt_p为纵波时差(μs/m)；

进而，可计算单井动态岩石力学参数，具体为：通过公式(2)计算动态弹性模量、通过公式(3)计算动态泊松比、通过公式(4)计算单轴抗压强度、通过公式(5)计算岩石内聚力以及通过公式(6)计算层理面的内摩擦角，公式(2)-(6)分别为：

S_c＝E_d0.008V_sh+0.0045(1-V_sh) (4)

式中：E_d为动态杨氏模量(GPa)，ρ_b为体积密度(g/cm3)，μ_d为动态泊松比，S_c为单轴抗压强度(MPa)，C为岩石内聚力(MPa)，为层理面的内摩擦角(°)，V_sh为泥质的体积含量(小数)；

步骤2.4：基于测井计算力学参数结果进行断裂韧性计算，优选出针对泥岩的I型和II型断裂韧性计算模型，分别为下述公式(7)和(8)；

K_ⅡC＝0.0466P+0.1674S_t-0.1851 (8)

式中：P为围压(MPa)，S_t为互层岩石的抗拉强度(MPa)。

进一步的，步骤3：计算脆性指数及恢复古力学参数，具体包括：

步骤3.1：由测井曲线计算出的动态力学参数和实验得到的静态参数会存在一定误差，静态力学参数更符合实际工程需要，因此，需要优选拟合公式分析动静转化关系，泥岩静态弹性模量采用公式(9)，泥岩静态泊松比采用公式(10)，用这两个公式对测井计算的力学参数进行连续校正，并制作单井静态力学参数图，公式(9)和(10)分别为：

E_S＝2.902E_d-14.92 (9)

μ_S＝4.59μ_d-1.07 (10)

式中：E_S为静态杨氏模量(GPa)，μ_S为静态泊松比；

步骤3.2：结合密度测井资料、阵列声波或偶极声波所测量的纵、横波数据和杨氏模量、泊松比参数，优选泊-杨法根据公式(11)计算泥岩段的脆性指数，公式(11)为：

式中：BI_e为弹性参数脆性指数，E_max和E_min分别为某一层段杨氏模量的最大和最小值，μ_max和μ_min分别为某一层段泊松比的最大和最小值，E为杨氏模量GPa。

步骤3.3：收集地震叠前反演得到的研究区P波阻抗、S波阻抗和密度数据体，在此基础上以单井测井解释结果为约束，计算现今应力条件下的弹性模量、泊松比，并运算得到泥岩段的脆性指数体；

步骤3.4：收集研究区深度域地震解释剖面，剖面尽量垂直于主体构造走向，剖面包含目的层和断裂解释结果，采用构造平衡剖面法恢复断距和二维褶皱形态，确定构造强烈活动期或关键造缝期；

步骤3.5：建立研究区泥岩层段构造地质模型，在划分断层相对级别的前提下，以一级、二级断层为控制，采用断层相关褶皱理论进行关键造缝期古构造形态恢复，具体地，采用弯曲去褶皱法、恢复到基准面法恢复伸展构造区古构造，采用“拼板恢复”法恢复走滑构造区古构造，采用斜剪切法、弯曲滑动法和断层平行流法恢复挤压区古构造；

步骤3.6：通过三轴力学实验拟合建立不同围压条件下的力学参数变化曲线，将围压等效为埋藏深度下的上覆地层压力，根据古构造恢复结果可以得到关键造缝期的泥岩层段古埋深，结合力学参数变化曲线计算出此时期的古力学参数和脆性指数三维体，进一步按照脆性指数的分布趋势和区间特征，将其划分为三个等级：强(脆性)、中等(脆-塑性)和弱(塑性)。

进一步的，步骤4：泥岩屈服-破坏复合强度准则建立，具体包括：

根据大量岩心观察结果，当泥岩段富含碳酸盐岩、硅质矿物时，即脆性指数较大时，常显示出明显的层理构造，将泥岩段的变形和破坏方式分为层理面滑移型和基质破坏型两种基本类型；

步骤4.1：首先判断三向应力状态下层理面的滑移、破坏情况，当最小主应力σ₃<0或张应力时与地层的夹角为45°-90°时，泥岩位于层面拉张破裂或剪切破裂可能范围内，但油气田地下的实际地层倾角极少能超过75°，根据平行四边形法则，泥岩层理面上既有张应力也有剪应力，更倾向于产生拉-剪复合破裂，采用公式(12)判断层理面是否发生滑移和破坏情况；

式中：F^t为岩石发生张破裂的应力条件，F^S为岩石发生剪破裂的应力条件，σ_t为层理面的抗张强度(MPa)，μ为岩石的泊松比，P为围压(MPa)，即σ₃；

步骤4.2：判断三向压应力状态下层理面的滑移、破坏情况，当水平挤压力σ₁与地层的夹角为30°-60°时，泥岩段位于层面剪切破裂可能范围内，此时采用莫尔－库仑剪切准则判别层理面是否发生滑移破坏，采用公式(13)判断层理面是否发生滑移和破坏；

式中：c_j为层理面的内聚力(MPa)，为层理面的内摩擦角(°)，θ_h为水平挤压力σ₁与地层的夹角(°)；

步骤4.3：判断基岩，即泥岩的屈服、破裂状态，三向挤压应力状态下，采用弹塑性屈服准则(德鲁克-普拉格准则)和塑性硬化准则判别岩石的屈服临界状态和破裂状态，德鲁克-普拉格屈服准则认为当材料的最大形状变形能(又称畸变能)达到某一常数时，材料就屈服并进入塑性状态，采用公式(14)判断岩石的屈服临界状态和破裂状态：

式中：I₁为应力第一不变量，J₂为应力第二不变量，a和K为与岩石内摩擦角和内聚力相关的实验常数， 为岩石的内摩擦角(°)；

步骤4.4：根据步骤4.3，当泥岩发生屈服后，仍然有可能发生破裂，依据脆性指数等级分两种情况进行判断，当泥岩脆性指数位于弱等级时，符合理想弹塑性本构材料，进入初始无限塑性状态即预示岩石很难发生规模性破坏，主要以发生持续塑性变形为主，而当泥岩脆性指数位于强或中等等级时，符合硬化弹塑性本构材料，破坏面是硬化面和软化面的分界面，或者说破坏是塑性变形过程发展的最终结果，这样，可以定义岩石峰值应力p-q塑性屈服曲线，参考附图2所示，根据公式(18)确定出适合泥岩的屈服-破坏准则；

式中：α_p为硬化函数，即塑性变形破坏过程中屈服面将随着硬化函数演变，C₀和A分别是岩石内聚力和内摩擦系数，通过以上力学实验获得，其中， 为岩石的内摩擦角(°)，这里为保证参数A无量纲，设置P₀＝1MPa，又σ_1p为不同围压σ₃下的泥岩峰值应力与前面砂岩的σ_c相当，σ_ij为加载过程中的应力偏量(MPa)，δ_ij为克罗内克符号；

步骤4.5：对于泥岩屈服-破坏准则中的硬化函数α_p，其公式(19)如下：

式中：为硬化函数初始值，可由初始屈服点予以确定，参考附图2所示，当α_p＝1表示岩石达到了宏观破坏时的硬化状态，故其取值范围为b为硬化参数，表示岩石的塑性硬化率，可基于硬化函数α_p和塑性内变量γ_p的关系来确定，参考附图3所示，b值越小，α_p变化越快，且逐渐趋近于1，为塑性应变增量偏张量，为塑性总应变量，ε^p为塑性体积塑性应变量，γ_p为等效剪切塑性应变量，这些应变和应力参数可在Ansys弹塑性应力场模拟过程中直接提取；

步骤4.6：判断基质(泥岩)中裂缝的倾角，在三维空间中，裂缝的倾角和走向要采取投影计算的方法来确定，Ansys坐标系的X轴与大地坐标的X轴(东)重合，Z轴与Y轴负方向(南)重合，Y和Z轴重合，因此若裂缝面法线方向向量在整体坐标系下的方向余玄已经确定为将投影到XOZ平台，其投影线与Z轴负方向的夹角为α_z＝arctan(-l/n)，走向-角α可确定：若0≤α_z<90°，α＝90°-α_z，若-90°<α_z<0，α＝(-90°-α_z)+360°，在Ansys三维坐标系中，裂缝倾角是裂缝表面与XY平面的夹角，也就是平面lx+my+nz＝0与平面y＝0之间的夹角α_dip(0°≤α_dip≤90°)，其公式(20)如下：

进一步的，步骤5：泥岩构造裂缝力学参数模型推导，具体包括：

步骤5.1：大量岩石力学实验结果表明：可释放弹性应变能是促使岩石破坏的根本原因，泥岩内主应力达到单轴抗压强度σ_c的0.85倍时，大量微裂缝开始集结并产生宏观裂缝，因此可以把该时刻的应变能称作要产生裂缝而必须克服的最大弹性应变能，于是新增裂缝表面积需要的能量就等于当前单元体内应变能与要产生裂缝而必须克服的最大弹性应变能之差，选取岩石内表征单元体进行分析(REV)，在适当简化的前提下，参考附图4所示，得到公式(15)；

式中：w_总为当前应力状态下所述单元体受力吸收的总应变能密度(J/m³)，w_d为岩石屈服后发生塑性变形所消耗的塑性应变能密度(J/m³)，w为存储的可释放弹性应变能密度(J/m³)，w_e为想产生裂缝必须克服的弹性应变能密度(J/m³)，V表征所述单元体的体积(m³)，S_f为新增裂缝的表面积(m²)，J为产生单位面积裂缝所需要的能量(J/m²)；

步骤5.2：由能量守恒原理和最大应变能密度理论，泥岩新增裂缝总表面积等于裂缝体密度，于是有公式(16)；

式中：D_vf为所述单元体内缝隙的体密度(m²/m³)，即裂缝体的总表面积与单元体体积之比，E为杨氏模量GPa，μ为泊松比；

步骤5.3：推导建立三向挤压应力状态下泥岩和层理面裂缝密度、开度的力学模型，采用拉-剪复合破裂、莫尔－库仑破裂准则和德鲁克-普拉格准则进行岩石破裂判别，若满足破裂条件，用公式(17)进行裂缝参数计算；

式中：J₀为零围压时裂缝表面能(J/m²)，ε₀为最大弹性张应变，对应于0.85σ_c时的应变，σ₃为最小张应变，对应0.45σ_c时的张应变，K₀为比例系数，w为应变能密度(J/m³)，w_f为裂缝应变能密度(J/m³)，w_e为新增单位裂缝表面积必须克服的应变能(J/m³)，σ_p为岩石破裂应力(MPa)，D_lf为所述单元体内缝隙的线密度(l/m)，b为裂缝开度m，为内摩擦角(°)，θ为破裂角(°)，其中，破裂压力可由公式(17)或公式(22)换算得到，对应0.85σ_c这一关键点。

进一步的，步骤6：地质模型构建及古应力状态恢复，具体包括：

步骤6.1：以Petrel地质建模软件为平台，收集地震解释、地质分层和测井解释结果，构建研究区断层模型以及储层地质模型，将层面和断层数据导入AutoCad软件进行面模型重构，并将面模型再次导出到Ansys软件中进行三维模型构建，划分网格后按照“形心”搜索法将井-震反演获得的古力学参数离散赋值到三维模型的每个网格内，完成三维地质力学模型的建立；

步骤6.2：对研究区钻井岩心进行声发射试验取样，钻井位置最好平均分布在研究区平面上或呈四边形，每组样品至少包括4个试样，其中1个取自垂直方向(平行井眼轴线)，另3个取自水平面内互成45°角的3个方向；

步骤6.3：声发射试验所用设备由伺服岩石刚性试验机和声发射测试系统组成，将加工好的试样在室内进行重复加载实验，加载速率为0.1MPa/s，测出岩样受载过程中的声发射信号随载荷的变化曲线，在第二次加载声发射载荷变化曲线上，找出声发射抹录不净点，参照抹录不尽点载荷值，在第一次加载声发射载荷曲线中，确定出凯塞尔点，取凯塞尔点和抹录不尽点载荷的平均值，即为岩心在地下所受的最大正应力，凯塞尔点的数量即代表了构造运动的期次；

步骤6.4：根据构造形迹力学方法，进行裂缝的分期配套，采用玫瑰花图进行裂缝走向统计分析，找出共轭裂缝的优势组系，以其共轭夹角平分线作为关键造缝期古应力场最大主应力方向。

进一步的，步骤7：三维古应力场模拟及裂缝定量预测，具体包括：

步骤7.1：基于Ansys中建立的三维地质力学模型和关键造缝期古应力测试结果，设置模型力学边界，设定加载速率为0.1MPa/s，进行弹塑性有限元三维应力场模拟，并反复尝试，直到古应力场的最大主应力分量与次级挤压断层分布趋势、最小主应力分量与次级伸展断层分布一致、剪应力分量与走滑断层分布趋势一致为止：

步骤7.2：将上述破裂准则、屈服准则、裂缝密度、开度和产状力学模型编写成APDL语言程序，植入到Ansys软件平台中，在古应力场模拟的基础上，提取每个节点的最大主应力、最下主应力、中间主应力、剪应力、最大主应变、最下主应变、中间主应变、屈服强度参数，经过计算后可得到每个节点裂缝的体密度、线密度和开度值；

步骤7.3：以岩石力学层为单位进行单井预测裂缝线密度、开度与实测值对比，如果平均吻合率低于90％，则修正裂缝力学参数、调整边界条件重新模拟三维应力场，直到平均吻合率达到90％以上为止，最终完成厚层泥岩段构造裂缝的三维定量预测。

根据上述所列，本发明实施例提出一种基于弹塑性力学定量预测泥岩构造裂缝的方法，能够解决脆性指数约束下的泥岩屈服-破裂复合强度准则建立的难题，适合于任何含泥岩隔层型碎屑岩储层的裂缝定量预测工作；将任何软、硬岩石中的裂缝产生都归结为可释放弹性应变能的释放，推导建立了泥岩裂缝参数的力学模型，依托有限元模拟平台，可有效预测碎屑岩储层中构造裂缝的空间分布特征，包括裂缝密度、开度和产状，为泥岩裂缝性油气藏的勘探开发提供了地质依据，为避免钻井过程中井漏和井壁失稳事故的发生提供了保障，减少了勘探开发的风险和成本。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (8)

1.一种基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，包括：

单井厚层泥岩段构造裂缝的分析和统计；

进行岩石力学实验及力学参数测井解释；

计算脆性指数及恢复古力学参数；

泥岩屈服-破坏复合强度准则建立；

泥岩构造裂缝力学参数模型推导；

地质模型构建及古应力状态恢复；

三维古应力场模拟及裂缝定量预测；

其中，所述泥岩屈服-破坏复合强度准则建立，具体为：

判断三向应力状态下层理面的滑移和破坏情况；

判断三向压应力状态下层理面的滑移和破坏情况；

判断基岩的屈服和破裂状态，在三向挤压应力状态下，采用弹塑性屈服准则和塑形硬化准则判断岩石的屈服临界状态和破裂状态；

判断基岩中的裂缝倾角；

其中，所述判断三向应力状态下层理面的滑移和破坏情况，具体为：

当最小主应力σ₃<0或σ₃为张应力时且与地层的夹角为45°-90°时，采用公式(12)判断层理面是否发生滑移和破坏情况；

式中：F^t为岩石发生张破裂的应力条件，F^S为岩石发生剪破裂的应力条件，σ_t为层理面的抗张强度MPa，μ为岩石的泊松比，P为围压MPa，为岩石的内摩擦角°，C为岩石内聚力MPa；

所述判断三向压应力状态下层理面的滑移和破坏情况，具体为：

当水平挤压力σ₁与地层的夹角为30°-60°时，采用公式(13)判断层理面是否发生滑移和破坏；

式中：c_j为层理面的内聚力MPa，为层理面的内摩擦角°，θ_h为水平挤压力σ₁与地层的夹角°；

所述采用弹塑性屈服准则和塑形硬化准则判断岩石的屈服临界状态和破裂状态，具体为：

采用公式(14)判断岩石的屈服临界状态和破裂状态；

式中：f为岩石屈服临街状态和破裂状态的应力条件，I₁为应力第一不变量，J₂为应力第二不变量，a和K为与岩石内摩擦角和内聚力相关的实验常数，为岩石的内摩擦角°，σ₂为中间主应力MPa，σ_x为x方向正应力MPa，σ_y为y方向正应力MPa，σ_z为z方向正应力MPa，τ_xy为平行于xy面方向的剪应力MPa，τ_yz为平行于yz面方向的剪应力MPa，τ_xz为平行于xz面方向的剪应力MPa。

2.根据权利要求1所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述单井厚层泥岩段构造裂缝的分析和统计，具体为：

收集现场成像测井资料和钻井取芯资料，统计泥岩段地层产状和构造裂缝参数；

绘制单井裂缝发育柱状图，明确厚层泥岩段裂缝分布的差异性，针对不同裂缝发育程度的深度点采样进行镜下薄片观察，识别矿物类型及含量，分析其对裂缝密度的影响；

其中，所述地层产状包括层理倾角和裂缝产状，所述裂缝参数包括裂缝线密度参数、裂缝体密度参数、裂缝开度参数及裂缝长度参数。

3.根据权利要求1所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述进行岩石力学实验，具体为：

在泥岩厚层段，钻取柱塞岩样，且定义岩心钻取方向与轴向垂直，进行巴西劈裂试验、三轴压缩试验和直剪试验，测试并计算泥岩的Ⅰ型和Ⅱ型断裂韧性，同时记录各类试验的应力-应变曲线；

其中，所述三轴压缩试验用于获取岩石力学参数，所述岩石力学参数包括密度、弹性模量、泊松比、抗压强度、屈服强度、内摩擦系数及内摩擦角参数，所述巴西劈裂试验用于获得岩石的抗压强度，所述直剪试验用于获取岩石的抗剪强度、内摩擦系数以及内摩擦角；

在具有层理面的泥岩段，钻取柱塞岩样，且定义岩心钻取方向与轴向平行，进行巴西劈裂试验和直剪试验，同时记录各类试验的应力-应变曲线；

其中，所述巴西劈裂试验用于获得所述层理面的抗拉强度，所述直剪试验用于获得所述层理面的抗剪强度、内摩擦系数和内摩擦角。

4.根据权利要求3所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述力学参数测井解释，具体为：

通过全波列测井和偶极横波测井计算得到横波时差曲线，并与常规测井资料进行相关性分析，拟合泥岩段的横波时差预测公式(1)，并通过公式(2)、(3)、(4)、(5)、(6)计算单井动态岩石力学参数；

Δt_s＝0.948Δt_p+57.42 (1)

式中：Δt_s为横波时差μs/m，Δt_p为纵波时差μs/m；

S_c＝E_d0.008V_sh+0.0045(1-V_sh) (4)

式中：E_d为动态杨氏模量GPa，ρ_b为体积密度g/cm³，μ_d为动态泊松比，S_c为单轴抗压强度MPa，C为岩石内聚力MPa，为层理面的内摩擦角°，V_sh为泥质的体积含量；

基于所述单井动态岩石力学参数结果进行断裂韧性计算，得出针对泥岩的Ⅰ型和Ⅱ型断裂韧性计算模型(7)(8)；

K_ⅡC＝0.0466P+0.1674S_t-0.1851 (8)

式中：P为围压MPa，S_t为互层岩石的抗拉强度MPa。

5.根据权利要求4所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述计算脆性指数及恢复古力学参数，具体为：

通过公式(9)(10)对单井动态岩石力学模型进行连续矫正，制作单井静态力学参数图；

E_S＝2.902E_d-14.92 (9)

μ_S＝4.59μ_d-1.07 (10)

式中：E_S为静态杨氏模量GPa，μ_S为静态泊松比；

结合密度测井资料、阵列声波或偶极声波所测量的纵、横波数据和杨氏模量、泊松比参数，选择泊-杨法根据公式(11)计算泥岩段的脆性指数；

式中：BI_e为弹性参数脆性指数，E_max和E_min分别为某一层段杨氏模量的最大和最小值，μ_max和μ_min分别为某一层段泊松比的最大和最小值，E为杨氏模量GPa；

收集地震叠前反演得到的研究区的P波阻抗、S波阻抗和密度数据体，并以单井测井解释结果为约束，计算现今应力条件下的杨氏模量、泊松比，并运算得到泥岩段的脆性指数体；

收集所述研究区的深度域地震解释剖面，所述剖面垂直于主体构造走向，且所述剖面包含目的层和断裂解释结果，采用构造平衡剖面法恢复断距和二维褶皱形态，确定构造强烈活动期或关键造缝期；

建立所述研究区的泥岩层段构造地质模型，在划分断层相对级别的前提下，以一级、二级断层为控制，采用断层相关褶皱理论进行关键造缝期的古构造形态恢复；

通过三轴力学实验拟合建立不同围压条件下的力学参数变化曲线，将围压等效为埋藏深度下的上覆地层压力，根据所述古构造形态恢复的结果得到关键造缝期的泥岩层段古埋深，结合力学参数变化曲线计算该时期的古力学参数和脆性指数三维体，并按照脆性指数的分布趋势和区间特征将其进行等级划分。

6.根据权利要求5所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述泥岩构造裂缝力学参数模型推导，具体为：

选取岩石内表征单元体进行分析并简化，得到公式(15)；

式中：w_总为当前应力状态下所述单元体受力吸收的总应变能密度J/m³，w_d为岩石屈服后发生塑性变形所消耗的塑性应变能密度J/m³，w为存储的可释放弹性应变能密度J/m³，w_e为想产生裂缝必须克服的弹性应变能密度J/m³，V表征所述单元体的体积m³，S_f为新增裂缝的表面积m²，J为产生单位面积裂缝所需要的能量J/m²；

根据能量守恒原理和最大应变能密度理论，得到公式(16)；

式中：D_vf为所述单元体内缝隙的体密度m²/m³，μ为泊松比；

推导建立三向挤压应力状态下泥岩和层理面裂缝密度、开度的力学模型，采用拉-剪复合破裂、莫尔－库仑破裂准则和德鲁克-普拉格准则进行岩石破裂判别，若满足破裂条件，通过公式(17)进行裂缝参数计算；

式中：J₀为零围压时裂缝表面能J/m²，ε₀为最大弹性张应变，σ₃为最小张应变，K₀为比例系数，w为应变能密度J/m³，w_f为裂缝应变能密度J/m³，w_e为新增单位裂缝表面积必须克服的应变能J/m³，σ_p为岩石破裂应力MPa，D_lf为所述单元体内缝隙的线密度J/m，为岩石的内摩擦角°，θ为破裂角°，C₀为岩石的内聚力MPa，L₁、L₃为最大主应力、最小主应力方向的单元网络长度m，b为裂缝开度m；

其中，所述最大弹性张应变为对应于0.85σ_c时的应变，所述最小张应变为对应0.45σ_c时的张应变，σ_c为压缩强度极限MPa。

7.根据权利要求6所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述地质模型构建及古应力状态恢复，具体为：

收集地震解释、地质分层和测井解释结果，构建所述研究区的断层模型以及储层地质模型，以进行面模型重构，并通过面模型重构进行三维模型构建，完成三维地质力学模型的建立；

对所述研究区的钻井岩心进行声发射试验取样，每组样品至少包括4个试样；

将加工好的所述试样进行重复加载实验，加载速率为0.1MPa/s，测出所述试样在受载过程中的声发射信号随载荷的变化曲线，在第二次加载声发射载荷变化曲线上，找出声发射抹录不尽点，参照所述抹录不尽点的载荷值，在第一次加载声发射载荷曲线中，确定出凯塞尔点，取所述凯塞尔点和抹录不尽点载荷的平均值，即为岩心在地下所受的最大正应力，所述凯塞尔点的数量代表构造运动的期次；

根据构造形迹力学方法，进行裂缝的分期配套，进行裂缝走向统计分析，找出共轭裂缝的优势组系，以其共轭夹角平分线作为关键造缝期古应力场最大主应力方向；

其中，在所述声发射试验取样中，当试样的数量为4个时，其中一个所述试样取自品行井眼轴线的方向，其余所述试样取自水平面内互成45°角的三个方向。

8.根据权利要求7所述的基于弹塑性力学定量预测泥岩构造裂缝的方法，其特征在于，

所述三维古应力场模拟及裂缝定量预测，具体为：

基于所述三维地质力学模型和关键造缝期古应力测试结果，设置模型力学边界，设定加载速率为0.1MPa/s，进行弹塑性有限元三维应力场模拟，并反复尝试，直到古应力场的最大主应力分量与次级挤压断层分布趋势、最小主应力分量与次级伸展断层分布一致、剪应力分量与走滑断层分布趋势一致为止；

基于得到的破裂准则、屈服准则、裂缝密度、开度和产状力学模型通过软件获取每个节点的最大主应力、最下主应力、中间主应力、剪应力、最大主应变、最下主应变、中间主应变、屈服强度参数，并经过计算得到每个节点裂缝的体密度、线密度和开度值；

以基岩力学层为单位进行单井预测裂缝线密度、开度与实测值对比，当平均吻合率低于90％时，修正裂缝力学参数、调整边界条件重新模拟三维应力场，直至平均吻合率达到90％以上为止，完成厚层泥岩段构造裂缝的三维定量预测。