三轴应力作用下颗粒离散元PFC矩张量声发射模拟

本文摘要（由AI生成）：

本文利用矩张量分析方法，深入探讨了岩石破坏产生的AE事件的震源机制。通过模拟颗粒间的接触，研究了不同岩石破坏类型，并定量分析了AE事件的震源行为。模拟试样的尺寸与实验试样相同，破坏模式、力-位移曲线和震级通过模拟得出。进一步，根据矩张量文件中的mag列数据，求得了b值。最后，利用T-k图和震源沙滩的MTLAB代码，绘制了相关图表，为理解岩石破坏的微裂纹产生机制提供了有力工具。

1. 引言

应力扰动，无论是自然的还是由人类活动引起的，都会引起地下地层的地震活动。当受到诱导应力时，完整的岩石或预先存在的裂缝可能会突然破坏。岩石灾难性破坏涉及以声发射（AE）、微震或诱发地震的形式发出不同震级地震活动。当进行实验室实验以分析断裂定位时，通过微裂纹的形成和相互作用可以产生声发射。

AE数据为研究岩石破坏的机理提供了宝贵的信息。AE事件的监测已用于研究导致岩石破坏的微裂纹聚结过程。由于完整岩石的破裂性质和断裂带的凹凸不平，因此对声发射数据进行了分析，以研究更大范围内的机制。在岩石力学试验（如三轴压缩试验）期间，可以使用岩石试样周围的一系列传感器记录实验室中的AE事件。AE震源的空间-时间分布能够表征完整岩石或具有预先存在裂缝的岩石的破坏。跟踪微裂纹位置和相互作用可以绘制剪切裂缝的成核、扩展和聚集图。此外，可以通过在不同围压下进行三轴压缩试验来研究不同带的形成（例如，压实或剪切）。

除了找到AE事件的位置外，对AE事件源机制的分析还提供了信息，以确定导致岩石破坏的破坏模式。近年来，大量研究工作者进行了不同的实验工作，以研究AE事件的震源机制。矩张量分析也被用作解释AE事件源机制的定量方法。

微尺度数值技术也可用于推断基本的岩石破坏机制。通过使用离散单元法（DEM），我们可以模拟完整和断裂岩石在压缩下的破坏，以研究断裂过程中涉及的主要机制。不同参数的影响，例如，颗粒形状和尺寸、晶粒互锁通过单个裂纹的能量释放诱导更多裂纹，颗粒的不规则性和角度对通过DEM模拟岩石力学行为和破坏模式进行了广泛研究。通过DEM研究晶粒位移和微裂纹发展，还可以模拟AE事件，并将其定位在微尺度数值样本中。

在这项工作中，我们通过确定岩石破坏产生的AE事件的矩张量、矩张量的分解来研究微裂纹产生的机制，以及解释每个AE事件的不同震源机制，重点对AE事件的震源行为进行定量分析。由于试验设计的复杂性，研究震源机制的实验工作可能限于特定场景（例如，岩石类型和围压）。使用机械地震耦合数值模型可以实现更广泛的场景。在这种情况下，我们可以更好地理解岩石破坏期间AE事件的震源机制。通过模拟颗粒和颗粒之间的接触，我们可以研究不同类型的岩石破坏。

2. 建模方法

AE模拟算法在PFC（particle Flow Code）中的平行粘结模型（BPM）中实现，PFC是Itasca开发的显式离散元方法（DEM）的一个分支。这里，假设粒子是刚性的（不可变形的），并通过模拟粒子之间相互作用行为的接触模型结合。假设每一个粘结断裂都是一个微裂纹，并且每个微裂纹都会在测试期间导致AE活动。通过在测试期间监测所有粘结断裂，可以收集AE信息，并计算矩张量（MT），以分析断裂类型和破坏机制。

2.1. 平行粘结模型

平行粘结模型（BPM）作为一种离散元方法，最初由Cundall引入，并由Cundal和Strack开发。Potyondy等人详细描述了BPM的原理。该模型模拟了由二维或三维粘结的非均匀尺寸圆形或球形颗粒表示的固体材料的运动和相互作用。模拟材料的特性通常由颗粒和粘结键的刚度和强度微参数决定。商业软件PFC中提供了两种类型的BPM，即接触粘结模型和平行粘结模型。接触粘结作用于接触点，只能传递法向力和剪切力。平行粘结作用于两个颗粒之间的横截面，可以传递力和力矩，这是岩石或类岩石材料的更现实的粘结模型。

2.2. AE模拟和矩张量分解

当粘结键断裂时，源粒子将移动，相邻触点将发生一些变形。因此，由于裂纹的形成，周围接触处的力将发生变化。矩张量的分量可以通过将接触点处的力变化乘以接触点到裂纹位置的距离来计算。值得注意的是，AE事件可以由多个微裂纹组成。在空间和时间上紧密发生的微裂纹被视为一个AE事件，事件质心被假设为事件的几何中心（图1（a））。力矩张量可以通过在源边界周围积分来计算，对于离散元方法，积分由以下公式计算：

其中ΔF_i是接触力变化的第i个分量，R_j是接触点和事件形心之间的距离。由于假设力矩张量是对称的，因此力矩张量的非对角分量将被假设为相等，并且将通过等式（1）计算的两个值进行平均来计算。此外，在AE事件持续期间的每个时间步计算矩张量，为每个AE事件存储的单个矩张量是在最大标量力矩的时间步计算的结果。

图1 示例AE事件及矩张量示意图

矩张量可以被认为是源的等效体力，矩张量表示将矩张量矩阵的主值描述为两组向量，其方向和长度分别表示方向和大小（图1（b））。

将公式（1）计算得到的矩张量分解为球张量与偏张量，则声发射破裂类型可根据球张量占比P_iso进行判定

2.3. 模型描述

数值试样的尺寸与实验试样的尺寸相同（如图2所示），直径为50 mm，高度100 mm。考虑到计算时间，颗粒尺寸遵循1 mm至1.4 mm的均匀分布，样品中总共产生15178个颗粒。

图2 数值模型

3. 模拟直接得到结果

模拟直接得出的破坏模式、力-位移曲线和震级分别如图3所示。

(a) 破坏模式

(b) 力-位移曲线

(c) 震级分布

（d）裂纹类型曲线

（e）体积应变

图3 模拟直接得到的结果

模拟得到的dat文件如图4所示。

(a) AE

（b）裂纹文件

（c）矩张量文件

图5 模拟得到的dat文件

3. 后处理得到结果

应用矩张量分析程序对试样在加载作用下的声发射信号进行分析，图6显示了声发射空间分布，利用图5（c）中ios和公式（2）可将破坏模式分为拉伸、剪切和混合。图6中青色球体表示剪切事件，蓝色球体表示张拉事件，红色球体表示混合事件，球体大小表示能力和震级的相对大小。此图为去掉仅有1条裂纹的声发射事件的定位图。

图6 不同类型AE事件空间分布图

在模拟过程中，通过矩张量文件中aenum记录了AE事件，进而可以做出应力应变、AE事件数和累计AE事件数之间的关系，如图7所示。

图7 AE事件和AE累积事件的演变过程

4. 通过AE b值评估裂缝网络的复杂性

断裂系统的产生伴随着现场的微震或实验室的AE，这有助于描述断裂分布。同时，我们注意到，Gutenberg和Richter提出的b值不仅可以用于描述AE震级的分形特征，还可以用于描述声发射源的空间分形特征。较大的b值意味着更复杂的裂缝网络，而较小的b值则意味着更简单的裂缝网络。描述地震震级和频率之间关系的古腾堡和里希特公式如下：

其中M为地震震级；N是震级大于M的地震的累积次数；a和b是两个统计常数。根据矩张量文件中的mag列数据可以得出下图，进而求出b值。

图8 AE事件数量和AE震级的关系（直方图显示频率-震级分布，散点显示大于某个震级的AE事件数量）。

5. T-k绘图法

Hudson（1989）提出了矩张量的两个参数T和k，以分析震源破裂类型。参数T用于测量源的偏差部分，其范围从位于−1的纯正CLVD到位于+1的纯负CLVD，与位于原点的纯DC交叉。参数k用于测量各向同性部分，其范围从位于+1的均匀爆炸类型到位于-1的均匀内爆类型。

根据矩张量文件里M11、M12 、M22、M31、M32和M33，利用T-k图MTLAB代码（包含在附件中）可做出图9。

图9 T-k图

6. 震源沙滩求

震源机制球是包围地震震源A的一个假想的球壳（震源球）的下半部分在水平面的投影。断层面与球壳相交处表现为一条线。断层破裂时应力场取向控制着断层面上滑动的方向，沙滩球也描述了应力的取向，灰色（深色、红色）区域包含着反映最小应力方向的张力轴（T），而白色域包含着反映最大应力方向的压力轴（P）。

根据矩张量文件里M11、M12 、M22、M31、M32和M33，利用震源沙滩求MTLAB代码（包含在附件中）可做出图10。此处挑选两个较大震级事件进行作图。

图10震源沙滩求

主要图形汇总