基于突變理論的采空區頂板安全厚度多因素預測模型

尚向凡，苗勝軍，于兆新，羅立民

(1.北京科技大學土木與資源工程學院，北京 100083；2.北京科技大學城市地下空間工程北京市重點實驗室，北京 100083)

作為采空區的相對薄弱部分，頂板的穩定性一直是影響礦山地下安全開采的主要問題之一，因此確定采空區頂板安全厚度至關重要。采空區頂板安全厚度會隨著采空區的跨度、高度及頂板巖體強度的變化而變化。對于采空區頂板安全厚度的確定，傳統的荷載傳遞線交匯法、厚跨比法、魯佩涅依特理論計算法、平板梁理論計算法等方法考慮的因素比較單一，對頂板的受力機制與破壞過程認識不足，應用范圍也受到一定的限制[1-3]。近年來，一些學者利用數值模擬得方法進行采空區頂板穩定性分析，該方法方便靈活、耗資少并且能進行反復試驗，可以廣泛應用于地下工程穩定性分析[4-6]。

采空區頂板失穩是一個不連續、不可逆的突變過程，具有突發性，因此一些學者將突變理論引入采空區頂板穩定性分析。趙延林等[7]提出采空區重疊頂板安全系數的概念，建立了頂板豎向位移序列與折減系數的尖點突變模型，以此作為采空區重疊頂板是否失穩的判據。張欽禮等[8]運用尖點突變理論分析采場突變過程中的能量釋放機理，構建采場破壞失穩的尖點突變模型，推導得出系統失穩的充要條件。徐恒等[9]根據頂板尖點突變模型失穩的充要條件，計算得到充填體下采空區頂板安全厚度的表達式。杜逢彬等[10]以我國東南某特大金銅礦山為背景，應用基于尖點突變的強度折減法建立了隔離頂柱安全系數與其厚度之間的函數關系。

本文運用有限差分軟件FLAC3D，應用基于尖點突變理論的強度折減法對西石門鐵礦采空區穩定性進行研究，探討了采空區縱深、跨度與高度、頂板黏聚力與抗拉強度對采空區頂板安全厚度的影響，并建立了綜合考慮采空區縱深與跨度、頂板黏聚力與抗拉強度的頂板安全厚度預測模型，為確定采空區頂板安全厚度提供了一個新的定量研究方法。

1 基于尖點突變理論的強度折減法

有限元強度折減法[11]是將頂板的黏聚力c與內摩擦角φ同時除以一個折減系數K，將折減后的黏聚力c′(式(1))和內摩擦角φ′(式(2))帶入模型進行計算，逐步增大K直至模型達到極限狀態發生破壞，將此時的折減系數K視為頂板的安全系數。

c′=c/K

(1)

φ′=φ/K

(2)

釆空區的失穩通常是瞬間完成的，可以認為是一個突變過程。因此采用強度折減法分析采空區穩定性時，采空區的失穩可以通過相應部位的位移突變來反映。將不同折減系數K下的采空區頂板中點豎向位移δ在坐標系中描出，再用4次項泰勒級數函數形式進行擬合[7]，計算過程見式(3)。

δ(K)=a0+a1K+a2K2+a3K3+a4K4

(3)

式中：a0、a1、a2、a3、a4為待定系數。

V(p)=c0+c1p+c2p2+c4p4

(4)

Q=p4+up2+vp+C

(5)

式中：u和v為控制變量；C為常數項，對突變無影響。

對式(5)進行求導，并求判別式，得式(6)。

Δ=8u3+27v2

(6)

根據突變理論，可知采空區頂板穩定性判據為：Δ>0，采空區頂板處于穩定狀態；Δ<0，采空區頂板處于失穩狀態；Δ=0，采空區頂板處于臨界狀態。

2 采空區頂板安全系數分析

根據西石門鐵礦中采區首采中段勘探資料，取采空區跨度L為55 m、高度h為45 m、縱深d為40 m、頂板厚度H為50 m建立計算模型(圖1)，通過基于尖點突變理論的強度折減法確定采空區的頂板安全厚度。為保證基礎數據的準確性，對室內試驗測得的巖石物理力學參數進行經驗折減(表1)。

根據西石門鐵礦的地應力資料[12]，結合“構造應力+重力”作用進行初始應力場計算，得到最大水平主應力、最小水平主應力和垂直主應力隨深度變化的回歸方程見式(7)。

σx=1.934+0.047 8D

σy=0.409+0.010 4D

σz=0.485+0.029 2D

(7)

式中：σx為水平x方向主應力，MPa；σy為水平y方向主應力，MPa；σz為垂直z方向主應力，MPa；D為埋深，m。

圖1 模型示意圖Fig.1 Model diagram

模型初始平衡后，通過不斷的折減頂板強度參數c和φ，得到不同折減系數K下頂板中點的位移值δ(K)，δ(K)與K的擬合曲線及擬合方程如圖2所示。將擬合曲線中各參數代入式(5)計算得到控制變量u和v，將u和v代入式(6)可得突變特征值Δ，進而可得Δ隨K的變化曲線(圖3)。由圖3可知，當K由1.65變為1.70時，Δ由正值變為負值，頂板由穩定變為失穩，可以認為臨界折減系數，即頂板安全系數在1.65～1.70之間。

為確定頂板安全系數的具體數值，提高模擬計算精度得到K分別為1.66、1.67、1.68、1.69時的頂板中點豎向位移，按照上述步驟擬合曲線并求突變特征值，最終判定西石門鐵礦采空區跨度為55 m、高度為45 m、縱深為40 m、頂板厚度為50 m時的頂板安全系數為1.67。

圖2 折減系數與頂板豎向位移擬合曲線Fig.2 Fitting curve between reduction coefficient androof vertical displacement

為研究頂板厚度對安全系數的影響，保持采空區大小與模型參數不變，分別建立頂板厚度為40 m、30 m、20 m、10 m的模型進行數值模擬，得到相應的安全系數分別為1.62、1.55、1.46、1.14，頂板安全系數與頂板厚度的關系曲線見圖4。

圖4 頂板厚度與頂板安全系數關系曲線Fig.4 The relationship between the thickness ofthe roof and the safety factor of the roof

采用對數函數對圖4曲線進行擬合，可得頂板安全系數與頂板厚度的關系式為式(8)。

K=1.08+0.16ln(H-8.56)

(8)

將頂板安全系數K=1.5作為評價頂板穩定性的判斷標準[7]，即：①K≥1.5時，采空區頂板穩定；②K<1.5時，采空區頂板失穩。

在式(8)中令K=1.5，求出H=22.36 m即為西石門鐵礦中采區首采中段的頂板安全厚度。

3 頂板安全厚度單因素試驗

通過上述計算得到了西石門鐵礦中采區首采中段的頂板安全厚度，為研究采空區縱深d、跨度L、高度h、頂板巖體黏聚力c、抗拉強度t對西石門鐵礦采空區頂板安全厚度的影響，設計了27種不同影響因素組合下的模擬方案(表2)，計算所得各因素與頂板安全厚度的相關關系曲線如圖5所示。

表2 模擬方案及結果Table 2 Simulation scheme and results

1) 開挖縱深。由試驗4和試驗26、試驗5和試驗27可知，采空區縱深d對頂板安全厚度H的影響受跨度比值d/L(L為采空區跨度)的限制。當d/L小于2時，d對H影響較大；當d/L大于2時，d對H影響相對較小。這也是采空區不斷向前延伸(d/L遠大于2)的原因，但只要跨度沒有變化，采空區頂板就可以保持穩定的原因。當d/L小于2時，縱深d與頂板安全厚度H的關系如圖5(a)所示，H隨著d的增大均勻增大，采用Origin對圖中曲線進行擬合可得H與d的關系式(式(9))。在實際工程計算中，若d/L大于2，可將d設為L的兩倍求出頂板安全厚度的近似值。

H=23.58+1.165d

R2=0.988 9

(9)

2) 采空區跨度。由圖5(b)可知，跨度L與頂板安全厚度H之間呈現非線性正相關的變化關系。對曲線進行擬合可得L與H的關系式，見式(10)。

H=80.42-3.74L+0.05L2

R2=0.993 2

(10)

3) 采空區高度。由圖5(c)可知，高度h對頂板安全厚度H的影響非常小，隨著h的增大，H幾乎不變。故在后續的預測模型中，h將不予考慮。

4) 頂板黏聚力。由圖5(d)可知，黏聚力c與頂板安全厚度H之間呈現非線性負相關的變化關系。對曲線進行擬合可得c與H的關系式，見式(11)。

R2=0.999 9

(11)

5) 頂板抗拉強度。由圖5(e)可知，抗拉強度t與頂板安全厚度H之間呈現非線性負相關的變化關系，隨著抗拉強度的減小，頂板安全厚度不斷增大。由試驗16和試驗17可知，當t小于一定值時，H會迅速增大。說明t存在一個臨界值，小于這個臨界值時，t的變化會對H造成很大的影響。而當t大于1 MPa時，隨著t的增大，H的變化幅度很小，說明此時的t已經大于模型最大拉應力，增大t對H的影響不大。

對曲線進行擬合可得t與H的關系式，見式(12)。

H=22.19+0.103t-5.4

R2=0.992 6

(12)

4 多因素頂板安全厚度線性預測模型

(13)

式中，α0、α1、α2、α3、α4均為待定系數。

采用多元回歸法求解待定系數，各因素按照顯著性的大小步進法輸入，具體求解過程如下所述。

1) 考慮t-5.4進行回歸，見式(14)。

H=23.232+0.104t-5.4

R2=0.477

(14)

2) 在式(14)的基礎上考慮d進行回歸，見式(15)。

H=-39.358+0.118t-5.4+1.389d

R2=0.909

(15)

圖5 各因素與頂板安全厚度的關系曲線Fig.5 Relationship between various factors and the thickness of the safety roof

R2=0.96

(16)

4) 在式(16)的基礎上考慮L2進行回歸，見式(17)。

H=-81.194+0.104t-5.4+1.158d+

R2=0.995

(17)

隨著擬合因素的增加，相關系數R2的數值逐漸增大，當4個因素全部考慮后，R2的數值達到了0.995，得到的預測公式可以為類似地質下鐵礦開采過程中的頂板厚度設計提供參考依據，該預測模型方法可為礦山頂板厚度設計提供思路。將表2中前21組模擬中所用采空區縱深、跨度、頂板巖體黏聚力、抗拉強度數據代入式(17)進行計算，所得結果與表2中模擬結果相差不大，相對誤差分布在0.3%～16.5%之間。

5 工程實例

為驗證預測模型的有效性，對與西石門鐵礦地質情況相似的某硫鐵礦采區瞬變電磁勘探結果進行分析。該硫鐵礦采區1 100水平、1 070水平、1 050水平、1 000水平分別呈低阻、中阻、中高阻和高阻，各水平視電阻率剖面疊加圖如圖6所示。

1 100水平視電阻率剖面圖中低阻異常區的位置和其余三張圖中高阻位置基本一致，范圍略有不同(圖6)，推測在1 050水平已進入采區空洞中，視電阻率急劇增大，表明頂板邊界位置大部位于1 050水平～1 070水平之間。

根據現場實踐情況，將采空區的尺寸簡化為跨度45 m、縱深80 m的矩形區域，對應于圖6中1 000水平視電阻率剖面圖的矩形方框部分。在考慮尺寸效應及地層結構面的影響后，巖體的黏聚力c和抗拉強度t通過經驗折減修正為c=11.40 MPa、t=4.37 MPa，將其帶入式(17)計算可得該硫鐵礦頂板安全厚度H=41.87 m。

根據瞬變電磁勘探數據，采空區上覆頂層的最小厚度為32 m。勘探結束一個月后，采空區頂板突然發生塌陷，因此可以認為勘探時頂板處于臨界狀態(即Δ=0)，對應的探測數據32 m為采空區的頂板安全厚度(即K=1.5)。勘探數據與預測模型所得值相差9.87 m，相對誤差為30.84%，說明該數學預測模型方法的有效性，其他礦山可根據自身巖石力學參數推導出相應的預測模型。

圖6 各水平視電阻率剖面疊加圖Fig.6 Overlay of each horizontal apparent resistivity profile

6 結 論

1) 應用基于尖點突變理論的強度折減法對西石門鐵礦采空區穩定性進行研究，確定了其頂板安全厚度，建立了綜合考慮采空區縱深與跨度、頂板黏聚力與抗拉強度的頂板安全厚度預測模型，為確定采空區頂板安全厚度提供了一個新的定量研究方法。

2) 西石門鐵礦頂板安全厚度H與各因素的關系為：d對H的影響受跨度比值d/L的限制，當d/L小于2時，d對H影響較大；當d/L大于2時，d對H影響較小；L與H之間呈現非線性正相關的變化關系；h對H的影響非常小；c與H之間呈現非線性負相關的變化關系；t存在一個臨界值，t小于這個臨界值時，t的變化會對H造成很大的影響。

3) 對與西石門鐵礦地質條件相似的某硫鐵礦采區瞬變電磁勘探結果進行分析，探測厚度與預測模型所得頂板安全厚度相對誤差為30.84%，驗證了預測模型方法的有效性。