王飛飛，鄒 平，孟中華，李?lèi)?ài)兵，劉正宇，虎萬(wàn)杰，馬 增

(1.長(zhǎng)沙礦山研究院有限責(zé)任公司，長(zhǎng)沙 410012；2.金屬礦山安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410012)

礦山采空區(qū)的穩(wěn)定性一直是礦山安全生產(chǎn)的關(guān)鍵性因素，關(guān)于礦山采空區(qū)穩(wěn)定性問(wèn)題，有眾多學(xué)者做了相關(guān)研究[1-8]，并提出了處理采空區(qū)的一些方法與技術(shù)，例如：崩落法處理采空區(qū)、充填法處理采空區(qū)及封閉處理采空區(qū)。由于地下采空區(qū)具有隱伏性強(qiáng)、空間分布特征規(guī)律性差、采空區(qū)頂板冒落以及采空區(qū)塌陷情況難以預(yù)測(cè)等特點(diǎn)。因此，準(zhǔn)確地分析采空區(qū)穩(wěn)定性是關(guān)乎礦山能否安全生產(chǎn)的重要因素。

關(guān)于采空區(qū)的穩(wěn)定性，大多學(xué)者主要研究的是影響采空區(qū)穩(wěn)定的因素，對(duì)采空區(qū)礦柱在上部巖體作用下?tīng)顟B(tài)的研究尚少。本文采用MIDAS-GTS/NX有限元數(shù)值模擬軟件，建立了采空區(qū)—礦柱—礦山三維力學(xué)分析模型，通過(guò)非線(xiàn)性靜力分析得到了礦柱剪應(yīng)力，礦柱軸向應(yīng)力與礦山塑性區(qū)分布。研究結(jié)果可以為相似礦山提供參考。

1 礦山概況

1.1 礦山地質(zhì)

大橋礦區(qū)磷礦礦體呈灰黑色薄層狀，產(chǎn)于震旦系上統(tǒng)燈影組麥地坪段含磷白云巖中部，含礦層位巖性特征為：上、下部灰白色薄至中厚層狀含磷砂屑白云巖-粉晶白云巖，中部為深灰—灰黑色膠磷礦磷塊巖。礦層呈層狀產(chǎn)出，層位穩(wěn)定，礦體沿大渡河?xùn)|岸陡坡地段分布，呈南北向出露，礦體長(zhǎng) 4 250 m，厚度 0.42～1.02 m，平均厚 0.57 m，傾向 111°～194°，傾角 5°～10°。礦層底板為灰白色薄層至中厚層狀含磷白云巖夾燧石條帶，頂板為灰白色薄層至中厚層狀含磷白云巖。礦層為深灰、灰黑色磷塊巖，地表及深部礦層厚度變化規(guī)律不明顯，礦層層位穩(wěn)定。礦層圍巖主要為含磷白云巖，與礦體界線(xiàn)清楚。

1.2 礦山采空區(qū)

大橋磷礦已經(jīng)過(guò)三十多年的有組織無(wú)規(guī)劃、采富棄貧等不同程度的開(kāi)采作業(yè)，形成了數(shù)量巨大的、不規(guī)則的采空區(qū)。形成采空區(qū)面積約為45萬(wàn)m2，留下了許多不規(guī)則礦柱，部分已經(jīng)變形，空區(qū)頂板出現(xiàn)少量冒落現(xiàn)象。經(jīng)現(xiàn)場(chǎng)勘查，近些年新增采空區(qū)48萬(wàn)m3左右，但通過(guò)人工電耙配合，使用前期廢石、廢碴充填接頂，后階段在荒石中攪拌入5%水泥通過(guò)人工電耙配合充填接頂，已治理約20萬(wàn)m3，目前尚存在大面積的采空區(qū)未治理。隨著礦山開(kāi)采工作的繼續(xù)，很難保證原有的老采空區(qū)穩(wěn)定性。因此，礦山老采空區(qū)亟待開(kāi)展礦柱穩(wěn)定性研究工作。

2 建立模型

2.1 非線(xiàn)性分析理論

非線(xiàn)性有限元解法是一種通過(guò)累積增量解收斂于正確解的方法，計(jì)算過(guò)程如圖1所示。

圖1 增量法Fig.1 Incremental method

圖中tfext和t+Δtfext為t時(shí)刻和t+Δt時(shí)刻的外力，兩個(gè)時(shí)刻的解和增量解關(guān)系如式(1)：

t+Δtμ=tμ+Δμ

(1)

式中：Δμ為時(shí)間增量Δt產(chǎn)生的增量解，累積增量解表達(dá)式如式(2)：

Δμi+1=Δμi+δμi+1

(2)

式中：Δμi為直到第i次迭代計(jì)算的累計(jì)增量解；δμi+1為第i+1次迭代計(jì)算時(shí)的增量解。δμi+1通過(guò)式(3)計(jì)算得到，Ki+1為切線(xiàn)剛度矩陣。

(3)

gi為殘余力、非平衡力。非平衡力gi表達(dá)式如下：

gi=t+Δtfext-fint,j

(4)

式中t+Δtfext為外力，fint,j為內(nèi)力。對(duì)方程(2～4)進(jìn)行重復(fù)迭代計(jì)算，直到滿(mǎn)足收斂準(zhǔn)則為止，收斂準(zhǔn)則包括力準(zhǔn)則、位移準(zhǔn)則和能量準(zhǔn)則等。

為了改善上述的迭代求解問(wèn)題，MIDAS -NX采用線(xiàn)性搜索功能來(lái)解決該問(wèn)題，線(xiàn)性搜索算法(如圖2所示)是在把增量δμi+1累計(jì)到增量解時(shí)采用一個(gè)比例值概念η來(lái)改善精度，此時(shí)累計(jì)增量解計(jì)算如式(5)：

Δμi+1=Δμi+ηδμi+1

(5)

假定Δμi+1滿(mǎn)足平衡狀態(tài)，利用總勢(shì)能不變?cè)?，可將線(xiàn)性搜索問(wèn)題歸結(jié)為解總勢(shì)能對(duì)η導(dǎo)數(shù)為零時(shí)的η值問(wèn)題，如式(6)所示：

(6)

圖2 線(xiàn)性搜索算法Fig.2 Linear search algorithm

假定能量對(duì)η的導(dǎo)數(shù)是線(xiàn)性變化，則滿(mǎn)足方程(6)的η計(jì)算見(jiàn)式(7)：

(7)

η=0或1時(shí)的斜率表達(dá)式如式(8)、(9)：

(8)

(9)

雖然線(xiàn)性搜索法不能十分準(zhǔn)確地滿(mǎn)足實(shí)際情況，按照式(8)、(9)計(jì)算s(η)一般不為零，但在MIDAS-NX中上面所述的過(guò)程會(huì)重復(fù)進(jìn)行直到s(ηi)/s(η=0)的值低于設(shè)定值為止，即達(dá)到預(yù)定的計(jì)算精度。

2.2 采空區(qū)礦柱模型建立

采用MIDAS-GTS/NX有限元模擬軟件建立采空區(qū)—礦柱—礦山三維力學(xué)分析模型，數(shù)值模型大小與原型一致。礦柱與巖體均采用實(shí)體單元模擬，采用摩爾—庫(kù)倫屈服準(zhǔn)則[9]，每個(gè)采空區(qū)的埋深均參照實(shí)際情況建立，礦柱埋深100～200 m，采空區(qū)高度1～2 m。采用四面體單元?jiǎng)澐秩S模型，礦柱與圍巖采用實(shí)體單元模擬。計(jì)算模型劃分節(jié)點(diǎn)總數(shù)有121 483個(gè)，單元總數(shù)有651 810個(gè)，模型為各項(xiàng)同性均質(zhì)理想的線(xiàn)彈塑性巖土體。模型底部與四周采用固定邊界，上部采用自由面邊界。三維模型分為三層，上層為白云巖圍巖，中部為含有采空區(qū)的礦體，下層為與上層同樣的巖體。礦山的巖體力學(xué)參數(shù)如表1所示。鑒于礦山采空區(qū)較多，限于篇幅，本文僅對(duì)大橋磷礦的5-2、9-2與320三個(gè)采空區(qū)進(jìn)行分析。所建立的三維礦山力學(xué)分析模型如圖3所示。

表1 材料物理力學(xué)參數(shù)

圖3 三維礦山數(shù)值模型Fig.3 Three-dimensional mine numerical model

3 數(shù)值模擬結(jié)果與分析

對(duì)建立的不同采空區(qū)三維數(shù)值模型進(jìn)行非線(xiàn)性靜力分析，得到了礦柱剪應(yīng)力云圖，礦柱軸應(yīng)力云圖與礦山塑性區(qū)分布圖。

3.1 礦柱剪應(yīng)力

礦柱剪應(yīng)力云圖如圖4所示。

在圖4(a)中，圖例中從上到下剪應(yīng)力由高到低，5-2采空區(qū)礦柱最大剪應(yīng)力為34.9 MPa，最小剪應(yīng)力為10 MPa。礦山5-2采空區(qū)埋深100～200 m。由圖可知不同礦柱的剪應(yīng)力大小差異較大。有的較大橫截面的礦柱剪應(yīng)力較大，而有的較小的橫截面的礦柱剪應(yīng)力較小，可以得到礦柱剪應(yīng)力大小與礦柱橫截面的大小沒(méi)有直接的關(guān)系。剪應(yīng)力較大或較小的礦柱集中在一個(gè)區(qū)域內(nèi)，說(shuō)明礦柱剪應(yīng)力的大小與上覆巖層厚度有直接聯(lián)系，較厚的上覆巖層的重力較大導(dǎo)致礦柱承受的剪切力較大。最大剪力產(chǎn)生在小礦柱中間位置處，這與現(xiàn)場(chǎng)礦柱的破壞形式非常吻合，驗(yàn)證了數(shù)值模擬計(jì)算得到的結(jié)果是可信的。

由圖4(b)可知，9-2采空區(qū)礦柱最大剪應(yīng)力為39.9 MPa，最小剪應(yīng)力為14 MPa。礦山9-2采空區(qū)埋深260～360 m，與5-2采空區(qū)相比，其埋深增加了160 m左右，故9-2采空區(qū)礦柱的剪應(yīng)力較大。9-2采空區(qū)中，有個(gè)別礦柱的剪應(yīng)力較大，而周邊的礦柱剪應(yīng)力較小，可能是因?yàn)榈V柱分布導(dǎo)致上覆壓力對(duì)個(gè)別礦柱產(chǎn)生集中力，最終極個(gè)別的礦柱剪應(yīng)力較大。9-2采空區(qū)礦柱剪應(yīng)力在整體采空區(qū)分布較為均勻，主要是因?yàn)樵摬煽諈^(qū)上部山體高低起伏度不大，對(duì)該采空區(qū)產(chǎn)生較為均勻的壓力。

由圖4(c)可知，320采空區(qū)礦柱最大剪應(yīng)力為38.9 MPa，最小剪應(yīng)力為14.8 MPa。礦山320采空區(qū)最大埋深在150～230 m范圍內(nèi)，其埋深比9-2采空區(qū)小，比5-2采空區(qū)大，故320采空區(qū)礦柱最大剪應(yīng)力比9-2采空區(qū)礦柱剪應(yīng)力小，比5-2采空區(qū)礦柱剪應(yīng)力大。

通過(guò)比較圖4中三個(gè)采空區(qū)礦柱剪應(yīng)力大小與分布可知，采空區(qū)剪應(yīng)力大小順序?yàn)椋?-2采空區(qū)>320采空區(qū)>5-2采空區(qū)，相比較5-2采空區(qū)礦柱較為穩(wěn)定。采空區(qū)礦柱剪應(yīng)力與礦柱的埋深有直接的聯(lián)系。

3.2 礦柱軸向應(yīng)力

采空區(qū)礦柱軸向應(yīng)力云圖如圖5所示。

由圖5(a)可知，5-2采空區(qū)礦柱的最大軸向應(yīng)力為84.44 MPa，最小軸向應(yīng)力為22.93 MPa。由圖5(b)可知，9-2采空區(qū)礦柱的最大軸力為100.75 MPa，最小軸力為39.56 MPa。由圖5(c)可知，320采空區(qū)礦柱的最大軸力為98.78 MPa，最小軸力為42.19 MPa。通過(guò)比較以上三個(gè)采空區(qū)礦柱的軸力，可知5-2采空區(qū)較為穩(wěn)定，其它兩個(gè)采空區(qū)具有不安全可能性，因此，要加強(qiáng)監(jiān)測(cè)與防護(hù)措施。

3.3 礦山塑性區(qū)

礦山塑性區(qū)分布如圖6所示。

圖4 采空區(qū)礦柱剪應(yīng)力Fig.4 Shear stress of pillar in goaf

圖5 礦柱軸向應(yīng)力云圖Fig.5 Axial stress cloud map of the pillar

圖6 礦山塑性區(qū)分布Fig.6 Distribution map of plastic zone in mine

大橋磷礦礦體賦存傾角5°～10°，礦體上盤(pán)有向山體內(nèi)側(cè)移動(dòng)的趨勢(shì)，不會(huì)向山體外側(cè)移動(dòng)，即不會(huì)發(fā)生山體垮塌堵塞大渡河的危險(xiǎn)。由圖6可知，礦山?jīng)]有產(chǎn)生大面積的塑性區(qū)分布，僅在礦體處有局部塑性區(qū)。整體而言，山體穩(wěn)定性良好。

4 結(jié)論

本文采用MIDAS-GTS/NX有限元數(shù)值模擬軟件，建立了采空區(qū)—礦柱—礦山三維力學(xué)分析模型，采用非線(xiàn)性靜力分析方法對(duì)大橋磷礦采空區(qū)開(kāi)展了三維數(shù)值模擬研究，得到以下結(jié)論：

1)礦柱剪應(yīng)力大小與礦柱橫截面的大小沒(méi)有直接的關(guān)系，礦柱剪應(yīng)力的大小與上覆巖層厚度有直接聯(lián)系。采空區(qū)剪應(yīng)力大小順序?yàn)椋?-2采空區(qū)>320采空區(qū)>5-2采空區(qū)，相比較5-2采空區(qū)礦柱較為穩(wěn)定。

2)采空區(qū)埋深對(duì)采空區(qū)礦柱的剪應(yīng)力與礦柱軸向應(yīng)力均有很大的影響，上部巖體的重力較小，礦柱產(chǎn)生的剪應(yīng)力與軸向應(yīng)力較小。

3)礦山?jīng)]有產(chǎn)生大面積的塑性區(qū)分布，僅在礦體處有局部塑性區(qū)。整體而言，山體穩(wěn)定性良好。