采空區(qū)穩(wěn)定性數(shù)值模擬研究

王永奇,戴 兵

(1.貴州開磷集團礦業(yè)總公司， 貴州 貴陽 550302；2.中南大學資源與安全工程學院, 湖南 長沙 410083)

0 引 言

礦產(chǎn)資源安全、高效開采對國民經(jīng)濟的發(fā)展至關重要。在上世紀中后期，房柱法、留礦法、階段礦房法等空場法[1]在地下礦山中廣泛運用，取得良好開采效果。受經(jīng)濟因素影響和空區(qū)安全隱患意識淡薄，部分采用空場法開采后殘留的空區(qū)未進行處理。由于部分采空區(qū)高度大、跨度大，且空區(qū)在受到長時間的風化、爆破振動的作用下，空區(qū)周圍巖體強度大幅降低，經(jīng)常出現(xiàn)采空區(qū)頂板坍塌、冒落現(xiàn)象，嚴重影響礦山后續(xù)開采中人員、設備的安全。

采空區(qū)穩(wěn)定性制約礦山的發(fā)展，因此對其進行研究具有重要意義。目前，對采空區(qū)穩(wěn)定性研究主要通過力學理論計算[2]、數(shù)學模型[3,4]和數(shù)值計算。由于巖體中存在大量節(jié)理、裂隙，通過理論計算所得結(jié)果往往存在較大誤差，且理論計算相對復雜，難以得到準確的計算結(jié)果。數(shù)學模型通過數(shù)學方法對采空區(qū)進行分類評價，雖能獲得分類結(jié)果，但不能直觀體現(xiàn)空區(qū)頂板受力狀態(tài)，且分類結(jié)果與實際存在一定誤差。而數(shù)值模擬計算能根據(jù)空區(qū)圍巖實際情況對力學參數(shù)進行調(diào)整，模擬精度相對較高。因此本文以開磷礦區(qū)內(nèi)空區(qū)為研究對象，擬通過數(shù)值模擬技術對采空區(qū)穩(wěn)定性進行研究。

1 礦山概況

開陽磷礦自1958年建礦以來，已連續(xù)生產(chǎn)了50多年。以往礦山使用的采礦方法為空場法，采后采空區(qū)未進行處理。礦區(qū)形成的采空區(qū)，導致山體塌陷和滑坡。礦體延伸15 km，從地表可見綿延15 km的山體塌陷區(qū)。礦體走向南北，傾向為280°～290°，傾角27°～30°。受F314斷層影響，斷層下盤礦體產(chǎn)狀由北向南逐漸變陡，在W11～W11+1線礦體出現(xiàn)倒轉(zhuǎn)。礦體厚度5.8～7.5 m，南薄北厚，厚度變化不大，平均6.0 m左右。

2 ANSYS數(shù)值軟件

ANSYS數(shù)值軟件采用有限元計算程序，集結(jié)構(gòu)、熱、流體、電磁場、聲場和耦合場分析于一體，能處理線性和非線性復雜力學計算問題，在開發(fā)之初便得到廣泛應用[5]。

有限單元法基本思路是將連續(xù)結(jié)構(gòu)離散成單元，單元上設定節(jié)點，使用變分原理建立有限元方程，求解離散域中的自由度問題。有限元分析求解過程主要可分為以下幾個步驟[6]。

(1)對結(jié)構(gòu)進行離散化。通過將求解的連續(xù)結(jié)構(gòu)化解與分割，轉(zhuǎn)換成有限單元和節(jié)點，保證相鄰單元力學性能參數(shù)連續(xù)。

(2)位移插值函數(shù)選擇。通過對單元位移函數(shù)的假設，利用節(jié)點位移來表示單元體位移、形變和應力，以對連續(xù)體問題進行分析。位移矩陣如下：

{f}=[N]{δ}e

(1)

式中，{f}—單元中任意一點的位移，

[N]—行函數(shù)，

{δ}—單元節(jié)點位移。

(3)單元力學特性。首先推導出單元應變(利用節(jié)點位移)表達式：

{ε}=[B]{δ}e

(2)

式中，{ε}—單元應變，

[B]—單元應變矩陣。

根據(jù)本構(gòu)方程，得出單元應力：

{σ}=[D][B]{δ}e

(3)

式中，[D]—彈性矩陣(與單元材料相關)。

最后根據(jù)變分原理，推導出單元節(jié)點力關于節(jié)點位移的關系式：

[F]=[k]e{δ}e

(4)

式中，[k]e—單元剛度矩陣：

[k]e=?[B]T[D][B]dxdydz

(5)

(4)根據(jù)(1)、(2)、(3)、(4)、(5)平衡方程，建立整體連續(xù)結(jié)構(gòu)的平衡方程：

[K]{σ}=[F]

(6)

式中：[K]為總剛度矩陣。

(5) 計算求解節(jié)點位移和單元應力。有限單元法求解程序的內(nèi)部過程如圖1所示。

圖1 有限單元法求解程序的內(nèi)部過程

3 采空區(qū)穩(wěn)定性的數(shù)值模擬

本文主要對開磷礦區(qū)3#采空區(qū)和6#采空區(qū)穩(wěn)定性進行數(shù)值模擬，3#空區(qū)尺寸為120 m×50 m×17.5 m；6#采空區(qū)主要由A、B、C三個小采空區(qū)構(gòu)成，三個空區(qū)尺寸分別為85 m×72 m×30 m,50 m×30 m×15 m、30 m×30 m×30 m，空區(qū)埋深均為230 m。

3.1 基本假定

假定礦體為理想彈塑性體，礦體和圍巖為局部均質(zhì)、各向同性的材料。考慮到采空區(qū)跨度極大，計算按平面應變問題進行。由于巖石具有脆性，分析中涉及到的所有物理量均與時間無關。

3.2 巖石物理力學參數(shù)

頂板巖性為白云巖，有輕微破碎，節(jié)理不太發(fā)育。頂板巖石密度為2.7 t/m3，抗拉強度為5.5 MPa; 單軸抗壓強度為159 MPa，彈性模量為30，泊松比為0.27， 內(nèi)聚力(抗剪強度)為37.49 MPa，內(nèi)磨擦角為32.27°，巖石較穩(wěn)定；礦石巖性為褐色磷礦石，較為致密堅硬。密度為3.3 t/m3， 抗拉強度為4.5 MPa，單軸抗壓強度為148 MPa，彈性模量為29，泊松比為0.25， 內(nèi)聚力(抗剪強度)為36.67 MPa，內(nèi)磨擦角為41.94°，巖石較穩(wěn)定；底板砂巖呈青灰色，有輕微破碎，節(jié)理不發(fā)育，其密度為2.7 t/m3，抗拉強度為3 MPa，單軸抗壓強度為110 MPa，彈性模量為18，泊松比為0.23，內(nèi)聚力(抗剪強度)為29.78 MPa，內(nèi)磨擦角為42.56°，巖石較穩(wěn)定。底板頁巖呈醬紅色，層理明顯、節(jié)理發(fā)育，巖石易風化，遇水易泥化，其密度為2.7 t/m3，抗拉強度為2.7 MPa，單軸抗壓強度為1 MPa，彈性模量為9.5，泊松比為0.39，內(nèi)聚力(抗剪強度)為14.09 MPa，內(nèi)磨擦角為42.56°，巖石不穩(wěn)定。

3.3 建立計算模型

采用ANSYS前處理分別對3#、6#采空區(qū)建模。為盡可能反映采空區(qū)對圍巖應力場分布的影響，所建模型的尺寸大小為5倍采空區(qū)尺寸。模型建立完成后進行網(wǎng)格劃分，3#空區(qū)和6#空區(qū)模型劃分網(wǎng)格后分別如圖2、圖3所示。

圖2 3#采空區(qū)模型單元

圖3 6#采空區(qū)模型單元

3.4 邊界條件與荷載

施加模型邊界條件主要考慮構(gòu)造應力的作用。首先對3#、6#采空區(qū)模型X軸方向的前后兩個面進行法向約束，然后再對Y軸方向的前后兩個面進行約束，最后對兩個采空區(qū)模型的底面進行約束。

本模型主要分析空區(qū)圍巖的穩(wěn)定性，因此主要考慮在重力場的作用下圍巖的應力應變特征。

原巖自重應力可分別按下式進行計算：

(7)

(8)

式中：γi為上覆第層巖體重度，kN/m3；為上覆巖體泊松比；Hi為上覆巖體分層厚度，m。

3.5 模擬結(jié)果分析

3#空區(qū)和6#空區(qū)X、Y、Z向的應力分布，分別如圖4、圖5所示。

圖4 3#采空區(qū)等效應力圖

模擬結(jié)果表明，3#采空區(qū)的頂板主要受到35.5 MPa的壓應力，并在空區(qū)中央出現(xiàn)拉應力集中現(xiàn)象，最大拉應力為13.6 MPa，遠大于頂板巖體抗拉強度，空區(qū)中央巖體發(fā)生拉伸破壞。從空區(qū)中心向空區(qū)兩端延伸，拉應力逐漸減小，直至變成壓應力，并在空區(qū)頂板與側(cè)幫相交位置出現(xiàn)壓應力集中，最大壓應力為84.7 MPa。空區(qū)底板與側(cè)幫相交位置拉應力作用較為明顯，最大拉應力為16.6 MPa，遠大于底板巖體的抗拉強度，發(fā)生拉伸破壞。對于采空區(qū)的凸出部位主要表現(xiàn)為拉應力，最大拉應力為11.6 MPa，大于巖體抗拉強度，也發(fā)生拉伸破壞。而采空區(qū)的凹進部位壓應力作用較為明顯，最大壓應力為84.7 MPa。

圖5 6#采空區(qū)等效應力圖

6#采空區(qū)由A、B、C三個采空區(qū)構(gòu)成，在三個采空區(qū)上方巖體主要受到拉應力作用，最大拉應力為20.136 MPa，均大于巖體抗拉強度，說明6#采空區(qū)上部及四周已發(fā)生拉伸破壞，這與實際觀測結(jié)果相符。從圖5中還可看出，A采空區(qū)頂板圍巖受到拉應力大于B、C空區(qū)，這是由于A采空區(qū)體積大、跨度大所致。

4 結(jié) 論

采用ANSYS有限元分析軟件，分別對用沙壩礦3#、6#采空區(qū)周圍巖體的穩(wěn)定性進行數(shù)值分析，得到如下結(jié)論：

(1) 3#采空區(qū)和6#采空區(qū)上方巖體受拉應力作用明顯，所受最大拉應力分別為13.581 MPa和20.136 MPa，大于巖石抗拉強度，說明3#、6#空區(qū)均不穩(wěn)定，易發(fā)生拉伸破壞；

(2) 采空區(qū)附近巖體主要受拉伸破壞，6#采空區(qū)圍巖所受拉應力大于3#采空區(qū)，說明采空區(qū)體積增大，采空區(qū)上部圍巖所受拉應力相對增大。

參考文獻：

[1]馮長根, 李俊平, 于文遠, 等. 東桐峪金礦空場處理機制研究[J]. 黃金, 2002, 23(10): 11-15.

[2]于學馥, 鄭穎人. 地下工程圍巖穩(wěn)定分析[M]. 北京: 煤炭工業(yè)出版社, 1983.

[3]唐勝利,唐 皓,郭 輝. 基于BP神經(jīng)網(wǎng)絡的空洞型采空區(qū)穩(wěn)定性評價研究[J]. 西安科技大學學報, 2012, 32(2): 234-238.

[4]程愛寶, 王新民, 劉洪強. 灰色層次分析法在地下采空區(qū)穩(wěn)定性評價中的應用[J]. 金屬礦山, 2011(2): 17-21.

[5]魏海波, 吳 敏. 邊坡的有限元分析及ANSYS軟件對邊坡開挖的模擬[J]. 云南水力發(fā)電, 2004, 20(4): 42-44.

[6]孫永剛. 有限元基礎教學的探索[J]. 科技信息, 2010(24): 517-518.