CN115859714B - 一种基于fem-dem联合仿真的岩石爆破全过程模拟方法 - Google Patents

Abstract

本发明公开了一种基于FEM‑DEM联合仿真的岩石爆破全过程模拟方法，涉及岩土工程爆破数值模拟技术领域，包括：基于有限元法和实际各段别炮孔爆破参数和装药结构，建立有限元模型；求解有限元模型，模拟炸药爆轰过程，获取各段别炮孔的孔壁压力荷载时程曲线和粉碎区边界压力荷载时程曲线，以及岩石爆破损伤分布和振动情况；基于离散元法和实际工程形状、尺寸和岩体类型，建立离散元模型；将各段别炮孔粉碎区边界荷载时程曲线分别施加至离散元模型各段别粉碎区边界，对离散元模型进行爆破动力响应分析，获取岩石爆破抛掷形态和爆破效果；综合来看，本发明可以提升数值模拟的预测精度和提高计算效率。

Description

一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法

技术领域

本发明涉及岩土工程爆破数值模拟技术领域，具体为一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法。

背景技术

钻孔爆破因其经济高效等优点，目前仍是我国水电施工、矿山开采、交通运输工程中最为常用的施工方法。在实施爆破之前，工程人员往往采用数值模拟方法初步预测爆破效果，可见，数值模拟结果在一定程度上可以指导现场的爆破方案。因此，有必要寻求一种较准确模拟岩石爆破的方法。

在岩石爆破中，目前主要的数值模拟方法有3类：一种是连续介质方法，如有限元和有限差分等方法；一种是非连续介质方法，如离散元和DDA等方法；还有一种是连续-非连续方法，如SPH-FEM等方法。

其中，连续介质方法如有限元法(FEM)是一种求解偏微分方程边值问题近似解的数值技术，求解时对整个问题域进行分解，每个子区域都成为简单的部分。以LS-DYNA软件在爆炸领域应用为例，在软件中建立炸药单元，定义炸药参数和状态方程，并通过流固耦合算法实现炸药-岩石的荷载传递，可以很好地模拟炸药起爆的爆破效果，如孔壁荷载、损伤和振动响应等。由于炸药爆炸是一种复杂的瞬态物理、化学反应，若想得到较准确的计算结果，对炮孔周围网格质量要求较高，进而造成模型网格较多，其计算时间较长。另外，FEM软件难以模拟岩石爆破的破碎、抛掷以及爆堆形态。

非连续介质方法如离散元法(DEM)是20世纪70年代由Cundall首先提出来的，起源于分子动力学，是为研究岩体等非连续介质的力学行为而发展起来的一种数值方法。在岩石爆破模拟方面，DEM软件常常被应用于的破碎、抛掷以及爆堆形态，其大致步骤主要是：首先，将爆源荷载简化为函数；再基于理论或经验公式，计算出粉碎区边界上的荷载；最后，将其施加至DEM模型粉碎区边界上。可知，DEM软件无法模拟炸药的爆轰过程，而往往通过在模型中施加简化荷载进行计算，其计算的准确性和可靠性难以保证。

现有的连续-非连续耦合方法中，对于模型中网格和接触面的设置均较为复杂，其计算时间较长，效率偏低。

综上所述，现有技术中对于岩石爆破的模拟，主要通过连续介质方法、非连续介质方法和连续-非连续耦合方法，但是三种方法均存在一定的缺点，导致数值模拟预测精度或计算效率较低，难以投入到实际的工程应用之中。

发明内容

为了克服现有数值模拟技术预测精度或计算效率低以及模拟方法的局限性，本发明提供了一种基于FEM-DEM联合仿真的岩石爆破全过程数值方法；该方法综合了FEM软件(本发明采用LS-DYNA软件)与DEM软件(本发明采用3DEC软件)各自的优点，改变了以往仅使用FEM或DEM软件模拟岩石爆破的局限性，结合了FEM软件模拟爆破动力响应和DEM软件模拟爆破抛掷情况及爆堆形态的各自优势，提出一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，不仅能够模拟岩石爆破全过程，而且提升了数值模拟的预测精度，还进一步提高了计算效率。总体上，本发明具有较高的实用性和信息化水平，能够应用于工程实际中。

本发明提出的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，包括以下步骤：

基于有限元法和实际各段别炮孔爆破参数和装药结构，建立有限元模型；

设置有限元模型的无反射边界和对称边界，并确定有限元模型中的炸药和空气参数及状态方程、岩石的本构模型和物理力学参数；

求解有限元模型，模拟炸药爆轰过程，获取各段别炮孔的孔壁压力荷载时程曲线和粉碎区边界压力荷载时程曲线，以及岩石爆破损伤分布和振动情况；

基于离散元法和实际工程形状、尺寸和岩体类型，建立离散元模型；

设置离散元模型中的围岩本构模型及物理力学参数和炮孔布置参数；

将有限元模型输出的各段别炮孔粉碎区边界荷载时程曲线分别施加至离散元模型各段别粉碎区边界，设置离散元模型的粘滞边界和固定边界，对离散元模型进行爆破动力响应分析，获取岩石爆破抛掷形态和爆破效果。

进一步的，所述有限元模型的建立是利用LS-DYNA软件；

所述离散元模型的建立是利用3DEC软件。

进一步的，通过现场或实验室获取实际各段别炮孔爆破参数、装药结构、实际工程形状、尺寸和岩体类型。

进一步的，岩土物理力学参数包括：密度ρ、弹性模量e、剪切模量G、泊松比ν、粘聚力c、内摩擦角μ、抗拉强度T_e、节理的切向和法向刚度比K_s/K_n。

进一步的，所述获取各段别炮孔的孔壁压力荷载时程曲线和粉碎区边界压力荷载时程曲线，具体包括：

求解有限元模型，得到各段别炮孔的爆炸荷载压力时程曲线和粉碎区爆炸荷载压力时程曲线，将孔壁压力荷载峰值和粉碎区压力荷载峰值与理论计算值进行对较，当孔壁压力荷载峰值和粉碎区压力荷载峰值与理论计算值的相对误差小于10％，则判定爆炸荷载压力曲线的荷载峰值符合理论公式；

其中，炮孔的孔壁荷载峰值理论计算公式为：

式中：P₀为炮孔的爆炸荷载压力，ρ₀为炸药密度，D为炸药爆轰速度，x为爆炸压力与起爆点的距离，l为药卷长度；

粉碎区外边界上径向应力峰值与距离的关系可表示为：

式中：P_b为粉碎区荷载压力，P₀为炮孔壁的爆炸荷载压力，r_c为粉碎区半径，r₀为炮孔半径，α为衰减指数。

进一步的，在对离散元模型进行爆破动力响应分析前，基于哈里斯爆破理论对离散元模型不同部位的网格尺寸进行预处理，使岩石爆破结果更符合实际；

其中，根据哈里斯爆破理论，半径为r的球形装药上，距球心为R的任意一点上的切应变为：

式中：β为应变波吸收系数；K为炮孔壁上的应变；

炮孔壁上的应变K计算公式为：

式中：P₀为爆压；ρ为岩石的密度；v_p为岩石纵波声速；μ为泊松比；γ为绝热指数。

在压缩应力波作用下岩石作径向位移，由径向位移衍生的切向应变值ε超过岩石的动态极限抗拉应变值T时，岩石形成径向裂缝；距炮孔距离为R处的径向裂缝的条数为：

设定破碎岩块的最大尺寸为两条最靠近的径向裂缝之间的距离，岩块的长L为：

因此，距离半径为b的炮孔距离R处的岩块的长度为：

基于上式得距爆心距离与预测爆破块度关系，选取不同块度的大小，计算出爆心距离并划分块度大小对应的区域范围，将模型按照距离炮孔中心的远近从小到大进行块度划分。

进一步的，所述设置有限元模型的无反射边界、对称边界、炸药和空气参数及状态方程、岩石的本构模型和物理力学参数，具体包括：

根据现场实际情况设置有限元模型的尺寸和参数，并对有限元模型网格进行精细划分，从炮孔近区开始由近到远采用渐变网格，网格最大控制在5cm以内。

进一步的，所述对离散元模型进行爆破动力响应分析，获取岩石爆破抛掷形态和爆破效果，具体包括：

通过设置对离散元模型的计算时间，获取岩石爆破不同时间的抛掷形态和爆破效果。

与现有技术相比，本发明提供的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其有益效果是：

本方法避免了单独使用有限元法无法模拟爆破过程中破碎块体的抛掷和爆破效果，以及单独使用离散元法无法准确获得爆炸荷载的局限性，综合了有限元法与离散元法各自的优势，不仅提升了数值模拟的预测精度，还进一步提高了计算效率，增强了实用性及信息化水平，便于应用于实际工程的模拟方法中。

附图说明

图1为本发明提供的基于FEM-DEM联合仿真的岩石爆破全过程模拟方法的流程图；

图2为本发明实施例1中提供的在FEM软件中建立的轴对称模型图；

图3为本发明实施例提供的爆破损伤云图；

图4为本发明实施例1中提供的孔壁荷载压力时程曲线图；

图5为本发明实施例1中提供的粉碎区荷载压力时程曲线图；

图6为本发明实施例1中提供的DEM软件三维真实模型示意图；

图7为本发明实施例1中提供的爆炸荷载加载示意图；

图8为本发明实施例1中提供的DEM软件计算0.1s、0.2s、0.3s、0.4s时的速度云图；

图9为本发明实施例1中提供的DEM软件最终模拟结果。

具体实施方式

下面结合附图1至图9，对本发明的具体实施方式作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例1：如图1-9所示，本发明提出了一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，具体包括以下步骤：

S1基于现场各段别炮孔爆破参数和装药结构，在FEM软件(本发明采用LS-DYNA软件)中建立爆破模型；

S2设置LS-DYNA软件中建立的爆破模型的边界条件；设置爆破模型中炸药和空气参数及状态方程、围岩的本构模型和物理力学参数；

S3对LS-DYNA软件中建立的爆破计算模型进行求解，得到各段别炮孔的孔壁爆炸荷载压力时程曲线和粉碎区爆炸压力荷载时程曲线，将孔壁和粉碎区荷载峰值与理论公式计算得到的结果进行对较；

S4在DEM软件(本发明采用3DEC软件)软件中根据现场情况建立模型，包括节理面数量的位置和地质构造；

S5设置3DEC软件中建立模型的岩石和节理的本构模型、物理力学参数和边界条件；

S6将LS-DYNA软件中计算得到的各段别炮孔的粉碎区爆炸荷载压力曲线分别施加至对应的3DEC软件建立的模型各段别炮孔粉碎区边界上；

S7在3DEC软件计算中可设置计算时间，由此可得到岩石爆破不同时间的爆破破坏形态；

S8在3DEC软件计算前，通过理论块度大小公式提前将岩体按照距离炮孔中心的远近划分为不同大小的块体，使模拟结果更符合实际，距离半径为b的炮孔距离R处的岩块的长度为：

式中：T为岩石的动态极限抗拉应变值；K为炮孔壁上的应变；β为应变波吸收系数；

具体实验数据和计算公式如下：

(1)岩土力学相关参数：现场或实验室获取；

(2)现场爆破相关参数：现场获取；

(3)炸药相关参数：从生产厂家获得；

(4)岩土力学相关参数包括：密度ρ、弹性模量e、剪切模量G、泊松比ν、粘聚力c、内摩擦角μ、抗拉强度T_e、节理的切向和法向刚度比K_s/K_n；

(5)对LS-DYNA软件中建立的模型参数进行调试求解，直至各段别炮孔的孔壁荷载和粉碎区荷载的峰值符合理论和经验公式，即炮孔的孔壁爆炸压力荷载峰值和粉碎区爆炸压力荷载峰值与理论计算值的相对误差小于10％，判定炮孔的爆炸荷载的荷载峰值符合理论和经验公式；

其中，炮孔的孔壁爆炸荷载压力理论计算公式为：

式中：P₀为炮孔的爆炸荷载压力，ρ₀为炸药密度，D为炸药爆轰速度，x为提取的荷载位置距起爆点的距离，l为药卷长度。

粉碎区外边界上径向应力峰值与距离的关系可表示为：

式中：P₀为炮孔的爆炸荷载压力，P_b为粉碎区荷载，r_c为粉碎区半径，r₀为炮孔半径，α为衰减指数。

(6)在3DEC软件模型计算前，基于哈里斯爆破理论对模型不同部位的网格尺寸进行预处理，使岩石爆破结果更符合实际；

其中，根据哈里斯等的研究，半径为r的球形装药上，距球心为R的任意一点上的切应变为：

式中：β为应变波吸收系数；K为炮孔壁上的应变并由下式表出：

式中：P₀为爆压；ρ为岩石的密度；v_p为岩石纵波声速；μ为泊松比；γ为绝热指数。

在压缩应力波作用下岩石作径向位移，由径向位移衍生的切向应变值ε超过岩石的动态极限抗拉应变值T时，岩石形成径向裂缝。距炮孔距离为R处的径向裂缝的条数为：

假设破碎岩块的最大尺寸为两条最靠近的径向裂缝之间的距离，岩块的长L为：

从式(3)(4)(5)(6)导出距离半径为b的炮孔距离R处的岩块的长度为：

实施例2：基于某现场台阶爆破试验,参阅图1，基于现场各段别炮孔爆破参数和装药结构，采用LS-DYNA软件建立轴对称模型并求解；设置边界条件(依据现场情况对模型施加无反射边界和对称边界条件)；设置模型中炸药和空气参数及状态方程、围岩的本构模型和物理力学参数，在本实施例中，围岩的本构模型选用RHT材料模型，炸药采用2号岩石乳化炸药，材料模型选用MAT_HIGH_EXPLOSIVE_BURN，该材料模型需与JWL状态方程联用。采用3DEC软件建立三维真实尺度模型，设置边界条件、岩石和节理的物理力学参数，然后将LS-DYNA软件计算得到的粉碎区荷载继承施加到3DEC软件中建立的模型对应的粉碎区边界，进行爆破动力响应分析。

参阅图2，图2为实施例在LS-DYNA软件中建立的轴对称模型，其爆破参数为：根据现场实际情况建立尺寸为17m×14m的模型，对模型精细划分得到水平方向尺寸最小为1cm，最大在5cm以内，在垂直方向尺寸为1cm的网格，装药直径为70mm，钻孔直径为115mm，药卷长度为10m，堵塞长度为4m，炸药密度为：1300kg/m³，爆速4000m/s，A为214.4GPa，B为0.182GPa，对LS-DYNA软件中建立的轴对称模型进行求解，直至各段别炮孔的爆炸荷载压力曲线的荷载峰值符合理论和经验公式，否则调整轴对称模型中炸药和空气参数及状态方程、围岩的本构模型和物理力学参数。

参阅图3，为LS-DYNA软件中建立的轴对称模型爆破后的损伤形状和范围；

参阅图4，为LS-DYNA软件中建立的轴对称模型爆破后在对应孔壁单元提取的爆破荷载曲线，在本实施例中，现场炮孔孔径为115mm，药径为70mm，为不耦合装药结构，通过公式(1)计算得孔壁压力峰值为2.6GPa，数值模拟得孔壁压力峰值为2.7GPa，相对误差为3.8％，满足要求。

参阅图5，为LS-DYNA软件中建立的轴对称模型爆破后在对应粉碎区单元提取的爆破荷载曲线，在本实施例中，将孔壁荷载按照理论公式(2)进行衰减计算，理论计算得粉碎区荷载峰值为517MPa，数值模拟得荷载峰值为550MPa，相对误差为6％，满足要求。

参阅图6，基于现场实际工程形状、尺寸和研究区域围岩类别，采用3DEC软件建立三维真实尺度模型进行求解。设置3DEC软件中建立的三维真实尺度模型的岩石、节理的物理力学参数和边界条件(设置固定边界和粘滞边界)。本实施例的岩石和节理力学参数为：岩石密度为2700kg/m³，体积模量为25GPa，剪切模量为18.75GPa，粘聚力为2MPa，内摩擦角为35°，抗拉强度为0.5MPa，泊松比为0.25，节理粘聚力为2MPa，节理内摩擦角为35°，节理抗拉强度为0.5MPa，节理法向刚度为18.75GPa，节理剪切刚度为25GPa。同时，在3DEC软件模型计算前，基于哈里斯爆破理论对模型不同部位的网格尺寸进行预处理，使岩石爆破结果更符合实际，将岩石参数带入公式进行计算，选取0.1m、0.15m、0.2m、0.3m、0.5m和1m为划分块度的依据，按照距爆心的距离由近到远进行划分，将距离炮孔0～1.4m划分为10cm的块体，距离炮孔1.4～2m划分为15cm的块体，距离炮孔2～2.8m划分为20cm的块体，距离炮孔2.8～3.5m划分为30cm的块体，距离炮孔3.5～4.5m划分为50cm的块体，距离炮孔4.5m以外划分为1m的块体；

参阅图7，将各段别炮孔的粉碎区爆炸压力荷载曲线分别施加至对应的3DEC软件中建立的三维真实尺度模型各段别粉碎区边界，某一粉碎区边界中，考虑爆轰传播的时间效应，P(t)表示随时间变化的炮孔爆炸荷载压力，t＝L/D表示距起爆点距离为L的炸药起爆滞后时间。将LS-DYNA模型中提取出的粉碎区爆炸荷载时程曲线时间效应依次施加于3DEC软件中的求解模型粉碎区边界上。

参阅图8，为3DEC软件中建立的三维真实尺度模型在0.1s、0.2s、0.3s、0.4s时的速度云图。

参阅图9，为3DEC软件中建立的三维真实尺度模型最终的爆堆形态。

综合来看，本发明提供的基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，具有以下优点：

本方法首先根据岩石爆破工程的爆破区域地形及工况条件，确定爆破区域的岩体参数和爆破参数；其中，岩体参数至少包括密度、弹性模量、剪切模量、泊松比、粘聚力、内摩擦角、抗拉强度、节理的切向和法向刚度，节理面位置和形态；爆破参数至少包括：多排炮孔的位置数据、炸药的尺寸、种类、采用的装药形式、不耦合装药时采用的介质、起爆雷管位置以及是否采用多段装药；然后根据获取的爆破参数通过FEM软件模拟计算得到爆破过程的孔壁荷载曲线、粉碎区荷载曲线与损伤情况；在DEM软件中建立现场爆破台阶模型，可在计算前提前将模型按照块度大小提前用节理面切割处理，使最后的计算结果更符合实际；最后将FEM软件中模拟计算得到的粉碎区荷载数据施加到DEM软件中，通过DEM软件来模拟爆破效果，包括爆破时的破碎块体抛掷速度、方向和抛掷距离以及爆堆形态。

本方法避免了单独使用有限元法无法模拟爆破过程中破碎块体的抛掷和爆破效果，以及单独使用离散元法无法准确获得爆炸荷载的局限性，综合了有限元法与离散元法各自的优势，不仅提升了数值模拟的预测精度，还进一步提高了计算效率，增强了实用性及信息化水平，便于应用于实际工程的模拟方法中。

以上所述实施例仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换，均属于本发明的保护范围。

Claims (7)

1.一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于，包括以下步骤：

基于有限元法和实际各段别炮孔爆破参数和装药结构，建立有限元模型；

设置有限元模型的无反射边界和对称边界，并确定有限元模型中的炸药和空气参数及状态方程、岩石的本构模型和物理力学参数；

求解有限元模型，模拟炸药爆轰过程，获取各段别炮孔的孔壁压力荷载时程曲线和粉碎区边界压力荷载时程曲线，以及岩石爆破损伤分布和振动情况；

基于离散元法和实际工程形状、尺寸和岩体类型，建立离散元模型；

设置离散元模型中的围岩本构模型及物理力学参数、炮孔布置参数；

将有限元模型输出的各段别粉碎区边界荷载时程曲线分别施加至离散元模型各段别粉碎区边界，设置离散元模型的粘滞边界和固定边界，对离散元模型进行爆破动力响应分析，获取岩石爆破抛掷形态和爆破效果；

其中，所述获取各段别炮孔的孔壁压力荷载时程曲线和粉碎区边界压力荷载时程曲线，具体包括：

求解有限元模型，得到各段别炮孔孔壁的爆炸荷载压力时程曲线和粉碎区爆炸荷载压力时程曲线，将孔壁压力荷载峰值和粉碎区压力荷载峰值与理论计算值进行对较，当孔壁压力荷载峰值和粉碎区压力荷载峰值与理论计算值的相对误差小于10％，则判定爆炸荷载压力曲线的荷载峰值符合理论公式；

其中，炮孔的孔壁荷载峰值理论计算公式为：

式中：P₀为炮孔的爆炸荷载压力，ρ₀为炸药密度，D为炸药爆轰速度，x为爆炸压力与起爆点的距离，l为药卷长度；

粉碎区外边界上径向应力峰值与距离的关系可表示为：

式中：P_b为粉碎区荷载压力，P₀为炮孔壁的爆炸荷载压力，r_c为粉碎区半径，r₀为炮孔半径，α为衰减指数。

2.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于：

所述有限元模型的建立是利用LS-DYNA软件；

所述离散元模型的建立是利用3DEC软件。

3.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于：

通过现场或实验室获取实际各段别炮孔爆破参数、装药结构、实际工程形状、尺寸和岩体类型。

4.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于：

岩土物理力学参数包括：密度ρ、弹性模量e、剪切模量G、泊松比ν、粘聚力c、内摩擦角μ、抗拉强度T_e、节理的切向和法向刚度比K_s/K_n。

5.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于：

在对离散元模型进行爆破动力响应分析前，基于哈里斯爆破理论对离散元模型不同部位的网格尺寸进行预处理，使岩石爆破结果更符合实际；

其中，根据哈里斯爆破理论，半径为r的球形装药上，距球心为R的任意一点上的切应变为：

式中：β为应变波吸收系数；K为炮孔壁上的应变；

炮孔壁上的应变K计算公式为：

式中：P₀为爆压；ρ为岩石的密度；v_p为岩石纵波声速；μ为泊松比；γ为绝热指数；

在压缩应力波作用下岩石作径向位移，由径向位移衍生的切向应变值ε超过岩石的动态极限抗拉应变值T时，岩石形成径向裂缝；距炮孔距离为R处的径向裂缝的条数为：

设定破碎岩块的最大尺寸为两条最靠近的径向裂缝之间的距离，岩块的长L为：

因此，距离半径为b的炮孔距离R处的岩块的长度为：

基于上式得距爆心距离与预测爆破块度关系，选取不同块度的大小，计算出爆心距离并划分块度大小对应的区域范围，将模型按照距离炮孔中心的远近从小到大进行块度划分。

6.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于，所述设置有限元模型的无反射边界、对称边界、炸药和空气参数及状态方程、岩石的本构模型和物理力学参数，具体包括：

根据现场实际情况设置有限元模型的尺寸和参数，并对有限元模型网格进行精细划分，从炮孔近区开始由近到远采用渐变网格，网格最大控制在5cm以内。

7.如权利要求1所述的一种基于FEM-DEM联合仿真的岩石爆破全过程模拟方法，其特征在于，所述对离散元模型进行爆破动力响应分析，获取岩石爆破抛掷形态和爆破效果，具体包括：

通过设置对离散元模型的计算时间，获取岩石爆破不同时间的抛掷形态和爆破效果。

Images