姜立春羅恩民

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

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

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

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

1　工程概況

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

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

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

2　動力響應(yīng)模型

2.1　基本假設(shè)

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

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

2.2　模型構(gòu)建

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

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

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

3　參數(shù)確定

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

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

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

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

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

式中，是振幅矩陣。

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

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

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

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

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

4　類動力響應(yīng)方程分析方法

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

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

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

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

5　數(shù)值分析模型

5.1　基本假設(shè)

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

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

5.2　數(shù)值模型構(gòu)建

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

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

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

6　動力響應(yīng)特征分析

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

6.1　位移響應(yīng)特征分析

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

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

6.2　速度響應(yīng)特征分析

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

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

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

7　結(jié)　論

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

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

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

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