基于離散元的高陡堆石邊坡失穩過程模擬

鄭成成,龍小剛,胡廣柱,袁 祥,馬春輝,李高超

(1.陜西鎮安抽水蓄能有限公司，陜西 西安 710061；2.西安理工大學水利水電學院，陜西 西安 710048；3.省部共建西北旱區生態水利國家重點實驗室，陜西 西安 710048)

西部大開發戰略推動了我國西部地區的水利、交通、市政等基礎設施建設的快速發展。由于西部地區多高山峽谷地貌，在工程建設中易形成高陡邊坡工程，如由棄渣料堆積形成的高陡堆石邊坡。這些高陡邊坡在自重、水壓力、作業荷載、風化和侵蝕等外部因素的影響下，其發生失穩的可能性將極大增加。邊坡失穩對工程運行、建筑物和人員安全都將會產生嚴重影響，高陡邊坡的穩定性已經成為影響和制約工程建設與運營的關鍵問題之一[1]。因此模擬高陡邊坡失穩過程，明確其影響范圍是工程領域中重要的研究課題。

邊坡穩定問題一直是工程建設中的重要研究內容，眾多學者采用各種穩定分析方法求解邊坡穩定相關問題，例如極限平衡法、有限元法、離散元法等。在傳統極限平衡法方面，陳祖煜等[2]將二維Spencer法擴展到三維條件下，提出了邊坡穩定分析的三維極限平衡法，并通過工程實例證明了其可行性。盧坤林等[3]將條柱間作用力等效為滑面正應力，提出了適用于空間形態滑面的邊坡三維極限平衡法。馮樹仁等[4]提出了能夠考慮多種滑面類型的邊坡三維極限平衡方法，并提供了計算程序和驗證算例。Jiang等[5]在極限平衡分析方法的框架下考慮了土體性質的二維空間變異性，提出了定量的邊坡破壞風險評估方法。此外，有限元法在邊坡穩定分析中也取得了良好的應用效果。史卜濤等[6]將坡頂豎直位移突變作為邊坡失穩判據，為邊坡穩定性分析提供了新思路。林姍等[7]提出了虛單元強度折減法，較真實地模擬滑裂面形態的多樣性，用于土石混合體邊坡的穩定性分析。蔣水華等[8]提出基于非侵入式隨機有限元法的邊坡可靠度分析方法，是解決復雜邊坡可靠度問題的有效手段。由于極限平衡法和有限元法分析方法具有連續介質特征，對非黏土、巖石混合體等離散介質的邊坡分析效果并不十分理想，Lu等[9]通過對比滑動穩定分析方法，認為離散元法更適合研究巖、土體破壞特性。Weng等[10]采用離散元法研究了邊坡角度、裂隙角度和材料劣化等影響下板巖邊坡的變形行為。毛佳等[11]根據滑坡問題的普遍性與重要性，探討了三維滑坡體的運動路徑、堆積范圍，提出了三維距離勢函數離散單元法。周健等[12]運用離散元方法分析土坡穩定性，并分析了其破壞機理。蔣明鏡等[13]提出了膠結尺寸離散元微觀接觸模型，研究了不同節理邊坡的破壞形式，并模擬了邊坡失穩演化過程。Xu[14]運用離散元方法模擬三軸試驗，研究了顆粒細觀參數對黏性土宏觀特性的影響。李世海等[15]通過對比滑坡穩定分析方法，認為將計算模型與現場監測結果相結合，是判斷邊坡穩定的最有效方法。Lu等[16]認為離散元方法能夠較好地模擬強降雨條件下徑流路徑、顆粒速度和滑動影響范圍，為滑動災害預警和決策支持提供有效信息。汪儒鴻等[17]結合邊坡災前位移變形程度、坡體解體現象等，采用顆粒離散元方法分析堆積體邊坡的突變失穩問題。趙蘭浩等[18]提出了適用于非均勻離散顆粒體系的接觸檢測算法，具有內存占用較小等優點，對于大規模、分布密集、大粒徑比的顆粒體系具有較高的計算效率。Wang等[19]建立了基于位移統計的離散元分析方法，并將其與抗剪強度折減法相結合，分析高陡順層巖質邊坡的穩定性。黃達等[20]為進一步探明軟硬互層反傾邊坡的傾倒變形機制，融合離心模型試驗與離散元模擬，研究了此類邊坡的破壞模式與影響因素。

通過上述分析可知，學者們采用極限平衡法、有限元法、離散元法等方法研究了水利、巖土工程中的巖質、土質邊坡穩定性問題，但對于特殊土體的邊坡穩定性及其失穩過程研究較少。為此，本文利用離散元方法對高陡堆石邊坡進行數值模擬，重點研究極限工況下的邊坡失穩過程，詳細分析其堆石運動、擋墻受力以及邊坡整體變化情況，旨在為高陡堆石邊坡穩定分析、失穩過程和影響范圍研究提供思路和方法。

1 離散元原理及模擬方法

離散元數值模擬方法不僅可用于砂石等非連續介質的細觀力學行為研究，還可用于求解實際工程中非連續變形問題[21]。離散元的基本思想是將非連續介質分成有限個單元的集合，并使每個剛性單元體滿足接觸本構關系。隨后，根據牛頓第二定律，建立離散元方法的運動方程并進行積分求解，得到單元體的速度、位移等物理量，經過多次循環計算后求得非連續介質的運動狀態。離散元的主要優勢體現在對非連續介質問題的處理上，允許單元間存在相對運動，可以將單個微觀介質的運動狀態進行宏觀放大，具有原理簡單、計算速度快等優勢。因此，離散元方法適用于模擬塊體屬性明顯的物體的靜、動力特性，如高堆石邊坡失穩過程模擬，其在重現堆石邊坡或塊石的位移、速度等方面具有很好的效果。離散元原理如圖1所示。

圖1 工程棄渣區域與辦公區域平面布置

在建立堆石邊坡離散元模型過程中，所采用的主要模擬技術與方法如下：①不同粒徑塊石生成：根據粒徑級配曲線，采用擠壓排斥法生成堆石料離散元三軸試驗模型與堆石邊坡離散元模型。②不平衡力計算：當邊坡無明顯位移，顆粒平均不平衡力小于0.1N，且最大不平衡力與平均不平衡力之比小于10時，可認為離散元模型已處于最終的穩定狀態。③能量計算：通過在每個時間步對模型中的墻體、球體進行能量分析，并將勢能、動能和摩擦耗能等進行分類合計，由此可得整個離散元模型的能量變化。④擋墻上作用力計算：通過監測球體作用在擋墻上的合力變化，實現對擋墻作用力的監測。

2 工程背景與離散元模型建立

某日調節抽水蓄能電站工程規模為Ⅰ等大(1)型工程，多年平均發電量23.41億kW·h。電站上水庫正常蓄水位1 392.00 m，死水位1 367.00 m，有效庫容856.00萬m3；下水庫正常蓄水位945.00 m，死水位910.00 m，有效庫容956.10萬m3。為修建上下庫連接路，在道路沿線的自然溝谷或地形低洼處修建了4個施工平臺，其中1號施工平臺與工程棄渣場的平面布置如圖1所示。建成后的1號施工平臺為大型堆石邊坡，最高處高程為1 116.40 m，最低高程為916.60 m，最大高差199.80 m。

在堆石邊坡初始建設階段，由于未能及時有效地平整棄渣料，導致降雨后棄渣體排水不暢，出現了局部泥石流現象。該次滑動最遠滑至邊坡下游與主溝道匯合處，存在影響工程棄渣場排洪洞泄洪安全的隱患，現場情況及匯合處位置關系如圖2所示。小規模泥石流發生后，在堆石邊坡的中部、底部分別修建了石籠擋墻和混凝土擋墻，以增強堆石邊坡的整體穩定性。因此，考慮到工程后期運行管理營地位于棄渣場主溝道的下游，堆石邊坡的可靠運行將直接影響運行人員生命安全及抽水蓄能工程的安全運行。

圖2 早期泥石流現場情況

采用傳統極限平衡法計算堆石邊坡穩定性，其安全系數略大于1，安全儲備較小。為進一步分析高陡堆石邊坡的穩定性，綜合考慮堆石邊坡所在區域的地形特征、邊坡體型特征和堆渣特點等，確定了最不利的邊坡二維剖面并建立了離散元模型，其現狀局部如圖3(a)所示，二維剖面如圖3(b)所示。相較于三維模型，二維模型忽略了實際地形對坡體、滾石等緩沖作用，將緩沖作用作為堆石邊坡的安全富裕，因此二維堆石邊坡穩定分析針對的是邊坡安全穩定的最不利情況。本文使用商業軟件進行離散元建模與計算，其水平與豎直方向的尺寸分別為945 m×312 m，共包含10 528個塊石顆粒，其底部基礎、石籠擋墻、混凝土擋墻均由wall構建[22]。關于堆石邊坡離散元中的塊石模擬，本文對堆石邊坡材料級配曲線進行適當地截取與放大，使模型在計算精度與計算效率間得到了較好的平衡，并通過模擬施工工況發生的小規?；瑒?，驗證了離散元模型的可靠性。

圖3 1號施工平臺(單位:m)

依據堆石料室內大三軸試驗結果、堆積邊坡級配曲線，建立了堆石邊坡離散元三軸試驗模型，對堆石料的細觀線性剛度模型參數進行標定?，F場取樣獲得的堆石料級配曲線如圖4(a)所示，據此建立的堆石料離散元三軸試驗模型如圖4(b)所示。采用基于本構模型的堆石料細觀參數標定方法進行計算[20]，確定顆粒的法向接觸剛度kn=3.5 MN/m，切向接觸剛度ks=2.6 MN/m以及摩擦因數μ=0.09，其標定結果如圖5所示。此外，在離散元計算中設定基巖面與塊石的摩擦因數均設定為0.30，考慮到系統中存在的其他阻尼，將局部阻尼設定為較小值0.10。通過查閱資料以及參考相關的工程經驗，確定石籠擋墻和混凝土的摩擦因數分別為0.35和0.20。

圖4 堆石邊坡離散元三軸試驗模型

圖5 堆石邊坡的室內三軸試驗與離散元模擬結果

3 堆石邊坡失穩過程模擬

在建立上述離散元模型的基礎上，模擬運行工況下堆石邊坡的失穩過程，整個失穩過程共持續418.38 s，其平均不平衡力比的變化過程如圖6所示，其體型變化如圖7所示。由圖可知：0～32.32 s為滑動啟動與快速滑動階段，通過對比圖7(a)與圖7(b)可知，在該階段堆石邊坡頂部、中部處于不穩定狀態的塊石迅速發生滑動。堆石邊坡頂部1 070 m高程以上的塊石滑向下游，直接導致石籠擋墻以上邊坡的塊石層厚度明顯增加。堆石邊坡中部的水平向300～400 m區間內塊石滑向下游，混凝土擋墻雖對塊石運動起到了一定的阻擋作用，但是仍有大量顆粒越過混凝土擋墻滑向下游；32.32～139.52 s堆石邊坡進入局部滑動階段，對比圖7(b)與圖7(c)可知，該階段堆石邊坡上部塊石越過石籠擋墻，堆積到堆石邊坡中下部的后緣部位；139.52～189.21 s堆石邊坡滑動緩慢，僅有個別塊石越過石籠擋墻滑向下游；隨著堆石邊坡中下部前端塊石的下滑，在后緣塊石的推動作用下，堆石邊坡中下部再次發生滑動。通過對比圖7(d)與圖7(e)可知，堆石邊坡中下部和后部變化較為明顯，且部分塊石已滾落至河床部位；在235.99～418.37 s時段內僅有個別塊石運動，滾落至河床的塊石數量有所增加，直至邊坡整體達到穩定狀態。

圖6 堆石邊坡失穩的平均不平衡力比過程線監測

圖7 運行工況下堆石邊坡各階段體型

在堆石邊坡失穩過程中，塊石平均水平向速度、平均豎直向速度和平均速度變化情況如圖8所示，水平向速度以順坡向為正，豎直向速度以向下為正。由圖8可知，塊石的速度變化情況與堆石邊坡失穩過程中各階段的劃分基本一致：①在邊坡失穩的初始階段，塊石平均速度急劇增大，隨后逐步減小，其平均速度最大值出現于3.73 s，為4.34 m/s。過程中，塊石的平均水平向速度明顯大于平均豎直向速度，這表明該階段塊石主要以水平向順坡運動為主；②在后續各階段中，塊石速度的變化規律與各階段的運動態勢基本相符，多以塊石的水平向運動為主，豎直向運動為輔。

圖8 堆石邊坡堆石速度變化

堆石邊坡初始狀態與失穩后狀態的塊石分布情況如圖9所示，水平方向上塊石分布范圍由初始狀態的0～480 m失穩后擴展至40～840 m。在堆石邊坡上部，峰值的塊石數量由383個增加至403個，表明石籠擋墻起到了明顯的阻滑作用，堆石邊坡上部的塊石堆積現象更加明顯。在堆石邊坡中下部，其水平方向的顆粒分布曲線明顯后移，表明中下部塊石滑向下游，塊石在水平向上的分布更加分散。對于塊石的豎直向分布，其變化規律、原因和水平向分布基本相同。此外，塊石所在位置的最大高程值下降，表明最高處的塊石出現下滑，同時堆石邊坡中下部塊石分布曲線出現下移，表明該處塊石出現明顯下滑。相比于石籠擋墻，受混凝土擋墻阻擋的堆石邊坡中下部塊石分布曲線改變更為明顯，表明堆石邊坡底部的混凝土擋墻更需要加強。

圖9 堆石邊坡運行工況模型

堆石邊坡失穩發生后，塊石的水平向、豎直向位移統計如圖10所示，由圖可知：對于水平向位移，運動距離為40 m的塊石數量最多，運動距離超過200 m后塊石數量明顯減??；對于豎直向位移，運動距離為10 m左右的塊石數量最多，塊石的最大下降高度為148.5 m。

圖10 運行工況邊坡失穩后堆石邊坡塊石位移統計

通過監測越過石籠擋墻、混凝土擋墻的塊石數量，能夠明確邊坡失穩過程以及擋墻所發揮的阻滑作用，其變化情況如圖11所示，由圖可知：對于石籠擋墻，其以上部位的塊石在滑動初期不斷越過石籠擋墻滑向下游，直至123 s后無塊石越過石籠擋墻。上述分析表明：石籠擋墻的修建有效地阻擋了堆石邊坡上部塊石向下游滾落，對維持其上部塊石的穩定起到了明顯的阻滑作用。對于混凝土擋墻，滑動啟動后，邊坡中下部的塊石不斷越過擋墻，直至264 s后無明顯增加。對于到達河床底部的塊石數量，其在213 s后幾乎不再發生變化，到達河床底部的塊石多數來源于邊坡中下部的堆石，因此加強邊坡底部工程措施將有助于減少塊石滾落至河床的可能性。因此，石籠擋墻、混凝土擋墻的修建起到了增加堆石邊坡穩定性及減少塊石滑落的重要作用。

圖11 堆石邊坡越過特定位置塊石數量變化監測

堆石邊坡滑動結束后，越過石籠擋墻、越過混凝土擋墻和到達河床的塊石粒徑情況如圖12所示。由圖可知：越過石籠擋墻、越過混凝土擋墻的塊石粒徑均為均勻，與堆石邊坡整體的粒徑基本一致。對于距離更遠的河床位置，粒徑大小表現出明顯的遞增趨勢，表明粒徑大的塊石更容易滾落至更遠的位置。

圖12 運行工況下越過特征位置的堆石邊坡粒徑統計

通過監測擋墻上所受壓力的變化，能夠掌握堆石邊坡塊石滑動下擋墻承受的沖擊力情況，以及擋墻所受到的最大壓力，從而為擋墻的設計提供依據。在堆石邊坡失穩過程中，石籠擋墻與混凝土擋墻的單寬受力變化如圖13所示。由圖可知，滑動啟動后，在大量塊石的沖擊下擋墻受力迅速增加，隨著滑動的發展擋墻受力處于波動狀態，當滑動逐漸停止時擋墻受力逐漸趨于穩定。對于石籠擋墻，其初始受力為1.98 MN/m，在塊石的沖擊下擋墻受力在57.83 s時達到峰值9.51 MN/m，因此擋墻在動力情況下的受力為靜力情況的4.8倍。對于混凝土擋墻，其初始受力為1.10 MN/m，在塊石的沖擊下擋墻受力在10.93 s時達到峰值4.15 MN/m，因此擋墻在動力情況下的受力為靜力情況的3.8倍。相比于石籠擋墻，混凝土擋墻在邊坡失穩過程中的受力變化不斷發生變化，表明堆石邊坡的中下部處于不斷調整中，其對混凝土擋墻受力亦有著持續的影響。

圖13 堆石邊坡擋墻上游面受力變化監測

在堆石邊坡失穩的過程中，堆石邊坡的能量變化如圖14所示?；瑒芋w的動能變化反映了在滑動過程中滑動體所攜帶的能量及其對防護工程造成破壞的能力。在t=3.51 s時滑動體動能出現峰值，與上述分析中塊石平均速度達到最大值的時刻相同，此時滑動體動能可達3.05×108J。峰值時刻坡體將攜帶巨大能量沖擊石籠擋墻與混凝土擋墻，對其安全穩定造成一定的影響。重力做功在滑動體啟動初期增長較快，后期增長速度變慢，表明堆石邊坡逐步趨于穩定。接觸累積耗能也隨時間逐步增長并趨于穩定，表明滑動體在運動過程中內部不斷發生碰撞、摩擦，消耗了大量能量。阻尼累積耗能主要為局部阻尼耗能。

圖14 堆石邊坡失穩過程中的能量變化

4 結 語

針對高陡堆石邊坡，采用離散元法對邊坡失穩過程進行定量分析。主要結論如下：①高堆石邊坡的失穩過程總體可分為啟動與快速滑動階段、局部滑動階段、平穩階段、再次局部滑動階段、最終平穩階段；②在整個邊坡失穩過程中，塊石平均速度初始階段較大后緩慢減小，主要以水平向順坡為主，豎直向運動為輔；③通過比較初始階段和最終狀態塊石分布、塊石越過特定位置的數量、粒徑統計和擋墻上游面的受力情況表明：采用石籠擋墻等工程措施能夠明顯起到了阻滑的作用，但混凝土擋墻迫切需要加強；④在堆石體邊坡失穩過程中，滑坡體將攜帶的大量動能，對防護措施造成的一定沖擊。通過離散元數值模擬，揭示了高陡堆石體邊坡整個失穩過程，并驗證了施工措施的阻滑效果，為特殊材料工程的失穩過程模擬、工程措施設計提供了模擬思路和方法。