王永增 王心泉 王 潤 馮 春 曹 洋 孟磊磊 吳恩澤
(1.鞍鋼集團(tuán)礦業(yè)有限公司齊大山分公司,遼寧鞍山114000;2.中國科學(xué)院力學(xué)研究所,北京100190;3.中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京100049)
地下礦產(chǎn)資源被開采出來后會遺留下采空區(qū),采空區(qū)的存在改變了巖體的原有應(yīng)力場,使得巖體的原有平衡被破壞[1-3],易產(chǎn)生地表塌陷、地面沉降等問題,對露天開采的安全性存在巨大的威脅[4-6]。采空區(qū)頂板的臨界厚度是開采過程安全性評價(jià)的重要指標(biāo),對采空區(qū)穩(wěn)定性評價(jià)有著重要的意義。近年來,不少學(xué)者對此進(jìn)行了深入研究。唐碩等[7]基于模糊物源理論,確定采空區(qū)穩(wěn)定性的影響因素,建立了采空區(qū)穩(wěn)定性評價(jià)模型,評價(jià)結(jié)果用關(guān)聯(lián)度表示。柳小波等[8]結(jié)合第一強(qiáng)度理論提出了頂板厚度計(jì)算的新方法,采用Abaqus模擬分析了各因素對采空區(qū)頂板破壞的影響規(guī)律。姜桂鵬等[9]基于量綱分析法建立了采空區(qū)頂板穩(wěn)定性判別的數(shù)學(xué)模型,結(jié)合有限元數(shù)值模擬分析,得到了穩(wěn)定的臨界厚跨比。唐勝利等[10]基于BP神經(jīng)網(wǎng)絡(luò)原理,結(jié)合采空區(qū)穩(wěn)定性影響因素構(gòu)建了BP神經(jīng)網(wǎng)絡(luò)模型,并對空洞型采空區(qū)進(jìn)行了穩(wěn)定性評價(jià)。王炳文等[11]綜合理論分析和數(shù)值模擬結(jié)果,對遺留采空區(qū)群進(jìn)行了穩(wěn)定性級別評定。張瑞等[12]利用三維探測掃描設(shè)備建立了采空區(qū)三維數(shù)值模型,綜合多種采空區(qū)穩(wěn)定性評價(jià)方法對采空區(qū)進(jìn)行安全分級。盧欣奇等[13]通過FLAC軟件對近地表采空區(qū)進(jìn)行數(shù)值模擬分析,提出了治理方案,并通過FLAC模擬分析對各方案的治理效果進(jìn)行了對比。張敏思等[14]采用多種理論方法計(jì)算了不同跨度采空區(qū)頂板的安全厚度,并采用RFPA數(shù)值模擬軟件對采空區(qū)進(jìn)行模擬得到頂板出現(xiàn)初始損傷及失穩(wěn)垮塌的厚度。周平等[15]依據(jù)尖點(diǎn)突變模型得到采空區(qū)在充填體荷載作用下的頂板厚度表達(dá)式。王志修等[16]采用理論計(jì)算和數(shù)值模擬手段研究了盧安夏露天銅礦地下采空區(qū)頂板的安全厚度。胡洪旺等[17]基于Ressiner厚板理論構(gòu)建了理論分析模型,推導(dǎo)出頂板撓曲線方程及頂板最大應(yīng)力表達(dá)式,并利用Matlab軟件計(jì)算得到頂板厚度隨坐標(biāo)軸變化的三維曲線圖。
現(xiàn)有的采空區(qū)頂板穩(wěn)定性研究理論中,多將采空區(qū)頂板視為彈性梁進(jìn)行近似處理,難以考慮巖石的材料性質(zhì)對采空區(qū)穩(wěn)定性的影響。同時多數(shù)研究只通過數(shù)值計(jì)算對采空區(qū)頂板進(jìn)行規(guī)律性分析,缺少對采空區(qū)頂板臨界厚度的定量分析。本研究從量綱分析出發(fā),探討與采空區(qū)臨界頂板厚度相關(guān)的因素,借助連續(xù)-非連續(xù)單元法(Continuous-Discontinuous Element Method,CDEM)通過數(shù)值分析,找出影響臨界頂板厚度的關(guān)鍵變量,得到臨界厚度與相關(guān)變量之間對應(yīng)公式,為評價(jià)采空區(qū)穩(wěn)定性提供定量性的參考指標(biāo)。
采空區(qū)臨界頂板厚度的影響因素復(fù)雜,根據(jù)以往工程經(jīng)驗(yàn)及采空區(qū)研究資料,綜合考慮各方面的作用,本研究將影響因素分為幾何形狀因素、地層巖性因素及外力荷載因素[18-19]。
影響采空區(qū)臨界頂板厚度Hc的自變量為:①幾何參數(shù)為采空區(qū)長度a、寬度b和高度h;②材料參數(shù)有密度ρ、彈性模量E、泊松比μ、黏聚力c、內(nèi)摩擦角?、抗拉強(qiáng)度σt和剪脹角ψ;③靜載荷F(假設(shè)載荷作用點(diǎn)位于采空區(qū)正中);④重力加速度g(圖1)。
臨界頂板厚度Hc與各影響因素之間的關(guān)系可表述為
以密度ρ、重力加速度g、采空區(qū)長度a為基本量,建立的無量綱表達(dá)式為
理論分析可知:彈性模量、泊松比對彈性變形有較大的影響,但對臨界頂板厚度Hc影響不大,因此可以忽略;此外,采空區(qū)高度h僅對頂板垮落后的下沉距離有影響,對臨界頂板厚度Hc無影響,因此也可以忽略。由此,式(2)可變?yōu)?/p>
本研究采用連續(xù)-非連續(xù)單元法(CDEM)建立二維采空區(qū)數(shù)值模型,重點(diǎn)探討采空區(qū)的幾何參數(shù)及巖體物理力學(xué)參數(shù)對采空區(qū)穩(wěn)定性的影響規(guī)律,并借助安全系數(shù)進(jìn)行定量分析。由于采空區(qū)臨界頂板厚度的影響因素眾多,本研究采用多參數(shù)分析理論與方法[20-21]對參數(shù)進(jìn)行靈敏度分析,進(jìn)而給出頂板安全系數(shù)的主要影響因素。
計(jì)算時根據(jù)采空區(qū)的常見特征設(shè)立一組參數(shù)基礎(chǔ)值,當(dāng)研究某一個參數(shù)的影響時,其他參數(shù)取基礎(chǔ)值,各參數(shù)基礎(chǔ)值取值如表1所示。
根據(jù)工程中常見的采空區(qū)幾何尺寸及采空區(qū)頂板力學(xué)參數(shù)的取值范圍,確定采空區(qū)特征尺寸及采空區(qū)頂板力學(xué)參數(shù)的下限值及上限值。根據(jù)以往研究成果,表2中的參數(shù)對采空區(qū)安全系數(shù)的影響都是單調(diào)的,故首先研究極限參數(shù)下采空區(qū)頂板的安全系數(shù),根據(jù)安全系數(shù)變化值分析參數(shù)的敏感性。采空區(qū)參數(shù)限值及對應(yīng)的安全系數(shù)如表2所示。
由表2可知:黏聚力對安全系數(shù)的影響最大,采空區(qū)跨度、巖體內(nèi)摩擦角、巖體密度次之。因此,本研究重點(diǎn)研究不同頂板厚度下,采空區(qū)跨度、黏聚力、內(nèi)摩擦角、巖體密度與空區(qū)頂板安全系數(shù)的對應(yīng)關(guān)系。當(dāng)研究某一參數(shù)對頂板穩(wěn)定性的影響時,剩余參數(shù)取值為基礎(chǔ)值。頂板穩(wěn)定性規(guī)律研究時,各影響因素的計(jì)算參數(shù)取值范圍在設(shè)定的基礎(chǔ)值附近,不超過理論及工程上的合理取值范圍,試驗(yàn)中各參數(shù)取值如表3所示。
在不同的采空區(qū)頂板厚度H下,頂板安全系數(shù)與采空區(qū)跨度的對應(yīng)關(guān)系如圖2所示。分析該圖可知:隨著采空區(qū)跨度增大,安全系數(shù)逐漸減小,但減小趨勢逐漸變緩;隨著頂板厚度增大,安全系數(shù)逐漸增大,但增大幅度也逐漸變緩。
巖體密度對采空區(qū)頂板穩(wěn)定性的影響如圖3所示。由圖3可知:隨著巖體密度增大,頂板安全系數(shù)基本呈線性減小趨勢;且隨著頂板厚度增大,安全系數(shù)逐漸增大,但增大幅度逐漸變緩。
巖體黏聚力對采空區(qū)頂板穩(wěn)定性的影響規(guī)律如圖4所示。由圖4可知:隨著巖體黏聚力增大,頂板安全系數(shù)基本呈線性增大趨勢;且隨著頂板厚度增大,安全系數(shù)與黏聚力之間直線關(guān)系的斜率逐漸增大。
巖體內(nèi)摩擦角對采空區(qū)頂板穩(wěn)定性的影響規(guī)律如圖5所示。由圖5可知:隨著巖體內(nèi)摩擦角增大,頂板安全系數(shù)呈現(xiàn)非線性增大趨勢;內(nèi)摩擦角小于65°時,基本呈直線變化;大于65°時,安全系數(shù)有突然增大現(xiàn)象??傮w上,隨著頂板厚度增大,安全系數(shù)逐漸增大,但增大趨勢逐漸變緩。
本研究借助連續(xù)-非連續(xù)單元法(CDEM)并采用二維模型及三維模型分別進(jìn)行采空區(qū)臨界頂板厚度分析。二維模型采用連續(xù)介質(zhì)模型進(jìn)行計(jì)算,單元模型為Mohr-Coulomb理想彈塑性本構(gòu)模型;三維模型采用非連續(xù)介質(zhì)模型進(jìn)行計(jì)算,單元模型為線彈性本構(gòu)模型,單元間的虛擬界面采用脆性斷裂的Mohr-Coulomb本構(gòu)模型。
首先探討采空區(qū)長度及寬度的影響,建立的三維模型外邊界尺寸為140 m×140 m×64 m(長×寬×高),采空區(qū)高度為24 m,本研究假設(shè)模型長度a和寬度b一致,如圖6所示。單元模型采用線彈性模型,單元間的邊界采用脆性斷裂模型。計(jì)算過程中的基礎(chǔ)參數(shù)取值如表1所示。
無量綱臨界厚度與無量綱強(qiáng)度間的關(guān)系如圖7所示。由圖7可知:隨著無量綱強(qiáng)度增加,無量綱臨界厚度迅速減小,兩者關(guān)系可用下式進(jìn)行表示:
本研究構(gòu)建了長100 m、高50 m的二維數(shù)值模型,研究無量綱臨界高度與無量綱強(qiáng)度間的對應(yīng)關(guān)系。由于模擬的是二維平面應(yīng)變問題,因此寬度b為無限大。本研究采用理想彈塑性模型進(jìn)行分析,計(jì)算過程中的基礎(chǔ)參數(shù)取值見表1。
共設(shè)計(jì)3類工況,工況1:設(shè)定臨界高度Hc為5 m,采空區(qū)高度h為5 m,采空區(qū)長度a分別為10、20、30、40、50 m,通過二分法求取對應(yīng)的強(qiáng)度值;工況2:設(shè)定臨界高度Hc為10 m,采空區(qū)高度h為5 m,取采空區(qū)長度a分別為10、20、30、40、50 m,通過二分法求取相應(yīng)的強(qiáng)度值;工況3:設(shè)定臨界高度Hc為15 m,采空區(qū)高度h為5 m,取采空區(qū)長度a分別為20、30、40、50 m,通過二分法求取相應(yīng)的強(qiáng)度值。
第1類工況下,臨界失穩(wěn)時的塑性體應(yīng)變云圖如圖8所示。由圖8可知,當(dāng)Hc/a較大時,采空區(qū)頂板主要發(fā)生沿著側(cè)壁切落;隨著Hc/a減小,采空區(qū)頂板側(cè)壁易發(fā)生拉剪破壞,采空區(qū)中部也易發(fā)生彎拉 破壞。
基于上述3種工況的模擬結(jié)果,得到了無量綱臨界高度與無量綱強(qiáng)度之間的關(guān)系曲線,如圖9所示。由圖9可知:盡管3種工況下,各類參數(shù)差異較大,但在無量綱的表述下,均可用一條曲線進(jìn)行表述,采用指數(shù)衰減型函數(shù)進(jìn)行擬合,得到的無量綱臨界高度與無量綱強(qiáng)度的關(guān)系式為
綜合分析式(4)、式(5),可得無量綱高度與無量綱強(qiáng)度之間的關(guān)系可表示為
其中,α、β、γ為待定系數(shù)。
需要指出的是,式(4)與式(5)的系數(shù)有所差別。一方面是由于兩個模型采用了不同的計(jì)算模型及計(jì)算本構(gòu)模型;另一方面是由于二維模型忽略了采空區(qū)第3個方向尺度對采空區(qū)頂板穩(wěn)定性的影響。
(1)借助量綱分析手段,分析了采空區(qū)臨界頂板穩(wěn)定性的主要影響因素,得到與采空區(qū)頂板厚度相關(guān)的無量綱參數(shù)。選取黏聚力、采空區(qū)跨度、巖體內(nèi)摩擦角、巖體密度為參數(shù),借助CDEM數(shù)值分析方法,通過改變各參數(shù)的取值計(jì)算了采空區(qū)頂板的安全系數(shù)。并對二維、三維兩種模型下無量綱臨界高度與無量綱強(qiáng)度間的對應(yīng)關(guān)系進(jìn)行了研究,得到了無量綱臨界高度與無量綱強(qiáng)度間的函數(shù)表達(dá)式。
(2)彈性模量、泊松比對臨界頂板厚度影響不大,采空區(qū)高度對臨界頂板厚度無影響;黏聚力對安全系數(shù)影響最大,采空區(qū)跨度、巖體內(nèi)摩擦角、巖體密度次之。在只改變某一參數(shù)取值的情況下進(jìn)行分析發(fā)現(xiàn),隨著采空區(qū)跨度增加、巖體密度增加、黏聚力減小、內(nèi)摩擦角減小,頂板安全系數(shù)呈減小趨勢;隨著頂板厚度增加,安全系數(shù)增加,但增加趨勢變緩;隨著無量綱強(qiáng)度增大,無量綱的臨界頂板高度基本呈指數(shù)衰減型的下降趨勢。
(3)在進(jìn)行二維及三維無量綱臨界高度與無量綱強(qiáng)度關(guān)系表達(dá)式計(jì)算中,由于兩種模型所采用的本構(gòu)模型有所差別,且二維模型忽略了采空區(qū)第3個方向尺度的影響,導(dǎo)致兩種模型雖然關(guān)系表達(dá)式一致,但在擬合公式的系數(shù)取值上存在一定的差異。在下一步工作中,可考慮改變模型本構(gòu),或?qū)ν荒P透淖冇?jì)算參數(shù),研究表達(dá)式系數(shù)是否發(fā)生變化,再做進(jìn)一步的研究。