基于FLAC3D的岩质边坡爆破动力响应规律研究

唐旭　方正峰　邹飞

摘要：边坡的动力稳定一直以来都是岩土工程及地震工程重点研究的问题。利用FLAC3D数值计算方法，对爆破作用下岩质边坡振动速度、位移、应力等进行模拟，并在与实测爆破振动峰值进行对比验证的基础上，分析边坡的动力响应规律。结果表明：地震波具有极强的方向性，主要朝施加荷载的正向传播;在边坡坡面和坡脚处因应力波在边坡自由面反射容易产生一定区域的振动加强场;应力的响应特征上，无论是拉应力还是压应力其最大值均出现在坡脚附近，在实际的边坡爆破开挖过程中，对边坡的坡面尤其是坡脚处应该加强保护措施。研究成果对岩质边坡开挖设计及施工具有一定的参考价值。

关键词：动力响应;岩质边坡;FLAC3D;爆破

中图法分类号：P642

文献标志码：A

DOI：l0.16232/j.cnki.1001-4179.2019.03.035

1研究背景

爆破作为一种最有效、最经济和最便捷的施工方法已经普遍应用到山地公路路堑岩质高边坡的形成中。巨大的工程规模伴随而来的是大量的爆破施工作业，而边坡的动力稳定一直以来都是岩土工程及地震工程重点研究的问题之一。

随着科技的发展，越来越多的研究人员借助计算机进行仿真模拟，提出了各类数值分析方法，并开发了相应的数值分析软件，如基于动力有限元法的Ansys/Ls-dyna模块[1]，离散元法的UDEC[2]，快速拉格朗日法的FLAC3D等[3-4]。

蒋楠用Ansys/Ls-dyna建立了不同坡度的边坡爆破振动数值计算模型，获得了不同坡度边坡在同一水平位置的各坡面监测点的质点爆破振动速度变化趋势，且同一爆源距的质点峰值振速随高程增大表现出高程放大效应;在边坡坡度相近时，质点振速随爆源距及高程的同时增加，高程效应不明显;不同坡度边坡，在同一水平质点振速随着边坡坡度增加以衰减为主，即质点振速在爆源距增加及高程增加的同时作用下，而不再表现为常规的振动速度随爆源距增加衰减的现象，而是以质点振速随高程增加而呈现增长的现象[5]。唐海，李海波等应用UDEC程序模拟了爆破振动波在不同地貌中的传播特征，得出凸形地貌下放大系数并不随台阶型凸形地貌的增高而单调增加，在孤立凸形地貌下，放大系数还与凸形地貌的高宽比有关[6-7]。言志信、祁生文等利用FLAC3D软件对不同坡高、坡角的边坡动力响应特征做了数值计算分析研究，获知当边坡坡高较低时，在一定的坡高范围内边坡内的质点振动速度随着高程的增大表现出高程放大效应，且速度放大系数等值线图呈近等间距分布;当边坡的坡高高出一定的范围时，边坡质点振动速度不再随高程的增大表现出线性的反放大，而是表现出节律性的变化，在坡肩附近时，才再一次表现出放大的迹象[8-9]何铮等结合西南某水电战顺层岩质高边坡工程实例，研究了该边坡地震后工程地质特性，恢复边坡震前的地形地貌，建立地震前未滑塌边坡的三维模型，而后利用FLAC3对模型进行动力-时程响应分析，绘制出位移、速度、加速度三量放大系数等值线图[10]。

岩质高边坡爆破开挖过程中引起的地震效应已成为影响边坡稳定性的重要诱发因素之一，越来越多的学者也开始关注岩质高边坡在爆破荷载作用下的动力响应特征问题。爆破荷载对边坡的稳定性影响主要是动荷载作用下岩体的响应特征问题[11-12]，因此掌握并深入分析岩质高边坡对爆破振动效应的响应特征规律，探索爆破开挖施工对边坡稳定性的影响，对于进一步控制和预测爆破振动效应对岩质边坡的危害，同时减少爆破作业对周边环境的影响具有重要的指导和实践意义[13-15]。

本文利用FLAC3D建立公路岩质边坡的数值分析模型，分析岩质边坡爆破振动下振速、位移、应力相应的响应特征，得到并分析岩质边坡爆破振动下振速、位移、应力相应的响应特征，细化研究传播路径中地质条件、地貌特征的影响。

2FLAC3D动力计算原理

FLAC3D的基本原理是采用线性显式拉格朗日差分算法和混合离散技术求解运动方程，将指定计算区域划分成若千个单元，单元之间由节点连接，在对某一节点施加荷载后，此节点的运动方程可以表示为时间不长的有限差分形式，在一微小时间步长内作用于节点的荷载只对相邻的一些节点产生影响，即根据某一t时刻的应力状态和△t时间步长应变增量，进而确定t+△t时刻的应力状态。根据高斯（Gauss）定律，由节点的速率求解单元的应变增量;再利用应力-应变关系，即本构方程计算单元应力，并对其进行积分，得到作用在节点上的应力矢量;最后结合平衡方程，由节点力进一步求解出节点速率及位移。如此反复循环迭代，这一过程将随时间步长正常扩展到整个计算范围直至计算收敛，从而准确地模拟岩体及其他材料的塑性破坏和流动。其计算流程如图1所示。

在速度和位移的动力计算中，FLAC3D以节点为计算对象，分别将质量和力都集中在节点上，再通过运动方程在时域内进行求解，节点运动方程可以表示为

公式

式中，Fli（t）为t时刻l节点在i方向的平衡力分量，可由虚功原理导出;ml为l节点的集中质量。

将式（1）左边用中心差分来近似，则可以得到节点的速度：

公式

此时，可以进一步得到节点位移：

公式

3數值计算模型的建立

3.1计算模型

图2为某高速公路K86+152～K86+600段地形地貌示意图，本文依托该工程对岩质边坡爆破振动响应特征进行数值研究。

本文采用三维计算模型，X正向为边坡右侧，Y正向为线路大里程方向，Z轴铅垂向上，模型计算区域为252mx230mx137m。。输入波动的最短波长控制着网格划分的尺寸，研究表明，当网格的最大尺寸为△l时，输入波动的最短波长为λ时，必须满足△lt<（1/10～1/8）λ。所建立模型及其网格划分如图3所示。

3.2本构模型及材料参数

本文模型选取莫尔-库伦（Mohr-Coulomb）弹塑性模型，该模型用于描述岩土体的力学行为。

该里程段岩性较为均一，仅为灰岩。模型中所需的岩石物理力学参数根据室内实验确定，详见表1。

3.3吸收边界

在FLA3D中可以设置吸收边界（也称静止边界、黏性边界）和自由场边界来减少边界中波的反射。

本文将X正、负向，Y正、负向4个侧面及Z负向地面均设置为黏性吸收边界，地表均为自由边界，详见图4。

3.4等效荷载输入

炸药在爆炸瞬间产生高温高压的爆轰气体，作用于炮孔壁，使其产生破碎区和塑性变形区，并且瞬间衰减为应力波沿岩体继续传播，整个过程仅仅持续数毫秒。其间影响因素甚多且极其复杂，目前量测技术有限，无法量化并准确地获知爆炸过程中的每一个具体细节，只能通过相关经验将爆破过程中荷载作用及变化情况等效处理，通常是将数值模拟中假定荷载以应力时程的形式输入。本文拟采用三角脉冲荷载，这种荷载的建立需要确定两个基本因素[16-19]。

3.4.1爆破振动峰值荷载

峰值荷载的确定即确定炸药爆炸过程中产生的爆轰气体作用在炮孔壁上的峰值压应力。根据爆轰波的C-J理论，炸药作用在炮孔壁上的平均炮轰压力为

公式

式中，Pj为瞬时炸药爆轰平均初始压力;ρe为炸药密度;D;为炸药爆轰速度;γ为炸药的等熵系数，y的取值与装药密度有关。研究表明[20]，当ρe<1.2g/cm3时，γ取2.1，当ρe≥1.2g/cm3时，y取3。

在耦合装药条件下，炸药作用在炮孔壁上的初始峰值压力P0为

P0=Pj

在不耦合装条件下，炸药作用在炮孔壁上的初始峰值压力为

公式

式中，k为装药径向不耦合系数，

公式

分别为炮孔直径和药卷直径。

研究表明，爆轰气体作用于炮孔壁时，其压力将明显增大，故得到不耦合系数较小的柱状药包爆炸对炮孔壁的冲击初始峰值压力为

公式

式中，n为爆炸产物膨胀作用炮孔壁时的压力增大系数，n=8～10，本文取10;l为轴向不耦合系数，

公式

;le、lb分别为装药长度和炮孔长度。

表2为爆破荷载峰值压力计算相关计算参数。

3.4.2爆破荷载作用时间

炸药爆炸作用时间非常短，通常认为炸药爆破时冲击波作用的持续时间约1.0x10-6～0.1s，爆轰气体压力作用时间约为1.0x10-3～0.1s。本次数值计算中爆破荷载作用时间假设为7ms，其中爆破荷载压力上升时间为1ms，'下降时间为6ms，且荷载按现场实际爆破设计情况单排爆破作用于K86+320～K86+350段第一级边坡坡面上，图5为爆破荷载施加示意图。

3.5力学阻尼的选取

岩体中的内摩擦以及可能存在的接触表面滑动是产生阻尼的主要原因。在动力计算中，动力反应的阻尼需要在数值模型计算分析时模拟重现自然系统中动荷载作用下的阻尼大小。

FLAC3D的动力计算中提供了3种阻尼形式，分别为滞后阻尼、局部阻尼及瑞利阻尼。实践证明，瑞利阻尼因其理论与常规动力分析方法类似，虽然采用瑞利阻尼时的计算时间步长较小，但是计算得到的加速度响应规律相较于前二者更加符合实际情况，故也是动力计算中通常使用的阻尼形式。

综合考虑模型计算时间及其模拟的准确性，经过反复调试，本文选取最小临界阻尼比为0.02，最小中心频率为10Hz。

4数值计算结果验证

为了验证数值计算结果，依托某在建高速公路爆破施工工程现场监测数据，实际爆破参数设计如下。

（1）钻孔直径d=90mm。

（2）台阶高度H=10m。

（3）最小抵抗线W光=K光d=20x120mm=2400mm，式中W光为光面爆破最小抵抗线，K光为光面爆破计算系数，一般取K=15～25，岩石越硬K取值越小。该工程取K=20。

（4）炮孔超深h=0.5～1.5m，孔深大和岩石坚硬完整者取大值，反之取小值。此工程取h=1.0m。

（5）孔距光=mW=（0.5～0.8）x2.4=（1.2～1.92）m，式中m為装药密集系数。该工程取光=1.2。

（6）炮孔长度L=（H+h）/sin。当坡度为1：0.75时，L=13.8m;当坡度为1：0.5时，L=12.3m。

在模型上布设相应监测点来跟踪监测其速度时程，位于K86+335断面的5个监测点位置，详见图6。因选取的数值计算荷载施加相关计算参数与现场监测时的实际爆破设计装药参数相近，且该工况实际测点布置与数值模型中测线位于同一断面，将二者峰值振动速度进行对比，得到数值计算各监测点时程曲线如图7所示，表3为数值计算各监测点与现场监测点峰值振动速度对比结果，并将其绘成折线图如图8所示。

由表3及图8可知，数值模拟结果与现场监测各测点峰值振动速度在爆源近1～3号测点区其相对误差均在10%以下，到了爆源远区误差相对较大，且数值模拟峰值振速普遍大于现场监测数据，可能原因是简化模型为均匀连续介质，忽略了实际岩体内部节理裂隙等不连续结构面，故而使爆破振动地震波的衰减速度相较与实际工程略有减小。总体而言，数值计算结果的精度尤其是在爆源近区的都在可接受范围内，说明借助FLAC3D研究该边坡爆破振动动力响应特征是可行的。

5动态响应分析

5.1速度响应分析

图9为爆炸后5，10ms时模型X向速度分布云图。由图可知爆破地震波的作用影响范围，因爆破振动地震波具有方向性，特别是对于这种深孔台阶开挖，随着时间的推进，虽说爆破荷载施加的反方向岩体也有一定的响应，但爆破振动地震波主要朝着施加荷载的右边坡传播。总的来说在这种边坡爆破中，地震波的传播可以概化成如图10所示的椭圆等振线。故可以得到的结论是，右边坡的爆破开挖对既有左边坡影响可以忽略不计。

图11为爆炸后5，10ms时K86+335断面X向速度分布云图。可清晰地观察到，随着爆破地震波的传播，在坡脚右侧某一半径内质点振动速度在某一区域得到加强，其原因是应力波在边坡自由面反射而产生了一定区域的振动加强场;又根据图11（b）所显示的各台阶周边速度大小分布情况可知，该时刻速度值在第2级台阶坡脚处形成一个极大中心，顺坡而上其速度值是先减小又增大，并在第4级台阶坡脚处再次形成一个极大中心，随后减小。

5.2位移响应分析

为进一步研究边坡在爆破振动荷载作用下相关参数的响应特征分布情况，引入无量纲位移放大系数，其定义为边坡爆破振动位移响应大小与第1级台阶边缘位移大小的比值。

根据数值计算结果，绘制出边坡地表及K86+335断面在10ms时Z向位移放大系数等值线图如图12所示。由图12（a）可知，位移响应特征与前述地震波椭圆等振线分布情况基本一致。图12（b）中的位移放大系数等值线分布情况再次反映了右边坡爆破开挖未对左边坡造成的影响。为进一步分析边坡坡体内的位移响应特征规律，将图12（b）右边坡局部放大显示如图13所示。

由图13可知，爆破荷载施加的1级边坡其位移在一定深度范围内随深度的增加而出现极大中心，超过某一深度后位移放大系数减小;再观察位移放大系数等值线随高程增加所表现出来的特征时发现，其分布规律为位移放大系数在1级边坡表面先减小，顺坡面而上在进入2级台阶时位移放大系数逐渐增大，随后在2级边坡中部位置顺坡面而上位移放大系数开始减小。此时位移放大系数所表现出的分布规律与速度响应大小的分布规律近乎一致。

5.3应力响应分析

图14为K86+335断面炸药爆炸后不同时刻最大主应力分布云图。

由图可知，各个时刻除爆源附近一定区域内应力在进一步不断地调整，边坡其他区域等值线变化不大，其原因是所施加的荷载相较于自重荷载比较小，尚不足以使整个边坡的应力重分布。因爆破地震波主要以压缩波在坡体内传播，爆破荷载作用下边坡受影响区域多出现正常的压应力，但在局部如边坡坡脚附近、个别坡面及台阶边缘，会出现一定区域的拉应力，而出现拉应力的区域基本上都是接近坡面，即自由面，其可能原因是爆破振动应力波在自由面发生发射，进而生成一定的拉伸波造成相应区域出现拉应力。

由图14可知，爆后随着应力的逐步调整，无论是拉应力还是压应力，最大值均出现在坡脚，即靠近爆源附近，且相应的应力值均小于岩石相应的临界强度值，在实际的边坡开挖爆破中，也是坡脚处最容易率先发生局部的失稳破坏。又由6～10ms的应力分布云图可知，压应力主要集中区域在该时段内逐渐向深处移动，而因应力波逐渐传播至坡面，同时在坡面发生反射，生成一定的拉伸波，该时段坡面逐渐开始产生拉应力，且逐渐地在坡面附近产生拉应力集中区域;10ms过后，坡面周边拉应力集中区域逐渐消散，仅在边坡坡脚及1，2级台阶突出物保持着较小拉应力集中区。同时，在10～1000ms，受爆破动荷载影响，较大的坡脚及1，2级台阶附近的应力状态几乎不发生变化，应力调整基本趋于稳定，其余部分逐渐向初始状态下的应力分布进行调整。

6结论

（1）边坡开挖爆破中，地震波的传播可以概化为椭圆等振线，地震波主要朝施加荷载的正向传播，对既有左边坡的影响较小。

（2）在岩质均一的边坡中，无论是速度还是位移，其爆破振动的响应规律较为一致，且在坡面和坡脚处因应力波在边坡自由面反射容易产生一定区域的振动加强场。

（3）在岩质边坡爆破作用影响下，无论是拉应力还是压应力其最大值均出现在坡脚附近。

（4）在实际的边坡爆破开挖过程中，对边坡的坡面，尤其是坡脚处应该加强保护措施。

参考文献：

[1]Aysilio E，Conte E，Dente G.Seismic stability analysis of reinforcedslopes[J].Soil dynamics and earthquake，2000，19（3）：159-172.

[2]Spyros Sklavounos，FotisRigas.Computer simulation of shock wavestransmission in obstructed Terrains[J].Journal of Loss Prevention inthe Process Industries，2004，（17）：407-417.

[3]Siyahi，Bilge G.Pseudo-static stability analysis innormally consolidated soil slopes subjected to earthquake[J].Teknik Dergi/Technicaljournal of Turksih chamber of civil enginers，1998，9（DEC）：457-461.

[4]UDEC（Universal Distinct Element Code）user' s manual Version 3.0[M].Lundon：Itasca Consulting Group，Inc.，1996.

[5]蔣楠，周传波，平雯，等.岩质边坡爆破振动高速度高程效应[J].中南大学学报：自然科学版，2014，45（1）：237-243.

[6]唐海，李海波，蒋鹏灿，等.地形地貌对爆破振动波传播的影响实验研究[J].岩石学与工程学报，2007（9）：1817-1823.

[7]唐海，李海波.地形地貌对爆破振动波传播的影响实验研究[J].岩石力学与工程学报，2007（09）：1817-1823

[8]言志信，蔡汉成，王群敏，等.纵波作用下边坡动力响应[J].工程地质学报，2010，18（4）：487-505.

[9]祁生文，伍法权.边坡动力响应规律研究[J].中国科学E辑技术科学，2003，33（s）：28-40.

[10]何铮，徐卫亚，石崇，等.顺层岩质高边坡地震变形破坏机制三维数值反演研究[J].岩土力学，2009，30（11）：3512-3518.

[11]陈明，卢文波，李鹏，等.岩质边坡爆破振动速度的高程放大效应研究[J].岩石力学与工程学报，2011（11）：2189-2195.

[12]李鹏，苏生瑞，王闫超，等.含软弱层岩质边坡的动力响应研究[J].岩土力学，2013（S1）：365-370，378.

[13]孟祥甜.露天边坡动载响应特性的数值模拟研究[D].包头：内蒙古科技大学，2015.

[14]宋光明，陈寿如，史秀志，等.露天矿边坡爆破振动监测与评价方法的研究[J].有色金属，2000（4）：24-27.

[15]王在泉，陆文兴.高边坡爆破开挖振动传统规律及质量控制[J].爆破，1994，11（3）：1-4.

[16]夏详，李俊如，李海波，等.爆破荷载作用下岩体振动特征的数值模拟[J].岩土力学，2005，26（1）：50-56.

[17]Kyhlemeyer R L，Lysmer J.Finite Element Method Accuracy forWave Propagation Problems[J].Soil Mech & Foundations Div，ASCE，1973，99（5）：417-421.

[18]Josef Henrych.The Dynamics of Explosion and Its Use[M].NewYork：Elsevier Scientific Publishing Company，1979.

[19]S.G.Chen，J.Zhao.UDEC modeling of a field explosion test[J].International Journal of Blasting and Fragmentation，2000（4）：149-163

[20]Josef Henrych.The Dynamics of Explosion and Its Use[M].NewYork：Elsevier Scientific Publishing Company，1979.

引用本文：唐旭，方正峰，鄒飞.基于FLAC3的岩质边坡爆破动力响应规律研究[J].人民长江，2019，50（3）：198-204.

Study on dynamic response laws of rock slope blasting based on FLAC3D

TANG Xu'，FANG Zhengfeng2，ZOU Fei，4

（1.Guizhou Province Quality and Safety Traffic Engineering Monitoring and Inspection Center Co.，Ltd.，Guiyang 550081，China;2.Guizhou Hongxin Chuangda Engineering Detection & Consultation Co，Ltd.，Guiyang 550014，China;3.China GuizhouTransportation Planning Survey & Design Academe，Guiyang 550081，China;4.School of Ciwil Engineering，Central South University，Changsha 410000，China）

Abstract：The dynamic stability of slope has always been the focus of geotechnical engineering and seismic engineering.Inthis paper，FLACED numerical calculation method is used to simulate the vibration velocity，displacement and stress of rocky slopeunder blasting action.And on the basis of comparing and verifying the peak value of blasting vibration，the dynamic responselaws of rocky slope was analyzed.The results show that the seismic wave has a strong directionality and propagates mainly to thepositive direction of the applied load.On the slope surface and at the foot of the slope，it is easy to produce the vibrationstrengthening field in a certain area due to the reflection of the force wave on the free surface of the slope.The maximum value ofstress，whether tensile stress or compressive stress，appears near the slope foot.In the actual blasting excavation，the slope surface should be protected，especially the slope foot.The research results have certain reference for the construction design of rockyslope excavation.

Key words：dynamic response;rock slope;FLAC3D;blasting