爆炸作用下围岩稳定性样条小波三维数值模拟

孙惠香, 胡 平, 梁丽娜, 刘远飞, 康 婷, 于 媛

(1.空军工程大学 航空航天工程学院，西安 710038；2.兰空房管处，兰州 730020)



爆炸作用下围岩稳定性样条小波三维数值模拟

孙惠香1, 胡 平2, 梁丽娜1, 刘远飞1, 康 婷1, 于 媛1

(1.空军工程大学 航空航天工程学院，西安 710038；2.兰空房管处，兰州 730020)

爆炸具有高压瞬时的特点，地下结构模拟属于半无限域，应用传统有限元模拟时单元划分较多，计算效率低。样条小波有限元具有良好的数值逼近性，一个单元具有多个节点。推导了三维小波转换矩阵，应用构造了三维区间B样条圆环、扇形小波单元。结合工程实例通过Matlab软件编程，对爆炸作用地下拱形结构围岩稳定性进行了三维数值模拟，模拟结果与现场采集数据基本一致，计算效率高于传统有限元。

爆炸作用；地下拱形结构；围岩稳定性；区间B样条小波；单元构造

岩石地下结构地质条件复杂，爆炸作用下地下结构围岩稳定性分析是地下工程开挖的重要研究课题，也是一个非常复杂的力学变化过程，其中包含了多种非线性。由于爆炸荷载具有高温高压瞬时的特点，爆炸应力波加载速度快，加载梯度大；混凝土与岩石的应力应变关系均有下降段，这些都是地下工程的奇异性，传统有限元在下降段计算困难，而小波有限元在计算中刚度矩阵非奇异，可以求逆，此外，小波函数具有多尺度、多分辨率与紧支性等特性，很好的解决了奇异性问题。

近年来小波有限元得到迅速发展，其尤其适合于解决工程中的奇异性问题[1]。区间B样条小波是样条函数和小波的结合，既具有样条函数高逼近精度和高效率的特点，又兼有小波多分辩等特性，一个区间B样条小波单元具有多个内部和边界节点，因此在计算中可以用较少的单元获得较高的精度。

近年来，小波有限元数值模拟取得了丰硕的研究成果。2003年，法国学者CHEN 首先构造了样条小波单元[1]，在此基础上，美国学者CHUI构造了[0，1]区间B样条小波，生成了有限区间上的多分辨分析，并给出了快速分解和重构算法。陈雪峰等[4-7]详细研究了区间B样条小波，构造了一系列小波单元，并在纸张大变形、温度场突变大梯度和裂纹奇异性等问题中进行了大量的应用。关履泰研究了B样条小波的性质及分解重构算法。

目前，三维小波单元构造及数值模拟研究还很少。本文通过构造三维区间B样条小波单元，并将其应用在爆炸作用下地下围岩稳定性动力分析中。

1 三维区间B样条小波单元构造

地下结构一般为直墙拱结构，对于围岩可以划分为三维区间B样条扇形小波单元，直墙部分划分为六面体实体单元。

u(ξ,η,ζ)=Φae

(1)

式中，Φ为小波插值基函数。

Φ=Φ1⊗Φ2⊗Φ3

(2)

式中，Φ1,Φ2,和Φ3分别为m阶j尺度下的一维小波尺度函数。ae为待求的小波插值系数列向量。令ξi=(i-1)/n；ηi=(i-1)/n；ζi=(i-1)/n；i=1,2,…,n+1为标准求解域中各节点的坐标值，设物理自由度列向量。

图1 三维区间B样条圆形扇形小波单元Fig.1 Three-way interval B-spline annulus sector wavelet element

(3)

可以得到

ue=Reae

(4)

式中

(5)

式中:

(6)

令

u(ξ,η,ζ)=Neue

(7)

由式(1)和式(4)，可得到三维C0型转换矩阵

Te=(Re)-1

(8)

形函数

Ne=ΦTe

(9)

1.1 能量泛函分析

三维柱坐标几何方程为[8]

ε=LV=

(10)

物理方程为

σ=[σrσθσzτrθτθzτzr]=Dε

(11)

式中:三维柱状坐标的本构矩阵为

(12)

设求解域上的体分布力为f，面分布力为q，集中力为Fi。由于爆炸荷载属于强动载，因此结构分析时，必须考虑其动力效应，设C为阻尼黏滞矩阵，c为黏滞系数；M为拱的质量矩阵，ρ为岩石密度，则某一瞬时单元的总能量泛函为

∫ΩuTfdv-∫sΩuTqds-∑Fiui=

(13)

1.2 运动微分方程和小波单元刚度推导

对于三维拱形结构来说，要对半径r，角度θ，轴线z方向位移分别独立插值，分别令

(14)

将式(7)、(9)代入(10)则

(15)

将单元求解域映射到标准求解域Ω(0,1)，将式(11)(12)、(15)代入式(13)，由瞬时变分原理[8]令δΠp=0可以得到运动微分方程

(16)

式中:

Ke=lerleθlez∫Ω(0,1)BTDBdξdηdζ

(17)

Fe=lerleθlez(∫Ω(0,1)NeTfdξdηdζ+

∫sΩqNeTdξdηdζ)+∑iNeTFi

(18)

Me=lerleθlez∫Ω(0,1)ρNeTNedξdηdζ

(19)

Ce=lerleθlez∫Ω(0,1)cNeTNdξddηdζ

(20)

式中:

(21)

经计算可以得到刚度矩阵每一元素的工程显式如下。

当本构矩阵选为非线性本构矩阵时，该方法可用于结构的非线性分析中。

2 围岩稳定性研究

2.1 数值模拟

数值模拟背景为某军事工程，该工程为直墙拱结构，以中等风化花岗岩为主，内摩擦角为30°，坚硬系数为6。模拟段岩体以Ⅱ级围岩为主体，隧道内轮廓跨度为14.5 m，拱矢高度为5.0 m，直墙高4m。埋深50 m，衬砌结构为1 m的钢筋混凝土结构，混凝土标号为C40，钢筋为HRB400级钢筋。参数见表1。炸药采用空气间隔装药，设计单位耗药量为0.68 kg/m3硝铵炸药，密度为1 630 kg/m3，爆速为6 717 m/s，试分析开挖过程中已支护和开挖围岩的稳定性。

表1 材料参数

区间B样条尺度函数为

(22)

图2 模拟模型Fig.2 The model of simulation

取硐室四周围岩各1倍跨度，长度为10 m的结构进行三维模拟，由于结构对称，取一半进行建模，见图2。用阶数m=2，尺度j=3的尺度函数作为插值函数，围岩和支护结构划分为6个圆环扇形小波单元，未开挖岩体划分为2个扇形小波单元，直墙部分划分为三维区间B样条实体小波单元，具体参见作者在文献[11]构造的实体单元。节点为等间距排列，阶数m=2，尺度j=3的尺度函数有9个，边界小波有2个，内部尺度函数为7个[1]。采用瑞雷阻尼，岩石中的爆炸荷载取突加线性荷载，为了简化模拟，将爆炸冲击波等效成节点荷载，在开挖部位施加，采用台阶法分两部开挖。岩石孔壁中的冲击峰值压力采用式(23)计算

(23)

式中:ρ0为炸药密度，Dv为炸药爆速，dc为炮孔直径，lc为炮孔长度，db为药卷直径，lb为装药长度，n为不耦合装药增大系数，一般取n=8～11。

荷载衰减时间按式(24)计算

τ=2i/ΔPm

(24)

式中：τ为等冲量作用时间，i为冲击波冲量。

选用Plastic-Kinematic硬化弹塑性本构模型，考虑围岩的非线性特征。钢筋混凝土单元采用整体式模型，钢筋配筋率为ρx=0.10，ρy=0.10，ρz=0.06，钢筋的本构矩阵Ds，按式(25)计算。

(25)

运动方程用增量法求解[9]，应用Matlab软件编程进行了数值计算。

2.2 现场测试

本工程爆破开挖时，利用应变仪和8通道动态数据采集仪组成数据采集系统，如图3，通过埋设压力传感器如图4，进行了现场爆破数据采集，埋设位置在拱顶，拱肩，直墙顶和直墙底部如图5。

图3 数据采集系统 图4 压力传感器

Fig.3 Data acquisition system Fig.4 Pressure sensor

图5 传感器埋设位置Fig.5 The embedding location of pressure sensor

爆破开挖时的震动对已支护结构围岩的稳定性产生威胁，图6为爆破开挖时，分别在现场采集和数值模拟得到的围岩压力时程曲线。有图可以看出，数值模拟最大围岩压力出现的时间和压力大小与现场采集到的数据基本相符。

图6 支护结构拱顶压力时程曲线Fig.6 The pressure-time curve ofsupporting structure roof

单位耗药量为0.68 kg/m3爆破开挖时，开挖处围岩各部位的围岩最大压力与现场采集数据对比如表2所示，可以看出，模拟结果与现场采集的最大围岩压力基本相符，误差均在10%以内。

表2 小波模拟和现场采集压力数据对比(MPa)

2.3 围岩稳定性分析

施工过程中，炸药的用量随意性很大，不同工程，同一工程不同开挖时间单位耗药量随意性较大，炸药多少会直接影响爆破效果和围岩稳定性[11-12]，本工程进行了多种耗药量模拟计算，得到了围岩不同部位的最大主应力如表3所示，由表可见，在炸药消耗量由0.68/m3分别增加10%、20%和40%。增大20%到0.816/m3时，根据Rankine强度准则，开挖后，进尺深度范围内岩石最大主拉应力均超过了花岗岩的抗拉强度，因此，围岩会产生受拉裂缝，但是最大主压强度均较低，竖向位移稳定，见图8所示，因此，围岩是稳定的，而炸药消耗量增大40%到0.952/m3时，围岩主应力呈非线性增长，主应力显著增大，拱顶主压应力超过岩石抗压强度，竖向位移持续增大，因此，拱顶围岩不稳定，会出现塌方或超挖现象。

表3 不同炸药消耗量主应力(MPa)

图7 拱顶竖向位移时程曲线Fig.7 The vertical displacement-time curve of roof

3 结 论

选用小波尺度函数作为插值函数，通过能量泛函分析和瞬时变分原理，构造了三维区间B样条圆环和扇形两种小波单元，应用构造的单元对爆破开挖时大跨地下拱形围岩稳定性进行数值模拟，同时进行了爆破开挖现场动态数据采集，与模拟结果对比可见：

(1) 小波有限元数值模拟结果与现场监测数据基本一致，表明单元构造是正确的，模拟结果是可信的。小波有限元可以应用于爆破开挖时地下围岩稳定性的数值模拟。

(2) 构造单元具有多节点，模拟中8个小波单元相当于传统有限元4 000余个单元。其刚度矩阵非奇异，计算效率高，模拟中计算效率要比传统有限元高40左右%。因此，小波有限元适用于爆炸强动载作用下的围岩稳定性动态分析。

(3) 岩体爆破开挖时，单位耗药量由0.68/m3提高40%到0.952/m3时，拱顶围岩最大主压应力和最大主拉应力均超过岩石的抗压和抗拉强度，竖向位移不断增大。表明随意提高单位炸药消耗量将威胁开挖洞室围岩稳定性，会出现超挖或塌方现象。

[1] 何正嘉，陈雪峰等.小波有限元理论及其工程应用[M].北京：科学出版社，2006，30.

[2] 秦荣.计算结构力学[M].北京：科学出版社，2003.

[3] A practical guide to splines[M].World Book Publishing Company, 2000.

[4] 陈雪峰，向家伟，何正嘉.区间B样条小波薄壳截锥单元构造[J].应用力学学报，2007 (12)：599-603.

CHEN Xuefeng,XIANG Jiawei, HE Zhengjia. Construction of thin truncated conical shell elements by inteval b-spline wavelet[J]. Cinese Journal of Applied Mechanics,2007(12):599-603.

[5] 向家伟，陈雪峰,李兵,等.一维区间B样条小波单元的构造研究[J].应用力学学报，2006(6)：222-227.

XIANG Jiawei, CHEN Xuefeng,LI Bing,et al. Construction of one-dimensional elements with B-Spline wavelet[J]. Cinese Journal of Applied Mechanics,2006(6):222-227.

[6] GOSWAMI J C, CHAN A K, CHUI C K. On solving first-kindintegral equations using wavelets on a bounded interval[J].IEEE Transactions on Antennas and Propagation, 1995, 43(6):614-622.

[7] JIANG Z W.Cubic spline wavelet bases of sobolev spaces and multilevel interpolat ion[J].Applied and Computational Harmonic Analysis，1996(3)：154-163.

[8] 张雄，王天舒.计算动力学[M].北京：清华大学出版社.2007：5-29

[9] 张波，盛和太.ANSYS有限元数值分析原理与工程应用[M].北京：清华大学出版社, 2005：216-270.

[10] 王冒晟，邵敏.有限单元法基本原理和数值方法[M].北京：清华大学出版社，1997.

[11] 孙惠香，许金余,朱国富,等. 爆炸作用下跨度对地下结构破坏形态的影响[J].空军工程大学学报，2013, 14(2)：90-94.

SUN Huixiang,XU Jinyu,ZHU Guofu,et al. The influence of span for deep underground arch structure on failure modes under blast loading [J]. Journal of Air Force Engineering University, 2013， 14(2)：90-94.

[12] 孔大庆，孙惠香,康婷，等.岩体特性对围岩与结构动力相互作用影响[J].空军工程大学学报，2014， 15(6)：77-81.

KONG Daqing,SUN Huixiang,KANG Ting,et al. The Influence of rock characteristics on dynamic interaction between adjoining rock and structure subjected to blast loading[J]. Journal of Air Force Engineering University, 2014,15(6)：77-81.

Spline wavelet 3D numerical simulation for adjoining rock stability under explosion

SUN Huixiang1, HU Ping2, LIANG Lina1, LIU Yuanfei1, KANG Ting1, YU Yuan1

(1. Department of Airport Construction, Air Force Engineering University, Xian 710038, China;2. Lanzhou Air Force Housing Management Office, Lanzhou 730020, China)

Blasting load has characteristics of high pressure and transient state. Underground structure blast simulation belongs to a semi-infinite domain simulation. When using FE for its blast simulation, the computation efficiency is very low due to too many elements. The wavelet finite element method has characteristics of good numerical approach property, one element has many nodes. The 3D wavelet conversion matrix was derived. New elements of 3D interval B-spline annulus wavelet element and sector wavelet element were constructed. They were used to study the adjoining rock stability of an underground arch structure under blasting load through programming with Matlab software. The simulation results were consistent with the data collected at the blasting site. It was shown that the computation efficiency of the proposed method is improved compared with the traditional finite element method.

blast action; underground arch structure; adjoining rock stability; interval B-spline wavelet; element construction

国家自然科学基金资助项目(51208506;51308540)

2015-07-20 修改稿收到日期：2015-09-27

孙惠香 女，博士，副教授，硕士生导师，1975年生

TU45

A

10.13465/j.cnki.jvs.2016.19.005