基于動態強度折減DDA法的邊坡多滑面穩定性分析

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

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

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

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

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

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

1 DSR-DDA法

1.1 DDA基本原理

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

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

(1)

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

(2)

基于最小勢能原理構建的系統整體平衡方程為

(3)

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

1.2 DSR-DDA法位移閾值

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

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

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

(4)

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

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

表2 安全系數解析解Table 2 Theoretical safety factor

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

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

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

表3 接觸狀態Table 3 Contact status

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

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

2 撫順西露天礦巖質邊坡模型及參數

2.1 撫順西露天礦巖質邊坡概況

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

2.2 計算分析模型及參數

圖5 UAV技術Fig.5 UAV technology

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

D=AX+BY+C

(5)

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

(6)

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

(7)

式中，α，β分別為結構面傾向和結構面傾角。

對結構面產狀進行計算，并對提取出的結構面產狀信息導出，部分計算結果見表4，結構面產狀統計云圖如圖6所示。

研究表明結構面產狀服從一定規律分布，對確定性結構面分組處理，并依照雙正態密度分布生成隨機結構面產狀，其標準變化形式為

(8)

式中，σα為確定性結構面傾向的標準差;σβ為確定性結構面傾角的標準差;ρ為結構面傾向、結構面傾角的相關系數;μα為確定性結構面傾向的均值;μβ為確定性結構面傾角的均值。

表4 產狀計算結果Table 4 Yield calculation results

圖6 結構面產狀統計云圖Fig.6 Statistical cloud chart of structure surface

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

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

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

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

3 撫順西露天礦巖質邊坡穩定分析

3.1 現場位移監測及數值計算分析

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

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

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

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

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

邊坡首次失穩結束后，塊體位置及內部應力重新分布。運用首次失穩過程結束后的邊坡模型繼續運行程序，旨在尋找次級滑落面。

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

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

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

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

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

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

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

圖11 關鍵塊體耦合Fig.11 Key block coupling

4 結 論

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

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

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

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