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

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

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

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

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

1 離散元原理及模擬方法

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

圖1 工程棄渣區(qū)域與辦公區(qū)域平面布置

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

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

某日調(diào)節(jié)抽水蓄能電站工程規(guī)模為Ⅰ等大(1)型工程，多年平均發(fā)電量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。

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

圖2 早期泥石流現(xiàn)場情況

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

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

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

圖4 堆石邊坡離散元三軸試驗?zāi)Ｐ?/p>

圖5 堆石邊坡的室內(nèi)三軸試驗與離散元模擬結(jié)果

3 堆石邊坡失穩(wěn)過程模擬

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

圖6 堆石邊坡失穩(wěn)的平均不平衡力比過程線監(jiān)測

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

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

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

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

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

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

圖10 運行工況邊坡失穩(wěn)后堆石邊坡塊石位移統(tǒng)計

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

圖11 堆石邊坡越過特定位置塊石數(shù)量變化監(jiān)測

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

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

通過監(jiān)測擋墻上所受壓力的變化，能夠掌握堆石邊坡塊石滑動下?lián)鯄Τ惺艿臎_擊力情況，以及擋墻所受到的最大壓力，從而為擋墻的設(shè)計提供依據(jù)。在堆石邊坡失穩(wěn)過程中，石籠擋墻與混凝土擋墻的單寬受力變化如圖13所示。由圖可知，滑動啟動后，在大量塊石的沖擊下?lián)鯄κ芰ρ杆僭黾樱S著滑動的發(fā)展擋墻受力處于波動狀態(tài)，當滑動逐漸停止時擋墻受力逐漸趨于穩(wěn)定。對于石籠擋墻，其初始受力為1.98 MN/m，在塊石的沖擊下?lián)鯄κ芰υ?7.83 s時達到峰值9.51 MN/m，因此擋墻在動力情況下的受力為靜力情況的4.8倍。對于混凝土擋墻，其初始受力為1.10 MN/m，在塊石的沖擊下?lián)鯄κ芰υ?0.93 s時達到峰值4.15 MN/m，因此擋墻在動力情況下的受力為靜力情況的3.8倍。相比于石籠擋墻，混凝土擋墻在邊坡失穩(wěn)過程中的受力變化不斷發(fā)生變化，表明堆石邊坡的中下部處于不斷調(diào)整中，其對混凝土擋墻受力亦有著持續(xù)的影響。

圖13 堆石邊坡?lián)鯄ι嫌蚊媸芰ψ兓O(jiān)測

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

圖14 堆石邊坡失穩(wěn)過程中的能量變化

4 結(jié) 語

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