一种快速确定岩石类脆性材料离散元模型参数的方法与流程

本发明属于脆性材料离散元模型参数确定领域，尤其涉及到一种岩石类脆性材料离散元模型参数获取方法。

背景技术

岩石是一种天然地质体，具有不连续性、非均质性、各向异性和非线性等特点。随着高铁的建设、隧道的开挖、防护工程建设和页岩气的开采，对岩石破坏准则的研究正在成为一个重要的热点方向。

离散元法是一种可靠而有效的数值模拟方法，已被广泛应用于研究岩石类脆性材料的力学性能，其基本思想就是把整个介质视为一系列离散的颗粒组成，以力-位移定律和牛顿第二运动定律为依据，对问题进行力学分析。采用离散元方法研究岩石问题时，首先必须确定与材料宏观特性相匹配的模型细观参数。然而，岩石类脆性材料离散元模型中所含参数众多，且不易确定。准确获取岩石类脆性材料离散元模型中的关键参数，对于岩石类脆性材料的力学性能研究具有重要的理论意义和实际作用。

目前，岩石类脆性材料离散元模型参数一般是通过实验法来确定，即通过对批量材料试样进行多工况下的压缩实验、巴西劈裂实验以及三点弯曲实验等来确定。脆性材料离散元模型中有些参数可通过这些实验直接测取，但有些关键参数难以通过实验直接测取，常常需要进行大量反复的实验，通过数据拟合的方式来获取。发明专利申请201610058702.4公开了一种母岩及其颗粒料的二维离散元模型构建方法，该方法基于颗粒流软件(pfc^2d)二维模块环境构建母岩试样以及颗粒料试样的离散元模型但，但该模型仅提供一种二维模型的建立方法，无法广泛应用于岩石类脆性材料离散元模型参数的确定。此外，由于材料试样的制备、实验设备的限制，进行大量有效的岩石类脆性材料实验是成本较高、效率较低的。因此，开发其他快捷高效的方法来确定这些难以直接获取的脆性材料离散元模型参数是十分必要的。

技术实现要素：

针对以往确定岩石类脆性材料离散元模型参数需大量材料实验的问题，本发明的目的在于提出一种有效的岩石类脆性材料离散元模型参数确定方法，仅需少量实验数据就能快速有效地获取材料离散元模型中的关键参数。

本发明的技术方案是提供了一种岩石类脆性材料离散元模型参数确定方法，用于确定bpm模型中参数σc、ec和σc表示黏结法向强度，ec表示颗粒接触模量，表示颗粒法向与切向刚度比；其特征在，包括以下步骤：

步骤1：利用万能材料实验装置，分别对上表面与水平面夹角α为0°和1°的石灰岩试样进行准静态单轴压缩实验5次，并通过传感器、位移控制采集系统及数据处理装置获得石灰岩试样的应力应变数据；

步骤2：通过离散元方法建立与步骤1中水平面夹角α为0°的试样准静态单轴压缩实验相对应的材料离散元数值模型，其中采用bpm模型描述岩石类脆性材料的力学行为；

步骤3：基于步骤2所建立的材料离散元数值模型，对α＝0°的石灰岩材料在准静态单轴压缩工况下开展材料测量应力应变响应关于材料bpm模型参数σc、ec和的敏感性分析；

步骤4：基于步骤3所得的敏感性分析结果确定待识别的模型参数，并确定模型参数的取值范围；

步骤5：利用拉丁超立方试验方法并结合步骤4中所确定参数取值范围进行样本点采样，基于样本点采用支持向量机方法构建近似模型；

步骤6：构建目标函数，将利用步骤2数值求解得到的计算响应与步骤1实验测量响应相对比，根据实际问题及其解的要求建立求解模型，即目标函数；

步骤7：为获取步骤6中所建立目标函数的最小值，选择蚁群算法作为优化方法进行求解，求解时调用步骤5中所建立的近似模型；

步骤8：判断反求目标函数值是否大于目标函数的阈值ξ，如果大于目标函数的阈值ξ，则转入步骤9，如果不大于目标函数的阈值ξ，则转入步骤10；

步骤9：在σc、ec和取值范围内增加新的样本点，然后转入步骤5；

步骤10：此时bpm模型参数σc、ec和的值作为石灰岩材料离散元模型参数的最佳取值输出。

步骤11：将步骤10得到的bpm模型参数σc、ec和的值应用于1°石灰岩试样单轴压缩离散元模拟中，对比模拟出的应力应变曲线结果与步骤1所得的实验结果，进行验证。

进一步地，建立的α＝0°的石灰岩材料离散元数值模型为二维轴对称离散元模型。

进一步地，σc的取值范围为[30mpa,40mpa]，ec的取值范围为[20gpa,30gpa]，的取值范围为[1,3]。

进一步地，步骤8中，目标函数的阈值ξ为0.1。

本发明的有益效果在于：在综合利用石灰岩的单轴压缩试验和相应的离散元模型基础上，本发明提出了一种确定岩石类脆性材料离散元模型参数方法，本发明仅需结合少数物理实验就能快速获取石灰岩材料的离散元模型参数，有效减少材料相关物理实验的次数和降低成本，有利于离散元方法在脆性材料性能研究的应用拓展。

附图说明

图1是一种岩石类脆性材料离散元模型参数确定方法的流程图；

图2是几何模型示意图；

图3是所示为本发明一个实施例中所建立的α＝0°的石灰岩材料离散元数值模型；

图4是离散元数值计算与实验测量的应力应变曲线；

图5是该实施例中α＝1°的材料离散元数值模型；

图6是基于确定的离散元模型参数。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

在综合利用石灰岩材料的单轴压缩试验和相应的离散元模型基础上，本发明提出了一种确定岩石类脆性材料离散元模型参数的方法，其用于黏结颗粒模型(bondedparticlemodel，bpm)离散元模型的参数确定，其中该模型为：

上述公式中σ为材料的单轴抗压强度、e为材料弹性模量、υ为泊松比；和分别表示单轴抗压强度、弹性模量和泊松比的函数；σc、τc和ec分别表示黏结法向强度、黏结切向强度和颗粒接触模量；μ为颗粒摩擦因子、为颗粒法向与切向刚度比、为黏结法向与切向刚度比、n为离散元模型孔隙率。

上述模型中描述了材料宏观的单轴抗压强度、弹性模量和泊松比与离散元模型中颗粒特性参数的关系，其中根据本领域技术人员的研究认为黏结法向与切向刚度比近似等于颗粒法向与切向刚度比由于τc、μ和n对于不同岩石类脆性材料其值变化不大，故可根据本领域技术人员的经验确定。其他离散元模型参数σc、ec和通常需要基于大量的实验数据利用试错法来确定，但是要获取大量的材料实验数据周期长、成本很高，因此有效获取这些参数是较为困难的。因此，该实施例将结合材料单轴压缩实验、相应离散元模型和优化算法的基础上获取这些较难确定的材料离散元模型参数。

如图1所示，该实施例提供了一种岩石类脆性材料离散元模型参数确定方法，用于确定bpm模型中参数σc、ec和该方法包括具体实施步骤如下：

步骤1：利用万能材料实验装置，分别对上表面与水平面夹角α为0°和1°的石灰岩试样(几何模型示意图如图2所示)进行准静态单轴压缩实验5次，并通过高精度传感器、位移控制采集系统及数据处理获得石灰岩试样的应力应变数据，结果如表1和表2所示。

表1α＝0°试样的单轴压缩实验数据

表2α＝1°试样的单轴压缩实验数据

现有技术中公知的万能材料实验装置，技术非常成熟，采用机电一体化技术，主要由伺服驱动器、测力传感器、微处理器等构成。采用该实验装置，较易获取在准静态单轴压缩工况下材料的应力应变数据。分别对两种材料试样进行5次实验，将5次实验测量结果的平均值作为材料的最终实验测量数据。

在本发明一个实施例中，材料实验中的加载条件以轴向位移作为控制指标，以0.1mm/min的速度加载，边界约束条件为圆柱体下端固定，轴向加载直至试件破坏将作为步骤2中离散元数值模型的初始条件。实验测量得到的材料应力应变数据将作为确定岩石材料离散元模型参数的已知信息，其中：

α＝0°的材料试样测量数据将用于确定岩石材料离散元模型参数σc、ec和

α＝1°的材料试样测量数据将用于对确定的模型参数进行验证。

步骤2：通过离散元方法建立与步骤1中水平面夹角α为0°的试样准静态单轴压缩实验相对应的材料离散元数值模型；

图3所示为本发明一个实施例中所建立的α＝0°的石灰岩材料离散元数值模型，该模型为二维轴对称离散元模型，通过公知的离散元方法建立。

该是实施例中，采用pfc软件进行建模，建模过程中，材料样本的尺寸为长53mm×宽126mm的矩形，共由12016个颗粒填充，颗粒平均半径为0.4mm，颗粒间接触模型采用bpm模型，如图3所示。该离散元模型是一个基本的数值模拟，能够实现有效的数值计算，其中石灰岩材料bpm模型参数还不准确，需要对其进行有效确定。

步骤3：基于步骤2所建立的材料离散元数值模型，对α＝0°的岩石类脆性材料在准静态单轴压缩工况下开展材料测量应力应变响应关于材料bpm模型参数σc、ec和的敏感性分析；

步骤3.1、采用正交分析法对这些模型参数的敏感性进行分析，将测量响应的应力应变曲线量化为3个宏观指标：单轴抗压强度、弹性模量和泊松比，基于σc、ec和三个模型参数对单轴抗压强度、弹性模量和泊松比进行敏感性分析。

步骤3.2、采用l9(3⁴)正交试验表进行了正交分析，因素水平如表3所示；结合极差分析法分析了准静态下单轴压缩石灰岩材料离散元模拟(采用步骤2所建立的模型)的正交试验数据，结果如表4-表7所示。

采用极差值r对结果进行分析可发现：表5中σc所对应的r最大(59.22)，这说明σc对单轴抗压强度的影响最大；表6中ec所对应的r最大(38.76)，这说明ec对材料弹性模量的影响最大；表7中所对应的r最大(0.21)，这说明对泊松比的影响最大。

因此，这三个模型参数对材料应力应变响应的敏感性较强，即可以根据所提出的方法，结合实验测量获取的材料应力应变响应对模型参数σc、ec和进行确定。

表3正交试验因素水平表

表4正交分析试验结果值

表5单轴抗压强度极差分析表

表6弹性模量极差分析表

表7泊松比极差分析表

注：(1)ki^a＝包含a因素的“i”水平的目标和，a＝a,b,c；

(2)

(3)

步骤4：基于步骤3所得的敏感性分析结果确定待识别的模型参数，并确定模型参数的取值范围；

通过步骤3可以得到这三个参数对模拟应力应变响应的具体影响，结合步骤1所得到的实验结果以及表4所得到的模拟结果，可以确定三个参数的取值范围。

在本发明一实施例中，σc的取值范围为[30mpa,40mpa]，ec的取值范围为[20gpa,30gpa]，的取值范围为[1,3]。

步骤5：利用拉丁超立方试验方法并结合步骤4中所确定参数取值范围进行样本点采样，基于样本点采用支持向量机方法构建近似模型；

基于步骤4中选定模型参数的取值范围内，采用已公知的拉丁超立方试验方法进行样本采样，用于构建近似模型，本发明中的样本点数为20个。然后，将20组模型参数值代入离散元模型中进行准静态下单轴压缩数值模拟，获取了相应的20组石灰岩材料应力应变数据，即单轴抗压强度、弹性模量和泊松比，如表8所示，这些数据将用于构建近似模型。在构建近似模型过程中，采用已公知的支持向量机模型建立近似模型，该模型中所采用的函数为gaussian核函数。

表8拉丁超立方试验采样及数值计算结果

步骤6：构建目标函数，将利用步骤2数值求解得到的计算响应与步骤1实验测量响应相对比，根据实际问题及其解的要求建立求解模型，即目标函数；

在进行目标函数最小值求解时，需要调用大量数值求解所得的计算响应，如果直接采用离散元进行模拟得到计算响应，那将耗费大量的人力物力。调用步骤5中所建立的近似模型则可以节省大量的时间和精力。

该实施例中以同一应变下的实验测量应力与离散元数值计算应力差的平方和作为目标函数，其中n为响应采样个数，这样有效提高求解的计算精度，同时避免由于差异的大小引起的计算误差。其中离散元数值计算应力基于步骤5中所建立的近似模型中获取，n为20。

步骤7：为获取步骤6中所建立目标函数的最小值，选择蚁群算法作为优化方法进行求解，求解时需调用步骤5中所建立的近似模型；

在本步骤中，利用蚁群算法求解获得目标函数的最小值，输出最小的目标值和相对应的3个模型参数值。

本领域技术人员可以选择公有技术中的某蚁群算法作为优化算法，在蚁群算法中采用ant-cyclesystem模型，其中蚂蚁循环一周所释放的总信息素量为150、蚁群数量为20、信息素启发因子为1、信息素残留系数为0.5、期望启发因子为5，最大迭代步数为300。

步骤8：判断反求目标函数值是否大于目标函数的阈值ξ，如果大于目标函数的阈值ξ，则转入步骤9，如果不大于目标函数的阈值ξ，则转入步骤10。

在该实施例一个实施例中，设定目标函数的阈值ξ为0.1。

步骤9：在σc、ec和取值范围内增加新的样本点，然后转入步骤5。

在该实施例一个实施例中，在步骤4中选定的ec、和σc参数取值范围内，额外增加采样点，用于重新构建近似模型，进而转入步骤5。

步骤10：将此时bpm模型参数σc、ec和的值作为石灰岩材料离散元模型参数的最佳取值输出。

此时如图4所示，离散元数值计算与实验测量的应力应变曲线基本吻合，则确定的模型参数值为石灰岩离散元本构参数的最佳取值，参数确定结果如下表9所示：

表9石灰岩离散元模型参数的确定结果

为进一步验证确定结果的正确性与适用性，将表9所确定的模型参数值代入α＝1°的石灰岩材料中进行准静态单轴压缩离散元数值计算。

图5是该实施例中α＝1°的材料离散元数值模型。图6是基于确定的离散元模型参数，α＝1°的材料离散元计算结果与实验结果对比图，从图中可以看出离散元数值计算响应与实验测量曲线较为吻合，这验证了石灰岩脆性材料离散元模型参数确定结果的正确性。

以上对本发明一种确定岩石类脆性材料离散元模型参数方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。