潘 青 張清照 李藝靈

(①同濟大學土木工程學院地下建筑與工程系，上海 200092，中國)(②同濟大學巖土及地下工程教育部重點實驗室，上海 200092，中國)

0 引 言

隨著國家大力推進公路建設(shè)發(fā)展進程，我國公路建設(shè)的重點也由東南沿海轉(zhuǎn)向了西南欠發(fā)展地區(qū)，該地區(qū)內(nèi)公路往往需要布設(shè)于沿河傍山地段，由于復雜的地質(zhì)環(huán)境條件及人類工程活動，泥石流災害十分發(fā)育，據(jù)不完全統(tǒng)計，我國西部幾大交通干線公路受泥石流災害影響多達3600余處(馬東濤等，2003)。特別是每年的雨季，由于泥石流造成的公路阻斷、堵車等問題層出不窮，以2019年8月20日四川阿壩州特大泥石流為例，泥石流的發(fā)生既對公路工程本身造成了損害，同時造成的公路中斷也對受災地區(qū)人民生命財產(chǎn)的及時搶救帶來了影響，一定程度上加重了災區(qū)災情。目前對于泥石流的工程防治措施可分為跨越工程、穿過工程、排導工程、攔擋工程等，為了合理進行此類工程的設(shè)計與施工，需要對泥石流的運動規(guī)律及其沖擊性能開展系統(tǒng)研究。

近年來，國內(nèi)外專家學者分別通過試驗研究、數(shù)值模擬以及理論分析等手段對泥石流的動力性能展開了大量研究，成果頗豐。何曉英等(2016)分別進行了清水流、水石流與泥石流工況下的泥石流沖擊特性模型試驗，模擬了溝谷的真實形態(tài)，考慮了漿體黏度與級配顆粒組合，采用小波理論對沖擊荷載時程曲線及沖擊荷載均值進行了分析；Scheidl et al.(2013)進行了小規(guī)模泥石流沖擊特性模型試驗，基于16組試驗數(shù)據(jù)得出了泥石流最大沖擊力的計算方法；劉興榮等(2018)通過開展室內(nèi)水槽試驗，研究了不同粒徑級配條件下工程棄渣泥石流啟動機理。相較于試驗研究方法，數(shù)值模擬由于其分析過程穩(wěn)定可靠、成本低廉，已成為泥石流災害研究的一種重要方式。李兆華等(2018)基于FEMLIP對泥石流固體、流體及固流轉(zhuǎn)化3個階段進行了模擬，并開展了泥石流對擋土墻的沖擊效應的啟發(fā)式研究；韓征等(2019)建立了整合HBP本構(gòu)關(guān)系的稀性泥石流動力過程三維數(shù)值模型；馬秋娟等(2015)，周健等(2015)均利用PFC軟件對泥石流啟動過程進行數(shù)值模擬研究；柳春等(2019)基于SPH-FEM耦合數(shù)值方法，考慮了大塊石-泥石流漿體-攔擋壩的相互作用，分析了泥石流沖擊被攔擋的過程、沖擊力時程等力學行為；湯碧輝等(2018)提出了黏性泥石流溝內(nèi)大顆粒的滾動啟動模型，并用FLUENT對多種工況下的顆粒啟動進行了數(shù)值模擬，驗證了模型的適用性；王秀麗等(2015)采用流固耦合分析方法，將CFX與ANSYS相結(jié)合，將常規(guī)重力壩與新型帶支撐攔擋壩進行了對比分析，得出了結(jié)構(gòu)極限抗沖擊承載特性；胡卸文等(2019)采用有限元與無限元耦合方法，系統(tǒng)分析了落石沖擊荷載下樁板攔石墻的抗沖擊特性。

從上述研究成果可以發(fā)現(xiàn)，泥石流的運動機理十分復雜，而其中固體顆粒間表現(xiàn)出的摩擦特性是泥石流-防治結(jié)構(gòu)相互作用過程中能量耗散的關(guān)鍵，因此，從固體顆粒出發(fā)開展泥石流基本運動機理及其沖擊性能的研究可為后續(xù)研究多相介質(zhì)運動提供堅實的基礎(chǔ)。本文基于離散單元法，通過EDEM 2018軟件進行碎屑流沖擊流槽試驗模擬，綜合考慮流槽坡度、防治結(jié)構(gòu)角度、顆粒級配三方面因素，探明各因素對碎屑流運動過程及沖擊性能的影響規(guī)律，為相關(guān)防治結(jié)構(gòu)設(shè)計與施工提供理論依據(jù)。

1 數(shù)值模擬方法

1.1 EDEM軟件介紹

EDEM作為一種DEM仿真商業(yè)軟件，最初被用于工業(yè)生產(chǎn)中的顆粒處理及其制造設(shè)備的生產(chǎn)過程的仿真和分析。用戶可利用EDEM輕松快速地創(chuàng)建顆粒實體的參數(shù)化模型，同時可以將力、材料特性和其他物理特性添加到EDEM中，形成完整的顆粒模型，為用戶建立個性化的模型處理環(huán)境。

因此，借助EDEM的顆粒工廠技術(shù)(Particle Factory)，可很好地模擬碎屑流的完整運動過程，并且該軟件具備高效的并行計算能力，可進行百萬級顆粒運算，能夠基本滿足碎屑流沖擊攔擋結(jié)構(gòu)的計算要求。

1.2 模型設(shè)計與參數(shù)選取

本文旨在通過數(shù)值模擬方法進行碎屑流沖擊流槽試驗模擬分析，進而開展碎屑流運動規(guī)律及其沖擊性能研究，因此，對比分析國內(nèi)外學者已經(jīng)開展的流槽試驗模型，綜合考慮其尺寸及參數(shù)完整性，最終選擇以Jiang et al.(2013)所完成的流槽試驗為參照進行數(shù)值模擬模型建立。其流槽試驗裝置如圖1所示，其中流槽整體長2.93im，寬0.3im，高0.35im，流槽坡度角α為40°，底部攔擋結(jié)構(gòu)角θ為50°，初始顆粒沉積區(qū)域長0.44im，寬0.3im，高0.150im。

圖1 流槽試驗裝置示意圖(Jiang et al.，2013)

圖2即本文參照上述參數(shù)建立的流槽試驗模型，主要包括槽底面、側(cè)墻、攔擋結(jié)構(gòu)、顆粒儲存盒、觸發(fā)門以及垂直與水平兩個方向的加壓面。設(shè)置加壓面是為了保證各組試驗的初始條件保持一致，即每組試驗開始前，均需通過垂直與水平加壓面對顆粒進行壓密。上述流槽試驗模型部件可由CAD導入，導入后需要對部件特性及部件運動進行指定，流槽試驗模型材料參數(shù)如表1所示。需要注意的是，該流槽試驗模型各部件除觸發(fā)門外均為靜態(tài)部件，即在模擬過程中不發(fā)生移動，因此只需要對觸發(fā)門進行運動指定，本文通過觸發(fā)門的線性旋轉(zhuǎn)模擬觸發(fā)門打開過程。

圖2 數(shù)值模擬模型組成

表1 數(shù)值模擬中材料參數(shù)取值

在顆粒級配選取方面，Jiang et al.(2013)所選取的物質(zhì)組成中粒徑為10～20imm的顆粒占比在95%左右(圖3)，已有研究表明對于粒徑10～20imm的顆粒占比在95%以上的情況，由于粒徑小于10imm的顆粒在初始顆粒樣品中的質(zhì)量百分比相對較低，忽略尺寸小于10imm的顆粒對計算結(jié)果的影響(Shen et al.，2018)，故本文直接選取10～20imm間的線性顆粒級配作為數(shù)值模擬中的物質(zhì)組成如圖3所示。

圖3 顆粒級配曲線

此外，固體顆粒參數(shù)及接觸參數(shù)如表1和表2所示，其中顆粒楊氏模量、泊松比、密度等參數(shù)主要根據(jù)顆粒介質(zhì)數(shù)值模擬中的常用值進行設(shè)置(Salciarini et al.，2010)。對于接觸參數(shù)，由于本文模擬的固體顆粒均為球形，而實際多為帶棱角且表面不光滑的粗糙顆粒，因此在Jiang et al.(2013)的試驗基礎(chǔ)上，本文通過增大靜摩擦系數(shù)，減小滾動摩擦系數(shù)來盡可能還原真實顆粒情況，具體取值如表2所示。

表2 數(shù)值模擬中有關(guān)接觸的系數(shù)選取

1.3 工況設(shè)計

本文運用控制變量法，選取對碎屑流運動過程及沖擊性能影響較大的3個因素作為自變量，即流槽坡度角α、底部攔擋結(jié)構(gòu)角θ以及顆粒級配，具體工況設(shè)計如表3所示。

表3 工況設(shè)計

1.4 模型驗證

為了驗證本文所建立模型的可靠度，本節(jié)將數(shù)值模擬結(jié)果與Jiang et al.(2013)的試驗數(shù)據(jù)進行了比較，圖4、圖5分別為Jiang et al.(2013)的試驗與數(shù)值模擬得到的泥石流運動過程圖像，可以看出兩組顆粒運動形態(tài)基本吻合，表明該模型可以很好地模擬顆粒運動過程，并且通過總結(jié)發(fā)現(xiàn)顆粒運動過程主要可分為如下4個階段：

圖4 流槽試驗中碎屑流運動過程照片(Jiang et al.，2013)

圖5 數(shù)值模擬中碎屑流運動過程圖像

(1)第1階段，顆粒儲存箱的觸發(fā)門被移除，顆粒開始崩塌(圖4a2和圖5b2)；

(2)第2階段，隨著顆粒的持續(xù)滑動，滑體后緣逐漸離開其原始位置，前部滑體轉(zhuǎn)變?yōu)楸?，而中段滑體擾動較小(圖4a3和圖5b3)；

(3)第3階段，滑體前端撞擊剛性攔擋結(jié)構(gòu)并發(fā)生反彈，同時，部分顆粒開始沉積在流槽末端的攔擋結(jié)構(gòu)處(圖4a4和圖5b4)；

(4)第4階段，越來越多的顆粒到達已沉積顆粒上，直到所有的顆粒運動停止(圖4a5和圖5b5)，最終顆粒沉積形態(tài)如圖4a6和圖5b6所示。

此外，圖6對比了試驗與數(shù)值模擬中的沖擊力時程曲線，基于離散單元法進行數(shù)值模擬時，當與攔擋結(jié)構(gòu)接觸的顆粒被困在相對長的力鏈中，該顆粒所產(chǎn)生的能量會顯著增加，進而產(chǎn)生力的峰值，同時導致沖擊力時程曲線發(fā)生明顯波動。因此，本節(jié)采用中值濾波方法對原始時程曲線進行處理，經(jīng)處理后的時程曲線與試驗結(jié)果吻合，即該模型可以很好地模擬顆粒流的沖擊性能。

圖6 流槽試驗與數(shù)值模擬中的沖擊力時程曲線

2 碎屑流運動規(guī)律分析

2.1 碎屑流基本運動規(guī)律

各工況下的數(shù)值模擬結(jié)果所表現(xiàn)出的碎屑流運動過程基本相似，以S1-2(α=40°，θ=50°)為例，運動過程可分為如下3個階段：

(1)第1階段(沖擊)：在顆粒撞擊攔擋結(jié)構(gòu)后，前端滑體形成薄層沉積體(圖7a～圖7b)，此時顆粒運動速度達到最大值；

圖7 碎屑流固體顆粒運動過程(S1-2)

(2)第2階段(爬升)：后端顆粒逐步向下運動，在前端既有沉積物上發(fā)生碰撞并逐步爬升，到達攔擋結(jié)構(gòu)的更高處(圖7c～圖7e)。顆粒的不斷爬升堆積使前端沉積高度逐步上升，而此時顆粒運動速度逐漸減小，靠近攔擋結(jié)構(gòu)的顆粒運動速度均趨于0im·s-1。此外，該階段的沉積形態(tài)具有一定的變化規(guī)律，即由起初的凹型(圖7c)變?yōu)橹本€型(圖7d)，再到最后的凸型(圖7e)；

(3)第3階段(堆積)：待所有顆粒完全靜止且穩(wěn)定后，其最終沉積形態(tài)如圖7f所示，由于上層流動顆粒的加入，以及顆粒向上的爬升運動，沉積形狀最終呈現(xiàn)為凸型。

圖8則進一步展示了碎屑流固體顆粒速度場的變化過程，其中紅色代表速度為3.08～3.96im·s-1，黃色代表速度為2.22～3.08im·s-1，綠色代表速度為1.35～2.22im·s-1，藍色代表速度為0.49～1.35im·s-1，灰色代表速度為0～0.49im·s-1。從圖9中可以發(fā)現(xiàn)靜態(tài)沉積區(qū)在不斷阻止顆粒的運動，導致顆粒減速向上爬升，同時，隨著沉積體積的增大，靜態(tài)沉積區(qū)逐漸向上游延伸。當既有沉積物坡度到達一個閾值后，隨后的顆粒無法越過沉積物(圖8a5)，只能沖擊并堆積在現(xiàn)有的沉積物上。

圖8 碎屑流固體顆粒速度場變化過程(S1-2)

圖9 不同流槽坡度角下碎屑流運動速度時程曲線(S1-1，S1-2，S1-3)

2.2 流槽坡度角α對碎屑流運動規(guī)律的影響

如2.1節(jié)中所述，各工況下運動規(guī)律大致相似，但從圖9中可以發(fā)現(xiàn)，隨著流槽坡度角的增大，固體顆粒越快達到速度峰值，且峰值越大，當流槽坡度角α由35°分別增加至40°和45°時，對應的速度峰值漲幅分別為8.37%和22.62%。此外，顆粒運動速度從峰值到停止運動所需時間隨著坡度的增大而減短，35°、40°和45°所對應的時間分別為1.50is、1.25is和1.10is。需要說明的是，為了保證各工況下初始顆粒運動狀態(tài)一致，t=0is至t=0.5is期間在料箱內(nèi)對顆粒進行預先壓密處理，因此該區(qū)間內(nèi)存在明顯的速度變化。

有效殘余力由最終沉積形狀的幾何形狀決定(Jiang et al.，2013)，在流槽下端形成的速度為0的三角形區(qū)域被廣大學者稱為dead zone，該區(qū)域的長度在很大程度上取決于坡度，同時也影響著下端結(jié)構(gòu)的平均受力。因此進一步觀察各工況下的沉積形態(tài)(圖10)可以發(fā)現(xiàn)，隨著流槽坡度的增大，沉積體高度越低而長度越寬。

圖10 碎屑流固體顆粒速度場變化過程(S1-1，S1-2，S1-3)

2.3 底部攔擋結(jié)構(gòu)角θ對泥石流運動規(guī)律的影響

本節(jié)對比分析了底部攔擋結(jié)構(gòu)角θ分別為30°、50°、70°和90°的工況，即S2-1，S3-1，S2-2和S2-3。數(shù)值模擬結(jié)果表明，在撞擊攔擋結(jié)構(gòu)前，各工況下固體顆粒的運動速度及運動規(guī)律大致相似，而撞擊后的顆粒堆積情況如圖11所示，其中a為側(cè)視圖，b為顆粒在攔擋結(jié)構(gòu)面上的分布。

圖11 碎屑流固體顆粒堆積形態(tài)(S2-1，S1-2，S2-2，S2-3)

圖11利用不同顏色表示不同粒徑的顆粒，其中紅色粒徑最大(18～20imm)，其次是黃色(16～18imm)、綠色(14～16imm)、藍色(12～14imm)、灰色(10～12imm)。可以發(fā)現(xiàn)隨著底部攔擋結(jié)構(gòu)角θ的增大，顆粒在攔擋結(jié)構(gòu)上的爬升距離逐漸降低，并且，當θ=90°時的固體顆粒的反序現(xiàn)象并不十分顯著，而隨著角度變緩反序現(xiàn)象愈發(fā)明顯，即細顆粒主要分布于沉積物的底部，粗顆粒分布于頂部。這是由于隨著固體顆粒爬升距離的增加，顆粒繼續(xù)前進，此時，細顆粒仍可沿著粗顆粒向下移動，而粗顆粒則只能繼續(xù)向前滑動，最終呈現(xiàn)為反序現(xiàn)象。

2.4 顆粒級配對碎屑流運動規(guī)律的影響

圖12為S3-2下不同粒徑顆粒的運動過程，其中2～5imm顯示為黃色，5～10imm顯示為藍色，10～15imm顯示為綠色，15～20imm顯示為紅色。當固體顆粒釋放后直至沉積停止運動后，可觀察到明顯的反序分離現(xiàn)象，即粗顆粒主要分布于頂部與前部，而細顆粒主要分布于底部與后部。這是由于大多數(shù)粗顆粒在運動過程相互碰撞，飛躍后分離至頂部與前部，而細顆粒由于比表面積較大，相應的能量損失較大，運動相對較慢，因此在粗顆?？p隙間穿過滑落至底部與后部。

圖12 不同粒徑顆粒的運動過程(S3-2)

圖13進一步對比了S3-1，S3-2，S3-3工況下固體顆粒運動速度變化曲線，分別對應了粗顆粒含量由低變高的情況，發(fā)現(xiàn)隨著粗顆粒含量的提升，固體顆粒峰值運動速度逐漸減小，表明細顆?？梢欢ǔ潭忍岣咚樾剂鞴腆w顆粒的運動速度。

圖13 不同顆粒級配下碎屑流運動速度時程曲線(S3-1，S3-2，S3-3)

3 碎屑流沖擊性能分析

3.1 流槽坡度角α對碎屑流沖擊性能的影響

圖14為不同流槽坡度角α下碎屑流固體顆粒對攔擋結(jié)構(gòu)的沖擊力時程曲線，觀察到流槽坡度角α越大，泥石流沖擊力則越快達到峰值，且峰值越大，同時后續(xù)也越快減小至殘余沖擊力，即沉積顆粒的重力分力穩(wěn)定值，這是由于重力加速度在沿流槽方向上的分量發(fā)生了改變，進而影響了沿流槽方向的沖擊力。

圖14 不同流槽坡度下碎屑流沖擊力時程曲線(S1-1，S1-2，S1-3)

3.2 底部攔擋結(jié)構(gòu)角θ對碎屑流沖擊性能的影響

圖15為不同底部攔擋結(jié)構(gòu)角θ下碎屑流固體顆粒對攔擋結(jié)構(gòu)的沖擊力時程曲線，發(fā)現(xiàn)各工況下，當泥石流到達攔擋結(jié)構(gòu)(t=1.4is)后，沖擊力迅速到達峰值，之后逐漸降低到殘余沖擊力，并且最大沖擊力隨攔擋結(jié)構(gòu)角θ的增大而增大，這是由于當攔擋結(jié)構(gòu)角θ增大后顆粒間的相互作用隨之降低，最終導致了顆粒對攔擋結(jié)構(gòu)沖擊力的提升。

圖15 不同攔擋結(jié)構(gòu)角下碎屑流沖擊力時程曲線(S2-1，S1-2，S2-2，S2-3)

3.3 顆粒級配對碎屑沖擊性能的影響

圖16為不同顆粒級配下碎屑流固體顆粒對攔擋結(jié)構(gòu)的沖擊力時程曲線，發(fā)現(xiàn)碎屑流最大沖擊力隨著粗顆粒占比的增大而顯著增大。其中：S3-3的粗顆粒占比最大，其沖擊力時程曲線的波動也更為劇烈，主要源于前端粗顆粒不斷撞擊攔擋結(jié)構(gòu)以及后續(xù)顆粒碰撞前緣堆積體而間接沖擊攔擋結(jié)構(gòu)的過程，相應地，該組泥石流沖擊力峰值最大；而S3-1中粗顆粒占比最小，主要以細顆粒為主，因此，顆粒間的摩擦耗能較為嚴重，此時泥石流對攔擋結(jié)構(gòu)的沖擊作用主要源于前緣堆積體的重力作用，同時細顆粒間的碰撞作用相對較小，導致相較于S3-3，其沖擊力時程曲線波動小且沖擊力峰值小。

圖16 不同顆粒級配下碎屑流沖擊力時程曲線(S3-1，S3-2，S3-3)

4 結(jié) 論

本文基于離散單元法，利用EDEM 2018軟件，通過數(shù)值模擬手段分析了碎屑流固體顆粒的運動過程及其與攔擋結(jié)構(gòu)間的相互作用，特別是碎屑流的沖擊性能。主要考慮了流槽坡度角、底部攔擋結(jié)構(gòu)角以及顆粒級配的影響，對碎屑流的運動規(guī)律及其沖擊性能進行了系統(tǒng)研究，得到以下主要結(jié)論：

(1)流槽坡度越大，碎屑流固體顆粒越快達到其速度與沖擊力峰值，且對應的速度與沖擊力的峰值越大。

(2)攔擋結(jié)構(gòu)角越大，與其相互作用的固體顆粒數(shù)量越多，因此碎屑流越快達到速度和沖擊力的峰值，同時也越快停止滑動，且對應的速度與沖擊力的峰值也越大。

(3)不同顆粒級配的碎屑流均呈現(xiàn)反序特性，即細顆粒主要分布于下部及后部，并且細顆粒的存在可顯著提高泥石流的運動速度，而粗顆粒對運動速度影響不大。

(4)粗顆粒占比越大，沖擊力主要源于前端粗顆粒對攔擋結(jié)構(gòu)的撞擊作用以及后續(xù)顆粒碰撞前緣堆積體而產(chǎn)生的能量，沖擊力時程曲線波動大且沖擊力峰值大，相反地，細顆粒占比越大，顆粒間的摩擦耗能越嚴重，沖擊力主要源于前緣堆積體的重力作用，其沖擊力時程曲線波動小且沖擊力峰值小。