任青云

(湖南有色金屬研究院，長(zhǎng)沙 410100)

地下礦山開(kāi)采的重點(diǎn)安全問(wèn)題是開(kāi)采形成的采空區(qū)對(duì)后續(xù)開(kāi)采產(chǎn)生的巖層移動(dòng)問(wèn)題，開(kāi)采過(guò)程中伴隨的地層沉降及地壓活動(dòng)[1-2]。礦山開(kāi)采面臨的重要安全問(wèn)題主要表現(xiàn)在，開(kāi)采區(qū)域的不斷擴(kuò)大，地層隨之進(jìn)行持續(xù)的運(yùn)動(dòng)，由于開(kāi)采現(xiàn)場(chǎng)會(huì)出現(xiàn)不同的地質(zhì)條件，例如溶洞、斷層、空區(qū)等可能會(huì)產(chǎn)生移動(dòng)變形的大范圍空間，此時(shí)會(huì)造成地下礦山的地壓顯現(xiàn)，導(dǎo)致頂板出現(xiàn)大變形，繼而影響圍巖發(fā)生冒落和片幫等現(xiàn)象。因此地下礦山長(zhǎng)期開(kāi)采導(dǎo)致的采空區(qū)穩(wěn)定性是影響礦山生命財(cái)產(chǎn)安全及后續(xù)開(kāi)采的重要因素[3-4]。

影響采空區(qū)的穩(wěn)定性及礦山地表位移沉降的主要原因有開(kāi)采礦體的埋深、采場(chǎng)的高度，礦體的厚度、傾角、巖石力學(xué)特性及礦山采用的采礦工藝。采用大型有限元數(shù)值模擬[5-6]前處理模式可以直觀地表示出地下采空區(qū)的分布情況及整改礦山的輪廓，其后處理模式能夠計(jì)算出礦山的地表位移沉降、采空區(qū)是否發(fā)生拉應(yīng)力破壞、礦柱是否產(chǎn)生主應(yīng)力集中及塑性區(qū)變形破壞情況。本文根據(jù)南川河石灰?guī)r礦礦區(qū)地質(zhì)資料，結(jié)合礦體情況及開(kāi)采設(shè)計(jì)，對(duì)南川河石灰?guī)r礦進(jìn)行取樣試驗(yàn)，采用數(shù)值模擬方法分析礦山開(kāi)采后巖層移動(dòng)變形情況，并結(jié)合傳統(tǒng)的理論計(jì)算方法對(duì)礦柱、頂板穩(wěn)定性進(jìn)行分析，以此確定無(wú)支護(hù)狀態(tài)下采空區(qū)跨度，從而評(píng)判采空區(qū)的穩(wěn)定性及地表位移沉降的影響范圍和程度。

1 工程概況

南川河石灰?guī)r礦位于瀏陽(yáng)市澄潭江鎮(zhèn)南川河畔，距瀏陽(yáng)市城區(qū)直距 21 km，公司始建于1972年，原采用露天開(kāi)采，1996年開(kāi)始地下開(kāi)采，主要采用斜坡道開(kāi)拓，對(duì)角式通風(fēng)。礦山長(zhǎng)期使用的采礦方法是房柱采礦法。礦山總的回采順序是沿礦體傾斜由上而下循階段進(jìn)行，按邊探邊采的原則進(jìn)行開(kāi)采。采場(chǎng)回采順序是后退式開(kāi)采，礦房?jī)?nèi)采用分層開(kāi)采。礦山2014年之前開(kāi)采130 m中段，2014年至今開(kāi)采100 m中段，目前礦山正在+100 m中段北部進(jìn)行生產(chǎn)。

由于石灰?guī)r的巖性較軟，并且礦山開(kāi)采歷史悠久，地下空區(qū)存在極大的不穩(wěn)定性及不確定性，為了掌握空區(qū)對(duì)下部礦體回采的影響，必須針對(duì)采空區(qū)進(jìn)行穩(wěn)定性研究和分析。根據(jù)礦山的實(shí)際生產(chǎn)現(xiàn)狀，針對(duì)空區(qū)進(jìn)行了現(xiàn)場(chǎng)調(diào)查，并通過(guò)統(tǒng)計(jì)分析獲得礦山不同開(kāi)采水平的采空區(qū)規(guī)模、尺寸。此次現(xiàn)場(chǎng)調(diào)查針對(duì)+130 m中段和+100 m中段共統(tǒng)計(jì)出130個(gè)采空區(qū)，+100 m中段采空區(qū)部分調(diào)查數(shù)據(jù)見(jiàn)表1?，F(xiàn)場(chǎng)調(diào)查發(fā)現(xiàn)，礦山部分采空區(qū)已形成采空區(qū)群；大部分空區(qū)寬度約12～14 m，采空區(qū)未發(fā)生冒頂片幫事故，側(cè)壁處理較好，有極少量的浮石[7]。溶洞構(gòu)造溶蝕形成，并有軟塑性、硬塑溶洞泥質(zhì)充填，與地表水源連通的構(gòu)造少。

表1 +100 m中段采空區(qū)特征統(tǒng)計(jì)表(部分)

2 礦巖力學(xué)參數(shù)室內(nèi)試驗(yàn)

2.1 礦巖力學(xué)參數(shù)試驗(yàn)

為了合理評(píng)價(jià)南川河石灰?guī)r礦工程地質(zhì)條件，為后續(xù)數(shù)值模擬提供依據(jù)，開(kāi)展了礦巖力學(xué)參數(shù)室內(nèi)試驗(yàn)，測(cè)試的內(nèi)容包括礦巖的容重、單軸抗壓強(qiáng)度、抗拉強(qiáng)度、彈性模量、泊松比等基本物理力學(xué)參數(shù)。針對(duì)礦山實(shí)際情況，選取了石灰?guī)r巖樣進(jìn)行加工，并對(duì)各巖樣進(jìn)行了相關(guān)參數(shù)測(cè)定。巖樣總計(jì)22個(gè)，其中5個(gè)用于開(kāi)展單軸壓縮試驗(yàn)、5個(gè)用于巴西劈裂試驗(yàn)、12個(gè)用于變角剪試驗(yàn)。

根據(jù)數(shù)值模擬計(jì)算分析需要的基本數(shù)據(jù)，主要進(jìn)行了如下四項(xiàng)內(nèi)容的實(shí)驗(yàn)室試驗(yàn)：1)礦、巖的天然容重；2)礦、巖的單軸抗壓力學(xué)特性；3)巖體的劈裂抗拉強(qiáng)度試驗(yàn)；4)巖體的泊松比與彈性模量試驗(yàn)。

試驗(yàn)開(kāi)始前，先對(duì)巖樣加工。首先，按照標(biāo)準(zhǔn)巖石力學(xué)試驗(yàn)所要求的規(guī)格，制成試樣尺寸為Φ50 mm×100 mm的標(biāo)準(zhǔn)巖石試樣；然后擦干試樣，在自然狀態(tài)下風(fēng)干；最后，試驗(yàn)在250 t全數(shù)字型液壓伺服剛性巖石力學(xué)試驗(yàn)系統(tǒng)(RMT150，圖1)及30 t萬(wàn)能材料試驗(yàn)機(jī)上完成。

圖1 RMT150型試驗(yàn)設(shè)備Fig.1 RMT150 test equipment

2.1.1 單軸抗壓試驗(yàn)

試樣尺寸為直徑Ф50 mm，高徑比約為2∶1，巖石單軸抗壓試驗(yàn)結(jié)果如表2所示。

表2 巖石單軸抗壓強(qiáng)度及靜力受壓彈性模量和泊松比測(cè)試

2.1.2 抗拉強(qiáng)度試驗(yàn)

巖石的抗拉強(qiáng)度試驗(yàn)如圖2所示，各種試驗(yàn)結(jié)果見(jiàn)表3和表4。

圖2 巴西劈裂試驗(yàn)Fig.2 Brazilian splitting experiment

表3 各試樣劈裂抗拉試驗(yàn)參數(shù)測(cè)定值

表4 變角剪試驗(yàn)結(jié)果

根據(jù)單軸抗壓試驗(yàn)和單軸抗拉試驗(yàn)，將室內(nèi)試驗(yàn)結(jié)果匯總?cè)绫?所示。

表5 礦巖力學(xué)參數(shù)試驗(yàn)結(jié)果匯總表

2.2 巖體強(qiáng)度參數(shù)分析

由于自然巖體內(nèi)部有很多缺陷，比如節(jié)理、裂隙和孔洞等，因此自然巖體與小試樣巖石的力學(xué)性質(zhì)是有差異的，導(dǎo)致實(shí)驗(yàn)室測(cè)得小試樣巖石力學(xué)參數(shù)不能直接作為數(shù)值模擬的力學(xué)參數(shù)。因此，巖體力學(xué)工程界，大多數(shù)依據(jù)室內(nèi)巖石力學(xué)實(shí)驗(yàn)，采用理論力學(xué)的方法，估算出巖體的力學(xué)參數(shù)，用于有限元數(shù)值模擬計(jì)算。

根據(jù)HOEK，CARRANZA-TORRES提出的巖體破壞經(jīng)驗(yàn)準(zhǔn)則：

(1)

巖體的mb、S和a由GSI和上式計(jì)算的mi確定。

(2)

(3)

(4)

根據(jù)確定的mb、S和a值，可以計(jì)算出巖體的單軸抗壓強(qiáng)度和抗拉強(qiáng)度：

σc=σci·sa

(5)

(6)

巖體的內(nèi)摩擦角φ和黏聚力c可以由以下公式進(jìn)行計(jì)算：

(7)

(8)

針對(duì)南川河石灰?guī)r礦工程地質(zhì)條件，取石灰?guī)rGSI=70，根據(jù)室內(nèi)巖石力學(xué)參數(shù)試驗(yàn)結(jié)果，利用HOEK和CARRANZA-TORRES提出的巖體破壞經(jīng)驗(yàn)準(zhǔn)則及巖體分類(lèi)RMR值，采用Roclab軟件進(jìn)行計(jì)算。選取巖體的力學(xué)參數(shù)如表6所示。

表6 巖體物理力學(xué)參數(shù)匯總表

3 采空區(qū)穩(wěn)定性數(shù)值分析

3.1 三維礦體模型的建立

制約礦山地表位移沉降及采空區(qū)穩(wěn)定性的重要因素，包括礦區(qū)地表的地形地質(zhì)條件及井下采空區(qū)的分布狀態(tài)，因此在三維有限元數(shù)值模擬計(jì)算過(guò)程中，缺少完整的地形模塊，不能真實(shí)地反映礦山的實(shí)際情況。本研究通過(guò)基于DIMINE強(qiáng)度的前處理功能，利用二維地質(zhì)剖面圖或者M(jìn)APGIS地形數(shù)據(jù)庫(kù)生成三維地質(zhì)模型，對(duì)勘探線(xiàn)剖面進(jìn)行初步的處理，僅保留勘探線(xiàn)與礦體解譯線(xiàn)[8]。將礦山剖面圖經(jīng)過(guò)CAD處理之后導(dǎo)入到DIMINE中，進(jìn)行線(xiàn)串運(yùn)算使勘探線(xiàn)“立起來(lái)”，對(duì)立起來(lái)的三維線(xiàn)文件進(jìn)行坐標(biāo)運(yùn)算，使各點(diǎn)坐標(biāo)為其真實(shí)坐標(biāo)值。通過(guò)上述基本工作的完成，可以聯(lián)接各礦體解譯線(xiàn)剖面，對(duì)于比較復(fù)雜的礦體解譯線(xiàn)，可以構(gòu)造中線(xiàn)，通過(guò)中線(xiàn)與剖面線(xiàn)來(lái)生成實(shí)體，礦體實(shí)體模型見(jiàn)圖3。在DIMINE中生成的礦山實(shí)體模型，無(wú)法直接用于有限元數(shù)值模擬的運(yùn)算，需要將此模型導(dǎo)入到MIDAS/GTS中，利用MIDAS/GTS的網(wǎng)格劃分功能，結(jié)合采礦工藝，進(jìn)一步建立采區(qū)模型，劃分網(wǎng)格單元，進(jìn)行數(shù)值模擬分析。

圖3 南川河石灰?guī)r礦三維實(shí)體模型Fig.3 3D solid model of Nanchuanhe limestone mine

3.2 理論基礎(chǔ)及相關(guān)參數(shù)設(shè)定

根據(jù)礦山的地層參數(shù)，確定有限元數(shù)值模擬的巖體力學(xué)參數(shù)見(jiàn)表6。計(jì)算過(guò)程礦巖及廢石均采用摩爾-庫(kù)倫((Mohr-Coulomb)屈服準(zhǔn)則，該屈服準(zhǔn)則的控制方程為：

(9)

最大拉應(yīng)力屈服準(zhǔn)則函數(shù)為：

ft=σ3-σt

(10)

式中：c—黏聚力，MPa；φ—內(nèi)摩擦角，(°)；σ1—最大主應(yīng)力，MPa；σ3—最小主應(yīng)力，MPa；σt—抗拉強(qiáng)度，MPa。MIDAS/GTS有限元計(jì)算程序，默認(rèn)當(dāng)fs大于0時(shí)，表明巖體已發(fā)生了剪切破壞；當(dāng)ft大于0時(shí)，表明巖體已發(fā)生了拉伸破壞。

為了模型建立與分析結(jié)果的準(zhǔn)確性，本數(shù)值模擬進(jìn)行了以下假設(shè)：模型的礦巖體為均質(zhì)、各向同性材料，模擬的三維網(wǎng)格模型底部為固定約束，四周為限制水平約束，只在豎直方向有位移變形。初始階段考慮地應(yīng)力與自重的影響，位移清零，表明礦山在原始狀態(tài)已處于穩(wěn)定的固結(jié)沉降。計(jì)算模擬不考慮地震、爆破等動(dòng)荷載的影響，不分析地下水的影響；礦山地質(zhì)條件較復(fù)雜，礦巖穩(wěn)固性在空間分布上具有較大的隨機(jī)性，礦體為斷層分割，且礦巖和圍巖中均有穿插，這些地質(zhì)現(xiàn)象的影響在巖體參數(shù)折減時(shí)已做考慮，在模擬過(guò)程中不予另行考慮。

根據(jù)地下采空區(qū)對(duì)礦山圍巖的影響范圍，設(shè)置有限元數(shù)模網(wǎng)格模型的邊界固定條件如下：1)模型X方向的邊界為X方向固定邊界，取u=0，v≠0，w≠0(u為X軸方向的位移，v為Y軸方向的位移，w為Z軸方向的位移)；2)模型Y方向邊界為Y方向固定邊界，取u≠0，v=0，w≠0；3)模型底部為全固定邊界，取u=0，v=0，w=0；4)模型地表部分為自由面，不限制位移約束。礦體邊界約束情況及三維網(wǎng)絡(luò)質(zhì)量情況見(jiàn)圖4和圖5。

圖4 三維網(wǎng)格模型+邊界條件Fig.4 3D grid model + boundary conditions

圖5 網(wǎng)格質(zhì)量檢查Fig.5 Grid quality inspection

3.3 數(shù)值模擬計(jì)算結(jié)果分析

3.3.1 位移場(chǎng)分析

采空區(qū)的穩(wěn)定性分析，位移是評(píng)判礦山穩(wěn)定性的一個(gè)重要指標(biāo)，因?yàn)榇罅康牡乇碜冃危瑫?huì)使礦山發(fā)生地表塌陷，且位移變形量超過(guò)巖體的變形限值，會(huì)發(fā)生采空區(qū)頂板冒頂片幫事故。從圖6～9的位移云圖得出，該礦山產(chǎn)生的位移變化情況，在礦體的頂板位移值小于0，說(shuō)明由于產(chǎn)生了采空區(qū)，發(fā)生了頂板下沉變形情況，但隨著礦體的開(kāi)挖過(guò)程進(jìn)行，在采空區(qū)的底部及巷道的底板，豎直方向位移大于0，是因?yàn)榈V體開(kāi)挖深度加大，采空區(qū)及巷道周?chē)l(fā)生了主應(yīng)力集中現(xiàn)象，水平方向的應(yīng)力大于垂直方向的應(yīng)力，使巷道底板產(chǎn)生上拱的位移變形。

從圖6豎向位移云圖可以看出，+130 m中段礦體開(kāi)采時(shí)，采空區(qū)頂板豎向位移最大值為-1.74 mm，礦體開(kāi)采底板巷道發(fā)生向上最大位移為+2.39 mm。

圖6 +130 m中段頂板位移云圖Fig.6 Displacement cloud image of the +130 m middle section roof

圖7為礦柱的位移云圖，+130 m中段采空區(qū)礦柱位移變形為受壓沉降，最大位移值為-1.8 mm，+100 m中段礦柱由于+130 m中段礦體已全部回采，+100 m礦柱所受側(cè)向壓力大于豎向壓力，變形為部分礦柱產(chǎn)生上拱位移，最大位移值為+2.5 mm。

圖7 礦柱位移云圖Fig.7 Cloud diagram of pillar displacement

圖8為+100 m中段礦體開(kāi)采后采空區(qū)頂板產(chǎn)生的位移沉降變形云圖，目前在+100 m中段回采時(shí)采空區(qū)頂板最大位移沉降為-2.4 mm，采空區(qū)底板巷道上拱位移值為+2.94 mm。

圖8 +100 m中段空區(qū)頂板位移云圖Fig.8 Displacement cloud image of the +100 m middle section roof

采空區(qū)形成的地表沉降有微弱的影響，地表沉降值不足1 mm，該位移沉降整體上可以忽略，采空區(qū)地表沉降等值線(xiàn)范圍見(jiàn)圖9，位移變形值數(shù)據(jù)見(jiàn)表7。

圖9 地表沉降位移云圖與沉降等值線(xiàn)圖Fig.9 Surface displacement cloud map and subsidence contour map

表7 采空區(qū)位移變化值計(jì)算結(jié)果

3.3.2 應(yīng)力場(chǎng)分析

礦體開(kāi)采后對(duì)圍巖造成的影響是一個(gè)非線(xiàn)性的不可逆的加載過(guò)程，它使得處于初始地應(yīng)力狀態(tài)下的圍巖進(jìn)行應(yīng)力重新分布，最后達(dá)到新的平衡，礦體開(kāi)挖引起徑向應(yīng)力釋放、切向應(yīng)力增加，導(dǎo)致礦體開(kāi)采頂板和底板產(chǎn)生拉應(yīng)力；從應(yīng)力云圖可以看出礦體開(kāi)采區(qū)域與礦柱接觸處發(fā)生較大的應(yīng)力集中，周邊最大。采空區(qū)影響范圍一定深度后，逐步變化為恢復(fù)到原巖應(yīng)力狀態(tài)。采空區(qū)形成后的主應(yīng)力值都小于0，表明主應(yīng)力(壓應(yīng)力)為受壓狀態(tài)。

圖10 最小主應(yīng)力云圖Fig.10 The minimum principal stress cloud diagram

圖11 最大主應(yīng)力云圖Fig.11 The maximum principal stress cloud diagram

根據(jù)計(jì)算結(jié)果圖10、圖11得出結(jié)論：在采空區(qū)的頂板和底板產(chǎn)生拉應(yīng)力破壞，其中+130 m中段采空區(qū)頂板產(chǎn)生的最大拉應(yīng)力約為2.27 MPa；+100 m中段采空區(qū)頂板產(chǎn)生的最大拉應(yīng)力為1.85 MPa，已接近巖體的抗拉強(qiáng)度，表明在采空區(qū)的頂板部分已發(fā)生拉應(yīng)力破壞，容易產(chǎn)生冒頂片幫事故。從最大應(yīng)力云圖可以得出：在礦柱與空區(qū)的邊界圍巖出現(xiàn)壓應(yīng)力集中現(xiàn)象，最大壓應(yīng)力集中發(fā)生在礦山+100 m中段東北方向的采空區(qū)與13號(hào)礦柱交界處，產(chǎn)生的最大壓應(yīng)力值為12.66 MPa，+130 m中段的(老采空區(qū))礦柱發(fā)生最大壓應(yīng)力集中值為8.88 MPa，表明礦柱發(fā)生受壓破壞。應(yīng)力結(jié)果見(jiàn)表8。

表8 最大、最小主應(yīng)力值計(jì)算結(jié)果

3.3.3 塑性區(qū)分析

圖12為分析得出的采空區(qū)塑性變形情況，從計(jì)算結(jié)果可以得出，在采空區(qū)的附近發(fā)生了塑性區(qū)變形破壞，塑性區(qū)面積占總面積約44.5%，最大塑性變形為2.06×10-3，礦體開(kāi)采采空區(qū)周邊礦柱僅有少量單元發(fā)生塑性破壞，塑性區(qū)變化主要出現(xiàn)在采空區(qū)頂板圍巖。

圖12 塑性區(qū)云圖Fig.12 Cloud image of plastic zone

4 結(jié)論

利用DIMINE強(qiáng)大的前處理功能，建立了南川河石灰?guī)r礦的三維地質(zhì)模型，導(dǎo)入到MIDAS/GTS，形成了包含采空區(qū)、礦體、地表地形的三維有限元數(shù)值模型。利用基于DIMINE-MIDAS/GTS的耦合技術(shù)，分析了礦體開(kāi)采后礦山的位移變形、應(yīng)力集中現(xiàn)象及塑性區(qū)破壞情況，分析論證該礦山采空區(qū)的穩(wěn)定性，結(jié)果表明：

1)+100 m中段采空區(qū)頂板產(chǎn)生的最大位移沉降為-2.4 mm，采空區(qū)的底板與巷道底板發(fā)生最大正向位移+2.94 mm。采空區(qū)形成的地表沉降有微弱的影響，沉降值不足1 mm，表明該礦山開(kāi)采對(duì)地表影響較小。

2)采空區(qū)頂板發(fā)生了拉應(yīng)力破壞。+130 m中段老采空區(qū)最大拉應(yīng)力為2.27 MPa，接近圍巖的抗拉強(qiáng)度值，采空區(qū)容易發(fā)生冒頂片幫事故。+100 m中段采空區(qū)與13號(hào)礦柱接觸處產(chǎn)生了最大壓應(yīng)力值集中12.66 MPa，礦柱處于壓應(yīng)力破壞。

3)采空區(qū)開(kāi)采后產(chǎn)生的塑性區(qū)變形主要分布在采空區(qū)的四周邊界，產(chǎn)生的塑性區(qū)變形范圍占總面積的44.5%，其中最大塑性區(qū)變形量為2.06×10-3，采空區(qū)周邊礦柱僅有少量單元發(fā)生塑性破壞，塑性區(qū)主要存在于采空區(qū)頂板。

4)綜合考慮巖體位移、應(yīng)力、塑性形變等計(jì)算結(jié)果，礦體開(kāi)采采空區(qū)頂板產(chǎn)生少量的豎直位移，地表發(fā)生微量的沉降，部分采空區(qū)頂板發(fā)生塑性變形，只在+130 m中段老采空區(qū)發(fā)生拉應(yīng)力破壞，表明礦山采空區(qū)基本處于穩(wěn)定狀態(tài)。