礦震擾動下立體采空區群動力響應研究

姜立春羅恩民

（1.華南理工大學土木與交通學院，廣東廣州510640；2.華南理工大學安全科學與工程研究所，廣東廣州510640）

地下金屬礦山在長期開采過程中遺留大量采空區，隨著開采進行將形成采空區群[1-6]。隨著開采的深入進行，頻繁的采礦開拓作業導致采動應力疊加嚴重，并導致巖體的承載應力逐漸增加和屈服損傷日益嚴重，易引發頂板冒落、圍巖片幫、沖擊地壓等地質災害，當受地質災害引發的礦震應力波擾動時，極易改變采空區群原始的應力狀態，空區群的整體穩定性會受到弱化，嚴重影響采空區群的安全。因此，開展礦震擾動下采空區群的動態響應規律研究具有非常重要的意義。

目前國內外學者對動荷載作用下地下金屬礦山采空區群的動力響應規律開展了研究。唐禮忠等[7]采用應力波方法分析了采空區圍巖受到側崩時的動力響應，采用有限差分軟件FLAC3D模擬了某金屬礦典型采空區圍巖受到動力擾動時的動力特征；賈楠等[8]采用三維動力有限差分法分析了動力擾動下地下金屬礦山采空區圍巖非線性演化的規律；付建新等[9]利用相似材料模型試驗，驗證了在外力擾動作用下采空區頂板和礦柱應變的演變規律；W C Zhu等[10]數值模擬了動態擾動波形對地下采空區圍巖破壞區發育的動態影響；SK REDDY等[11]采用有限差分的方法在爆破作用下研究了障礙柱的應力分布規律。總體來說，目前金屬礦山采空區動力響應問題研究尚處于起步階段，主要研究方法為數值模擬，動態擾動主要指的是爆破荷載的擾動。然而，針對礦震應力波擾動作用下采空群動力響應的研究尚不多見，亟待開展該領域研究方法，減少巨大工程災害發生的概率。

以往關于采空區受到礦震作用的研究，主要從單元采空區和礦震機理的角度進行研究[12-13]，關于采空區群在礦震作用下的動態響應規律研究目前較少，而采空區群又在礦山地質中頻頻出現。為此，筆者在前期工作中借助結構動力學的離散化手段，提出水平和縱向采空區群動力響應模型分析方法，研究單行和單列采空區群的動力響應問題，取得良好成效[13-15]。在此基礎上，以某金屬礦山8號礦體1 100～1 250 m中段36～42勘探線為工程背景，構建動力響應模型，運用動力學理論，引入剪切力作用系數矩陣B，采用Newmark-β分析方法，利用Matlab工具包和自編程，研究礦震擾動下多自由度系統的動力響應規律，并與數值分析方法結果進行對比，分析該方法的合理性，進而為采空區群礦震防治提供理論指導。

1　工程概況

某金礦位于秦嶺褶皺系南秦嶺印支褶皺帶鳳縣-鎮安褶皺束的北緣，自上而下分為泥盆系中統的王家楞組（D2W）和古道嶺組（D2g），礦床賦存于古道嶺組（D2g1）下段的含金角礫巖地帶，上覆巖層和圍巖為王家楞組（D2W）和古道嶺組（D2g）的鈉長石化板巖、絹云板巖和結晶灰巖。礦體與圍巖為過渡接觸，礦體邊界與角礫巖邊界一致時，礦體與角礫巖的邊界線容易辨認[14]。礦山由KT5、KT8、KT9等6個礦體組成。

選取8號礦體1 100～1 250 m中段36～42勘探線間的3行、3列的復雜采空區群作為研究對象，如圖1所示。單元采空區跨度40 m；間柱寬10 m、高35 m；頂板厚5 m、底板厚10 m；圍巖厚度20 m。巖體物理力學參數見表1。

隨著采礦開拓作業的深入或殘采的進行，礦區內在區域應力場原始應力發生改變處于失調不穩的異常狀態，若在空區頂板處積累了一定能量后其極限應力區發生脆性破壞，導致頂板覆巖大面積垮落至空區22頂板上產生了巖層震動，造成了采空區群的礦震災害，圖2為立體采空區礦震示意圖。

2　動力響應模型

2.1　基本假設

（1）立體采空區群巖體構造完整，不考慮斷層、節理和水的滲流的影響；

（2）立體采空區群頂（底）板、間柱和圍巖未發生嚴重變形。

2.2　模型構建

將采空區群及圍巖作為一個系統進行考慮，采空區群系統離散化處理后，考慮到上覆巖層和圍巖水平地應力的影響，上覆巖層等效為自重應力。圍巖和間柱共3行，每行有4個元素；頂板共2行，每行有3個元素；系統共含18個元素。在平面剖面方向上，1個質點為1個自由度，系統共有18個自由度單元（圖3）。

根據結構動力學的D’Alembert原理，采空區系統動力響應的運動微分方程為：

式中，M、C、K、Q分別為質量矩陣、阻尼矩陣、剛度矩陣、剪切力矩陣，分別為加速度矩陣、速度矩陣和位移矩陣，λ為剪切力作用因子矩陣，B為剪切力作用系數矩陣。

3　參數確定

（1）質量矩陣M。根據文獻和現場工況，3行、3列立體采空區群各質點的質量矩陣為

（2）剛度矩陣K。根據文獻[14-16]可得剛度系數k為：

式中，E為巖體彈性模量，；A為單元采空區圍巖的水平截面積，A=1 050 m2；L為單元采空區的高度。

利用求得的剛度系數借助于Matlab自編程可求得剛度矩陣K。

（3）頻率系數ω。在Rayleigh阻尼的計算過程中需要考慮結構的振動快慢，即自振頻率ω，ω通過求解特征方程式求得：

式中，是振幅矩陣。

（4）阻尼矩陣C。計算式（6）求得模型的18階頻率系數，選取第3階頻率系數78.96 Hz和第7階頻率系數147.95 Hz，計算巖體的Rayleigh阻尼，最終獲取系統的阻尼矩陣C。Rayleigh阻尼計算式如下：

式中，和分別是質量阻尼比例系數和剛度阻尼比例系數。

工程中通常取巖體阻尼比，則和為

（5）剪切力矩陣Q。外部荷載擾動下，采空區群內部的頂（底）板、間柱與周邊圍巖接觸面的剪切應力，可通過廣義胡克定律進行求解，由文獻[14-16]可得剪切力矩陣Q。剪切應力的大小受多種因素影響，受既有工程技術水平發展的限制，精確計算剪切力存在一定的困難。本研究引入剪切力作用因子矩陣進行修正，的取值參考文獻[14]。修正后的剪切力矩陣為：

（6）礦震等效沖擊荷載矩陣P。若在空區頂板處積累了一定能量后其極限應力區發生脆性破壞，導致頂板覆巖大面積垮落至空區22頂板上產生了沖擊荷載，沖擊荷載的研究參考文獻[14]，本研究中沖擊荷載即為礦震等效沖擊荷載矩陣P，P的表達式為：

4　類動力響應方程分析方法

對于3行、3列采空區群動力的響應問題，剪切力矩陣和礦震等效沖擊荷載矩陣構成復雜，類動力響應運動響應方程計算困難，而且它是一個對時間函數非常敏感的剛性方程，采用顯示方法計算時，必須取非常小的時間步長[17]。相比之下，隱式方法對時間步長要求較低，且計算結果更具穩定性和精確性[17]。因此，本文擬采用隱式方法Newmark-β法進行[18]。計算程序如下。

（1）Newmark-β法的穩定性分析表明[17-18]，當β≥0.5，γ≥0.25（0.5+β）2時，它是穩定的（β和γ為算法參數）。設爆破荷載步長Δt=0.001，參數γ=0.5，β=0.25，Newmark-β法中各積分常數為

（4）求解時刻的位移、加速度和速度

選取計算總時間T=0.8 s，循環式（12）至式（16）計算步驟，借助于Matlab自編程對式（3）進行迭代求解，獲取爆破荷載作用下各頂（底）板、間柱和圍巖的動力響應參數位移矩陣和速度矩陣。

5　數值分析模型

5.1　基本假設

（1）立體采空區群巖體構造完整，不考慮斷層、節理和水的滲流的影響

（2）立體采空區群巖體具有均質各向同性的半無限體。

5.2　數值模型構建

根據礦山現場實際情況結合圣維南定理合理簡化邊界，利用Midas/GTS軟件構建采空區3行3列的區域數值模型（如圖4）。該模型長為380 m，寬為30 m，左側高為400 m，右側高為300 m。

在數值模擬過程中，巖體破壞準則服從摩爾-庫侖（Mohr-Coulomb）準則，僅考慮巖體的初始地應力，巖體力學參數參考表1。數值模擬過程中，為避免爆炸振動波在模型邊界反射回來影響結果，在模型左右、底部設置粘性邊界，前后設置位移條件約束。為了提高數值模擬的精度，對模型進行了網格細分（圖4）。

礦震等效荷載施加于如圖2在空區22頂板的上向臨空面上。

6　動力響應特征分析

利用動力響應模型法計算時，由于動力響應模型是將采空區的頂（底）板和間柱經過離散化處理得到的18個自由度單元質體，每一條質體對應一個動力響應時程曲線，得出動力響應模型有18條動力響應時程曲線。而數值模型法計算時，要經過特征值分析和時程分析獲得了動力響應特征參數，但是，數值模型的頂（底）板和間柱是連續分布的，不同節點對應不同的動力響應時程曲線。為了與動力響應模型對比分析，需要對數值模型的頂（底）板和間柱的數值模擬結果進行加權平均處理。

6.1　位移響應特征分析

圖5所示動態響應模型法和數值模型法得到的位移響應時程曲線對比分析可以看出它們的變化趨勢趨于一致，峰值都處于同一數量級。礦震沖擊荷載激勵瞬間，采空區的頂（底）板和間柱的位移瞬時到達峰值（見圖5（a）～（f）），直接作用巖體的位移是在塑性區域內發生有規律的振動（見圖5（c））。相鄰中段的頂（底）板和間柱的位移極值隨著與礦震沖擊荷載的距離增大而減小。由于巖體內部阻尼的能量耗散和剪切力的共同作用，位移的幅值隨時間增大而逐漸衰減，最后實現平衡。

圖5（c）中，在t=0.0～0.1 s時，巖體振動強烈，振幅較大；t=0.1～0.2 s時巖體位移的振幅忽高忽低，位移時程曲線出現跳躍式的衰減；t≥0.2 s時，位移振幅逐漸衰減，采空區群最終恢復平衡。

6.2　速度響應特征分析

從圖6中動力響應模型法與數值模擬法的速度時程曲線對比分析可以看出它們的變化趨勢出現相似的規律。首先達到采空區的頂（底）板和間柱的位移迅速到達峰值，由于巖體內部的阻尼和剪切力的共同作用，位移的幅值隨時間增大而逐漸衰減，最后實現平衡。直接作用巖體的位移是在塑性區域內發生有規律的振動，相鄰中段的頂（底）板和間柱的位移極值隨著與礦震沖擊荷載的距離增大而減小。

圖6中，在t=0.0～0.1 s時巖體動力響應強烈；t=0.1～0.2 s時巖體速度的振幅忽高忽低，速度時程曲線出現跳躍式的衰減；t≥0.3 s速度響應振幅逐漸衰減，采空區群最終恢復平衡。

綜合圖5和圖6分析可知，2種方法存在一定的時間差，主要原因歸結于動力響應模型法的能量傳遞瞬間完成，而數值模擬法的能量是逐步完成的。也可看出直接作用巖體的位移是在塑性區域內發生有規律的振動。同時，動力響應模型法和數值模型法計算得出的頂（底）板和間柱的位移和速度時程曲線變化趨勢趨于一致，具有相類似的變化規律。從而在一定程度上揭示了礦震擾動下立體采空區的動力響應規律。

7　結　論

（1）提出了3行3列立體采空區群礦震擾動下動力響應模型，研究了礦震擾動下立體采空區群動力響應規律。

（2）動力響應模型法利用動力學理論推導出的類動力響應方程，其分析過程中引入剪切力作用系數矩陣B表征了剪切力，采用了Newmark-β法，并借助于Matlab自編程，提高了分析過程的效率。

（3）礦震沖擊荷載作用下，巖體振動能量主要來源于自身和相鄰近的外力荷載激勵，巖體位移增幅最大的部位是承載垮落物的頂板巖體，直接作用巖體的位移在塑性區域內發生有規律的振動。

（4）數值模型法用時5 h，而動力響應模型法只需要0.8 s，兩種方法得到的結果一致，得出動力響應模型法運算效率更加高效。該方法為研究礦震條件下立體采空區群的動力響應特征提供了一種新思路。