本发明涉及声发射模拟领域，特别是涉及一种基于离散单元法的声发射模拟方法及系统。
背景技术：
理解岩石的破裂过程对很多工程具有重要意义。然而，很少有方法能直观看到这一过程，研究者们常常借助其他手段来进行研究。这其中最主要的技术手段就是声发射，其主要优点在于能实现岩石破裂过程中的无损监测和连续监测。
室内对岩石样品的声发射试验能得到可靠的声发射结果，但它往往有很多限制条件。如isrm(internationalsocietyforrockmechanics,国际岩石力学学会)对岩石样品的尺寸和端面不平整度做了严格规定，而野外采样难以保证岩石样品的完整性；岩石样品从野外运输至室内不确定性因素很多，要保证岩石样品的完好需要大量的人力物力财力；此外，该实验要求周围低噪音的实验环境。实验过程中也难以保证对单一因素的控制，降低了实验结果的可对比性。
数值模拟可以克服实验技术上的不足，是对岩石破裂过程研究的有力手段。离散单元法在模拟岩石宏观破裂、损伤演化方面具有显著优势，是模拟岩石破裂的主流方法之一。在离散单元法中，岩石材料被表示为离散化的颗粒的组合体，他们之间通过粘结键相连接。这些颗粒的运动和颗粒间的相互作用实现了受力岩石的破裂过程。当外界作用力超过了粘结键的强度时，这些粘结键将发生断裂，对应了实际岩石样品的裂纹扩展。在模拟声发射这一问题上，基于离散单元法的颗粒流程序pfc(particleflowcode)，hazzard和young在2000年率先提出了一种方法，是通过监测粘结键断裂时周围的颗粒动能，并将其用于直接量化震源辐射的声发射能量。通过在空间和时间上对多个粘结键的断裂进行聚类分析，最终得到了花岗岩的声发射特征。2002年，hazzard和young引入了粘结键断裂时变化的力矩张量这一参数，改进了上述方法，成功模拟了岩石破裂过程中的声发射实验。这一方法巧妙利用了离散单元法的特点，能够较好地模拟声发射，因此被广泛应用至今。然而许多学者指出该方法存在许多问题：①该模拟技术并非根据真实声发射实验的原理。真实的声发射实验是通过在岩石试样表面布置监测点，在监测点处安装压电传感器，直接监测试样表面的法向运动分量，并对得到的运动波形进行分析，最终得到声发射特征参数。而该模拟技术是通过监测粘结键断裂时的矩张量变化。因此二者的原理存在明显差异；②该模拟技术只能监测新裂纹扩展，不能监测裂纹扩展后的重新活动，如摩擦导致的声发射不能被监测到。由于模拟技术直接监测的矩张量变化，一个或多个临近的粘结键断裂时视为一个新的声发射事件，并没有考虑裂纹的重新活动，因此很多声发射活动不能被记录下来，不能真实反映声发射；③计算矩张量的变化将占用较多的计算资源，延缓模拟速度，导致计算效率偏低。
针对以上矩张量法模拟声发射技术的原理不符、真实性差、计算资源占用较多的问题，我们提出了一种基于离散单元法的声发射模拟方法及系统。
技术实现要素：
本发明的目的是提供一种基于离散单元法的声发射模拟方法及系统，解决现有技术中模拟声发射技术的原理不符、真实性差、计算资源占用较多的问题。
为实现上述目的，本发明提供了如下方案：
一种基于离散单元法的声发射模拟方法，所述模拟方法包括：
在计算模型表面布置监测点；
监测所述监测点的速度波形；
对所述速度波形进行分析，得到声发射数、能量以及b值。其中，声发射数被定义为超过门槛值并导致系统通道积累的信号数据，它与裂纹增长成正比；能量被定义为信号波形的平方，它与破裂的震级成正比，第i个数据的能量e(i)计算公式如下：
其中vx指速度的x方向分量，vy指速度的y方向分量；b值代表声发射的规模，被广泛用于表征裂纹扩展的尺度，常用gutenberg–richter公式进行线性拟合得出：
logn＝a-b×m
其中n指的是震级大于等于m的数据数量，a和b为常数，m为震级，常用以下公式得出：
可选的，对所述速度波形进行分析，得到声发射数具体包括：
a1：导入监测点的速度波形；
a2：依次对每个采样数据进行分析；
a3：判断第i个数据是否小于或等于采样总数据数量；
a4：若是，则结束；
a5：若否，则执行下一步；
a6：判断第i个数据是否大于或等于预设的门槛值；
a7：若否，则i＝i 1，返回步骤a2，对下一个数据进行判断；
a8：若是，则执行下一步骤；
a9：判断一次声发射事件是否结束；
a10：若否，则强行切断，该声发射事件结束，执行下一步骤；
a11：若是，则记录本次声发射事件的特征参数，所述特征参数包括到达时间、幅值、信号持续时间、振铃计数；
a12：对下一个声发射信号进行分析，直到所有采样数据被分析完毕。
可选的，对所述速度波形进行分析，得到能量具体包括：
b1：判断第i个数据是否小于或等于采样总数据数量；
b2：若否，则结束；
b3：若是，则计算e(i)＝[vx2(i) vy2(i)]；
b4：计算
b5：i＝i 1，返回步骤b1，直到所有采样数据分析完毕。
可选的，对所述速度波形进行分析，得到b值具体包括：
c1：输入震级范围[mmin,mmax]，并输入m＝mmin；
c2：判断第i个数据是否小于或等于采样总数据数量；
c3：若否，则结束；
c4：若是，则执行下一步骤；
c5：获取m(1)～m(i)中大于或等于mmin的数据个数，计为n1，获取m(1)～m(i)中大于或等于mmin 0.1的数据个数，计为n2，直到获取完m(1)～m(i)中大于或等于mmax的数据个数，记为nm；
c6：得到m组(mmin～max)～n(1～m)数据，根据公式logn＝a-bm进行拟合得b(i)值，i＝i 1，返回步骤c2。
本发明另外提供一种基于离散单元法的声发射模拟系统，所述模拟系统包括：
监测点布置模块，用于在计算模型表面布置监测点；
速度波形监测模块，用于监测所述监测点的速度波形；
速度波形分析模块，用于对所述速度波形进行分析，得到声发射数、能量以及b值。
可选的，所述速度波形分析模块包括：声发射数分析单元、能量分析单元以及b值分析单元。
可选的，所述声发射数分析单元具体包括以下步骤：
a1：导入监测点的速度波形；
a2：依次对每个采样数据进行分析；
a3：判断第i个数据是否小于或等于采样总数据数量；
a4：若是，则结束；
a5：若否，则执行下一步；
a6：判断第i个数据是否大于或等于预设的门槛值；
a7：若否，则i＝i 1，返回步骤a2，对下一个数据进行判断；
a8：若是，则执行下一步骤；
a9：判断一次声发射事件是否结束；
a10：若否，则强行切断，该声发射事件结束，执行下一步骤；
a11：若是，则记录声发射特征参数与到达时间、幅值、信号持续时间、振铃计数；
a12：对下一个声发射信号进行分析，直到所有采样数据被分析完毕。
可选的，所述能量分析单元具体包括以下步骤：
b1：判断第i个数据是否小于或等于采样总数据数量；
b2：若否，则结束；
b3：若是，则计算e(i)＝[vx2(i) vy2(i)]；其中，e(i)表示第i个数据的能量，vx指速度的x方向分量，vy指速度的y方向分量；
b4：计算
b5：i＝i 1，返回步骤b1，直到所有采样数据分析完毕。
可选的，所述b值分析单元具体包括以下步骤：
c1：输入震级范围[mmin,mmax]，并输入m＝mmin；
c2：判断第i个数据是否小于或等于采样总数据数量；
c3：若否，则结束；
c4：若是，则执行下一步骤；
c5：获取m(1)～m(i)中大于或等于mmin的数据个数，计为n1，获取m(1)～m(i)中大于或等于mmin 0.1的数据个数，计为n2，直到获取完m(1)～m(i)中大于或等于mmax的数据个数，记为nm；
c6：得到m组(mmin～max)～n(1～m)数据，根据公式logn＝a-bm进行拟合得b(i)值，i＝i 1，返回步骤c2。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
1.原理与实际声发射实验一致。传统的矩张量法虽然能够通过矩张量的变化模拟声发射，但与实际声发射实验原理相差较大。实际的声发射实验是通过在岩石试样表面布置监测点，在监测点处安装压电传感器，直接监测试样表面的法向运动分量，并对得到的运动波形进行分析，进而得到声发射特征参数。本方法通过在模型表面布置监测点，直接监测模型在模拟过程中x方向的速度分量，得到速度波形图，并依据速度波形图，编写了声发射特征参数的计算单元，得到了声发射特征参数。故本方法解决了利用真实声发射实验原理来模拟的难题；
2.真实性更高。在与实际声发射实验原理相一致的基础上，该方法相比于传统的矩张量法，真实性更高。算例表明，通过该方法得到的声发射能量、事件数、b值与矩张量法接近，并更接近于真实的声发射实验。
3.计算速度快。矩张量法需要计算每步粘结键的力与位移的变化，进而构建矩阵进行运算。而新方法只需要监测监测点的速度分量即可。相同条件下，本方法相比于矩张量法运算速度提高了4.9倍，大大节约了运算时间和内存占用。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例基于离散单元法的声发射模拟方法流程图；
图2为本发明实施例获取声发射数的方法流程图；
图3为本发明实施例获取声发射能量的方法流程图；
图4为本发明实施例获取声发射b值的方法流程图；
图5为本发明实施例应力应变曲线的拟合示意图；
图6为本发明实施例监测点布置示意图；
图7为本发明实施例声发射数对比图；
图8为本发明实施例声发射能量对比图；
图9为本发明实施例b值对比图；
图10为本发明实施例基于离散单元法的声发射模拟系统结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种基于离散单元法的声发射模拟方法及系统，解决现有技术中模拟声发射技术的原理不符、真实性差、计算资源占用较多的问题。
本发明所述的方法及系统适用于离散单元法(包含pfc、udec、3dec等多个软件)，只是以pfc软件为例进行分析说明。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明所述方法完全从真实声发射实验的原理出发，通过直接在pfc中数值计算模型的表面布置监测点，并监测这些点在x方向的速度波形，通过对波形的分析，得到模型破裂过程中的声发射信息。具体技术方案如下：
1.根据室内物理模型实验的结果，得到岩石样品的单轴抗压强度、杨氏模量等参数。
2.用pfc对该物理模型实验进行再现，得到满意的模拟效果。对该岩石样品进行模拟时，由于无法得到岩石矿物颗粒间的胶结力学参数，故pfc中的微观参数只能通过试错法确定，即通过不断调整微观参数直至模拟的宏观力学属性与实际岩石样品的宏观力学属性相一致。
3.在pfc中的计算模型表面布置监测点。真实声发射实验中的声发射探头布置位置即为声发射的波形采样点，也叫监测点。为尽量贴近真实的声发射实验，在pfc中可将监测点布置位置与实验中的声发射探头布置位置相一致，即都位于模型中部。
4.得到监测点的速度波形。真实声发射实验是通过声发射探头监测采样点的运动波形。在pfc中，可直接通过内置fish语言的“ballhistoryvelocity”命令直接监测预设监测点的运动波形。
5.对速度波形进行分析，得到声发射数、能量、b值等声发射参数。具体步骤如下：
图1为本发明实施例基于离散单元法的声发射模拟方法流程图，如图1所示，所述方法包括：
步骤101：在计算模型表面布置监测点。
步骤102：监测所述监测点的速度波形。
步骤103：对所述速度波形进行分析，得到声发射数、能量以及b值。
具体的，步骤103中：
①发射数的获取。在pfc模拟过程中，系统默认每10步得到一个数据，即得到的速度波形其采样时间为10步。分析开始前，需预设门槛值、结束定义时间和声发射最大持续时间。导入监测点处的速度波形，从头至尾依次对每个采样数据进行分析。当第i个数据大于或等于预设的门槛值时，标志着一次声发射事件的开始。如果小于预设门槛值则对下一个数据进行判断。下一步判断该声发射事件是否结束，即在最大持续时间内是否出现采样数据连续低于门槛值达到结束定义时间的情况。如果未结束则认定该声发射信号属于连续型，该信号在达到最大持续时间时将被强行截断，该声发射事件结束。此时，记录本次声发射事件的特征参数，如到达时间，幅值，持续时间，振铃计数等，并准备对下一个声发射信号进行分析，直至所有采样数据被分析完毕。
具体流程如图2，详细步骤如下：
a1：导入监测点的速度波形；
a2：依次对每个采样数据进行分析；
a3：判断第i个数据是否小于或等于采样总数据数量；
a4：若是，则结束；
a5：若否，则执行下一步；
a6：判断第i个数据是否大于或等于预设的门槛值；
a7：若否，则i＝i 1，返回步骤a2，对下一个数据进行判断；
a8：若是，则执行下一步骤；
a9：判断一次声发射事件是否结束；
a10：若否，则强行切断，该声发射事件结束，执行下一步骤；
a11：若是，则记录声发射特征参数与到达时间、幅值、信号持续时间、振铃计数；
a12：对下一个声发射信号进行分析，直到所有采样数据被分析完毕。
②能量的获取。实际声发射实验中声发射能量e常以振幅平方计数。对于得到的速度波形平方即可得到采样点处的能量。运用经验公式
b1：判断第i个数据是否小于或等于采样总数据数量；
b2：若否，则结束；
b3：若是，则计算e(i)＝[vx2(i) vy2(i)]；
b4：计算
b5：i＝i 1，返回步骤b1，直到所有采样数据分析完毕。
③b值的获取。首先预设合理的震级范围，构建合理的震级梯度。其次判断1～i个采样数据中大于等于各震级m的个数n，并用公式logn＝a-bm进行线性拟合，其中a，b为常数，进而得到第i个采样数据对应的b值。具体流程如图4，详细步骤如下：
c1：输入震级范围[mmin,mmax]，并输入m＝mmin；
c2：判断第i个数据是否小于或等于采样总数据数量；
c3：若否，则结束；
c4：若是，则执行下一步骤；
c5：获取m(1)～m(i)中大于或等于mmin的数据个数，计为n1，获取m(1)～m(i)中大于或等于mmin 0.1的数据个数，计为n2，直到获取完m(1)～m(i)中大于或等于mmax的数据个数，记为nm；
具体的，对于1～i个采样点，预先设定一个震级范围，比如[1,2,3,4…]，判断m(1)～m(i)有多少个大于等于1，记录为n(1)，再判断有多少大于等于2，记录为n(2)，依次类推，得到震级范围[1,2,3,4…]所对应的[n(1)，n(2)，n(3)，n(4)]，再根据公式进行拟合，得到b(i)值。
c6：得到m组(mmin～max)～n(1～m)数据，根据公式logn＝a-bm进行拟合得b(i)值，i＝i 1，返回步骤c2。
实施例：
首先进行了页岩在单轴压缩下的声发射试验，得到岩石样品的单轴抗压强度、杨氏模量等宏观力学参数；
通过离散元软件pfc对该实验进行了再现。数值计算模型尺寸与岩石样品一致，所得应力应变曲线如图5所示，可得所建数值计算模型能较好再现实验过程；
在模型表面的轴向和侧向分别均匀布置三个监测点，间距分别为20mm和25mm，并依次标号为1-9。如图6所示。
监测9个监测点的速度波形；
对速度波形进行分析，得到声发射计数、能量和b值，以2号监测点为例进行对比分析；
用传统的矩张量法对模拟过程中的声发射进行了模拟，得到声发射计数、能量和b值；
将实验、矩张量法、本方法所得声发射结果进行了对比，如图7、图8、图9所示。
图10为本发明实施例基于离散单元法的声发射模拟系统结构示意图，如图10所示，所述系统包括：
监测点布置模块201，用于在计算模型表面布置监测点。
速度波形监测模块202，用于监测所述监测点的速度波形。
速度波形分析模块203，用于对所述速度波形进行分析，得到声发射数、能量以及b值。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
1.一种基于离散单元法的声发射模拟方法，其特征在于，所述模拟方法包括：
在计算模型表面布置监测点；
监测所述监测点的速度波形；
对所述速度波形进行分析，得到声发射数、能量以及b值。
2.根据权利要求1所述的基于离散单元法的声发射模拟方法，其特征在于，对所述速度波形进行分析，得到声发射数具体包括：
a1：导入监测点的速度波形；
a2：依次对每个采样数据进行分析；
a3：判断第i个数据是否小于或等于采样总数据数量；
a4：若是，则结束；
a5：若否，则执行下一步；
a6：判断第i个数据是否大于或等于预设的门槛值；
a7：若否，则i＝i 1，返回步骤a2，对下一个数据进行判断；
a8：若是，则执行下一步骤；
a9：判断一次声发射事件是否结束；
a10：若否，则强行切断，该声发射事件结束，执行下一步骤；
a11：若是，则记录本次声发射事件的特征参数，所述特征参数包括到达时间、幅值、信号持续时间、振铃计数；
a12：对下一个声发射信号进行分析，直到所有采样数据被分析完毕。
3.根据权利要求1所述的基于离散单元法的声发射模拟方法，其特征在于，对所述速度波形进行分析，得到能量具体包括：
b1：判断第i个数据是否小于或等于采样总数据数量；
b2：若否，则结束；
b3：若是，则计算e(i)＝[vx2(i) vy2(i)]；其中，e(i)表示第i个数据的能量，vx指速度的x方向分量，vy指速度的y方向分量；
b4：计算
b5：i＝i 1，返回步骤b1，直到所有采样数据分析完毕。
4.根据权利要求1所述的基于离散单元法的声发射模拟方法，其特征在于，对所述速度波形进行分析，得到b值具体包括：
c1：输入震级范围[mmin,mmax]，并输入m＝mmin；
c2：判断第i个数据是否小于或等于采样总数据数量；
c3：若否，则结束；
c4：若是，则执行下一步骤；
c5：获取m(1)～m(i)中大于或等于mmin的数据个数，计为n1，获取m(1)～m(i)中大于或等于mmin 0.1的数据个数，计为n2，直到获取完m(1)～m(i)中大于或等于mmax的数据个数，记为nm；
c6：得到m组(mmin～max)～n(1～m)数据，根据公式logn＝a-bm进行拟合得b(i)值，i＝i 1，返回步骤c2。
5.一种基于离散单元法的声发射模拟系统，其特征在于，所述模拟系统包括：
监测点布置模块，用于在计算模型表面布置监测点；
速度波形监测模块，用于监测所述监测点的速度波形；
速度波形分析模块，用于对所述速度波形进行分析，得到声发射数、能量以及b值。
6.根据权利要求5所述的基于离散单元法的声发射模拟系统，其特征在于，所述速度波形分析模块包括：声发射数分析单元、能量分析单元以及b值分析单元。
7.根据权利要求6所述的基于离散单元法的声发射模拟系统，其特征在于，所述声发射数分析单元具体包括以下步骤：
a1：导入监测点的速度波形；
a2：依次对每个采样数据进行分析；
a3：判断第i个数据是否小于或等于采样总数据数量；
a4：若是，则结束；
a5：若否，则执行下一步；
a6：判断第i个数据是否大于或等于预设的门槛值；
a7：若否，则i＝i 1，返回步骤a2，对下一个数据进行判断；
a8：若是，则执行下一步骤；
a9：判断一次声发射事件是否结束；
a10：若否，则强行切断，该声发射事件结束，执行下一步骤；
a11：若是，则记录声发射特征参数与到达时间、幅值、信号持续时间、振铃计数；
a12：对下一个声发射信号进行分析，直到所有采样数据被分析完毕。
8.根据权利要求6所述的基于离散单元法的声发射模拟系统，其特征在于，所述能量分析单元具体包括以下步骤：
b1：判断第i个数据是否小于或等于采样总数据数量；
b2：若否，则结束；
b3：若是，则计算e(i)＝[vx2(i) vy2(i)]；
b4：计算
b5：i＝i 1，返回步骤b1，直到所有采样数据分析完毕。
9.根据权利要求6所述的基于离散单元法的声发射模拟系统，其特征在于，所述b值分析单元具体包括以下步骤：
c1：输入震级范围[mmin,mmax]，并输入m＝mmin；
c2：判断第i个数据是否小于或等于采样总数据数量；
c3：若否，则结束；
c4：若是，则执行下一步骤；
c5：获取m(1)～m(i)中大于或等于mmin的数据个数，计为n1，获取m1～mi中大于或等于mmin 0.1的数据个数，计为n2，直到获取完m(1)～m(i)中大于或等于mmax的数据个数，记为nm；
c6：得到m组(mmin～max)～n(1～m)数据，根据公式logn＝a-bm进行拟合得b(i)值，i＝i 1，返回步骤c2。
技术总结