非均质岩体中一维应力波演化过程分析

摘要

地下球形药包爆轰后将产生很强的地震波，其传播受到多种因素的影响。把岩石介质的弹塑性帽盖模型耦合到有限差分法程序中，对一维爆轰波在非均质岩石体中的传播、演化机理进行了数值模拟与分析；着重探讨了一维应力波从“硬”岩石层进入“软”岩石层、从“软”岩石层进入“硬”岩石层时，波在形状和强度等方面的变化规律；同时探讨了岩石介质中存在的“软弱夹层”对一维应力波传播的影响。

正文

《1 前言》

1 前言

地下爆炸是民用与国防工程中常用的破岩手段。矿山建设、道路施工、地下工事的修建等均需要开展地下爆炸效应的研究, 地下爆炸波的传播规律和演化机理是地下爆炸力学的核心问题 ^[1]。

20世纪60年代初, 美国的Livermore国立实验室为配合地下核爆炸试验, 组织研究了球对称爆破效应的一维计算;我国几乎在同一时期开始这方面的研究。张景山等曾对我国的地下核爆炸应力波传播、演化进行了一维数值计算 1;工程兵三所也做过这方面的工作。

岩石爆轰波的传播、演化过程是极其复杂的, 其影响因素较多, 如地形、地质条件以及物理力学参数、爆破的种类和方法、爆源的大小和形状等。笔者重点探讨来自岩石强度 (软、硬) 等方面的因素的影响, 这对人防工程的防震、减震等都具有积极意义。

《2 有限差分法基本方程》

2 有限差分法基本方程

地下强爆炸的初始反应阶段十分复杂, 超出了力学分析的范围。因此, 在数值计算中, 一般以该阶段结束的瞬间作为初始时刻, 把引起的力学效应等效为一个简化的爆炸源模型, 然后基于合理的岩石本构模型, 对第二阶段 (波动力学阶段) 开展数值模拟。其中, 球活塞爆炸源模型计算起来比较简单, 其假定药室空腔内的爆生气体作等熵膨胀, 推动岩石介质向外运动, 腔壁上的压力就是计算的边界条件 ^[1]。

从严格意义上说, 地下爆炸是二维或三维波传播。如果暂不考虑埋深和地表反射的影响, 可以简化成一维球对称问题, 用Lagrange坐标表示的岩石介质的基本控制方程为 ^[1,2,3]

质量守恒:

$v = v_{0} \frac{r^{2}}{R^{2}} \frac{\partial r}{\partial R} (1)$ $v = v_{0} \frac{r^{2}}{R^{2}} \frac{\partial r}{\partial R} (1)$

动量守恒:

$\begin{array}{l} ρ_{0} \frac{\partial u}{\partial t} = \frac{r^{2}}{R^{2}} \frac{\partial [- (p + q) + s_{r}]}{\partial R} + \\ \frac{(s_{r} - s_{θ})}{R} \frac{r}{R} \frac{\partial r}{\partial R} (2) \end{array}$ $\begin{array}{l} ρ_{0} \frac{\partial u}{\partial t} = \frac{r^{2}}{R^{2}} \frac{\partial [- (p + q) + s_{r}]}{\partial R} + \\ \frac{(s_{r} - s_{θ})}{R} \frac{r}{R} \frac{\partial r}{\partial R} (2) \end{array}$

速度定义:

$u = \frac{\partial r}{\partial t} (3)$ $u = \frac{\partial r}{\partial t} (3)$

人工粘性:

$q = {\begin{matrix} \frac{b^{2}}{v} (Δ u)^{2} & \frac{\partial u}{\partial R} < 0 \\ 0 & \frac{\partial u}{\partial R} \geq 0 \end{matrix} (4)$ $q = {\begin{matrix} \frac{b^{2}}{v} (Δ u)^{2} & \frac{\partial u}{\partial R} < 0 \\ 0 & \frac{\partial u}{\partial R} \geq 0 \end{matrix} (4)$

式中u表示质点速度;ρ₀代表初始密度;v₀, v分别表示介质的初始与当前比容;r, R分别为欧氏坐标与L氏坐标;b是人工粘性系数;s_r, s_θ分别表示径向和环向的应力偏量;p为静水压力。

数值模拟中采用的是Lagrange坐标下中心交错差分格式 ^[2,3], 其中u定义在空间整格点、时间半格点上;r定义在时间与空间整格点上;其他量均定义在时间整格点、空间半格点上。

《3 岩石帽盖模型》

3 岩石帽盖模型

岩石爆破模型研究是爆破理论和爆破技术发展的关键, 因此人们十分关注岩石本构的建立。一般认为, 帽子模型属于经典的加工硬化、塑性增量理论的范畴, 在每一个荷载增量作用下, 可以承受较小的弹性和塑性应变。在数学上, 总应变增量可以表示为弹性应变和塑性应变之和 ^[4]。从整体而言, 该模型能较全面地反映岩土类特殊介质的应力-应变关系。在动荷载的数值模拟中, 它是基于典型的塑性增量理论的连续介质模型, 已被用来表示许多岩土材料在高、低应力状态下的力学特性, 可较好地考虑到岩体地下强爆炸的弹塑性特点。

图1是一种典型的帽盖示意图, 由两部分组成:第一部分是理想的塑性固定失效包络面f₁;第二部分是帽盖屈服面f₂, 随着硬化参量的增加而增大。该模型是一种塑性硬化效应的本构模型。设图1中椭圆与包络面的切点坐标为 $(J_{1 s}, \sqrt{J_{2 s}})$ $(J_{1 s}, \sqrt{J_{2 s}})$ 。此模型屈服函数的一般形式为 ^[1,5]:

${\begin{cases} f_{1} (J_{1}, \sqrt{J_{0}}) = 0 包络面 \\ f_{2} (J_{1}, \sqrt{J_{2}}, Η) = 0 帽盖 \end{cases} (5)$ ${\begin{cases} f_{1} (J_{1}, \sqrt{J_{0}}) = 0 包络面 \\ f_{2} (J_{1}, \sqrt{J_{2}}, Η) = 0 帽盖 \end{cases} (5)$

《图1》

图1 帽盖模型示意图

Fig.1 Plot of cap model

其中, J₁, J₂为应力偏量第一、第二不变量;H是硬化参量。为了简化起见, 取帽盖为一簇 $(J_{1}, \sqrt{J_{2}})$ $(J_{1}, \sqrt{J_{2}})$ 平面上与包络面相切的椭圆;若椭圆的中心 (J_c, 0) 在J₁轴上, 竖轴长度为 $\sqrt{Q}$ $\sqrt{Q}$ , 水平轴与竖轴之比为R, 则帽盖椭圆的方程式为:

$\frac{J_{2}}{Q} + \frac{(J_{1} - J_{c})^{2}}{Q R^{2}} = 1 (6)$ $\frac{J_{2}}{Q} + \frac{(J_{1} - J_{c})^{2}}{Q R^{2}} = 1 (6)$

笔者拟采用Dimaggio和Sandler等提出的具体帽盖模型屈服面方程 ^[5,6]:

$\begin{array}{l} f_{1} = J_{2} - (A - C e^{- B J_{1}})^{2} (7) \\ f_{2} = J_{2} - J_{2 s} + \frac{(J_{1} - J_{1 s})^{2}}{R^{2}} - \\ 2 B C (J_{1} - J_{1 s}) e^{- B J_{1 s}} (A - C e^{- B J_{1 s}}) (8) \end{array}$ $\begin{array}{l} f_{1} = J_{2} - (A - C e^{- B J_{1}})^{2} (7) \\ f_{2} = J_{2} - J_{2 s} + \frac{(J_{1} - J_{1 s})^{2}}{R^{2}} - \\ 2 B C (J_{1} - J_{1 s}) e^{- B J_{1 s}} (A - C e^{- B J_{1 s}}) (8) \end{array}$

帽盖模型长、短轴半径之比为:

$R = R_{0} e^{R_{1} D Η} (9)$ $R = R_{0} e^{R_{1} D Η} (9)$

假定体模量与静水压的关系为:

$Κ = Κ_{0} (1 - Κ_{1} e^{- Κ_{2} J_{1}}) (10)$ $Κ = Κ_{0} (1 - Κ_{1} e^{- Κ_{2} J_{1}}) (10)$

剪切模量G被看成是常数;包络面与帽盖切点的横坐标J_1s是硬化参数H的线性函数J_1s=DH。如此, 模型参数共有10个:A, B, C, D, K₀, K₁, K₂, G₀, R₀和R₁。当帽盖模型取以上形式时, 可推导出包络面与帽盖上的塑性流动因子dλ具体形式 (略) ^[1,5,6], 据此即可计算出塑性应变、应力、偏应力等。

《4 数值计算过程》

4 数值计算过程

假定炸药完全填塞药室, 炸药的爆轰过程是稳定的, 并忽略起爆过程的化学反应, 如此可对从球心开始的起爆过程进行数值计算。第一步采用C-J爆轰条件, 爆炸产物采用BKW方程, 以形成初始条件;第二步将起爆过程结束时的爆炸产物的压力、速度分布作为初始条件, 对爆轰波的传播进行计算。第一步生成初始条件的过程参见文献[1,7,8]。

图2中带充填的小圆代表药室, 里面填满炸药, 外围是岩石层, 虚线表示岩石分界线。这样, 图2中共有2层岩石Ⅰ和Ⅱ。炸药从球心起爆后, 波阵面后方是已完全反应后的气体产物。计算时认为爆炸产物的外边界是以恒定的速度D_cj向外传播的爆轰波波阵面, 因此有波阵面条件:

《图2》

;波阵面后方的气体产物则按一维球对称的气体基本方程进行计算。装药密度ρ₀等于1 680.0 kg/m³, 爆速D_cj为7 830.0 m/s, 爆压p_cj为27.0 GPa。爆轰波后方气体产物的BKW方程参数如表1所示 ^[1,2]。

《图3》

图2 一维球面波传播示意图

Fig.2 Propagation model of 1-D spherical wave

表1 爆炸产物的BKW方程参数

Table 1 Parameters for BKW equation of explosive products

《表1》

a₀	a₁	a₂	a₃	a₄
-3.469	-2.488	0.246 7	0.0294 7	-0.012 01

药室中充填的炸药为1 350 kg, 若充填密实且无空隙, 则可算出球形药室的半径为0.58 m。岩石帽盖模型参数列于表2 ^[5,6]。

表2 岩石帽盖模型参数

Table 2 Parameters of cap model for rock mass

《表2》

A/Pa	B/Pa^-1	C/Pa	D	G/Pa	ρ₀/kg·m^-3	R₀	R₁/Pa^-1	K₀/Pa	K₁	K₂/Pa^-1
1.047×10⁹	4.209×10^-10	9.991×10⁸	450	2.274×10¹⁰	2 500.0	4.0	7.257×10^-9	5.618×10¹⁰	0.75	2.903×10^-9

笔者主要讨论非均匀岩石介质对波的影响。当应力波从“硬”岩石层进入“软”岩石层时, 为了进行模拟, 在图2中设Ⅰ代表“硬”岩层, Ⅱ代表“软”岩层。在模型参数的选取上, 对于“硬”岩层, 其帽盖参数取表2提供的数值;对于 “软”岩石介质, 由于缺乏试验数据, 做了近似处理——把表2中的参数G, K₀分别取原值的0.1倍, 密度ρ₀取原值的0.8倍, 其他参数保持不变, 以此作为“软"岩层的帽盖模型参数。反之, 当应力波从“软”岩层进入“硬”岩石层时, Ⅰ代表“软”岩层, Ⅱ代表“硬”岩层, 其参数取值与上类似。第三种情形是考虑岩石中的“软弱夹层”对爆轰波传播的影响。中间的“软弱夹层”模型10个参数的具体取值与上面两种情形的“软”岩层参数取值可作同样考虑, “软弱夹层”两边的岩石设为相同的“硬”岩层, 模型参数即为表2中提供的数据。计算中每个独立岩层都假定为连续的、各向同性介质。

程序中一维计算域的最大半径取13.5 m。对于第一、第二种情形, 岩层Ⅰ、Ⅱ分界面的半径为5.0 m;对于第三种有“软弱夹层”的情形, 设定“软弱夹层”厚0.6 m (“软弱夹层”圆环的内、外径分别为3.5和4.1 m) 。图3至图7给出了3种不同情形的计算结果曲线, 其中虚线代表药室外围为同一均质岩层 (模型参数取表2数据) 的计算结果, 它是用来进行对比的, 并且与非均质岩层计算曲线一一对应。

第一种情形——应力波从“软”岩层到“硬”岩层 (即在图2中Ⅰ代表“软”岩层、Ⅱ代表“硬”岩层) 。从图3的应力波形和图4的速度波形均能看出, 当爆轰波从“软”岩层传播到“硬”岩层时, 应力波和速度波的波阵面峰值与均质的岩层相比, 均有所提高。同时, 由于“软”岩层的存在, 应力波形和速度波形在径向出现的位置都要小于均质岩层, 这是由于“软”岩层的剪切模量与体积模量较小、波速也小的缘故。

《图4》

图3 从软岩层到硬岩层的应力波形

Fig.3 Stress-waves from soft rock to hard rock

《图5》

图4 从软岩层到硬岩层的速度波形

Fig.4 Velocity-waves from soft rock to hard rock

《图6》

图5 从硬岩层到软岩层的应力波形

Fig.5 Stress-waves from hard rock to soft rock

第二种情形, 即爆轰波从“硬”岩层到“软”岩层的径向应力波形, 虚线是均质岩层的, 实线与虚线在时刻上相对应, 图5中共有3个不同的时刻。当应力波碰到“软”岩层时, 岩石材料的强度指标偏低, 于是变形增大, 爆生气体能量被释放出来, 与相同时刻的均质岩层相比, 压力有所减小。而对于非均质岩层, 相同时刻应力波的传播距离则要小于均质岩层的。因为在“软”岩层中波速偏小。

存在“软弱夹层”的情形。“软弱夹层”位于3.5～4.1 m之间, 设“软弱夹层”两侧为“硬”岩层, “软弱夹层”模型参数取值与上述“软”岩层相似。从图6可见, “软弱夹层”的存在使得波阵面上的应力波峰值有所降低, 传播的速度也较均质岩层的慢。其次, “软弱夹层”还能影响波阵面后方的波形曲线, 在“软弱夹层”附近, 应力波形发生跃动, 经反射波的干涉、叠加, 应力值幅度也会加大。可见“软弱夹层”的存在能起到削弱波峰、吸收能量的作用。

《图7》

图6 存在软夹层的不同时刻应力波形

Fig.6 Stress-waves of different time when soft interlayer existed

图7是为径向距离4.0 m处 (位于“软弱夹层”内) 的定点波形, 可以看出, “软弱夹层”能大大降低应力峰值, 从而达到削弱爆轰波所携带的能量。

《图8》

图7 存在软夹层的定点应力波形

Fig.7 Stress-waves at special point when soft interlayer existed

《5 结语》

5 结语

岩石是一种复杂性的材料, 爆破地震波在岩石介质中传播的影响因素众多, 笔者主要针对岩石强度指标上的差异对应力波传播的影响进行了初步计算与分析。

爆轰波从“硬”岩层进入“软”岩层时, 应力波的峰值明显降低;从“软”岩层进入“硬”岩层时, 应力、速度波的波阵面峰值有所提高;“软弱夹层”的存在或岩石分化带的存在, 使应力波的传播速度明显降低, 应力波的峰值被大大地削弱, 波形也被改变, 在一定程度上可以起到缓冲、吸能与减震的作用。

在弹塑性帽盖模型中, G, K₀是反应岩石强度的直接指标, 在程序计算中的“软”岩层参数选取上主要调整该2个参数, 同时适当降低密度ρ₀, 尽管B, R₀和R₁等也有重要影响, 暂不计入。

岩石属于典型的复杂介质, 对于其细观结构 (裂纹、微孔洞) 引起的材料初始损伤和爆炸过程中产生的损伤, 以及岩石在高温、高压下的粘-弹-塑性效应还有待深入研究。

《注释》

注释

1 张景山. 地下核爆炸一维球对称数值计算总结, 西北核技术研究所.1976

展示更多