基于DIMINE-MIDAS/GTS的采空區(qū)穩(wěn)定性分析

任青云

(湖南有色金屬研究院，長沙 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)定性及地表位移沉降的影響范圍和程度。

1 工程概況

南川河石灰?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)計表(部分)

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

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

為了合理評價南川河石灰?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é)果匯總表

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

由于自然巖體內(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ù)匯總表

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

3.1 三維礦體模型的建立

制約礦山地表位移沉降及采空區(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

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

根據(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 數(shù)值模擬計算結(jié)果分析

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

4 結(jié)論

利用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)。