馬宗源,焦 貝,黨發(fā)寧,成玉祥

(1.貴州交通職業(yè)技術(shù)學(xué)院,貴州 貴陽 551400;2.西安理工大學(xué)巖土工程研究所,陜西 西安 710048;3.長安大學(xué)地質(zhì)工程與測繪學(xué)院,陜西 西安 710054)

0 引言

我國是煤炭開采量較大的國家,許多地區(qū)由于煤層埋深較淺,地下采煤誘發(fā)地表變形并造成嚴(yán)重塌陷。近年來，由于地下煤層開采誘發(fā)的黃土邊坡破壞及滑坡失穩(wěn)情況日趨嚴(yán)重,嚴(yán)重影響了煤礦生產(chǎn)安全及當(dāng)?shù)孛癖姷纳睢?991年8月,陜西銅川金華山煤礦西側(cè)黃土梁邊由于采煤引起大規(guī)模的崩塌性黃土滑坡,土體達(dá)1 050萬m3,滑坡體將坡腳處的西龍村埋沒,摧毀大片農(nóng)田。2004年4月4日1時20分,由于多個煤礦開挖導(dǎo)致陜西彬縣城關(guān)鎮(zhèn)下溝村發(fā)生3次大面積山體滑坡,滑坡體黃土涌入水簾小河,將河流截斷,數(shù)百名群眾緊急撤離。2005年甘肅華亭煤礦區(qū)硯峽村受采空區(qū)塌陷影響誘發(fā)了大型滑坡,體積達(dá)395萬m3,致使硯峽村部分房屋出現(xiàn)裂縫及倒塌,學(xué)校、醫(yī)院、信用社等公用設(shè)施遭到不同程度的破壞,直接經(jīng)濟損失達(dá)2 000多萬元。2010年7月,陜西省延安市某礦區(qū)發(fā)生大規(guī)模黃土滑坡,滑坡體積達(dá)10萬m3,滑坡及坡體變形如圖1所示。從現(xiàn)場滑體特征可以看出,滑坡速度極快,滑體沖向?qū)γ嫘逼?所幸未造成人員傷亡。2011年7月,呂梁市方山縣山西焦煤集團(tuán)霍州煤電有限責(zé)任公司方山店坪煤礦發(fā)生一起黃土滑坡,滑體掩埋多個建筑及設(shè)備,造成5名人員死亡。由此可見,地下采煤引發(fā)的滑坡災(zāi)害已經(jīng)成為煤礦開采區(qū)除地面塌陷災(zāi)害外最嚴(yán)重的地質(zhì)災(zāi)害形式[1]。根據(jù)張銀洲等[2]的研究,地下采煤形成采空區(qū),當(dāng)開采煤層的采深介于30～100 m時,就會產(chǎn)生較為強烈的地表變形,從而誘發(fā)黃土邊坡失穩(wěn)。

圖1 礦區(qū)采礦誘發(fā)黃土滑坡及開裂Fig.1 Loess landslide and cracking in slope triggered by underground coal mining

目前國內(nèi)外對地下開采誘發(fā)的滑坡已經(jīng)進(jìn)行了大量研究,但缺乏地下采煤誘發(fā)邊坡失穩(wěn)的有效分析方法及評價標(biāo)準(zhǔn)?，F(xiàn)有研究大多是基于靜力學(xué)理論進(jìn)行分析[3-9]，而地下采煤對地表邊坡的影響是一個動態(tài)變化的過程,開采過程中對地表邊坡的擾動作用屬于動荷載范疇,因此需使用動力學(xué)的理論及方法對地下采煤誘發(fā)邊坡失穩(wěn)問題進(jìn)行計算分析。

本文以陜西省延安市黃陵礦區(qū)某煤礦發(fā)生的典型黃土滑坡為研究對象,使用顯式有限元大變形計算方法分析不同煤層開挖高度情況下,開采擾動對地表黃土邊坡安全系數(shù)的動態(tài)影響,提出一種適應(yīng)于地下采煤誘發(fā)黃土邊坡失穩(wěn)及漸進(jìn)破壞過程的分析方法,以期為黃土地區(qū)地下采礦誘發(fā)的邊坡及滑坡災(zāi)害防治工程提供一定的理論指導(dǎo)。

1 研究方法及理論

1.1 土動力學(xué)基本理論

動三軸試驗數(shù)據(jù)說明Hardin-Drnevich骨干加載曲線模型(下文簡稱HD模型)更適于預(yù)測飽和黃土的動應(yīng)力-應(yīng)變關(guān)系[10],因此采用土動力學(xué)模型分析地下采煤過程中開挖擾動對黃土邊坡的影響,選取HD模型預(yù)測土的動模量和阻尼比隨動應(yīng)變的衰減關(guān)系。HD模型中土的動剪切模量和阻尼比可寫為動剪應(yīng)變的函數(shù)[11-12]:

(1)

(2)

式中:Gmax為土體動剪切模量的最大值;Dmax為土體的最大阻尼比;γ為土體的動剪應(yīng)變幅值;γref為參考剪應(yīng)變，γref=τmax/Gmax,即土體所受最大剪應(yīng)力與剪切模量最大值的比值,可換算為G/Gmax=0.5時動剪切模量對應(yīng)的動剪應(yīng)變幅值γ,γref越大,土的動剪切模量隨動應(yīng)變衰減越快,土的阻尼越大。

基于HD模型建立黃土的動應(yīng)力-動應(yīng)變本構(gòu)模型,由式(1)可推導(dǎo)出剪應(yīng)力與剪應(yīng)變的關(guān)系:

(3)

動力情況下,土體處于加、卸載循環(huán)狀態(tài),卸載時剪應(yīng)力可寫為:

(4)

再加載時有:

(5)

式中:γc為加、卸循環(huán)荷載轉(zhuǎn)換時刻土體的剪應(yīng)變;τc為加、卸循環(huán)荷載轉(zhuǎn)換時刻土體的剪應(yīng)力,τc可寫為:

(6)

圖2為不同土的動剪切模量、阻尼比與動剪應(yīng)變的關(guān)系以及HD模型計算結(jié)果對比,其中砂土和黏土數(shù)據(jù)取自文獻(xiàn)[13]和文獻(xiàn)[14],陜北黃土數(shù)據(jù)取自文獻(xiàn)[15]。采用Visual Fortran語言編制黃土動本構(gòu)關(guān)系計算程序,并使用VUMAT子程序接口將模型導(dǎo)入ABAQUS有限元計算軟件。用四節(jié)點減縮積分單元計算剪應(yīng)力及剪應(yīng)變關(guān)系,并與黃土動三軸試驗測試數(shù)據(jù)進(jìn)行對比。對該單元的頂部約束水平自由度,單元底部水平方向上加載正弦加速度時程曲線,正弦時程函數(shù)如式(7)所示,峰值加速度為0.04g:

(7)

黃土動力學(xué)參數(shù)如下:剪切模量最大值Gmax=75 MPa,泊松比ν=0.25,參考剪應(yīng)變γref=0.03。黃土的動應(yīng)力與動應(yīng)變關(guān)系理論預(yù)測與試驗數(shù)據(jù)的對比如圖3所示,可以看出根據(jù)HD模型建立的黃土動力本構(gòu)模型能夠很好地模擬黃土加卸載過程中呈滯回曲線形狀的應(yīng)力-應(yīng)變關(guān)系。

圖3 黃土動應(yīng)力與動應(yīng)變關(guān)系理論預(yù)測與試驗 數(shù)據(jù)的對比Fig.3 Comparison between theoretical prediction and experimental data of the relationship between dynamic stress and dynamic strain of loess

1.2 強度折減法計算邊坡安全系數(shù)

使用強度折減法計算邊坡的安全系數(shù)FOS(Factor of Safety)。邊坡抗剪強度參數(shù)(黏聚力及內(nèi)摩擦角)折減方式可寫為[16]:

(8)

式中:SRF(Strength Reduction Factor)為折減系數(shù);c和φ分別為邊坡土體原始的黏聚力及內(nèi)摩擦角;cf和φf分別為經(jīng)過折減之后的邊坡土體的黏聚力及內(nèi)摩擦角。將邊坡位移突然增大時刻(折減系數(shù)與邊坡最大位移關(guān)系曲線拐點)的折減系數(shù)作為邊坡安全系數(shù)。

以文獻(xiàn)[16]中無下臥地層邊坡算例作為對比,分別使用顯式及隱式有限元方法計算邊坡的安全系數(shù),其中顯式有限元方法采用大變形,隱式有限元方法采用小變形模式進(jìn)行計算。算例邊坡坡高10 m,坡比為1∶2。土體計算參數(shù)為:彈性模量E=100 MPa,泊松比υ=0.25,密度γ=2 000 kg/m3,黏聚力為c=10 kPa,內(nèi)摩擦角φ=20°,不考慮剪脹角的影響,重力加速度g=9.81 m/s2。圖4為本文使用顯式有限元方法及文獻(xiàn)[16]使用隱式有限元方法通過強度折減計算出的邊坡最大位移結(jié)果的對比。由圖4可以看出去,使用顯式有限元和隱式有限元方法都可以通過強度折減系數(shù)與邊坡最大位移關(guān)系曲線的突變確定邊坡安全系數(shù),但由于顯式有限元方法采用大變形模式,因此邊坡位移比隱式方法要大。圖5為顯式及隱式有限元方法得出的網(wǎng)格變形圖,可看出二者差別不大。顯式有限元方法采用顯式時間積分方案求解動力學(xué)方程,材料及幾何非線性問題不需要進(jìn)行迭代,因此更適用于動力學(xué)及大變形問題的計算分析[17-18]。

圖4 通過強度折減計算出的邊坡最大位移Fig.4 Maximum displacement of slope using the strength reduction method

圖5 顯式及隱式方法計算得出的邊坡極限狀態(tài)時的網(wǎng)格變形圖(變形放大15倍)Fig.5 Mesh deformation of slope in limit state calculated by explicit and implicit methods (Scale factor is 15)

2 數(shù)值模擬研究

2.1 計算模型及邊界條件

將實際問題簡化為二維平面應(yīng)變下均質(zhì)土邊坡問題,使用ABAQUS有限元軟件的顯式動力分析模塊及大變形模式進(jìn)行分析研究。根據(jù)實際工程概況,建立黃土邊坡下煤層開挖模型,模型中黃土邊坡的坡比為1∶2,坡高為30 m,計算模型尺寸及網(wǎng)格剖分情況如圖6所示。其中,H為煤層開挖頂板埋深,h為煤層開挖高度,煤層向邊坡方向開挖。使用二維四邊形減縮積分單元進(jìn)行網(wǎng)格劃分,共劃分14 076個單元,28 658個節(jié)點。由于將采煤擾動過程視為動力學(xué)問題進(jìn)行計算分析,模型邊界條件采用動力邊界條件進(jìn)行設(shè)置，將模型底部及兩側(cè)邊界設(shè)置為零加速度邊界,以保證計算域內(nèi)產(chǎn)生的加速度、速度及位移以波的形式傳播至邊界處時不會反射回計算域內(nèi)。但該邊界條件會導(dǎo)致邊界節(jié)點的位移及速度變化也為零,所以要盡量遠(yuǎn)離計算分析區(qū)域設(shè)置,以免影響位移及應(yīng)力分布。煤層開挖頂板埋深H為100 m,開挖速度為10 m/天,在煤層開挖高度h為2.5 m、5 m、7.5 m的情況下分別計算邊坡的穩(wěn)定性。計算分析中黃土邊坡土體強度準(zhǔn)則使用Mohr-Coulomb準(zhǔn)則,黃土體物理力學(xué)參數(shù)列于表1。

表1 計算參數(shù)匯總表Table 1 Summary of calculation parameters

圖6 地下采煤邊坡計算模型尺寸及網(wǎng)格 劃分示意圖(單位：m)Fig.6 Dimension and meshing of the slope model with underground mining (Unit:m)

2.2 計算結(jié)果及分析

圖7、圖8及圖9分別為不同煤層開挖高度情況下,采煤過程中邊坡及地層的位移云圖。從采煤過程中邊坡及地層的位移及邊坡安全系數(shù)變化可以發(fā)現(xiàn),邊坡安全系數(shù)隨開挖面與邊坡坡腳距離L的減少而逐漸降低(負(fù)值代表挖過邊坡坡腳的距離),同時邊坡位移逐漸增大,邊坡穩(wěn)定性逐漸降低。

圖7 不同開采位置邊坡及地層的位移云圖及安全系數(shù)(h=2.5 m)Fig.7 Displacement contours of slope and FOS with different mining positions (h=2.5 m)

圖8 不同開采位置時邊坡的位移云圖及安全系數(shù)(h=5 m)Fig.8 Displacement contours of slope and FOS with different mining positions (h=5 m)

圖9 不同開采位置時邊坡的位移云圖及安全系數(shù)(h=7.5 m)Fig.9 Displacement contours of slope and FOS with different mining positions (h=7.5 m)

圖10為折減系數(shù)SRF=1.33時,地下采煤開挖過程中邊坡坡頂水平向位移的變化情況,可以看出邊坡位移隨著開挖進(jìn)尺的增加逐漸增大,在接近坡腳正下方時邊坡坡體滑動破壞。圖11為不同煤層開挖高度下邊坡安全系數(shù)變化曲線,可以看出開挖前邊坡處于靜力情況下的安全系數(shù)(FOS=1.41),隨著采煤開挖進(jìn)尺的不斷加大,邊坡安全系數(shù)逐漸降低。此外,還可以看出地下采煤過程中邊坡安全系數(shù)是動態(tài)變化的,因此地下采煤誘發(fā)邊坡失穩(wěn)也是一個漸進(jìn)破壞的過程。不同煤層開挖高度下的計算結(jié)果(圖11)顯示,邊坡安全系數(shù)變化的拐點相似,都是與坡腳水平距離40 m處,而煤層開挖高度7.5 m情況下邊坡安全系數(shù)數(shù)值降幅明顯,說明煤層開挖高度超過5 m將對地表邊坡穩(wěn)定性產(chǎn)生顯著影響。

圖10 地下采煤過程中邊坡坡頂水平向位移變化Fig.10 Horizontal displacement variation of slope crest during underground coal mining

圖11 地下采煤過程中邊坡安全系數(shù)變化Fig.11 Variation of slope SOF during underground mining

3 結(jié)論

本文采用顯式有限元及大變形計算方法分析了地下采煤對地表黃土邊坡的穩(wěn)定性影響問題,首先使用黃土動三軸試驗結(jié)果對黃土動力本構(gòu)模型進(jìn)行驗證,再對顯式和隱式有限元計算方法進(jìn)行對比分析,最后對地下采煤情況下地表邊坡穩(wěn)定性問題進(jìn)行計算分析，并得出如下主要結(jié)論:

(1) 本文探討了一種新的邊坡穩(wěn)定性分析方法,用來分析地下采煤等動態(tài)擾動情況下地表邊坡的穩(wěn)定性,該方法計算效率高且能夠分析邊坡的漸進(jìn)破壞過程。

(2) 通過本文分析研究可知,地下采煤對地表邊坡穩(wěn)定性的影響是一個動態(tài)的過程,地表邊坡失穩(wěn)也是一個漸進(jìn)破壞的過程,因此必須使用動力學(xué)方法及土動力學(xué)理論分析該類型邊坡的穩(wěn)定性問題。

(3) 在本文研究成果的基礎(chǔ)上,下一步將考慮三維、煤層頂板塌落及地表裂縫等因素的影響,開展進(jìn)一步研究。