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

姜立春羅恩民

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

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

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

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

1　工程概況

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

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

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

2　動力響應模型

2.1　基本假設

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

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

2.2　模型構建

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

根據(jù)結構動力學的D’Alembert原理，采空區(qū)系統(tǒng)動力響應的運動微分方程為：

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

3　參數(shù)確定

（1）質量矩陣M。根據(jù)文獻和現(xiàn)場工況，3行、3列立體采空區(qū)群各質點的質量矩陣為

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

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

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

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

式中，是振幅矩陣。

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

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

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

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

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

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

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

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

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

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

5　數(shù)值分析模型

5.1　基本假設

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

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

5.2　數(shù)值模型構建

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

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

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

6　動力響應特征分析

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

6.1　位移響應特征分析

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

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

6.2　速度響應特征分析

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

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

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

7　結　論

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

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

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

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