王玉鎖， 李俊杰， 李正輝， 馮高飛， 吳 浩， 何俊男

（1．西南交通大學土木工程學院，四川成都610031；2．西南交通大學地球科學與環(huán)境工程學院，四川成都610031）

落石沖擊力評定的離散元顆粒流數(shù)值模擬

王玉鎖1， 李俊杰2， 李正輝2， 馮高飛1， 吳 浩1， 何俊男2

（1．西南交通大學土木工程學院，四川成都610031；2．西南交通大學地球科學與環(huán)境工程學院，四川成都610031）

為使落石沖擊力評定更合理、有效，以離散元模型（DEM）為理論基礎，利用顆粒流（PFC）單元法研究垂直下落條件下落石對結構的沖擊力，分析了落石高度、重力及回填土厚度對沖擊力的影響規(guī)律．將PFC數(shù)值模擬結果與其他幾種落石沖擊力計算方法所得結果進行比較，結果表明：落石沖擊力與落石高度及重力均呈明顯的線性關系，落石沖擊力隨落石高度線性變化幅度與落石重力有關，而落石沖擊力隨落石重力線性變化幅度與回填土厚度有關；回填土厚度為0．6m時的沖擊力比無回填土時減小50%～60%，回填土厚度大于2．0m時，緩沖幅度趨于穩(wěn)定，說明回填土厚度有一個合理范圍，超出時，結構總荷載會增大；離散元顆粒流數(shù)值模擬方法通過阻尼設置來模擬空氣在落石下落過程中產生的作用，使得在評定大體積落石從較高高度下落時產生的沖擊力更加合理．

落石沖擊力；回填土；離散元模型（DEM）；顆粒流（PFC）

落石災害在我國山區(qū)鐵路、公路修建及運營過程中經常出現(xiàn)，例如寶成線K307＋100－K308＋720段20年間發(fā)生46起較大落石災害，損壞軌道、列車和輸電網，引發(fā)停車事故．5．12地震時，邊坡崩塌和落石引起寶成線109隧道洞口結構嚴重破壞．成昆鐵路北段，1971至1992年間發(fā)生落石災害214處共218次．據(jù)不完全統(tǒng)計，全國鐵路每年落石災害數(shù)千次．因此，落石災害已成為山區(qū)交通運輸安全的重要影響因素，落石對鐵路、公路工程結構的沖擊效應日益引起各方關注并展開深入研究［1-3］．

在落石沖擊力評估或計算方面，西南交通大學楊其新、關寶樹［4］在較早時期利用鐵球垂直下落沖擊結構頂板的方法，得出了垂直下落條件下落石沖擊力的計算經驗公式．我國路基規(guī)范、隧道手冊等根據(jù)沖量定理提出了相應的落石沖擊力的計算方法［5-6］．國外包括日本道路協(xié)會（Japan Road Association）方法［7］、Vincent Labiouse博士等提出的計算方法［8］．

由于試驗條件及計算參數(shù)選取等因素的限制，以上這些落石沖擊力評定方法在同樣條件下所得到的結果差異較大，說明各種方法都有一定的適用性和局限性［9］．防落石沖擊措施尤其是被動防護結構如隧道棚洞、拱形明洞等的設計目前仍主要采用經驗方法，缺乏設計理論指導．因此有必要深入研究落石沖擊力的確定方法．

近年來，采用離散單元模型（discrete element modeling，DEM）模擬落石災害的研究成為一種趨勢［10-13］．本文采用三維離散元顆粒流（particle flow code in 3 dimensions，PFC3D）方法，研究了在垂直下落條件下，落石高度、重力及回填土厚度等因素對落石沖擊力的影響規(guī)律，對合理回填土厚度進行了探討，介紹了離散元顆粒流方法的基本原理及模型設置，并與其他落石沖擊力計算方法得到的結果進行了分析比較．

1 離散元顆粒流理論基礎及模型設置

1．1 PFC3D簡介及離散單元模型基本理論

基于離散單元方法，PFC3D軟件由“Itasca公司”發(fā)行．模型是由可獨立運行的球體顆粒單元及墻單元組成，球體單元之間以及球體與墻體單元之間的相互作用是以牛頓第二運動定律為理論基礎［12，14-15］．本次所用到的本構關系包括接觸剛度、粘結、滑動等可參見文獻［16］．

1．2 落石沖擊模型及計算工況

初速度為零、垂直下落條件下，研究不同高度、重力的落石，對具有不同厚度的回填土的結構（頂面）引起的沖擊力．其中，落石采用單個剛性球體單元（ball），根據(jù)落石物理性質設定其容重或天然密度，通過改變球體顆粒的大?。ò霃剑┇@得不同重力的落石．本次研究沒有考慮落石下落及沖擊過程中的碎裂．利用多個球體單元（ball）的集合模擬回填土，通過接觸粘結模型，模擬回填土與落石的相互作用，以及回填土受力的分布及傳遞．回填土底部為半平面無限剛性體，用墻單元（wall）模擬．

應用顆粒流PFC3D建模軟件模擬落石在垂直下落條件下，對不同回填土厚度的結構沖擊過程．落石下落的過程中，考慮重力、黏滯阻力和碰撞力，其中阻尼采用黏滯阻尼，法向和切向黏滯阻尼系數(shù)都取0．2［17］．

回填土范圍為10 m×10 m（長×寬），用球體單元模擬，通過改變球體單元數(shù)量來獲得不同回填土厚度，在四周用墻單元將回填土約束，底部wall用來模擬結構．計算模型見圖1．

圖1 落石沖擊計算模型Fig．1 PFC3D rockfall model of impact progress on base structure with backfill soil

落石高度H、重力W及回填土厚度h取值見表1．

表1 落石沖擊數(shù)值模擬工況Tab．1 Working conditions of rockfall impact numerical simulation

為全面評價各因素對落石沖擊力的影響，本次采用全面組合，故共有5×6×7＝210種計算組合情況．

1．3 材料單元參數(shù)設置

1．3．1 PFC3D模型材料細觀參數(shù)

對所選用的各PFC3D單元（ball，wall）需要設置對應的細觀參數(shù)．根據(jù)顆粒流理論及相關資料，對落石沖擊模型的細觀參數(shù)設置如表2所示．

1．3．2 回填土顆粒模型宏觀力學參數(shù)的校準

由于PFC模型中，顆粒單元間的細觀力學參數(shù)并不是所模擬材料的實際宏觀力學參數(shù)，故需要對所采用的單元顆粒的細觀力學參數(shù)進行校準，使模型材料的微觀力學與實際工程材料的力學參數(shù)一致．

表2 模型材料單元細觀參數(shù)Tab．2 Micro-properties of the model

采用PFC3D模擬回填土的三軸試驗［18］．試驗樣本采取與回填土模型相同的細觀特性參數(shù)及排列方式，單元半徑為0．200 3 m．試驗尺寸高4 m，直徑2 m，圍壓分別取25、50、100、150和200 kPa．三軸試驗試樣見圖2．

圖2 PFC3D三軸數(shù)值模擬試驗Fig．2 Numerical simulation of triaxial test

由三軸數(shù)值模型試驗可得到相應的應力-應變、體積應變-軸向應變曲線，根據(jù)彈性理論，推導出彈性模量與壓縮模量間的換算關系，用應力為100～200 kPa段割線計算壓縮模量E12，用壓縮模量Es來表示土的壓縮性［19-20］，由此得出回填土模型材料的壓縮模量Es及泊松比ν．同時，根據(jù)不同圍壓下的試驗結果，可繪出回填土模型材料的莫爾圓，進而得出回填土模型材料抗剪強度的指標黏聚力c及內摩擦角φ．

根據(jù)土的彈性模量Eo與壓縮模量之間的關系如式（1），可得到彈性模量Eo．

由此，得出回填土模型材料的宏觀力學參數(shù)如表3．

表3 三軸試驗數(shù)值模擬宏觀參數(shù)結果Tab．3 Micro-properties results from numerical simulation of triaxial test

將該結果與工程地質手冊（2006版）［21］中土的平均物理力學指標對比，基本接近有關回填土的宏觀物理力學指標，故所設計的回填土顆粒模型及細觀參數(shù)能夠反映實際回填土的物理力學性能．

2 數(shù)值模擬結果及分析

2．1 落石沖擊力評定

由PFC3D計算得到回填土底部結構所受到的接觸力記錄情況如圖3所示．

圖3 回填土底部結構接觸力Fig．3 Impact force of wall

根據(jù)所得到的回填土底部結構接觸力（圖3）隨時間變化的記錄曲線，通過讀取接觸力突變處的峰值，該峰值減去由回填土自重引起的平均接觸力（圖3中呈水平分布曲線部分為回填土自重引起的接觸力），即可得到最大落石沖擊力，也就是本文中所稱的落石沖擊力．

圖4 落石沖擊力與落石高度的關系Fig．4 Relationship between impact force and height of rockfall

2．2 落石沖擊力影響因素及規(guī)律

根據(jù)計算結果，落石沖擊力與落石高度、重力及回填土厚度的關系如圖4～圖6所示（只列出部分結果）．

圖5 落石沖擊力與落石重力的關系Fig．5 Relationship between impact force and gravity of rockfall

由圖4～圖6可知，在垂直下落條件下：

（1）落石沖擊力隨落石高度增大而增大，且線性關系較為明顯（圖4），當回填土厚度一定時，沖擊力增大的幅度或趨勢線的斜率與落石重力有關，重力越大，沖擊力變化幅度也越大；從圖4可看出，各系列的曲線的垂直分離程度明顯，說明落石重力是影響沖擊力的重要因素．

（2）落石沖擊力隨落石重力增大而增大，且線性關系較為明顯，而無回填土時，結構受到的落石沖擊力明顯大于有回填土情況（圖5），落石高度一定時，沖擊力增大的幅度或趨勢線的斜率與回填土厚度有關，厚度越小，沖擊力變化幅度越大；從圖4曲線的上下分離看，無回填土趨勢線明顯高于有回填土的各趨勢線，說明回填土緩沖作用明顯，即使厚度為0．6 m也可緩沖至少一倍的沖擊力；同時從圖5可以看出，回填土厚度越大，沖擊力隨落石重力的變化趨勢越不明顯．

（3）從圖6可明顯看出回填土對結構所受落石沖擊力的緩沖作用，同時從圖5可明顯看出，在本次計算范圍內（回填土厚度0．6～5 m），回填土可緩沖至少一半的落石沖擊力．但當回填土達到一定厚度時，落石沖擊力變化會趨于穩(wěn)定，說明回填土厚度有一個合理值或范圍．在本次計算條件下，即落石垂直下落、最大高度為90 m及落石最大重力為50 kN（體積大約為2．5 m3）條件下，回填土厚度2～3 m為合理厚度，如再增大，則會使結構所受的總荷載效應（沖擊力＋回填土自重）增大．

2．3 數(shù)值模擬結果與已有方法計算結果對比

國內外現(xiàn)有的5種沖擊力計算方法分別為鐵路工程設計技術手冊隧道版［6］（后面簡稱“隧道手冊方法”）、西南交通大學楊其新、關寶樹方法［4］（后面簡稱“交大方法”）、公路路基設計規(guī)范［5］（后簡稱“路基規(guī)范方法”）、日本道路協(xié)會（Japan Road Association）方法［7］（后簡稱“日本方法”）和Vincent Labiouse博士等人提出的計算公式［8］（考慮到該公式主要建立在瑞士學者的研究基礎上，故后面簡稱“瑞士方法”）．前述方法所用的公式見所引文獻．

圖6 落石沖擊力與回填土厚度的關系Fig．6 Relationship between impact force of rockfall and thickness of backfill soil

日本方法中的拉梅常數(shù)λ建議取1 000 kN／m2，實際工程中此值變化范圍較大［7］，故本次研究對于拉梅系數(shù)的取值分別用建議值和彈性理論換算值［22］．

而瑞士方法中，ME是通過標準荷載板試驗得到的基床反力系數(shù)，在原文中作者按砂礫墊層取ME＝3 200 kN／m2［8］，本次研究中，ME分別按原文的取值和彈性換算［23］取值ME＝96 023．8 kN／m2．

現(xiàn)選取數(shù)值模擬工況（表1）中落石重力為5、 10、20、30、40和50 kN，回填土厚度h為1 m，落石下落高度H為10、30、50、70和90 m共30種工況，按以上5種計算方法分別計算落石沖擊力，并與顆粒流數(shù)值模擬結果進行對比．5種方法所需基本計算參數(shù)見表4，計算中涉及到的其他參數(shù)均是通過所給基本參數(shù)計算得到．表4及以上各方法中涉及到的回填土及落石的物理力學參數(shù)均按數(shù)值模擬方法中的參數(shù)選取．

表4 計算參數(shù)Tab．4 Caculate parameters

另外，為便于對比分析，也比較了無回填土即h＝0時的數(shù)值模擬結果．

由此，當給定落石重力時，由不同計算方法及參數(shù)，共得到9個系列的落石沖擊力隨高度變化的趨勢線，分別為：

①PFC數(shù)值模擬1 m厚回填土；

②PFC數(shù)值模擬無回填土；

③隧道手冊方法；

④路基規(guī)范方法；

⑤交大方法；

⑥日本方法按建議值（λ＝1 000 kN／m2）；

⑦日本方法按計算值（λ＝21 994．7 kN／m2）；

⑧瑞士方法，ME＝3 200 kN／m2；

⑨瑞士方法，ME＝96 023．8 kN／m2．

落石沖擊力比較情況見圖7，為便于比較分析，圖中的趨勢線對應的系列名稱總體上從上到下按落石沖擊力由大到小排列，例如圖7（a）中，PFC數(shù)值模擬無回填土（h＝0）時，對應落石沖擊力最大，瑞士方法（ME＝96 023 kN／m2）次之，以此類推，而PFC數(shù)值模擬結果（回填土厚h＝1m）結果最?。?/p>

圖7 不同方法沖擊力計算結果Fig．7 Results of impact force using different methods

從圖7可以看出，9個系列的趨勢線具有如下特征：

（1）PFC數(shù)值模擬無回填土、日本方法及瑞士方法涉及的λ和ME分別按彈性理論公式求得．3者結果（對應的趨勢線系列名稱分別為“PFC h＝0”、“日本（λ取值為21 994 kN／m2）”及“瑞士（ME取值為96 023 kN／m2）”較為接近，所得的落石沖擊力最大；當落石重力越大且落石高度也增大時，PFC數(shù)值模擬無回填土即“PFC h＝0”落石沖擊力與其他兩種方法差距越來越顯著．

（2）當落石重力較小時（如5 kN）時，隧道手冊方法得到的落石沖擊力變化趨勢與路基方法、日本及瑞士方法（λ和ME分別為1 000和3 200 kN／m2）三者結果接近，當落石重力及高度增大時，結果越來越大，也越接近PFC無回填土的結果．

（3）路基方法、日本及瑞士方法（λ和ME分別為1 000和3 200 kN／m2）三者結果接近，而路基方法偏大．

（4）當落石重力較?。ㄈ?，10 kN）時，PFC數(shù)值模擬有回填土即“PFC h＝1 m”結果與交大方法接近，而當落石重力較大時，則與路基方法、日本及瑞士方法（λ和ME分別按1 000和3 200 kN／m2）三者結果接近，稍大于交大方法的結果．

通過以上各結果的相互驗證、比較，可有如下分析或推論：

（1）日本方法與瑞士方法中，拉梅常數(shù)λ及基床反力系數(shù)ME的取值對結果影響較大，原方法中建議取λ＝1 000 kN／m2和ME＝3 200 kN／m2，所得結果較為合理，與我國路基規(guī)范方法所得結果較為接近；而這兩個參數(shù)如果按彈性理論公式換算后，所得沖擊力結果偏大，甚至比由PFC數(shù)值模擬無回填土結果還大，說明回填土緩沖效果不能按彈性理論求解，這也同時從側面反映出PFC數(shù)值模擬方法的合理性．

（2）由于瑞士、日本方法中，并未涉及落石下落過程中空氣等的阻力影響，而PFC數(shù)值模擬中則通過阻尼系數(shù)的設置（取推薦值0．2）來反映空氣阻力．當落石高度不高或重力（體積）不大時，影響不明顯；但當落石高度和體積增大，空氣的阻力影響程度會加大，這可從圖7中各圖的最上面3條趨勢線（“PFC h＝0”、“日本（λ＝21 994 kN／m2）”及“瑞士（ME＝96 023 kN／m2）”3個系列）隨落石高度增大而上下分離越來越明顯看出．

（3）PFC數(shù)值模擬1 m厚回填土（“PFC h＝1 m”）結果，當落石重力較?。?0 kN以下）或落石高度較?。?0 m以下）時，與交大方法結果接近，而當落石重力較大或高度較大時，則與路基方法、日本及瑞士方法（λ和ME分別按1 000和3 200 kN／m2）三者結果接近，比交大方法偏大而小于隧道手冊方法．

通過以上綜合比較分析，可認為顆粒流落石沖擊數(shù)值模擬方法能更加全面的考慮各種落石沖擊影響因素，所得到的結果也較為合理．

3 結 論

基于離散元模型理論，通過顆粒流PFC3D模擬分析，對垂直下落條件下落石引起的地基或結構的沖擊力進行了研究，分析了落石高度、重力及回填土厚度對沖擊力的影響規(guī)律，得出如下結論：

（1）落石沖擊力隨落石重力及高度增大而增大，且線性關系明顯，在其他條件一定時，可認為沖擊力分別與落石重力及高度成正比關系；當落石高度一定時，沖擊力隨落石重力的變化幅度與回填土厚度有關，厚度越小，變化幅度越大；當回填土厚度一定時，沖擊力隨落石高度的變化幅度與落石重力有關，重力越大，變化幅度越大．

（2）回填土厚度對落石沖擊力緩沖作用明顯，當回填土厚度為0．6m時，即可緩沖一半以上的落石沖擊力；當回填土厚度為2～3 m時，落石沖擊力變化已趨于穩(wěn)定，當進一步增大厚度，緩沖作用已不明顯，反而會使結構所受的總荷載（沖擊力＋回填土自重）增大．因此，從防落石角度出發(fā)，結構合理回填土厚度為2～3 m．

（3）由于試驗條件及參數(shù)選取問題，現(xiàn)有的國內外幾種落石沖擊力評定方法都有一定的局限性和適用性，而PFC數(shù)值模擬方法能方便的考慮各種影響因素．通過與國內外不同落石沖擊力方法的對比，用離散元顆粒流數(shù)值模擬方法得到的落石沖擊力較為合理．

［1］ 胡厚田．崩塌與落石［M］．北京：中國鐵道出版社，1998：71．

［2］ 王玉鎖，楊國柱．隧道洞口段危巖落石風險評估［J］．現(xiàn)代隧道技術，2010，47（6）：33-39．WANG Yusuo，YANG Guozhu．Rockfall risk assessment for a tunnel portal section［J］．Modern Tunnelling Technology，2010，47（6）：33-39．

［3］ 王玉鎖．高速鐵路隧道洞口段危巖落石運動軌跡及沖擊特性研究［J］．前言動態(tài)，2011（2）：16-21．

［4］ 楊其新，關寶樹．落石沖擊力計算方法的試驗研究［J］．鐵道學報，1996，18（1）：101-106．YANG Qixin，GUAN Baoshu．Test and research on calculating method of falling stone impulsive force［J］．Journal of the China Railway Society，1996，18（1）：101-106．

［5］ 交通部第二公路勘察設計研究院．JTJ013—95，公路路基設計規(guī)范［S］．北京：人民交通出版社，1995．

［6］ 鐵道第二勘測設計院．鐵路工程設計技術手冊．隧道［M］．2版．北京：中國鐵道出版社，1995：148-150．

［7］ KAWAHARA S，MURO T．Effects of dry density and thickness of sandy soil on impact response due to rockfall［J］．Journal of Terramechanics，2006，43（3）：329-340．

［8］ LABIOUSE V，DESCOEUDRES F，MONTANI S．Experimental study of rock sheds impacted by rock blocks［J］．Sruct．Eng．Int．，1996，3（1）：171-175．

［9］ 唐建輝．落石沖擊對隧道明洞結構的影響研究［D］．成都：西南交通大學，2013．

［10］ ASHAYER P．Application of rigid body impact mechanics and discrete element modeling to rockfall simulation［D］．Toronto：Department of Civil Engineering University of Toronto，2007．

［11］ THOENIK，GIACOMINIA，LAMBERT C，et al．A 3D discrete element modelling approach for rockfall analysis with drapery systems［J］．International Journal of Rock Mechanics＆Mining Sciences，2014，68：107-119．

［12］ CUNDALL P A，HART R D．Numerical modeling of discontinua［J］．Engineering Computations，1992，9（2）：101-113．

［13］ CUNDALL P A，STRACK O D L．A discrete numerical model for granular assemblies［J］．Geotechnique，1979，29（1）：47-65．

［14］ SHIUW J，DONZéF V．Numerical study of rockfalls on covered galleries by the discrete elements method［c］∥Annual Report 2004．Grenoble cedex 9 France：［s n］，2005：67-82．

［15］ DONZE F V，RICHEFEU V，MAGNIER S A．Advances in discrete element method applied in soil，rock，concrete mechanics［J］．Electronic Journal of Geotechnical Engineering，2009，8：1-44．

［16］ 鄭智能，張永興，董強，等．邊坡落石災害的顆粒流模擬方法［J］．中國地質災害與防治學報，2008，19（3）：46-49．ZENG Zhineng，ZHANG Yongxing，DONG Qiang，et al．Visual simulation of rock-fall of slope based on particle flow theory［J］．Journal of Geological Hazard and Control，2008，19（3）：46-49．

［17］ ItascaConsulting Group．Inc．PFC3D（particle flow code in three dimensions）version 4．0 manual［M］．Minneapolis：Itasca Consulting Group Inc．，2008：1-5

［18］ 李廣信．高等土力學［M］．北京：清華大學出版社，2004：114-130．

［19］ 王玉鎖，王明年，童建軍，等．砂類土體隧道圍巖壓縮模量的試驗研究［J］．巖土力學，2008，29（6）：1607-1617．WANG Yusuo，WANG Mingnian，TONG Jianjun，et al．Test research on compression modulus of sandy soil tunnel surrounding rock［J］．Rock and Soil Mechanics，2008，29（6）：1607-1617．

［20］ 王玉鎖．砂質土隧道圍巖力學參數(shù)及分級方法研究［D］．成都：西南交通大學，2008．

［21］ 工程地質手冊［M］．4版．中國建筑工業(yè)出版社，2006：192．

［22］ 王光欽．彈性力學［M］．北京：中國鐵道出版社，2008：131．

［23］ 黃斌，羅菊，何曉民等．基床系數(shù)的試驗方法及合理取值［J］．四川大學學報：工程科學版，2007，39（增刊）：61-65．HUANG Bin，LUO Ju，HE Xiaomin，et al．Research on test method and reasonably valuing of coefficient of subgrade reaction［J］．Journal of Sichuan University：Engineering Science Edition，2007，39（Sup．）：61-65．

（中、英文編輯：徐 萍）

Assessment of Rock fall Impact Force by Particle Flow Code Numerical Simulation Based on Discrete Element Model

WANG Yusuo1， LI Junjie2， LI Zhenghui2， FENG Gaofei1， WU Hao1， HE Junnan2
（1．School of Civil Engineering，Southwest Jiaotong University，Chengdu 610031，China；2．Faculty of Geosciences and Environmental Engineering，Southwest Jiaotong University，Chengdu 610031，China）

To make a reasonable and effective assessment of rockfall impact force，a study on the impact force of falling down rockfall to a structure was carried out by using the method of PFC（particle flow code）based on the theory of DEM（discrete element model）．The effects of rockfall height and gravity，and backfill soil thickness on the impact force were analyzed．The results of PFC numerical simulation were compared with those reached by other kinds of calculation methods．The results show that there is an obvious linear relationship between the impact force and the rockfall height or gravity．The linear variation range of the impact force with the changing of the rockfall height relates to the rockfall gravity，while the linear variation range of the impact force with the changing of the rockfall gravity relates to the backfill soil thickness．The impact force reduces by 50%-60%when the backfill soil thickness is 0．6 meter compared with no backfill soil．it tends to be stable when the backfill soil thickness is more than 2meters，showing that there is a reasonable range of the backfill soil thickness．The total load of structure will increase if the backfill soil thickness exceeds the range．The method of PFC numerical simulation can simulate the effect of air acting on rockfall by setting damping coefficient，which makes the assessment of bigger size rockfall impact force more reasonable．

impact force of rockfall；backfill soil；discrete element model；particle flow code

P642．21

A

0258-2724（2016）01-0022-08

10．3969／j．issn．0258-2724．2016．01．004

2015-11-16

四川省科技計劃資助項目（2013GZ0047）

王玉鎖（1974—），男，副教授，博士，研究方向為地下工程結構及工程材料，E-mail：wangysuo＠home．swjtu．edu．cn

王玉鎖，李俊杰，李正輝，等．落石沖擊力評定的離散元顆粒流數(shù)值模擬［J］．西南交通大學學報，2016，51（1）：022-029．