孫惠香, 胡 平, 梁麗娜, 劉遠飛, 康 婷, 于 媛

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

?

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

孫惠香1, 胡 平2, 梁麗娜1, 劉遠飛1, 康 婷1, 于 媛1

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

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

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

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

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

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

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

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

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

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為標準求解域中各節(jié)點的坐標值，設(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 能量泛函分析

三維柱坐標幾何方程為[8]

ε=LV=

(10)

物理方程為

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

(11)

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

(12)

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

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

(13)

1.2 運動微分方程和小波單元剛度推導

對于三維拱形結(jié)構(gòu)來說，要對半徑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)

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

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

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

2.1 數(shù)值模擬

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

表1 材料參數(shù)

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

(22)

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

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

(23)

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

荷載衰減時間按式(24)計算

τ=2i/ΔPm

(24)

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

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

(25)

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

2.2 現(xiàn)場測試

本工程爆破開挖時，利用應(yīng)變儀和8通道動態(tài)數(shù)據(jù)采集儀組成數(shù)據(jù)采集系統(tǒng)，如圖3，通過埋設(shè)壓力傳感器如圖4，進行了現(xiàn)場爆破數(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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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

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

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

[4] 陳雪峰，向家偉，何正嘉.區(qū)間B樣條小波薄殼截錐單元構(gòu)造[J].應(yīng)用力學學報，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)用力學學報，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有限元數(shù)值分析原理與工程應(yīng)用[M].北京：清華大學出版社, 2005：216-270.

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

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