炮孔近區巖體振動主頻及能量的數值模擬

羅 憶，鄭博聞,2，劉飛香，李新平 ，張雪屏,2

(1.武漢理工大學道路橋梁與結構工程湖北省重點實驗室，武漢430070；2.武漢理工大學土木工程與建筑學院，武漢430070；3.中國鐵建重工集團股份有限公司，長沙 410100)

在國家“一帶一路”的戰略背景下，巖體爆破技術在水利、交通、礦山等行業的應用越來越廣泛，由于巖體結構的特殊性以及爆破振動信號的復雜性、非平穩性，如何有效精準地預測和控制爆破振動對保留巖體所造成的破壞效應仍是當前研究的一大難題[1]。

單一利用爆破振動速度來控制巖體穩定性存在一些問題，現行爆破振動安全規程[2]綜合考慮振動速度和振動頻率來監測和控制巖體穩定性，而在爆破振動速度研究方面現有研究成果較多，所以近些年爆破振動頻率逐漸成為研究的熱點問題。夏祥等[3]模擬得出了質點振動主頻的衰減規律，以及振動主頻與炸藥藥量之間的關系。盧文波等[4]發現爆破地震波幅頻分布曲線存在多峰結構，爆破振動的主頻并不是隨著爆心距的增大嚴格遞減，而是會在局部出現跳躍現象。周俊汝等[5]結合實測資料發現爆破振動平均頻率隨爆心距的增大規則地衰減。楊潤強等[6]發現隨著地應力的增加，爆破振動的主頻會減小，低頻振動成分會增加。李洪濤等[7-8]對比研究了地下洞室梯段爆破與露天梯段爆破等不同類型爆源形式下爆破地震能量的分布特征。

單孔爆破作為最簡單最基礎的情況，一直是學者們研究的重點，而研究側重點主要在于不耦合系數，徐穎等[9]發現不耦合系數為1.67時的爆破裂紋長度最長，并且不耦合系數的增加會導致裂紋數目減少。王志亮等[10]發現當不耦合系數在3.0以內時，孔壁壓力、加速度、速度的下降速率較快，而大于3.0時衰減速率很小。王偉等[11]發現合理的不耦合系數，可使巖石不形成粉碎區，并且能量耗散大幅度減少。姜鵬飛等[12]發現當不耦合系數大于2.0時，不耦合系數對巖石應力及峰值速度影響不明顯。鐘明壽等[13]發現當不耦合系數為1.5時，采集的地震波信號主頻低、信噪比最高。梁為民等[14]發現耦合裝藥條件下粉碎區更大，而不耦合裝藥條件下形成的裂紋更稀疏、更長。王夏南[15]發現徑向不耦合系數在1.15～1.5內時，波峰波谷處的應力強度最強。

相對于爆破中遠區，爆破近區振動響應更加復雜，傅洪賢等[16]通過現場實測對薩道夫斯基公式進行了完善和補充。胡英國等[17]利用數值模擬手段，得出了爆破近區巖體臨界損傷峰值質點速度。謝烽等[18]得出適用于近區振動預測模型，隨著離爆源距離的增加，振速的衰減速率逐漸變小。

在爆破振動頻率方面的研究相對較少，特別是針對爆破近區方面研究更少，需要盡快揭秘爆破近區巖體爆破振動頻率傳播和衰減特性。筆者從最簡單的單孔爆破模型入手，采用ANSYS/LS-DYNA模擬研究了爆破近區的振動響應，通過改變不耦合系數，研究不耦合系數與爆破振動主頻以及能量之間的變化規律。

1 理論分析

1.1 振動主頻

在振動頻率衰減機理研究方面，國內學者[19-20]參照薩道夫斯基公式，通過量綱分析得到爆破振動主頻預測公式：

(1)

(2)

式中：f為爆破振動主頻預測值；Q為裝藥量；R為爆心距；PPV為質點峰值振動速度；K、α為與介質特性、爆破方式、傳播途徑有關的系數。

另有學者[21-22]推得的主頻預測公式分別為

(3)

(4)

式中：f為爆破振動主頻預測值；a1、a2為待定常數；Q為裝藥量；R為爆心距；Cs為巖體橫波波速,cm/s；kf為頻率系數，kf=0.01～0.03。

(5)

式中：f為頻率；i為虛數單位。

爆炸荷載的頻譜可以表示為

(6)

文獻[23]根據三角形爆炸荷載推得爆炸荷載的幅值譜為

(7)

式中：pmax為爆炸荷載形成的峰值壓力；tr為爆炸荷載達到峰值所需的上升時間；td為爆炸荷載全歷程持續時間。

徑向不耦合裝藥時，爆炸峰值荷載p0為

(8)

式中：ρ0為炸藥的裝藥密度；D為爆轟波沿炮孔方向的傳播速度；γ為等熵指數，對于一般的巖石炸藥起爆，取γ=3.0；a為炮孔的裝藥直徑；b為炮孔直徑。

1.2 振動能量

爆炸過程中，爆生氣體在炮孔中膨脹壓縮間隙中的空氣產生沖擊波，而后又沖擊炮孔壁，可認為爆生氣體膨脹充滿整個炮孔。在該過程中，爆生氣體的體積增大，密度減小并且音速也降低，從而導致其波阻抗發生變化。

文獻[24-25]認為爆生氣體膨脹充滿整個炮孔過程視為等熵絕熱膨脹，其密度和音速分別為

(9)

(10)

(11)

將式(9)和式(11)代入式(10)中，得

(12)

傳遞給巖石的能量可表示為

(13)

由式(12)可知，巖體中爆炸能量的傳遞與炸藥和巖石的波阻抗、裝藥不耦合系數有關。

用|F(ω)|2表示爆破振動信號中不同頻率成分所對應的功率譜密度，將其定義為功率譜函數，功率譜密度PSD表示了一定頻率對應諧波分量的能量的相對大小。

爆破發生時，對于空間中質量為Δm的質元在t時刻的動能可以表示為：

(14)

式中：E(t)為爆破振動在t時刻的振動能量;v(t)為質元在t時刻的振動速度。

對質元質量進行歸一化處理，在振動時程t1→t2內對質元振動速度進行積分，得到爆破振動信號的總能量Eq：

(15)

式中：Eq為爆破振動的總能量；t1，t2分別為爆破振動信號記錄的起止時刻。

由于爆破振動采集到的是一系列離散值，所以上式可以表示為

(16)

式中：v(ti)為離散的爆破振動速度采樣序列；Δt為采樣時間間隔。

2 單孔爆破數值模擬

2.1 數值模型

采用ANSYS/LS-DYNA有限元動力分析軟件，建立了爆破二維平面有限元模型。模型長10 m、寬10 m，在模型四周設置無反射邊界，模型左側施加位移約束，炮孔位于模型左下角(見圖1)，炸藥直徑32 mm，采用不耦合裝藥，不耦合系數K分別取1.3、2、2.5、3、3.5五種不同工況對比分析，在模型對角線上間隔1 m依次選取10個測點。

圖1 有限元數值模型Fig.1 Finite element numerical model

在模擬計算中，采用流固耦合的計算方法，用ALE算法將炸藥、空氣、巖體三者相耦合，采用JWL狀態方程[26]描述爆炸過程，其狀態方程表達式如下：

(17)

式中：V為爆生氣體相對體積；A、B、R1、R2、ω、E均為狀態方程材料參數。

JWL狀態方程對應參數如表1所示，其中：A、B、R1、R2、ω為與炸藥相關的材料參數；E0為初始內能比，V0為初始相對體積[27]。

表1 炸藥和JWL狀態方程參數

空氣介質模型狀態方程采用Gruneisen方程描述，具體材料參數如表2所示。

(γ0+αμ)E2

(18)

式中：ρ2為空氣材料密度；γ0為Gruneisen參數；C為曲線截距；α為γ0的一階體積修正；S1、S2、S3均為曲線斜率的系數；μ為體積修正量。

表2 空氣材料參數

巖體采用隨動強化本構模型，結合河南南陽天池抽水蓄能電站現場巖體參數，具體巖石參數如表3所示。

表3 巖石參數

2.2 爆心距與主頻及能量的變化規律

根據文獻[3]的研究，質點振動速度的主頻按下式定義：

(19)

式中：F1，F2分別為水平直線y=Fmax/2與速度頻譜曲線交點橫坐標的最小值和最大值；Fmax為傅里葉頻譜曲線上的極大值。

模擬計算爆心距1 m處的振速頻譜極大值Fmax為4 600 Hz，取Fmax/2為2 300 Hz，與頻譜圖得到的交點F1，F2分別為3 560、8 850 Hz，進而求得爆心距1 m處的振動主頻為6 205 Hz(見圖2)。

圖2 爆心距1 m處振動主頻Fig.2 Dominant frequency of vibration at 1 m from blast center

不同爆心距的振動速度頻譜如圖3所示。對比爆心距分別為2、3、4 m處的振速頻譜發現，振速頻譜圖在2、3、4 m處分別呈現5峰、6峰、7峰(極大值點)的結構，可見隨著爆心距的增加振速頻譜的峰數目是逐漸增大的，即頻率由某個區段集中分布逐漸變為多個區段均勻分布的形態。對比爆心距1、5、10 m三處測點，在1 m處振速頻譜呈現4峰結構，頻率主要集中在4 000～10 000 Hz，但在5、10 m處頻譜呈現多峰結構，頻率主要集中在3 000～8 000 Hz。

圖3 不同爆心距的測點頻譜Fig.3 Frequency spectrum of measuring point at different distance from blast center

根據式(19)，求出在不耦合系數K=1.3的工況下，爆心距1～10 m范圍內爆破振動速度主頻，K=1.3代表工程中常用的炮孔直徑42 mm，藥卷直徑32 mm的工況。為了便于分析主頻的衰減規律，選取10 m處的最小主頻作為基準，將主頻歸一化整理，得出主頻衰減變化規律(見圖4)，主頻在2 m處最大，3 m處主頻迅速衰減變小，3～6 m主頻基本保持穩定，在6 m處還有小幅度增大，但在6 m之后主頻又迅速直線衰減，說明了爆破振動速度主頻并不是隨著爆心距的增大而嚴格遞減的。

圖4 主頻衰減Fig.4 Dominant frequency attenuation

爆破振動速度時程曲線經傅里葉變換，獲得離散化的頻率值系列和對應的功率譜密度PSDi系列。分別對不同頻率段的功率譜進行積分，得出在一定頻帶范圍內的能量分布情況[28]。頻率范圍(fm,fn)內的爆破振動能量占總能量的比例：

(20)

式中：PEi為頻率范圍(fm,fn)內的爆破振動能量占比。

將0～15 000 Hz按照2 500 Hz為界限劃分為6個區段，選取爆心距分別為1、5、10 m處測點的能量占比分析(見圖5)，在爆心距1 m處5 000～7 500 Hz段的占比最多，達到了29.7%，但與其他頻段占比相差不大，總體呈現正態分布。在爆心距5 m處，占比在2 500～5 000 Hz和5 000～7 500 Hz較大，兩個頻段相差不大，僅相差1.8%，并且能量主要集中在2 500～10 000 Hz這個范圍內，約占總能量的80%以上，對比1 m處，該占比約僅有70%，說明隨著爆心距的增大，頻帶能量逐漸往主頻靠攏。在10 m處，能量最大占比則在2 500～5 000 Hz頻段，并且達到了37.9%，進一步說明了隨著爆心距的增大，能量往主頻段附近集中，能量分布由高頻均勻分布變為低頻集中分布。

圖5 各頻率段能量占比Fig.5 Energy ratio of each frequency segment

爆心距1、5、10 m處的峰值振速、主頻、總能量如表4所示，從1 m增至5 m處峰值振速和總能量分別衰減了73.0%和71.7%，5 m到10 m處則分別衰減了50.9%和43.1%，可見隨著爆心距的增大，峰值速度和能量迅速減小并且衰減速率是逐漸變小的。但從1 m到5 m處，振動主頻衰減了17.0%，5 m處到10 m處則衰減了18.8%，可見振動主頻的衰減速率是基本穩定的，與峰值振速和能量的衰減速率有所不同。

表4 速度、主頻、能量衰減對比

2.3 不耦合系數與主頻及能量的變化規律

選取最小主頻作為基準值，將主頻歸一化處理，同樣選取不同不耦合系數下1、5、10 m處的振動主頻(見圖6)，隨著不耦合系數的增大，振動主頻總體上呈現先增大后減小再增大的趨勢。當K<2時，振動主頻隨K的增大而增大；當22.5時，振動主頻隨K的增大而增大，并且增大速率逐漸變快。

圖6 主頻隨不耦合系數變化規律Fig.6 Variation of dominant frequency with decouple ratio

同樣將能量歸一化處理(見圖7)得出，能量隨著不耦合系數的增大先增大后減小，變化規律與峰值速度相同，當不耦合系數K=2.5時，在爆心距1、5、10 m處能量均達到了最大值，并且在1 m處能量的增減幅度更大。說明不耦合系數的改變對離炮孔較近處能量的變化影響更加明顯。

圖7 能量隨不耦合系數變化規律Fig.7 Variation law of energy with decouple ratio

各頻率段能量占比如表5所示，主頻所在的2 500～5 000 Hz和5 000～7 500 Hz兩個區段，爆心距1 m測點，隨著不耦合系數的增大主頻先增大后減小再增大，同樣主頻所在5 000～7 500 Hz區段能量也呈現相同的變化趨勢，與2 500～5 000 Hz占比的差值同樣先從5%增大到15%后減小至10%再增大到15%。在5 m測點，K=1.3時兩段能量占比相差不大，隨著不耦合系數的增大到2時，5 000～7 500 Hz段能量占比要多15%，但當K增大到2.5、3時，兩段能量占比又相差不大，當K=3.5時，兩段能量占比又拉大到15%左右，總體與主頻一樣呈現先增大后減小再增大的變化趨勢。在10 m測點，主頻位于2 500～5 000 Hz區段，K=1.3時，2 500～5 000 Hz區段能量占比比5 000～7 500 Hz多9%，不耦合系數增大到2.0時，5 000～7 500 Hz段能量占比比2 500～5 000 Hz段多約20%，發現不同于1、5 m測點，兩段的能量占比差距變得更大，總體規律與主頻變化規律一致，進一步說明了主頻與能量之間的密切關系。

表5 各頻率段能量占比

3 結論

1)在爆心距10 m范圍內，隨著爆心距的增大，爆破振動主頻在爆心距2 m處達到峰值，隨后迅速衰減，在3～6 m范圍內保持穩定，6 m之后再次迅速衰減。

2)在同一爆心距位置，隨著不耦合系數的增大，振動主頻呈現先增大后減小再增大的趨勢，而能量呈現先增大后減小的趨勢。當不耦合系數為2.5時，爆破振動主頻最小，爆破振動能量最大。

3)爆破振動能量在距離炮孔較近處在各個頻率段均勻分布，但在炮孔較遠處各頻段能量逐漸由分散變為相對集中，并且向主頻所在的頻率段集中。

本研究數值模擬計算基于各項同性假設，為了排除其他因素干擾和提高計算效率，沒有考慮巖體中的節理、裂隙等結構面的影響，在該方面還需要進一步的研究完善。所得研究結論將為后續多孔爆破研究提供理論依據和參考。