任紅崗，汪旭光，譚卓英，夏志遠(yuǎn)

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

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

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

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

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

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

1.2 力學(xué)模型構(gòu)建

1.2.1 采場(chǎng)力學(xué)模型

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

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

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

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

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

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

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

1.2.2 頂板撓度算法

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

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

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

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

頂板邊緣破壞條件為

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

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

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

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

頂板整體破壞時(shí)的條件為

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

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

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

1.2.3 采場(chǎng)蠕變模型

礦柱變形具有隨時(shí)間變化的蠕變性，采用Burgers體模型來(lái)計(jì)算，如圖5所示，本構(gòu)方程如下[11]：

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

圖5 礦柱Burgers物理本構(gòu)模型Fig.5 Burgers physical constitutive model of pillar

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

1.3 撓度曲線方程

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

令

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

階段II：當(dāng)頂板處于簡(jiǎn)支模型時(shí)，將式(4)代入式(14)求得

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

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

式(19)微分方程通解為

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

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

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

2.1 突變理論分析

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

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

通過(guò)對(duì)勢(shì)函數(shù)求導(dǎo)可得相空間的平衡曲面方程：

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

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

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

圖6 尖點(diǎn)突變模型Fig.6 Cusp catastrophe model

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

2.2 勢(shì)函數(shù)構(gòu)建

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

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

階段I：當(dāng)頂板處于固支模型時(shí)，應(yīng)用薄板理論的頂板應(yīng)變能表達(dá)式為

利用Matlab軟件解算出

階段Ⅱ：當(dāng)頂板處于固支模型時(shí)，頂板應(yīng)變能表達(dá)式為

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

階段Ⅲ：當(dāng)頂板處于自由邊模型時(shí)，頂板應(yīng)變能為零。

礦柱壓縮應(yīng)變能Us為

外力做的功W為

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

對(duì)式(31)求偏導(dǎo)數(shù)，得

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

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

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

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

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

表1所示為雙向分層條帶式采空區(qū)巖體參數(shù)，采場(chǎng)長(zhǎng)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)計(jì)算結(jié)果，分別列出在固支模型、簡(jiǎn)支模型下的頂板位移和破壞條件，表2中所列出的采空區(qū)群均會(huì)發(fā)生突變，采礦過(guò)程中要及時(shí)對(duì)采空區(qū)進(jìn)行充填處理。

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

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

2.3.2 數(shù)值模擬分析

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 影響因子對(duì)頂板下沉位移影響

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

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

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

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

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

3.2 影響因子對(duì)臨界點(diǎn)的影響

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

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

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

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

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

圖14 階段Ⅰ影響因子對(duì)臨界點(diǎn)的影響Fig.14 Influence of various factors on critical points in stageⅠ

圖15 階段Ⅱ影響因子對(duì)臨界點(diǎn)的影響Fig.15 Influence of various factors on critical points in stageⅡ

4 結(jié)論

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

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

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