深埋直墙拱形地下洞室爆破振动数值模拟分析

Numerical Simulation Analysis of Blasting Vibration in Deep Buried Straight Wall Arch Underground Cavern

摘要: 地下水封油库的地下洞室通常采用钻爆法开挖，爆炸释放出的巨大能量会引起一系列危及工程安全的问题会影响到地下水封油库的正常使用，严重时还会造成人员伤亡和重大的经济损失。本文借助于ANSYS/LS-DYNA动力有限元分析软件来模拟爆破开挖对某地下水封石油洞库1#洞室围岩动力响应的影响，研究内容主要包括爆破开挖引起的洞室围岩质点振动速度分布情况，沿洞室轴向的传播规律。在相同时刻内地下洞室围岩质点各方向的振速呈对称性分布；X方向的最大振速出现在直立墙中部位置，而Y和Z方向最大振速出现在底板中部位置；随着爆心距R的增大，地下洞室截面上X、Y方向的峰值振速分布图形状变化较小，而Z方向的峰值振速分布图形状逐渐趋于圆形。

Abstract: The underground caverns of underground sealed oil depots are usually excavated using the drilling and blasting method. The huge energy released by the explosion can cause a series of problems that endanger engineering safety and affect the normal use of underground sealed oil depots. In severe cases, it can also cause casualties and significant economic losses. This paper uses ANSYS/LS-DYNA dynamic finite element analysis software to simulate the impact of blasting excavation on the dynamic response of the surrounding rock of 1# cavern of an underground water sealed oil cavern. The research content mainly includes the distribution of particle vibration velocity of surrounding rock caused by blasting excavation, and the propagation law along the axial direction of the cavern. At the same time, the vibration velocities of the surrounding rock particles in the underground cavern exhibit a symmetrical distribution in all directions. The maximum vibration velocity in the X direction appears in the middle of the vertical wall, while the maximum vibration velocity in the Y and Z directions appears in the middle of the bottom plate. As the distance between the explosion center R increases, the shape of the peak vibration velocity distribution map in the X and Y directions on the cross-section of the underground chamber changes less, while the shape of the peak vibration velocity distribution map in the Z direction gradually tends to be circular.

1. 引言

在地下工程的爆破施工中，如果爆破工艺设计不当，则爆炸释放出的巨大能量会引起爆心附近一定范围内硐室围岩的松动失稳、原裂隙扩张、新裂隙产生、渗水率增加、岩爆等一系列危及工程安全的问题，严重时还会造成人员伤亡和重大的经济损失。近年来，有关学者对爆破地震波传播过程的振动衰减规律展开了研究，并取得了进展。李新平等将损伤变量引入弹塑性本构方程中，通过数值模拟方法及现场爆破振动测试，研究了爆破损伤范围及判据 [1] ；曹峰等采用LS-DYNA有限元软件建立三维模型，基于HJC本构模型，通过引入损伤变量，模拟循环荷载作用下的累积损伤演化过程，得到累积损伤程度、损伤范围与爆破次数之间的关系 [2] ；王登科等基于叠加原理和应力波衰减理论，优化爆破等效荷载在隧道群孔中的施加方法，提出大型地下洞室爆破等效荷载计算模型，实现地下洞室群爆破开挖时围岩动力响应的快速计算 [3] ；L. F. Trivino等在研究多种爆源不同起爆条件下的能量和频率的变化规律时引入对平均频率的分析 [4] ；卢文波等通过引入介质阻尼项，建立了球形药包爆破条件下实际岩体介质中爆破振动的频谱表达式，进而分析爆破振动频率的衰减机制和影响因素 [5] ；周俊汝等通过数值分析和现场监测的方法分析黏性岩体中地震波的衰减机制，得到球状药包振动主频及平均频率衰减规律 [6] ；费鸿禄等对振动信号进行小波包能量谱分析，得到爆破振动信号能量在各频带上的分布，在沿隧道掘进方向，随着与掌子面之间距离的增大地震波的主频下降 [7] 。

本文借助于ANSYS/LS-DYNA动力有限元分析软件来模拟爆破开挖对某地下水封石油洞库1^#洞室围岩动力响应的影响。研究内容主要包括：爆破开挖引起的洞室围岩质点振动速度分布情况，沿洞室轴向的传播规律等。

2. 有限元模型

为分析爆破开挖对本洞室围岩质点动力响应的影响，假设地下洞室围岩为各向同性介质。以炸药中心位置为原心，X方向为水平方向指向洞室右侧边墙，Y方向为竖直向上方向，Z方向为洞室轴线方向指向掌子面方向；数值模型采用m-kg-s单位制。模型整体尺寸为100 m × 100 m × 250 m，由岩体、空气和炸药三种材料组成。有限元数值模型如图1所示。

Figure 1. Finite element model

图1. 计算模型

在计算过程中对模型的外边界面均施加相应法向位移约束；另外，由于实际洞室为深埋洞室，而且该模型只是无限岩体的一部分，为消除人为边界面的反射波对结构动力响应的影响，将模型的外边界面均设为无反射边界面。

根据现场工程地质勘察报告资料，场地内深部岩体较完整，裂隙稀疏，选取LS-DYNA中的塑性随动硬化材料模型作为围岩本构模型。塑性随动硬化材料模型可考虑与材料的应变率相关的失效作用，其应变率用Cowper-Symonds模型来考虑，其本构关系和材料屈服应力σ_y (Pa)如公式(1)所示。

${\sigma }_y=\left[1+{\left(\frac{\epsilon }{C}\right)}^{\frac{1}{P}}\right]\left({\sigma }_0+\beta \frac{E_{\tan }E}{E-E_{\tan }}{\epsilon }_p^{eff}\right)$ (1)

式中：ε为岩体的应变率，无量纲；

C和P为Cowper-Symonds模型中应变率参数，无量纲；

σ₀为岩体的初始屈服应力，Pa；

E为岩体的弹性模量，Pa；

E_tan为岩体的切线模量，Pa；

${\epsilon }_p^{eff}$ 为岩体的有效塑性应变，无量纲；

β为模型的硬化参数，β = 0~1，0表示模型仅考虑随动硬化、1表示仅模型考虑各向同性硬化。

岩体模型的物理力学参数根据原位岩体试验和室内岩样试验计算求得，具体数值如表1所示。

Table 1. Parameters of rock

表1. 岩体材料参数

在爆炸问题的数值模拟中，由于炸药爆炸产生的爆轰产物压力P波动范围很大，一般在0.1~1000 MPa，一般材料模型的状态方程很难适合这个范围。本位炸药材料选用LS-DYNA中的内嵌模型来模拟，其中爆轰压力P与单位体积内能E₀、相对体积V的关系采用JWL状态方程来描述。JWL状态方程是通过圆桶爆破实验，得到爆轰产物等熵线随炸药初始密度和爆速而变化的实验数据，通过热力学原理计算得到相应的函数关系。JWL状态方程具有精确的物理意义，可以对爆炸过程中产生的压力做出与实验结果非常相近的预测，因而在爆炸问题的数值模拟中得到了广泛应用。该模型在模拟炸药爆炸过程时，爆轰产物的压力P (Pa)与相对体积V和单位体积初始能量E₀的函数方程如式(2)所示：

$P=A\left(1-\frac{\omega }{R_{1}V}\right){\text{e}}^{-R_{1}V}+B\left(1-\frac{\omega }{R_{2}V}\right){\text{e}}^{-R_{2}V}+\frac{\omega E_0}{V}$ (2)

式中： $A,B$ 均为炸药材料的有关参数，Pa；

$R_1,R_2,\omega$ 表示炸药材料的相关常数，无量纲；

V为单位体积炸药产生的爆轰产物体积，无量纲；

E₀为炸药的初始单位体积内能，J/m³。

实际过程中所用的岩石乳化炸药的材料参数和状态方程参数见表2所示。

Table 2. Parameters of rock emulsion explosive

表2. 岩石乳化炸药材料参数

为能够真实反映爆破振动作用下洞室围岩质点的动力响应，需要在地下洞室开挖部分建立空气材料模型。在LS-DYNA中，通常采用MAT_NULL材料模型来模拟空气材料，空气材料的压力值P (Pa)表示为空气的单位初始体积内能E₀的线性关系，由公式(3)给出：

$P=C_0+C_1\mu +C_2{\mu }^2+C_3{\mu }^3+\left(C_4+C_5\mu +C_6{\mu }^2\right)E_0$ (3)

式中： $\mu =\rho /{\rho }_0$ ；

ρ为当前密度，kg/m³；

ρ₀为初始密度，kg/m³；

为材料的单位初始体积内能，J/m³；

C₀~C₆为状态方程参数，无量纲。

该模型中的空气材料参数见表3所示。

Table 3. Parameters of air

表3. 空气材料参数

在ANSYS/LS-DYNA中可供模拟岩土体中炸药爆炸的数值计算方法有：共用节点算法、面面接触算法、侵彻接触算法和流固耦合算法等。文献在对比不同算法的计算结果之后，得出采用共用节点算法与滑动接触算法获得的结果与爆破实际比较符合，同时这两种算法获得的计算结果也比较一致。事实上每一种计算方法中都有多个控制选项，其参数的不同取值直接影响到计算结果。因此对于一个具体问题，应该以实验为基础，微调参数使计算结果与理论解和试验值等相吻合。为节省计算时间，本文选用共用节点算法。

3. 数值模拟准确性验证

部分现场监测结果与对应位置的数值模拟结果见表4所示。数值模拟与监测试验的平均绝对误差为0.03 cm/s、平均相对误差为9.69% < 10.00%。由于实测V_max不是在同次爆破振动作用下测试得到的，离散性较大，导致平均相对误差较高，但从平均绝对误差的角度考虑，该数值模拟的准确度是可以接受的。

Table 4. Result comparison of numerical simulation and monitoring test

表4. 数值模拟与监测试验结果对比表

4. 结果分析

4.1. 轴向峰值振速分布规律

数值模拟得到的在单响炸药量Q = 64.8 kg的爆破振动作用下，不同隧道断面不同位置不同方向质点的峰值振速如表5所示。

Table 5. Peak velocity of surrounding rock masses

表5. 质点峰值振速

利用数值模拟结果，根据《爆破安全规程》中规定的M.A.萨道夫斯基经验公式，回归得到围岩质点峰值振速在洞室轴向上的传播衰减规律如表6所示。

Table 6. Attenuation law of peak velocity in longitudinal direction

表6. 峰值振速洞室轴向衰减规律

由衰减指数α的大小可以看出：① 在拱顶和底板中部位置，Z方向的峰值振速衰减速率最大，衰减指数α分别为1.25和1.20；而X方向的峰值振速衰减速率最小，其衰减指数α分别为0.88和0.99。② 在左拱角位置，Y方向的峰值振速衰减速率最大，其衰减指数α为1.22；而X方向的峰值振速衰减速率最小，其衰减指数α为0.98。③ 在左直立墙中部位置，Z方向的峰值振速衰减速率最大，其衰减指数α为1.19；而Y方向的峰值振速衰减速率最小，其衰减指数α为0.89。④ 在墙角位置，X方向的峰值振速衰减速率最大，其衰减指数α为1.20；而Y方向的峰值振速衰减速率最小，其衰减指数α为1.01。

4.2. 截面峰值振速分布规律

取爆心距R = 10.0 m处不同位置不同方向的截面峰值振速，作隧洞截面上水平方向(X方向)、铅直方向(Y方向)、洞轴线方向(Z方向)、隧道面法向(N方向)的峰值振速分布图，如图2所示。截面峰值振速分布图均采用极坐标系，坐标原点为直墙拱形隧道断面的矩形区域中心。可以看出：在爆心距R = 10.0 m处，X方向上的最大振速为17.76 cm/s，出现在地下洞室的直立墙中部位置；Y方向上的最大振速为23.91 cm/s，出现在地下洞室的底板中部位置；Z方向上的最大振速为27.69 cm/s，同样出现在地下洞室的底板中部位置，这一点的基本规律与文献是一致的。

Figure 2. Peak velocity distribution at R = 10.0 m section

图2. R = 10.0 m截面峰值振速分布图

隧道面法向(N方向)峰值振速定义为：在拱顶和地板处为质点的铅直方向(Y方向)峰值振速；在直立墙部位为质点的水平方向(X方向)峰值振速；在拱脚和墙角处为质点X方向和Y方向峰值振速的算数平方根。可以看出：拱脚位置的法向振速小于拱顶位置，而墙角位置的法向振速小于直立墙中部以及底板中部位置的法向振速。

Figure 3. Peak velocity distribution at R = 20.0 m section

图3. R = 20.0 m截面峰值振速分布图

爆心距R = 20.0 m处不同位置不同方向的峰值振速分布图如图3所示。与爆心距R = 10.0 m处的截面峰值振速分布图对比可知：随着爆心距R的增大，地下洞室截面上X、Y和N方向的峰值振速分布图形状变化较小，而Z方向的峰值振速分布图形状逐渐趋于圆形。

5. 结论

本文通过数值模拟方法，分析了某地下水封石油洞库1^#洞室围岩质点在爆破振动作用下的动力响应规律。具体如下所述：

① 从振动速度的分布云图来看，在相同时刻内地下洞室围岩质点X、Y、Z方向的振速对X = 0平面呈对称性分布；X方向的最大振速出现在直立墙中部位置，而Y和Z方向最大振速出现在底板中部位置。

② 根据M.A.萨道夫斯基经验公式，回归得到地下洞室围岩质点峰值振速在洞室轴向上的传播规律，由衰减指数α的大小可以看出：在拱顶和底板中部位置，Z方向的峰值振速衰减速率最大，衰减指数α分别为1.25和1.20；而X方向的峰值振速衰减速率最小，其衰减指数α分别为0.88和0.99。在左拱角位置，Y方向的峰值振速衰减速率最大，其衰减指数α为1.22；而X方向的峰值振速衰减速率最小，其衰减指数α为0.98。在左直立墙中部位置，Z方向的峰值振速衰减速率最大，其衰减指数α为1.19；而Y方向的峰值振速衰减速率最小，其衰减指数α为0.89。在墙角位置，X方向的峰值振速衰减速率最大，其衰减指数α为1.20；而Y方向的峰值振速衰减速率最小，其衰减指数α为1.01。

③ 在单响炸药量Q = 64.8 kg的爆破振动作用下，1^#洞室围岩峰值振速在爆心距R = 10.0 m截面处，X方向上的最大振速为17.76 cm/s，出现在地下洞室的直立墙中部位置；Y方向上的最大振速为23.91 cm/s，出现在地下洞室的底板中部位置；Z方向上的最大振速为27.69 cm/s，同样出现在地下洞室的底板中部位置。另外随着爆心距R的增大，地下洞室截面上X、Y和N方向的峰值振速分布图形状变化较小，而Z方向的峰值振速分布图形状逐渐趋于圆形。