基于SPH-DEM的泥石流沖擊剛性防護結構的動力過程研究

趙德博,郝夢潔,熊 昊

(1. 深圳大學 土木與交通工程學院,廣東 深圳 518060; 2. 深圳大學 濱海城市韌性基礎設施教育部重點實驗室,廣東 深圳 518060)

0 引言

泥石流作為一種山區地帶多發的地質災害,包含泥砂、塊石的固液兩相流體,呈黏性層流或稀性紊流狀態[1]。由于其往往突然暴發、來勢兇猛、破壞力強,泥石流對人類的生命財產安全構成巨大威脅。我國泥石流幾乎遍布每個省、市和自治區,據自然資源部統計,僅2022年上半年,特別是5月底以來,南方多省遭遇多輪強降雨,地質災害多發群發,特別是浙江、福建、江西、湖南、廣東、廣西6省(區)發生4000多處崩塌滑坡泥石流,已超去年全年的數量[2]。2019年6月27日20時30分,四川省金川縣曾達鄉曾達溝(流域面積125.53 km2)暴發特大型泥石流災害,造成房屋、基礎設施、農業、公益設施損失約1.51億元[3]。泥石流災害后果極為嚴重,因此對于泥石流易發區域迫切需要采取治理措施。沿著泥石流路徑建造防護結構是目前泥石流治理采用的主要工程措施。防護結構主要包括剛性[4]和柔性[5-6]2種,它的設計可以改變泥石流的運動狀態,減弱甚至避免泥石流對下游造成破壞。然而目前在工程實踐中,主要還是采用經驗方法對防護結構進行設計與施工。究其原因,主要是缺乏對泥石流運動特點以及泥石流與防護結構相互作用機理的研究。

許多研究人員都通過現場實測、室內水槽試驗和數值模擬[7-9]等手段對泥石流對防護結構的沖擊開展過探究。劉道川等[10]通過小型水槽試驗開展了系列泥石流的沖擊試驗研究了兩相泥石流沖擊力的時空分布特征。而NG等[11]將泥石流簡化為均一的流體,借助水槽試驗研究了玻璃球和砂子分別作用于剛性防護結構的沖擊特性和機理,揭示了玻璃球作用下的爬升機制和砂顆粒作用下的堆積機制。此外,其他學者利用水槽試驗裝置研究了不同組成成分的泥石流與防護結構的作用機制,發現泥石流的堆積形態、爬升高度和泥石流的組成密切相關[12-13]。相比于原型場試驗,小尺度水槽試驗存在明顯的尺寸效應,無法完全還原泥石流的動力學特征。但大尺度的原型試驗成本較高,工作量大,不確定因素多,不適合泥石流對防護結構沖擊機理的探究。因此,近年來數值模擬在巖土工程方面[13-15],尤其是泥石流運動機理方面得到廣泛的應用[16]。甘建軍等[8]采用RAMMAS3D數值軟件準確模擬了泥石流的運動變形過程并確定了其影響區域。HE等[17]利用光滑粒子流體動力學(SPH)方法模擬了干顆粒流和剛性防護結構的作用,而CUI等[18]利用離散元法研究了不同粒徑的顆粒流撞擊防護系統的過程,并與試驗結果進行比對,取得了較好的效果。CALVETTI等[19]則從微觀角度解釋了干顆粒流對防護結構的沖擊。但以上數值模擬方法都傾向于將泥石流模擬為單一相的流體,這種簡化限制了對泥石流這種多相流體的深入理解。

因此,本研究采用SPH-DEM耦合方法開展水槽試驗的數值模擬,其中SPH模擬泥石流中的黏性流體,DEM模擬泥石流中的石塊,探究含石塊的泥石流的運動過程及其對剛性防護結構的沖擊效應,并研究了不同水槽傾角對泥石流的運動形態以及對剛性防護結構的沖擊力的影響,為剛性防護結構強度的工程設計標準提供參考。

1 SPH-DEM數值模擬方法

1.1 SPH數值方法

光滑粒子流體動力學(smooth particle hydrodynamics,SPH)是一種純拉格朗日無網格粒子方法,是當前數值計算的熱點話題之一。該方法基于連續介質假設以一系列任意分布的粒子來模擬流體運動規律,其中每一個粒子都有獨立的影響域和插值域,計算過程中粒子的速度、壓強和位置等性質以及它們的梯度分布都是根據插值域內粒子構造的核近似函數插值得到的。

在SPH方法中,位于r處的流體粒子特性的函數可由以下的積分表達式定義:

(1)

式中:h為光滑長度,控制積分域的大小;Ω為支持域;r′為支持域內其他粒子的位置;W為核函數,它不僅決定了SPH近似的內值精度,而且與數值穩定性有關。本文采用Wendland核函數[20]進行平滑,其表達式為:

(2)

式中:αD在二維問題中等于7/(4πh2),在三維問題中等于21/(16πh3);q=|r-r′|/h。

1.2 DEM數值方法

離散單元法(discrete element method,DEM)最早是由CUNDALL[21]于1971年提出的,已被廣泛應用于顆粒材料[22]、巖土工程[23]等領域。DEM的基本思想是利用一系列離散的粒子來模擬固體介質,其中每個粒子作為一個單元,具有固定的大小和形狀,單元之間可以在小范圍內相互擠壓重疊從而產生相互作用力。當單元i與單元j接觸時,根據牛頓第二定律,單元i的運動方程可以寫為:

(3)

(4)

式中:mi為流體粒子i的質量;u為單元i的運動速度;g為重力加速度;Fn,ij和Ft,ij分別是單元j對單元i的法向和切向作用力;Ii和ωi分別是單元i的轉動慣量和角速度;而Tt,ij為Fn,ij和Ft,ij對單元i的力矩。

為了簡化接觸模型,單元均用球形粒子來表示,那么對于發生碰撞的球形粒子,粒子間的法向接觸力如式(5)所示:

Fn=knδn+cnvn

(5)

式中:kn為法向剛度系數;δn為法向疊合量;cn為法向黏性阻尼系數;vn為法向相對速度。

粒子間的切向力大小與相接觸的粒子間的相對速度有關,滿足庫侖定律,如式(6)所示:

Ft=min(ktvt,μFn)

(6)

式中:kt為切向剛度系數;vt為切向相對速度。

1.3 SPH-DEM耦合方法

SPH和DEM的耦合計算實際上是將流體作用于固體顆粒表面上的作用力進行積分得到流體對固體的作用力,而固體顆粒表面單元內插到與其相鄰的流體粒子的SPH方程中,從而可以計算固體對流體的作用力。

在SPH-DEM耦合方法中,流體粒子的動量方程式為:

(7)

固體DEM顆粒所受的作用力主要包括固體顆粒間的接觸力、潤滑力、流體施加的力(浮力和阻力)和重力等,其控制方程式為:

(8)

(9)

2 SPH-DEM方法驗證

為了驗證本文采用的SPH-DEM耦合方法在模擬流體與固體之間復雜相互作用的可靠性,模擬了單個非旋轉球體的入水測試,并結合ARISTOFF等[24]開展的室內試驗和XU等[25]進行的模擬校驗結果。

試驗模型的裝置如圖1所示,水箱的長為0.20 m,寬為0.20 m,高為0.14 m,水深0.11 m。直徑為25.4 mm的球體以2.17 m/s的初始速度與水面相切釋放。數值模擬中使用的球體和水的參數如表1所示。模擬時考慮重力,重力加速度為9.81 m/s2。SPH流體粒子初始間距dp為球體半徑的1/10,采用Verlet時間積分算法,模擬1 s內小球進入水中的運動狀態。

圖1 球體入水的數值模型Fig. 1 Numerical model for water-entry of a sphere表1 球體入水模擬的關鍵參數Table 1 Parameters used in the simulation for water-entry of a sphere材料參數取值球體密度/(kg/m3)860楊氏模量/(MPa)100泊松比0.3水密度/(kg/m3)1000黏度/(×10-6 Pa?s)1

不同時刻小球進入水中的運動過程如圖2所示,與ARISTOFF等[24]的試驗在下降趨勢、水面形狀等方面一致。球體與水面接觸并開始向下運動,沖擊并擠壓水,致使球體周圍的水流向側向運動,并在球體上方形成一個空腔。隨著小球繼續向下運動,空腔塌陷,形成向上的水射流。而球體由于密度較小,速度逐漸減小為0,隨后在浮力的作用下被抬升。球體的整個運動軌跡都是沿一條垂直直線,沒有任何的旋轉。這是因為球體在耦合過程中只受來自垂直方向的力,這與理論情況也是相符的。

圖2 小球入水過程模擬圖Fig. 2 Simulation snapshots of the sphere entering water

球體入水深度隨時間的變化情況如圖3所示,并與試驗數據、理論數據和XU等[25]的模擬結果做了對比,發現在0.07 s之前,隨著時間的不斷增大,球體在水中的深度逐漸變大,但深度的增大率,即球體的速度卻逐漸變小。模擬效果與理論結果較為接近,且由于所采用的粒子間距大于XU等[25]的模擬但結果相差不大,模擬效率大大提升。由此可見,本文采用的SPH-DEM耦合方法具有可行性,能夠較好地模擬流體與固體之間的相互作用。

圖3 球體入水的深度隨時間的演變Fig. 3 Evolution of the depth of sphere into the water with time

考慮到多數泥石流中的流體均為黏性流體,本文進一步研究了球體進入黏性流體的運動,并提取其深度,如圖4所示?？梢钥吹皆陴ば粤黧w(黏度=1)中,球體的深度隨時間變化得慢,速度減小得快,在0.8 s左右,速度基本為0。因此,耦合模型模擬黏性流體與固體的相互作用依舊有效且可靠。

3 泥石流沖擊剛性防護結構

長期以來,泥石流頻發區防護結構的設計一直是人們關注的問題。目前許多學者已經針對泥石流與防護結構的作用過程開展了大量試驗研究。其中CHOI等[26]利用小型水槽裝置研究了不同傾角下干性泥石流沖擊下,防護結構的動態受荷過程,揭示了水槽傾角與沖擊力對防護結構設計的重要性。本文采用SPH-DEM的耦合方法模擬水槽試驗中包含石塊的泥石流在水槽內的流動以及對剛性防護結構的沖擊,在此基礎上改變水槽傾角、石塊的粒徑、泥石流中流體的體積,從而分析了不同條件下防護結構所受沖擊力以及其上水面高程的演變。

3.1 計算模型

水槽試驗模擬的數值模型尺寸如圖5所示。水槽的長、寬和高分別為2.3、0.6、0.7 m,防護結構距離料箱門1.6 m,料箱的容量為0.168 m3。防護結構的寬度、高度和厚度分別為600、700、5 mm。水槽底部與水平面的傾角為15°。

圖5 水槽試驗模擬的模型示意圖Fig. 5 Schematic diagram of the model for the simulation of the flume experiments表2 水槽試驗模擬的關鍵參數Table 2 Parameters used in the simulation of the flume experiment材料參數取值流體密度/(kg/m3)1000黏度/(Pa·s)1玻璃珠密度/(kg/m3)2000楊氏模量/MPa200泊松比0.2動態摩擦系數0.2恢復系數0.6

在初始狀態時,泥石流位于水槽頂部,由固液兩相組成。其中由SPH模擬泥石流中的黏性流體,流體的體積為0.072 m3,流體的黏度為1 Pa·s,流體粒子之間的間距為0.005 m,流體總粒子數為573934,而采用直徑為4 cm 的玻璃珠模擬泥石流中的石塊,本模擬中包含27個玻璃珠,即模擬的泥石流的固體體積分數為1.26%。整個模擬過程中考慮重力,重力加速度被設置為9.81 m/s2。模擬需要用到的其它參數如表2所示。

3.2 泥石流運動過程分析

泥石流在水槽中運動過程的模擬圖如圖6所示,SPH粒子的顏色代表流體的流速。在模擬之初,泥石流被生成在水槽頂部的料箱中,料箱門可以有效地阻止泥石流的運動。當料箱中的泥石流達到穩定狀態(T=1.5 s)后,料箱門向上抽出,模擬泥石流的啟動。在重力的作用下,兩相流沿著水槽底部運動。隨著運動速度的不斷增大,泥石流中的流體前部逐漸靠近剛性防護結構,當其一接觸到防護結構,流體的速度迅速減小,并沿著防護結構向上爬升(T=2.5 s),在其動能減少為0后,流體達到最大爬升高度,開始向后回滾(T=2.8 s),一部分堆積在防護結構的底部,形成一個穩定的堆積區,一部分與后續泥石流相撞,使其速度減少(T=3.1 s),而后續的泥石流不會直接撞擊到防護結構,而是直接匯入堆積體,堆積區的大小逐漸增大,直至達到穩定,而泥石流中的玻璃珠由于與水槽底部的摩擦以及與流體的相互作用的影響,運動速度要小于流體,故其到達防護結構的時間更晚,在流體已經形成堆積體時,玻璃珠會沿著堆積體向上運動(T=5.0 s)。

圖6 泥石流的運動過程Fig. 6 Runout process of the debris flow

為了進一步量化分析泥石流對防護結構的沖擊效應,模擬中還提取了泥石流中的黏性流體對防護結構的沖擊力隨時間的演變,如圖8中的綠色曲線所示。在t=2.36 s時,流體前端接觸到剛性防護結構底部,動態沖擊力由0迅速增大,流體沿著防護結構向上爬升時,沖擊力增加速度放緩,隨后隨著流體爬升到最大高度后開始回流以及后續泥石流的沖擊,沖擊力存在波動并在2.8 s左右達到峰值。之后的泥石流由于被堆積在防護結構后部的流體阻擋,并不能直接接觸到防護結構,故對防護結構的沖擊力逐漸減弱并趨于平緩,造成這種現象的原因主要是泥石流體積是一定的,防護結構后的沖擊力流量不再隨著時間繼續增加。而穩定前曲線仍存在的微小波動,是由于防護結構后堆積體的回流以及后續流體匯入堆積體對防護結構存在擠壓間接造成的。

泥石流與防護結構的相互作用受多種因素的共同影響,包括泥石流的屬性(如泥石流的體積、固體體積分數、密度和流體黏度等)、防護結構的位置、坡度以及水槽傾角等,這些都是設計防護結構需要考慮的重要因素。因此,本研究基于以上水槽試驗模型進行多次模擬,探究不同影響因素下泥石流對防護結構的沖擊效應。

3.3 水槽傾角的影響

為了比較不同水槽傾角對泥石流與防護結構相互作用的影響,分別設置θ=0°、8°、15°、24°、30°進行模擬。圖7為在模擬最終時刻t=6.0 s時,不同水槽傾角下泥石流的最終形態,對比發現水槽傾角越大,泥石流形成的堆積區體積就越大(θ=30°時除外),且只有θ=30°時,泥石流撞擊防護結構發生了溢流。這是因為水槽傾角越大,泥石流的動能越大,前端接觸到防護結構的泥石流越不易發生回流。同時,可以觀察到傾角越大,懸浮在流體上的玻璃珠數量越少,這是因為流體對顆粒向上的浮力減小導致的。

圖7 不同水槽傾角下泥石流的最終形態Fig. 7 Final shape of debris flow at different inclinations of the flume

通過比較不同水槽傾角下,泥石流對防護結構的沖擊力隨時間的變化情況,進一步分析了水平傾角對于泥石流對防護結構的沖擊效應的影響,如圖8所示。顯然,沖擊力的變化趨勢基本相似,先迅速增加,達到峰值沖擊力后再慢慢下降,并逐漸趨于穩定。但θ>15°時,沖擊力曲線會出現2個明顯的波峰,且后者為峰值沖擊力,這種現象的發生是因為第一個波峰主要是由流體前部對防護結構的撞擊導致,隨著流體的回流沖擊力減弱,后部流體的動能仍較大,超過前部流體對其的阻礙,故出現了第二個波峰。而θ<15°時,后部流體的動能小于前部流體回流的阻礙,所以只存在一個波峰。同時,注意到若以θ=0°時的峰值沖擊力50.26 N為基準,則θ=8°時峰值沖擊力為其的2.59倍,隨著傾角的增大,其倍數依次為4.75、8.88和11.45,傾角越大其峰值沖擊力的增長率越大。且通過圖像可以看出,在趨于穩定前,傾角越大,沖擊力波動越大,這是因為流體的反復沖擊與回流。穩定后的沖擊力差異也很大,這是由于防護結構底部堆積體的體積不一致導致的。提取防護結構中心軸處水面高程隨時間的變化情況如圖9所示。由圖9可知,傾角越大,到達防護結構的時間越短,并且沿防護結構爬升的高度越大,最后穩定時的水面高程也越大,但穩定后水面高程的增加率是隨傾角的增大而逐漸減小的。同時水面高程達到穩定的時間隨傾角的增大而變長,即泥石流與防護結構的動態作用時間隨傾角的增大而增大。

圖8 不同水槽傾角下沖擊力隨時間的變化Fig. 8 Relationship of impact force-time at different inclinations of the flume

圖9 不同水槽傾角下水面高程隨時間的變化Fig. 9 Relationship of water surface elevation-time at different inclinations of the flume

綜上所述,改變水槽傾角對于泥石流對防護結構的沖擊效應的影響是顯著的,不僅沖擊力會成倍地變化,當防護結構的高度低于泥石流沿防護結構爬升的高度時,還會有溢流現象發生。因此在工程實踐中,設計防護結構前考慮泥石流易發處的地形傾角,模擬泥石流對防護結構的沖擊特性是非常必要的。

3.4 泥石流體積的影響

在自然災害中,泥石流的暴發規模不同,其造成的危害程度是不一致的,那么需要采取的預防措施便也不同。本文為了探究不同規模泥石流對防護結構沖擊的影響,將泥石流中流體部分的體積V分別設置為0.024、0.048、0.072、0.096、0.120 m3進行模擬。圖10比較t=6.0 s時刻不同流體體積的泥石流運動形態,可見流體體積的增多導致固體體積分數的減小,會促進流體與顆粒間的相互作用,從而一定程度上增加玻璃珠的動能。圖11進一步比較不同流體體積下沖擊力隨時間的變化,可以看出,流體體積越大其接觸到防護結構所需的時間會越短,達到穩定沖擊力的時間越長。同時,流體規模越大,其對防護結構的峰值沖擊力和靜態沖擊力都越大。隨著流體體積的增大,其峰值沖擊力增長率也逐漸變大,相較V=0.024 m3分別增長了0.79倍、1.64倍、3.20倍和4.56倍,但是穩定下來的靜態沖擊力基本與體積的變化成正比。結合圖12水面高程隨時間的變化可以看到,雖然水面高程隨流體體積的增大而增高,但流體達到最大水面高程所需時間基本相同。而且最大水面高程的增長率隨體積的增大表現為先增加后減小的趨勢,分別為27%、47%、33%和21%;穩定后水面高程的增長率隨體積的增大則依次遞減的,分別為76%、21%、15%和13%。也就是說,在實際工程中為能更加經濟環保地設計防護結構,可以考慮在保持一定高度時,適當增強防護結構的強度以阻擋較大規模的泥石流。

圖10 不同流體體積下泥石流的最終形態Fig. 10 Final morphology of debris flow with different fluid volumes

圖11 不同流體體積下沖擊力隨時間的變化Fig. 11 Relationship of impact force-time with different fluid volumes

圖12 不同流體體積下水面高程隨時間的變化Fig. 12 Relationship of water surface elevation-time with different fluid volumes

3.5 石塊粒徑的影響

此外,本研究還考慮了其他條件一定的情況下,泥石流中石塊粒徑的改變對泥石流與防護結構相互作用的影響。在確保石塊總體積不變時,模擬球形石塊半徑R分別為10、20、30 mm時,泥石流對防護結構的沖擊。由圖13中t=6.0 s時刻不同石塊粒徑下泥石流的運動形態可以看出,相同體積下,石塊粒徑越小,其與流體的接觸面積越大,對流體的拖曳力的反作用力越強,則越易阻礙流體的運動。而圖14展現的防護結構所受沖擊力的變化則表明,由于固體體積分數較小,石塊粒徑的改變對峰值沖擊力的影響不大,但是石塊粒徑越大,泥石流對防護結構的靜態沖擊力越易波動,這可能是由于大塊石對流體回流的阻礙更弱導致的。同時,由圖15可知,石塊粒徑越大,流體的最大水面高程越大,但石塊粒徑對穩定后的水面高程基本沒有影響。因此,對于含有大塊石的泥石流,在設計防護結構時可以適當提高防護結構的高度來對泥石流起到更好的阻擋效果。

圖13 不同石塊粒徑下泥石流的最終形態Fig. 13 Final morphology of the debris flow with different rock radii

圖14 不同石塊粒徑下沖擊力隨時間的變化Fig. 14 Relationship of impact force-time with different rock radii

圖15 不同石塊粒徑下水面高程隨時間的變化Fig. 15 Relationship of water surface elevation-time with different rock radii

4 結論

本文利用球體入水試驗驗證了SPH-DEM耦合方法在模擬流固兩相相互作用時的可靠性?；隍炞C方法進一步建立了水槽試驗中含石塊的泥石流沖擊剛性防護結構的模型,并研究了泥石流的運動過程及對剛性防護結構的沖擊效應,分析了不同水槽傾角、不同泥石流流體體積和不同石塊粒徑的情況下沖擊力和水面高程的變化,得出以下結論:

1)通過與試驗、理論結果的對比分析表明,采用耦合的SPH-DEM方法可以較好地模擬固液兩相的相互作用,其中SPH模擬流體,而DEM模擬固體。而且該方法對模擬含石塊泥石流的運動過程依舊可靠。

2)含石塊泥石流與防護結構的相互作用可以分為3個階段:沖擊、爬升和堆積。當水槽傾角較大時,也可能會有溢流發生。泥石流中的黏性流體對防護結構的沖擊力先迅速上升達到峰值沖擊力后逐漸下降并趨于穩定,穩定后的沖擊力不為0。而其中的石塊由于與流體的相互作用以及水槽底部的摩擦的影響,運動速度小于流體且會沿著堆積體向上攀爬。

3)水槽傾角、流體體積和石塊粒徑對于泥石流與防護結構的相互作用影響顯著。隨著傾角由0°依次增長為8°、15°、24°和30°,峰值沖擊力增長倍數依次為2.59、4.75、8.88和11.45,傾角越大其峰值沖擊力的增長率越大。但穩定后水面高程的增加率是隨傾角的增大而逐漸減小的;隨著流體體積的增大,其峰值沖擊力增長率也逐漸變大,相較V=0.024 m3分別增長了0.79倍、1.64倍、3.20倍和4.56倍,但是穩定下來的靜態沖擊力基本與體積的增長成正比;石塊粒徑的改變對沖擊力的影響不大,但會影響峰值水面高程的形成,即石塊粒徑越大,峰值水面高程越高。