摘 要：岩体爆炸是一个复杂的瞬时过程，确定爆炸荷载更是一个难点。本文应用FLAC3D软件，模拟了两种不同形式的爆破荷载，并将计算结果与现场监测结果进行了比较。数模拟结果表明，不同的爆破荷载施加方法对应的数值模拟结果在近区存在一定差异，但在中远区吻合得较好。最后给出了两者的适用范围。

　　关键词：爆破荷载;爆破震动;岩质边坡;数值模拟

　　Abstract: The rock explosion is a complex transient process，determine the explosion load was a more difficult point. In this paper, use FLAC3D simulate two different forms of blast loading, and monitoring the results with field results. The Imposed different methods of blastingt load，the corresponding numerical simulation results show that in the near field there are some differences, but coincide better in further area. Finally，give the application scope of both methods.

　　Key words: blast loading; blasting vibration; rock slope; numerical simulation 1 引言

　　动力输入荷载的确定是分析爆破荷载作用下岩质边坡动态响应的基础，也一直是难点之一。岩体爆炸是一个复杂的瞬时过程，爆破荷载作用下岩石边坡的响应特征研究主要集中在现场监测和数值模拟2个方面。数值模拟的目的是研究距爆源中远距离处岩体质点的振动情况，以便对爆破荷载作

　　用下岩质边坡的速度、位移和应力状态等动

　　态响应特征进行正确的描述。爆破荷载计算方法有DDA法，离散元法，遗传算法，复合遗传算法和体积平衡法。由以上的爆破荷载计算方法,可以进行爆破模拟。目前爆破震动传播的数值模拟通常采用三种方法:(1)文献[1 2]在炮孔壁上直接施加半理论半经验的爆破荷载压力曲线(如三角型脉冲荷载或指数衰减型荷载等)(以下称模拟1);(2)文献[3 4]在弹性震动区，爆破动力荷载直接采用现场实测的爆破振动速度时程曲线输入(以下称模拟2);(3)文献[5 6]使用LS-DYNA模块中内建的高性能炸药材料及其状态方程来模拟炸药的整个爆轰过程及爆轰产物与周围岩石介质的相互作用过程。由于方法(3)只能较好地模拟爆源近区的冲击响应问题，若用来研究中远区的爆破震动传播问题，尚存在一定缺陷[7]。因此，本文应用FLAC3D程序模拟前两种爆破荷载作用下岩质边坡的动态响应，分析爆破荷载作用下边坡质点的振动速度，并与现场监测结果进行了比较。

　　2 建立计算模型

　　2.1 模型网格

　　本次数值模拟的目的在于模拟爆破荷载作用下岩质边坡的动力响应规律，模拟对象为广贺高速某段路堑左侧边坡爆破震动监测的一次实验。本次模拟的是单孔爆破，装药量为25kg。模型中边坡上底宽25 m，下底宽84 m，高58 m，单台阶高10 m。边坡地形较陡峭,节理裂隙发育，边坡地质条件较好。完整岩石材料的力学参数见表1。 表1 岩石物理力学指标 岩层 密 度(kg/m3) 粘聚力(MPa) 内摩擦角(。) 泊松比μ 弹性模量(GPa) 体积模量(GPa) 剪切模量(GPa) 抗拉强度(MPa) 砾岩 2200 0.12 45 0.25 9.12 6.08 3.65 0.67 在有限元计算中，网格单元的尺寸对计算结果有很大影响。一般而言，网格单元的尺寸越小，计算精度越高，但这样会占用更多的内存和消耗更多的计算时间。在FLAC3D模拟中，一般将块体切割成大量的有限差分元六面体网格单元。为精确模拟波在模型中的传播，设置网格尺寸时，应使得网格尺寸Δl满足下式[8]：

(1)

式中,为最高频率对应的波长。

　　在FLAC3D动力计算中，其时步需要遍历模型中所有的单元。时步是由几何尺寸最小、模量最大的单元确定的。要减少计算时间和获得较为精确的计算结果，就要尽量避免较小的单元，同时保证单元尺寸满足式

　　(1)。模型网格如图1所示。

　　图1有限元模型

　　2.2 边界条件的确定

　　进行岩土工程问题的模拟时，相对于分析尺寸而言，进行模型计算时，所涉及的介质在分析范围内的表述是无边界的。由于模型是有限的，人工提取的边界也是有限的。在动力分析中，由于计算模型的区域边界有可能造成外传波的反射，给数值分析结果的准确性造成一定影响。为了解决这个问题，本文采用了文[9]提出的粘滞(不反射)边界。所谓无反射边界是一种人工边界，它是通过在边界上施加切向和法向的阻尼筒来实现。边坡的左右侧面均设为粘性(不反射)边界，其x和z方向均不施加约束。边坡的下底也设为粘性边界，且其底部z向固定。根据实际装药情况，均布爆炸荷载施加于第一级边坡的炮孔壁上(图2)。

图2 模型计算简图

　　2.3 输入荷载的确定

　　对于目前常用的三角形脉冲波，它是根据相关的经验公式和质点振动的监测结果推测爆炸冲击波的大小及衰减情况。该荷载形式需要确定荷载峰值和爆破荷载作用时间。当采用的装药形式为不耦合装药时，爆破荷载峰值由下式确定[10]。

(2)

对于不耦合装药，根据研究，爆轰产物撞击药室壁时，孔壁压力将明显增大，增大倍数n=8-11。因此，可以得到柱状不耦合装药时，药室受到的冲击压力表达式为 (3)

　　式中，P0为爆破荷载峰值压力;ρe为炸药密度，在现场使用的是2号岩石乳化炸药，密度为1300kg/m3;D为炸药爆轰速度,取4000m/s;γ为炸药的等熵指数，与装药密度有关;dc为装药直径，为50mm;db为炮孔直径，为90mm;lc为装药长度，lb为炮孔长度;n为增大倍数，一般取10。

　　对于爆破荷载的作用时间，一般认为，炸药爆炸时的冲击波作用持续时间约为1.0×10-6-0.1s,爆生气体压力作用时间更长，约为1.0×10-3-0.1s。而卢文波[11]与张建华[12]的研究显示，炮孔压力持续时间约为数百微秒。因此可以认为，取几百微秒至毫秒级的持续时间较为合理。本文中，假定爆破荷载作用总时间为7ms，其中爆破荷载压力上升时间为2ms。图3为爆破荷载示意图，图中横轴为时间(ms)，纵轴为作用在炮孔壁上的峰值压力400(MPa)。

　　图3 药室爆破压力曲线

　　在FLAC3D模拟中，动载荷的输入可以采用加速度时程，速度时程，位移时程和应力(力)时程四种方式。FLAC3D通过设置动态多步可以保持数值的稳定性。需要注意的是，以速度时程曲线表示的动力荷载不能直接加在粘性边界面上，故须运用下述公式将速度曲线转化为应力曲线然后加在边界上[13]。而对于三角型爆破冲击荷载可以施加在炮孔壁的法向。粘滞边界通过在边界的法线方向和水平方向上设置独立的黏壶以便吸收来自模型内部的入射波。法向黏滞力fn和剪切黏滞力fs的计算公式为

(4)

(5)

　　式中：vn、vs分别为边界上的法向和切向速度分量;ρ为密度;Cp、Cs分别为P波和S波的波速。

　　图4现场实测水平径向与垂直向振动速度时程曲线

　　2.4 阻尼的选取

　　进行动力计算时，FLAC3D提供了3种可供选择的阻尼形式，分别为瑞利阻尼、局部阻尼和滞后阻尼。由于瑞利阻尼的理论与常规动力分析方法类似，故本文阻尼形式选取瑞利阻尼。设置瑞利阻尼时，需要确定两个参数，分别为最小临界阻尼比ξmin和最小中心频率fmin。因为数值计算中精确地模拟岩体振动的阻尼效应是很困难的，根据试算，本次模拟采用的最小临界阻尼比为0.03，最小中心频率为10Hz。

　　2.5 本构模型的选取

　　FLAC3D的非线性动力分析计算可以遵循任何指定的非线性本构模型。其程序内置10多种岩土本构模型。

　　莫尔-库伦模型和德鲁克-普拉格模型是最常用的岩土本构关系。前者适用于剪切应力作用下达到屈服状态的材料。后者适用于摩擦角较小的粘性土。

FLAC3D计算采用的是弹性增量定律。设主应力增量为，，(，拉应力为正，压应力为负)，对应的弹性应变增量为，，，其本构关系表示如下：

(6)

式中和分别为材料的体积模量和剪切模量。

　　模拟时选取的莫尔-库伦弹塑性模型。模型的破坏包络线和莫尔-库伦强度准则(剪切屈服函数)以及拉破坏准则(拉屈服函数)相对应。图5为莫尔-库伦强度准则的示意图。

图5 Mohr-Coulumb强度准则示意图

　　3 结果分析

　　3.1 振动速度峰值

　　图6～9分别为距爆源12m和58m处的水平径向和垂直方向FLAC3D模拟计算得到的质点振动速度与现场监测结果的比较。

　　图6距爆源12 m质点水平振动速度

　　图7距爆源12 m质点垂直振动速度

　　图8距爆源58 m质点水平振动速度

　　图9距爆源58 m质点垂直振动速度

　　图6-9分别为由FLAC3D模拟得到的垂直向与水平径向的振动速度时程曲线与现场实测结果的对比，从图中可以看出：

　　(1)本文简化了爆炸能量的分布，但由于其它参数的取值比较合理，因此模拟结果与实测结果吻合得比较好。在距爆源较近的区域爆破地震波急剧衰减，到远区趋于平缓。由FLAC3D模拟得到的垂直向和水平向振动速度峰值，与监测得到的峰值的误差均在工程允许误差范围内。

　　(2)从图中可以它们达到各自的峰值的时间是不一样的。距离越远，达到峰值所需要的时间就越长。即使是同一个点，达到水平径向峰值与垂直向峰值的时间也是不一样的。

　　(3)在近区，采用三角波加载所得的上升速度比实测曲线快。因此，采用半经验公式的压力曲线比现场实测的压力曲线更早达到振速峰值。在远区，由于实测曲线的应力衰减很慢，因此，在远区由半经验公式的压力引起的速度比由现场实测的压力引起的速度要大。

　　3.2 衰减特性分析

图10竖直向衰减规律

图11水平径向衰减规律

　　从图10-11可以看出

　　(1)随着爆心距的增大，不论是哪种载荷施加方式，岩坡的爆破震动响应速度逐渐衰减。后面出现了振速衰减变缓的现象，这可以看作是出现了不同程度的高程效应。

　　(2)数值模拟得到的质点振速普遍大于现场实测振速。这可能是由于数值模拟建模时没有考虑到边坡的节理裂隙的影响。

　　通过以上分析，可以得到两种方法各自的适用范围：

　　(1)如果要考虑应力波对爆破近区的作用，就需要采用半经验的压力公式。它能模拟爆破作用下，爆破近区的变化。并且要现在爆破现场测得爆破近区的振动情况是很相当困难的。

　　(2)由于应力波很快会衰减成地震波，而地震波只会导致裂隙的扩展，不会对岩石造成破坏。如果只是考虑振速对边坡的影响，两种方法均适用。

　　4 结论

　　本文以广贺高速某段路堑边坡爆破工程为背景，通过施加两种不同的等效爆破荷载进行爆破震动数值模拟，并与实测波形进行对比。在对振速衰减规律及进行比较之后，发现两种荷载所得结果在近区的差异较大，但在量级上没有差别。在中远区差别较小，而且随着爆心距的增加，差异也越来越小。总的来说，对于爆破近区，半经验公式所确定加载方法比较合适;而中远区，则更适合采用近区所得的速度曲线进行加载的方法。而关于爆破振动问题的研究，考察范围着重于中远区。因此, 实际中可以根据近区所得的速度曲线进行数值模拟，模拟结果的误差在工程应用的允许范围内，可以用来指导爆破施工。

　　参考文献

　　[1]刘国华,王振宇.爆破荷载作用下隧道的动态响应和抗爆分析[J].浙江大学学报(工学版), 2004, 38(2):15～18.

　　[2]夏祥，李俊如，李海波等。爆破荷载作用下岩体振动特征的数值模拟[J].岩土力学，2005，26(1)，50～56

　　[3]刘亚群,李海波.爆破荷载作用下黄麦岭磷矿岩质边坡动态响应的UDEC模拟研究[J].岩石力学与工程学报, 2004, 23(11): 42～45.

　　[4]陈占军,朱传云,周小恒.爆破荷载作用下岩石边坡动态响应的FLAC3D模拟研究[J].爆破，2005，22(4)，8～13

　　[5]Livermore Software Technology Corporation. LS-DYNA keyword user’s manual [R]. Version 970.2003.

　　[6]赵海鸥. LS-DYNA动力分析指南[M].北京:兵器工业出版社,2003.

　　[7]丁桦,郑哲敏.爆破震动等效载荷模型[J].中国科学(E辑),2003,33(1):82～90.

　　[8]祁生文,伍法权,严福章等著.岩质边坡动力反应分析，科学出版社[M]，2007

　　[9] Lysmer J，Kublmeyer R L.Finite dynamics model for infinitemedia[J].J.Eng.Mech.，&Foundations Div.，ASCE，1969，95(4)：859～877

　　[10]戴 俊，王树仁.爆破工程[M].机械工业出版社，2005

　　[11]卢文波，陶振宇.预裂爆破中炮孔压力变化历程的理论分析[J].爆炸与冲击，1994, 14( 2) : 140- 147.

　　[12]张建华，李世禄，王玉杰等.爆炸扩腔数值模拟及分析[J].武汉科技大学学报(自然科学版)，2001, 24(2):174-177

　　[13]陈育明，徐鼎平.FLAC/FLAC3D基础与工程实例[M]，中国水利水电出版社，2008