爆炸实验水池防护性能及动力学响应分析

胡亚峰; 金建峰; 顾文彬; 陈亮; 张瑞江

doi:10.11883/1001-1455(2017)06-1001-09

摘要: 针对爆炸实验水池强度设计问题，利用非线性动力学程序LS-DYNA对10 kg TNT爆炸后的水中冲击波传播规律及爆炸水池结构动态响应情况进行了数值模拟, 对空气桶和气泡帷幕削弱水中冲击波的能力进行了定量计算，结果表明：空气桶对冲击波峰值压力削弱作用接近50%，对比冲量削弱作用达到16.2%；气泡帷幕对壁面反射冲击波的削弱作用高达86.2%，对比冲量削弱作用达到75.6%。在此基础上进一步分析了冲击载荷作用以及结构响应机理，指出了内衬钢板和混凝土围堰的危险区域，对爆炸能量在各物质间的分配和传递规律进行了初步探索，为相关爆炸水池的工程设计提供参考。

关键词:

Abstract: In this work, in view of the design issues concerning the strength of the explosion testing pool, we simulated numerically the shock wave propagation in water and structural dynamic responses of the testing pool subjected to a 10 kg TNT explosion impact loading using the nonlinear dynamics program LS-DYNA. We also quantitatively calculated the capability of the air tank and the bubble curtain to weaken the shock wave in water. The results show that the weakening effect of the air tank on the shock wave peak pressure and the specific impulse is close to 50% and 16.2%, that of the bubble curtain on the shock wave reflection and the specific impulse is as high as 86.2% and 75.6%. Based on this, we further analyzed the mechanism of the impact loading and the structural response and carried out a preliminary investigation of the distribution and transmission of the explosion energy between each substance. Our work can be used as reference for the engineering design of similar explosion testing pools.

Key words:

构造地震很少由新剪切断裂的产生而传播，而是沿着已有断层边界面突然滑动扩散，地震是板块粘滑摩擦失稳的结果。此过程中脆性断裂只起次要作用，起主要作用的是断裂面和摩擦磨损的交替发展^[1-3]。地震过程由表征断层短暂运动状态的“滑”和表征断层长期准静状态的“粘”组成，后者是地震事件间弹性应变积累的时期。地震过程中岩石界面动摩擦的变形机理主要表现为^[4]在临界温度之上的晶体塑性变形和碎裂流动。虽然都显示出较强的塑性特征，但两种变形有不同的流变学机理。从动力学角度，岩石界面动摩擦过程大致可以分为波致摩擦、惯性摩擦以及滑移三个阶段^[5-7]。本文中进行岩石界面动摩擦过程中剪切失稳行为的研究，主要关注与局部大变形相关的地震事件的“粘”这一过程的形成和发展，贯穿于波致摩擦和惯性摩擦两个阶段，但仍属于界面动摩擦的初期。相关研究可增强对孕震过程的认识。

对岩石界面动摩擦性能的实验研究较多。在低速冲击下，Di Toro等^[8]采用一种伺服控制压-扭装置进行了一系列的快速滑动实验；Rice^[9]基于旋转圆盘研究了动摩擦因数与滑动速度之间的关系，发现在滑动速度为4 m/s时界面出现了零摩擦强度的现象；张磊等^[10-11]基于SHPB（split Hopkinson pressure bar）杆束技术研究了含斜节理面岩石的动摩擦行为；Zhou等^[12]基于SHPB斜截面垫块技术研究了岩石的压剪动力学特性，发现在此过程中，岩石界面滑移速度可以控制在1 m/s量级，相对滑移位移可以控制在毫米量级。Okada等^[13]采用平板斜撞击动摩擦实验，得到了压剪联合高速冲击下界面滑移时的临界滑动速度；Xu等^[14]利用轻气炮压剪联合冲击加载装置对具有斜节理的花岗岩进行了界面动滑移特性研究，此过程的岩石界面滑移速度可以控制在10 m/s量级，相对滑移位移可以控制在0.1 mm量级。当岩石界面发生动态滑移时，岩石界面细微观结构必然发生不可逆变化，如局部损伤、局部断裂等，会产生一定频率的声发射信息和红外辐射信息，对这些信息的追踪分析可得到局部结构发生变化的位置和规模。因此，研究这些物理信息可以推测伴随岩石界面动摩擦过程发展的物理现象，主要包括：(1) 波致摩擦过程中界面Ⅱ型裂纹表面的亚瑞利波、超剪切波的传播^[15-18]，粗糙界面上跨声速的波、亚瑞利波，以及一种慢波的传播^[19-20]，这两组波有较大的差异，但在各自物理现象的描述中均有较好的效果，其机制目前尚不明确；(2) 波致摩擦过程中的红外温升和声发射信息^[6,21-23]，前者反映岩石界面滑移破裂时的规模，后者可以对破裂区域进行定位。岩石界面动摩擦是一个典型的失稳过程，界面作用非常复杂。Huang等^[24]采用多尺度PFC离散元方法跟踪了岩石界面动摩擦因数的演化；Gao等^[25]采用有限元与离散元结合的方法分析了岩石界面的破碎特性。以上研究关注的都是界面破碎带来动摩擦因数的降低及相关机理，但是此过程是如何形成的，则需要探讨岩石界面动态变形发展与界面失稳的联系，尤其是与剪切失稳扩散的关系。

对岩石失稳现象的研究已有较多^[26]。基于线弹性和弹塑性小应变的理论一般不会产生本构失稳现象^[27]。产生失稳现象的情况主要基于下述三种理论：(1)非关联流动理论^[28]，此理论认为体积应变与偏量应变具有不同的流动规律，假定界面失稳，应用弹塑性理论进行分析，由此研究平面应变条件下失稳特性与岩石扩容角和内摩擦角的关系，此理论可以较好地应用于脆性剪切带分析；(2)应变梯度塑性理论^[29]，此理论认为材料力学行为不仅与应变、应变历史有关，而且与应变梯度有关，由此采用屈服函数的应变梯度研究岩石的局部化变形失稳问题；(3)有限变形理论^[30-33]，此理论采用有限变形表述，导出一个简洁的四阶多项式判别方程，以确定产生局部化变形失稳的临界应力状态，结合有限群表征，可以得到平面压缩^[34]和剪切^[35]条件下岩石失稳发展过程，结合扩散型流函数分析，进一步则可以分析失稳时速率扰动的局部分布形状^[36]，主要适用于中等或较低强度岩石在围压作用下剪切带局部化变形失稳分析。

上述失稳分析主要基于准静态加载情况，产生剪切变形带等特征失稳变形的边界条件仍采用Hill等^[30]的零应力边界，这对于动力加载情况不再适合。在岩石界面动态剪切变形的过程中，局部存在复杂的变形特征^[11]。因此，本文中在有限变形理论的基础上，引入惯性力和剪切边界条件的影响，研究岩石动摩擦过程中的剪切扩散行为，结合实验观测和有限元分析探讨岩石界面动态剪切失稳的变形特征，以及其对波动传播的影响。

1. 岩石动态剪切失稳理论分析

1.1 基于有限变形的动态剪切失稳控制方程

在体积为V(t)、边界为∂V的岩石界面系统中，总能量为 $\bar E$ 的位移表达式为^[37-38]：

$\bar E({u_i}) = {E_{\text{int} }} + {E_{{\text{kin}}}} + {E_{{\text{ext}}}}$

(1)

式中：E_int为系统内能，E_kin为系统动能，E_ext为系统外力势能，其表达式分别为：

${E_{\rm{int} }} = \int {W({\varepsilon _{ij}})} {\text{d}}V$

(2)

${E_{{\text{kin}}}} = \int {\rho v_i^2} {\text{d}}V$

(3)

${E_{{\text{ext}}}} = - \int_V {\rho {b_i}} {u_i}{\text{d}}V - \int_{\partial V} {{f_i}} {u_i}{\text{d}}{\varGamma }$

(4)

式中： $W({\varepsilon _{ij}})$ 为应变能密度， ${v_i}$ 为系统变形过程中的速度， ${b_i}$ 为系统的体力分量，u_i为系统变形过程中的位移， ${f_i}$ 为作用于系统边界的面力分量，i=1, 2。一般情况下，忽略体力分量 ${b_i}$ 。图1所示为剪切加载下，岩石局部剪切变形带发展示意图，图1(b)中n为剪切失稳面的外法线方向，α为x₁轴与n方向的夹角，Q代表局部失稳位置。图1(b)中 $\rho {a_n}$ 为法向惯性力， $\rho {a_\tau }$ 为切向惯性力。由于岩石微结构的影响，局部变形不会沿着加载路径，而是会产生分叉变形，形成剪切变形带。

图 1 剪切局部化失稳模型

Figure 1. Shear localized instability model

下载: 全尺寸图片幻灯片

应变能密度W的表达式为：

$W = \frac{1}{2}{L_{ijkl}}{\varepsilon _{ij}}{\varepsilon _{kl}}$

(5)

式中： ${L_{ijkl}}$ 为刚度张量，相应的本构关系满足 ${\sigma _{ij}} = {L_{ijkl}}{\varepsilon _{kl}}$ 。 ${\varepsilon _{ij}}$ 为拉格朗日应变的有限变形表达，具体形式为：

${\varepsilon _{ij}} = \frac{1}{2}\left(\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}} + \frac{{\partial {u_k}}}{{\partial {x_i}}}\frac{{\partial {u_k}}}{{\partial {x_j}}}\right)$

(6)

对式(1)取一阶变分可以得到系统的运动方程，即：

${\text{δ}}\bar E = {\text{δ}}{\smallint} _V{W_\varepsilon }{\text{d}}V = {\smallint }_V{\sigma _{ij}}{\text{δ}}{\varepsilon _{ij}}{\text{d}}V - \int_{\partial V} {{f_i}{\text{δ}}{u_i}{\text{d}}{\varGamma }} + {\smallint} _V\rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}}{\text{δ}}{u_i}{\text{d}}V$

(7)

在稳定系统中，总势能不但满足式(7)，而且它的二阶变分应大于零^[39]，即：

${{\text{δ}}^2}\bar E = [{\bar E_{,uu}}(u({\varLambda })){\text{δ}}u]{\text{δ}}u {\text{＞}} 0$

式中： ${\varLambda }$ 为载荷参数。

当系统达到失稳临界状态时，载荷参数 ${\varLambda }$ 达到临界值 $\mathop {\varLambda }\limits^{\mathrm{o}}$ ，位移场 $u(\mathop {\varLambda }\limits^{\mathrm{o}} )$ 和对应的应力场 $\sigma (\mathop {\varLambda }\limits^{\mathrm{o}} )$ 均达到临界状态，此时：

${{\text{δ}}^2}\bar E = [{\bar E_{,uu}}(u(\mathop {\varLambda }\limits^{\mathrm{o}} )){\text{δ}}u]{\text{δ}}u = 0$

(8)

将位移变分 ${\text{δ}}u$ 记为给定特征位移模式的 $\Delta u$ ，考虑到应力场 $\sigma _{ij}^{(0)} = {L_{ijkl}}{\varepsilon _{kl}}$ (上角标⁽⁰⁾表示远场)，则由系统失稳的临界条件可得到系统的稳定性函数为：

$S({\varLambda }) = \int_V {[{L_{ijkl}}} \Delta {\varepsilon _{kl}}{\text{δ}}{\varepsilon _{ij}} + \sigma _{ij}^{(0)}\Delta {\text{δ}}{\varepsilon _{ij}}]{\text{d}}V - \int_{\partial V} {{f_i}{\text{δ}}\Delta {u_i}{\text{d}}{\varGamma }} - {\int} _V\rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}}{\text{δ}}\Delta {u_i}{\text{d}}V$

(9)

式中：应变变分的表达式为：

${\text{δ}}{\varepsilon _{ij}} = \frac{1}{2}{\text{δ}}\left(\frac{{\partial u_i^0}}{{\partial {x_j}}} + \frac{{\partial u_j^0}}{{\partial {x_i}}} + \frac{{\partial u_k^0}}{{\partial {x_i}}}\frac{{\partial u_k^0}}{{\partial {x_j}}}\right) = {({\text{δ}}{u_{i,j}} + u^0_{k,i}{\text{δ}}{u_{k,j}})_{\text{Q}}}$

(10)

式中：上角标⁰表示临界位移，Q代表局部失稳位置。

系统局部剪切带状失稳分析如图1所示。远场应力 $\sigma _{ij}^{(0)}$ 作用下产生局部剪切变形区，Hill等^[30]采用零应力条件来处理此区的边界，但是，在动力学过程中，在局部剪切变形区的边界显然存在惯性力。因此，本文中的分析将考虑边界惯性力的影响，如图1(c)所示。基于上述分析，剪切失稳系统的运动方程为：

$\sigma _{ij,j}^{(0)} = \rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}}\qquad {i,j = 1,2}$

(11)

根据高斯散度定理，局部运动方程可进一步表达为：

${[{L_{ijkl}}\Delta {u_{k,l}} + \sigma _{pj}^{(0)}\Delta {u_{i,p}}]_{,j}} = {\left( {\rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}}} \right)_{\text{Q}}}$

(12)

下面将基于式(12)和剪切边界条件，即：x₂=0时， $\sigma _{21}^{(0)} = {\tau _0}$ ， $\sigma _{1i}^{(0)} = 0$ (i=1, 2)进行分析。

1.2 动态剪切局部化失稳的扩散分析

由式(12)，表征界面失稳临界状态的运动方程可进一步表达为：

${L_{1212}}\Delta {u_{1,22}} + {L_{1111}}\Delta {u_{1,11}} + \sigma _{21}^{(0)}\Delta {u_{1,21}} + {L_{1122}}\Delta {u_{2,21}} + {L_{1221}}\Delta {u_{2,12}} = \rho \frac{{{\partial ^2}\Delta {u_1}}}{{\partial {t^2}}}$

(13)

${L_{2211}}\Delta {u_{1,12}} + {L_{2112}}\Delta {u_{1,21}} + {L_{2222}}\Delta {u_{2,22}} + {L_{2121}}\Delta {u_{2,11}} + \sigma _{21}^{(0)}\Delta {u_{2,21}} = \rho \frac{{{\partial ^2}\Delta {u_2}}}{{\partial {t^2}}}$

(14)

式中：本构关系 ${\sigma _{ij}} = {L_{ijkl}}{\varepsilon _{kl}}$ 采用Jaumann率下塑性可膨胀本构关系^[40]，则刚度张量 ${L_{ijkl}}$ 表达为:

${L_{ijkl}} = \dfrac{E}{{1 + \mu }}\left[ {\dfrac{1}{2}\left( {{g_{ik}}{g_{jl}} + {g_{il}}{g_{jk}}} \right) + {g_{ij}}{g_{kl}}\left( {\dfrac{{\mu - {E/ {3E_{{\text{tm}}}^{\text{p}}}}}}{{1 - 2\mu + {E / {E_{{\text{tm}}}^{\text{p}}}}}}} \right) - \dfrac{3}{{2\sigma _{\text{e}}^{\text{2}}}}\dfrac{E}{{E_{{\text{te}}}^{\text{p}}}}\dfrac{{{S_{ij}}{S_{kl}}}}{{\dfrac{2}{3}\left( {1 + \mu } \right) + \dfrac{E}{{E_{{\text{te}}}^{\text{p}}}}}}} \right]$

(15)

式中：等效塑性切线模量 $E_{{\text{te}}}^{\text{p}} = {{{\text{d}}{\sigma _{\text{e}}}} /{{\text{d}}\varepsilon _{\text{e}}^{\text{p}}}}$ ，无损情况下， $\dfrac{1}{{E_{{\text{te}}}^{\text{p}}}} = \dfrac{1}{{{E_{\text{t}}}}} - \dfrac{1}{E}$ ；

单轴加载下的切线模量 ${E_{\text{t}}} = {{{\text{d}}{\sigma _1}}/ {{\text{d}}{\varepsilon _1}}}$ ， ${\text{d}}{\sigma _{\text{e}}} = \dfrac{3}{{2{\sigma _{\text{e}}}}}{S_{kl}}{\text{d}}{\sigma _{kl}}$ ；

体积模量 $E_{{\text{tm}}}^{\text{p}} = {{{\text{d}}{\sigma _{kk}}} / {{\text{d}}\varepsilon _{kk}^{\text{p}}}} = {{{\text{d}}{\sigma _m}} / {{\text{d}}\varepsilon _m^{\text{p}}}}$ ；

等效应力 ${\sigma _{\text{e}}} = \sqrt {\dfrac{3}{2}{S_{ij}}{S_{ji}}}$ ；等效塑性应变 $\varepsilon _{\text{e}}^{\text{p}} = \sqrt {\dfrac{3}{2}e_{ij}^{\text{p}}e_{ji}^{\text{p}}}$ 。

基于实验对岩石应力应变全过程^[34-35]进行拟合，塑性阶段： ${E/{E_{{\text{te}}}^{\text{p}}}}$ =0～3， ${E /{E_{{\text{tm}}}^{\text{p}}}}$ =0～120。结合文献[26,34-35]中的实验过程分析和局部化分叉的现象选取材料参数，如表1所示，其中 $\lambda$ 和 $\mu$ 为拉梅系数。取塑性应变 $\varepsilon _{\text{e}}^{\text{p}}$ =0.001，参数 ${E/ {E_{{\text{tm}}}^{\text{p}}}}$ 和 ${E/ {E_{{\text{te}}}^{\text{p}}}}$ 采用文献[40]推荐的方法进行试算，反复迭代确定其内部应力状态而得到。

表 1 岩石材料参数

Table 1. Rock material parameters

ρ/（kg·m⁻³）	E/GPa	ν	λ/GPa	μ/GPa	${E/{E_{{\text{tm}}}^{\text{p}}}}$	${E / {E_{{\text{te}}}^{\text{p}}}}$
2600	30	0.23	10.4	12.2	50	2

下载: 导出CSV

| 显示表格

求解剪切局部化失稳的扩散行为。Hill等^[30]给出了直角坐标系下速率变分的表达式为 ${\text{δ}}{V_i} = {V_i}\left( {{n_0}} \right)$ ， ${n_0} = {n_{0i}}{x_i}$ ， $i = 1,2$ 。其梯度为（图1）： ${\text{δ}}{V_{i,j}} = \dfrac{{{\text{d}}{V_i}}}{{{\text{d}}n}}{n_j} = {c_i}{n_j}$ ， ${n_1} = \cos \alpha$ ， ${n_2} = \sin \alpha$ 。进一步，采用流函数 $\psi = v\left( {{x_2}} \right)\cos \left( {{c_1}{x_1}} \right)$ 来求解速度扰动。李国琛等^[40]对平面应变问题采用的速率扰动函数表达式为： ${\text{δ}}{V_\alpha } = \sin \left(\dfrac{{{\text{π}}{x_\alpha }}}{L}\right)\varphi {{\varDelta }_2}$ ， ${\text{δ}}{V_\alpha } = \sin \left(\dfrac{{{\text{π}}{x_\alpha }}}{L}\right)\varphi {{\varDelta }_4}$ ，其中： ${{\varDelta }_2}$ 、 ${{\varDelta }_4}$ 为待求的广义速率参量， $\varphi = 1 + \left( {1 - \varpi } \right)\cos \dfrac{{2{\text{π}}{x_3}}}{h} - \varpi \cos \dfrac{{4{\text{π}}{x_3}}}{h}$ ， $\varpi$ 为获取最低临界值指标的调节变量。Miles等^[41]、Chau^[33]则采用一种对称形位移函数 $\chi \left( {r, {\textit{z}} } \right) = X\left( r \right)\cos \left( {\eta {\textit{z}} } \right)$ ，其中：η=kπ/L，k=1, 2, ···。对于剪切问题，Chau^[42]采用一种反对称势函数形式： $\varphi = {C_1}{\psi _1}(r)\sin \left( {\eta {\textit{z}} } \right)\cos \left( {n\theta } \right)$ ， $\psi = {C_2}{\psi _2}(r)\cos \left( {\eta {\textit{z}} } \right)\sin \left( {n\theta } \right)$ 。其中 $\eta = {{m{\text{π}}} \mathord{\left/ {\vphantom {{m{\text{π}}} L}} \right. } L}$ ， C₁、C₂为常数, m, n=0, 1, 2, ···。

对于剪切局部化失稳的扩散行为，参照Chau^[42]的反对称势函数形式，可设本征位移为：

$\Delta {u_1} = {U_1}({x_2})\sin (\omega {x_1})\sinh \left( { - {k_1}t} \right)$

(16)

$\Delta {u_2} = {U_2}({x_2})\cos (\omega {x_1})\sinh \left( { - {k_2}t} \right)$

(17)

式中：U_i为关于x₂的位移函数，k₁、k₂为参数。

将式(16)～(17)代入(13)～(14)，为简化计算近似认为两个方向的动力影响相同，可以得到描述局部化失稳的扩散行为的控制方程组为：

${L_{1212}}{U''_1} + \frac{{\sigma _{21}^{\left( 0 \right)}\omega }}{{\tan (\omega {x_1})}}{U'_1} - \left( {{\omega ^2}{L_{1111}} + \rho k_1^2} \right){U_1} - \omega \left( {{L_{1122}} + {L_{1221}}} \right){U'_2} = 0$

(18)

${L_{2222}}{U''_2} - \sigma _{21}^{\left( 0 \right)}\omega \tan (\omega {x_1}){U'_2} - \left( {{\omega ^2}{L_{2121}} + \rho k_2^2} \right){U_2} + \omega \left( {{L_{2211}} + {L_{2112}}} \right){U'_1} = 0$

(19)

进一步整理，可得到：

${a_1}U_1^{(4)} + {b_1}{U'''_1} + {c_1}{U''_1} + {d_1}{U'_1} + {e_1} = 0$

(20)

${a_2}U_2^{(4)} + {b_2}{U'''_2} + {c_2}{U''_2} + {d_2}{U'_2} + {e_2} = 0$

(21)

式中：

${a_1} = \dfrac{{{L_{2222}}{L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}} {\text{，}}\;\;\;\; {a_2} = - \dfrac{{{L_{2222}}{L_{1212}}}}{{\omega \left( {{L_{2211}} + {L_{2112}}} \right)}}$

${b_1} = \frac{{{{\sigma _{21}^{\left( 0 \right)}}/ {\tan (\omega {x_1})}}{L_{2222}}}}{{{L_{1122}} + {L_{1221}}}} - \frac{{\sigma _{21}^{\left( 0 \right)}\tan (\omega {x_1}){L_{1212}}}}{{{L_{1122}} + {L_{1221}}}} {\text{，}}\quad {b_2} = - \frac{{{{\sigma _{21}^{\left( 0 \right)}} / {\tan (\omega {x_1})}}{L_{2222}}}}{{{L_{2211}} + {L_{2112}}}} + \frac{{\sigma _{21}^{\left( 0 \right)}\tan (\omega {x_1}){L_{1212}}}}{{{L_{2211}} + {L_{2112}}}}$

${c_1} = - \frac{{{L_{2222}}\left( {{\omega ^2}{L_{1111}} + \rho {k^2_1}} \right)}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}} - \frac{{{{\left( {\sigma _{21}^{\left( 0 \right)}} \right)}^2}\omega }}{{{L_{1122}} + {L_{1221}}}} - \frac{{\left( {{\omega ^2}{L_{2121}} + \rho {k^2_2}} \right){L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}} + \omega \left( {{L_{2211}} + {L_{2112}}} \right)$

${c_2} = \frac{{{L_{2222}}\left( {{\omega ^2}{L_{1111}} + \rho {k^2_1}} \right)}}{{\omega \left( {{L_{2211}} + {L_{2112}}} \right)}} + \frac{{{{\left( {\sigma _{21}^{\left( 0 \right)}} \right)}^2}\omega }}{{{L_{2211}} + {L_{2112}}}} + \frac{{\left( {{\omega ^2}{L_{2121}} + \rho {k^2_2}} \right){L_{1212}}}}{{\omega \left( {{L_{2211}} + {L_{2112}}} \right)}} - \omega \left( {{L_{1122}} + {L_{1221}}} \right)$

${d_1} = \frac{{\sigma _{21}^{\left( 0 \right)}\tan (\omega {x_1})\left( {{\omega ^2}{L_{1111}} + \rho {k^2_1}} \right)}}{{{L_{1122}} + {L_{1221}}}} - \frac{{{{\sigma _{21}^{\left( 0 \right)}} / {\tan (\omega {x_1})}}\left( {{\omega ^2}{L_{2121}} + \rho {k^2_2}} \right)}}{{{L_{1122}} + {L_{1221}}}}$

${d_2} = - \frac{{\sigma _{21}^{\left( 0 \right)}\tan (\omega {x_1})\left( {{\omega ^2}{L_{1111}} + \rho {k_1^2}} \right)}}{{{L_{2211}} + {L_{2112}}}} + \frac{{{{\sigma _{21}^{\left( 0 \right)}} / {\tan (\omega {x_1})}}\left( {{\omega ^2}{L_{2121}} + \rho {k_2^2}} \right)}}{{{L_{2211}} + {L_{2112}}}}$

${e_1} = \frac{{\left( {{\omega ^2}{L_{1111}} + \rho {k^2_1}} \right)\left( {{\omega ^2}{L_{2121}} + \rho {k^2_2}} \right)}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}} {\text{，}}\quad{e_2} = - \frac{{\left( {{\omega ^2}{L_{1111}} + \rho {k^2_1}} \right)\left( {{\omega ^2}{L_{2121}} + \rho {k^2_2}} \right)}}{{\omega \left( {{L_{2211}} + {L_{2112}}} \right)}}$

此控制方程组即为动力学过程中局部失稳存在的判断方程。若 $\tan( \omega {x_1})$ 存在非平凡解，则系统中存在剪切局部化变形带，剪切带内的变形分布将由 $\tan (\omega {x_1})$ 决定。对其直接求解比较困难，下面将进行简化分析。

2. 岩石动态剪切失稳求解

2.1 基于剪切变形带的求解

岩石动摩擦滑移主要发生在界面两侧1～5 mm范围内^[9]，剪切变形带作为界面破坏滑移的先导过程，一般小于此范围，岩石界面失稳主要发生在非常狭窄的区域。因此，直角坐标系下速率^[30,40]变分的平面应力分量可采用Hill等^[30]的表达式，即为：

${\text{δ}}{V_{i,j}} = {V_i}\left( n \right)\qquad {n = {n_i}{x_i}}{\text{，}} i=1,2$

(22)

式中：V_i为关于x_i的速度函数，n_i为参数。

其梯度为：

${\text{δ}}{V_{i,j}} = \frac{{{\text{d}}{V_i}}}{{{\text{d}}n}}{n_j} = {c_i}{n_j}{\text{，}}\qquad {n_1} = \cos \alpha {\text{，}}\qquad {{n_2} = \sin \alpha }$

(23)

类似地，对于剪切局部化失稳的扩散行为，引入动力学项影响，可设本征位移为：

$\Delta {u_i} = {X_i}\left( {n,k} \right)\qquad {n = {n_i}x}{\text{，}} {k = kt} {\text{，}} i=1,2$

(24)

式中：X_i为关于x_i的位移函数，k为局部动力参数。

将式(24)代入式(12)，可得到描述剪切局部化失稳的扩散行为的控制方程组：

${L_{1212}}{X''_1}n_2^2 + {L_{1111}}{X''_1}n_1^2 + \sigma _{21}^{\left( 0 \right)}{X''_1}{n_1}{n_2} + \left( {{L_{1122}} + {L_{1221}}} \right){X''_2}{n_1}{n_2} = \rho {X''_1}{k^2}$

(25)

$\left( {{L_{2211}} + {L_{2112}}} \right){X''_1}{n_1}{n_2} + {L_{2222}}{X''_2}n_2^2 + {L_{2121}}{X''_2}n_1^2 + \sigma _{21}^{\left( 0 \right)}{X''_2}{n_1}{n_2} = \rho {X''_2}{k^2}$

(26)

设 $r = \dfrac{{{{X_1''}}}}{{{{X_2''}}}}$ ， $s = \dfrac{{{n_2}}}{{{n_1}}}$ ， ${\textit{z}} = \dfrac{{{k^2}}}{{{n_1}{n_2}}} = {k^2}\left( {s + \dfrac{1}{s}} \right)$ 代入式(25）～（26)，可得：

$\begin{split} {L_{1212}}{L_{2222}}{s^4} +& \left( {\sigma _{21}^{\left( 0 \right)} - \rho {\textit{z}} } \right)\left( {{L_{1212}} + {L_{2222}}} \right){s^3} + \left( {L_{1111}}{L_{2222}} + {L_{2121}}{L_{1212}} + {{\left( {\sigma _{21}^{\left( 0 \right)} - \rho {\textit{z}} } \right)}^2} - \right.\\ & \left.\left( {{L_{1122}} + {L_{1221}}} \right)\left( {{L_{2211}} + {L_{2112}}} \right) \right){s^2} + \left( {\sigma _{21}^{\left( 0 \right)} - \rho {\textit{z}} } \right)\left( {{L_{1111}} + {L_{2121}}} \right)s + {L_{1111}}{L_{2121}} = 0 \end{split}$

(27)

此方程即为冲击过程中剪切局部化失稳存在的判断方程。其中，局部动力参数k体现了时间效应的影响。此方程中，若 $s = \tan \alpha$ （α为剪切带角度）存在非平凡解，则系统中存在剪切带。

为研究冲击过程惯性力和载荷对局部化剪切变形带的影响特征，对式(27)无量纲化，取 ${\sigma ^ * } = {{\sigma _{21}^{\left( 0 \right)}}/ G}$ ， ${a^ * } = {{\rho {\textit{z}} } /{\sigma _{21}^{\left( 0 \right)}}}$ ，其中G为剪切模量。则式(27)可以改写为：

$\begin{split} \left( {{{{L_{1212}}{L_{2222}}} / {{G^2}}}} \right){s^4} + &{\sigma ^*}\left( {1 - {a^*}} \right)\left[ {{{\left( {{L_{1212}} + {L_{2222}}} \right)} / G}} \right]{s^3} + \left[ {{{L_{1111}}{L_{2222}}} / {{G^2}}} + {{{L_{2121}}{L_{1212}}} / {{G^2}}} + {{\left( {{\sigma ^*}\left( {1 - {a^*}} \right)} \right)}^2} -\right.\\ &\left.{{\left( {{L_{1122}} + {L_{1221}}} \right)\left( {{L_{2211}} + {L_{2112}}} \right)} / {{G^2}}} \right]{s^2} + {\sigma ^*}\left( {1 - {a^*}} \right)\left[ {{{\left( {{L_{1111}} + {L_{2121}}} \right)} / G}} \right]s + {{{L_{1111}}{L_{2121}}} / {{G^2}}} = 0 \end{split}$

(28)

式(28)即为剪切局部化失稳存在的判别式。此式存在两组对称解（也称为共轭解），取对称部分换算为角度进行分析。图2所示为剪切带角度 $\alpha$ 随无量纲剪切力σ^*和局部动力参数k的变化特征。剪切变形带角度 $\alpha$ 为5°～9°，与岩石直剪实验结果^[35]相近。由此可见：局部动力系数k的影响比外部剪切作用力的影响更为显著；随着外部剪切作用力的增大，剪切变形带角度 $\alpha$ 有一定程度的增大，但增幅不是很大；随着局部动力系数k的增大，即局部惯性作用力的增大，剪切变形带角度 $\alpha$ 明显减小。较强的动态作用过程，惯性效应明显，局部化程度加大，此时的剪切失效也集中于局部区域。

图 2 剪切带角度的动力学影响

Figure 2. Dynamic influence of shear band angle

下载: 全尺寸图片幻灯片

2.2 控制方程的近似求解

对含惯性力作用的式(18)～(21)进行 $\tan( \omega {x_1})$ 的直接求解比较困难，因此，可寻求其近似解，即：先对式(20)～(21)求解其静力学解，然后在静力学解的基础上进行动力学修正。基于此，对式(13）～（14)，令加速度项等于零，即：

${L_{1212}}\Delta {u_{1,22}} + {L_{1111}}\Delta {u_{1,11}} + \sigma _{21}^{\left( 0 \right)}\Delta {u_{1,21}} + \left( {{L_{1122}} + {L_{1221}}} \right)\Delta {u_{2,21}} = 0$

(29)

$\left( {{L_{2211}} + {L_{2112}}} \right)\Delta {u_{1,21}} + {L_{2222}}\Delta {u_{2,22}} + {L_{2121}}\Delta {u_{2,11}} + \sigma _{21}^{\left( 0 \right)}\Delta {u_{2,21}} = 0$

(30)

将式(24)代入式(29）～（30)，可得：

${L_{1212}}{X''_1}n_2^2 + {L_{1111}}{X''_1}n_1^2 + \sigma _{21}^{\left( 0 \right)}{X''_1}{n_1}{n_2} + \left( {{L_{1122}} + {L_{1221}}} \right){X''_2}{n_1}{n_2} = 0$

(31)

$\left( {{L_{2211}} + {L_{2112}}} \right){X''_1}{n_1}{n_2} + {L_{2222}}{X''_2}n_2^2 + {L_{2121}}{X''_2}n_1^2 + \sigma _{21}^{\left( 0 \right)}{X''_2}{n_1}{n_2} = 0$

(32)

设 $r = \dfrac{{{X''_1}}}{{{X''_2}}}$ ， $s = \dfrac{{{n_2}}}{{{n_1}}}$ ，代入(31）～（32)，可以得到：

$\begin{split} {L_{1212}}{L_{2222}}{s^4} + &\sigma _{21}^{\left( 0 \right)}\left( {{L_{1212}} + {L_{2222}}} \right){s^3} + \left[ {L_{1111}}{L_{2222}} + {L_{2121}}{L_{1212}} + {{\left( {\sigma _{21}^{\left( 0 \right)}} \right)}^2} - \right.\\ &\left.\left( {{L_{1122}} + {L_{1221}}} \right)\left( {{L_{2211}} + {L_{2112}}} \right) \right]{s^2} + \sigma _{21}^{\left( 0 \right)}\left( {{L_{1111}} + {L_{2121}}} \right)s + {L_{1111}}{L_{2121}} = 0 \end{split}$

(33)

同理，取 ${\sigma ^ * } = {{\sigma _{21}^{\left( 0 \right)}} \mathord{\left/ {\vphantom {{\sigma _{21}^{\left( 0 \right)}} G}} \right. } G}$ 对方程无量纲化，可得到：

$\begin{split} \left( {{{{L_{1212}}{L_{2222}}} / {{G^2}}}} \right){s^4} + &{\sigma ^*}\left[ {{{\left( {{L_{1212}} + {L_{2222}}} \right)} / G}} \right]{s^3} + \left[ {{{L_{1111}}{L_{2222}}}/ {{G^2}}} + {{{L_{2121}}{L_{1212}}} / {{G^2}}} + {{\left( {{\sigma ^*}} \right)}^2} - \right.\\ &\left. {{\left( {{L_{1122}} + {L_{1221}}} \right)\left( {{L_{2211}} + {L_{2112}}} \right)} / {{G^2}}} \right]{s^2} + {\sigma ^*}\left[ {{{\left( {{L_{1111}} + {L_{2121}}} \right)} / G}} \right]s + {{{L_{1111}}{L_{2121}}}/ {{G^2}}} = 0 \end{split}$

(34)

上式为关于s的多项式，一般有4个根(包括实根和虚根)。一旦得到实根s₀，则表明确实存在剪切带失稳，此时得到静力学的解为： ${\left( {\tan \alpha } \right)_{{\text{static}}}}$ =s₀。然后，令s=s₀+s₁，将其回代到式(28)，并忽略高阶小量，可得到：

${{s_1} = \frac{{{\sigma ^*}{a^*}\left[ {{{\left( {{L_{1212}} + {L_{2222}}} \right)} / G}} \right]s_0^3 - {{\left( {{\sigma ^*}} \right)}^2}\left[ {{{\left( {{a^*}} \right)}^2} - 2\left( {{a^*}} \right)} \right]s_0^2 + {\sigma ^*}{a^*}\left[ {{{\left( {{L_{1111}} + {L_{2121}}} \right)} / G}} \right]{s_0}}}{{4\left( {{{{L_{1212}}{L_{2222}}} / {{G^2}}}} \right)s_0^3 + 3{\sigma ^*}\left[ {{{\left( {{L_{1212}} + {L_{2222}}} \right)} / G}} \right]\left( {1 - {a^*}} \right)s_0^2 + 2\left[ {{{{L_{1111}}{L_{2222}}} / {{G^2}}} + {{{L_{2121}}{L_{1212}}}/ {{G^2}}} - {{\left( {{L_{1122}} + {L_{1221}}} \right)\left( {{L_{2211}} + {L_{2112}}} \right)} / {{G^2}}} + {{\left( {{\sigma ^*}} \right)}^2}{{\left( {1 - {a^*}} \right)}^2}} \right]{s_0} + {\sigma ^*}\left( {1 - {a^*}} \right)\left[ {{{\left( {{L_{1111}} + {L_{2121}}} \right)} / G}} \right]}} }$

(35)

此即为一阶近似解。再次令 $s = {s_0} + {s_1} + {s_2}$ ，将其回代到式(28)，并忽略高阶小量，可得到 ${s_2}$ 的表达式，即为二阶近似解。以此类推，反复迭代，可以得到s在动力条件下的n阶近似表达式，即： $s = {s_0} + \displaystyle\sum\limits_{i = 1}^n {{s_i}}$ 。其二阶近似解计算过程如图3所示。由此可见，惯性效应的存在使得剪切带的角度减小，反映出冲击过程变形高度局部化的趋势。

图 3 冲击作用下剪切局部变形带角度的二阶近似求解

Figure 3. Second-order approximate solution of shear local deformation zone angle under impact

下载: 全尺寸图片幻灯片

3. 岩石动态剪切失稳变形的扩散分析

3.1 失稳变形的扩散分析

若方程组存在剪切局部化变形解，对于剪切变形区域，基于式(20）～（21)的扩散型位移解中的U₁和U₂具有如下形式：

${U_1} = A{}_1{{\text{e}}^{ - {p_1}{x_2}}} + A{}_2{{\text{e}}^{ - {p_2}{x_2}}} + A{}_3{{\text{e}}^{ - {p_3}{x_2}}} + A{}_4{{\text{e}}^{{p_4}{x_2}}}$

(36)

${U_2} = {B_1}{{\text{e}}^{ - {p_1}{x_2}}} + {B_2}{{\text{e}}^{ - {p_2}{x_2}}} + {B_3}{{\text{e}}^{ - {p_3}{x_2}}} + {B_4}{{\text{e}}^{{p_4}{x_2}}}$

(37)

式中：A₁、A₂、A₃、A₄为待定系数；

${p_1} = {Z_3} + \frac{{{b_1}}}{{4{a_1}}} + {Z_2} \text{，} {p_2} = {Z_3} + \frac{{{b_1}}}{{4{a_1}}} - {Z_2} \text{，} {p_3} = \frac{{{b_1}}}{{4{a_1}}} - {Z_3} + {Z_1} \text{，} {p_4} = {Z_3} - \frac{{{b_1}}}{{4{a_1}}} + {Z_1} \text{；}$

${Z_1} = \frac{{\sqrt { - 9{Z^{2/3}_6}\sqrt {{Z_5}} - 12{Z_8}\sqrt {{Z_5}} - {Z^2_9}\sqrt {{Z_5}} - 3\sqrt 6 {Z_{10}}\sqrt {27{Z^2_{10}} - 72{Z_8}{Z_9} + 2{Z^3_9} + 3\sqrt 3 {Z_7}} - 12{Z_9}{Z^{1/3}_6}\sqrt {{Z_5}} } }}{{{Z_4}}} \text{，}$

${Z_2} = \frac{{\sqrt {3\sqrt 6 {Z_{10}}\sqrt {27{Z^2_{10}} - 72{Z_8}{Z_9} + 2{Z^3_9} + 3\sqrt 3 {Z_7}} - 12{Z_8}\sqrt {{Z_5}} - {Z^2_9}\sqrt {{Z_5}} - 9{Z^{2/3}_6}\sqrt {{Z_5}} - 12{Z_9}{Z^{1/3}_6}\sqrt {{Z_5}} } }}{{{Z_4}}} \text{，}$

${Z_3} = \frac{{\sqrt {{Z_5}} }}{{6{Z^{1/6}_6}}} \text{，} {Z_4} = 6{Z^{1/6}_6}{Z^{1/4}_5} \text{，} {Z_5} = {Z^2_9} + 9{Z^{2/3}_6} - 6{Z_9}{Z^{1/3}_6} + \frac{{12{e_1}}}{{{a_1}}} - \frac{{9b_1^4}}{{64a_1^4}} + \frac{{3b_1^2{c_1}}}{{4a_1^3}} - \frac{{3{b_1}{d_1}}}{{a_1^2}} \text{，}$

${Z_6} = \frac{{{Z^2_{10}}}}{2} - \frac{{4{Z_8}{Z_9}}}{3} + \frac{{{Z^3_9}}}{{27}} + \frac{{\sqrt 3 {Z_7}}}{{18}} \text{，} {Z_7} = \sqrt {27{Z^4_{10}} - 256{Z^3_8} - 16{Z_8}{Z^4_9} + 4{Z^3_9}{Z^2_{10}} + 128{Z^2_8}{Z^2_9} - 144{Z_8}{Z_9}{Z^2_{10}}} \text{，}$

${Z_8} = \frac{{{e_1}}}{{{a_1}}} - \frac{{3b_1^4}}{{256a_1^4}} + \frac{{b_1^2{c_1}}}{{16a_1^3}} - \frac{{{b_1}{d_1}}}{{4a_1^2}}\text{，}{Z_9} = \frac{{{c_1}}}{{{a_1}}} - \frac{{3b_1^2}}{{8a_1^2}}\text{，} {Z_{10}} = \frac{{{d_1}}}{{{a_1}}} + \frac{{b_1^3}}{{8a_1^3}} - \frac{{{b_1}{c_1}}}{{2a_1^2}} \text{；}$

$\begin{split}& {B_1} = \left[ {\frac{{{L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}}{p_1} + \frac{{{{\sigma _{21}^{\left( 0 \right)}} / s}}}{{{L_{1122}} + {L_{1221}}}} - \frac{{\omega \left( {{L_{1111}} + \rho {k^2_1}} \right)}}{{{L_{1122}} + {L_{1221}}}}\frac{1}{{{p_1}}}} \right]{A_1} \text{，} \\ & {B_2} = \left[ {\frac{{{L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}}{p_2} + \frac{{{{\sigma _{21}^{\left( 0 \right)}} / s}}}{{{L_{1122}} + {L_{1221}}}} - \frac{{\omega \left( {{L_{1111}} + \rho {k^2_1}} \right)}}{{{L_{1122}} + {L_{1221}}}}\frac{1}{{{p_2}}}} \right]{A_2} \text{，} \end{split}$

$\begin{split}&{B_3} = \left[ {\frac{{{L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}}{p_3} + \frac{{{{\sigma _{21}^{\left( 0 \right)}} / s}}}{{{L_{1122}} + {L_{1221}}}} - \frac{{\omega \left( {{L_{1111}} + \rho {k^2_1}} \right)}}{{{L_{1122}} + {L_{1221}}}}\frac{1}{{{p_3}}}} \right]{A_3} \text{，}\\ &{B_4} = \left[ {\frac{{{L_{1212}}}}{{\omega \left( {{L_{1122}} + {L_{1221}}} \right)}}{p_4} + \frac{{{{\sigma _{21}^{\left( 0 \right)}} / s}}}{{{L_{1122}} + {L_{1221}}}} - \frac{{\omega \left( {{L_{1111}} + \rho {k^2_1}} \right)}}{{{L_{1122}} + {L_{1221}}}}\frac{1}{{{p_4}}}} \right]{A_4} \end{split}$

剪切作用下，即：x₂=0时， $\sigma _{21}^{\left( 0 \right)} = {\tau _0}$ ， $\sigma _{1i}^{\left( 0 \right)} = 0$ (i=1, 2)。此问题的解以剪切变形为主，因此，为简化问题起见，在剪切变形带的边界，忽略法向惯性力，即令 $\rho \dfrac{{{\partial ^2}\Delta {u_1}}}{{\partial {t^2}}}$ =0。同时，考虑其切向惯性力，取 $\rho \dfrac{{{\partial ^2}\Delta {u_2}}}{{\partial {t^2}}}$ = $\tau$ 。则有：

${L_{2112}}\Delta {u_{1,2}} + {L_{2121}}\Delta {u_{2,1}} + \sigma _{21}^{\left( 0 \right)}\Delta {u_{2,2}} = \tau$

(38)

${L_{1111}}\Delta {u_{1,1}} + {L_{1122}}\Delta {u_{2,2}} + \sigma _{21}^{\left( 0 \right)}\Delta {u_{1,2}} = 0$

(39)

此过程中， ${\tau / {\sigma _{21}^{\left( 0 \right)}}} = \beta$ 。其中，β为无量纲的切向惯性力。为简化讨论过程，β取0.1。

由于岩石动摩擦滑移主要发生在界面两侧1～5 mm范围内，此过程为破坏滑移的先导过程，失稳界面变形的位移小于滑移位移。剪切变形带通常处于一个狭窄区域，靠近基体侧的位移更小。因此，可增加位移边界条件：

${x_2} = h,\qquad {\Delta {u_i} = 0}\qquad i = 1,2$

(40)

式中：h为岩石动摩擦影响厚度。式(36)的待定系数A₁、A₂、A₃、A₄可由上述边界条件确定。

另外，参数k₁、k₂由下式所列的初始条件确定：

$t = 0 ,\qquad {\Delta {u_i} = 0},\qquad {\frac{{\partial \Delta {u_i}}}{{\partial t}} = 0} \qquad i = 1,2$

(41)

定义式(33)表达的等效位移u表征岩石界面失稳的位移。

$u = \sqrt {{{\left( {\Delta {u_1}} \right)}^2} + {{\left( {\Delta {u_2}} \right)}^2}}$

(42)

基于此，由式(36）～（42)，结合表1所列岩石参数，可以进行岩石界面的动态剪切失稳变形的扩散分析。图4所示为局部化失稳变形扩散过程的初步计算结果，失稳界面的厚度h取为0.3 mm，此时， $\tan \alpha$ =0.1659， ${\sigma ^*}$ =0.002，β=0.1。由此可见：随着深度的增加，剪切扩散作用下的失稳位移越来越小；随着加载时间的增加，局部位移越来越大，直到局部产生破坏形成滑移系。计算结果与实验观察现象^[10-11]基本一致。

图 4 局部化失稳变形扩散过程

Figure 4. Local instability deformation diffusion process

下载: 全尺寸图片幻灯片

3.2 剪切力的影响

图5(a)所示为失稳表面的位移特征随着剪应力增大而变化的规律。在所加载的范围内，随着剪切应力的增大，剪切变形带的角度变化不大，与剪切失稳对应的5个峰的峰位比较接近；但是最大局部位移随剪切应力的增大增加较快。图5(b)所示为将失稳界面的最大位移进行投影后的图，图中剪切力较大时用蓝色和淡蓝色线表示。由此可见：与剪切失稳对应的5个峰的峰位比较接近；界面失稳位移增加幅度随着界面上剪应力的增大而增大。这也表明：当岩石所受剪切作用力过大时，材料会在局部产生失稳的基础上逐渐发展到材料破坏，例如产生岩石类材料常见的局部的剪切破裂。

图 5 界面失稳形貌随剪应力的变化规律

Figure 5. Variation of interface instability morphology with shear stress

下载: 全尺寸图片幻灯片

3.3 界面厚度的影响

图6(a)为失稳表面位移随界面厚度变化的情况，界面厚度分别设置为0.3、0.5和1.0 mm，此时，令 ${\sigma ^*}$ =0.002。随着界面厚度的增大，可变性位移增大，为了研究界面厚度对变形的影响，应用下式定义的相对变形来分析：

图 6 随着厚度变化的界面失稳形貌

Figure 6. Interface instability morphology with thickness variation

下载: 全尺寸图片幻灯片

${\varepsilon _{\text{h}}} = \frac{{{u_{{\text{max}}}}}}{h}$

(43)

等效位移和应变随界面厚度变化的关系如图6(b)所示，随着界面厚度的增大，可变形区域增大；由于变形区的累积和几何改变，界面表面褶皱加深，即界面越厚，失稳位移越大。

4. 岩石界面剪切失稳数值分析

4.1 有限元模型和材料参数

应用ABAQUS软件建立如图7所示的双层结构模型进行有限元模拟计算。模型尺寸为12.5 mm×40 mm×100 mm。模型中岩石失稳界面等效成一个厚度较小的弹性层，其上层材料与下层基底完全粘结，右端施加完全约束条件，左端施加压应力，底部与侧面限制离面位移。动态载荷在ABAQUS软件载荷界面使用表（tabular）添加，如图7所示，上表面局部区域施加恒值剪力载荷作用，以模拟岩石受压过程中的局部剪切作用力。有限元计算网格类型采用C3D8R，同时对受局部剪力部分的网格进行加密。由此数值分析局部剪应力处界面失稳特征。

图 7 有限元计算模型

Figure 7. FEM model

下载: 全尺寸图片幻灯片

基体和界面层的材料参数选取主要基于张磊等^[10-11]对含斜节理岩石试样进行SHPB杆束冲击的实验结果。图8(a)所示为冲击前后岩石表面粗糙系数的变化规律，图8(b)所示为岩石界面摩擦应力与表面粗糙系数对比度的关系。冲击后岩石表面的粗糙系数有一定程度的升高，冲击后粗糙表面的高程有0.2～0.4 mm的降低，均反映出一定程度的局部破碎的加剧。由表面的粗糙系数的变化幅度估计摩擦界面区材料性能的折减系数约0.2～0.4，为讨论界面影响，有限元计算参数选取如下：基体密度为2600 kg/m³，弹性模量为30 GPa，泊松比为0.23；界面区按照基体材料属性进行折减，折减系数为0.1、0.2、0.3、0.4和0.5。

图 8 不同角度下界面的动摩擦性能

Figure 8. Interface dynamic friction properties at different angles

下载: 全尺寸图片幻灯片

4.2 有限元计算结果

图9所示为局部剪切力为20 MPa、折减系数为0.1时的位移形貌。计算结果表明：当试样受到外部载荷时，若局部薄弱处存在集中剪力，则会产生褶皱行为；随着剪应力的增大，局部变形增加；试样将会沿着失稳界面处产生滑移等大变形。利用数值模拟得到了类似现象，印证了上述理论的合理性。在有限变形的条件下，随着加载进一步施加，即会产生大变形，进而产生破坏。局部应力过大导致岩石材料分叉，为材料破坏的先导阶段。

图 9 有限元计算结果

Figure 9. Result of FEM

下载: 全尺寸图片幻灯片

4.3 界面参数的影响

对于小尺度的界面薄层，确定其具体的力学参数是十分困难的，类似于混凝土界面过渡区（interfacial transition zone，ITZ），大部分是根据模量进行折减，并经过有限元反复试算得到^[43]。近年来，随着纳米压痕技术的发展，对中微尺度薄弱界面的力学性能参数有一定程度的评估，但获得其准确参数仍然是个挑战。下面将基于杆束实验^[10-11]前后粗糙度的变化程度，讨论参数对界面起皱形貌的影响，为详细讨论，按照折减系数0.1、0.2、0.3、0.4和0.5给予界面层参数。图10所示为最大位移（u_max）与折减系数的关系，此时保持界面厚度为0.3 mm, 局部剪切力为20 MPa。

图 10 界面参数的影响

Figure 10. Effect of interface parameters

下载: 全尺寸图片幻灯片

从图10可以看出，界面参数对数值计算结果影响较大，因此，在数值计算中应选择合适的计算参数来讨论其影响。当界面参数与基体参数差异较小时不易引起起皱，差异较大时，在局部受载时容易形成失稳扰动，影响表面位移形貌。为了方便讨论此现象，以下分析中的折减系数取为0.1。

4.4 剪切力的影响

图11所示为界面厚度为0.3 mm，局部剪切力分别为10、20、30 MPa条件下的有限元计算结果。图11结果表明，当试样受到外部载荷时，试样局部薄弱位置产生分叉变形，表现为起皱现象；随着剪应力的增大，形貌变化，局部变形加剧；随着加载时间的延长，局部变形逐渐由小变形向大变形转变，导致试样将沿着失稳截面处产生滑移，从而产生损伤破坏。

图 11 局部剪切力的影响

Figure 11. Effect of local shear force

下载: 全尺寸图片幻灯片

4.5 界面厚度的影响

图12所示为剪切力20 MPa，界面厚度分别为0.3、0.5、1 mm条件下的有限元计算结果。由图12可以看出，随着失稳界面的增厚，局部力所产生的响应将被失稳界面层的变形进一步消耗掉，峰值整体位移表现为增大。

图 12 界面厚度对位移和应变的影响

Figure 12. Influences of interface thickness on displacement and strain

下载: 全尺寸图片幻灯片

岩石SHPB杆束实验^[10-11]中，由于岩石的不均匀性，发生分叉时各处界面厚度变化不一，因此裂纹萌生发展并不仅仅产生在一处，而是会形成多处裂纹路径，致密性较好的岩石可能呈规律的片状破坏（沿着分叉处滑移破坏），非均匀程度较大的岩石或呈无规律性的剪切破坏（分叉处较多整体表现较为无序）^[44-45]。

4.6 局部化变形对波传播的影响

图13上半部分所示为剪切力为20 MPa，界面厚度分别为0.3、0.5、1.0 mm条件下，界面中心位置上层和下层应力变化的计算结果。结果表明，局部剪切力的大小随着界面层的起皱，载荷输入能量会被消耗，体现出衰减现象；弱界面越厚，相对衰减幅度越大。图13下部分所示为对失稳界面上层和下层两条波进行傅里叶变换得到的频域图，发现相关频域的幅值发生了改变，即界面前后不仅波的幅值改变了，波的频率也改变了。

图 13 界面层上下的应力传递

Figure 13. Stress transfer above and below the interface layer

下载: 全尺寸图片幻灯片

岩石SHPB杆束实验^[10-11]中，由于岩石的不均匀性，发生分叉时各处界面厚度变化不一，这种变形会消耗加载波的能量使之产生衰减，若局部化失稳持续产生，则会演化成为断裂行为；若仅局部化到一定程度继续加载，宏观应力-应变曲线上表现为：在塑性阶段，随着应变的增大，应力先微小下降再上升。此现象为断裂破坏的先导过程，表现出剪切变形带耗能的特点。局部化剪切变形带的形成和发展，为断裂滑移的先导，研究此现象及其动力学特性，可为评估局部化失稳对岩石局部化断裂、地震波传播过程中滑移界面等的影响提供参考。

5. 结　论

利用有限变形理论，对岩石界面进行了失稳分析，建立了岩石材料界面剪切局部化失稳的理论分析模型；通过实验和数值分析，说明了模型的可靠性；基于此分析了岩石界面失稳的影响因素，发现变形局部化对波传播具有衰减作用；得到的主要结论如下。

(1)基于能量法和扰动分析，提出了界面受剪力影响下局部化失稳的判别式和局部变形扩散方程。

(2)较强的动态作用过程中，惯性效应的影响明显，岩石界面局部化程度较深。剪切变形带角度随着外部剪切作用力的增大而增大，随着局部动力系数的增大明显减小。

(3)剪切力越大，局部位移越大，岩石越容易发生破坏；界面厚度越小，局部位移越大，更容易沿界面破坏。界面剪切扩散行为极大降低了透射波的幅值，同时也改变了透射波的频率。

图 1 爆炸实验水池示意图

Figure 1. Sketch of explosion testing pool

下载: 全尺寸图片幻灯片

图 2 爆炸水池模型及其网格划分情况

Figure 2. Sketch of finite element structure and meshing condition

下载: 全尺寸图片幻灯片

图 3 数值计算拟合曲线与P.Cole压力曲线的对比

Figure 3. Contrast curves of numerical calculation and P. Cole empirical formulas

下载: 全尺寸图片幻灯片

图 4 模型Ⅰ水中冲击波传播过程

Figure 4. Propagation process of chock wave in water (model Ⅰ)

下载: 全尺寸图片幻灯片

图 5 观察点A、C两点比冲量时程曲线

Figure 5. Contrast curves of specific impulse in positions A, C

下载: 全尺寸图片幻灯片

图 6 模型Ⅱ与模型Ⅰ冲击波压力对比曲线

Figure 6. Contrast curves of shock wave pressure (models Ⅰ & Ⅱ)

下载: 全尺寸图片幻灯片

图 7 模型Ⅲ与模型Ⅰ冲击波压力对比曲线

Figure 7. Contrast curves of shock wave pressure (models Ⅰ & Ⅲ)

下载: 全尺寸图片幻灯片

图 8 不同模型位置A比冲量时程曲线

Figure 8. Contrast curves of specific impulse of steel in position A

下载: 全尺寸图片幻灯片

图 9 不同模型位置A等效应力时程曲线

Figure 9. Contrast curves of effective stress of steel in position A

下载: 全尺寸图片幻灯片

图 10 不同模型混凝土围堰的塑性应变区域

Figure 10. Plastic strain zone of concrete cofferdam in each model

下载: 全尺寸图片幻灯片

图 11 不同模型位置B的等效应力时程曲线

Figure 11. Contrast curves of effective stress in position B

下载: 全尺寸图片幻灯片

图 12 不同模型C位置的等效应力时程曲线

Figure 12. Contrast curves of effective stress in position C

下载: 全尺寸图片幻灯片

图 13 不同模型钢板内衬的内能变化曲线

Figure 13. Contrast curves of internal energy in steel lining

下载: 全尺寸图片幻灯片

图 14 不同模型混凝土的内能变化曲线

Figure 14. Contrast curves of internal energy in concrete cofferdam

下载: 全尺寸图片幻灯片

表 1 TNT炸药材料参数

Table 1. Material parameters of TNT explosive

ρ/(g·cm^-3)	D₀/(km·s^-1)	P_CJ/GPa	A/GPa	B/GPa	R₁	R2	ω	E₀/MJ
1.63	6.93	27	371	7.43	4.15	0.95	0.3	0.007

下载: 导出CSV

表 2 内衬钢板材料参数

Table 2. Material parameters of steel inner-lined plate

ρ/(g·cm^-3)	E/GPa	λ	σ_s/GPa	E_t/GPa	C	P	β
7 850	206	0.3	0.337	9.4	40.4	5	1

下载: 导出CSV

表 3 混凝土材料参数

Table 3. Material parameters of concrete

ρ/(g·cm^-3)	E/GPa	λ	σ_s/GPa	E_t/GPa	C	P	β
2 400	30	0.18	0.035	12.7	0	0	0.8

下载: 导出CSV

表 4 土壤压力与体应变关系

Table 4. Relation between pressure and volumetric strain

体应变	0.0	0.05	0.09	0.15	0.19	0.21	0.22	0.25
压力/MPa	0.0	0.34	0.45	1.27	2.08	2.71	3.92	5.66

下载: 导出CSV

表 5 水的状态参数

Table 5. State parameters of water

ρ/(g·cm^-3)	c/(m·s^-1)	S₁	S₂	S₃	γ₀	α	E₀/MJ
1 020	1 650	1.92	-1.98	0.22	0.5	0	0

下载: 导出CSV

表 6 水中冲击波峰值压力计算值与理论值对比

Table 6. Contrast of numerical calculation and empirical formulas

距离/m	峰值压力/MPa		计算与理论值的相对误差/%
距离/m	P. Cole理论计算值	数值计算结果平均值	计算与理论值的相对误差/%
2.0	57.97	64.1	10.6
2.5	45.05	47.0	4.3
3.0	36.66	37.4	2.0
3.5	30.80	29.8	-3.2
4.0	26.49	25.3	-4.5
4.5	23.19	21.6	-6.8
5.0	20.59	31.2^*
注：由于池壁位置入射压力峰值被反射压力峰值覆盖，因此此值仅为反射压力峰值。

下载: 导出CSV

表 7 空气桶和气泡帷幕对钢板内衬典型位置应力峰值的影响

Table 7. Effect of air tank and bubble curtain on effective stress in typical positions of steel lining

典型位置	应力峰值/MPa				衰减率/%
典型位置	模型Ⅰ	模型Ⅱ	模型Ⅲ	模型Ⅳ	模型Ⅱ	模型Ⅲ	模型Ⅳ
A	328.1	231.5	53.5	24.2	29.4	83.6	92.6
B	76.9	44.9	31.7	25.2	43.9	58.8	67.2
C	201.5	108.2	132.0	83.1	46.3	34.5	58.8

下载: 导出CSV

[1]	马俊海, 施力琳.美军常规武器试验靶场手册[M].北京:国防工业出版社, 2014.
[2]	王益森, 刘徽伯.弹药试验鉴定技术[M].北京:国防工业出版社, 2001.
[3]	Saito T, Marumoto M, Yamashita H. Experimental and numerical studies of underwater shock wave attenuation[J].Shock Waves, 2003, 13(2):139-148. doi: 10.1007/s00193-003-0201-6
[4]	Kira A, Fujita M, Itoh S. Underwater explosion of spherical explosive[J]. Journal of Materials Processing Technology, 1999, 85(1):64-71. https://www.sciencedirect.com/science/article/pii/S092401369800257X
[5]	樊自建, 沈兆武, 马宏昊, 等.空气隔层对水中冲击波衰减效果的实验研究[J].中国科技大学学报, 2007, 37(10):1306-1311. http://d.old.wanfangdata.com.cn/Periodical/zgkxjsdxxb200710025 Fan Zijian, Shen Zhaowu, Ma Honghao, et al. Experimental study on attenuation of underwater shock wave by air interlayer[J]. Journal of University of Science and Technology of China, 2007, 37(10):1306-1311. http://d.old.wanfangdata.com.cn/Periodical/zgkxjsdxxb200710025
[6]	周睿, 冯顺山.气泡帷幕对水中冲击波峰值压力衰减特征的研究[J].工程爆破, 2001, 7(2):13-17. doi: 10.3969/j.issn.1006-7051.2001.02.004 Zhou Rui, Feng Shunshan. Study on weaking peck pressure of underwater shock wave by bubble curtain[J]. Engineer Blasting, 2001, 7(2):13-17. doi: 10.3969/j.issn.1006-7051.2001.02.004
[7]	张志波, 李春军, 李红勇, 等.气泡帷幕在水下爆破减震工程中的应用[J].爆破, 2003, 20(2):75-76. doi: 10.3963/j.issn.1001-487X.2003.02.028 Zhang Zhibo, Li Chunjun, Li Hongyong, et al. Application of air bubble purdah in the damping measure in the underwater blasting[J]. Blasting, 2003, 20(2):75-76. doi: 10.3963/j.issn.1001-487X.2003.02.028
[8]	顾文彬, 胡亚峰, 徐浩铭, 等.复合结构防爆罐抗爆特性的数值模拟[J].含能材料, 2014, 22(3):325-331. doi: 10.3969/j.issn.1006-9941.2014.03.010 Gu Wenbin, Hu Yafeng, Xu Haoming, et al. Numerical simulation of blast resistant characteristics for the composite structure anti-explosion container[J]. Chinese Journal of Energetic Materials, 2014, 22(3):325-331. doi: 10.3969/j.issn.1006-9941.2014.03.010
[9]	张秀华, 张春巍, 段忠东.爆炸荷载作用下钢框架柱冲击响应与破坏模式的数值模拟研究[J].沈阳建筑大学学报(自然科学版), 2009, 25(4):656-661. http://www.cnki.com.cn/Article/CJFDTOTAL-SYJZ200904009.htm Zhang Xiuhua, Zhang Chunwei, Duan Zhongdong, et al. Numerical simulation on impact responses and failure modes of steel frame structural columns subject to blast loads[J]. Journal of Shenyang Jianzhu University (Natural Science), 2009, 25(4):656-661. http://www.cnki.com.cn/Article/CJFDTOTAL-SYJZ200904009.htm
[10]	时党勇, 李裕春, 张国民.基于ANSYS/LS-DYNA8.1进行显示动力分析[M].北京:清华大学出版社, 2002.
[11]	LSTC. LS-DYNA keyword user's manual[R]. California: Livermore Software Technology Corporation, 2003.
[12]	P.库尔.水下爆炸[M].北京:国防工业出版社, 1960.
[13]	黄芬.压力动态标准方法研究[D].南京: 南京理工大学, 2003. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y624741
[14]	伍俊, 庄铁栓, 闰鹏, 等.水中爆炸实验装置结构设计与防护研究[J].振动与冲击, 2013, 32(11):131-136. doi: 10.3969/j.issn.1000-3835.2013.11.026 Wu Jun, Zhuang Tieshuan, Yan Peng, et al. Structural design of a test facility for underwater explosion and its protection measure to reduce shock wave[J]. Journal of Vibration and Shock, 2013, 32(11):131-136. doi: 10.3969/j.issn.1000-3835.2013.11.026
[15]	徐秉业, 刘信声.应用弹塑性力学[M].北京:清华大学出版社, 1995.

期刊类型引用(1)

1.

姚梦圆，董春亮，郝建平，谢雅. 基于Poyting-Thomson的岩石黏弹塑性时效损伤蠕变研究. 河南城建学院学报. 2024(06): 25-34 .

百度学术

其他类型引用(0)