基于數(shù)值模擬正交試驗(yàn)的采空區(qū)穩(wěn)定性因素的敏感性分析

何榮興，韓智勇，劉 洋，張 晶，鄒開(kāi)華

(1.東北大學(xué)資源與土木工程學(xué)院，遼寧 沈陽(yáng) 110819；2.內(nèi)蒙古科技大學(xué)礦業(yè)研究院，內(nèi)蒙古 包頭 014010)

地下礦產(chǎn)資源采出后形成采空區(qū)，采空區(qū)周?chē)鷰r體的力學(xué)平衡被打破，在重力和次生應(yīng)力作用下，會(huì)產(chǎn)生裂隙和變形。采空區(qū)的長(zhǎng)期存在使一些礦區(qū)的地壓活動(dòng)加劇，相繼從局部崩落發(fā)展到大范圍的突然崩落，引起地表大范圍的巖移，造成了嚴(yán)重的破壞。目前，采空區(qū)已成為金屬非金屬礦山安全生產(chǎn)的主要危險(xiǎn)源[1]，因此準(zhǔn)確評(píng)估采空區(qū)穩(wěn)定性狀態(tài)，對(duì)于確定采空區(qū)處理方法和處理時(shí)機(jī)具有重要的指導(dǎo)意義。

目前，采空區(qū)穩(wěn)定性評(píng)價(jià)方法主要有定性分析法、定量分析法、現(xiàn)場(chǎng)監(jiān)測(cè)分析法等。采空區(qū)變形破壞總是在若干特定影響因素作用下產(chǎn)生的，而每一個(gè)影響因素又由若干因子組成，采空區(qū)穩(wěn)定性受多種因素制約，各因素間又相互關(guān)聯(lián)，決定了采空區(qū)穩(wěn)定性評(píng)判是一個(gè)復(fù)雜的系統(tǒng)工程。因此，越來(lái)越多的學(xué)者采用數(shù)學(xué)方法對(duì)采空區(qū)穩(wěn)定性進(jìn)行綜合評(píng)價(jià)[2-9]。在建立評(píng)價(jià)體系中，采空區(qū)穩(wěn)定性因素的影響權(quán)重確定方法主要是專家打分法和1～9標(biāo)度法，但二者均帶有一定的主觀性，影響采空區(qū)的評(píng)價(jià)結(jié)果。確定各影響因素的敏感性對(duì)確定復(fù)雜采空區(qū)失穩(wěn)主控因素及相應(yīng)的臨界值并設(shè)計(jì)采場(chǎng)結(jié)構(gòu)參數(shù)、為采空區(qū)穩(wěn)定性綜合評(píng)判體系構(gòu)建過(guò)程中權(quán)重確定具有十分重要意義。國(guó)內(nèi)外學(xué)者在此領(lǐng)域開(kāi)展了大量研究工作。宋衛(wèi)東等[10]從礦柱載荷、強(qiáng)度、失穩(wěn)形式及影響因素四個(gè)方面推導(dǎo)出2種礦柱的安全系數(shù)計(jì)算公式，采用6因素、5水平正交試驗(yàn)分析礦柱穩(wěn)定性影響因素的敏感性，并研究主要影響因素與礦柱安全系數(shù)之間的關(guān)系。 尹升華等[11]建立計(jì)算矩形礦柱安全系數(shù)的簡(jiǎn)化公式，采用正交極差分析對(duì)礦柱穩(wěn)定性影響因素的敏感性進(jìn)行評(píng)價(jià)，得到礦房寬度、礦柱寬度與礦柱安全系數(shù)間的關(guān)系曲線。隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬在多方案理論研究方面具有明顯優(yōu)勢(shì)，數(shù)值模擬已成為研究采空區(qū)穩(wěn)定性的重要手段之一。王曉軍等[12]利用FLAC3D軟件完成多因素組合影響階段礦柱上采過(guò)程數(shù)值模擬正交試驗(yàn)，分析了單一影響因素與頂板臨界厚度的關(guān)系，利用多元非線性回歸的數(shù)學(xué)方法，建立了上采過(guò)程3因素組合影響下頂板臨界厚度數(shù)學(xué)預(yù)測(cè)公式。何標(biāo)慶[13]利用FLAC3D軟件對(duì)福建省某礦區(qū)大型采空區(qū)群的穩(wěn)定性進(jìn)行了模擬，結(jié)果表明采空區(qū)間柱多發(fā)生塑性破壞，采空區(qū)上下盤(pán)圍巖存在冒落風(fēng)險(xiǎn)。彭超等[14]利用FLAC3D軟件對(duì)某大理巖礦山采空區(qū)穩(wěn)定性進(jìn)行了模擬分析，結(jié)果表明采空區(qū)塑性破壞主要發(fā)生在頂板淺層圍巖，采空區(qū)不會(huì)發(fā)生大規(guī)模失穩(wěn)，并根據(jù)不同區(qū)域的穩(wěn)定性制定相應(yīng)的處理方案。

鑒于采空區(qū)穩(wěn)定性受多因素影響的特點(diǎn)，本文綜合各類巖體分級(jí)標(biāo)準(zhǔn)評(píng)價(jià)指標(biāo)，能夠反映巖體力學(xué)參數(shù)且容易量化的影響采空區(qū)穩(wěn)定性因素，采用正交試驗(yàn)確定研究方案，利用霍克布朗準(zhǔn)則，將各方案中的巖體質(zhì)量指標(biāo)轉(zhuǎn)化了巖體力學(xué)參數(shù)，利用FLAC3D軟件模擬分析影響采空區(qū)穩(wěn)定性各因素的敏感性，對(duì)確定采空區(qū)穩(wěn)定性評(píng)價(jià)的權(quán)重及礦山采場(chǎng)結(jié)構(gòu)參數(shù)提供參考。

1 采空區(qū)穩(wěn)定性影響因素分析及選取

影響采空區(qū)穩(wěn)定性的因素較多，一般可概括為巖體質(zhì)量因素、采空區(qū)規(guī)格參數(shù)和誘發(fā)因素三類。巖體質(zhì)量是影響巖體穩(wěn)定性的最基本因素，目前國(guó)際巖體分級(jí)標(biāo)準(zhǔn)包括Q分級(jí)[15-16]、RMR分級(jí)[17-18]、MRMR分級(jí)[19-20]以及《工程巖體分級(jí)標(biāo)準(zhǔn)》(GB/T 50218—2014)，各分級(jí)標(biāo)準(zhǔn)所需的巖體質(zhì)量因素見(jiàn)表1。不同分級(jí)標(biāo)準(zhǔn)的側(cè)重點(diǎn)和選擇巖體質(zhì)量因素不盡相同，因此，本文選定的分析采空區(qū)穩(wěn)定性影響因素能基本反映巖體質(zhì)量，容易量化，而且通過(guò)這些因素可直接或間接估算出巖體的力學(xué)參數(shù)。點(diǎn)荷載強(qiáng)度或單軸抗壓強(qiáng)度是反映巖體質(zhì)量的重要指標(biāo)，是大部分分級(jí)標(biāo)準(zhǔn)都不可缺少的，而且也是估算巖體力學(xué)參數(shù)的基礎(chǔ)，二者具有良好的相關(guān)性，考慮到單軸抗壓強(qiáng)度更加直接，因此選擇巖石單軸抗壓強(qiáng)度作為采空區(qū)穩(wěn)定性的分析指標(biāo)。巖石質(zhì)量指標(biāo)RQD、節(jié)理間距、巖體完整性系數(shù)、節(jié)理面產(chǎn)狀、節(jié)理面條件等表征結(jié)構(gòu)面特征的因素，存在相關(guān)性或是同一參數(shù)在不同方面的反映程度，根據(jù)上述原則選擇巖石質(zhì)量指標(biāo)RQD、節(jié)理間距作為采空區(qū)穩(wěn)定性因素進(jìn)行分析。地下水對(duì)巖體工程穩(wěn)定性和巖體本身強(qiáng)度具有重要影響，因此地下水因素在采空區(qū)穩(wěn)定性分析中是不能忽略的。在已知節(jié)理面條件和產(chǎn)狀的條件下，通過(guò)選定的上述因素即可獲得巖體的RMR值即反映巖體質(zhì)量，再利用Hoek-Brown強(qiáng)度準(zhǔn)則便可估算出巖體的抗拉強(qiáng)度、內(nèi)聚力、內(nèi)摩擦角、彈性模量等力學(xué)參數(shù)。在模擬計(jì)算中還需要泊松比來(lái)確定水平應(yīng)力大小和體積模量等參數(shù)，因此最終選擇的影響采空區(qū)穩(wěn)定性巖體質(zhì)量因素包括單軸抗壓強(qiáng)度、RQD、節(jié)理面間距、地下水、泊松比等5個(gè)參數(shù)。

表1 不同分級(jí)標(biāo)準(zhǔn)包含的巖體質(zhì)量因素Table 1 Rock mass quality factors included in different grading standards

采空區(qū)穩(wěn)定性很大程度上受采空區(qū)本身規(guī)格參數(shù)影響，主要包括采空區(qū)的高度、面積、埋藏深度、形狀。生產(chǎn)實(shí)踐中形成的采空區(qū)形狀較為復(fù)雜，很難量化，為了便于分析量化，本文簡(jiǎn)化為規(guī)則采空區(qū)，用采空區(qū)的長(zhǎng)寬比來(lái)定量表征采空區(qū)的形狀。因此，在采空區(qū)規(guī)格參數(shù)中選擇采空區(qū)的高度、面積、埋深、長(zhǎng)寬比作為采空區(qū)穩(wěn)定性評(píng)價(jià)因素。

實(shí)踐表明，眾多采空區(qū)災(zāi)害的發(fā)生，多是由于連續(xù)強(qiáng)降雨、山體滑坡、大規(guī)模爆破振動(dòng)、采動(dòng)擾動(dòng)影響以及相鄰采空區(qū)的影響等因素誘發(fā)的，可稱為誘發(fā)因素，這些因素和采空區(qū)本身的穩(wěn)定性無(wú)直接關(guān)系，且難于量化，因此本文不予考慮。

2 正交試驗(yàn)方案及力學(xué)參數(shù)估算

2.1 正交試驗(yàn)方案設(shè)計(jì)及數(shù)值模型建立

根據(jù)上述分析，本文選取9個(gè)因素，分別為單軸抗壓強(qiáng)度、RQD、節(jié)理面間距、地下水、泊松比、采空區(qū)高度、采空區(qū)面積、采空區(qū)埋深和采空區(qū)的長(zhǎng)寬比，并對(duì)各因素對(duì)空區(qū)穩(wěn)定性影響程度進(jìn)行分析。根據(jù)正交原則設(shè)計(jì)了“9因素4水平”的正交試驗(yàn)方案(表2)。利用FLAC3D軟件建立各方案采空區(qū)模型，由于后續(xù)正交分析中要比較各方案的塑性區(qū)體積，為使結(jié)果具有可比性，建立模型時(shí)盡量使各模型采空區(qū)周?chē)膰鷰r網(wǎng)格尺寸和密度保持一致，由于篇幅有限，只列出前4個(gè)方案的數(shù)值模型，如圖1所示。

表2 正交試驗(yàn)方案Table 2 Orthogonal test scheme

圖1 前四個(gè)方案模型圖Fig.1 The simulation model of the four test schemes

對(duì)模型施加邊界條件，模型底部采取鉸支固定邊界法，四周邊界約束其橫向位移，初始垂直應(yīng)力為自重應(yīng)力，根據(jù)模型尺寸以及方案中采空區(qū)的埋深，在模型頂部分別施加1.08 MPa、6.48 MPa、11.87 MPa和17.28 MPa的垂直應(yīng)力，水平應(yīng)力根據(jù)凈水壓力假設(shè)，計(jì)算公式見(jiàn)式(1)。

(1)

式中：σh、σv分別為圍巖水平應(yīng)力和垂直應(yīng)力；v圍巖泊松比。

本次數(shù)值計(jì)算將礦巖視作非線性彈塑性介質(zhì)，采用莫爾-庫(kù)倫準(zhǔn)則作為屈服準(zhǔn)則，其控制方程見(jiàn)式(2)。

(2)

式中：σ1、σ3分別為最大主應(yīng)力和最小主應(yīng)力；c、φ分別為材料黏聚力和內(nèi)摩擦角；F為破壞判斷系數(shù)，當(dāng)F≥0時(shí)，材料將發(fā)生剪切破壞。

另外，材料在拉應(yīng)力狀態(tài)下，采用抗拉破壞強(qiáng)度準(zhǔn)則。其力學(xué)模型見(jiàn)式(3)。

F=σt-Rt

(3)

式中：σt為材料所受拉應(yīng)力；Rt為材料抗拉強(qiáng)度。如果拉應(yīng)力超過(guò)材料抗拉強(qiáng)度(F≥0)，材料將發(fā)生拉破壞。

2.2 巖體質(zhì)量參數(shù)量化和力學(xué)參數(shù)估算

巖體質(zhì)量參數(shù)對(duì)采空區(qū)穩(wěn)定性的影響通過(guò)對(duì)圍巖力學(xué)參數(shù)反映出來(lái)，因此各方案的巖體質(zhì)量參數(shù)估算圍巖的力學(xué)參數(shù)，不僅可將巖體質(zhì)量參數(shù)量化，還可以為數(shù)值模擬采空區(qū)穩(wěn)定性提供基礎(chǔ)。目前，最常用的巖體參數(shù)估算方法是經(jīng)驗(yàn)公式法，被廣泛采用的方法主要有完整性系數(shù)(Kv)修正法、費(fèi)森科法[21]、格吉法[22]以及霍克-布朗法[23]。 眾多研究[24-25]認(rèn)為Hoek-Brown強(qiáng)度準(zhǔn)則比較全面地反映了巖體結(jié)構(gòu)等特征對(duì)巖體強(qiáng)度的影響，是發(fā)展最完善的方法。Hoek-Brown巖體破壞準(zhǔn)則表達(dá)式見(jiàn)式(4)。

(4)

式中：σ1、σ3分別為巖體破壞時(shí)的最大主應(yīng)力和最小主應(yīng)力；m為巖體的Hoek-Brown常數(shù)；s為巖體材料性質(zhì)的無(wú)量綱系數(shù)，m、s可用式(5)和式(6)估算。

(5)

(6)

式中：mi可根據(jù)不同類型巖體的Hoek-Brown常數(shù)表[20-21]來(lái)選取；RMR可根據(jù)巖石的強(qiáng)度、RQD值、節(jié)理間距以及地下水影響的綜合評(píng)價(jià)巖體分類指標(biāo)，其值可通過(guò)節(jié)理巖體的巖石力學(xué)分類(RMR)表計(jì)算得出[26-27]。

當(dāng)σ1=0時(shí)，由式(4)可以得到巖體的抗拉強(qiáng)度，見(jiàn)式(7)。

(7)

Hoek-Brown給出巖體彈性模量Em與地質(zhì)強(qiáng)度指標(biāo)GSI的關(guān)系見(jiàn)式(8)。

(8)

式(8)可以用巖體分類指標(biāo)RMR替代GSI值進(jìn)行彈性模量估算。

Mohr-Coulomb強(qiáng)度準(zhǔn)則認(rèn)為巖體強(qiáng)度主要與巖體的內(nèi)聚力c和內(nèi)摩擦角φ有關(guān)，而最大主應(yīng)力和最小主應(yīng)力σ1、σ3之間存在線性關(guān)系，見(jiàn)式(9)。

σ1=σmc+kσ3

(9)

當(dāng)0<σ3<0.25σc時(shí)，所做的三軸實(shí)驗(yàn)結(jié)果同時(shí)符合Hoek-Brown準(zhǔn)則與Mohr-Coulomb準(zhǔn)則，可由式(4)確定的不同σ1、σ3對(duì)式(9)回歸得到巖體抗壓強(qiáng)度σmc和k。

另外，巖體內(nèi)摩擦角和內(nèi)聚力可由式(10)和式(11)求解。

(10)

(11)

假設(shè)采空區(qū)圍巖為石灰?guī)r，節(jié)理面稍粗糙，寬度<1 mm，節(jié)理面巖石軟弱，根據(jù)不同方案的巖體質(zhì)量數(shù)據(jù)獲得相應(yīng)的RMR等相關(guān)參數(shù)，再利用Hoek-Brown準(zhǔn)則便可估算出相應(yīng)的力學(xué)參數(shù)，見(jiàn)表3。

表3 各試驗(yàn)方案的計(jì)算參數(shù)Table 3 Parameters of each test scheme

3 模擬結(jié)果分析

為了對(duì)表1的試驗(yàn)方案進(jìn)行評(píng)價(jià)，需要制定評(píng)估指標(biāo)。采空區(qū)圍巖的塑性變形，可表征采空區(qū)失穩(wěn)的開(kāi)始，塑性區(qū)總體積的大小可評(píng)估采空區(qū)最終的穩(wěn)定程度。不同采空區(qū)參數(shù)條件下，采空區(qū)頂板和邊壁的塑性變形程度是不同的，有些采空區(qū)邊壁巖體不斷片幫導(dǎo)致最終采空區(qū)失穩(wěn)，有些采空區(qū)失穩(wěn)形式是頂板直接冒落，為評(píng)估采空區(qū)參數(shù)和失穩(wěn)形式關(guān)系，頂板塑性區(qū)和邊壁巖體塑性區(qū)可作為模擬結(jié)果的評(píng)估指標(biāo)。因此，選擇采空區(qū)塑性區(qū)總體積，頂板塑性區(qū)體積和邊壁巖體塑性區(qū)體積作為各方案評(píng)估的指標(biāo)。利用Fish語(yǔ)言編輯的塑性區(qū)體積查詢小程序，可獲得各方案的塑性體積，結(jié)果見(jiàn)表4。

對(duì)表4中的采空區(qū)模擬結(jié)果進(jìn)行極差分析，結(jié)果見(jiàn)表5。其中，K1i、K2i和K3i分別為指標(biāo)采空區(qū)頂板塑性區(qū)體積、采空區(qū)邊壁塑性區(qū)體積和采空區(qū)圍巖塑性區(qū)總體積在水平i下指標(biāo)的偏差平均值(i=1、 2、 3、 4)，Rj為對(duì)應(yīng)三個(gè)指標(biāo)的極差(j=1、2、3)。

表4 各試驗(yàn)方案的計(jì)算結(jié)果Table 4 Statistics of calculation results of each test scheme

通過(guò)極差計(jì)算得到的采空區(qū)頂板塑性區(qū)體積的影響因素敏感性主次順序依次為：?jiǎn)屋S抗壓強(qiáng)度>采空區(qū)埋深>采空區(qū)面積>泊松比>涌水量>RQD>節(jié)理面間距>采空區(qū)長(zhǎng)寬比>采空區(qū)高度；采空區(qū)邊壁塑性區(qū)體積的影響因素敏感性主次順序依次為：采空區(qū)埋深>采空區(qū)長(zhǎng)寬比>單軸抗壓強(qiáng)度>采空區(qū)高度>采空區(qū)面積>RQD>節(jié)理面間距>泊松比>涌水量；采空區(qū)塑性區(qū)總體積的影響因素敏感性主次順序依次為：?jiǎn)屋S抗壓強(qiáng)度>采空區(qū)埋深>采空區(qū)面積>采空區(qū)高度>采空區(qū)長(zhǎng)寬比>泊松比>節(jié)理面間距>RQD>涌水量。綜合這3個(gè)指標(biāo)可知，采空區(qū)埋深、采空區(qū)面積、單軸抗壓強(qiáng)度、長(zhǎng)寬比均排在前列，為進(jìn)一步評(píng)估這4個(gè)因素隨水平變化對(duì)采空區(qū)穩(wěn)定性的影響趨勢(shì)，將表5中各指標(biāo)偏差的平均值，繪制成與對(duì)采空區(qū)敏感性最為敏感的4個(gè)因素、水平的關(guān)系趨勢(shì)圖，如圖2～圖4所示。各指標(biāo)的偏差平均值大小與單軸抗壓強(qiáng)度、采空區(qū)面積、采空區(qū)埋深和采空區(qū)長(zhǎng)寬比因素水平變化呈曲線變化關(guān)系，采空區(qū)頂板塑性區(qū)體積偏差平均值隨水平變化呈快速降低-緩慢減小或緩慢升高-快速上升的變化關(guān)系，邊壁塑性區(qū)體積偏差平均值隨水平變化呈先增加后降低的變化關(guān)系，圍巖塑性區(qū)總體積偏差平均值隨水平先變化近似呈線性增加或降低變化關(guān)系?？傊?，采空區(qū)埋深、采空區(qū)面積、單軸抗壓強(qiáng)度、長(zhǎng)寬比4個(gè)因素是對(duì)采空區(qū)穩(wěn)定性影響最為敏感的，在評(píng)估采空區(qū)穩(wěn)定性和設(shè)計(jì)礦房尺寸時(shí)，應(yīng)將這4個(gè)因素作為重點(diǎn)因素。

續(xù)表4

表5 各試驗(yàn)方案的極差分析Table 5 Extract analysis results of each test scheme

圖2 頂板塑性區(qū)指標(biāo)偏差平均值隨因素和水平的變化趨勢(shì)Fig.2 Trend of the deviation average value of volume of roof plastic zone along with factors and levels

圖3 邊壁塑性區(qū)指標(biāo)偏差平均值隨因素和水平的變化趨勢(shì)Fig.3 Trend of the deviation average value of volume of sidewall plastic zone along with factors and levels

圖4 總塑性區(qū)指標(biāo)偏差平均值隨因素和水平的變化趨勢(shì)Fig.4 Trend of the deviation average value of volume of total plastic zone along with factors and levels

4 結(jié) 論

1) 在頂板塑性區(qū)體積，邊壁塑性區(qū)體積和塑性區(qū)總體積評(píng)估指標(biāo)下，各因素對(duì)采空區(qū)穩(wěn)定性影響的敏感性不盡相同，根據(jù)不同指標(biāo)下的因素敏感性順序，采空區(qū)埋深、采空區(qū)面積、單軸抗壓強(qiáng)度、長(zhǎng)寬比4個(gè)因素是對(duì)采空區(qū)穩(wěn)定性影響主要因素。

2) 采空區(qū)頂板塑性區(qū)體積偏差平均值隨水平變化呈快速降低-緩慢減小或緩慢升高-快速上升的變化關(guān)系，邊壁塑性區(qū)體積偏差平均值隨水平變化呈先增加后降低的變化關(guān)系，圍巖塑性區(qū)總體積偏差平均值隨水平先變化近似呈線性增加或降低變化關(guān)系。

3) 對(duì)于具體礦山而言，巖體質(zhì)量參數(shù)已固定，某一采深條件下的采場(chǎng)穩(wěn)定性主要取決于采空區(qū)面積和長(zhǎng)寬比，采場(chǎng)結(jié)構(gòu)參數(shù)和穩(wěn)定性也需要根據(jù)采深的變化進(jìn)行動(dòng)態(tài)評(píng)估。