应力波在各向异性介质中传播及损伤数值模拟★

鞠冬冬，王 清，陈 鹏，胡少斌，蔡余康

(1.南京市江北新区公共工程建设中心,江苏 南京 210098; 2.河海大学土木与交通学院,江苏 南京 210098)

在人类的生活中经常有诸如冲击、爆炸、碰撞等强荷载作用下材料发生高应变的现象。比如地震、洪水、交通事故这样的灾难,还有矿山开采、钻爆法开挖、结构抗震等造福民生的技术方法。这些现象中诸如地震、矿山开采以及钻爆法开挖都可纳入岩石力学范畴,更具体地讲,属于岩石动力学范畴,即应变率在1×10-1～1×104s-1之间的情况[1]。而岩石动力学以应力波作为其理论基础,主要研究涵盖应力波在岩体中的传播与衰减[2-4],应力波与节理裂隙的相互作用[5-7],应力波与结构面的透、折、反射关系[8-10],应力波入射角[11-12]等多个方面。除此之外,岩石动力学所研究的另一个基本问题是岩石或岩体受冲击后的动态响应,包括变形、损伤、断裂的结果以及演化过程。目前,对应力波传播及岩石的动态响应的研究多种多样,从一开始常规的应力波的传播和衰减规律[13],到后来探究循环荷载下的岩石力学响应[14-15],进一步地研究其裂纹扩展和损伤演化的规律[16-18],并考虑不同的工况如高地应力爆破和水下深孔爆破等。模拟研究的对象也从各向同性岩石扩展到正交各向异性[19]和各向异性岩石。具体到工程应用领域,如典型的岩石动力学应用范畴——矿山开采及钻爆法破岩。这两个领域经典的核心施工流程是钻孔、装药以及爆破。然而从目前来看,还很少有人基于应力波对该过程的岩体,尤其是各向异性岩体损伤演化以及能量耗散进行研究,从而确定更合理的布孔间距和布孔方式。本文将采用数值模拟的手段对无限元边界下各向异性岩体在孔周脉冲波作用情况下岩体的损伤演化展开研究,并根据损伤程度来确定各向异性条件下合理的布孔间距和布孔角。

1 数值建模

本章将基于COMSOL MULTIPHYSICS建立围压状态下爆破花岗岩层的数值模拟模型,并给出本模型的相关控制方程、几何建模、网格划分、边界条件,以及模型所涉及的材料参数。

1.1 控制方程

本模型采用相场损伤法来描述线弹性岩石的损伤应变软化过程,其本构方程采用胡克定律:

σ=D×ε

(1)

其中,σ为弹性应力;D为考虑损伤的弹性矩阵;ε为弹性应变。

本文研究应力波在各向异性介质中的传播,故弹性矩阵D自定义为:

(2)

由于考虑了应变软化过程,引入关于损伤的裂纹相场因子,其相场损伤演化函数fd为:

fd=1-(1-φ)2

(3)

其中,φ为裂纹相场因子,与fd间存在弱函数关系约束,弱表达式为:

(4)

其中,Dd为裂纹驱动力;test()为试函数;l为长度尺寸;φx和φy分别为φ在x和y方向上的导数。同时,为确定固体的应力应变,还需要加入平衡微分方程,考虑惯性项后为:

(5)

其中,ρs为固体密度;x为固体位移;f为单位体积力;t为时间。此外,还需通过变形协调方程来根据应变确定固体位移,表达式为:

(6)

1.2 几何模型

本模型模拟爆破孔处由爆破或高能气体压裂产生的冲击应力波的岩层中的传播,故采用长宽为8 m×4 m的模型尺寸,取模型尺寸宽度的1/8为无限元域层以表示有限域内的无限空间,其物理宽度为8 944 m,爆破孔直径为160 mm,间距为Δd,爆破孔圆心相连的轴线与竖向轴的夹角为α,所加荷载为指数函数与正弦函数的乘积形式,其表达式为:

(7)

其中,a为放大系数,取10;P0为作用力大小取7 000 kN;γ为衰减率,取1.86;R为爆破孔直径;l为爆破孔深度;t0为加载时间,取2.7×10-8s。

具体模型如图1所示。

在二维平面的网格划分当中,主要有三角形网格单元和四边形网格单元两种,相比于三角形网格单元,四边形网格具有更为平顺,更适用于高度非线性问题,同等密度下网格数更少,收敛性更强等优点[20-22],因此本模型中,周边无限元域采用四边形网格,中间采用三角形网格,在接近圆孔处进行加密,如图2所示,网格数共11 336,一定程度上保证计算结果的精度。

本模型边界的水平载荷与竖向荷载皆取ρ=5 MPa,所采用的花岗岩参数如表1所示。

表1 花岗岩参数表

2 模型结果与分析

本章展现了无限元边界下应力波在花岗岩层中的传播过程以及伴随发生的损伤演化过程。通过对爆破孔间距和夹角α进行参数化扫描得到了不同工况下的模拟结果。同时令损伤大于0.9的面积作为有效破坏面积,将有效破坏面积作为爆破效果的评价指标,展示各向异性条件下孔间距和轴线转角对应力波冲击破坏岩石的影响。

2.1 应力波在各向异性岩层中的传播及损伤规律

本文研究爆破孔间距和轴线转角对应力波传播以及损伤演化的影响,因此对间距和转角进行参数化扫描,即在完全相同的模型条件及求解器设置下计算仅改变指定参数下的所有结果,具体模拟组设置如表2所示,其中峰值应力考虑到爆破时的放大效应,设为435 MPa。

表2 模拟组参数

本文将A组作为对照组,其应力波传播过程如图3所示,可以看出应力波的横向传播速度显著大于竖向传播速度,由弹性波波速计算公式v=√(E/P)可得弹性波在水平方向的相速度为7 817.4 m/s,在竖直方向上的相速度为4 792.0 m/s,符合图3中t=2×10-5s的结果。同时可以观察到水平方向的应力远大于竖直方向上的应力,联系各向同性条件下的计算结果(如图4所示)可知,这是由弹性模量在方向上的差异性所导致的。

A组的损伤演化过程如图5所示,由于本模型采用的是理想的线弹性模型,没有考虑到岩石的脆性,因此裂隙区较不明显,对于破碎区,将损伤面积(d>0.9)作为等效的破碎面积,其等效直径为爆破孔的7倍,符合其他学者得出的破碎区为孔径2倍～10倍的结论[23],从图5中还可以看出两个爆破孔的连线中点处损伤较中点两侧更高,从时间在0.296 ms的损伤云图中可以发现中点处更早发生破碎,这是由于应力波在中点处产生了叠加,这点在下文关于孔间距的影响中会进一步说明。

2.2 轴线转角对应力波传播及损伤的影响

在各向异性的均质岩石中,应力波向不同的方向传播会产生有规律性的差异,从各组的应力结果中(如图6所示)可以看出不同轴线转角下的峰值应力相对近似,且由于距离水平边界和竖直边界的距离不等导致没有明显规律。

而观察模型中心点的峰值应力曲线(如图7所示)可以观察到不同角度下第一个波峰的峰值和出现时间不同,其出现时间随着转角的减小逐个延后,而峰值应力也逐渐减小。但之后的波峰由于次波的干扰并未展现清晰的规律,其中D组相对其他组较大是由于传播过程中产生的次波共振恰好处于中心位置。

如上所述本文将损伤d>0.9的区域作为破碎区,并得出A组—E组的损伤破碎面积,可以看出在不同转角下,一开始的破碎面积扩展是一致的,这时双孔处加载的应力波使各自临近处发生破碎,但尚未传播到邻孔处,但在一定时间后,应力波传播到邻孔位置,对彼此的应力场造成扰动,同时由于各向异性条件下在不同方向上应力波的传播速率和衰减率的差异导致破碎区的扩展产生了较为复杂的变化,但从8×10-4s的结果可以看出破碎区面积随着角度的减小而增大,值得注意的是该时间点B组面积较C组更大(见图8),然而从该点的切线斜率可以判断在更大的时间尺度上C组的最终损伤面积更大。本文认为是角度更小的情况下,双孔间弹性模量更小,其间得到损伤演化速率更快,两孔间更容易贯通。

2.3 孔间距对应力波传播及损伤的影响

爆破孔间距向来是爆破孔布置方案中的关键参数,对爆破效果有着重要影响。由图9可以看出F,G两组峰值应力较大,且出现多个波峰,而H,I两组应力较小,且波动不明显。表明孔间距较小时产生的峰值应力较大,且由于距离较近,相应传播时间更短,因此应力波的叠加效应更加明显。而从损伤演化曲线(如图10所示)可以看出破碎区大致随着距离的增加而增加,本文认为其原因为双孔应力波的相对冲击增强了能量的衰减,导致近孔间距爆破的能量利用率降低。

3 结语

本文通过数值模拟对应力波在各向异性岩层中传播展开研究与分析,得到如下结论:

1)各向异性条件下,应力波在弹性模量更大的方向上传播速率更快,且应力幅值较高。2)当轴线转角不同时,双孔中心点处应力的第一个波峰的出现时间随着转角的减小逐个延后,且峰值应力也逐渐减小。同时破碎区面积随着角度的减小而增大。3)孔间距较小时产生的峰值应力较大,应力波的叠加效应更加明显。破碎区面积大致随着距离的增加而增加。