CN116611248A - 基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法 - Google Patents

Abstract

本发明公开了基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，具体为：收集深部各向异性储层岩石的物理力学参数、工程地质环境参数和双边缺口巴西圆盘几何参数；根据几何参数，绘制模拟不同层理角度真实剪切断裂的岩体计算域；对计算域添加固体力学场、历史应变场以及相场等多个物理场，建立相互耦合的各向异性真实剪切相场模型，并采用交错方案求解耦合，保持位移加载速率不变，重复上述步骤，直至岩石完全剪切断裂。该方法充分考虑岩石的各向异性特征以及工程地质环境特征，定量地表征了各向异性岩石真实剪切断裂模式，为评价深部地下工程施工过程中岩石的力学参数提供了重要参考。

Description

基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算 方法

技术领域

本发明属于岩土工程技术领域，具体涉及基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法。

背景技术

随着经济的不断发展，地下工程建设逐渐向更大深度扩展，高围压下的工程地质环境导致储层岩石更容易发生剪切破坏。例如，在非常规油气钻探和深部隧道开挖过程中，井壁以及洞室附近的围岩容易因应力释放而发生剪切破坏和失稳。同时，由于大部分地层是由不同岩性岩石组成的层状结构，储层岩石普遍具有各向异性特征，造成剪切裂缝扩展行为复杂。因此，了解各向异性岩石的剪切破裂机制对于地下能源的有效开发和地下岩石工程的安全至关重要。

大量文献表明，各向异性岩石的剪切破裂机制在很大程度上尚未被探索，主要是由于难以进行基于剪切的断裂试验。现有标准室内剪切裂纹扩展试验中通常会观察到裂纹出现一定角度的扭结，这种扭结裂纹的形成是由拉应力驱动的。因此，大部分实验无法准确了解岩石的真实剪切破裂机制。尽管已有报道提出的穿通剪切(PTS)试验和剪切箱试验可以捕捉到岩石的真实剪切断裂形态，但由于其装置的复杂型和额外要求，仅能获取有限的实验数据。因此，受限于室内试验的苛刻条件，开展各向异性岩石真实剪切断裂扩展的数值模拟是有必要的。

双边缺口巴西圆盘(DNBD)的径向压缩作为一种测试真实剪切断裂的新兴方法，不仅操作简单，还可以确定任何方向上相对于各向异性方向的真实剪切断裂韧性。然而，目前基于双边缺口巴西圆盘的剪切断裂扩展模型尚未考虑岩石的各向异性和深部工程地质环境对岩石真实剪切裂纹扩展所带来的影响，存在一定局限性。

发明内容

本发明的目的是提供基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，综合考虑了深部各向异性储层岩石的工程环境特征，从而准确预测岩石真实剪切断裂形态。

本发明所采用的技术方案是，基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，具体按照以下步骤实施：

步骤1、收集深部各向异性储层岩石的物理力学参数、地下工程地质环境参数和双边缺口巴西圆盘几何参数；

步骤2、根据物理力学参数中获得的不同主方向下的剪切断裂韧性，计算出不同层理角度下岩石纯剪切断裂的临界能量释放率G_c，公式如式(1)及式(2)所示：

K_IIc＝K_IIc,1 sin²α₁+K_IIc,2 cos²α₁ (1)；

式中，K_IIc,1为垂直于层理方向剪切断裂韧性；K_IIc,2为沿层理方向剪切断裂韧性；E₁为垂直于层理方向弹性模量；ζ＝E₁/E₂，E₂为沿层理方向弹性模量；α₁为层理角度；v₁₂和v₂₃分别为沿层理平面内和垂直层理平面内的泊松比；

步骤3、利用COMSOL软件绘制双边缺口巴西圆盘几何计算域，并利用COMSOL内置基矢坐标系表征岩石不同层理角度下的横观各向同性特征；

步骤4、对计算域分别添加固体力学场、历史应变场以及相场，建立相互之间的耦合关系，从而获得各向异性真实剪切相场模型；

步骤5、对各向异性真实剪切相场模型采用交错格式在分离耦合求解器上独立求解；

步骤6、保持位移加载速率不变，重复步骤5，直至岩石完全断裂后停止计算，即可获得不同层理角度时岩石的真实剪切断裂形态。

本发明的特点还在于，

步骤1中，物理力学参数包括：岩石沿层理方向和垂直层理方向的弹性模量、剪切模量、泊松比、剪切断裂韧性以及岩石的比奥系数；工程地质环境参数包括：岩石所受围压、孔隙压力以及原位应力；双边缺口巴西圆盘几何参数包括：圆盘的半径、预制缺口尺寸、加载角度。

步骤3中，绘制双边缺口巴西圆盘几何计算域的具体步骤为：首先以坐标原点为中心生成半径R为37.5mm的圆形以及94mm×1.8mm的矩形，随后以坐标形式建立两条插值曲线用以确定预制缺口长度，将所有几何对象通过布尔并集合并后删除两侧预制缺口区域，最后获得双边缺口巴西圆盘模型计算域；

建立内置基矢坐标系，如式(3)所示：

式中，x和y分别表示全局坐标系x轴方向和y轴方向；x′和y′分别表示垂直于层理方向和沿层理方向；α表示沿层理方向与y轴方向的夹角，分别选取五个不同层理角度α＝0°、22.5°、45°、67.5°以及90°用以真实剪切模拟分析。

步骤4中，各向异性真实剪切相场模型包括固体力学方程、相场控制方程和历史应变方程；

其中，固体力学方程如式(4)所示：

式中，σ＝[(1-κ)(1-φ)²+κ]￡:ε，0＜κ≤1用以避免数值模型的奇异性，φ为裂纹萌生以及扩展的相场，￡为未退化的四阶刚度张量，ε为应变张量；f_b和f_s分别为岩石的体力和面力；n为方向矢量；▽为梯度算子；

相场控制方程如式(5)所示：

式中，l₀为长度尺寸参数；G_c为临界能量释放率；B为各向异性梯度；弹性应变能的拉伸部分

其中，ε_ii为线应变，ε_ij(i≠j)为切应变；i＝1,2,3、j＝1,2,3；

▽φ为相场梯度；C^l为未退化的四阶刚度张量￡的Voight表示法，其表示为：

式中，为刚度分量，其不同表达式如下：

式中，E_i(i＝1,2,3)表示沿三个主要材料方向的弹性模量(e₁,e₂,e₃)；G_ij(i＝1,2,3、j＝1,2,3)是由e_i和e_j定义的正交平面内的剪切模量；v_ij＝v_jiE_i/E_j为正交平面中的泊松比；

历史应变方程如式(6)所示：

式中，t为时间；x为二维域任意点的位置矢量；σ₀为初始应力张量。

步骤5中，具体为：

将步骤1的深部各向异性储层岩石的物理力学参数、地下工程地质环境参数代入各向异性真实剪切相场模型中，依据步骤3给定模型不同的层理角度，将圆盘计算域上方设置为指定位移边界用以模拟压裂过程，下方设置为辊支撑，其余为自由边界，指定位移边界输入位移增量Δu＝0.05mm/min，按照位移场-历史应变场-相场的顺序进行分离耦合瞬态计算，分别得到应力、应变和相场分布。

本发明的有益效果是：针对岩石在深部地下工程地质环境中易发生剪切破坏，建立各向异性纯固体相场模型，提供预测岩石真实剪切断裂的计算方法，由于该方法充分考虑岩石的各向异性特征以及工程地质环境特征，可以准确预测岩石真实剪切断裂形态，定量评估剪切断裂面的形态参数差异。其计算结果不仅为评价深部地下工程施工过程中岩石的力学参数提供了重要参考，而且对准确选择断裂面研究不同工况下的强度和变形也具有重要意义。

附图说明

图1是本发明基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法的流程图；

图2为本发明方法中双边缺口巴西圆盘计算域及边界条件示意图；

图3为本发明各向异性纯固体相场模型中多物理场耦合方式图；

图4为本发明实施例中考虑不同层理角度下的剪切裂缝形与实验结果对比图；

图5为本发明实施例中最大主应变的数值结果与实验结果对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

实施例1

本发明基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，具体为：收集深部各向异性储层岩石的物理力学参数、工程地质环境参数和双边缺口巴西圆盘几何参数；根据几何参数，绘制模拟不同层理角度真实剪切断裂的岩体计算域；对计算域添加固体力学场、历史应变场以及相场等多个物理场，建立相互耦合的各向异性真实剪切相场模型；采用交错方案求解耦合系统，通过历史应变能更新相场从而驱动裂纹萌生及扩展；保持位移加载速率不变，重复上述步骤，直至岩石完全剪切断裂。

本发明方法充分考虑岩石的各向异性特征以及工程地质环境特征，定量地表征了各向异性岩石真实剪切断裂模式，不仅为评价深部地下工程施工过程中岩石的力学参数提供了重要参考，而且对准确选择断裂面研究不同工况下的强度和变形也具有重要意义。

实施例2

本发明基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，如图1所示，具体按照以下步骤实施：

步骤1、收集深部各向异性储层岩石的物理力学参数、地下工程地质环境参数和双边缺口巴西圆盘几何参数；

其中物理力学参数包括：岩石两个主方向(沿层理方向和垂直层理方向)的弹性模量、剪切模量、泊松比、剪切断裂韧性以及岩石的比奥系数；

工程地质环境参数包括：岩石所受围压、孔隙压力以及原位应力；

双边缺口巴西圆盘几何参数包括：圆盘的半径、预制缺口尺寸、加载角度；

步骤2、根据物理力学参数中获得的不同主方向下的剪切断裂韧性，计算出不同层理角度下岩石纯剪切断裂的临界能量释放率G_c，公式如式(1)及式(2)所示：

K_IIc＝K_IIc,1 sin²α₁+K_IIc,2 cos²α₁ (1)；

式中，K_IIc,1为垂直于层理方向剪切断裂韧性；K_IIc,2为沿层理方向剪切断裂韧性；E₁为垂直于层理方向弹性模量；ζ＝E₁/E₂，E₂为沿层理方向弹性模量；α₁为层理角度；v₁₂和v₂₃分别为沿层理平面内和垂直层理平面内的泊松比；

步骤3、根据几何参数，利用COMSOL软件绘制双边缺口巴西圆盘几何计算域，并利用COMSOL内置基矢坐标系表征岩石不同层理角度下的横观各向同性特征；

绘制双边缺口巴西圆盘几何计算域的具体步骤为：首先以坐标原点为中心生成半径R为37.5mm的圆形以及94mm×1.8mm的矩形，随后以坐标形式建立两条插值曲线用以确定预制缺口长度，将所有几何对象通过布尔并集合并后删除两侧预制缺口区域，最后获得双边缺口巴西圆盘模型计算域，如图2所示。

建立内置基矢坐标系，如式(3)所示：

式中，x和y分别表示全局坐标系x轴方向和y轴方向；x′和y′分别表示垂直于层理方向和沿层理方向；α表示沿层理方向与y轴方向的夹角，分别选取五个不同层理角度α＝0°、22.5°、45°、67.5°以及90°用以真实剪切模拟分析；

步骤4、对计算域分别添加固体力学场、历史应变场以及相场等多个物理场，建立相互之间的耦合关系，从而获得各向异性真实剪切相场模型，如图3所示；

使用COMSOL中的固体力学方程、历史应变方程和相场方程组成了各向异性岩石真实剪切断裂的相场模型。首先从固体力学场中提取拉伸弹性能量密度储存在变量列表中。随后，通过历史应变场调用拉伸弹性能量密度/>进而更新局部历史应变场H。采用更新后的H求解相场，更新后的相场用于修改固体力学场中的刚度矩阵。通过使用历史字段H，可以确保相位字段在压缩或卸载期间单调增加，从而避免裂纹愈合。

各向异性真实剪切相场模型包括固体力学方程、相场控制方程和历史应变方程；

其中，固体力学方程如式(4)所示：

式中，0＜κ≤1用以避免数值模型的奇异性，φ为裂纹萌生以及扩展的相场，￡为未退化的四阶刚度张量，ε为应变张量；f_b和f_s分别为岩石的体力和面力；n为方向矢量；▽为梯度算子；

相场控制方程如式(5)所示：

式中，l₀为长度尺寸参数；G_c为临界能量释放率；B为各向异性梯度；弹性应变能的拉伸部分

其中，ε_ii为线应变，ε_ij(i≠j)为切应变；i＝1,2,3、j＝1,2,3；

▽φ为相场梯度；C^l为未退化的四阶刚度张量￡的Voight表示法，其表示为：

式中，为刚度分量，其不同表达式如下：

式中，E_i(i＝1,2,3)表示沿三个主要材料方向的弹性模量(e₁,e₂,e₃)；G_ij(i＝1,2,3、j＝1,2,3)是由e_i和e_j定义的正交平面内的剪切模量；v_ij＝v_jiE_i/E_j为正交平面中的泊松比。

历史应变方程：

为了真实模拟地下工程地质环境，必须充分考虑岩石所受围压、孔隙压力以及原位应力的影响，则相应的历史应变方程如式(6)所示：

式中，t为时间；x为二维域任意点的位置矢量；σ₀为初始应力张量。

步骤5、对各向异性真实剪切相场模型采用交错格式在分离耦合求解器上独立求解；

具体为：将步骤1的深部各向异性储层岩石的物理力学参数、地下工程地质环境参数代入各向异性真实剪切相场模型中，依据步骤3给定模型不同的层理角度(即α＝0°、22.5°、45°、67.5°以及90°)，将圆盘计算域上方设置为指定位移边界用以模拟压裂过程，下方设置为辊支撑，其余为自由边界，指定位移边界输入位移增量Δu＝0.05mm/min，按照位移场-历史应变场-相场的顺序进行分离耦合瞬态计算，分别得到应力、应变和相场分布。

步骤6、保持位移加载速率不变，重复步骤5，直至岩石完全断裂后停止计算，即可获得不同层理角度时岩石的真实剪切断裂形态。

实施例3

以某双边缺口巴西圆盘(DNBD)实验为例，利用表1所示的各向异性花岗岩的实验参数模拟岩石真实剪切断裂形态。按照上面步骤开展模拟计算，最终获得岩石的真实剪切断裂形态。随后，将模拟结果与实验结果进行对照，如图4所示，可以看出，数值计算结果与实验结果吻合较好。同时，还将最大主应变的数值结果与实验结果进行比较，如图5所示，这与实验结果非常一致。基于以上对花岗岩的破坏形态和最大主应变的比较结果，表明本计算方法在模拟真实剪切断裂方面是有效的。

表1某各向异性花岗岩的实验参数

(注：储层系数中“1”和“2”分别代表沿层理方向以及垂直层理方向；“12”和“23”分别代表各向异性平面以及各向同性平面。)

本发明针对岩石在深部地下工程地质环境中易发生剪切破坏，建立了各向异性纯固体相场模型，提供了一种预测岩石真实剪切断裂的计算方法。由于该方法充分考虑岩石的各向异性特征以及工程地质环境特征，可以准确预测岩石真实剪切断裂形态，定量地表征了各向异性岩石真实剪切断裂模式。计算结果不仅为评价深部地下工程施工过程中岩石的力学参数提供了重要参考，而且对准确选择断裂面研究不同工况下的强度和变形也具有重要意义。

Claims (6)

1.基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，具体为：收集深部各向异性储层岩石的物理力学参数、工程地质环境参数和双边缺口巴西圆盘几何参数；根据几何参数，绘制模拟不同层理角度真实剪切断裂的岩体计算域；对计算域添加多个物理场，建立相互耦合的各向异性真实剪切相场模型；采用交错格式在分离耦合求解器上独立求解；保持位移加载速率不变，重复上述步骤，直至岩石完全剪切断裂，即可获得不同层理角度时岩石的真实剪切断裂形态。

2.根据权利要求1所述的基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，具体按照以下步骤实施：

步骤1、收集深部各向异性储层岩石的物理力学参数、地下工程地质环境参数和双边缺口巴西圆盘几何参数；

步骤2、根据物理力学参数中获得的不同主方向下的剪切断裂韧性，计算出不同层理角度下岩石纯剪切断裂的临界能量释放率G_c，公式如式(1)及式(2)所示：

K_IIc＝K_IIc,1 sin²α₁+K_IIc,2 cos²α₁ (1)；

式中，K_IIc,1为垂直于层理方向剪切断裂韧性；K_IIc,2为沿层理方向剪切断裂韧性；E₁为垂直于层理方向弹性模量；ζ＝E₁/E₂，E₂为沿层理方向弹性模量；α₁为层理角度；v₁₂和v₂₃分别为沿层理平面内和垂直层理平面内的泊松比；

步骤3、利用COMSOL软件绘制双边缺口巴西圆盘几何计算域，并利用COMSOL内置基矢坐标系表征岩石不同层理角度下的横观各向同性特征；

步骤4、对计算域分别添加固体力学场、历史应变场以及相场，建立相互之间的耦合关系，从而获得各向异性真实剪切相场模型；

步骤5、对各向异性真实剪切相场模型采用交错格式在分离耦合求解器上独立求解；

步骤6、保持位移加载速率不变，重复步骤5，直至岩石完全断裂后停止计算，即可获得不同层理角度时岩石的真实剪切断裂形态。

3.根据权利要求2所述的基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，所述步骤1中，物理力学参数包括：岩石沿层理方向和垂直层理方向的弹性模量、剪切模量、泊松比、剪切断裂韧性以及岩石的比奥系数；工程地质环境参数包括：岩石所受围压、孔隙压力以及原位应力；双边缺口巴西圆盘几何参数包括：圆盘的半径、预制缺口尺寸、加载角度。

4.根据权利要求2所述的基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，所述步骤3中，绘制双边缺口巴西圆盘几何计算域的具体步骤为：首先以坐标原点为中心生成半径R为37.5mm的圆形以及94mm×1.8mm的矩形，随后以坐标形式建立两条插值曲线用以确定预制缺口长度，将所有几何对象通过布尔并集合并后删除两侧预制缺口区域，最后获得双边缺口巴西圆盘模型计算域；

建立内置基矢坐标系，如式(3)所示：

式中，x和y分别表示全局坐标系x轴方向和y轴方向；x′和y′分别表示垂直于层理方向和沿层理方向；α表示沿层理方向与y轴方向的夹角，分别选取五个不同层理角度α＝0°、22.5°、45°、67.5°以及90°用以真实剪切模拟分析。

5.根据权利要求2所述的基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，所述步骤4中，各向异性真实剪切相场模型包括固体力学方程、相场控制方程和历史应变方程；

其中，固体力学方程如式(4)所示：

式中，σ＝[(1-κ)(1-φ)²+κ]￡:ε，0＜κ≤1用以避免数值模型的奇异性，φ为裂纹萌生以及扩展的相场，￡为未退化的四阶刚度张量，ε为应变张量；f_b和f_s分别为岩石的体力和面力；n为方向矢量；▽为梯度算子；

相场控制方程如式(5)所示：

式中，l₀为长度尺寸参数；G_c为临界能量释放率；B为各向异性梯度；弹性应变能的拉伸部分

其中，ε_ii为线应变，ε_ij(i≠j)为切应变；i＝1,2,3、j＝1,2,3；

▽φ为相场梯度；C^l为未退化的四阶刚度张量￡的Voight表示法，其表示为：

式中，为刚度分量，其不同表达式如下：

式中，E_i(i＝1,2,3)表示沿三个主要材料方向的弹性模量(e₁,e₂,e₃)；G_ij(i＝1,2,3、j＝1,2,3)是由e_i和e_j定义的正交平面内的剪切模量；v_ij＝v_jiE_i/E_j为正交平面中的泊松比；

历史应变方程如式(6)所示：

式中，t为时间；x为二维域任意点的位置矢量；σ₀为初始应力张量。

6.根据权利要求5所述的基于巴西圆盘实验考虑的各向异性岩石真实剪切断裂计算方法，其特征在于，所述步骤5中，具体为：

将步骤1的深部各向异性储层岩石的物理力学参数、地下工程地质环境参数代入各向异性真实剪切相场模型中，依据步骤3给定模型不同的层理角度，将圆盘计算域上方设置为指定位移边界用以模拟压裂过程，下方设置为辊支撑，其余为自由边界，指定位移边界输入位移增量Δu＝0.05mm/min，按照位移场-历史应变场-相场的顺序进行分离耦合瞬态计算，分别得到应力、应变和相场分布。