基于CDEM的采空區(qū)穩(wěn)定性及臨界頂板厚度分析

王永增 王心泉 王 潤 馮 春 曹 洋 孟磊磊 吳恩澤

（1.鞍鋼集團礦業(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ū)頂板的臨界厚度是開采過程安全性評價的重要指標，對采空區(qū)穩(wěn)定性評價有著重要的意義。近年來，不少學(xué)者對此進行了深入研究。唐碩等［7］基于模糊物源理論，確定采空區(qū)穩(wěn)定性的影響因素，建立了采空區(qū)穩(wěn)定性評價模型，評價結(jié)果用關(guān)聯(lián)度表示。柳小波等［8］結(jié)合第一強度理論提出了頂板厚度計算的新方法，采用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ū)進行了穩(wěn)定性評價。王炳文等［11］綜合理論分析和數(shù)值模擬結(jié)果，對遺留采空區(qū)群進行了穩(wěn)定性級別評定。張瑞等［12］利用三維探測掃描設(shè)備建立了采空區(qū)三維數(shù)值模型，綜合多種采空區(qū)穩(wěn)定性評價方法對采空區(qū)進行安全分級。盧欣奇等［13］通過FLAC軟件對近地表采空區(qū)進行數(shù)值模擬分析，提出了治理方案，并通過FLAC模擬分析對各方案的治理效果進行了對比。張敏思等［14］采用多種理論方法計算了不同跨度采空區(qū)頂板的安全厚度，并采用RFPA數(shù)值模擬軟件對采空區(qū)進行模擬得到頂板出現(xiàn)初始損傷及失穩(wěn)垮塌的厚度。周平等［15］依據(jù)尖點突變模型得到采空區(qū)在充填體荷載作用下的頂板厚度表達式。王志修等［16］采用理論計算和數(shù)值模擬手段研究了盧安夏露天銅礦地下采空區(qū)頂板的安全厚度。胡洪旺等［17］基于Ressiner厚板理論構(gòu)建了理論分析模型，推導(dǎo)出頂板撓曲線方程及頂板最大應(yīng)力表達式，并利用Matlab軟件計算得到頂板厚度隨坐標軸變化的三維曲線圖。

現(xiàn)有的采空區(qū)頂板穩(wěn)定性研究理論中，多將采空區(qū)頂板視為彈性梁進行近似處理，難以考慮巖石的材料性質(zhì)對采空區(qū)穩(wěn)定性的影響。同時多數(shù)研究只通過數(shù)值計算對采空區(qū)頂板進行規(guī)律性分析，缺少對采空區(qū)頂板臨界厚度的定量分析。本研究從量綱分析出發(fā)，探討與采空區(qū)臨界頂板厚度相關(guān)的因素，借助連續(xù)-非連續(xù)單元法（Continuous-Discontinuous Element Method，CDEM）通過數(shù)值分析，找出影響臨界頂板厚度的關(guān)鍵變量，得到臨界厚度與相關(guān)變量之間對應(yīng)公式，為評價采空區(qū)穩(wěn)定性提供定量性的參考指標。

1 量綱分析

1.1 采空區(qū)臨界頂板厚度影響因素分析

采空區(qū)臨界頂板厚度的影響因素復(fù)雜，根據(jù)以往工程經(jīng)驗及采空區(qū)研究資料，綜合考慮各方面的作用，本研究將影響因素分為幾何形狀因素、地層巖性因素及外力荷載因素［18-19］。

影響采空區(qū)臨界頂板厚度Hc的自變量為：①幾何參數(shù)為采空區(qū)長度a、寬度b和高度h；②材料參數(shù)有密度ρ、彈性模量E、泊松比μ、黏聚力c、內(nèi)摩擦角?、抗拉強度σt和剪脹角ψ；③靜載荷F（假設(shè)載荷作用點位于采空區(qū)正中）；④重力加速度g（圖1）。

1.2 無量綱表達式建立

臨界頂板厚度Hc與各影響因素之間的關(guān)系可表述為

以密度ρ、重力加速度g、采空區(qū)長度a為基本量，建立的無量綱表達式為

理論分析可知：彈性模量、泊松比對彈性變形有較大的影響，但對臨界頂板厚度Hc影響不大，因此可以忽略；此外，采空區(qū)高度h僅對頂板垮落后的下沉距離有影響，對臨界頂板厚度Hc無影響，因此也可以忽略。由此，式（2）可變?yōu)?/p>

2 采空區(qū)穩(wěn)定性影響因素

2.1 頂板安全系數(shù)主要影響因素

本研究采用連續(xù)-非連續(xù)單元法（CDEM）建立二維采空區(qū)數(shù)值模型，重點探討采空區(qū)的幾何參數(shù)及巖體物理力學(xué)參數(shù)對采空區(qū)穩(wěn)定性的影響規(guī)律，并借助安全系數(shù)進行定量分析。由于采空區(qū)臨界頂板厚度的影響因素眾多，本研究采用多參數(shù)分析理論與方法［20-21］對參數(shù)進行靈敏度分析，進而給出頂板安全系數(shù)的主要影響因素。

計算時根據(jù)采空區(qū)的常見特征設(shè)立一組參數(shù)基礎(chǔ)值，當研究某一個參數(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)摩擦角、巖體密度次之。因此，本研究重點研究不同頂板厚度下，采空區(qū)跨度、黏聚力、內(nèi)摩擦角、巖體密度與空區(qū)頂板安全系數(shù)的對應(yīng)關(guān)系。當研究某一參數(shù)對頂板穩(wěn)定性的影響時，剩余參數(shù)取值為基礎(chǔ)值。頂板穩(wěn)定性規(guī)律研究時，各影響因素的計算參數(shù)取值范圍在設(shè)定的基礎(chǔ)值附近，不超過理論及工程上的合理取值范圍，試驗中各參數(shù)取值如表3所示。

2.2 頂板安全系數(shù)變化規(guī)律

在不同的采空區(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ù)逐漸增大，但增大趨勢逐漸變緩。

3 采空區(qū)臨界頂板厚度分析

本研究借助連續(xù)-非連續(xù)單元法（CDEM）并采用二維模型及三維模型分別進行采空區(qū)臨界頂板厚度分析。二維模型采用連續(xù)介質(zhì)模型進行計算，單元模型為Mohr-Coulomb理想彈塑性本構(gòu)模型；三維模型采用非連續(xù)介質(zhì)模型進行計算，單元模型為線彈性本構(gòu)模型，單元間的虛擬界面采用脆性斷裂的Mohr-Coulomb本構(gòu)模型。

3.1 三維模型下無量綱參數(shù)間關(guān)系

首先探討采空區(qū)長度及寬度的影響，建立的三維模型外邊界尺寸為140 m×140 m×64 m（長×寬×高），采空區(qū)高度為24 m，本研究假設(shè)模型長度a和寬度b一致，如圖6所示。單元模型采用線彈性模型，單元間的邊界采用脆性斷裂模型。計算過程中的基礎(chǔ)參數(shù)取值如表1所示。

無量綱臨界厚度與無量綱強度間的關(guān)系如圖7所示。由圖7可知：隨著無量綱強度增加，無量綱臨界厚度迅速減小，兩者關(guān)系可用下式進行表示：

3.2 二維模型下無量綱參數(shù)間關(guān)系

本研究構(gòu)建了長100 m、高50 m的二維數(shù)值模型，研究無量綱臨界高度與無量綱強度間的對應(yīng)關(guān)系。由于模擬的是二維平面應(yīng)變問題，因此寬度b為無限大。本研究采用理想彈塑性模型進行分析，計算過程中的基礎(chǔ)參數(shù)取值見表1。

共設(shè)計3類工況，工況1：設(shè)定臨界高度Hc為5 m，采空區(qū)高度h為5 m，采空區(qū)長度a分別為10、20、30、40、50 m，通過二分法求取對應(yīng)的強度值；工況2：設(shè)定臨界高度Hc為10 m，采空區(qū)高度h為5 m，取采空區(qū)長度a分別為10、20、30、40、50 m，通過二分法求取相應(yīng)的強度值；工況3：設(shè)定臨界高度Hc為15 m，采空區(qū)高度h為5 m，取采空區(qū)長度a分別為20、30、40、50 m，通過二分法求取相應(yīng)的強度值。

第1類工況下，臨界失穩(wěn)時的塑性體應(yīng)變云圖如圖8所示。由圖8可知，當Hc/a較大時，采空區(qū)頂板主要發(fā)生沿著側(cè)壁切落；隨著Hc/a減小，采空區(qū)頂板側(cè)壁易發(fā)生拉剪破壞，采空區(qū)中部也易發(fā)生彎拉 破壞。

基于上述3種工況的模擬結(jié)果，得到了無量綱臨界高度與無量綱強度之間的關(guān)系曲線，如圖9所示。由圖9可知：盡管3種工況下，各類參數(shù)差異較大，但在無量綱的表述下，均可用一條曲線進行表述，采用指數(shù)衰減型函數(shù)進行擬合，得到的無量綱臨界高度與無量綱強度的關(guān)系式為

綜合分析式（4）、式（5），可得無量綱高度與無量綱強度之間的關(guān)系可表示為

其中，α、β、γ為待定系數(shù)。

需要指出的是，式（4）與式（5）的系數(shù)有所差別。一方面是由于兩個模型采用了不同的計算模型及計算本構(gòu)模型；另一方面是由于二維模型忽略了采空區(qū)第3個方向尺度對采空區(qū)頂板穩(wěn)定性的影響。

4 結(jié)論

（1）借助量綱分析手段，分析了采空區(qū)臨界頂板穩(wěn)定性的主要影響因素，得到與采空區(qū)頂板厚度相關(guān)的無量綱參數(shù)。選取黏聚力、采空區(qū)跨度、巖體內(nèi)摩擦角、巖體密度為參數(shù)，借助CDEM數(shù)值分析方法，通過改變各參數(shù)的取值計算了采空區(qū)頂板的安全系數(shù)。并對二維、三維兩種模型下無量綱臨界高度與無量綱強度間的對應(yīng)關(guān)系進行了研究，得到了無量綱臨界高度與無量綱強度間的函數(shù)表達式。

（2）彈性模量、泊松比對臨界頂板厚度影響不大，采空區(qū)高度對臨界頂板厚度無影響；黏聚力對安全系數(shù)影響最大，采空區(qū)跨度、巖體內(nèi)摩擦角、巖體密度次之。在只改變某一參數(shù)取值的情況下進行分析發(fā)現(xiàn)，隨著采空區(qū)跨度增加、巖體密度增加、黏聚力減小、內(nèi)摩擦角減小，頂板安全系數(shù)呈減小趨勢；隨著頂板厚度增加，安全系數(shù)增加，但增加趨勢變緩；隨著無量綱強度增大，無量綱的臨界頂板高度基本呈指數(shù)衰減型的下降趨勢。

（3）在進行二維及三維無量綱臨界高度與無量綱強度關(guān)系表達式計算中，由于兩種模型所采用的本構(gòu)模型有所差別，且二維模型忽略了采空區(qū)第3個方向尺度的影響，導(dǎo)致兩種模型雖然關(guān)系表達式一致，但在擬合公式的系數(shù)取值上存在一定的差異。在下一步工作中，可考慮改變模型本構(gòu)，或?qū)ν荒Ｐ透淖冇嬎銋?shù)，研究表達式系數(shù)是否發(fā)生變化，再做進一步的研究。