王永奇,戴 兵

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

0 引 言

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

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

1 礦山概況

開陽(yáng)磷礦自1958年建礦以來，已連續(xù)生產(chǎn)了50多年。以往礦山使用的采礦方法為空?qǐng)龇?，采后采空區(qū)未進(jìn)行處理。礦區(qū)形成的采空區(qū)，導(dǎo)致山體塌陷和滑坡。礦體延伸15 km，從地表可見綿延15 km的山體塌陷區(qū)。礦體走向南北，傾向?yàn)?80°～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ù)值軟件采用有限元計(jì)算程序，集結(jié)構(gòu)、熱、流體、電磁場(chǎng)、聲場(chǎng)和耦合場(chǎng)分析于一體，能處理線性和非線性復(fù)雜力學(xué)計(jì)算問題，在開發(fā)之初便得到廣泛應(yīng)用[5]。

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

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

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

{f}=[N]{δ}e

(1)

式中，{f}—單元中任意一點(diǎn)的位移，

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

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

(3)單元力學(xué)特性。首先推導(dǎo)出單元應(yīng)變(利用節(jié)點(diǎn)位移)表達(dá)式：

{ε}=[B]{δ}e

(2)

式中，{ε}—單元應(yīng)變，

[B]—單元應(yīng)變矩陣。

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

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

(3)

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

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

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

(4)

式中，[k]e—單元?jiǎng)偠染仃嚕?/p>

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

(5)

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

[K]{σ}=[F]

(6)

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

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

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

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

本文主要對(duì)開磷礦區(qū)3#采空區(qū)和6#采空區(qū)穩(wěn)定性進(jìn)行數(shù)值模擬，3#空區(qū)尺寸為120 m×50 m×17.5 m；6#采空區(qū)主要由A、B、C三個(gè)小采空區(qū)構(gòu)成，三個(gè)空區(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ū)跨度極大，計(jì)算按平面應(yīng)變問題進(jìn)行。由于巖石具有脆性，分析中涉及到的所有物理量均與時(shí)間無關(guān)。

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

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

3.3 建立計(jì)算模型

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

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

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

3.4 邊界條件與荷載

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

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

原巖自重應(yīng)力可分別按下式進(jìn)行計(jì)算：

(7)

(8)

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

3.5 模擬結(jié)果分析

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

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

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

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

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

4 結(jié) 論

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

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

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

參考文獻(xiàn)：

[1]馮長(zhǎng)根, 李俊平, 于文遠(yuǎn), 等. 東桐峪金礦空?qǐng)鎏幚頇C(jī)制研究[J]. 黃金, 2002, 23(10): 11-15.

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

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

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

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

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