基于PFC3D和Rockfall的中武山危巖體數值模擬研究

戎澤鵬， 范宣梅， 董遠峰， 楊 帆， 戴嵐欣， 馮澤濤

(成都理工大學地質災害防治與地質環境保護國家重點實驗室， 成都 610059)

中國西南地區地質條件復雜，高陡岸坡和深切峽谷隨處可見，崩塌、滑坡、泥石流等典型地質災害特別發育，不僅嚴重影響這一地區的重大工程建設，同時給當地群眾生命財產帶來嚴重威脅[1]。中武山位于烏東德水電站上游地區，在烏東德水電站蓄水后，庫水位變動將引起兩岸地質災害易發[2]。

利用數值模擬方法對危巖體失穩運動過程及運動路徑的研究及預測相比普通經驗方法有更大的優勢，目前已被廣泛應用于地質災害及水利水電工程建設的相關研究中[3-8]。

1971年Cundall[3]提出適用于巖石力學的離散元法。常規的數值模擬方法可分為連續介質方法和非連續介質方法兩大類。連續介質方法包括有限元(FEM)、有限差分法(FDM)、邊界元(BEM)、無單元(EFM)等；非連續介質方法包括離散元(DEM)、顆粒元、流形元(MM)、非連續變形分析(DDA)等。PFC3D程序是美國ITASCA咨詢集團開發的顆粒流分析程序(particle flow code,PFC)，屬于離散元范疇，許多學者利用此方法對相關災害問題展開了研究[4-8]。Heim[4]通過模型試驗和典型滑坡分析，提出了顆粒流模型。施鳳根[5]采用顆粒流離散單元法對文家溝滑坡的運動堆積特征進行了PFC模擬分析；Scaringi等[6]采用PFC等4種軟件模型將新磨村滑坡事件進行了數值模擬再現還原；Robert等[7]對節理裂隙在巖體破壞中的控制作用進行了塊體單元法分析；Wang等[8]對邊坡穩定性問題進行了數值模擬。

Rockfall軟件是中仿科技公司開發的一款專門針對危巖體崩塌數值模擬的軟件[9]。該軟件以概率論為理論基礎，根據滾石隨機下落路徑統計分析的結果，最終給出可靠性較高的結論[10]。代婧瑜[11]利用概率分析中的正交試驗設計思路，結合Rockfall軟件通過大量的數值模擬收集防護結構控制指標，從而揭示出規律；胡聰[12]采用Rockfall軟件進行數值模擬，結合室內試驗的手段分析了高陡邊坡危巖體失穩的特征以及單體滾石的運動特征。

目前中國在危巖體崩塌災害預測數值模擬研究方面，大部分學者均采用PFC2D或者Rockfall等單個數值模擬軟件進行二維崩塌運動過程模擬分析研究。本研究對比以往的危巖體數值模擬研究，將PFC3D三維地質模型數值模擬方法結合Rockfall二維運動過程計算，以三維數值模擬方法代替了傳統的物理試驗方法，減少了繁瑣的試驗過程，使危巖體數值模擬具有了更好的可視化特性，為后續災害防治提供更多的參考。

1 中武山陡立危巖體地質概況

中武山位于四川省涼山彝族自治州會理縣，現皎平渡大橋上游約270 m，金沙江左岸。中武山中部位于高程1 016～1 070 m處，為皎平渡山洞遺址保護區，其上部山體存在一處板狀板裂式危巖體，巖層陡立，結構面發育，已發生局部崩塌，危巖體估算體積約15 000 m3，嚴重威脅皎平渡山洞遺址的安全運營。

1.1 崩塌區地質條件

中武山出露地層主要有下元古界康定群、中元古界會理群變質巖系，震旦系下統砂礫巖、上統碳酸鹽巖，三疊系上統、侏羅系、白至系碎屑巖類，不同侵入時期的巖漿巖。由于地處川滇山地，屬強烈構造侵蝕作用為主的中山地貌，地貌結構多是以丘狀高原面或分割山頂面為“基面”，中武山巖體節理裂隙發育，易形成崩塌危巖體。

中武山危巖體位于山體頂部，新田溝右側，高程范圍為1 070～1 210 m，平面基本呈長270 m、寬100 m的舌型分布，崩塌區面積約2.7×104m2；坡度較陡為35°～45°，坡面無明顯地形起伏，僅少量順坡面的溝槽，寬為1～2 m；臨空條件良好，共存在3處潛在不穩定危巖體塊體，最低處高程約為1 336 m，最高處高程約為1 470 m，如圖1所示。

1.2 崩塌區物源特征

崩塌區物源主要有強風化卸荷巖體、脫離母巖的危石等，均位于西北側后山上方坡體，如圖2所示。后山山頂為一聳立突出的尖棱狀山峰，由強風化卸荷巖體構成，推測為周圍山體錯斷下切形成，東南面臨空，高程范圍為1 230～1 210 m，陡崖與斜坡交界處高程為1 210 m，水平距離為270～300 m。

后山山峰處強風化卸荷巖體為前震旦系中元古界落雪組(Pt2l)白云巖，淺表層強風化，局部塊體突出懸空瀕臨失穩狀態。臨空面即坡面產狀約120°∠75°，坡頂面較平緩；白云巖巖性產狀45°∠77°；除此之外還發育一組與崖面和地層斜交的節理裂隙及一組平行崖面的卸荷裂隙。

L1：構造裂隙。產狀180°∠20°，壓剪性，與崖面近垂直，緩傾發育，裂隙面平直，可視延伸長1 m，張開性，局部充填植被根系和腐殖土，淺表層貫通性較好，間距為1～2 m。

L2：卸荷裂隙。高陡的臨空面及局部掉塊形成凹巖腔，引起強卸荷作用，形成卸荷裂隙，或加寬加大原構造裂隙。卸荷裂隙產狀120°∠85°，基本平行于危巖體臨空面，在崖壁淺表層豎向發育，延伸長1～2 m，張開性，較貫通，裂隙面平直，大多無充填。

1.3 崩塌堆積體特征

崩塌碎石堆積情況如圖3所示，崩落物以塊石為主，次為碎石。塊石巖性主要為白云巖、輝長巖、花崗巖、礫巖等。崩塌堆積碎塊石分選一般，塊徑不均勻，主要以20～50 cm為主，碎石以3～8 cm為主，最大崩落塊石體積約2 m×1 m×3 m。

圖3 崩塌碎石堆積Fig.3 Collapse and gravel accumulation

崩塌區中下部為緩傾斜坡，高程范圍為1 039～1 070 m，平均坡度為0°～30°，平緩的坡度和少量坡表植被為崩塌落石提供了良好的消能條件，形成崩塌堆積區。

2 基于PFC3D程序的危巖體三維運動軌跡分析

基于顆粒離散元方法對危巖體崩塌過程進行數值模擬分析，方法流程如圖4所示。

圖4 危巖體快速評價方法流程Fig.4 Flow chart of rapid evaluation method for dangerous rock mass

2.1 中武山危巖體三維模型建立

本次模擬使用了557個球形顆粒代替危巖體，模擬中主要使用線性接觸模型[13]，為模擬塊體及表面的相互作用，使用聯合剛度計算顆粒、墻之間的接觸剛度[14]；危巖體崩塌模擬過程中的能力損失[15]，通過加有阻尼器的塊體接觸模型，以其滯后阻尼的特性來展現[16]。選取了危巖體不同的3個部分進行崩塌運動過程的數值模擬研究，預測崩塌滾石其滑動距離以及滑動路徑及其范圍[17]。

通過采用隨機球形顆粒生成替代真實的危巖體顆粒進行崩塌三維運動過程模擬，危巖體顆粒數量為557個，球形顆粒半徑范圍為0.51～2.08 m，崩塌體量約13 300 m3。根據前期勘察報告資料中的巖性剛度及黏結強度，對球形顆粒進行合理參數取值(表1)，建成危巖體三維模型如圖5所示。

表1 PFC中使用的最佳擬合模型參數

圖5 危巖體三維模型Fig.5 3-D model of dangerous rock mass

2.2 基于離散元模型危巖體三維運動軌跡分析

通過PFC軟件對危巖體崩塌運動過程進行模擬，研究危巖體根據實際地形滾動過程中的運動路徑，并為后續Rockfall軟件對危巖體滾動、彈跳、沖擊能量等進行計算提供實際運動路徑剖面。

PFC3D軟件可較好地模擬其運動路徑，并可還原危巖體崩塌運動的全過程。

3 基于Rockfall軟件崩塌過程分析

Rockfall程序基于危巖體基本參數進行大量模擬和概率統計[18]，模擬出危巖體崩塌過程的運動路徑，并得出彈跳高度、沖擊能量、巖塊速度等參數[19]，其軟件具有運算過程快、運算結果準確、數據結果全面等特點，在科研領域或者生產實踐領域均被廣泛應用[20]。

根據PFC三維模型中的危巖體運動軌跡預測，確定了3個危險塊體的6個滾石運動剖面(圖7)。區內巖體破壞模式以傾倒變形破壞為主，根據Rockfall軟件對不同剖面的落石運動軌跡、彈跳高度、水平速度、沖擊能量等數據進行計算，從而判斷不同剖面的落石運動情況，根據以上模擬結果可進行高效預測并采取有效治理措施。

圖7 危巖體剖面布置Fig.7 Section arrangement of dangerous rock mass

3.1 危巖體運動計算算法

如前所述，區內危巖體破壞模式主要以傾倒變形破壞為主，危巖體傾倒脫離母巖后，多為順地形或溝槽發生運動，主要以滑落、滾動、彈跳為主，運動分區如圖8所示。

圖8 危巖體運動分區Fig.8 Zoning of dangerous rock mass movement

由于Rockfall軟件中的運動模型計算公式均為理想模型，故建立以下假定。

(1)在自由下落以及跳躍的過程中，為忽略運動路徑中的空氣阻力，故將危巖體產生的落石視為質點。

(2)將危巖體產生的落石在斜坡滾動過程中視為球體，并且具有定量的質量以及慣性。

(3)危巖體產生的落石在運動過程中無視其質量及半徑，質量及半徑只在判斷其運動模式中運用。

(4)將崩塌運動過程中的動力性能視為定量。

3.1.1 自由墜落速度

危巖體崩塌自由墜落過程中，設落石自某一時刻從(x1，y1)處以初速度脫離母巖或坡面后，根據運動學原理，失穩后的落石運動軌跡為

xi+1=xi+Δtvix

(1)

(2)

墜地瞬間的速度v(i+1)x和v(i+1)y為

v(i+1)x=Vix，v(i+1)y=viy

(3)

式中：(xi，yi)為i時刻的危巖體位置坐標；(xi+1，yi+1)為墜落點坐標；Δt為墜落時間，s。g為重力加速度。

3.1.2 碰撞彈跳高度

危巖體崩塌過程中的碰撞屬于非彈性碰撞，因此須把落石的運動速度分解為斜面的法向和切向，即

vin=vicosα，vit=visinα

(4)

危巖體崩塌滾落過程中，由于坡面巖土體性質與能量損失有關，故引入回彈系數表征碰撞過程中的能量損失，則法向與切線的速度分別為

(5)

(6)

式中：vi為落石碰撞坡面的入射速度，m/s;v′i為落石碰撞坡面后的反彈速度，m/s；其中vin和vit分別為落石碰撞坡面后的入射速度沿斜坡與切線的分速度，m/s；v′in和v′in分別為落石碰撞坡面后的反彈速度沿斜坡法向與切線的分速度,m/s；en和et分別為法向回彈系數與切向回彈系數。

則落石反彈速度v′i及其與坡面夾角θ為

(7)

(8)

落石落地速度分解：

β=θ-α

(9)

為了計算方便，引入高度h，H為反彈點至墜落點豎直高度(m)，存在

(10)

式(10)中：s為落石反彈拋程，m；smax為落石反彈最大拋程，m。

因此，基于運動學原理，落石彈跳最高點距離起跳點的距離：

(11)

最大彈跳高度

(12)

3.1.3 沖擊力計算

危巖體沖擊力的計算有助于設計防護治理措施，以及對威脅對象的破壞程度可有一個較好的預估。建議使用基于實驗數據和Hertz彈性碰撞理論的日本道路公團法，沖擊力公式為

p=2.1082/3λ2/5H3/5

(13)

式(13)中：p為危巖沖擊墻體時的沖擊力，kN；λ為拉梅系數，建議取1 000 kPa；Η為危巖自由下落的高度，m。

3.2 運動計算結果

運用Rockfall軟件對6條剖面分別進行崩塌運動模擬，落石隨機模擬50次危巖體運動軌跡(圖9)，巖石容重取26.1 kN/m3，分別得出危巖體崩塌在不同剖面運動產生的騰躍高度、沖擊能量以及水平運動速度，由于剖面3—3為主要威脅剖面，且威脅路徑下多為民房，故取剖面3—3數據處理結果圖進行運動計算結果說明。其統計分析結果如圖10～圖12所示，所有剖面的計算結果如表2所示。

由以上崩塌運動模擬結果可知，此危巖體未在計算剖面運動路徑中停止運動，將直接對危巖體下方民房區等產生威脅。

圖9 滾石運動軌跡計算Fig.9 Calculating chart of falling stone trajectory

圖10 滾石彈跳高度統計分析Fig.10 Statistical analysis of bounce height

圖11 滾石沖擊能量統計分析Fig.11 Statistical analysis of impact energy of rockfall

圖12 滾石水平速度統計分析Fig.12 Statistical analysis of horizontal velocity of falling rocks

表2 所有剖面計算結果

Rockfall崩塌落石運動計算結果表明，最大騰躍高度為1.2 m，最大沖擊能量為960～1 524 kJ。可確定危巖體崩塌運動將對居民區造成危害，根據危巖體崩塌運動過程中的騰躍高度、沖擊能量和巖塊速度等計算結果可布置有效的工程治理防御措施。

4 結論

傳統治理工程處理巖體邊坡崩塌等問題一般選用Rockfall軟件對所選剖面直接進行模擬計算，但是存在剖面選擇是否合理的問題。因此運用顆粒離散元PFC3D三維模型對危巖體崩塌整體運動過程進行模擬，根據模擬結果選取典型剖面，結合Rockfall軟件進行運動軌跡計算，得出其運動過程中彈跳、沖擊能量等數據，形成危巖體崩塌過程數值模擬快速預測方法，為治理設計提供準確、有效的預測數據。研究得出危巖體快速預測方法如下：

(1)通過無人機航拍影像獲取危巖體崩塌地區準確的地表信息，再通過前期勘察資料獲取等高線等信息建立崩塌三維地質模型。

(2)利用真實反映地表信息的三維模型通過離散元PFC軟件進行崩塌模擬，判斷出不同危巖體在不同剖面的運動情況，獲取危巖體墜落路徑以及危害范圍，最終選取典型剖面進行后續計算。

(3)通過Rockfall軟件對已經選取好的剖面進行運動軌跡計算，計算落石的彈跳高度、沖擊能量、水平速度等參數，根據以上結果確定危巖體崩塌危害范圍及強度大小，可有針對性地對危巖體進行高效危險性分析并采取有效治理措施。

通過以上的步驟，形成一種通過數值模擬對危巖體危害性進行快速預測的方法，與傳統方法相比，彌補了其在三維空間運動分析的缺失，可根據PFC3D三維模型準確預測崩塌運動路徑以及危害范圍，同時輔以Rockfall二維剖面沖擊能量等計算結果更為準確的布置治理預防措施，從而將危巖體崩塌造成的災害影響降至最低，形成一套針對于危巖體的數值模擬快速預測方法。