冲击载荷作用下岩体拉-压损伤破坏的边坡抛掷爆破模拟

梁瑞; 周文海; 余建平; 李珍宝; 杜超飞; 王敦繁

doi:10.11858/gywlxb.20180535

冲击载荷作用下岩体拉-压损伤破坏的边坡抛掷爆破模拟

doi: 10.11858/gywlxb.20180535

1.
兰州理工大学石油化工学院，甘肃兰州　730050
2.
浙江大学海洋学院，浙江舟山　316021

基金项目: 国家自然科学基金（51566010，51076061）；甘肃省自然科学基金（B061709）

详细信息

作者简介:
梁　瑞（1968－），男，博士，教授，主要从事安全工程与工程爆破研究. E-mail: liangr@lut.cn

通讯作者:
周文海（1989－），男，硕士，助教，主要从事安全工程与工程爆破研究. E-mail: 18394499554@139.com

中图分类号: O347.1; TU457
计量
- 文章访问数: 8241
- HTML全文浏览量: 2993
- PDF下载量: 37
出版历程
- 收稿日期: 2018-04-03
- 修回日期: 2018-04-19
- 刊出日期: 2019-02-25

Numerical Simulation of Rock Tension-Compression Fracture Caused by Impact Load during Slope Casting Blast

1.
School of Petrochemical Technology, Lanzhou University of Technology, Lanzhou 730050, China
2.
Ocean College, Zhejiang University, Zhoushan 316021, China

摘要

摘要: 为了获得边坡台阶爆破时岩体在冲击载荷作用下破碎和抛掷过程中各物理参量的变化规律，将构建的岩体动态拉-压损伤本构关系嵌入模拟软件进行数值分析。结果表明：振动波三轴合成速率曲线与质点振动位移曲线的变化趋势所对应的时间节点和步长基本吻合，可作为降振减灾的判定指标；岩体最早于0.6 ms左右从坡脚位置产生裂纹，直至12.5 ms左右裂纹裂隙完成扩展，所形成的炮孔粉碎区半径约28 cm；抛掷块状分离现象从炮孔药包的中间部位开始，最大抛掷速度集中于该部位至边坡自由面之间的垂直区域内，边坡自由面抛掷速度小于炮孔周围岩块抛掷速度，导致抛掷过程中形成二次挤压破碎现象；破碎大块主要来源于边坡坡脚、炸药与堵塞物接触面两侧围岩以及台阶顶部自由面处，大块块体直径分布于1.6～2.7 m范围。
- 冲击载荷 /
- 拉-压损伤破坏 /
- 损伤耦合 /
- 抛掷爆破
Abstract: In order to study the change law of physical parameters in the process of rock fracturing and throwing during slope bench blasting, the equation of rock damage under dynamic tension-compression effect was established and applied to numerical analysis. The results showed that the tendency of time node and step size in simulation was basically identical with the triaxial synthetic rate curve of vibration wave and particle vibration displacement, which can be used as a criterion for reducing vibration and decreasing disaster. Cracks were formed in the foot of slope at about 0.6 ms and completely extended at about 12.5 ms. The pulverizing area radius around the blast hole was 28 cm. The rock separation phenomenon was preliminarily observed at the middle part of the blast hole. The maximum throwing velocity was distributed in the vertical region between this part and the free surface of the slope. The throwing velocity at the free surface was less than that of the rocks around the blast hole, which results in the secondary crushing phenomenon during the throwing process. The large bulk rocks were mainly produced in the toe of slope, the surrounding rock on both sides of the contact surface between explosive and plug, and the free-surface at the top of the step. The range of large rock diameter in the process of blasting was 1.6–2.7 m.
- impact load /
- tension-compression fracture /
- damage coupling /
- casting blast

HTML全文

工程爆破是一个极为复杂的瞬态动力破坏过程，很难用具体的理论去全面解释，所以借助软件运用数值法预估冲击载荷作用下脆性岩体爆破的整个过程成为当前工程爆破的研究热点。如Budiansky等^[1]为研究岩体爆破过程中裂纹裂隙的扩展情况，运用细观力学对岩体中存在的裂隙群进行分析，在此基础上建立了爆破冲击载荷作用下岩体损伤本构关系；朱传云等^[2]将非连续变形分析法（Discontinuous Deformation Analysis，DDA）运用到台阶抛掷爆破模拟研究中，通过分析给出非均质岩体在爆轰气体准静态作用下岩块抛掷轨迹以及形成爆堆的堆积情况；Munjiza^[3]采用FEM/DEM（有限元法/离散元法）模拟爆破过程中裂纹扩展情况，考虑了岩体本身存在的节理面，将网格划分为三角形单元，相邻单元边上嵌入粘结的节理单元，通过节理断裂模拟裂隙拓展；黄永辉等^[4]采用RHT材料本构模型对台阶抛掷爆破进行数值模拟，获得了坡面岩块抛掷轨迹和爆堆形状，并运用岩体塑性拉伸损伤理论，给出了岩块抛掷速率传播规律。

学者们虽然通过建立岩体在冲击载荷作用下单纯受拉时的不同损伤本构关系研究了岩质边坡爆破情况，但是却未考虑岩体在受压状态下的硬化效应，而合理的屈服破坏本构关系对于岩体爆破模拟结果的准确性至关重要。在岩体受拉状态下，若用Taylor-Chen-Kuszmaul（TCK）损伤演化方程表征岩体破坏规律，并设定材料响应服从弹性胡克定律，则简单明了；在受压状态下，设定其应力-应变行为服从双线性、弹塑性本构，并遵循Mises屈服准则，解决了以往该本构关系不考虑岩体在受压状态下的硬化效应的问题，从而较合理地反映爆炸载荷下岩体的弹塑性变形规律。

1. 岩体宏细观缺陷耦合的节理拉压损伤本构模型

1.1 拉裂损伤计算模型

岩体本身存在各式各样的裂纹裂隙，也称原始损伤，在爆破载荷拉应力作用下该损伤进一步扩展。Budiansky等^[1]基于“自洽法”求出了损伤介质的弹性模量。TCK模型提出在爆炸载荷拉应力作用下岩体介质中的裂纹将被激活，并且该裂纹相互挤压碰撞形成相互作用力。存在原始损伤的同性岩体有效体积模量可表示为

$\frac{K}{{{K_0}}} = 1 - \frac{{16}}{9}{f_1}(\nu ){C_{\rm{d}}} = 1 - D$

(1)

式中：K、 ${K_0}$ 分别为损伤介质体积模量和初始体积模量； ${f_1}(\nu ) = {{(1 - {\nu ^2})} / {(1 - 2\nu )}}$ ，其中 $\nu$ 为泊松比；C_d为介质裂纹密度；D为损伤标量。

Taylor等^[5]提出岩体介质在冲击载荷作用下的泊松比变化关系式

$\nu = {v_0}\left(1 - \frac{{16}}{9}{C_{\rm{d}}}\right)$

(2)

式中： ${v_0}$ 为介质的初始泊松比。

Kipp等^[6]提出可用Weibull函数精确计算在冲击拉应力作用下单位介质的原始损伤裂纹被激活数目N，即

$N = \xi {\left(\frac{p}{{3k}}\right)^m}$

(3)

式中：p为冲击应力； $\xi$ 、m为Weibull函数的分布参数。

Taylor等^[5]曾提出岩体等脆性介质在动态冲击载荷作用下的碎块半径公式

$r = \frac{1}{2}{\left( {\frac{{\sqrt {20} {K_{{\rm{IC}}}}}}{{\rho C{{\dot \varepsilon }_{v,\,{\rm max}}}}}} \right)^{{2 / 3}}}$

(4)

式中: ${K_{{\rm{IC}}}}$ 、 $\rho$ 分别为介质的断裂韧性和密度，C为纵波波速， ${\dot \varepsilon _{v,\max }}$ 为介质断裂时的体应变率。碎块平均半径r和裂纹密度C_d之间的函数关系可表示为

${C_{\rm{d}}} = N\varpi {r^3}$

(5)

式中： $\varpi$ 为比例系数。

把（3）式和（4）式代入（5）式可得

${C_{\rm{d}}} = \xi {\left(\frac{p}{{3k}}\right)^m}\left[ {\varpi \cdot \frac{1}{8}{{\left( {\frac{{\sqrt {20} {K_{{\rm{IC}}}}}}{{\rho C{{\dot \varepsilon }_{v,\max }}}}} \right)}^2}} \right] = \frac{5}{2}\xi \varpi {\left(\frac{p}{{3k}}\right)^m}{\left(\frac{{{K_{{\rm{IC}}}}}}{{\rho C{{\dot \varepsilon }_{v,\max }}}}\right)^2}$

(6)

设 $\kappa = \xi \varpi$ ，则

${C_{\rm{d}}} = \frac{5}{2}\kappa {\left(\frac{p}{{3k}}\right)^m}{\left(\frac{{{K_{{\rm{IC}}}}}}{{\rho C}}\right)^2}{\dot \varepsilon ^{ - 2}}_{v,\max }$

(7)

由（1）式可求得拉裂损伤函数式

$D = \frac{{16}}{9}{f_1}(v){C_{\rm{d}}}$

(8)

根据弹性力学理论中以应力表示应变的广义胡克定律，有

${\varepsilon _x} = \frac{{{\sigma _x}}}{E} - \frac{\nu }{E}({\sigma _y} + {\sigma _z})$

(9)

式中： ${\varepsilon _x}$ 为x方向应变， ${\sigma _x}$ 、 ${\sigma _y}$ 、 ${\sigma _z}$ 分别为x、y、z轴方向应力，E为弹性模量。

（9）式可改写为

${\varepsilon _x} = \frac{{1 + \nu }}{E}{\sigma _x} - \frac{\nu }{E}({\sigma _x} + {\sigma _y} + {\sigma _z})$

(10)

式中： ${\sigma _x} = {s_x} + \sigma _{\rm{m}}$ ， $3{\sigma _{\rm{m}}}{\rm{ = }}{\sigma _x} + {\sigma _y} + {\sigma _z}$ ； ${s_x}$ 为x方向偏主应力， $\sigma _{\rm{m}}$ 为平均正应力。将该等式代入（10）式可得

${\varepsilon _x} = \frac{{1 + \nu }}{E}({s_x} + {\sigma _{\rm{m}}}) - 3\frac{\nu }{E}{\sigma _{\rm{m}}} = \frac{{1 + \nu }}{E}{s_x} + \frac{{1 - 2\nu }}{E}{\sigma _{\rm{m}}}$

(11)

令 $K = \dfrac{E}{{3(1 - 2\nu )}}$ 为体积模量， $G = \dfrac{E}{{2(1 + \nu )}}$ 为剪切弹性模量，将K、G代入（11）式可得

${\varepsilon _x} = \frac{{s_x}}{{2G}} + \frac{{{\sigma _{\rm{m}}}}}{{3K}}$

(12)

在弹性力学中，以应力表示应变的广义胡克定律为

$\left\{ \begin{aligned} & {\varepsilon _x} = \frac{{{\sigma _x}}}{E} - \frac{\nu }{E}({\sigma _y} + {\sigma _z})\\ & {\varepsilon _y} = \frac{{{\sigma _y}}}{E} - \frac{\nu }{E}({\sigma _x} + {\sigma _z})\\ & {\varepsilon _z} = \frac{{{\sigma _z}}}{E} - \frac{\nu }{E}({\sigma _y} + {\sigma _x})\\ & {\gamma _{xy}} = \frac{{{\tau _{xy}}}}{G}\\ & {\gamma _{yz}} = \frac{{{\tau _{yz}}}}{G}\\ & {\gamma _{zx}} = \frac{{{\tau _{zx}}}}{G} \end{aligned} \right.$

(13)

对用应力表示应变的广义胡克定律（即（13）式）中其他5式按照（9）式～（12）式的方法改写，最后可将K、G形式的本构关系张量形式写成如下统一形式

${\varepsilon _{ij}} = \frac{{s_{ij}}}{{2G}} + \frac{{{\sigma _{\rm{m}}}}}{{3K}}{{{\text{δ}}} _{ij}}$

(14)

式中： ${\varepsilon _{ij}}$ 为应变张量， $s_{ij}$ 为偏应力张量， ${{{\text{δ}}} _{ij}}$ 为Kronecker记号。以应变表示应力的本构关系张量下标式为

${\sigma _{ij}} = 2G{e_{ij}} + 3K{\varepsilon _{\rm{m}}}{{{\text{δ}}} _{ij}}$

(15)

由弹性体积变形定律可知 ${\sigma _{\rm{m}}} = K{\varepsilon _v} = 3K{\varepsilon _{\rm{m}}}$ ，故岩体拉应力损伤-应变本构关系可表示为

$\Delta {\sigma _{ij}} = {K_{\rm d}}{{{\text{δ}}}_{ij}}\Delta {\varepsilon _v} + 2{G_{\rm d}}\Delta {e_{ij}}$

(16)

式中： ${K_{\rm {d} }}$ 和 ${G_{\rm {d} }}$ 代表损伤介质的模量； $\Delta {\varepsilon _v}$ 和 $\Delta {e_{ij}}$ 为体应变增量和偏应变增量， $\Delta {\varepsilon _v} = 3{\varepsilon _{\rm{m}}}$ 。

1.2 压应力作用下硬化效应的本构关系

假设岩体介质在爆轰气体压缩作用下服从各向同性硬化的弹塑性本构关系，该屈服方程可表示为^[7]

$\varPhi = \sigma _i^2 - \sigma _y^2$

(17)

式中： ${\sigma _i} = \sqrt {\dfrac{3}{2}{s_{ij}}}$ 。根据广义胡克定律，用张量形式表示偏应力 ${s_{ij}}$

${s_{ij}} = {\sigma _{ij}} - \frac{1}{3}{\sigma _{kk}}{{{\text{δ}}} _{ij}}$

(18)

爆炸过程中，在动力载荷作用下岩体介质有效塑性应变 $\varepsilon _{{\rm{eff}}}^{\rm{P}}$ 和达到断裂破坏的屈服极限 ${\sigma _{\rm{y}}}$ 的函数关系式为

$\varepsilon _{{\rm{eff}}}^{\rm{P}} = \int_0^t {\sqrt {\frac{3}{2}} \dot \varepsilon _{ij}^{\rm{P}}{\rm{d}}t}$

(19)

${\sigma _{\rm{y}}} = {\sigma _0} + \beta {E_{\rm{p}}}\varepsilon _{{\rm{eff}}}^{\rm{P}}$

(20)

式中： ${E_{\rm{p}}}$ 为塑性硬化模量，通常可取0.1E； $\beta$ 为硬化参数；而 $\dot \varepsilon _{ij}^{\rm{P}} = {\dot \varepsilon _{ij}} - \dot \varepsilon _{ij}^{\rm{e}}$ 。（17）式中，当 $\varPhi \leqslant 0$ 时为弹性或中性加载； $\varPhi > 0$ 为塑性硬化。

1.3 拉-压损伤本构模型嵌入模拟软件

大型有限元分析软件LS-DYNA对于瞬时动力和大变形破坏模拟自带丰富的材料库，同时还为用户提供方便的模型导入窗口。为此本研究将存在宏细观缺陷的边坡岩体拉裂损伤计算模型和压应力作用下硬化效应的本构关系耦合，导入该模拟软件进行台阶爆破计算，模拟爆炸冲击作用下振动波速度传播规律、裂纹裂隙扩展状态、爆轰压力分布、抛掷运动轨迹以及岩体破碎块度情况。具体流程如图1所示。

图 1 损伤模型嵌入流程

Figure 1. Flowchart for damage model

下载: 全尺寸图片幻灯片

2. 模型及材料参数选取

2.1 模型构建

以某一实际高陡边坡台阶抛掷爆破为背景，建立单孔起爆三维实体模型。设置无反射边界条件约束台阶底面以及左右边界面，整个过程对所有材料施加重力约束，起爆方式为孔底起爆，抛掷方向设定在垂直于边坡自由面法线方向；台阶高度取12 m，坡面角70°，最小抵抗线设置3.5 m，孔径150 mm，孔深14.5 m，堵塞长度4 m。参照Ma等^[8]的研究结果确定模型整体规格，台阶下部平盘与台阶上部平盘均取10 m长，坡面高度取12 m，下层台阶高8 m。为节省运算时间和空间，并且动态显示爆破过程中岩体裂纹裂隙扩展和抛掷趋势，将模型厚度定义为一个单元厚度，具体模型如图2所示。

图 2 台阶平面模型

Figure 2. Model of the bench plane

下载: 全尺寸图片幻灯片

岩石材料选取典型的弹塑性材料作为研究对象，即采用MAT_PLASTIC_KINEMATIC模型。岩石的主要物理力学参数见表1，其中 $\rho$ 为密度。

表 1 岩石主要物理力学参数

Table 1. Physical and mechanical parameters of rock

$\rho$ /(g·cm^–3)	E/GPa	$\nu$	${\sigma _0}$ /MPa	$\beta$	G/GPa	K/GPa	${K_{{\rm{IC}}}}$ /(10⁵ N·m^-3/2)
2.7	18.23	0.23	102	1.0	7.41	11.25	5.32

下载: 导出CSV

| 显示表格

2.2 炸药材料及模型

炸药选取2^#岩石乳化炸药，其材料类型为HIGH_EXPLOSIVE_BURN。炸药状态方程选取不考虑炸药产物成分的JWL状态方程，具体参数如表2所示，其中：A、B、R₁、R₂、ω为材料参数，E₀为初始比内能，p为压力，V₀为爆轰前炸药的初始比容。

表 2 炸药主要参数

Table 2. Parameters of explosive

$\rho$ /(g∙cm^-3)	p/GPa	A/GPa	B/GPa	${R_1}$	${R_2}$	$\omega$	E₀/GPa	${V_0}$
1.2	27	52.1	0.182	4.2	0.9	1.1	4.192	1.0

下载: 导出CSV

| 显示表格

炸药爆炸时JWL状态方程的具体形式为^[9]

$p = A\left(1 - \frac{\omega }{{R_1V}}\right){{\rm e}^{ - {R_1}V}} + B\left(1 - \frac{\omega }{{R_2V}}\right){{\rm e}^{ - {R_2}V}} + \frac{{\omega E}}{V}$

(21)

式中：V为爆轰产物的相对比容。

3. 数值模拟结果分析

3.1 合成速度分析

爆破降振作为当下工程爆破研究中最广泛的课题之一，研究人员通常将装药量、孔排间距、起爆方式、延时、振动速率、振动频谱等作为研究指标。依据试验的可行性和理论的成熟度，可将通过测振试验得到的爆破振动速率，特别是垂直方向（y轴）的振动速率作为振动大小（质点位移）的评判标准^[10-14]。然而，分析现场试验实测的振动速率-时间曲线和振动位移-时间曲线变化趋势发现，振动波三轴合成速率 $\left(v = \sqrt {v_x^2 + v_y^2 + v_z^2}\right)$ -时间曲线与介质振动位移-时间曲线变化趋势基本吻合，并且三轴合成速率与介质振动位移之间的相关性优于y轴振动速率^[15]。故将三轴合成速率作为爆破振动评判标准更为合理。通过数值模拟得到的振动波三轴合成速率-时间云图如图3所示。

图 3 三轴合成速率-时间云图

Figure 3. Nephograms of synthetic velocity

下载: 全尺寸图片幻灯片

在靠近边坡面位置上等距离选取6个监测点，网格介质序列依次为：A（178 321）、B（175 881）、C（173 217）、D（170 627）、E（168 332）、F（165 890）。通过对该6个监测质点的合成速度-时间曲线、振动位移-时间曲线、合成加速度-时间曲线进行分析，得到以下结论。（1） 0～4 ms内，所选介质的合成速度迅速增大，且变化率最大，说明加速度数值较大；通过观测振动波合成速度-时间曲线与介质合成加速度-时间曲线得到一致结论；该时间段介质振动位移数值较小，变化率较大，观测结果符合理论实际。其原因是：在该时间段内冲击波主要对岩体起动力粉碎作用，起爆初期冲击波使岩体产生微小裂纹；随着冲击波的集聚、传播、转化为应力波并进入微裂隙，应力波反射拉伸作用将微裂隙进一步扩展，形成炮孔粉碎区和显著的裂纹裂隙区。（2）在4～10 ms时间段内，合成速率继续增大，且增长速率减缓，加速度逐渐降低，振动位移不断增大，且增长速率逐渐变小；其主要原因为在该时间段内应力波基本消失，爆轰气体压力超过岩石断裂强度，爆轰气体进入裂纹裂隙对岩体进行二次拉裂破坏。（3） 10 ms以后，合成速度基本趋于稳定，达到最大值7.9 km/s，加速度趋于零，振动位移呈线性增长；原因为该时间段内炸药能量释放结束，炮孔内集聚的爆轰气体和破碎岩块运动过程中相互挤压碰撞作用促使岩体进一步破裂分离。合成速度-时间曲线、振动位移-时间曲线、合成加速度-时间曲线分别如图4、图5和图6所示。

图 4 合成速度-时间曲线

Figure 4. Resultant velocity vs. time

下载: 全尺寸图片幻灯片

图 5 振动位移-时间曲线

Figure 5. Vibration displacement vs. time

下载: 全尺寸图片幻灯片

图 6 合成加速度-时间曲线

Figure 6. Resultant acceleration vs. time

下载: 全尺寸图片幻灯片

3.2 裂纹裂隙扩展分析

裂纹裂隙扩展主要是爆破初期冲击波（应力波）的动力破坏和爆轰气体准静态破坏综合作用的结果。起爆瞬间冲击波的动态应力使炮孔周围岩层产生径向和环向微裂隙；随着冲击波传播能量递减，部分冲击波逐渐转化为应力波，微裂隙在后期传播过来的冲击波以及已经转换成的应力波的叠加作用下再次受到反射拉伸破坏，加剧裂隙进一步扩展；随着时间推移，冲击波（应力波）能量减小，而爆轰气体压力集聚超过岩体的断裂强度，其高能高压作用对微裂隙岩块施加准静态破坏，使得岩体形成贯穿性裂隙。

分析裂纹裂隙扩展过程的数值模拟结果可知：应力波最早于0.6 ms到达坡脚位置，此时波速v=5.78 km/s；在0.9 ms左右，炮孔孔底开始产生裂纹；随着应力波继续传播，边坡自由面发生反射拉伸破坏，直至1.2 ms左右坡脚处产生贯穿性裂隙；在2.8 ms左右，应力波传至坡顶位置，扩腔现象表现明显，此时粉碎区基本形成，岩体破裂呈块状趋势，不再发生粉碎现象；6.1 ms左右，裂纹裂隙首次贯穿整个岩体和自由面，直至12.5 ms左右，裂纹裂隙扩展基本完成。裂纹裂隙的动态扩展过程如图7所示。

图 7 不同时间裂纹裂隙扩展形态分布

Figure 7. Variety of crack development of rock at different time

下载: 全尺寸图片幻灯片

3.3 爆轰压力分布

同样对边坡面附近6个监测点能量和所受压力-时间分布结果进行分析，发现：整个爆破过程中，边坡自由面附近所受压力的最大值为1.36 GPa，并且在1～4 ms时间段内压力峰值集中出现在药包周围介质点；平行于装药段边坡区域的峰值压力明显高于平行于堵塞段边坡区域的峰值压力，平行于装药段边坡区域压力分布在0.70～1.36 GPa范围，平行于堵塞段边坡区域的压力分布在0.25～0.55 GPa范围，说明炮孔装药段位置应该最早发生岩块分离的抛掷现象；峰值压力过后，边坡面所受压力急剧下降，在6 ms左右降至最小，约52 MPa，此时炮孔围岩未完全断裂形成块体分离的抛掷运动；7.5 ms左右冲击波（应力波）基于整个边坡区域的动力破坏结束，炮孔粉碎区早已形成，半径约为28 cm，约等于装药半径的4倍。而冷振东等^[16]研究发现炮孔粉碎区半径约为装药半径的1.2～5.0倍，说明本研究所得结果与该理论结果吻合。6个监测点的压力-时间曲线如图8所示。

图 8 边坡压力-时间曲线

Figure 8. Slope pressure vs. time

下载: 全尺寸图片幻灯片

3.4 抛掷运动分析

随着爆轰气体的集聚，当气体压力超过岩石裂隙尖端断裂因子强度时，贯穿裂隙在气体压力作用下产生沿边坡自由面方向的抛掷现象。不同时刻岩体的抛掷形态分布如图9所示。

图 9 起爆后不同时刻岩体的抛掷形态

Figure 9. Variety of damage and deformation of rock at different time

下载: 全尺寸图片幻灯片

黄永辉等^[4]通过高速摄影机对工程爆破抛掷现象进行拍摄，通过试验获得了抛掷时间、块体初速度以及抛掷运动轨迹等，然而由于爆破过程中爆轰波引起的振动和大量粉尘影响拍照，加之抛掷过程中岩块之间的相互挤压和碰撞造成二次破碎现象，增大了物理量测量的误差范围；本研究通过数值模拟可以直观地动态反演整个抛掷现象。

分析抛掷过程模拟结果可得：最早于7 ms左右，坡脚最先出现抛掷速度4.21 m/s；随着时间的推移，抛掷现象由坡脚向自由面演化，通过模型抛掷单元监测发现，最大抛掷速度出现在炮孔药包中间位置，为21.6 m/s，该值与黄永辉等^[4]提出的最大抛掷速度为18～28 m/s的结论一致；边坡自由面抛掷速度小于炮孔周围的岩块抛掷速度，内外岩块之间存在速度差，正好解释了二次挤压碰撞的破碎现象；抛掷块状分离现象从炮孔药包的中间部位开始，并且块体最大抛掷速度集中在该部位至边坡自由面之间的垂直区域。

3.5 岩体破碎块度分析

基于岩体破碎理论，材料单一的均质岩体在明确的失效和损伤准则规范下，爆破形成的破碎岩块当属均质岩块，块体尺寸分布在0.1～0.5 m。而实际爆破工程中由于岩层结构、原始裂纹裂隙发育状况、装药量、起爆方式、最小抵抗线等因素不同而导致岩体破碎块度大小分布各异。本研究通过数值模拟分析爆破过程中的大块分布，故忽略岩层结构对块度的影响。

分析发现，大块主要来源于边坡坡脚、炸药与堵塞物接触面两侧以及台阶顶部自由面处，如图10所示。坡底出现大块主要取决于底盘抵抗线数值的选取，抵抗线取值不同导致块度分布不同，如图10中5号岩块所示，直径约为2.3 m。炸药与堵塞物接触面两侧出现大块主要是由于该接触面两侧材料的力学性质和物理参数不同，应力波在该接触面发生反射和透射现象，反射拉伸作用于贯穿性微裂隙，使该处岩体在拉应力作用下脱离母岩形成大块，如图10中1、2号和3、4号岩块所示，直径分布在1.6～2.7 m范围内。台阶顶部自由面出现大块是由于炸药和堵塞物接触面的应力波透射传至该顶部自由面，透射波同样形成反射拉伸作用，使该区域岩体沿内部微裂隙断裂破碎形成大块，如图10中6号岩体所示，其直径约为1.7 m。

图 10 大块分布形态

Figure 10. Distribution of rock blocks

下载: 全尺寸图片幻灯片

4. 结　论

（1）假定岩体在冲击载荷作用下的受压状态服从Mises屈服破坏，且应力-应变关系满足双线性弹塑性本构关系，则可构建压应力作用下硬化效应的本构损伤模型；而在受拉状态下满足TCK关系，且材料响应服从弹性胡克定律，从而构建拉裂损伤计算模型；将存在宏细观缺陷的岩体拉裂损伤计算模型和压应力硬化效应的本构关系耦合，能够更加准确地反映岩体在冲击载荷作用下的破裂情况。这样既考虑了岩体受压状态下的硬化效应，又合理地解释了冲击载荷作用下岩体的弹塑性变形规律。

（2）模拟发现振动波三轴合成速率曲线与介质振动位移曲线变化对应的时间节点和步长基本吻合。结合文献[14]分析实测振动波三轴合成速率峰值-位移分布特征，可得三轴振动速率与振动位移之间的相关性优于垂直方向（y轴）振动速率与对应振动位移之间的相关性，故三轴振动速率可作为降振减灾的判定指标。

（3）岩质边坡最早于0.6 ms左右从坡脚位置产生裂纹，0.9 ms抵达炮孔孔底，2.8 ms时应力波传至坡顶位置形成粉碎区，半径为28 cm，约等于装药半径的4倍，6.1 ms时整个岩体和自由面形成贯穿裂隙，直至12.5 ms左右裂纹裂隙扩展完成，1～4 ms时边坡面质点集中出现峰值压力，最大值达到1.36 GPa；炮孔装药段位置最早发生块状分离的抛掷现象，最大抛掷速度集中于该部位至边坡自由面之间的垂直区域，最大值达到21.6 m/s；破碎大块主要来源于边坡坡脚、炸药与堵塞物接触面两侧围岩以及台阶顶部自由面处，大块块体直径分布在1.6～2.7 m范围内。

（4）爆炸冲击载荷作用下岩体不仅存在拉压损伤和硬化效应，还受到塑性压剪损伤以及应变软化和损伤软化效应，本研究所构建的损伤模型未计入塑性压剪损伤以及应变软化和损伤软化效应，相关后续工作将另文探讨。

图 1 损伤模型嵌入流程

Figure 1. Flowchart for damage model

下载: 全尺寸图片幻灯片

图 2 台阶平面模型

Figure 2. Model of the bench plane

下载: 全尺寸图片幻灯片

图 3 三轴合成速率-时间云图

Figure 3. Nephograms of synthetic velocity

下载: 全尺寸图片幻灯片

图 4 合成速度-时间曲线

Figure 4. Resultant velocity vs. time

下载: 全尺寸图片幻灯片

图 5 振动位移-时间曲线

Figure 5. Vibration displacement vs. time

下载: 全尺寸图片幻灯片

图 6 合成加速度-时间曲线

Figure 6. Resultant acceleration vs. time

下载: 全尺寸图片幻灯片

图 7 不同时间裂纹裂隙扩展形态分布

Figure 7. Variety of crack development of rock at different time

下载: 全尺寸图片幻灯片

图 8 边坡压力-时间曲线

Figure 8. Slope pressure vs. time

下载: 全尺寸图片幻灯片

图 9 起爆后不同时刻岩体的抛掷形态

Figure 9. Variety of damage and deformation of rock at different time

下载: 全尺寸图片幻灯片

图 10 大块分布形态

Figure 10. Distribution of rock blocks

下载: 全尺寸图片幻灯片

表 1 岩石主要物理力学参数

Table 1. Physical and mechanical parameters of rock

$\rho$ /(g·cm^–3)	E/GPa	$\nu$	${\sigma _0}$ /MPa	$\beta$	G/GPa	K/GPa	${K_{{\rm{IC}}}}$ /(10⁵ N·m^-3/2)
2.7	18.23	0.23	102	1.0	7.41	11.25	5.32

下载: 导出CSV

表 2 炸药主要参数

Table 2. Parameters of explosive

$\rho$ /(g∙cm^-3)	p/GPa	A/GPa	B/GPa	${R_1}$	${R_2}$	$\omega$	E₀/GPa	${V_0}$
1.2	27	52.1	0.182	4.2	0.9	1.1	4.192	1.0

下载: 导出CSV

参考文献(16)

[1]	BUDIANSKY B, O'CONNELL R J. Elastic moduli of a cracked solid [J]. International Journal of Solids and Structures, 1976, 12(2): 81–97. doi: 10.1016/0020-7683(76)90044-5
[2]	朱传云, 戴晨, 姜清辉. DDA方法在台阶爆破仿真模拟中的应用 [J]. 岩石力学与工程学报, 2002, 21(Suppl 2): 2461–2464 ZHU C Y, DAI C, JIANG Q H. Numerical simulation of bench blasting by discontinuous deformation analysis method [J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(Suppl 2): 2461–2464
[3]	MUNJIZA A. The combined finite-discrete element method [M]. London: John Wiley & Sons Ltd, 2004: 29-32.
[4]	黄永辉, 刘殿书, 李胜林, 等. 高台阶抛掷爆破速度规律的数值模拟 [J]. 爆炸与冲击, 2014, 34(4): 495–500 doi: 10.11883/1001-1455(2014)04-0495-06 HUANG Y H, LIU D S, LI S L, et al. Numerical simulation on pin-point blasting of sloping surface [J]. Explosion and Shock Waves, 2014, 34(4): 495–500 doi: 10.11883/1001-1455(2014)04-0495-06
[5]	TAYLOR L M, CHEN E P, KUSZMAUL J S. Microcrack-induced damage accumulation in brittle rock under dynamic loading [J]. Computer Methods in Applied Mechanics and Engineering, 1986, 55(3): 301–320. doi: 10.1016/0045-7825(86)90057-5
[6]	KIPP M E, GRADY D E. Numerical studies of rock fragmentation [M]. Albuquerque: Sandia National Laboratories, 1978.
[7]	Livermore Software Technology Corporation. LS-DYNA keyword user's manual [M]. California: Livermore Software Technology Corporation, 2003.
[8]	MA G W, AN X M. Numerical simulation of blasting-induced rock fractures [J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(6): 966–975. doi: 10.1016/j.ijrmms.2007.12.002
[9]	王肖钧, 刘文韬, 张刚明, 等. 爆炸载荷下钢板层裂的二维数值计算 [J]. 爆炸与冲击, 1999, 19(2): 97–102 doi: 10.3321/j.issn:1001-1455.1999.02.001 WANG X J, LIU W T, ZHANG G M, et al. 2D numerical simulation of spallation in a steel plate due to explosive loading [J]. Explosion and Shock Waves, 1999, 19(2): 97–102 doi: 10.3321/j.issn:1001-1455.1999.02.001
[10]	赵铮, 陶钢, 杜长星. 爆轰产物JWL状态方程应用研究 [J]. 高压物理学报, 2009, 23(4): 277–282 doi: 10.3969/j.issn.1000-5773.2009.04.007 ZHAO Z, TAO G, DU C X. Application research on JWL equation of state of detonation products [J]. Chinese Journal of High Pressure Physics, 2009, 23(4): 277–282 doi: 10.3969/j.issn.1000-5773.2009.04.007
[11]	李夕兵, 古德生, 赖海辉. 冲击载荷下岩石动态应力-应变全图测试中的合理加载波形 [J]. 爆炸与冲击, 1993, 13(2): 125–130 LI X B, GU D S, LAI H H. On the reasonable loading stress waveforms determined by dynamic stress-strain curves of rocks by SHPB [J]. Explosion and Shock Waves, 1993, 13(2): 125–130
[12]	吕涛, 石永强, 黄诚, 等. 非线性回归法求解爆破振动速度衰减公式参数 [J]. 岩土力学, 2007, 28(9): 1871–1878 doi: 10.3969/j.issn.1000-7598.2007.09.019 LÜ T, SHI Y Q, HUANG C, et al. Study on attenuation parameters of blasting vibration by nonlinear regression analysis [J]. Rock and Soil Mechanics, 2007, 28(9): 1871–1878 doi: 10.3969/j.issn.1000-7598.2007.09.019
[13]	许名标, 彭德红. 某水电站边坡开挖爆破震动动力响应有限元分析 [J]. 岩土工程学报, 2006, 28(6): 770–775 doi: 10.3321/j.issn:1000-4548.2006.06.018 XU M B, PENG D H. Finite element analysis of dynamic response on blasting vibration in slope excavation of a hydroelectric power station [J]. Chinese Journal of Geotechnical Engineering, 2006, 28(6): 770–775 doi: 10.3321/j.issn:1000-4548.2006.06.018
[14]	傅洪贤, 赵勇, 谢晋水, 等. 隧道爆破近区爆破振动测试研究 [J]. 岩石力学与工程学报, 2011, 30(2): 335–340 FU H X, ZHAO Y, XIE J S, et al. Study of blasting vibration test of area near tunnel blasting source [J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(2): 335–340
[15]	楼晓明, 周文海, 简文彬, 等. 微差爆破振动波速度峰值-位移分布特征的延时控制 [J]. 爆炸与冲击, 2016, 36(6): 839–846 doi: 10.11883/1001-1455(2016)06-0839-08 LOU X M, ZHOU W H, JIAN W B, et al. Control of delay time characterized by distribution of peak velocity-displacement vibration of millisecond blasting [J]. Explosion and Shock Waves, 2016, 36(6): 839–846 doi: 10.11883/1001-1455(2016)06-0839-08
[16]	冷振东, 卢文波, 陈明, 等. 岩石钻孔爆破粉碎区计算模型的改进 [J]. 爆炸与冲击, 2015, 35(1): 101–107 doi: 10.11883/1001-1455(2015)01-0101-07 LENG Z D, LU W B, CHEN M, et al. Improved calculation model for the size of crushed zone around blasthole [J]. Explosion and Shock Waves, 2015, 35(1): 101–107 doi: 10.11883/1001-1455(2015)01-0101-07