随着城市地铁的发展，隧道爆破掘进穿越坚硬岩土层的案例屡见不鲜^[1-2]，但由于地铁沿线建筑林立、地下管线密布，隧道施工过程中爆破诱发的环境振动效应成为了该工程领域研究的重点和难点。爆破振动强度过大导致围岩破坏^[3-4]、爆破引起建筑物共振以及高频爆破振动迫使刚性结构产生局部应力集中^[5]，均可能导致爆破影响区内的构筑物产生不同程度的损伤。为减振降害，通过电子雷管或非电雷管精确控制多孔毫秒延时爆破，既能显著改善岩石破碎效果，又能有效减小爆破振动^[6]。不少学者从大量的爆破资料中总结出毫秒延时爆破的工程参数^[7-8]，如毫秒延时爆破最佳间隔时间^[9]、最优单孔药量以及爆破顺序^[10]等。然而，由于实际工程地质条件千差万别^[11]，此类现场爆破试验往往具有较强的地域性，且现场试验前单孔药量及爆破顺序等设计参数尚不明确，容易导致相邻建筑物损伤。因此，合理有效的数值模拟成了较理想的预测爆破效果和确定爆破参数的方法^[12]。然而，由于延时爆破中各起爆点的单孔装药量不同，爆破荷载不规则，而常规的数值模拟较难反映此类不规则爆破荷载作用下岩土材料的动力响应特征。

针对上述问题，本文中基于南京地铁隧道毫秒延时爆破环境振动特性现场试验，首先建立有限元-无限元耦合边界的地层-爆源体系的三维精细化有限元模型。其中，考虑爆破荷载的不规则特性，采用基于非对称加卸载准则修正的Davidenkov本构模型描述场地土体的动力非线性特性；通过改进的Friedlander方程来模拟内源爆炸在圆柱形炮孔表面产生的瞬态空气冲击波。通过与已开展现场爆破试验的实测数据进行对比，验证该模型方法的有效性。最后，对齐爆和毫秒延时爆破的峰值振速及频谱特性进行对比，探讨毫秒延时爆破环境振动特性及减振机理。

1.   现场爆破试验回顾

1.1   工程背景

南京地铁四号线云南路～鼓楼站区间全长约660 m，拟采用矿山法施工，包括区间隧道、风道、鼓楼小里程暗挖段和三联拱部分。鼓楼站东侧为既有地铁一号线鼓楼站，东侧有已通车的鼓楼隧道，西侧为鼓楼公园，公园环岛中间为文物古建筑。区间隧道距鼓楼城阙仅约1.3 m，竖向距离14.5 m；1#风道距城阙约为19.8 m，竖向距离13.5 m；三联拱距城阙约为18.1 m，竖向距离11.0 m，暗挖结构与文物鼓楼的相对位置如图1所示。

图  1  地铁隧道与文物鼓楼的相对位置

Figure  1.  Relative position between the metro tunnel and the cultural relic drum tower

下载: 全尺寸图片 幻灯片

1.2   现场爆破试验方案

由调查可知，文物鼓楼的病害较严重，且区间隧道距鼓楼城阙水平距离最近仅约1.3 m，为避免直接近区爆破对文物鼓楼的安全性产生影响，在文物鼓楼周边选取了爆破试验区域。爆破试验前，采用环境随机激振法（也称脉动法）对鼓楼的固有频率进行测试，结果列于表1。

表  1  鼓楼固有频率测试结果

Table  1.  Test results of natural frequencies of the drum tower

模态阶数	固有频率/Hz
	水平向	竖直向
1	1.37	2.00
2	2.79	2.98

下载: 导出CSV 

| 显示表格

结合现场实际情况和勘查资料（见图2）可知，1#风道位于K-3-2及K-3-3中风化砂砾岩层，该区域地质情况与文物鼓楼处较接近，因此将隧道1#风道作为模拟爆破点。此外，区域场地构成上软下硬软硬不均的工作面，相对于单纯的岩体，岩石-土软硬交界面的存在将放大爆破引起的地表峰值振速，也将间接影响文物鼓楼的峰值振速。预先在掌子面掏槽，形成临空面，掘进进尺为0.6 m，钻孔孔深取0.7 m，炮孔孔径为42 mm，炮孔距取0.5 m。试爆分为二类，共进行了6次爆破测试。一类为单孔按顺序爆破，相互之间没有干扰，1～5号孔单孔药量依次为100、200、300、400、500 g，共5种工况；二类为六孔毫秒延时爆破，其中6～10号孔每孔装药200 g，11号孔装药100 g，共1 100 g，各孔间延时50 ms。模拟爆破作业在K-3-2层的中风化砂砾岩中，围岩等级为Ⅱ～Ⅳ的坚硬岩石，如图2所示。

图  2  鼓楼与模拟爆破点的相对位置及地层分布

Figure  2.  Relative position between the drum tower and the blasting site as well as stratum distribution

下载: 全尺寸图片 幻灯片

2.   土体动力本构模型

采用非线性滞回本构模型描述岩土体对爆破振动波的能量耗散作用^[13]。为更好地描述并简化不规则荷载作用下Davidenkov模型的剪应力-剪应变关系曲线，赵丁凤等^[14]引入Pyke的n倍法对加卸载准则进行修正（见图3），当施加的应力转向后，后续的剪应力-剪应变曲线沿当前拐点指向历史上最大（小）点的方向前行，此时循环剪应力$\tau$-剪应变$\gamma$关系由以下公式来描述：

图  3  不规则加卸载准则修正的Davidenkov模型^[14]

Figure  3.  The Davidenkov model modified by the irregular loading-unloading rules^[14]

下载: 全尺寸图片 幻灯片

式中：${G_{\max }} = \rho v_{\rm{s}}^2$为最大剪切模量，$\rho$为土体密度，${v_{\rm{s}}}$为土体剪切波速；${\tau _{\rm{c}}}$和${\gamma _{\rm{c}}}$分别为加卸载转折点处的剪应力和剪应变。

通过对式（1）中的变量($\gamma - {\gamma _{\rm c}}$)求导，得到应力-应变滞回曲线段t时刻的时变切线剪切模量${G^{(t)}}$的表达式：

式中：A、B和${\gamma _0}$为土的试验参数；${\tau _{{\rm{ex}}}}$和${\gamma _{{\rm{ex}}}}$分别为历史上的最大剪应力和剪应变；符号“±”在加载时取“−”，卸载时取“+”。

根据上述修正，除记录历史上最值点外，当施加的应力转向后，只需要记忆当前转折点处的应力和应变，就可以确定应力-应变曲线的走向，通过这样的方法可以有效解决转折点信息记忆量大的问题^[14]。此时，不规则加卸载条件下的应力-应变路径为曲线段 0→1→2→3→4→5→6→7，修正2倍Masing法则中的曲线段2→1→3或曲线段6→5→7，模型参数见表2。

表  2  土体剖面及Davidenkov模型参数

Table  2.  Parameters for the soil profiles and the Davidenkov model

土层描述	厚度/m	重度/(kN·m−3)	剪切波速/(m·s−1)	Davidenkov模型参数
				A	2B	${\gamma _0}$/10−4
杂填土	1.3	18.5	139.5	1.05	0.84	5.5
粉质黏土	6.3	20.2	250.4	1.09	0.82	6.2
砂砾岩	12.2	22.6	558.3	1.30	0.40	21.0

下载: 导出CSV 

| 显示表格

应力-应变滞回曲线段t+Δt时刻的时变体积模量${K^{(t)}}$可由下式确定：

式中：$\mu$为泊松比。

由此，对t时刻的应力张量$\sigma _{ij}^{(t)}$进行更新：

3.   三维数值模型的建立

3.1   毫秒延时爆破荷载输入

采用ABAQUS计算平台中的ConWep模型，可实现基于改进Friedlander方程的爆破冲击波模拟。如图4所示，随着空气冲击波的传播，超压降低，使波振面后的压力衰减至环境压力以下。由改进的Friedlander方程表示爆炸冲击波产生的入射压强${p_{\rm{i}}}(t)$和反射压强${p_{\rm{r}}}(t)$的传播规律^[15]：

图  4  爆炸引起的超压时程曲线

Figure  4.  Time history of overpressure induced by blasting

下载: 全尺寸图片 幻灯片

式中：${p_{{\rm{atm}}}}$为环境压力，${p_{{\rm{i}},\max }} - {p_{{\rm{atm}}}}$为入射超压，t_a为激波前沿到达炮孔表面的时间，t_d为正相超压的持时，a为衰减因子，I为给定位置处的爆破冲量。根据${p_{{\rm{i}},\max }} - {p_{{\rm{atm}}}}$、t_d和I，衰减因子a可根据式(10)进行迭代计算。将炮孔内空气近似描述为理想气体，取$\gamma =1.4$。

入射冲击波作用在炮孔表面时将发生反射，从而导致爆炸冲击波压力和冲量增强。此时，炮孔表面压强$p(t)$为入射压强${p_{\rm{i}}}(t)$和反射压强${p_{\rm{r}}}(t)$的耦合效应，为介质颗粒动量变化速率产生的单位面积上的反作用力，并与冲击波至炮孔表面的入射角度$\theta$有关。$\theta$为炮孔内单元面中心点与爆源连接射线与单元面外法线之间的夹角（见图5），则$p(t)$与入射压强${p_{\rm{i}}}(t)$和反射压强${p_{\rm{r}}}(t)$满足如下经验关系。

图  5  爆炸冲击波至炮孔表面的入射角度

Figure  5.  Illustration of the incident angle of the blast wave to the surface of the blast hole

下载: 全尺寸图片 幻灯片

当$\cos \theta {\text{≥}} 0$时：

当$\cos \theta {\text{＜}}0$时：

ConWep模型参数涉及爆源至迎爆面的距离、装药量和起爆时刻等。相比于直接施加半理论半经验的爆破荷载压力时程曲线，ConWep模型中考虑了爆破冲击波入射压强和反射压强在炮孔表面的空间耦合效应，增加了空间维度。

3.2   地层-爆源体系三维有限元模型

爆破是一个非常复杂的力学过程，现有的有限元、离散元等数值方法难以准确模拟炮孔起爆条件下从炸药爆轰、岩石破碎到弹性地震波激发这一复杂过程，需寻求一种实用的等效模拟方法。本文中采用等效弹性边界的概念^[16-19]，把炮孔周围岩体根据其破坏程度划分为粉碎区、破碎区和弹性振动区，将整个非弹性区(粉碎区和破碎区)等效为爆炸源，由ConWep模块激发的爆炸冲击波直接作用在等效弹性边界上，产生随时间和空间变化的压强。

考虑到圆柱形炮孔能量传播的三维特性，基于ABAQUS平行计算平台，建立了三维有限元精细化模型。炮孔等效弹性边界按照文献[16-19]近似取值r=10r₀ (r₀为实际孔径)。爆源设置在炮孔中心处，由ConWep模块施加在炮孔等效弹性边界上。计算区域为75 m×83 m×20 m。为弥补有限元模型中场地尺寸过大对计算效率的影响，模型边界采用有限元-无限元耦合的方法来模拟，以考虑爆破地震波在边界的透射问题。模型除自由表面及炮孔所在面外，均设置为无限单元，炮孔所在面约束法向位移。计算采用单元总数8 009 728个，具体三维有限元模型如图6所示。

图  6  地层-爆源体系三维有限元模型

Figure  6.  A three-dimensional finite element model of the ground-blasting-source system

下载: 全尺寸图片 幻灯片

4.   数值模拟与实测结果的对比

为初步验证此次建立的模型对大型场地爆破地震波传播规律的预测能力，现对比现场爆破试验与数值模拟得到的地表振动响应。图7对比了地表试验测点（横坐标）与数值模拟（纵坐标）所得的竖直向/水平向峰值振动速度。由图7可知，无论是竖直向峰值振动速度，还是水平向峰值振动速度，两者的实测值与模拟值较均匀地分布在45°线附近。这表明，基于修正的Davidenkov 本构模型以及改进的Friedlander方程的三维精细化有限元模型能较真实地反映爆破试验场地的地表振动特性。

图  7  实测地表测点峰值振速与模拟结果的对比

Figure  7.  Comparison of peak vibration velocities between monitored and simulated surface measurement points

下载: 全尺寸图片 幻灯片

图8对比了测点5 (距爆源水平距离8 m处) 和测点9 (距爆源水平距离40 m处) 地表水平振动速度时程和傅里叶谱的模拟与实测结果。整体而言，模拟与实测的振动速度衰减规律较一致，远处测点的实测速度时程稍有滞后，这是由于实际场地中土的阻尼特性更显著。此外，对比这两个测点速度响应的傅里叶谱可知，模拟与实测速度响应的频谱成分、谱形及谱值基本相同。综上可知，本文模型方法能较理想地模拟爆破荷载作用下场地的环境振动特性，反映爆破地震波在场地传播过程中的耗能特性。

图  8  地表振动速度时程和傅里叶谱对比

Figure  8.  Comparison of vibration waveforms and Fourier spectra bewteen different monitoring points

下载: 全尺寸图片 幻灯片

5.   爆破模式对地表速度响应的影响

由于岩土体的阻尼特性，随着传播距离的增大，爆破地震动的能量在岩土介质不断耗散。此外，场地反应也与爆破地震动的特性密切相关。因此，了解不同爆破模式对场地地表响应的影响机理，有利于爆破方案的合理设计。对此，开展了两种爆破方案的数值计算：(1)六孔齐爆，单孔装药量为500 g；(2) 六孔毫秒延时爆破，单孔装药量为500 g，延时间隔为5 ms。值得一提的是，本文数值模拟为精确多孔齐爆，而目前瞬发工业雷管或同段毫秒雷管难以实现精确的多孔齐爆，即使精度较高的数码电子雷管也存在±1 ms的误差。

地表峰值振速通常随爆心距呈指数衰减：

式中：v为爆破振动质点最大速度（cm/s)；Q为炸药量（kg)，齐爆时取总药量，延时起爆时取单段最大药量；R为爆源距(m)；d为比例距离；K和$\alpha$分别为爆破点至测点间地形、地质条件的相关系数和衰减指数。

图9中给出了两种爆破模式下地表峰值振速随比例距离的变化趋势。由图9可知，在总药量相同的情况下，在相同测点处，延时爆破引起的地表峰值振动速度较齐爆更低，表明延时爆破的降振效果明显。此外，对比图9(a)和图9(b)可以看出，相较于两种爆破模式下的竖向地表峰值振动速度衰减曲线，水平向地表峰值振动速度衰减曲线的差异较大。这表明，延时爆破显著降低了爆破施工引起邻近建筑物水平向的振动响应。

图  9  不同爆破方式下地表振动速度衰减规律对比

Figure  9.  Comparison of ground vibration attenuation under different blasting modes

下载: 全尺寸图片 幻灯片

先起爆炮孔的速度波形和后起爆炮孔的波形相互干扰，导致波形峰值相互错开进而影响频谱特性^[20]。因此，进一步对比两种工况下地表速度响应(距爆源水平距离8 m)的频谱特性。图10所示：齐爆产生的频带分布较广，并在0～200 Hz内出现多个峰值；延时爆破产生的频率成分较单一，主频也较高，约为200 Hz。从图10可以看出，地表土层具有显著吸收爆破振动中高频成分的特性，且对总能量有衰减作用。综上可知，延时爆破对分散爆破振动能量的作用显著，地表速度响应的主频较高，更远离建筑结构自振频率(见表1)，可显著降低爆破施工引起的邻近建筑物的结构振动水平。

图  10  不同爆破方式下地表测点的频谱对比

Figure  10.  Comparison of frequency spectra at surface monitoring points under different blasting modes

下载: 全尺寸图片 幻灯片

6.   结　论

依托南京地铁隧道的爆破振动现场试验，考虑毫秒延时爆破的不规则荷载作用，采用非线性土体本构模型和改进的Friedlander方程，建立了地层-爆源体系三维精细化有限元模型，对比分析了50 ms延时爆破和齐爆引起的环境振动特性，得到的主要结论如下：

（1）采用不规则加卸载准则构造的基于Davidenkov骨架曲线的应力-应变滞回圈描述场地土动力非线性特性，能够较好地模拟毫秒延时爆破不规则荷载作用下的土体动力响应。

（2）通过改进的Friedlander方程来模拟内源爆炸在圆柱形炮孔表面产生的瞬态空气冲击波，实现了包含爆破荷载输入和有限元-无限元耦合边界的地层-爆破荷载体系的三维精细化有限元模型，与现场实测数据对比验证了该模型方法的有效性。

（3）两种爆破方式中，地表峰值振速均随比例距离的增大而呈指数衰减。50 ms延时爆破引起的地表振动速度显著低于齐爆引起的地表振动速度，在地表土层传播中衰减较快，且地震波主频集中在高频段，远离建筑结构自振频率，能够显著降低邻近建筑物的结构振动水平。