爆炸作用下圍巖穩(wěn)定性樣條小波三維數(shù)值模擬

孫惠香, 胡 平, 梁麗娜, 劉遠(yuǎn)飛, 康 婷, 于 媛

(1.空軍工程大學(xué) 航空航天工程學(xué)院，西安 710038；2.蘭空房管處，蘭州 730020)

?

爆炸作用下圍巖穩(wěn)定性樣條小波三維數(shù)值模擬

孫惠香1, 胡 平2, 梁麗娜1, 劉遠(yuǎn)飛1, 康 婷1, 于 媛1

(1.空軍工程大學(xué) 航空航天工程學(xué)院，西安 710038；2.蘭空房管處，蘭州 730020)

爆炸具有高壓瞬時(shí)的特點(diǎn)，地下結(jié)構(gòu)模擬屬于半無(wú)限域，應(yīng)用傳統(tǒng)有限元模擬時(shí)單元?jiǎng)澐州^多，計(jì)算效率低。樣條小波有限元具有良好的數(shù)值逼近性，一個(gè)單元具有多個(gè)節(jié)點(diǎn)。推導(dǎo)了三維小波轉(zhuǎn)換矩陣，應(yīng)用構(gòu)造了三維區(qū)間B樣條圓環(huán)、扇形小波單元。結(jié)合工程實(shí)例通過(guò)Matlab軟件編程，對(duì)爆炸作用地下拱形結(jié)構(gòu)圍巖穩(wěn)定性進(jìn)行了三維數(shù)值模擬，模擬結(jié)果與現(xiàn)場(chǎng)采集數(shù)據(jù)基本一致，計(jì)算效率高于傳統(tǒng)有限元。

爆炸作用；地下拱形結(jié)構(gòu)；圍巖穩(wěn)定性；區(qū)間B樣條小波；單元構(gòu)造

巖石地下結(jié)構(gòu)地質(zhì)條件復(fù)雜，爆炸作用下地下結(jié)構(gòu)圍巖穩(wěn)定性分析是地下工程開(kāi)挖的重要研究課題，也是一個(gè)非常復(fù)雜的力學(xué)變化過(guò)程，其中包含了多種非線性。由于爆炸荷載具有高溫高壓瞬時(shí)的特點(diǎn)，爆炸應(yīng)力波加載速度快，加載梯度大；混凝土與巖石的應(yīng)力應(yīng)變關(guān)系均有下降段，這些都是地下工程的奇異性，傳統(tǒng)有限元在下降段計(jì)算困難，而小波有限元在計(jì)算中剛度矩陣非奇異，可以求逆，此外，小波函數(shù)具有多尺度、多分辨率與緊支性等特性，很好的解決了奇異性問(wèn)題。

近年來(lái)小波有限元得到迅速發(fā)展，其尤其適合于解決工程中的奇異性問(wèn)題[1]。區(qū)間B樣條小波是樣條函數(shù)和小波的結(jié)合，既具有樣條函數(shù)高逼近精度和高效率的特點(diǎn)，又兼有小波多分辯等特性，一個(gè)區(qū)間B樣條小波單元具有多個(gè)內(nèi)部和邊界節(jié)點(diǎn)，因此在計(jì)算中可以用較少的單元獲得較高的精度。

近年來(lái)，小波有限元數(shù)值模擬取得了豐碩的研究成果。2003年，法國(guó)學(xué)者CHEN 首先構(gòu)造了樣條小波單元[1]，在此基礎(chǔ)上，美國(guó)學(xué)者CHUI構(gòu)造了[0，1]區(qū)間B樣條小波，生成了有限區(qū)間上的多分辨分析，并給出了快速分解和重構(gòu)算法。陳雪峰等[4-7]詳細(xì)研究了區(qū)間B樣條小波，構(gòu)造了一系列小波單元，并在紙張大變形、溫度場(chǎng)突變大梯度和裂紋奇異性等問(wèn)題中進(jìn)行了大量的應(yīng)用。關(guān)履泰研究了B樣條小波的性質(zhì)及分解重構(gòu)算法。

目前，三維小波單元構(gòu)造及數(shù)值模擬研究還很少。本文通過(guò)構(gòu)造三維區(qū)間B樣條小波單元，并將其應(yīng)用在爆炸作用下地下圍巖穩(wěn)定性動(dòng)力分析中。

1 三維區(qū)間B樣條小波單元構(gòu)造

地下結(jié)構(gòu)一般為直墻拱結(jié)構(gòu)，對(duì)于圍巖可以劃分為三維區(qū)間B樣條扇形小波單元，直墻部分劃分為六面體實(shí)體單元。

u(ξ,η,ζ)=Φae

(1)

式中，Φ為小波插值基函數(shù)。

Φ=Φ1?Φ2?Φ3

(2)

式中，Φ1,Φ2,和Φ3分別為m階j尺度下的一維小波尺度函數(shù)。ae為待求的小波插值系數(shù)列向量。令ξi=(i-1)/n；ηi=(i-1)/n；ζi=(i-1)/n；i=1,2,…,n+1為標(biāo)準(zhǔn)求解域中各節(jié)點(diǎn)的坐標(biāo)值，設(shè)物理自由度列向量。

圖1 三維區(qū)間B樣條圓形扇形小波單元Fig.1 Three-way interval B-spline annulus sector wavelet element

(3)

可以得到

ue=Reae

(4)

式中

(5)

式中:

(6)

令

u(ξ,η,ζ)=Neue

(7)

由式(1)和式(4)，可得到三維C0型轉(zhuǎn)換矩陣

Te=(Re)-1

(8)

形函數(shù)

Ne=ΦTe

(9)

1.1 能量泛函分析

三維柱坐標(biāo)幾何方程為[8]

ε=LV=

(10)

物理方程為

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

(11)

式中:三維柱狀坐標(biāo)的本構(gòu)矩陣為

(12)

設(shè)求解域上的體分布力為f，面分布力為q，集中力為Fi。由于爆炸荷載屬于強(qiáng)動(dòng)載，因此結(jié)構(gòu)分析時(shí)，必須考慮其動(dòng)力效應(yīng)，設(shè)C為阻尼黏滯矩陣，c為黏滯系數(shù)；M為拱的質(zhì)量矩陣，ρ為巖石密度，則某一瞬時(shí)單元的總能量泛函為

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

(13)

1.2 運(yùn)動(dòng)微分方程和小波單元?jiǎng)偠韧茖?dǎo)

對(duì)于三維拱形結(jié)構(gòu)來(lái)說(shuō)，要對(duì)半徑r，角度θ，軸線z方向位移分別獨(dú)立插值，分別令

(14)

將式(7)、(9)代入(10)則

(15)

將單元求解域映射到標(biāo)準(zhǔn)求解域Ω(0,1)，將式(11)(12)、(15)代入式(13)，由瞬時(shí)變分原理[8]令δΠp=0可以得到運(yùn)動(dòng)微分方程

(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)

經(jīng)計(jì)算可以得到剛度矩陣每一元素的工程顯式如下。

當(dāng)本構(gòu)矩陣選為非線性本構(gòu)矩陣時(shí)，該方法可用于結(jié)構(gòu)的非線性分析中。

2 圍巖穩(wěn)定性研究

2.1 數(shù)值模擬

數(shù)值模擬背景為某軍事工程，該工程為直墻拱結(jié)構(gòu)，以中等風(fēng)化花崗巖為主，內(nèi)摩擦角為30°，堅(jiān)硬系數(shù)為6。模擬段巖體以Ⅱ級(jí)圍巖為主體，隧道內(nèi)輪廓跨度為14.5 m，拱矢高度為5.0 m，直墻高4m。埋深50 m，襯砌結(jié)構(gòu)為1 m的鋼筋混凝土結(jié)構(gòu)，混凝土標(biāo)號(hào)為C40，鋼筋為HRB400級(jí)鋼筋。參數(shù)見(jiàn)表1。炸藥采用空氣間隔裝藥，設(shè)計(jì)單位耗藥量為0.68 kg/m3硝銨炸藥，密度為1 630 kg/m3，爆速為6 717 m/s，試分析開(kāi)挖過(guò)程中已支護(hù)和開(kāi)挖圍巖的穩(wěn)定性。

表1 材料參數(shù)

區(qū)間B樣條尺度函數(shù)為

(22)

圖2 模擬模型Fig.2 The model of simulation

取硐室四周?chē)鷰r各1倍跨度，長(zhǎng)度為10 m的結(jié)構(gòu)進(jìn)行三維模擬，由于結(jié)構(gòu)對(duì)稱(chēng)，取一半進(jìn)行建模，見(jiàn)圖2。用階數(shù)m=2，尺度j=3的尺度函數(shù)作為插值函數(shù)，圍巖和支護(hù)結(jié)構(gòu)劃分為6個(gè)圓環(huán)扇形小波單元，未開(kāi)挖巖體劃分為2個(gè)扇形小波單元，直墻部分劃分為三維區(qū)間B樣條實(shí)體小波單元，具體參見(jiàn)作者在文獻(xiàn)[11]構(gòu)造的實(shí)體單元。節(jié)點(diǎn)為等間距排列，階數(shù)m=2，尺度j=3的尺度函數(shù)有9個(gè)，邊界小波有2個(gè)，內(nèi)部尺度函數(shù)為7個(gè)[1]。采用瑞雷阻尼，巖石中的爆炸荷載取突加線性荷載，為了簡(jiǎn)化模擬，將爆炸沖擊波等效成節(jié)點(diǎn)荷載，在開(kāi)挖部位施加，采用臺(tái)階法分兩部開(kāi)挖。巖石孔壁中的沖擊峰值壓力采用式(23)計(jì)算

(23)

式中:ρ0為炸藥密度，Dv為炸藥爆速，dc為炮孔直徑，lc為炮孔長(zhǎng)度，db為藥卷直徑，lb為裝藥長(zhǎng)度，n為不耦合裝藥增大系數(shù)，一般取n=8～11。

荷載衰減時(shí)間按式(24)計(jì)算

τ=2i/ΔPm

(24)

式中：τ為等沖量作用時(shí)間，i為沖擊波沖量。

選用Plastic-Kinematic硬化彈塑性本構(gòu)模型，考慮圍巖的非線性特征。鋼筋混凝土單元采用整體式模型，鋼筋配筋率為ρx=0.10，ρy=0.10，ρz=0.06，鋼筋的本構(gòu)矩陣Ds，按式(25)計(jì)算。

(25)

運(yùn)動(dòng)方程用增量法求解[9]，應(yīng)用Matlab軟件編程進(jìn)行了數(shù)值計(jì)算。

2.2 現(xiàn)場(chǎng)測(cè)試

本工程爆破開(kāi)挖時(shí)，利用應(yīng)變儀和8通道動(dòng)態(tài)數(shù)據(jù)采集儀組成數(shù)據(jù)采集系統(tǒng)，如圖3，通過(guò)埋設(shè)壓力傳感器如圖4，進(jìn)行了現(xiàn)場(chǎng)爆破數(shù)據(jù)采集，埋設(shè)位置在拱頂，拱肩，直墻頂和直墻底部如圖5。

圖3 數(shù)據(jù)采集系統(tǒng) 圖4 壓力傳感器

Fig.3 Data acquisition system Fig.4 Pressure sensor

圖5 傳感器埋設(shè)位置Fig.5 The embedding location of pressure sensor

爆破開(kāi)挖時(shí)的震動(dòng)對(duì)已支護(hù)結(jié)構(gòu)圍巖的穩(wěn)定性產(chǎn)生威脅，圖6為爆破開(kāi)挖時(shí)，分別在現(xiàn)場(chǎng)采集和數(shù)值模擬得到的圍巖壓力時(shí)程曲線。有圖可以看出，數(shù)值模擬最大圍巖壓力出現(xiàn)的時(shí)間和壓力大小與現(xiàn)場(chǎng)采集到的數(shù)據(jù)基本相符。

圖6 支護(hù)結(jié)構(gòu)拱頂壓力時(shí)程曲線Fig.6 The pressure-time curve ofsupporting structure roof

單位耗藥量為0.68 kg/m3爆破開(kāi)挖時(shí)，開(kāi)挖處圍巖各部位的圍巖最大壓力與現(xiàn)場(chǎng)采集數(shù)據(jù)對(duì)比如表2所示，可以看出，模擬結(jié)果與現(xiàn)場(chǎng)采集的最大圍巖壓力基本相符，誤差均在10%以內(nèi)。

表2 小波模擬和現(xiàn)場(chǎng)采集壓力數(shù)據(jù)對(duì)比(MPa)

2.3 圍巖穩(wěn)定性分析

施工過(guò)程中，炸藥的用量隨意性很大，不同工程，同一工程不同開(kāi)挖時(shí)間單位耗藥量隨意性較大，炸藥多少會(huì)直接影響爆破效果和圍巖穩(wěn)定性[11-12]，本工程進(jìn)行了多種耗藥量模擬計(jì)算，得到了圍巖不同部位的最大主應(yīng)力如表3所示，由表可見(jiàn)，在炸藥消耗量由0.68/m3分別增加10%、20%和40%。增大20%到0.816/m3時(shí)，根據(jù)Rankine強(qiáng)度準(zhǔn)則，開(kāi)挖后，進(jìn)尺深度范圍內(nèi)巖石最大主拉應(yīng)力均超過(guò)了花崗巖的抗拉強(qiáng)度，因此，圍巖會(huì)產(chǎn)生受拉裂縫，但是最大主壓強(qiáng)度均較低，豎向位移穩(wěn)定，見(jiàn)圖8所示，因此，圍巖是穩(wěn)定的，而炸藥消耗量增大40%到0.952/m3時(shí)，圍巖主應(yīng)力呈非線性增長(zhǎng)，主應(yīng)力顯著增大，拱頂主壓應(yīng)力超過(guò)巖石抗壓強(qiáng)度，豎向位移持續(xù)增大，因此，拱頂圍巖不穩(wěn)定，會(huì)出現(xiàn)塌方或超挖現(xiàn)象。

表3 不同炸藥消耗量主應(yīng)力(MPa)

圖7 拱頂豎向位移時(shí)程曲線Fig.7 The vertical displacement-time curve of roof

3 結(jié) 論

選用小波尺度函數(shù)作為插值函數(shù)，通過(guò)能量泛函分析和瞬時(shí)變分原理，構(gòu)造了三維區(qū)間B樣條圓環(huán)和扇形兩種小波單元，應(yīng)用構(gòu)造的單元對(duì)爆破開(kāi)挖時(shí)大跨地下拱形圍巖穩(wěn)定性進(jìn)行數(shù)值模擬，同時(shí)進(jìn)行了爆破開(kāi)挖現(xiàn)場(chǎng)動(dòng)態(tài)數(shù)據(jù)采集，與模擬結(jié)果對(duì)比可見(jiàn)：

(1) 小波有限元數(shù)值模擬結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)基本一致，表明單元構(gòu)造是正確的，模擬結(jié)果是可信的。小波有限元可以應(yīng)用于爆破開(kāi)挖時(shí)地下圍巖穩(wěn)定性的數(shù)值模擬。

(2) 構(gòu)造單元具有多節(jié)點(diǎn)，模擬中8個(gè)小波單元相當(dāng)于傳統(tǒng)有限元4 000余個(gè)單元。其剛度矩陣非奇異，計(jì)算效率高，模擬中計(jì)算效率要比傳統(tǒng)有限元高40左右%。因此，小波有限元適用于爆炸強(qiáng)動(dòng)載作用下的圍巖穩(wěn)定性動(dòng)態(tài)分析。

(3) 巖體爆破開(kāi)挖時(shí)，單位耗藥量由0.68/m3提高40%到0.952/m3時(shí)，拱頂圍巖最大主壓應(yīng)力和最大主拉應(yīng)力均超過(guò)巖石的抗壓和抗拉強(qiáng)度，豎向位移不斷增大。表明隨意提高單位炸藥消耗量將威脅開(kāi)挖洞室圍巖穩(wěn)定性，會(huì)出現(xiàn)超挖或塌方現(xiàn)象。

[1] 何正嘉，陳雪峰等.小波有限元理論及其工程應(yīng)用[M].北京：科學(xué)出版社，2006，30.

[2] 秦榮.計(jì)算結(jié)構(gòu)力學(xué)[M].北京：科學(xué)出版社，2003.

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

[4] 陳雪峰，向家偉，何正嘉.區(qū)間B樣條小波薄殼截錐單元構(gòu)造[J].應(yīng)用力學(xué)學(xué)報(bào)，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] 向家偉，陳雪峰,李兵,等.一維區(qū)間B樣條小波單元的構(gòu)造研究[J].應(yīng)用力學(xué)學(xué)報(bào)，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] 張雄，王天舒.計(jì)算動(dòng)力學(xué)[M].北京：清華大學(xué)出版社.2007：5-29

[9] 張波，盛和太.ANSYS有限元數(shù)值分析原理與工程應(yīng)用[M].北京：清華大學(xué)出版社, 2005：216-270.

[10] 王冒晟，邵敏.有限單元法基本原理和數(shù)值方法[M].北京：清華大學(xué)出版社，1997.

[11] 孫惠香，許金余,朱國(guó)富,等. 爆炸作用下跨度對(duì)地下結(jié)構(gòu)破壞形態(tài)的影響[J].空軍工程大學(xué)學(xué)報(bào)，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] 孔大慶，孫惠香,康婷，等.巖體特性對(duì)圍巖與結(jié)構(gòu)動(dòng)力相互作用影響[J].空軍工程大學(xué)學(xué)報(bào)，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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51208506;51308540)

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

孫惠香 女，博士，副教授，碩士生導(dǎo)師，1975年生

TU45

A

10.13465/j.cnki.jvs.2016.19.005