一种通用爆破数值模拟方法与流程

1.本发明涉及数值模拟方法，具体涉及一种通用爆破数值模拟方法。


背景技术：

2.爆破技术在矿山、铁路、交通、水利水电工程以及城市各种基础设施建设及改造中被广泛采用。爆破工程中的炸药是迄今为止人类所用工程中最危险、最难驾驭的器材之一，如何科学又安全地运用爆破技术来达到特定的工程目的，是目前仍需关注的问题。爆破效果的好坏与爆破方案、爆破参数的选取密切相关。因此，在实施现场爆破施工前，需要一种有效的手段，来对爆破方案和爆破参数的选取进行评价，并根据评价结果优化爆破方案及爆破参数。随着计算机硬件及计算仿真技术的发展，数值模拟方法成为评价爆破方案及爆破参数的有效手段，数值模拟的结果可以为优化爆破方案提供指导。由于爆破过程中涉及材料的断裂、破碎、碰撞，常规的数值方法如有限元法、颗粒离散元法很难对爆破问题进行很好地建模。例如有限元方法，通常采用生死单元法来对爆破损伤进行建模，但是这会导致模拟过程中质量不守恒，另外较难对爆破过程产生的块体飞溅等进行建模。块体离散元虽然对于模拟块体的碰撞、运动、堆积等有较好的效果，但在模拟块体断裂及破碎方面能力不足。颗粒离散元虽然能够较好地对材料断裂、破碎，运动、碰撞、堆积等过程进行建模，但输入参数无法和试验测得的材料参数对应，需要进行繁琐的参数标定。因此，上述这些方法在模拟岩土材料及混凝土材料爆破方面仍然存在很多问题。


技术实现要素：

3.本发明的目的在于，针对现有技术的不足，提供一种通用的爆破模拟方法，能够对土、岩石、混凝土或其他材料中的爆破进行数值模拟。
4.本发明采用的技术方案为：一种通用爆破数值模拟方法，包括如下步骤：
5.s1：建立爆破几何模型，并在模型中以炮孔边界划分炮孔区域和炮孔外区域；
6.s2：划分实体单元有限元网格：分别对炮孔区域和炮孔外区域划分实体单元有限元网格，两区域模型边界两侧的实体单元不共用节点；
7.s3：在炮孔外区域相邻的实体单元之间插入节理单元；
8.s4：在装药后炮孔区域内，对药包的传爆过程进行建模，并计算爆轰波首次传至炮孔区域内各实体单元的时刻；
9.s5：以爆轰波首次传至炮孔区域内第l个实体单元的时刻t
l
为起始时刻，对炮孔内相应的第l个实体单元施加一个随时间变化的应变增量；
10.s6：根据应力-应变本构关系，将应变增量转换成应力增量；
11.s7:将应力增量施加相应的第l个实体单元上进行力学断裂计算，计算炮孔区域外的所有节理单元的拉应力或剪应力，由此获得应力场，并判断该拉应力或剪应力是否超过爆破外区域结构的抗拉强度或者抗剪强度，由此推断哪些节理单元断裂，即完成一个时间步内爆破模拟计算；
12.s8：循环重复步骤s3-s7，对整个炮孔区域的实体单元进行力学断裂计算，每计算一次判断节理单元是否断裂，获得整个炮孔区域外的应力场，实现爆破全过程的模拟。
13.按上述方案，在s2中，实体单元有限元网格可为三角形、四边形、五边形，或其他任意多边形中的一种或多种。
14.按上述方案，在s4中，爆轰波传至第l个实体单元的时刻t
l
的计算方法：
15.设起爆点的坐标为(x0,y0)，炮孔区域内第l个实体单元中心点的坐标为(x
l
,y
l
)，爆轰波的传波速度为v(单位为m/s)，沿着炮孔轴向且与爆轰波传播方向一致的单位向量为n＝(e1,e2)，此为已知量，则t
l
为：
[0016][0017]
按上述方案，在s5中，应变增量δε
ij
与时间的关系可表示为：
[0018][0019]
式中：l为炮孔区域内实体单元的编号；δ
ij
为置换张量；ε
max
为实体单元发生的最大膨胀应变；f(x)为最大值为1、最小值为0的无量纲函数。
[0020]
按上述方案，在s6中，施加在爆孔内任意三角形单元l的应力增量δσ
ij
为：
[0021][0022]
式(3)中：k
*
为体积模量，对于平面应力问题，k
*
＝6kg/(3k 4g)，对于平面应变问题k
*
＝k，k为体积模量。
[0023]
按上述方案，在s7中，力学断裂计算及应力场的获得可采用有限元方法、混合有限-离散元方法、离散元方法、非连续变形分析方法或数值流行方法计算。
[0024]
本发明的有益效果为：本发明建立模型、划分实体单元网格，并对爆破过程进行模拟，炮孔区域内实体单元进行断裂计算并获得应力场，实现对爆破过程的模拟；本发明利用裂缝扩展沿着实体单元边界这一特征，在炮孔外区域相邻的实体单元之间插入起粘结作用的节理单元或界面单元，在模拟裂缝起裂和扩展和交汇时，不用追踪裂缝扩展和网格重划分，可以非常方便地模拟任意复杂的断裂及破碎问题。
附图说明
[0025]
图1为本实施例中几何模型示意图。
[0026]
图2为模型中炮孔区域和炮孔外区域的实体单元网格划分示意图。
[0027]
图3为节理单元的示意图。
[0028]
图4为本实施例中炮孔区域内实体单元起始时刻t
l
的计算。
[0029]
图5为本实施例中函数f(t)的典型曲线一。
[0030]
图6为本实施例中函数f(t)的典型曲线二。
具体实施方式
[0031]
为了更好地理解本发明，下面结合附图和具体实施例对本发明作进一步地描述。
[0032]
实施例
[0033]
一种通用爆破数值模拟方法，包括如下步骤：
[0034]
s1：建立爆破二维几何模型，并在模型中以炮孔边界划分炮孔区域和炮孔外区域，其中，位于炮孔边界内部的区域为炮孔区域，其余为炮孔外区域。一个几何模型中，可以有任意多个炮孔区域，为了简化起见，此处仅仅展示一个炮孔区域，如图1所示。
[0035]
s2：划分实体单元有限元网格：分别对炮孔区域和炮孔外区域划分实体单元有限元网格。
[0036]
实体单元有限元网格的形状和尺寸可不同，可为三角形、四边形、五边形、二维voronoi单元，或其他任意多边形中的一种或多种；炮孔区域的实体单元有限元网格完全独立于炮孔外区域的实体单元有限元网格，两区域模型边界两侧的实体单元不共用节点。本实施例中，如图2所示，实体单元有限元网格为三角形。
[0037]
s3：在炮孔外区域相邻的实体单元之间插入起粘结作用的节理单元或界面单元，如图3所示。若相邻的实体单元中插入节理单元或界面单元，则可通过节理单元或界面单元的断裂来模拟裂纹的萌生及扩展(若相邻的实体单元不插入节理单元和界面单元，则通过分离公共面两侧的实体单元来模拟裂缝的萌生及扩展，此为现有技术，这里不再赘述)。
[0038]
s4：在装药后的炮孔区域内，对药包的传爆过程进行建模，根据传爆速度、实体单元中心点的坐标及起爆点的坐标确定激活哪些实体单元需要施加应力增量(具体为：通过实体单元中心点离起爆点的轴向距离确定爆轰波传至对应实体单元的起始时刻，再结合爆轰持续时间，即可确定特定时刻哪些实体单元需要时间应力增量)，计算爆轰波首次传至炮孔区域内各个实体单元的时刻。
[0039]
在爆速已知的情况下，爆轰波自起爆点沿着炮孔轴向传播。如图4所示，设起爆点的坐标为(x0,y0)，炮孔区域内第l个实体单元中心点的坐标为(x
l
,y
l
)，爆轰波的传波速度为v(单位为m/s)，沿着炮孔轴向且与爆轰波传播方向一致的单位向量为n＝(e1,e2)，此为已知量，则爆轰波传至第l个实体单元的时刻t
l
为：
[0040][0041]
本实施例中，炮孔区域内第l个实体单元的中心点坐标可通过该实体单元三角形网格三个节点的坐标取平均值求得。
[0042]
s5：以爆轰波首次传至炮孔区域内第l个实体单元的时刻t
l
为起始时刻，对炮孔内相应的第l个实体单元施加一个随时间变化的应变增量δε
ij
(单位为ij表示是一个张量)，该应变增量δε
ij
可表示成时间的函数。
[0043]
应变增量δε
ij
与时间的关系可表示为：
[0044][0045]
式中：l为炮孔区域内实体单元的编号；δ
ij
为置换张量；ε
max
为实体单元发生的最大膨胀应变；f(x)为最大值为1、最小值为0的无量纲函数，可以为实体单元的分段函数或其他任意函数形式，只要这种函数形式能够表征实际的爆破作用一致即可，f(x)的典型函数形式如图5和图6所示。
[0046]
s6：根据应力-应变本构关系，将应变增量转换成应力增量δσ
ij
。
[0047]
本实施例中，施加在爆孔内任意三角形实体单元l的应力增量δσ
ij
为：
[0048][0049]
式(3)中，k
*
为体积模量，对于平面应力问题，k
*
＝6kg/(3k 4g)，对于平面应变问题k
*
＝k，k是体积模量。本实施例属于平面应变问题。
[0050]
s7：将(3)式的应力增量施加到相应的第l个实体单元上进行力学断裂计算，计算炮孔区域外的所有节理单元的拉应力或剪应力，由此获得应力场，并判断该拉应力或剪应力是否超过爆破外区域结构的抗拉强度或者抗剪强度，由此推断哪些节理单元断裂，即完成一个时间步内爆破模拟计算。
[0051]
炮孔区域内实体单元在上述压力增量的作用下，会向炮孔周边运动并与炮孔边界外的实体单元(炮孔外区域的实体单元)接触，使得炮孔外区域内的节理单元或界面单元产生断裂，并将断裂的节理单元或界面单元删去，或者公共边两侧的实体单元发生分离，从而模拟爆孔外区域产生裂缝扩展甚至破碎。
[0052]
力学断裂计算及应力场的获得可采用有限元(fem)方法、混合有限-离散元(fdem)方法、离散元(dem)方法、非连续变形分析(dda)方法或数值流行(nmm)方法计算，上述方法均为现有技术，这里不再赘述。本实施例中以fdem方法进行力学断裂计算为例，将相应应力增量施加在炮孔区域内的实体形单元上，然后进行fdem力学断裂计算(fdem可以处理断裂及接触)，判断节理单元是否断裂，即完成一个时间步内爆破模拟计算。
[0053]
s8：循环重复步骤s3-s7，对整个炮孔区域的所有实体单元进行力学断裂计算，每计算一次判断炮孔区域外的节理单元是否断裂，获得整个炮孔区域外的应力场。
[0054]
本发明中，不同材料的抗拉强度/者抗剪强度不同，将计算获得的节理单元的拉应力或剪应力与不同材料的抗拉强度或者抗剪强度进行对比并判断，如此循环往复，即可实现对土体、岩石、混凝土等不同材料爆破全过程的模拟。
[0055]
在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。