基于突變理論的采空區(qū)頂板安全厚度多因素預(yù)測模型

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

(1.北京科技大學(xué)土木與資源工程學(xué)院，北京 100083；2.北京科技大學(xué)城市地下空間工程北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083)

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

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

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

1 基于尖點(diǎn)突變理論的強(qiáng)度折減法

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

c′=c/K

(1)

φ′=φ/K

(2)

釆空區(qū)的失穩(wěn)通常是瞬間完成的，可以認(rèn)為是一個突變過程。因此采用強(qiáng)度折減法分析采空區(qū)穩(wěn)定性時，采空區(qū)的失穩(wěn)可以通過相應(yīng)部位的位移突變來反映。將不同折減系數(shù)K下的采空區(qū)頂板中點(diǎn)豎向位移δ在坐標(biāo)系中描出，再用4次項(xiàng)泰勒級數(shù)函數(shù)形式進(jìn)行擬合[7]，計算過程見式(3)。

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

(3)

式中：a0、a1、a2、a3、a4為待定系數(shù)。

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

(4)

Q=p4+up2+vp+C

(5)

式中：u和v為控制變量；C為常數(shù)項(xiàng)，對突變無影響。

對式(5)進(jìn)行求導(dǎo)，并求判別式，得式(6)。

Δ=8u3+27v2

(6)

根據(jù)突變理論，可知采空區(qū)頂板穩(wěn)定性判據(jù)為：Δ>0，采空區(qū)頂板處于穩(wěn)定狀態(tài)；Δ<0，采空區(qū)頂板處于失穩(wěn)狀態(tài)；Δ=0，采空區(qū)頂板處于臨界狀態(tài)。

2 采空區(qū)頂板安全系數(shù)分析

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

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

σx=1.934+0.047 8D

σy=0.409+0.010 4D

σz=0.485+0.029 2D

(7)

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

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

模型初始平衡后，通過不斷的折減頂板強(qiáng)度參數(shù)c和φ，得到不同折減系數(shù)K下頂板中點(diǎn)的位移值δ(K)，δ(K)與K的擬合曲線及擬合方程如圖2所示。將擬合曲線中各參數(shù)代入式(5)計算得到控制變量u和v，將u和v代入式(6)可得突變特征值Δ，進(jìn)而可得Δ隨K的變化曲線(圖3)。由圖3可知，當(dāng)K由1.65變?yōu)?.70時，Δ由正值變?yōu)樨?fù)值，頂板由穩(wěn)定變?yōu)槭Х€(wěn)，可以認(rèn)為臨界折減系數(shù)，即頂板安全系數(shù)在1.65～1.70之間。

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

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

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

圖4 頂板厚度與頂板安全系數(shù)關(guān)系曲線Fig.4 The relationship between the thickness ofthe roof and the safety factor of the roof

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

K=1.08+0.16ln(H-8.56)

(8)

將頂板安全系數(shù)K=1.5作為評價頂板穩(wěn)定性的判斷標(biāo)準(zhǔn)[7]，即：①K≥1.5時，采空區(qū)頂板穩(wěn)定；②K<1.5時，采空區(qū)頂板失穩(wěn)。

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

3 頂板安全厚度單因素試驗(yàn)

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

表2 模擬方案及結(jié)果Table 2 Simulation scheme and results

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

H=23.58+1.165d

R2=0.988 9

(9)

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

H=80.42-3.74L+0.05L2

R2=0.993 2

(10)

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

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

R2=0.999 9

(11)

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

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

H=22.19+0.103t-5.4

R2=0.992 6

(12)

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

(13)

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

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

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

H=23.232+0.104t-5.4

R2=0.477

(14)

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

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

R2=0.909

(15)

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

R2=0.96

(16)

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

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

R2=0.995

(17)

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

5 工程實(shí)例

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

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

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

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

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

6 結(jié) 論

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

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

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