爆破地震的数值模拟及爆破振动规律分析

张智超，刘汉龙，陈育民，王维国

（1.河海大学 岩土力学与堤坝工程教育部重点实验室，江苏 南京210098；2.河海大学 土木与交通学院，江苏 南京210098）

0 引言

土体地震动力响应一直是岩土领域的研究重点.以往所开展的地震动力试验，主要是通过动三轴等试验的结果来分析土体的地震液化特性和动力响应［1］，以及振动台、剪切模型箱和离心机等配合使用的手段来还原土工建筑或土体结构遭受天然地震时的受力状态［2－5］.然而，常规1 g室内试验难以满足土体原型应力水平，而离心机试验虽然能够还原土体的实际应力状态，但也存在一定缺陷，如离心应力场不均匀，在沿着离心力作用的半径方向上，离心力与半径大小成正比，这就给原位应力场的模拟带来一定误差，且离心机造价昂贵，技术复杂，很大程度上制约了其推广使用.

因此，王兰民［6］采用现场微差爆破试验的方法，创造出与天然地震相似的持续振动状态，以此来分析场地的地震震陷特征，取得了不错的效果；此外，通过爆破手段来模拟天然地震的试验研究也不断地为学者们所尝试［7－10］，以试图避免前述室内实验和模型试验的缺陷，从原位的角度来研究土体力学行为.

目前，对此类利用微差爆破模拟天然地震的研究，主要集中在现场试验方面.然而，现场试验的一些弊端，如对场地要求较高，数据监测范围有限，较为耗费人力、物力，难以进行大量重复等，都在一定程度上制约了爆破模拟天然地震的研究.

作者利用数值手段，对文献［6］中的微差爆破试验进行了模拟，将计算结果与实测值进行比较，结合理论分析，为数值结果的正确性提供保障；在此基础上，对地层中的爆破振动规律进行了分析，丰富了爆破地震动的研究，对将来利用现场爆破试验来模拟天然地震的研究工作具有一定指导意义.

1 现场试验

试验［6］位于甘肃某黄土场地，场地深度0～1 m为植耕土，1～15 m为风成黄土，15～28 m为次生黄土，28 m以下为卵砾石层.

炮点布局如图1所示，作为爆源的30个药包埋设在以场地中心为圆心、半径为15 m的圆周上，等间距分布，埋深23 m.每个药包装有40 kg二号岩石炸药，轴对称上的两个药包作为一对同时起爆.如图1所示，炮点边上的标序表示引爆顺序，而炮点边上无标序则表示此处药包引爆失败，模拟中将忽略之.

图1 药包布局及引爆顺序［6］Fig.1 Distribution of explosive charges and their detonation or der［6］

场地中心的地表布置了三分相标准强震仪，用于捕捉试验过程中的爆破振动加速度时程.

2 数值模型

2.1 几何模型

根据试验现场情况，数值模型采取直径50 m、高度30 m的圆柱体计算区域，如图2所示.分别将深度0～1 m、1～15 m、15～28 m和28～30 m的土体划分为不同材料，用不同颜色表示.x向对应现场的WE（东西）向，y向对应于现场的NS（南北）向，z向对应于现场的UD（竖直）向.30个药包在深度23 m、半径15 m的圆环上等间距分布，每个药包质量为40 kg.

图2 有限元模型Fig.2 FEMmodel

2.1.1 有限元网格划分

模型中的所有单元均为六面体八节点实体单元，并在药包附近进行了网格细化，共计527 085个单元和545 184个节点；同时，毗邻单元间的网格尺寸比都小于1∶4，并对单元大小进行限制，通过适当的网格渐变，将单元尺寸由炸药单元的0.2 m逐渐过度到最外围土体单元的1 m，即远离炸药的单元网格逐渐变得稀疏.这在一定程度上保证了计算精度，防止单元退化.

2.1.2 边界条件

模型上表面设置为自由边界，以示与空气接触，四周与底部设置为无反射边界，以真实地反映波在此的透射情况.此外，药包与底部边界距离为7 m，与圆周边界距离10 m，在一定程度上也保证了足够的计算范围，消除了应力波的反射.

2.1.3 求解时间

孙军杰［10］对实测加速度结果进行了处理，合理地截掉了6对未爆爆源所导致的振动中断时段，得到了叠续合成的振动时程.因此，数值计算中也将对这6段未爆爆源的微差时段进行忽略，只考虑成功引爆的微差时段，总计算时间为6 s.根据试验中的微差时间，计算爆破时刻如表1所示.

表1 数值计算的爆破顺序Tab.1 Blasting sequence in the numerical simulation

式中：P为压力；φ为内摩擦角；J2为偏应力张量第二不变量；K（θ）为张量平面角的函数；c为黏聚力；A HYP为决定修正后的Mohr－Coulo mb屈服面和标准Mohr－Coulomb屈服面相似度的参数.

当AHYP＝0时，上式表示标准 Mohr－Coulomb屈服面，当AHYP采用较大值时，修正后的屈服面明显偏离标准屈服面，对于数值模拟来讲，A HYP的取值应该小于c cotφ，一般按照下式进行选取：

2.2 物理模型

2.2.1 土体本构

LS－DYNA中的 MAT＿FH WA＿SOIL［11］本构可以模拟诸多土体特性，如爆炸荷载下的应变率效应、应变硬（软）化以及孔压等，是较为合适的计算土体爆炸响应的模型.该模型在爆炸分析方面的优越性已在文献［7］中得到充分验证.因此，作者选取其来描述现场爆破试验中土体的力学行为.

该模型将标准Mohr－Coulo mb屈服面修正为光滑的曲面，并且垂直于压力轴，提高了数值计算的有效性和稳定性.屈服函数为：

图3为标准Mohr－Coulomb屈服面与修正的屈服面对比.除了在低应力区外，二者几乎一致.

深度0～28 m的黄土层采用MAT＿FH WA＿SOIL本构.根据现场探测的剪切波速和压缩波速［6］，按 照 公 式算出各土层的剪切模量和体积模量；由于黄土含水率较低，将孔压系数Ksk设置为0，即不考虑孔压发展，进行总应力分析；其余参数根据现场实测数据并结合相关经验得出［11，13］，参数值见表2.

图3 标准与修正后的Mohr－Coulomb屈服面对比Fig.3 The standard and modified Mohr－Coulomb yield surface

表2 MAT＿FHWA＿SOIL模型参数Tab.2 Parameters of MAT＿FHWA＿SOIL

对于第28～30 m的卵砾石层，则采用线弹性本构，密度2 000 kg／m3，杨氏模量5e8 Pa，泊松比0.25.

2.2.2 炸药材料模型

采用高能炸药燃烧模型和J WL状态方程来模拟现场试验中所使用的二号岩石炸药［13］.JWL状态方程能够精确描述爆炸过程中爆轰产物的压力、体积、能量特性，其关系式为

式中：A，B，R1，R2，ω 为材料参数；E01为爆轰产物单位体积的内能；V为单位体积装药产生的爆轰产物的体积；P为爆炸爆轰压力.炸药材料参数取值分别为［14］：密度1 000 kg／m3，爆轰速度3 600 m／s，CJ压力3.24 GPa，A＝214 GPa，B＝0.182 GPa，R1＝4.15，R2＝0.95，ω＝0.38，E01＝4.5 GPa，V＝1.

2.3 结果分析

在爆破计算前，利用LS－DYNA本身的动力松弛算法对土体初始自重静应力进行计算，结果如图4所示.可以看出，自重结果合乎规律，为其后的爆破分析提供了准确的初始应力条件.

图4 初始自重应力云图侧视图Fig.4 Side view of initial gravity stress nephogram

2.3.2 模拟与实测地表加速度时程的比较

为了验证模拟结果的准确性，分别将模拟得到的x、y和z向加速度时程与所对应E W、SN和UD向加速度时程实测值［10］进行比较，图5是对比图.由图可以看出，模拟的加速度时程出现了9个较为明显的峰值群，与现场9对药包依次爆破后激发的地震动相对应.由于地震本身就是一种杂乱无章的随机振动现象，加之现场监测数据的离散性，以及数值计算的近似性，无法做到计算曲线与实测曲线的完全一致，但由图仍可看出，模拟的x、y向加速度时程与现场实测的EW、SN向加速度时程在幅值、曲线形状和峰值出现的时刻都吻合得较好，从一定程度上说明了数值手段的准确性；同时还可以看出，微差爆破在场地中心所产生的连续的爆破地震波时程曲线在宏观形态上与典型的天然地震波曲线具有比较好的相似性，因此，可以通过爆破手段来对原位土体场地的地震响应进行分析，其反映出的总体振动趋势大体上可以被用来模拟天然地震.

值得一提的是，对于模拟的z向加速度时程与实测竖向加速度时程，虽然在曲线形状和峰值出现时刻上较为吻合，但峰值大小差别相对较大，模拟的峰值普遍接近2 g，且在第一组峰值群有着不合常理的跃迁；而实测峰值普遍不到1 g，模拟值几乎是实测峰值的2倍.分析原因可能是由于计算参数不够精确以及有限元计算方法的近似性等造成的，加之实际环境的干扰，如土质不均匀、监测设备精度有限等，与数值计算中完全理想的情况有所不同，就更容易导致实测值与模拟值的差异，但我们仍可以通过数值手段预测爆破震动的大致规律.

2.3.3 爆破振动随地层深度的变化规律

基于上述结论，选取了从深度23 m——即炸药埋深至地表共10个不同深度的节点A～J作为参考，分析爆破振动变化规律.参考节点位置见图6.

图6 参考节点示意图Fig.6 Location of reference nodes

图7 为不同深度x、y和z向加速度时程，子图中的上图为时程全纪录，下图将加速度特征局部放大，以便更为清晰地展示其变化规律.由图7可以看出，随着埋深的减小、爆心距的增加，各节点x、y向（水平向）加速度幅值逐渐减小，而z向（竖向）加速度幅值则先增大、后减小，到了近地表又趋于稳定.之所以产生这种现象，原因分析如下：爆炸引起的球形应力波在土体空间中扩张，爆炸振动最为剧烈的方向近似在爆心与参考质点的连线方向上，即球形波阵面的外法线方向，因此，对于位于一对同时引爆的药包连线中点上的质点来说，z向为球形波阵面的切线方向，振动强度趋于最小，所以，当参考节点埋深减小，即节点上移、离开了药包的连心线中点后，z向的振动强度分量增大，这导致了z向加速度幅值的首次增大；进而，随着节点深度的继续减小，爆心距也逐渐增大，由于爆破振动强度与爆心距呈反关系，这就导致了z向（竖向）加速度幅值产生了减小的趋势；但药包附近存在自由面时，土体中的应力状态将变得较为复杂［15］，自由面的存在，会增强爆破效果，较之自由面切线方向（x、y向）来说，自由面反射的增强作用对自由面法线方向（z向）的影响更大，因此，这就导致临近自由面时，反射波的放大作用和爆心距增加而导致的减弱作用大致相互抵消，所以虽然爆心距增加了，但质点竖向振动强度并不随着爆心距的增加而减小，而是基本维持稳定；另一方面，自由面切向的质点运动特征受到自由面反射的影响相对较小，因而切向的振动幅值并未因为临近自由面而增强，而仅是随着爆心距的增加呈减小的趋势.

质点速度与加速度的变化规律呈现一致性，即从炸药所在的埋置深度向上统计，随着深度的减小、爆心距的增大，质点x、y向速度幅值呈现逐渐减小的趋势，而质点z向速度幅值则先增大、后减小，最后在临近自由面处又趋于稳定.分析原因与前述类似，在此不再赘述.

3 结论

对微差爆破模拟天然地震的现场试验进行了数值分析，可以得到以下几点结论.

（1）数值手段可以很好地还原微差爆破作用下的土体质点振动特征.数值手段简便易行，可灵活地针对不同的研究目标和内容进行调整，因此，开展爆破地震动的数值分析工作，对将来利用现场爆破来模拟天然地震的试验设计、施工及爆破振动预测都具有一定指导意义.

（2）从数值计算的角度验证了采用爆破手段来模拟天然地震在一定程度上是可行的，可以利用爆破地震波创造出与天然地震相似的震动环境，进而较为符合实际地进行一系列原位土体地震响应研究，避免其他室内试验和模型试验的缺陷.

（3）在本文的算例中，自炸药埋深向上，随着深度的减小、爆心距的增加，爆破水平向加速度幅值和速度幅值逐渐减小，而竖直向加速度幅值和速度幅值则先增大、后减小，到了近地表又因自由面反射效应的影响而基本维持稳定.

［1］ 何昌荣.动模量和阻尼的动三轴试验研究［J］.岩土工程学报，1997，19（2）：39－48.

［2］ 吴俊贤，倪至宽，高汉棪.土石坝的动态反应：离心机模型试验与数值模拟［J］.岩石力学与工程学报，2007，26（1）：1－14.

［3］ 汪明武，Susu mu IAI，Test uo TOBITA.液化场地堤坝地震响应动态土工离心试验及模拟［J］.水利学报，2008，39（12）：1346－1352.

［4］ 黄春霞，张鸿儒.大型叠层剪切变形模型箱的研制［J］.岩石力学与工程学报，2006，25（10）：2128－2134.

［5］ 李敏霞，杨泽群，陈建秋.地震模拟振动台技术的开发与应用［J］.世界地震工程，1996，2（1）：49－54.

［6］ 王兰民，孙军杰，徐舜华，等.爆破模拟地震动条件下黄土场地震陷研究［J］.岩石力学与工程学报，2008，27（5）：913－921.

［7］ Lee，WAYNE Y.NUMERICAL MODELING OF BLAST－INDUCED LIQUEFACTION ［D］.Brigham Young University，Ph D Thesis.August 2006.

［8］ PAT HIRAGE K S.Critical assessment of the CANLEX blast experi ment to facilitate a development of an in－situ liquefaction methodology using explosives［D］.Depart ment of Civil Engineering，Degree of Master of Science，2000.

［9］ ASHFORD S A，ROLLINS K M.Blast－Induced Liquefaction for Full－Scale Foundation Testing［J］.Jour nal of Geotechnical and Geoenviron mental Engineering.2004，130（8）：798－806.

［10］孙军杰.黄土场地震陷与桩基负摩阻力现场试验研究 ［D］.兰州：兰州大学土木工程与力学学院博士学位论文，2010.

［11］Lewis B A.Manual for LS－DYNA Soil Material Model 147［R］.Federal Highway Administration，2004.

［12］徐芝纶.弹性力学 ［M］.北京：人民教育出版社，1983.

［13］Anon.Ls－dyna key word users’menual（Version 971／Release 4）［M］.s.l.Liver more Soft ware Technology Cor poration，2009.

［14］刘红岩，王新生.地表垂直爆破震动速度的数值计算［J］.爆炸与冲击，2008，28（3）：255－270.

［15］亨利奇.爆炸动力学及其应用 ［M］.熊建国，译.北京：科学出版社，1987.