任青云
(湖南有色金屬研究院,長沙 410100)
地下礦山開采的重點安全問題是開采形成的采空區(qū)對后續(xù)開采產(chǎn)生的巖層移動問題,開采過程中伴隨的地層沉降及地壓活動[1-2]。礦山開采面臨的重要安全問題主要表現(xiàn)在,開采區(qū)域的不斷擴大,地層隨之進(jìn)行持續(xù)的運動,由于開采現(xiàn)場會出現(xiàn)不同的地質(zhì)條件,例如溶洞、斷層、空區(qū)等可能會產(chǎn)生移動變形的大范圍空間,此時會造成地下礦山的地壓顯現(xiàn),導(dǎo)致頂板出現(xiàn)大變形,繼而影響圍巖發(fā)生冒落和片幫等現(xiàn)象。因此地下礦山長期開采導(dǎo)致的采空區(qū)穩(wěn)定性是影響礦山生命財產(chǎn)安全及后續(xù)開采的重要因素[3-4]。
影響采空區(qū)的穩(wěn)定性及礦山地表位移沉降的主要原因有開采礦體的埋深、采場的高度,礦體的厚度、傾角、巖石力學(xué)特性及礦山采用的采礦工藝。采用大型有限元數(shù)值模擬[5-6]前處理模式可以直觀地表示出地下采空區(qū)的分布情況及整改礦山的輪廓,其后處理模式能夠計算出礦山的地表位移沉降、采空區(qū)是否發(fā)生拉應(yīng)力破壞、礦柱是否產(chǎn)生主應(yīng)力集中及塑性區(qū)變形破壞情況。本文根據(jù)南川河石灰?guī)r礦礦區(qū)地質(zhì)資料,結(jié)合礦體情況及開采設(shè)計,對南川河石灰?guī)r礦進(jìn)行取樣試驗,采用數(shù)值模擬方法分析礦山開采后巖層移動變形情況,并結(jié)合傳統(tǒng)的理論計算方法對礦柱、頂板穩(wěn)定性進(jìn)行分析,以此確定無支護(hù)狀態(tài)下采空區(qū)跨度,從而評判采空區(qū)的穩(wěn)定性及地表位移沉降的影響范圍和程度。
南川河石灰?guī)r礦位于瀏陽市澄潭江鎮(zhèn)南川河畔,距瀏陽市城區(qū)直距 21 km,公司始建于1972年,原采用露天開采,1996年開始地下開采,主要采用斜坡道開拓,對角式通風(fēng)。礦山長期使用的采礦方法是房柱采礦法。礦山總的回采順序是沿礦體傾斜由上而下循階段進(jìn)行,按邊探邊采的原則進(jìn)行開采。采場回采順序是后退式開采,礦房內(nèi)采用分層開采。礦山2014年之前開采130 m中段,2014年至今開采100 m中段,目前礦山正在+100 m中段北部進(jìn)行生產(chǎn)。
由于石灰?guī)r的巖性較軟,并且礦山開采歷史悠久,地下空區(qū)存在極大的不穩(wěn)定性及不確定性,為了掌握空區(qū)對下部礦體回采的影響,必須針對采空區(qū)進(jìn)行穩(wěn)定性研究和分析。根據(jù)礦山的實際生產(chǎn)現(xiàn)狀,針對空區(qū)進(jìn)行了現(xiàn)場調(diào)查,并通過統(tǒng)計分析獲得礦山不同開采水平的采空區(qū)規(guī)模、尺寸。此次現(xiàn)場調(diào)查針對+130 m中段和+100 m中段共統(tǒng)計出130個采空區(qū),+100 m中段采空區(qū)部分調(diào)查數(shù)據(jù)見表1?,F(xiàn)場調(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)計表(部分)
為了合理評價南川河石灰?guī)r礦工程地質(zhì)條件,為后續(xù)數(shù)值模擬提供依據(jù),開展了礦巖力學(xué)參數(shù)室內(nèi)試驗,測試的內(nèi)容包括礦巖的容重、單軸抗壓強度、抗拉強度、彈性模量、泊松比等基本物理力學(xué)參數(shù)。針對礦山實際情況,選取了石灰?guī)r巖樣進(jìn)行加工,并對各巖樣進(jìn)行了相關(guān)參數(shù)測定。巖樣總計22個,其中5個用于開展單軸壓縮試驗、5個用于巴西劈裂試驗、12個用于變角剪試驗。
根據(jù)數(shù)值模擬計算分析需要的基本數(shù)據(jù),主要進(jìn)行了如下四項內(nèi)容的實驗室試驗:1)礦、巖的天然容重;2)礦、巖的單軸抗壓力學(xué)特性;3)巖體的劈裂抗拉強度試驗;4)巖體的泊松比與彈性模量試驗。
試驗開始前,先對巖樣加工。首先,按照標(biāo)準(zhǔn)巖石力學(xué)試驗所要求的規(guī)格,制成試樣尺寸為Φ50 mm×100 mm的標(biāo)準(zhǔn)巖石試樣;然后擦干試樣,在自然狀態(tài)下風(fēng)干;最后,試驗在250 t全數(shù)字型液壓伺服剛性巖石力學(xué)試驗系統(tǒng)(RMT150,圖1)及30 t萬能材料試驗機上完成。
圖1 RMT150型試驗設(shè)備Fig.1 RMT150 test equipment
2.1.1 單軸抗壓試驗
試樣尺寸為直徑Ф50 mm,高徑比約為2∶1,巖石單軸抗壓試驗結(jié)果如表2所示。
表2 巖石單軸抗壓強度及靜力受壓彈性模量和泊松比測試
2.1.2 抗拉強度試驗
巖石的抗拉強度試驗如圖2所示,各種試驗結(jié)果見表3和表4。
圖2 巴西劈裂試驗Fig.2 Brazilian splitting experiment
表3 各試樣劈裂抗拉試驗參數(shù)測定值
表4 變角剪試驗結(jié)果
根據(jù)單軸抗壓試驗和單軸抗拉試驗,將室內(nèi)試驗結(jié)果匯總?cè)绫?所示。
表5 礦巖力學(xué)參數(shù)試驗結(jié)果匯總表
由于自然巖體內(nèi)部有很多缺陷,比如節(jié)理、裂隙和孔洞等,因此自然巖體與小試樣巖石的力學(xué)性質(zhì)是有差異的,導(dǎo)致實驗室測得小試樣巖石力學(xué)參數(shù)不能直接作為數(shù)值模擬的力學(xué)參數(shù)。因此,巖體力學(xué)工程界,大多數(shù)依據(jù)室內(nèi)巖石力學(xué)實驗,采用理論力學(xué)的方法,估算出巖體的力學(xué)參數(shù),用于有限元數(shù)值模擬計算。
根據(jù)HOEK,CARRANZA-TORRES提出的巖體破壞經(jīng)驗準(zhǔn)則:
(1)
巖體的mb、S和a由GSI和上式計算的mi確定。
(2)
(3)
(4)
根據(jù)確定的mb、S和a值,可以計算出巖體的單軸抗壓強度和抗拉強度:
σc=σci·sa
(5)
(6)
巖體的內(nèi)摩擦角φ和黏聚力c可以由以下公式進(jìn)行計算:
(7)
(8)
針對南川河石灰?guī)r礦工程地質(zhì)條件,取石灰?guī)rGSI=70,根據(jù)室內(nèi)巖石力學(xué)參數(shù)試驗結(jié)果,利用HOEK和CARRANZA-TORRES提出的巖體破壞經(jīng)驗準(zhǔn)則及巖體分類RMR值,采用Roclab軟件進(jìn)行計算。選取巖體的力學(xué)參數(shù)如表6所示。
表6 巖體物理力學(xué)參數(shù)匯總表
制約礦山地表位移沉降及采空區(qū)穩(wěn)定性的重要因素,包括礦區(qū)地表的地形地質(zhì)條件及井下采空區(qū)的分布狀態(tài),因此在三維有限元數(shù)值模擬計算過程中,缺少完整的地形模塊,不能真實地反映礦山的實際情況。本研究通過基于DIMINE強度的前處理功能,利用二維地質(zhì)剖面圖或者M(jìn)APGIS地形數(shù)據(jù)庫生成三維地質(zhì)模型,對勘探線剖面進(jìn)行初步的處理,僅保留勘探線與礦體解譯線[8]。將礦山剖面圖經(jīng)過CAD處理之后導(dǎo)入到DIMINE中,進(jìn)行線串運算使勘探線“立起來”,對立起來的三維線文件進(jìn)行坐標(biāo)運算,使各點坐標(biāo)為其真實坐標(biāo)值。通過上述基本工作的完成,可以聯(lián)接各礦體解譯線剖面,對于比較復(fù)雜的礦體解譯線,可以構(gòu)造中線,通過中線與剖面線來生成實體,礦體實體模型見圖3。在DIMINE中生成的礦山實體模型,無法直接用于有限元數(shù)值模擬的運算,需要將此模型導(dǎo)入到MIDAS/GTS中,利用MIDAS/GTS的網(wǎng)格劃分功能,結(jié)合采礦工藝,進(jìn)一步建立采區(qū)模型,劃分網(wǎng)格單元,進(jìn)行數(shù)值模擬分析。
圖3 南川河石灰?guī)r礦三維實體模型Fig.3 3D solid model of Nanchuanhe limestone mine
根據(jù)礦山的地層參數(shù),確定有限元數(shù)值模擬的巖體力學(xué)參數(shù)見表6。計算過程礦巖及廢石均采用摩爾-庫倫((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—抗拉強度,MPa。MIDAS/GTS有限元計算程序,默認(rèn)當(dāng)fs大于0時,表明巖體已發(fā)生了剪切破壞;當(dāng)ft大于0時,表明巖體已發(fā)生了拉伸破壞。
為了模型建立與分析結(jié)果的準(zhǔn)確性,本數(shù)值模擬進(jìn)行了以下假設(shè):模型的礦巖體為均質(zhì)、各向同性材料,模擬的三維網(wǎng)格模型底部為固定約束,四周為限制水平約束,只在豎直方向有位移變形。初始階段考慮地應(yīng)力與自重的影響,位移清零,表明礦山在原始狀態(tài)已處于穩(wěn)定的固結(jié)沉降。計算模擬不考慮地震、爆破等動荷載的影響,不分析地下水的影響;礦山地質(zhì)條件較復(fù)雜,礦巖穩(wěn)固性在空間分布上具有較大的隨機性,礦體為斷層分割,且礦巖和圍巖中均有穿插,這些地質(zhì)現(xiàn)象的影響在巖體參數(shù)折減時已做考慮,在模擬過程中不予另行考慮。
根據(jù)地下采空區(qū)對礦山圍巖的影響范圍,設(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ì)量情況見圖4和圖5。
圖4 三維網(wǎng)格模型+邊界條件Fig.4 3D grid model + boundary conditions
圖5 網(wǎng)格質(zhì)量檢查Fig.5 Grid quality inspection
3.3.1 位移場分析
采空區(qū)的穩(wěn)定性分析,位移是評判礦山穩(wěn)定性的一個重要指標(biāo),因為大量的地表變形,會使礦山發(fā)生地表塌陷,且位移變形量超過巖體的變形限值,會發(fā)生采空區(qū)頂板冒頂片幫事故。從圖6~9的位移云圖得出,該礦山產(chǎn)生的位移變化情況,在礦體的頂板位移值小于0,說明由于產(chǎn)生了采空區(qū),發(fā)生了頂板下沉變形情況,但隨著礦體的開挖過程進(jìn)行,在采空區(qū)的底部及巷道的底板,豎直方向位移大于0,是因為礦體開挖深度加大,采空區(qū)及巷道周圍發(fā)生了主應(yīng)力集中現(xiàn)象,水平方向的應(yīng)力大于垂直方向的應(yīng)力,使巷道底板產(chǎn)生上拱的位移變形。
從圖6豎向位移云圖可以看出,+130 m中段礦體開采時,采空區(qū)頂板豎向位移最大值為-1.74 mm,礦體開采底板巷道發(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中段礦體開采后采空區(qū)頂板產(chǎn)生的位移沉降變形云圖,目前在+100 m中段回采時采空區(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ū)地表沉降等值線范圍見圖9,位移變形值數(shù)據(jù)見表7。
圖9 地表沉降位移云圖與沉降等值線圖Fig.9 Surface displacement cloud map and subsidence contour map
表7 采空區(qū)位移變化值計算結(jié)果
3.3.2 應(yīng)力場分析
礦體開采后對圍巖造成的影響是一個非線性的不可逆的加載過程,它使得處于初始地應(yīng)力狀態(tài)下的圍巖進(jìn)行應(yīng)力重新分布,最后達(dá)到新的平衡,礦體開挖引起徑向應(yīng)力釋放、切向應(yīng)力增加,導(dǎo)致礦體開采頂板和底板產(chǎn)生拉應(yīng)力;從應(yīng)力云圖可以看出礦體開采區(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ù)計算結(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,已接近巖體的抗拉強度,表明在采空區(qū)的頂板部分已發(fā)生拉應(yīng)力破壞,容易產(chǎn)生冒頂片幫事故。從最大應(yīng)力云圖可以得出:在礦柱與空區(qū)的邊界圍巖出現(xiàn)壓應(yīng)力集中現(xiàn)象,最大壓應(yīng)力集中發(fā)生在礦山+100 m中段東北方向的采空區(qū)與13號礦柱交界處,產(chǎn)生的最大壓應(yīng)力值為12.66 MPa,+130 m中段的(老采空區(qū))礦柱發(fā)生最大壓應(yīng)力集中值為8.88 MPa,表明礦柱發(fā)生受壓破壞。應(yīng)力結(jié)果見表8。
表8 最大、最小主應(yīng)力值計算結(jié)果
3.3.3 塑性區(qū)分析
圖12為分析得出的采空區(qū)塑性變形情況,從計算結(jié)果可以得出,在采空區(qū)的附近發(fā)生了塑性區(qū)變形破壞,塑性區(qū)面積占總面積約44.5%,最大塑性變形為2.06×10-3,礦體開采采空區(qū)周邊礦柱僅有少量單元發(fā)生塑性破壞,塑性區(qū)變化主要出現(xiàn)在采空區(qū)頂板圍巖。
圖12 塑性區(qū)云圖Fig.12 Cloud image of plastic zone
利用DIMINE強大的前處理功能,建立了南川河石灰?guī)r礦的三維地質(zhì)模型,導(dǎo)入到MIDAS/GTS,形成了包含采空區(qū)、礦體、地表地形的三維有限元數(shù)值模型。利用基于DIMINE-MIDAS/GTS的耦合技術(shù),分析了礦體開采后礦山的位移變形、應(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,表明該礦山開采對地表影響較小。
2)采空區(qū)頂板發(fā)生了拉應(yīng)力破壞。+130 m中段老采空區(qū)最大拉應(yīng)力為2.27 MPa,接近圍巖的抗拉強度值,采空區(qū)容易發(fā)生冒頂片幫事故。+100 m中段采空區(qū)與13號礦柱接觸處產(chǎn)生了最大壓應(yīng)力值集中12.66 MPa,礦柱處于壓應(yīng)力破壞。
3)采空區(qū)開采后產(chǎn)生的塑性區(qū)變形主要分布在采空區(qū)的四周邊界,產(chǎn)生的塑性區(qū)變形范圍占總面積的44.5%,其中最大塑性區(qū)變形量為2.06×10-3,采空區(qū)周邊礦柱僅有少量單元發(fā)生塑性破壞,塑性區(qū)主要存在于采空區(qū)頂板。
4)綜合考慮巖體位移、應(yīng)力、塑性形變等計算結(jié)果,礦體開采采空區(qū)頂板產(chǎn)生少量的豎直位移,地表發(fā)生微量的沉降,部分采空區(qū)頂板發(fā)生塑性變形,只在+130 m中段老采空區(qū)發(fā)生拉應(yīng)力破壞,表明礦山采空區(qū)基本處于穩(wěn)定狀態(tài)。