基于非连续变形分析方法的爆破荷载下不同类别岩体的动力响应

江 巍，程圣国，冯 强

(三峡大学三峡库区地质灾害教育部重点实验室，湖北 宜昌 443002)

爆破是工程岩体开挖的常见手段，20世纪60年代以来爆破技术已经得到了长足的发展，但爆破带来的不利影响问题(如洞室安全[1]、边坡稳定[2]、近区民居破坏[3]等)也愈发明显。随着爆破理论和计算机技术的发展，数值模拟已经成为分析爆破产生不利影响的重要手段，动力有限元[4]、离散元[5]、有限差分法[2]等多种计算方法被引入该领域并取得了一定的研究成果。

采用连续数值计算方法分析爆炸荷载的效应时将工程岩体视为均一介质，这与岩体中广泛存在裂隙和节理的实际情况并不相符[6]，而且对于爆炸荷载作用过程中，岩体结构面的张开拉裂、块体飞散等动态过程，采用不连续数值计算方法更为合适。非连续变形分析(discontinuous deformation analysis，DDA)方法是石根华先生提出的用于统一分析块体系统静力和动力行为的不连续数值方法[7]，关于其有效性和正确性Mac Laughlin等[8]已经对1986—2006年的相关工作做了全面的综述。该方法在岩土工程领域已经得到了广泛的研究与应用，但是在模拟工程岩体中的爆破问题方面则刚刚起步，如朱传云等[9]进行了台阶爆破的仿真模拟，刘红岩等[10]进行了爆破漏斗形成过程模拟等。

对于工程岩体爆破振动引起的不利影响，主要通过考察工程区岩体的动力响应进行衡量。工程岩体根据岩体质量常划分为不同类别，本文以非连续变形分析方法为工具，对不同类别的岩体在爆破荷载作用下的动力响应进行数值研究，分析岩体振动的相关参数(起振时间、峰值振速)随距离的变化特征，为爆破的控制提供技术依据。

1 爆破荷载的等效

用数值方法模拟爆破后岩体的动力响应，首要前提就是确定爆破荷载。目前爆破荷载模型可以粗略划分为经验公式法的爆破荷载模型和考虑爆生气体状态的爆破荷载模型两类。

经验公式法建立的爆破荷载模型将作用在爆孔壁上的爆破荷载压力简化为随时间变化的函数，基本形式为

式中：Pb为压力峰值;f(t)通常为指数型的时间滞后函数。考虑药包在爆孔中爆炸时，爆孔周围岩体在瞬间高压作用下变成似流体介质，形成粉碎区，该类方法直接将压力作用于爆孔壁上，则如何定义粉碎区的岩石力学特性将成为一个难题。

考虑爆生气体状态的爆破荷载模型将爆孔内的爆生气体压力P 与爆腔体积V 对应起来：

式中的爆腔体积V 包含了爆孔体积和孔壁裂隙扩展产生的体积。李宁等[11-12]均对这一模型做了优秀的研究工作，并取得了丰富的成果，但该方法要求将爆孔周围网格划分得十分细小，在面对大型的计算模型时，会产生过多的计算单元和庞大的计算量。

鉴于上述两类方法的局限性，寻求一种适当的爆破荷载简化等效方法是有必要的，本文采用徐全军等[13]提出的粉碎区边缘施加爆炸压力方法，该方法通过分析水孔法测量得到的靠近炮孔岩体中的压力波形，获取粉碎区边缘的压力曲线，其表达式如下：

式中：Pe为作用在粉碎区边缘的压力峰值，MPa;α为衰减系数，与压力持续时间有关。对于完全耦合装药，作用在粉碎区边缘的压力峰值为

式中：γ为绝热指数;ρ0为炸药初始密度，kg/m3;Qv为单位质量爆热，J/kg;D为爆速，m/s;d0为爆孔直径，m;de为粉碎区直径，m。爆炸荷载的持续时间为

式中：R1为离爆源的距离，m;C为爆生气体膨胀后的音速，取370m/s。以铵油炸药为对象[13]，可得Pe=5.87MPa，T=0.00216 s，α=4500，粉碎区边缘压力曲线如图1 所示。

图1 粉碎区边缘压力曲线

该方法的基本思想为：当药包爆炸后，以爆孔为中心周围的岩体按照由近至远形成粉碎区、破碎区和振动区，将爆破荷载直接施加在粉碎区的边缘，则可以避开粉碎区的岩石力学特性定义难题。基于此方法徐全军等[13]采用非线性有限元程序对爆破振动效应进行了分析，取得了一定的成果，但在计算中为使用有限元技术将破碎区和振动区的岩体均近似为均匀连续介质。本文以非连续变形分析方法为手段，视破碎区和振动区的岩体由线弹性体的岩块和节理组成，破碎区岩体中除岩块发生弹性变形外更主要的是节理新生、扩张和延伸等导致岩体发生破碎行为，振动区的岩体中节理受力过小变形细微，两侧岩块之间力的传递良好，其力学行为表现接近连续介质。

2 非连续变形分析基本理论与爆破荷载的施加

2.1 非连续变形分析方法基本理论简介

非连续变形分析方法中的基本单元是块体。块体的位移和变形是由分时步的位移和变形累加而成的，在每一时步所有单元的位移都是小量。考虑每个块体为常应力和常应变，块体中一点(x，y)的位移(u，v)可用位移不变量表示为

式中：(u0，v0)为块体重心(x0，y0)的刚体平动;r0为块体绕点(x0，y0)的转动角;εx，εy，γxy为该块体的法向和切向应变。块体单元任意点的变形均可由上述位移不变量乘形函数矩阵得到。

块体系统方程组的建立根据系统势能最小化原理获得，系统的总势能包括块体单元的应变能、初始应力的势能、点荷载和线荷载的势能、体荷载势能、锚杆作用力的势能、惯性力势能以及块体相互接触力的势能。块体系统总势能Π的一般表达形式

式中：D，K，F 分别为整个块体系统的位移不变量、总刚度矩阵和总荷载向量。 Π取极值的条件为其对D的每一分量求偏导均为零，于是可得支配方程：

施加边界条件和块体系统的运动学条件，则可对方程(8)进行求解，得到每一块体的位移和变形状态。

2.2 爆破荷载的施加及边界条件

在第1 节中爆破荷载已经被等效为随时间变化的作用于粉碎区边缘的压力，因此在非连续变形分析计算模型中该荷载为一分布力，考虑非连续变形分析块体单元几何形状一般为多边形，则该分布力体现为非连续变形分析方法中的线荷载形式。因此，爆破荷载的作用能够非常简便地在程序中得到体现，在非连续变形分析计算过程中的每一时步开始时，将该时刻对应的粉碎区边缘压力施加于边缘块体即可。

在模拟爆破荷载在岩体中的动力响应时，必然需要采用有限区域代替无限区域或者半无限区域，这样建立的模型边界需要进行一定的处理。经典非连续变形分析中对于边界通常设置为固定边界，其不发生任何变形和移动，这样的边界对于应力波起着完全反射的作用，因此必须对其做一定的修改。借鉴FLAC3D动力分析中的黏滞边界来实现应力波的无反射，该边界可以吸收有限域向外界传递的能量，具体操作为在边界块体的沿边界侧两个端点分别施加独立的法线方向黏滞力f n和切线方向黏滞力fs：

式中：ρ为块体密度，kg/m3;cp和cs分别为质点的纵波和横波波速，m/s;vn和vs分别为该端点的法向速度和切向速度，m/s。该黏滞力体现为非连续变形分析方法中的点荷载形式，很容易添加至方程(8)中。

3 不同类别岩体在爆破荷载下的动力响应

在实际工程中，通常根据一定的标准将岩体进行质量分级，GB 50218—1994《工程岩体分级标准》[14]在我国应用广泛，该标准先考虑岩石的坚硬程度和岩体的完整程度进行基本质量分级：

式中：BQ为岩体质量基本指标;Rc为岩块饱和单轴抗压强度，MPa;Kv为岩体完整性系数。然后再根据地下水、岩体软弱结构面和地应力的情况进行修正，岩体按照计算得到的BQ 值共分为5 级：大于550时为Ⅰ级;450和550 之间为Ⅱ级;350和450 之间为Ⅲ级;250和350 之间为Ⅳ级;小于250 时为Ⅴ级。本文将该问题进行适当简化，仅考虑由于岩石坚硬程度和完整程度变化而导致的不同级别岩体。由于实际工程中，对于动力响应主要关注点为爆炸荷载在岩体中的传递时间和造成的不利振动效应，因此在计算中着重考察距离爆孔一定距离处岩石的起始振动时间及其峰值振速并进行分析。

3.1 坚硬程度差异导致的不同类别岩石的动力响应

采用同样结构面分布的厚层状大理岩[15]和厚层状细粒石英砂岩[16]作为样本，两种岩体的相关数据如表1 所示，建立的计算模型如图2 所示。

表1 坚硬程度不同的岩石样本质量分级参数

图2 计算模型

模型长、高均为25m，上下左右四侧边界均为无反射的黏滞边界。模型中的节理分为天然节理和虚拟节理两种，设置虚拟节理后岩石在爆破荷载下的断裂行为能够被模拟，另外由于非连续变形分析块体单元均为常应力常应变单元，虚拟节理的添加可以更加精确地分析块体的应力场和速度场。天然节理不考虑其抗拉强度，内摩擦角取45°，黏聚力取0.05MPa。虚拟节理在计算过程中达到破坏准则后转变为真实节理，其破坏有拉伸破坏和剪切破坏两种情况。当虚拟节理面上的拉应力达到节理面抗拉强度值时，拉伸破坏发生;当虚拟节理面上的剪应力达到摩尔-库伦准则容许的最大剪应力时，剪切破坏发生。虚拟节理的内摩擦角取45°，黏聚力取2.00MPa，抗拉强度取0.5MPa。

计算时采用的两种岩石块体材料参数为：大理岩密度2750kg/m3，弹性模量70GPa，泊松比0.10;石英砂岩密度2650 kg/m3，弹性模量55 GPa，泊松比0.20。将第1 节获得的压力曲线施加于粉碎区边缘，分别对距离爆孔中心3m，5m，10 m和20m 处的块体振动速度进行分析，计算得到的不同坚硬程度的两种岩体各距离处的块体振速随时间变化曲线如图3 所示，为了便于观察，仅取了速度曲线的首个半振幅。

图3 不同坚硬程度岩石振速随时间变化曲线

观察图3可知，随着距离的增加，块体的振速峰值变小，起始振动时间变长，但是两种坚硬程度不同的岩石确存在一定差异。同样距离的情况下，石英砂岩的起始振动时间晚于大理岩，振速峰值小于大理岩。距离爆孔3 m 处，石英砂岩的起始振动时间比大理岩晚0.3 ms，振速峰值比大理岩小0.2 m/s。随着距离的增长，这种趋势愈发明显，距离爆孔5m处，石英砂岩的起始振动时间比大理岩晚0.5 ms，振速峰值比大理岩小0.4m/s。至距离爆孔较远处时，由于此时两种岩体的振速均已较小，所以差别缩小，但起始振动时间的差别仍在继续增大，距离20 m 时石英砂岩的起始振动时间比大理岩晚2 ms。这说明，在同样的结构面情况下，岩石的坚硬程度对爆破荷载的传递影响主要体现为振速峰值与起始振动时间，同样距离情况下，岩石愈坚硬，振速峰值愈大，起始振动时间愈短。将计算结果与文献[16]的结果进行比较，发现坚硬程度对爆破荷载的传递影响趋势基本一致，但是本文计算的起始振动时间与振速峰值均偏低，其原因在于文献[16]中岩体为均质连续介质，而本文中节理的存在在一定程度上阻碍了爆炸荷载的迅速响应。

3.2 完整程度差异导致的不同类别岩石的动力响应

为分析完整程度差异导致的不同类别岩石动力响应，将前面所做的厚层状大理岩的结构面进行改变，以获得一种新的大理岩样本，并将其与原样本的计算结果进行比较分析，两种样本的相关数据见表2，计算模型的区域尺寸和边界条件也保持不变，区别在于大理岩2的结构面分布更加稀疏，间距增大。

表2 完整程度不同的岩石样本质量分级参数

为便于分析比较，仍然采用3.1 节中的大理岩材料参数和节理参数，对距离爆孔中心3m，5m，10m和20 m 处的块体振速进行计算分析，得到的不同完整程度的两种岩体各距离处的振速随时间变化曲线的首个半振幅如图4 所示。

图4 不同完整程度岩石振速随时间变化曲线

根据计算结果，随着距离的增加，块体的振速峰值变小和起始振动时间变长这一基本趋势没有改变，但是结构面的增加对于爆破荷载的传递产生很大的影响。同样距离的情况下，大理岩2的起始振动时间早于大理岩1，块体振速峰值大于大理岩1。距离爆孔3 m 处，大理岩2的起始振动时间比大理岩1 早0.2ms，振速峰值比大理岩1 高0.3 m/s。这种差异也有随距离增长而增大的倾向，距离爆孔5m处，大理岩2的起始振动时间比大理岩1 早0.5 ms，振速峰值比大理岩高0.4 m/s。图4 显示，在同样坚硬的岩石组成情况下，岩石的完整程度也会影响爆破荷载的传递，同样距离时岩石愈完整，振速峰值愈大，起始振动时间愈短。值得关注的是，结构面的增多对于爆破荷载的传递时间影响比岩石坚硬程度更加明显，距离20 m 处，大理岩2的起始振动时间比大理岩1 早3 ms。将计算结果与文献[6]的结果进行比较，爆破荷载的传递受完整程度影响的变化趋势也基本相似，而且在数值上完整程度高的大理岩2的计算结果更接近于采用动力有限元计算得到的结果，这说明在完整程度较优的情况下，采用均质连续介质近似是可以接受的。

值得注意的是，本文采用的节理仍然是一种非常理想化的力学节理，真实的岩体节理是非常复杂的，其特性及分布状态同样也会对岩体中爆破荷载的动力响应造成一定的影响，如朱哲明等[17]针对节理中充填水或者不同材料等情况下应力波传递过程进行的分析等。目前非连续变形分析方法仍处于发展之中，对于节理面的不同充填材料和不同充填度等的处理手段仍然不够成熟，因此对这个问题进行进一步深入的研究是很有价值的。

4 结 论

a.在分析爆破荷载引起的工程岩体中的动力响应问题时非连续变形分析方法能够较好处理岩体的不连续变形行为，但是在处理岩体中多种复杂的节理特性方面还值得进一步研究，如节理中的充填物质性质和充填程度如何反映等。

b.不同坚硬程度和完整程度的岩体对于爆破荷载的动力响应存在明显差异。岩石的坚硬程度愈高，完整程度愈好，岩体对爆破荷载的动力响应愈迅速，能量损耗愈小，振速峰值愈大，完整程度对于爆破荷载传递时间的影响更为明显。

[1]高涣涣.爆破荷载对领近洞室的及支护结构的的影响机理分析[D].西安：西安理工大学，2008.

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

[3]史秀志，周健，崔松，等.露天采矿爆破振动对民房危害预测的DDA 模型及应用[J].中南大学学报：自然科学版，2011，42(2)：441-448.

[4]许红涛，卢文波，周小恒.爆破振动场动力有限元模拟中爆破荷载的等效施加方法[J].武汉大学学报：工学版，2008，41(1)：67-71.

[5]刘亚群，李海波，李俊如，等.爆破荷载作用下黄麦岭磷矿岩质边坡动态响应的UDEC 模拟研究[J].岩石力学与工程学报，2004，23(21)：3659-3663.

[6]张秀丽，焦玉勇，刘泉声，等.节理对爆炸波影响的数值研究[J].岩土力学，2008，29(3)：712-721.

[7]SHI Gen-hua.Discontinuous deformation analysis：a new numerical model for the statics and dynamics of block systems[D].Berkeley，California：Department of Civil Engineering University of California，1988.

[8]Mac LAUGHLIN M M，DOOLIN D M.Review of validation of the discontinuous deformation analysis(DDA)method[J].International Journal for Numerical and AnalyticalMethodsin Geomechanics，2006，30(4)：271-305.

[9]朱传云，戴晨，姜清辉.DDA 方法在台阶爆破仿真模拟中的应用[J].岩石力学与工程学报，2002，21(增刊2)：2461-2464.

[10]刘红岩，杨军，陈鹏万.爆破漏斗形成过程的DDA 模拟分析[J].工程爆破，2004，10(2)：17-20.

[11]李宁，SWOBODA G.爆破荷载的数值模拟与应用[J].岩石力学与工程学报，1994，13(4)：357-364.

[12]甯尤军，杨军，陈鹏万.节理岩体爆破的DDA 方法模拟[J].岩土力学，2010，31(7)：2259-2263.

[13]徐全军，龙源，张庆明，等.微差爆破震动叠加起始位置数值模拟[J].力学与实践，2000，22(5)：45-48.

[14]GB50218—1994 工程岩体分级标准[S].

[15]田昌贵，陈世华.工程岩体分级标准在地下采矿工程中的应用[J].采矿技术，2005，5(4)：89-93.

[16]蔡斌.国标《工程岩体分级标准》应用中的几个问题[J].岩土力学，2003，24(增刊1)：74-76.

[17]朱哲明，李元信，周志荣，等.爆炸荷载下缺陷岩体的动力响应[J].岩石力学与工程学报，2011，30(6)：1157-1167.