基于動態(tài)強度折減DDA法的邊坡多滑面穩(wěn)定性分析

王述紅，朱承金，張紫杉，任藝鵬,王鵬宇，邱 偉

(東北大學 資源與土木工程學院，遼寧 沈陽 110819)

高邊坡由于具有規(guī)模大、地質(zhì)條件復(fù)雜、易受周邊環(huán)境干擾而破壞失穩(wěn)的特點，在初次失穩(wěn)后易發(fā)生二次失穩(wěn)，形成次級滑動面。巖質(zhì)高邊坡受層理結(jié)構(gòu)的影響，使得邊坡的穩(wěn)定問題更加復(fù)雜和突出，故對潛在多滑面的巖質(zhì)高邊坡穩(wěn)定性進行深入探究任重而道遠[1-3]。

邊坡穩(wěn)定性分析中有2種常用的計算方法:極限平衡法(LEM)[4]和有限元法(FEM)[5]。與有限元相比，極限平衡法計算效率高，通過幾何假定并基于已知或假定的規(guī)則滑面求解安全系數(shù)，但其忽視了邊坡滑落過程中巖土體的本構(gòu)關(guān)系[6]。近年來有限元快速發(fā)展并將強度折減法結(jié)合，該技術(shù)不必預(yù)先假設(shè)滑動面的位置，且易得安全系數(shù)及相應(yīng)潛在滑動面，其局限性在于安全系數(shù)及滑動面求解的單一性。CALA等[7]開創(chuàng)了新型強度折減法，在常規(guī)方法基礎(chǔ)上得出最危險失穩(wěn)面，進而確定首次滑落的步數(shù)Nr，繼續(xù)增大折減系數(shù)，隨即1.1Nr步將會被計算，相繼求得對應(yīng)次級潛在失穩(wěn)面及安全系數(shù)。李小春等[8]利用檢索邊坡體單元的方法，對不同折減系數(shù)下各范疇邊坡單元進行統(tǒng)籌，以此得到多個滑落面及相應(yīng)安全系數(shù)。

非連續(xù)性是巖體固有的屬性，由SHI[9]首次提出的非連續(xù)變形分析(DDA)在模擬塊體大位移、大變形方面獨具優(yōu)勢，在邊坡工程領(lǐng)域也得到迅速發(fā)展。MACLAUGHLIN等[10]模擬了傾斜邊坡的平面和弧形破壞模式，發(fā)現(xiàn)其結(jié)果比常規(guī)分析方法具有更好的精度。FU等[11]將矢量和方法(VSM)應(yīng)用到DDA中，基于實際應(yīng)力狀態(tài)和矢量和算法計算邊坡穩(wěn)定安全系數(shù)。

發(fā)展邊坡穩(wěn)定性評價核心問題在于安全系數(shù)及相應(yīng)滑落面的求解。過去研究中的強度折減法大都是對巖體統(tǒng)一折減，在實際情況下，主要是結(jié)構(gòu)面的力學性質(zhì)決定了巖體的力學性質(zhì)，且?guī)r體結(jié)構(gòu)面破壞時并不是同時損傷同等程度;再之邊坡首次發(fā)生失穩(wěn)后塊體位置及應(yīng)力已重新分布，故應(yīng)使用首次失穩(wěn)結(jié)束后的邊坡模型對次級滑落面進行分析，而以往對邊坡次級滑動面的研究中鮮有考慮該實際狀況。為了解決該問題，筆者借鑒DDA計算位移的優(yōu)勢，考慮塊體間相對位移，提出基于DSR-DDA法的邊坡多滑面搜索并進行邊坡穩(wěn)定性評價，利用傾斜平面滑塊經(jīng)典案例測試位移閾值，驗證方法的可行性及計算精度，并嘗試將其應(yīng)用到撫順西露天礦巖質(zhì)高邊坡的穩(wěn)定性分析中。

1 DSR-DDA法

1.1 DDA基本原理

DDA利用一階位移函數(shù)來表述每個塊體的運動參量，并假設(shè)應(yīng)力和應(yīng)變是恒定的[12]。每個塊體的位移矢量包含6個變量:

D=(u0,v0,r0,εx,εy,γxy)T

(1)

其中，(u0,v0)為塊體內(nèi)特定點(x0,y0)的剛性位移;r0為塊體繞特定點(x0,y0)的旋轉(zhuǎn)角度;εx,εy,γxy分別為該塊體的正應(yīng)變與切應(yīng)變。塊體內(nèi)任意點(x,y)的位移(u,v)為

(2)

基于最小勢能原理構(gòu)建的系統(tǒng)整體平衡方程為

(3)

式中，系數(shù)矩陣中Kij為6×6子矩陣;Kii由塊體單元的材料屬性和幾何參數(shù)決定;Kij(i≠j)則由塊體i和塊體j間的接觸條件而決定;[Di]和[Fi]為6×1子矩陣;Di為塊體i的變形變量(d1i,d2i,d3i,d4i,d5i,d6i);Fi為塊體i上分配給6個變形變量的荷載。

1.2 DSR-DDA法位移閾值

為了考慮真實情況中塊體結(jié)構(gòu)面損傷的程度不一性，提出對位移變化不小于閾值的塊體結(jié)構(gòu)面剪切強度進行折減變化，利用經(jīng)典斜坡滑塊案例測試剪切強度折減的位移閾值，模型如圖1所示，巖體參數(shù)見表1。在程序計算中，塊體結(jié)構(gòu)面抗剪強度參數(shù)每隔一定時步(10 000步，每一時步為0.001 s)按照塊體相對位移變化與閾值的關(guān)系只針對位移變化符合要求的結(jié)構(gòu)面進行折減，折減系數(shù)從1.000開始，依次增加0.001，直至發(fā)生急劇變形，則此時其取值為該邊坡滑落面安全系數(shù)。該程序計算所模擬的折減過程在一定程度上與實際情況下巖體材料的損傷退化相符。

圖1 滑塊沿斜坡滑動模型Fig.1 Model of sliding block along slope

該簡易模型運動問題安全系數(shù)解析解為

(4)

該算例結(jié)構(gòu)面選取各內(nèi)摩擦角對應(yīng)解析解值見表2。圖2為結(jié)構(gòu)面內(nèi)摩擦角φ取不同值下，模型安全系數(shù)及安全系數(shù)與解析解的誤差百分比隨位移閾值的變化曲線，由圖2可得，不論內(nèi)摩擦角φ取何值，當位移閾值接近1 mm時，安全系數(shù)取值接近最大值，且與解析解的誤差百分比有一個明顯的凹槽，說明這時的誤差與解析解相對最小且在0.5%以內(nèi)，所以位移閾值暫時取為1 mm，與解析解的結(jié)果保持一致。從而說明了DSR-DDA法是滿足計算要求的，并且能保證一定的計算準確性。

表1 邊坡模型巖體參數(shù)Table 1 Mechanical parameters of the slope model

表2 安全系數(shù)解析解Table 2 Theoretical safety factor

圖2 安全系數(shù)及其誤差隨位移閾值變化曲線Fig.2 Curves of safety factor and its error with displacement threshold

塊體間接觸力的精確求解是解決總平衡方程的核心步驟。在每一時步內(nèi)，都要重新確定彈簧的施加與否及彈簧的位置，需要反復(fù)生成求解總剛度矩陣，剛性彈簧的施加與去除過程稱之為開-合迭代。接觸有3種狀態(tài):張開、滑動和鎖定。模式變化的判定準則見表3，其中，N為法向位移，N>0為張開;t為剪切位移矢量;f為摩擦力矢量；T為剪切位移;‖為兩個矢量方向相同。

基于以上分析，DSR-DDA法在程序中實現(xiàn)如圖3所示。

表3 接觸狀態(tài)Table 3 Contact status

圖3 DDA程序步驟流程Fig.3 Flow chart outlining the procedures of the DDA

圖4 邊坡全貌及地質(zhì)剖面Fig.4 Slope topography and geological section

2 撫順西露天礦巖質(zhì)邊坡模型及參數(shù)

2.1 撫順西露天礦巖質(zhì)邊坡概況

露天開采在煤炭資源開采中占據(jù)至關(guān)重要的地位[13-14]。撫順西露天礦地處撫順市中心，其邊坡穩(wěn)定性直接決定了礦坑與城市交界處的安全與否。由于受到渾河斷裂的牽引作用,向斜北翼地層遭受了強烈改造和破壞。向斜軸部呈圓弧型褶曲,軸面呈現(xiàn)曲面狀,向北傾斜,傾角30°左右軸跡方向為NE60°,長為20 km,兩翼間距為3 km。北翼產(chǎn)狀各異,傾角15°～60°,間或伴有小型褶曲。目前,采坑北幫邊坡滑坡、崩塌、地裂縫等地質(zhì)災(zāi)害仍持續(xù)不斷，不時產(chǎn)生大型滑坡。故選取北幫W區(qū)某標段巖質(zhì)高邊坡作為案例。圖4為該標段邊坡全貌及地質(zhì)剖面圖，為監(jiān)測邊坡實時位移，在鉆孔中布置深部位移測斜儀。

2.2 計算分析模型及參數(shù)

圖5 UAV技術(shù)Fig.5 UAV technology

節(jié)理、裂隙等不連續(xù)結(jié)構(gòu)面對巖體的變形起到控制作用，結(jié)構(gòu)面信息采集精度對數(shù)值模擬分析準確性至關(guān)重要。運用高精度、高效率的UAV(圖5)技術(shù)對確定性結(jié)構(gòu)面信息進行非接觸精準獲取，生成點云模型。假定結(jié)構(gòu)面平面方程為

D=AX+BY+C

(5)

其中，A，B，C，D為結(jié)構(gòu)面平面參數(shù)，提取結(jié)構(gòu)面所在平面的法向量n=(-A,-B,1)，并對結(jié)構(gòu)面進行點的拾取，可得

(6)

借助最小二乘法解(A，B，C)，結(jié)構(gòu)面傾向及傾角分別為

(7)

式中，α，β分別為結(jié)構(gòu)面傾向和結(jié)構(gòu)面傾角。

對結(jié)構(gòu)面產(chǎn)狀進行計算，并對提取出的結(jié)構(gòu)面產(chǎn)狀信息導出，部分計算結(jié)果見表4，結(jié)構(gòu)面產(chǎn)狀統(tǒng)計云圖如圖6所示。

研究表明結(jié)構(gòu)面產(chǎn)狀服從一定規(guī)律分布，對確定性結(jié)構(gòu)面分組處理，并依照雙正態(tài)密度分布生成隨機結(jié)構(gòu)面產(chǎn)狀，其標準變化形式為

(8)

式中，σα為確定性結(jié)構(gòu)面傾向的標準差;σβ為確定性結(jié)構(gòu)面傾角的標準差;ρ為結(jié)構(gòu)面傾向、結(jié)構(gòu)面傾角的相關(guān)系數(shù);μα為確定性結(jié)構(gòu)面傾向的均值;μβ為確定性結(jié)構(gòu)面傾角的均值。

表4 產(chǎn)狀計算結(jié)果Table 4 Yield calculation results

圖6 結(jié)構(gòu)面產(chǎn)狀統(tǒng)計云圖Fig.6 Statistical cloud chart of structure surface

將獲取的結(jié)構(gòu)面信息導入DDA程序中，結(jié)構(gòu)面網(wǎng)格切割巖體，形成DDA的塊體單元，即算例數(shù)值計算分析模型(圖7)。

圖7 撫順西露天礦北幫邊坡DDA數(shù)值計算分析模型Fig.7 DDA numerical calculation and analysis model of slope in West Fushun Open-Pit Mine

天然應(yīng)力場僅將重力作用納入范圍，不計算區(qū)域構(gòu)造應(yīng)力影響。根據(jù)地質(zhì)勘探和試驗成果，得表5所示各項參數(shù)。程序計算中，塊體結(jié)構(gòu)面抗剪強度參數(shù)每隔一定時步(10 000步，每一時步約為0.001 s)針對塊體相對位移變化達到或超過位移閾值的結(jié)構(gòu)面進行折減，折減系數(shù)從1.000開始，依次增加0.001，在一定程度上與實際情況下巖體材料的損傷退化相符。

表5 邊坡模型巖體參數(shù)Table 5 Mechanical parameters of the slope model

3 撫順西露天礦巖質(zhì)邊坡穩(wěn)定分析

3.1 現(xiàn)場位移監(jiān)測及數(shù)值計算分析

測斜孔布置參考圖4，測點間隔為1 m，近2個月內(nèi)各測斜孔不同深度測得的累計位移值如圖8所示。

圖8 測斜孔測點累積總位移值Fig.8 Cumulative displacements of monitoring point in surveying slant hole

分析圖8可得，測斜孔69002內(nèi)前100 m左右測點總位移值隨時間增長變化極大，且在深度為100 m附近存在位移值劇增，說明測斜孔附近巖體已發(fā)生相對滑移，并有繼續(xù)滑動的趨勢，且滑動面在距坡面100 m附近。測斜孔55026內(nèi)前20 m左右測點總位移值增加相對迅速，且存在繼續(xù)增長的趨勢，后續(xù)測點位移變化效果不明顯，且趨于穩(wěn)定，推測該測斜孔上部巖體存在個別不穩(wěn)定塊體。測斜孔74003內(nèi)測點位移隨時間增長較慢，無明顯規(guī)律且趨于穩(wěn)定，說明暫時該區(qū)域巖體穩(wěn)定性較好。

為進一步確定邊坡潛在滑落面位置，根據(jù)DSR-DDA法，動態(tài)折減抗剪強度參數(shù)，對邊坡漸進失穩(wěn)過程進行表征，并給出滑落面與相應(yīng)的安全系數(shù)。選取該邊坡最危險滑落面與次級滑面驗證了方法的可行性，圖9～10為邊坡首次失穩(wěn)與二次失穩(wěn)過程。

第1.360×106時步前，各塊體監(jiān)測點位移及邊坡整體基本保持平穩(wěn)狀態(tài)，當運算至第1.360×106時步時，即折減系數(shù)達到1.136時，邊坡發(fā)生急劇變形，相應(yīng)塊體位移激增，且后續(xù)變形亦呈逐步增加趨勢，并且從邊坡整體變形中可以明顯看出下部滑坡體沿相應(yīng)滑落面發(fā)生滑移，說明邊坡發(fā)生首次失穩(wěn)，最危險滑落面所對應(yīng)安全系數(shù)為1.136。

邊坡首次失穩(wěn)結(jié)束后，塊體位置及內(nèi)部應(yīng)力重新分布。運用首次失穩(wěn)過程結(jié)束后的邊坡模型繼續(xù)運行程序，旨在尋找次級滑落面。

圖9 撫順西露天礦北幫邊坡首次失穩(wěn)過程Fig.9 First failure process of the north slope in West Fushun Open-Pit Mine

圖10 撫順西露天礦北幫邊坡二次失穩(wěn)過程Fig.10 Second failure process of the north slope in West Fushun Open-Pit Mine

第1.890×106時步前，邊坡整體趨于穩(wěn)定，當運算至第1.890×106時步時，即折減系數(shù)達到1.189時，邊坡再次發(fā)生急劇變形，形成次級滑落面，其對應(yīng)安全系數(shù)為1.189，且較首次邊坡失穩(wěn)其滑落面范圍更大，而且是在首次失穩(wěn)發(fā)生后邊坡強度略微減小后，與常規(guī)分析得出的結(jié)果大相徑庭(常規(guī)分析中不能剝離首次失穩(wěn)滑坡體，導致其對二次滑坡體起到保護作用，從而計算得到的安全系數(shù)偏高)。究其主要原因是邊坡首次發(fā)生失穩(wěn)后塊體位置及應(yīng)力已重新分布，故應(yīng)使用首次失穩(wěn)結(jié)束后的邊坡模型對次級滑落面進行分析，而以往研究中并未充分考慮這一因素。

撫順西露天礦邊坡案例破壞模式為牽引式滑動破壞，中下部油母頁巖對邊坡穩(wěn)定性起關(guān)鍵作用，故不可繼續(xù)開采剩余油母頁巖，以免引起邊坡上部大面積滑坡及礦坑-城市邊界滑坡。

3.2 DSR-DDA法與GeoSMA-3D耦合

GeoSMA-3D(巖土工程結(jié)構(gòu)與模型分析系統(tǒng))是一款以Key Block Theory為基礎(chǔ)，在C++框架下編譯源程序，團隊研發(fā)的三維關(guān)鍵塊體可視化分析軟件。該軟件能實現(xiàn)工程巖體空間結(jié)構(gòu)的建模、結(jié)構(gòu)面的空間表征、模型表面跡線顯示等多項功能，并且其精度已經(jīng)過大量實例驗證[15-17]。

將結(jié)構(gòu)面產(chǎn)狀導入GeoSMA-3D軟件，完成關(guān)鍵塊體搜索過程，并實現(xiàn)其三維表征;DSR-DDA法中最先發(fā)生滑落的塊體即為關(guān)鍵塊體，關(guān)鍵塊體耦合如圖11所示，且耦合效果良好。

圖11 關(guān)鍵塊體耦合Fig.11 Key block coupling

4 結(jié) 論

(1)提出了受位移閾值控制的DSR-DDA法，給出了DSR-DDA法計算邊坡多滑面安全系數(shù)的基本過程，為邊坡多滑面穩(wěn)定性分析貢獻了一種新手段。

(2)利用經(jīng)典斜坡滑塊案例驗證了DSR-DDA法的可行性與計算精度，并給出了位移閾值為1mm，對不同破壞程度下巖體結(jié)構(gòu)面動態(tài)折減不同的強度，解決了數(shù)值分析中巖體結(jié)構(gòu)面損傷的程度不一性問題。

(3)基于DSR-DDA法分析了撫順西露天礦北幫邊坡多滑面穩(wěn)定性，得到了其最危險滑落面及次級滑落面位置，與常規(guī)分析結(jié)果大相徑庭，佐證了邊坡多滑面穩(wěn)定性分析中模型需動態(tài)變化的必要性。該邊坡破壞模式為牽引式滑動破壞，中下部油母頁巖對邊坡穩(wěn)定性起關(guān)鍵作用，故不可繼續(xù)開采剩余油母頁巖，以免引起邊坡上部大面積滑坡及礦坑-城市邊界滑坡。

致謝衷心感謝石根華博士對筆者的熱忱幫助。