任紅崗，汪旭光，譚卓英，夏志遠

(1.北京礦冶科技集團有限公司，北京，100160；2.北京科技大學土木與資源工程學院，北京，100083)

在地下礦山開采中，采用空場法開采后形成了大量采空區(qū)群，其穩(wěn)定性控制是礦山生產(chǎn)過程中面臨的關鍵技術難題之一。當?shù)V柱寬度、礦柱高度、頂板荷載和巖體強度等參數(shù)處于臨界范圍時，稍有變化將可能引發(fā)采空區(qū)失穩(wěn)和坍塌，進而誘發(fā)礦區(qū)大規(guī)?？逅瑢Φ刭|(zhì)環(huán)境、礦山安全生產(chǎn)造成極大破壞和威脅，因此，研究采空區(qū)群的失穩(wěn)機理具有重要意義[1-4]。巖石是一種能夠儲備高應變能的材料，地下被開采后，采空區(qū)為巖石能量的釋放提供了條件。巖石所積累的應變能和勢能隨著開采不斷變化，當這種能量平衡被打破時，將導致突發(fā)性失穩(wěn)[5-7]。在開采三維空間、開采時間和應力變化的五維環(huán)境下，由頂板和礦柱構成的采空區(qū)群系統(tǒng)的破壞涉及蠕變和突變，是一個復雜的非線性力學問題。研究巖體力學為主的采空區(qū)穩(wěn)定性問題，可以用多種多樣的力學理論和方法分析，如流變力學、斷裂力學、非連續(xù)介質(zhì)力學、彈塑性理論、灰色理論、耗散結(jié)構理論和系統(tǒng)工程理論等。這些理論中突變理論應用相對較多，如譚毅等[8]基于損傷力學和突變理論，建立了條帶煤柱開采形成的采空區(qū)群失穩(wěn)模型；徐恒等[9]建立了頂板結(jié)構的尖點突變力學模型，計算充填體下的采空區(qū)穩(wěn)定性及失穩(wěn)機制。夏開宗等[10]基于非線性力學，構建了礦柱-護頂層支撐體系發(fā)生破壞的力學模型；王金安等[11-12]基于巖體流變力學理論，建立了采空區(qū)流變力學模型，揭示了采空區(qū)突變和失穩(wěn)的力學機理與過程。以往研究采空區(qū)穩(wěn)定性時，較多應用二維平面力學，不能直觀表達三維真實情況，且多集中于研究單個或幾個采空區(qū)，缺少對采空區(qū)群的整體穩(wěn)定性研究，未能全面、精準反映采空區(qū)群的動態(tài)破壞過程。為此，本文在前人基礎上，依據(jù)彈性力學和蠕變突變力學理論，建立礦柱頂板為系統(tǒng)的采空區(qū)群模型，在三維環(huán)境下，將礦柱-頂板系統(tǒng)等效為一系列Burger體蠕變力學模型，應用突變理論分析失穩(wěn)機制，得到在不同采空區(qū)群參數(shù)下，礦柱寬度、礦柱高度、頂板荷載和巖體強度等主控因子影響特性，并總結(jié)各主控因子臨界值分布規(guī)律，在礦柱蠕變作用下，可以對采空區(qū)穩(wěn)定時限進行預測，為采空區(qū)群穩(wěn)定性分析、識別和控制提供了量化參考，是實現(xiàn)礦山安全、高效、經(jīng)濟開采不可缺少的關鍵數(shù)據(jù)。

1 采空區(qū)群力學模型

1.1 雙向分層條式充填采礦法

雙向分層條式充填采礦法是開采中厚、層狀、緩傾斜和破碎礦體的一種新的方法，圖1所示為雙向分層條式充填采礦示意圖。由圖1可見：礦塊在水平方向上可劃分成條帶采場，縱向可劃分成2層或若干層，每個采場分上、下步驟開采，上層礦體開采完之后要對頂板進行支護，礦石由鏟運機經(jīng)位于采場兩端的上部巷道運出，隨后開采下層礦體，礦石由鏟運機經(jīng)位于采場中間的下部巷道運出。待2層礦開采完之后，采用廢石膠結(jié)充填采場并接頂。采場間開采順序通常按照隔一采一的方式，其開采工藝可簡單地概括為切頂、護頂、降底和充填4個工序。采用雙向分層條式充填采礦法在第1步驟開采后，采場和礦柱在空間上呈長條狀分布，形成由一系列礦柱支撐頂板的系統(tǒng)，在縱向上形成采場、礦柱交錯布置矩形結(jié)構。

圖1 雙向分層條式充填采礦法示意圖Fig.1 Sketch of slicing and strip filling mining method

1.2 力學模型構建

1.2.1 采場力學模型

礦房礦柱呈條形依次間隔布置，開采時，通常將若干個條形礦房礦柱組合成1個開采單元，開采后，形成由若干個礦柱支撐頂板的系統(tǒng)。在三維空間上，礦柱可等效為彈性體，采場頂板可等效為彈性薄板，系統(tǒng)受力主要考慮頂板上方的垂直應力。設采場頂板長和寬分別為2a和2b，采場高度為H，礦柱頂板力學模型體系如圖2所示。

假設條件：1)采場呈規(guī)整的矩形，由若干個采場和礦柱組成1個系統(tǒng)，且采場礦柱分布均勻；2)采場頂板長度大于其寬度，即b＞a；3)采場頂板所受荷載相對恒定；4)采空區(qū)底板足夠穩(wěn)定，忽略其底板的微小位移。

圖2 頂板-礦柱力學模型體系Fig.2 Mechanical model system of roof and pillar

Burgers體蠕變模型能夠描述彈性應變、衰減和穩(wěn)態(tài)蠕變的力學模型，將礦柱系統(tǒng)等效為一系列Burgers體支撐的模型，建立頂板-礦柱彈性模型如圖3所示，頂板上覆巖層的均布荷載為q0，頂板密度為ρ，抗拉強度為σT，Burgers體元件簡寫為B。

圖3 頂板-礦柱彈性模型Fig.3 Elastic model of roof and pillar

頂板-礦柱系統(tǒng)模型的控制方程為[11]

式中：D為頂板剛度，v為泊松比；w為頂板的下沉撓度；ξ為礦柱支撐頂板的面積與頂板總面積的比率，稱為面積比率；σ為礦柱應力。

1.2.2 頂板撓度算法

礦房開采后，采空區(qū)群發(fā)生破壞一般要經(jīng)歷3個階段[13]，每個階段支撐頂板模型如下。

階段I：固支邊界模型。在開采初始條件下，礦柱和頂板均未發(fā)生破壞，此時可視為固支模型，頂板的撓度曲方程為[6]：

式中：w(x,y)為頂板的下沉撓度；w0為頂板中心下沉撓度；φ(x,y)為頂板下沉撓度的比例函數(shù)。

階段Ⅱ：簡支邊界模型。隨著礦柱承載力逐漸下降，頂板下沉位移逐漸變大，進而導致頂板邊緣發(fā)生塑性變形。當頂板邊緣開始發(fā)生破壞時，固支模型轉(zhuǎn)化為簡支模型，此時有

頂板邊緣破壞條件為

式中：u1為系統(tǒng)進入第I階段后頂板中心下沉位移；h為頂板的厚度；σT為頂板抗拉強度；λ1為常數(shù)。

利用Matlab軟件將式(3)和式(4)函數(shù)方程繪制成三維曲面，如圖4所示。由圖4可見：固支模型的頂板下沉撓度呈“尖細”狀，簡支模型的頂板下沉撓度呈“寬粗”狀，且中心下沉撓度范圍也增大。

圖4 頂板中心下沉撓度三維曲面Fig.4 3D surface of roof center subsidence deflection

階段Ⅲ：自由邊界模型。當頂板邊緣整體破壞時，四周支撐相當于自由邊，礦柱完全支撐著頂板及其上覆巖層的荷載，此時，頂板下沉撓度與其剛度無關，主要受礦柱的自身強度影響，此時，D=0，w=u，且有

頂板整體破壞時的條件為

式中：u2為第Ⅲ階段頂板中心下沉位移；λ2為常數(shù)。

根據(jù)損傷力學理論，礦柱應力應變關系可用Weibull分布描述，其應力與應變關系如下[14]：

式中：ε為礦柱應變；ε0為曲線σ-ε峰值點縱坐標的平均數(shù)；m為試驗擬合的均勻性指標。Weibull分布σ-ε峰值前的曲線與一元三次函數(shù)曲線具有較高的相似性，σ-ε峰值前曲線也可用下列表達式表示：

1.2.3 采場蠕變模型

礦柱變形具有隨時間變化的蠕變性，采用Burgers體模型來計算，如圖5所示，本構方程如下[11]：

式中：E1和E2為彈性系數(shù)；η1和η2為黏性系數(shù)。

圖5 礦柱Burgers物理本構模型Fig.5 Burgers physical constitutive model of pillar

聯(lián)立式(1)和式(11)，消去礦柱應力σ得

1.3 撓度曲線方程

將式(1)與式(2)聯(lián)立，通過伽遼金法求得

令

階段I：當頂板處于固支模型時，將式(3)代入式(14)求得

階段II：當頂板處于簡支模型時，將式(4)代入式(14)求得

階段Ⅲ：當頂板處于自由邊模型時，將式(6)代入式(14)求得

聯(lián)立式(2)與式(12)，采用伽遼金法求得

式(19)微分方程通解為

根據(jù)王金安等[13]的研究成果，礦柱在最初受壓時會產(chǎn)生瞬時變形，階段I中式(20)初始位移和下沉速度為

式中：σ0=w0ξE2，為礦柱初始應力。階段II和階段Ⅲ位移和速度的初始值分別為上一階段末期的對應值。

2 采空區(qū)群突變及失穩(wěn)機制

2.1 突變理論分析

突變理論用于研究系統(tǒng)從一種穩(wěn)定狀態(tài)向另一種穩(wěn)定狀態(tài)轉(zhuǎn)化的規(guī)律[15-16]，采空區(qū)群失穩(wěn)是一個漸進的、非線性過程，突變理論是計算這種失穩(wěn)過程的有效方法。對于復雜的內(nèi)部系統(tǒng)奇點附近不連續(xù)問題，可利用突變理論構建的數(shù)學模型，找到在系統(tǒng)臨界點的突變規(guī)律，能有效預測和判別巖體工程穩(wěn)定性。尖點突變的勢函數(shù)模型如下：

式中：V(x)為尖點突變勢函數(shù)；x為狀態(tài)變量；p和q為控制變量；這3個變量構成系統(tǒng)的相空間。

通過對勢函數(shù)求導可得相空間的平衡曲面方程：

對勢函數(shù)求二階倒數(shù)得系統(tǒng)的奇點集方程：

聯(lián)立式(23)和(24)得分叉集方程為

式(25)中，p為非正數(shù)是方程有實數(shù)解的前提條件，這也是突變產(chǎn)生的必要條件。

圖6 尖點突變模型Fig.6 Cusp catastrophe model

圖6所示為尖點突變模型的示意圖。相空間的圖形可以看作一個褶皺的曲面，該曲面的折疊圖形上分別由上葉、中葉和下葉構成。當平衡曲面位于上葉或下葉時，系統(tǒng)處于發(fā)展過程的準穩(wěn)定狀態(tài)；當平衡曲面位于中葉時，系統(tǒng)處于不穩(wěn)定狀態(tài)，即將發(fā)生突變，進入一個新的平衡狀態(tài)[17-18]。

2.2 勢函數(shù)構建

從能量轉(zhuǎn)化角度分析，采空區(qū)群勢函數(shù)表達式為

式中：U為采空區(qū)勢能；Ue和Us分別為頂板和礦柱的應變能；W為外力作用所做的功。按照薄板的彎曲理論，在固支、簡支和自由邊模型下，頂板應變能Ue表達式不盡相同。

階段I：當頂板處于固支模型時，應用薄板理論的頂板應變能表達式為

利用Matlab軟件解算出

階段Ⅱ：當頂板處于固支模型時，頂板應變能表達式為

式中：u為頂板中心下沉位移；μ為頂板等效剛度。

階段Ⅲ：當頂板處于自由邊模型時，頂板應變能為零。

礦柱壓縮應變能Us為

外力做的功W為

由式(22)～(25)得勢函數(shù)表達式：

對式(31)求偏導數(shù)，得

對式(32)整理成如下格式：

2.3 采空區(qū)群系統(tǒng)突變失穩(wěn)分析

2.3.1 突變失穩(wěn)計算

圖7所示為采空區(qū)群系統(tǒng)突變失穩(wěn)計算流程。由式(27)求出μ1，將其代入式(34)并聯(lián)立式(25)，可求出面積率臨界值ξ1，將其代入(33)求出頂板下沉撓度u1。若u1＜λ1，則突變點系統(tǒng)處于階段I，否則依次計算μ2，ξ2和u2；若u2＜λ2，則突變點系統(tǒng)處于階段Ⅱ，否則，處于階段Ⅲ。

圖7 采空區(qū)群系統(tǒng)突變失穩(wěn)計算流程圖Fig.7 Flow chart for calculating catastrophic instability of goaf group system

表1所示為雙向分層條帶式采空區(qū)巖體參數(shù)，采場長180 m，寬120 m，代入數(shù)據(jù)求得ξ1=0.471 3，u1=0.124 1 m，λ1=0.045 7 m，故u1>λ1，系統(tǒng)不處于階段Ⅰ；繼而求得ξ2=0.491 2，u2=0.130 6 m，λ2=0.088 6 m，故u2>λ2，系統(tǒng)不處于階段Ⅱ，而處于階段Ⅲ。

表2所示為不同采空區(qū)群系統(tǒng)突變失穩(wěn)計算結(jié)果，分別列出在固支模型、簡支模型下的頂板位移和破壞條件，表2中所列出的采空區(qū)群均會發(fā)生突變，采礦過程中要及時對采空區(qū)進行充填處理。

表1 采空區(qū)群頂板-礦柱巖體力學參數(shù)Table 1 Mechanical parameters of roof and pillar rock mass in goaf group

表2 不同采空區(qū)群系統(tǒng)突變失穩(wěn)計算結(jié)果Table 2 Calculation results of catastrophic instability of different goaf groups

2.3.2 數(shù)值模擬分析

采用數(shù)值模擬研究分析[20]，將3Dmine，MIDAS和FLAC3D這3種軟件耦合，充分利用軟件在圖形處理、復雜建模和定制計算等方面優(yōu)勢，分別研究采場空區(qū)群在不同邊界狀態(tài)下的穩(wěn)定性。建立的礦區(qū)及采場數(shù)值計算三維模型如圖8所示。

圖8 礦區(qū)及采場三維數(shù)值計算模型Fig.8 3D numerical calculation model of mining area and stope

圖9所示為采空區(qū)群下沉位移云圖，由圖9可見：采場頂板中心位移及其范圍隨采場參數(shù)減小而變小。圖9(a)中頂板中心最大位移區(qū)域占采場頂板區(qū)域的1/2，表現(xiàn)出明顯的不穩(wěn)定性，引發(fā)頂板大面積垮塌和突變失穩(wěn)；圖9(c)中頂板相對較穩(wěn)定；圖9(b)中頂板穩(wěn)定性居于圖9(a)和圖9(c)中頂板穩(wěn)定性之間。

圖10所示為從三維角度分析同一采空區(qū)群在階段Ⅰ和階段Ⅱ的頂板下沉位移。由圖10可見：階段Ⅰ和階段Ⅱ采場頂板下沉中心位移分別為0.10 m和0.12 m，即簡支模型下的頂板位移要大于固支模型的頂板位移，說明采場頂板邊緣發(fā)生塑性變形后，頂板整體剛度減小，進而導致頂板中心下沉撓度增加，這與理論計算結(jié)果是一致的。

2.4 采空區(qū)群系統(tǒng)蠕變失穩(wěn)分析

2.4.1 系統(tǒng)階段Ⅰ破壞時間

礦柱蠕變參數(shù)按照已有研究成果選取[19]，E1=8.0×105MPa，E2=5.5×104MPa，η1=1.05×104MPa·h，η2=1.98×109MPa·h。由式(21)可計算初始位移和下沉速度分別為w=0.000 6 m，w=v0=0.001 6 m/h。t=0時，代入式(20)求解C1和C2，根據(jù)這些初始值，確定微分方程系數(shù)：

根據(jù)式(5)破壞條件，求解上述線性方程，得系統(tǒng)在階段Ⅰ破壞時間為688 h。

2.4.2 系統(tǒng)階段Ⅱ破壞時間

階段Ⅱ的頂板初始位移和速度為階段Ⅰ結(jié)束時的位移和速度，分別為w=0.045 7 m，w=6.552×10-5m/h，代入式(33)得t=0時，微分方程系數(shù)為

圖9 采空區(qū)群頂板下沉位移Fig.9 Roof subsidence displacement of goaf group

根據(jù)式(7)破壞條件，求解上述線性方程，得系統(tǒng)在階段Ⅱ破壞時間為693 h。系統(tǒng)在進入階段Ⅱ之后，頂板出現(xiàn)破裂，雖然沒有整體垮塌，頂板已經(jīng)處于破裂狀態(tài)，認為已失去穩(wěn)定性，采空區(qū)穩(wěn)定預測時間按照階段Ⅰ和階段Ⅱ破壞時間之和計算，因此，頂板-礦柱系統(tǒng)破壞時間為1 381 h，約為2月。

3 采空區(qū)群系統(tǒng)失穩(wěn)因素分析

3.1 影響因子對頂板下沉位移影響

頂板下沉位移為采空區(qū)群系統(tǒng)穩(wěn)定性的關鍵指標，從式(33)～(34)可以看出影響采空區(qū)群穩(wěn)定性的主要因子有：采場參數(shù)a和b，礦柱支撐面積率ξ，均布荷載q0，礦柱峰值抗壓強度σm和礦柱彈性模量E0。下面選取幾組不同采場參數(shù)，分別分析各影響因子對頂板下沉位移的影響。分別將式(15)和(16)代入式(33)和(34)，可得到以ξ，q0，H，σm和E0為自變量，以u為因變量的隱函數(shù)。階段Ⅰ和Ⅱ在蠕變突變發(fā)生機制方面具有代表性、相關性和可比性，因此，分析采空區(qū)群系統(tǒng)失穩(wěn)因素，主要針對這2個階段進行對比。當其他自變量值一定時，繪制單個自變量對因變量u的影響曲線，如圖11和圖12所示。

圖11(a)和圖12(a)可見：對于參數(shù)確定的采場，隨著ξ減小，頂板位移逐漸增大，直至曲線的拐點處(臨界值)產(chǎn)生突變失穩(wěn)；采場參數(shù)越大，保持系統(tǒng)穩(wěn)定所需的ξ越大；系統(tǒng)由階段Ⅰ至變化至階段Ⅱ時，保持系統(tǒng)穩(wěn)定所需ξ均變大，采場參數(shù)越小時變化越明顯；在階段Ⅱ，不同采場ξ的差值在縮小。從圖11(b)和圖12(b)可見：采場參數(shù)越大，能承受的q0越小。此外，階段Ⅰ與階段Ⅱq0變化幅度與采場大小呈反比。在圖11(a)和(b)以及圖12(a)和(b)中，曲線拐點處以下部分為有效研究內(nèi)容，拐點處以上部分無意義，研究時不予考慮。

從圖11(c)和圖12(c)可見：頂板位移隨礦柱峰值抗壓強度σm增大而呈“L”型曲線減小，在接近σm的臨界值時，頂板位移變化幅度最大；當σm超過某一特定值后，頂板位移接近恒定值。從圖11(d)和圖12(d)可見：頂板位移隨礦柱彈性模量E0增大而呈“C”型曲線減小。

從圖11和圖12中可以總結(jié)出如下規(guī)律：頂板位移與礦柱支撐面積率、礦柱峰值抗壓強度和礦柱彈性模量呈負相關，與均布荷載呈正相關；當頂板-礦柱系統(tǒng)處在固支模式時，采場參數(shù)對系統(tǒng)穩(wěn)定性影響顯著，隨著采場四周出現(xiàn)塑性變形，采場參數(shù)對系統(tǒng)的穩(wěn)定性差異逐漸縮小。

圖13所示為不同采場礦柱穩(wěn)定支撐的臨界值。從圖13可見：當采場參數(shù)a＞90 m，b＞60 m時，階段Ⅰ與階段Ⅱ中礦柱面積支撐率ξ趨于相近，這說明當采場參數(shù)較大時，固支模型和簡支模型對采場頂板位移影響不明顯。

3.2 影響因子對臨界點的影響

令f=4p3+27q2，與式(34)相結(jié)合，可組成分別以ξ，q0，H，σm和E0為自變量，f(ξ)，f(q0)，f(σm)和f(E0)為因變量的函數(shù)，將其稱之為分叉集函數(shù)。只有當f<0時，系統(tǒng)才發(fā)生突變；當f＞0，系統(tǒng)沒有穩(wěn)定可能性。在礦山生產(chǎn)開采中，可通過調(diào)節(jié)影響因子參數(shù)維持系統(tǒng)的穩(wěn)定性。

圖12 階段Ⅱ影響因子對頂板下沉位移的影響Fig.12 Influence of various factors on roof sinking displacement in stageⅡ

圖13 不同采場礦柱穩(wěn)定支撐率的臨界值Fig.13 Critical value of stability support ratio of pillars in different stopes

圖14所示為階段Ⅰ影響因子對臨界點的影響情況，圖15所示為階段Ⅱ影響因子對臨界點的影響情況。圖14(a)和圖15(a)顯示礦柱支撐面積率ξ對突變判別函數(shù)影響。在階段Ⅰ與階段Ⅱ中，不同采場結(jié)構參數(shù)下f(ξ)差距不大。圖15(a)顯示，當采場參數(shù)a=90 m，b=60 m時，系統(tǒng)在進入階段Ⅱ中f(ξ)＞0，即在階段Ⅱ中沒有穩(wěn)定可能性。圖14(b)和圖15(b)中顯示均布荷載q0對突變判別函數(shù)的影響特征，采場幾何尺寸越大，q0的臨界值越小。此外，階段Ⅰ與階段Ⅱq0的臨界值變化幅度與采場幾何尺寸呈反比。由圖14(c)和(d)以及圖15(c)和(d)可見：突變判別函數(shù)f(σm)隨礦柱峰值抗壓強度σm增大逐漸遠離臨界點，而f(E0)隨礦柱彈性模量E0的增大逐漸接近恒定值。

從圖14和圖15總結(jié)如下規(guī)律：分叉集函數(shù)f與礦柱支撐面積率和礦柱峰值抗壓強度呈負相關，與均布荷載和礦柱彈性模量呈正相關；當f<0時，f越接近零越容易發(fā)生突變。在分叉集函數(shù)f曲線中，當?shù)V柱支撐率和峰值抗壓強度大于臨界值時，系統(tǒng)將有可能發(fā)生突變；當其小于臨界值時，系統(tǒng)在初始狀態(tài)就不夠穩(wěn)定，均布荷載與之相反。系統(tǒng)在階段Ⅱ時，不同采場參數(shù)間均布荷載對臨界點的影響趨于一致。當柱彈性模量在超過臨界值后，f保持恒定，此時將不會發(fā)生突變。

圖14 階段Ⅰ影響因子對臨界點的影響Fig.14 Influence of various factors on critical points in stageⅠ

圖15 階段Ⅱ影響因子對臨界點的影響Fig.15 Influence of various factors on critical points in stageⅡ

4 結(jié)論

1)將采空區(qū)勢函數(shù)與尖點突變理論相結(jié)合，通過實例計算和數(shù)值模擬分析，得出多種采場參數(shù)在不同支撐邊界模式下的變形特性和力學響應，從而判別出采空區(qū)群突變失穩(wěn)狀態(tài)，計算出采空區(qū)失穩(wěn)定時間，為優(yōu)化采場參數(shù)、尋找最佳充填方案、穩(wěn)定性控制和安全隱患管理等方面提供依據(jù)。

2)頂板位移與礦柱支撐面積率、礦柱峰值抗壓強度和礦柱彈性模量呈負相關，與均布荷載呈正相關；分叉集函數(shù)f與礦柱支撐面積率、礦柱峰值抗壓強度呈負相關。與均布荷載和礦柱彈性模量呈正相關，通過調(diào)控影響因子參數(shù)，可以改變采空區(qū)的穩(wěn)定狀態(tài)。

3)當頂板-礦柱系統(tǒng)處在固支模式時，采場參數(shù)對系統(tǒng)穩(wěn)定性影響最顯著；隨著采場四周出現(xiàn)塑性變形，不同采場參數(shù)對系統(tǒng)的穩(wěn)定性差異在逐漸縮小；當采場參數(shù)較大時，固支模型和簡支模型對系統(tǒng)穩(wěn)定性影響差距不明顯。