CN113569442A - 一种基于rkpm-pd耦合算法的岩石裂纹扩展预测方法 - Google Patents

Abstract

本发明公开了一种基于RKPM‑PD耦合算法的岩石裂纹扩展预测方法。与原有算法相比，RKPM‑PD耦合算法继承了基于PD的非局部特性，不需要对裂纹尖端进行动态追踪；另一方面，该耦合算法有效解决了PD方法中的数值震荡现象，具备较好的数值稳定性。该算法在位移连续处采用RKPM方法，并应用PD中键的思想处理裂纹，有效结合了两方法的优势，非常适用于求解岩石材料内复杂三维裂纹扩展问题。

Description

一种基于RKPM-PD耦合算法的岩石裂纹扩展预测方法

技术领域

本发明专利属于岩土工程数值计算领域，尤其涉及一种基于重构核函数(Reproducing kernel particle method,RKPM)与近场动力学(Peridynamic,PD)耦合 算法的裂纹扩展预测方法。

背景技术

在岩石材料中存在许多的微裂隙，在工程尺度，岩体中又含有大量的节理与 裂隙。岩体中的这些软弱结构面，是岩石工程失稳破坏的主要影响因素。因此， 预测岩体中裂纹的扩展过程对于揭示岩体的变形和破坏规律以及评价岩土工程 的安全可靠性具有非常重要的意义。同时，高效可控的岩石破裂技术，也是煤层 气及页岩气开采中的关键技术。

对于岩石裂纹扩展问题，由于其强烈的几何非线性与材料非线性，往往无法 通过解析的形式求解，只能通过数值方法进行模拟研究。传统的有限元方法在模 拟裂纹扩展问题时，由于新裂纹与原始单元界面不一致，需要不断地重新划分网 格，极大地增加了计算工作量。对于复杂的三维裂纹模型，网格的生成也极具挑 战性。

鉴于有限元方法中网格重构的缺陷，基于节点离散的无网格重构核函数方法(Reproducing kernel particle method,RKPM)被广泛用于计算裂纹扩展问题。该方 案不再需要对网格进行划分与重构，降低了求解裂纹扩展问题的难度。然而，目 前RKPM方法对裂纹的处理方案仍基于断裂力学准则，需要对裂纹前端进行动 态的捕捉。但数值计算中如何实现对裂纹，特别是三维裂纹前端的动态追踪，至 今仍未有较好的解决方案。

近年来，一种基于非局部思想的近场动力学(Peridynamic,PD)方法开始在国 内外受到关注。该方法以非局部作用的积分模型代替传统理论的微分模型，避免 了传统方法在位移不连续处求导产生的应力奇异性问题，因而特别适用于裂纹等 不连续问题的求解。

在PD方法中，裂纹不再需要显式的追踪，而是通过键的断裂表示。定义距 节点小于δ范围的区域为该节点的支撑域，则节点与其支撑域内其他节点的相互 作用即被称为键。如图1所示，圆形区域代表节点x_i的支撑域，处于该支撑域内 的节点x_j与节点x_i间的黑色实线即代表两节点间的键。当键的临界伸长率超过一 定极限时，两点间的键发生断裂，若干键的断裂即可勾勒出裂纹的扩展路径。

但另一方面，目前的PD方法在应用时会遇到零能模式(zero-energy mode)这 一问题，该问题会导致PD方法在计算时产生数值震荡现象，阻碍近场动力学方 法更进一步的应用。

综上所述，目前的岩石裂纹扩展预测方法，均存在一定的缺陷。RKPM方法 计算稳定，但需要对裂纹前沿不断追踪，难以应用于三维问题。PD方法摆脱了 传统断裂力学的束缚，但又因零能模式问题而饱受数值震荡现象的困扰。因此， 目前缺乏一种基于RKPM-PD耦合计算的新型裂纹扩展预测方法。

发明内容

本发明的目的是针对现有技术存在的上述不足，提出一种基于RKPM-PD耦 合算法的裂纹扩展预测方法。其中，无裂纹部分采用RKPM方法进行计算，裂 纹部分采用PD中键的思路进行处理。

本发明通过以下技术方案实现：

一种基于RKPM-PD耦合算法的岩石裂纹扩展预测方法，其特征在于：

(1)首先将岩石物理模型离散为带背景网格的节点模型，生成带自由度的 节点与带权重的高斯点，确定各点支撑域内的节点信息，所述各点包括所述节点 和高斯点；

(2)依据输入的初始裂纹空间坐标信息，建立初始裂纹的PD模型；

(3)采用RKPM近似函数形成刚度矩阵，计算出每一时步各点的位移与应 力状态；

(4)施加边界条件，基于增量格式的Newmark算法进行迭代求解；对于每 一时步计算出的应力值，判断各点间键的断裂状态，如果有断裂发生，更新该点 处的刚度矩阵，并重新迭代计算，直至该时步内外力平衡；在所有时步施加完毕 后，得到基于RKPM-PD耦合算法进行裂纹扩展预测的结果。

包括如下具体步骤：

步骤一、首先在Abaqus内建立模型的有限元背景网格，并通过inp文件导入Matlab内，通过背景网格上的一阶有限元插值函数计算出节点上的体积，高斯点 处的坐标与权重信息；

步骤二、给定各节点与高斯点的支撑域半径，通过kdtree及rangesearch函 数得到各点支撑域内的节点信息；

步骤三、输入模型初始裂纹面信息，将各点支撑域内穿过初始裂纹面的键预 先断掉；与传统PD法不同，这里的键不仅指支撑域节点与节点间的相互作用， 也指支撑域内高斯点与节点间的相互作用；

步骤四，依据步骤二与步骤三得到的前处理信息，计算各点形函数及其导数， 并依据各高斯点的权重，得到该处的节点刚度矩阵，并组装为全局刚度矩阵；

步骤五，将模型总的位移荷载或面力荷载划分为若干增量步，进行增量格式 的Newmark时间迭代计算；

步骤六，在每一时步中，依据计算得到的应力，通过莫尔-库伦等岩石强度准 则，判断节点与节点及高斯点与节点间的键是否有新的断裂；如果有新键断裂， 更新该点处刚度矩阵，并重新迭代计算，直至该时步计算平衡；

步骤七，迭代完成后，计算各点损伤值，得到岩石裂纹扩展的预测结果。

本发明中，裂纹处理方案可以采用：

相较于PD方法，RKPM-PD耦合算法的裂纹处理方案虽然同样基于节点间 的“键”，但形式上更为复杂。如图2所示，RKPM-PD算法中既存在节点与节 点之间的键(浅黑色实线)，也存在高斯点与节点之间的键(浅色实线)。因此， 在采用该耦合算法计算裂纹时，不仅需要处理节点与节点间的键(深黑色间断线)， 还需处理高斯点与节点间的键(浅色间断线)。

本发明中，键的断裂准则可以采用：

对于岩石类材料，为充分引入岩石类材料的本构模型，本发明基于应力准则 来判断键的断裂情况。定义键上的最大主应力为其两端点最大主应力的平均值， 即：

当键上的最大主应力超过材料容许的最大强度f_t时，该键发生断裂。对于岩 石材料的受拉破坏情况，f_t可通过最大拉应力准则来确定，对于岩石材料的剪切 破坏情况，可选用莫尔-库伦准则作为岩石材料的剪切破坏准则。

根据上述键的破坏准则，节点x的损伤定义为：

物质点损伤程度

的范围为

0代表完整无损，1代表物质点影响域内 的键完全断裂。根据物质点的损伤程度，即可表示出裂纹的扩展路径。

本发明中，时间迭代格式可以采用：

为提高计算时间迭代的计算精度，本发明采用增量格式Newmark算法进行 时间迭代。增量格式下的Newmark方程可写为：

其中：

根据全量格式下Newmark法中速度与加速度的关系，可得：

其中

式中，α与δ分别为Newmark法中参数。当δ≥0.5，α≥0.25(0.5+δ)²时， 算法是无条件稳定的。

将式(5)带入式(3)中，可得：

即

其中：

在计算出ΔU后，即可通过U_t+Δt的值得到出新的σ_t+Δt。

本发明中，裂纹更新方案：

利用公式(8)中得到的位移U_t+Δt可进一步计算出σ_t+Δt，之后可通过公式(1)判 断是否有新的键断掉。根据是否有新键断裂，后续计算可分为：

a)如果当前步中没有新的断键出现，说明t+Δt时刻方程达到平衡，可进入 下一步计算。在通过式(4)与式(5)得到t+Δt时刻的速度及加速度后，进行变量回 代，并开启下一时步的计算。

b)如果当前步中有新键断裂，说明t+Δt状态并未达到平衡，方程(8)仍需进 行迭代计算。此时，对于支撑域内出现新断键的物质点，需要更新其状态，并重 新计算该点上的形函数[N]及对应的包含一阶导数的矩阵[B]。在新的断裂出现后， 体系内的不平衡力R为：

其中，[σ]为断键后的应力，可通过更新后的[B]与U_t+Δt得到。同时，新的 有效刚度矩阵

为：

[K]为新的全局刚度矩阵。新的位移增量[ΔU]为：

t+Δt时刻的位移为：

U_t+Δt＝U_t+ΔU+[ΔU] (13)

根据上式的U_t+Δt重新计算出σ_t+Δt，并继续判断是否依然有新的断键出现， 如果继续出现，则返回b)中继续迭代，直至系统内不再出现新的断键。此时， 结束当前时步内的循环，转至a)中。

本发明中，刚度矩阵更新说明：

在更新刚度矩阵[K]时，为降低计算量，新的刚度矩阵[K]可通过前一步刚度 度矩阵K与当前所有待更新物质点x_a上节点刚度矩阵的变化求得，可写为：

式中，上标′old′代表上一迭代步中的变量，K_a代表当前待更新物质点x_a上 的节点刚度矩阵，m代表需要更新的物质点数目。通过上式，全局刚度矩阵更新 时的计算量可大幅降低。

本发明基于RKPM-PD耦合算法，提出了一种新型的裂纹扩展预测方法。与传 统RKPM方法相比，本发明能够有效处理三维裂纹扩展问题；与近场动力学相比， 本发明采用的高斯积分方案计算精度更高，不存在零能模式引起的数值震荡现象。

附图说明

图1是PD方法键与裂纹示意图。

图2是RKPM-PD键与裂纹示意图。

图3是本发明的算法流程图。

图4是四点弯曲算例示意图。

图5是四点弯曲算例裂纹扩展数值及试验结果图，其中(a)为RKPM-PD 算法的裂纹扩展结果，(b)为DH-PD算法的裂纹扩展结果，(c)为裂纹扩展的 室内试验结果。

图6是带预制裂纹的三维巴西圆盘示意图。

图7是含完全贯穿裂纹三维巴西圆盘的裂纹扩展路径，其中，(a)裂纹初始 位置(Step300)，(b)裂纹扩展图(Step340)，(c)裂纹扩展图(Step400)，(d)裂纹扩展 最终位置(Step420)。

图8是含完全贯穿裂纹三维巴西圆盘在NMM法及室内试验下的裂纹扩展路 径，其中，(a)为3D NMM方法下的裂纹扩展结果，(b)三维巴西圆盘裂纹扩 展的室内试验结果。

图9是含未贯穿裂纹(嵌入深度6mm)三维巴西圆盘的裂纹扩展路径；其中， (a1)：Step470,正面图(z＝0)；(a2)：Step470,背面图(z＝10mm)；(b1)：Step500,正 面图(z＝0)；(b2)：Step500,背面图(z＝10mm)；(c1)：Step550,正面图(z＝0)；(c2)： Step550,背面图(z＝10mm)。

图10是三维巴西圆盘算例中裂纹扩展过程的应力-应变曲线。

具体实施方式

本发明提供的一种基于RKPM-PD耦合算法预测岩石裂纹扩展路径的计算 方法，其计算流程如图3所示，具体包括如下步骤：

1)、离散模型，确定节点与高斯点位置坐标，节点体积与高斯点权重等信息。

本发明采用Matlab语言进行编程，前处理工作采用Abaqus有限元软件进行 辅助。具体来说，计算模型可通过Abaqus等有限元商业软件得到，并通过inp 文件导入Matlab内。随后，通过背景网格上的一阶插值函数计算出节点与高斯 点的坐标、体积或权重等信息。

2)、确定各点(节点+高斯点)支撑域内的节点信息；

RKPM-PD耦合算法要求找到节点与高斯点(积分点)支撑域内的节点信息。 该步可采用MATLAB内置kdtree函数进行简化计算。

3)、输入初始裂纹信息；

当模型存在初始裂纹时，输入初始裂纹信息，将穿过初始裂纹面的键全部断 掉。

4)确定各点的形函数及其导数，并形成全局刚度矩阵K；

a.计算高斯点处的RKPM形函数N及其导数；

b.根据公式，计算单元刚度矩阵；

c.组装刚度矩阵；

在该步中，可采用MATLAB的稀疏矩阵函数sparse提高计算速度。

5)施加边界条件，设置位移或外力加载条件；为便于实现，位移边界条件 可采用罚函数法进行施加。

6)计算位移，得到各点处的位移，应力与应变。

7)检验各点处键的断裂情况；

如果有断裂发生，更新断裂点处形函数及其导数，形成新的全局刚度矩阵， 再转入6)；如果没有，转到5)，进入下一增量步。

8)完成规定增量步，退出程序。

本发明通过如下两个例子进一步展示了本发明预测结果的准确性，特别是本 发明对于三维裂纹扩展问题的有效性和适用性。

(1)四点弯曲试验

四点弯曲模型经常被作为裂纹扩展问题的基准算例。如图4所示，该算例模 型参数为：L＝800mm，a＝40m，v＝0.18，E＝27Pa。

该算例采用平面应力假设，抗拉强度f＝2MPa，粘聚力c＝1.4MPa，内摩擦角 θ＝35°.模型采用间距h＝1mm的均匀网格离散，并采用Δt＝5×10^-5s的增量Newmark隐式算法，每个时步内新增压力为100N.图5(a)给出了RKPM-PD算法 下四点弯曲算例的最终裂纹扩展路径。从图中可以看到，两条初裂纹互相远离， 并近似地沿弧形路径向载荷点或固定位移端延伸。对比图5(b)中态基PD得到的 结果与图5(c)给出的试验结果，可以看到RKPM-PD算法得到的最终裂纹路径， 与其他数值方法及物理实验的结果相吻合。需要指出的是，图5(b)中的态基PD 方法在模拟该算例时首先需要在模型边界上布置一定厚度的虚拟节点来消除边 界效应的影响，这种处理方式在模拟不规则模型边界时会存在一定的困难。

(2)含预制裂纹的三维巴西圆盘试验

实际工程中遇到的岩石裂纹，大都无法简单地化为平面二维问题，其本质都 是三维的。因此，三维裂纹的扩展模拟更具有实际意义。但由于三维裂纹扩展问 题的复杂性，目前对于该问题仍未有较好的解决方案。本发明提出了RKPM-PD 算法，可有效预测三维裂纹的裂纹扩展路径。

基于RKPM-PD算法，本发明对图6所示的含预制裂纹的三维巴西圆盘试验 了进行模拟。模型在xoy平面内的直径为100mm，z方向的厚度为10mm(厚度方 向为垂直于纸面的方向)。圆盘上下两端受压平台长度为6.2mm，采用增量 Newmark隐式积分方案，时间步长取为Δt＝5×10^-7s，平均加载速度为 1×10^-8m/step。岩石试样力学参数为：弹性模量E＝10Pa，泊松比v＝0.25，密 度ρ＝2500kg/m³，抗拉强度T＝0.5MPa，粘聚力C＝5MPa以及内摩擦角φ＝ 40°。离散模型中平均节点间距h＝2mm，共计15288个节点，98840个高斯点。 在本算例中，预制裂纹长度为30mm，与水平方向呈45°夹角。为探究预制裂纹 深度对结果的影响，预制裂纹深度分别取为10mm，6mm及2mm。

当预制裂纹深度为10mm时，预制裂纹在厚度方向完全贯穿平板，模型在z 方向完全一致。该算例的裂纹扩展路径如图7所示，可以看到，裂纹扩展路径与 二维情况下相一致，裂纹尖端分别向圆盘上下两端的压缩载荷施加处扩展，直至 最后裂纹完全贯穿整个圆盘。该模型采用RKPM-PD的计算结果与在三维NMM 方法下的模拟结果及试验结果(图8)均保持一致。

通过对裂纹完全贯穿试件的分析，证明了该方法在三维裂纹问题中的有效性。 之后，考虑裂纹未完全贯穿的情况，设置裂纹嵌入深度分别为6mm与2mm，预 制裂纹从正面嵌入，终止于圆盘内部。为更好的显示裂纹未完全贯穿圆盘情况下 的裂纹扩展预测路径，对于同一时刻的裂纹扩展结果分别给出其正面图(z＝0)与 背面图(z＝10mm)。基于RKPM-PD预测的嵌入深度为6mm时的裂纹扩展路径如 图9所示。

通过对未穿透裂纹的结果分析可知，裂纹首先沿厚度方向扩展，从背面可逐 步的看到损伤区域的形成；之后，正面与背面的裂纹均沿竖直方向朝两端扩展。 对比两个面的裂纹最终扩展轨迹，可以看到两者大体一致，细微处有所差别。同 时，可以看到，裂纹切割深度越浅，启裂所需计算步越多。图10给出了三种不 同切割深度所对应的应力应变曲线，从图中可以清晰的看出，巴西圆盘在表面裂 纹(2mm)时所能承受的载荷远远大于深度裂纹(6mm)及完全贯穿裂纹(10mm)。

Claims (5)

1.一种基于RKPM-PD耦合算法的岩石裂纹扩展预测方法，其特征在于：

(1)首先将岩石模型离散为带背景网格的节点模型，生成带自由度的节点与带权重的高斯点，确定各点支撑域内的节点信息，所述各点包括所述节点和高斯点；

(2)依据输入的初始裂纹空间坐标信息，建立初始裂纹的PD模型；

(3)采用RKPM近似函数形成刚度矩阵，计算出每一时步各点的位移与应力状态；

(4)施加边界条件，基于增量格式的Newmark算法进行迭代求解；对于每一时步计算出的应力值，判断各点间键的断裂状态，如果有断裂发生，更新该点处的刚度矩阵，并重新迭代计算，直至该时步内外力平衡；在所有时步施加完毕后，得到基于RKPM-PD耦合算法进行裂纹扩展预测的结果。

2.如权利要求1所述的一种基于RKPM-PD耦合算法的岩石裂纹扩展预测方法，其特征在于，包括如下具体步骤：

步骤一、首先在Abaqus内建立模型的有限元背景网格，并通过inp文件导入Matlab内，通过背景网格上的一阶有限元插值函数计算出节点上的体积，高斯点处的坐标与权重信息；

步骤二、给定各节点与高斯点的支撑域半径，通过kdtree及rangesearch函数得到各点支撑域内的节点信息；

步骤三、输入模型初始裂纹面信息，将各点支撑域内穿过初始裂纹面的键预先断掉；这里的键不仅指支撑域节点与节点间的相互作用，也指支撑域内高斯点与节点间的相互作用；

步骤四，依据步骤二与步骤三得到的前处理信息，计算各点形函数及其导数，并依据各高斯点的权重，得到该处的节点刚度矩阵，并组装为全局刚度矩阵；

步骤五，将模型总的位移荷载或面力荷载划分为若干增量步，进行增量格式的Newmark时间迭代计算；

步骤六，在每一时步中，依据计算得到的应力，通过莫尔-库伦岩石强度准则，判断节点与节点及高斯点与节点间的键是否有新的断裂；如果有新键断裂，更新该点处刚度矩阵，并重新迭代计算，直至该时步计算平衡；

步骤七，迭代完成后，计算各点损伤值，得到岩石裂纹扩展的预测结果。

3.根据权利要求1所述的方法，其特征在于，通过RKPM形函数与高斯积分计算全局刚度矩阵；通过节点与节点间的键，高斯点与节点间的键的断裂情况表示裂纹。

4.根据权利要求1所述的方法，其特征在于，基于增量格式Newmark算法，具体方程如下：

ΔU为位移增量，

与

分别为等效刚度矩阵与等效增量载荷矩阵，具体为：

K，M，ΔQ分别为全局刚度矩阵，质量矩阵及增量载荷矩阵，

分别为t时刻的速度与加速度。上式各参数为：

α与δ分别为Newmark法中参数。当δ≥0.5，α≥0.25(0.5+δ)2时，算法是无条件稳定的。增量格式与全量格式间转化关系如下：

5.依据权利要求1所述的方法，其特征在于所述“对于每一时步计算出的应力值，判断各点间键的断裂状态，如果有断裂发生，更新该点处的刚度矩阵，并重新迭代计算，直至该时步内外力平衡”的隐式迭代流程，具体内容为：

在每一迭代步中，当有新的键断裂时，需要更新该键的状态，并重新计算与该键相关联点上的形函数[N]及对应的计算矩阵[B]；更新完成后，重新计算体系内的不平衡力R：

其中，[σ]为断键后的应力，通过[B]与U_t+Δt得到。新的有效刚度矩阵

为：

通过不平衡力R与有效刚度矩阵

得到新的位移增量[ΔU]：

t+Δt时刻新的位移为：

U_t+Δt＝U_t+ΔU+[ΔU]

重复上述过程，直至系统内不再出现新的断键，进入下一时步；

上述过程中，全局刚度矩阵K在裂纹扩展隐式迭代过程中的快速计算方法的公式为

式中，上标′old′代表上一迭代步，K_a代表当前待更新物质点x_a上的节点刚度矩阵，m代表需要更新的物质点数目。

Images