基于沿程坐標積分模式顆粒流與結構物陣列相互作用的數值模擬1)

楊 肅 張會琴 余王昕 程鵬達 劉青泉 王曉亮 ,2)

* (北京理工大學宇航學院,北京 100081)

? (中國科學院力學研究所,北京 100190)

** (中國科學院大學工程科學學院,北京 100049)

引言

中國山地面積廣袤,超過2/3 的面積屬于山地和高原,各種大型災害時有發生[1].其中,以顆粒物質流動主導的大尺度災害如滑坡、泥石流,對人類的生命安全和各類基礎設施有著巨大的威脅[2-3].2010 年8 月8 日,甘肅舟曲縣發生特大規模山洪泥石流,造成1467 人遇難,298 人失蹤[4].2017 年6 月24 日,四川省阿壩州茂縣發生特大滑坡災害,造成83 人死亡,沖毀民房數十座,給當地民眾帶來了巨大損失[5].為了有效應對這類含大量顆粒物質的地質災害,通常在重要的聚居區或建筑物前設置一些防護結構,改變泥石流等災害的運動路徑、消耗能量降低沖擊作用,以達到防治或減輕災害的效果[6].近年來這類顆粒流動與結構物的相互作用越來越受到關注,主要集中在建筑物抵抗災害沖毀作用評估以及設置障礙物系統防災減災兩方面.因此,研究顆粒流與結構物相互作用的動力過程有助于防災減災,具有重要的工程意義.

顆粒材料的本構模型與通常的水體介質有顯著不同,導致顆粒流與結構物相互作用的模式和以水流為主的災害(如洪水)差別很大[7-8].顆粒流與結構物作用包含沖擊、繞流、爬升等行為,近10 年來顆粒流與結構物作用的動力過程吸引了大量學者進行研究[9-10].Pudasaini 等[11]進行了干顆粒流沖擊擋土墻的室內斜槽實驗,擋土墻會導致顆粒流從快速運動的超臨界薄層變化到厚度逐漸增加的靜止堆積體,最后的堆積面由顆粒材料的休止角決定.他們利用粒子圖像測速技術(PIV)進行觀測,分析了激波的形成和傳播過程、超臨界流體深度的演化、沖擊速度和動量的變化,并進行了初步的一維數值模擬.Pudasaini 和Kroener[12]進一步對快速密集顆粒流中的強激波進行了數值研究,預測得到了大部分Pudasaini等[11]的實驗結果.Teufelsbauer 等[13]發展了一種模擬干顆粒在斜面上流動的離散元模型(DEM),研究發現考慮旋轉約束的DEM 更準確地反映干顆粒物質的行為,同時還模擬了流通、沉積和撞擊過程.Zhou 等[14]利用DEM 研究了泥石流沖擊狹縫壩的爬升高度,并通過動量定理給出了解析表達式.Bi等[15]利用DEM 模擬研究了不同行間距和列間距的結構物陣列對顆粒流沖擊力的影響.孫新坡等[16]采用DEM 研究了都江堰汶川高速公路龍洞子隧道出口的崩塌體運動過程,對攔石墻進行了優化設計.Dai 等[17]建立了一個泥石流流過攔沙壩的光滑粒子流體動力學(SPH)模型,其中泥石流模型使用了黏性流體模型,攔沙壩被處理為彈性固體,并研究了文家溝和紅椿溝地區的泥石流沖擊攔沙壩的沖擊力演化過程.Gray 等[18]采用深度積分模型,計算得到了小尺度顆粒流快速流過四面體結構物形成的激波和真空結構.Cui 和Gray[19]研究了顆粒流流過圓柱體結構物所形成的弓形激波以及真空區等流態結構.

在顆粒流與多個結構物相互作用的研究方面,Juez 等[20]建立了一個類淺水波模型模擬顆粒流通過3 個半球結構物的過程,采用修正重力項表達陡峭地形效應,模擬得到的障礙物附近的顆粒流流態與實驗觀測基本吻合,但在結構物周圍附近的深度演化存在較大的偏差.張睿驍等[21-23]通過DEM對顆粒流和結構物相互作用的影響因素進行了研究,分別改變結構物尺寸、攔擋距離主要參數,得到不同形式的攔擋結構對顆粒流沖擊效應的影響.Fei 等[24]在模擬過程中采用了3 種不同的方式處理結構物,并比較了3 種方法在研究顆粒流與結構物相互作用問題上的優缺點.

綜上所示,當前對顆粒流與結構物的相互作用的研究主要集中在兩個方面: 第一是結構物作用下顆粒流動的流態及演化;第二是顆粒流對結構物的沖擊力作用.目前看起來,顆粒流和結構物相互作用的模型主要采用垂向積分模式和DEM 模型,垂向積分模式在處理復雜陡峭地形、障礙物處理以及顆粒流本構方面并沒有形成共識,由此在數值計算方法層面也帶來了一定的問題,而DEM 模型的研究對顆粒的尺寸、參數校核等方面提出的要求很多,造成當前的研究主要集中在實驗室尺度的顆粒流和結構物相互作用,對復雜結構陣列情形下的研究還比較缺乏.本文建立了基于沿程坐標積分模式的數值模型,模擬在陡峭地形上顆粒流與結構物的相互作用,特別是結構物作用下顆粒流的流態演化和堆積特征.在此基礎上深入研究了結構物陣列對顆粒流流態演化和堆積行為的影響,提出了一個新型的指標偏轉效率,和流通效率一起,定量評價了四面體結構物陣列對顆粒流阻礙作用和偏轉作用.

1 模型與驗證

1.1 控制方程

自然環境中的顆粒流災害在運動過程中,相對于流動方向,垂向的尺度很小,所以可以對其進行垂向積分,采取沿程坐標建立曲線坐標系,如圖1 所示.其中紅色四面體代表結構物,紫色曲線代表處在一個陡峭坡面上的顆粒流.綠色線條代表基礎地形,右側的水平流通區通過一個圓弧過渡區域與左側的斜坡區域連接在一起.(x,y,z)為建立的沿程坐標系.選取沿程坐標的主要目的是為了精確處理陡峭地形效應和顆粒流床面作用力的封閉問題.

圖1 陡峭地形下顆粒流與結構物陣列相互作用模型示意圖Fig.1 Illustration of granular flow past array of obstacles on steep curved terrain

建立基于沿程坐標系的控制方程[25-26]

其中,t代表時間,x,y是坐標系的兩個方向.h為顆粒物質的深度.u和v分別是顆粒流在x和y方向的深度平均速度.g=(gx,gy,gz)T為重力在坐標系3 個方向上的分量,zb為相對于曲面地形的局部偏移量,rx和ry為曲率半徑,δ 為床面摩擦角.

顆粒流是一種特殊的復雜流體,既表現出流體的性質,也有一定的固體特性,所以需要構建能夠準確描述這種特性的本構關系.從土力學角度出發,考慮應力的各向異性,kx和ky是由摩爾-庫倫定律確定的沿坡和跨坡方向的土壓力系數,表征運動土體水平和垂向壓力的比值,控制顆粒流沿流向和橫向的鋪展速率.對于主要運動方向為沿坡面向下的顆粒流.

方程主要參數是混合流體的內摩擦角 φ 和床面摩擦角 δ .其中,內摩擦角主要控制顆粒流體的變形程度,床面摩擦角主要控制其流動距離[25,27].

1.2 地形與結構物的參數化處理

在模擬過程中,顆粒流流體模型的準確性和結構物的處理方法都是問題的難點.在本文的研究中,采用參數化方案生成基礎地形和結構物陣列.

首先建立了一個如圖1 所示的沿程坐標系.基礎地形由3 部分組成,即藍色虛線代表的坡面區域,綠色虛線代表的圓弧過渡區域和藍色實線代表的水平流通區.其中,坡面區域與水平面的角度為 θ,圓弧過渡區與坡面區域和水平流通區的邊界分別為xl和xr.通過對這些關鍵參數,以及x方向和y方向的范圍和網格數等參數進行簡單設置,就可以得到基礎地形.

對于結構物(圖1 中紅色凸起區域),將其處理為一個相對于基礎地形的局部偏移量,用局部床面變化的場函數表示[10-11].同樣采用參數化的生成方案來生成結構物,只需輸入不同的結構物種類參數就可以生成長方體、四面體、圓柱體等形式的結構物,再對結構物的幾何參數,位置以及數量等參數進行設置,就可以得到不同形式的結構物陣列.生成初始顆粒堆積體所需要的參數與生成結構物陣列類似.此時,結構物被處理為地形,本文模型可以在統一框架下處理顆粒流在結構物作用下的反射、爬升以及繞射等動力過程.

1.3 數值方法

顆粒流的控制方程是帶源項的雙曲系統,這種非線性雙曲系統的典型特征是會產生激波.一階格式往往會把激波抹平,所以需要一種能夠很好地捕捉激波的高分辨率格式求解控制方程.除此以外,還需要正確處理干濕邊界條件,控制啟停行為,以達到在統一的歐拉框架下計算整個顆粒流與結構物相互作用動力過程.擴展Kurganov 和Petrova[28]在求解淺水方程時提出的激波捕捉格式(Kurganov 格式)處理對流項,采用算子分裂方法計算源項.

采用基于結構網格的有限體積法的網格劃分如圖2 所示.對于單元格 (i,j),定義場變量平均值為而方程1 左側部分的時間推進格式為

圖2 有限體積法網格劃分和數值通量示意圖Fig.2 Index of the mesh and numerical fluxes for finite volume methods

其中,采用MUSCL 方法插值,結合Minmod 限制器,構造單元格 (i,j) 邊界上的左右兩側的向量和使數值格式的空間精度在光滑區達到二階[28],分別是單元格 (i,j) 各邊界上最大和最小的波速.通量的計算方式與類似.

為了在統一計算框架下處理顆粒流的移動邊界或顆粒流前緣,仿照淺水波模型的術語,將顆粒流的前緣也稱之為“干濕邊界”,這里應用了一種改進的“干濕邊界”處理方法[29].把干單元格定義為單元格內顆粒平均堆積深度小于一個非常小深度hmin的單元.當所有相鄰單元的平均深度都小于hmin時,認為流動停止;當一些相鄰的單元內的平均深度大于hmin時,只計算質量守恒方程,速度從最深的相鄰單元外推.這樣可以避免非物理結果,比如負深度.先前關于快速顆粒流運動[29]和淺水波方程的數值求解[30]都表明,hmin值取為1 0-6時,界面前緣的捕捉較為理想,可以避免虛假振蕩,防止薄層水流出現超大速度.

采用Strang 分裂方法[31]對源項進行計算,由于源項不涉及微分運算,因此直接采用時間向前積分,即在對流算子推進的基礎上再進行源項直接積分更新動量和速度.為保障在單個時間步內不出現回流現象,單個單元的速度如果出現回流,則直接在本時間步的計算中將其設置為零.

1.4 模型驗證

為了驗證本文所建立的基于沿程坐標積分模式的數值模型,模擬Caviedes-Voullième 等[32]開展的典型顆粒流與結構物作用實驗,即顆粒流與一大二小共計3 個半球相互作用的實驗.首先將文獻[32]中的全局坐標轉變為本文的沿程坐標,得到計算區域為[0,1.88 m] × [0,1 m].初始顆粒堆積體為一個半徑為104 mm 的半球形堆積體,位于(162.4 mm,500 mm)處.大的半球形結構物的半徑為59 mm,位于(1.083 4 m,0.5 m)處,兩個小的半球形結構物的半徑為29 mm,分別位于(1.024 8 m,0.41 m)處和(1.024 8 m,0.59 m)處,如圖3 所示.主要地形相關參數和物理參數如表1 所示.

圖3 顆粒初始時刻及4 個固定測點位置示意圖,虛線之間為過渡區,虛線左側為斜坡,虛線右側為流通區Fig.3 Initial positions of granular flow and four probes,region in between the two dashed lines is transition area,the slope is to the left of the two dashed lines,and the runout zone is to the right of the two dashed lines

表1 “三半球實驗”模擬算例基礎地形生成關鍵參數數值表Table 1 Parameters involved in generating steep curved terrain of “three hemi-spheres experiment”

設置了5 組不同的網格密度對“顆粒流沖擊三半球問題”進行模擬,對數值模型進行網格無關性測試.提取了障礙物附近的兩個典型測點PSL 和PS0L(見圖3),在5 次模擬結果中的顆粒流深度并進行對比,結果如圖4 所示.當網格密度設置為50 × 50 和100 × 100 時,固體測點的時程曲線變化較大,并沒有收斂;而網格密度達到200 × 200 及以上時,計算結果達到收斂.在后續數值算例中,為保障收斂性,如無特殊說明,所有的算例中,網格數量默認為400 ×400.

圖4 不同網格數典型測點的顆粒流深度時程曲線Fig.4 Depth of granular flow over time at fixed measuring points under different grid densities

文獻[32] 并沒有給出床面摩擦角 δ 具體數值,本文對其進行了多次測試,發現當 δ=34°時,可以更好地描述顆粒堆積與鋪展的狀態,與實驗結果也更加匹配.而原實驗中,結構物所在區域的下游地形有一個逐漸彎曲向上的曲率,但由于那部分區域幾乎沒有顆粒流經過,也并不影響顆粒流與結構物相互作用的過程,所以本文的模型將其處理為水平面,并不會對最終結果造成顯著差異.

將Caviedes-Voullième 等[32]的實驗結果和基于全局坐標積分模式的數值模擬結果與本文的模擬結果進行了比對.本文對幾個典型時刻進行了細致的分析,如圖5 和圖6 所示.在t=460 ms 時,顆粒流前緣即將與較小半球發生接觸;在t=500 ms 時,顆粒流前緣與較大半球開始接觸,仔細觀察可以發現,Juez 等[20]的模擬結果在前緣部分有一個狹長的弧形區域集中了大量顆粒物質,而本文模型模擬結果與實驗數據則更為發散.在t=640 ms 時,顆粒流與結構物發生沖擊,形成了3 個強弱不等的激波.當t=740 ms 時,3 個激波向上游運動并開始接觸;在t=900 ms 時,3 個激波合并成新的波狀激波并向上游運動,隨后慢慢松弛,形成最終的堆積體.觀察圖5 結果,可以發現Juez 等[20]的數值模擬結果和本文的結果均較實驗結果略微落后.本文模型的模擬結果顯示,大量的顆粒物質被截留在結構物所處區域及上游區域,與實驗結果更為吻合,體現了結構物對顆粒流的截留作用和對下游的保護作用.

圖5 幾個典型時刻顆粒流深度分布文獻結果[20,32]與本文模擬結果對比.(a) 第1 列為實驗結果,第2 列為文獻中的模擬結果(圖片使用得到Elsevier 的許可),(b) 本文數值模型模擬結果Fig.5 Comparisons between the results from literature[20,32] and the simulation results of this paper at different times.(a) The first column is the experimental results,and the second column is simulation results from literature (permitted by Elsevier).(b) The numerical simulation results of this paper

圖6 幾個典型時刻顆粒流深度分布文獻結果[20,32]與本文模擬結果對比.(a) 第1 列為實驗結果,第2 列為文獻模擬結果(圖片使用得到Elsevier 的許可),(b) 本文數值模型模擬結果Fig.6 Comparisons between the results from literature [20,32] and the simulation results of this paper at different times.(a) The first column is the experimental results,and the second column is simulation results from literature (permitted by Elsevier).(b) The numerical simulation results of this paper

圖7 給出了4 個特定測點的顆粒流深度與實驗結果和基于全局坐標垂向模式的模擬結果的比較.4 個測點分別是PSL,PSR,PS0L 和PS0R,它們分別位于大球的肩部和側面,如圖3 所示.圖7 中的3 條曲線分別是多次實驗的平均深度時程演化[32]、Juez等[20]的模擬結果和本文的模擬結果.在大球肩部的兩個測點PSL 和PSR,本文模擬的峰值深度小于Juez 等[20]的模擬結果,更加貼近實驗結果,而最后達到的穩定深度基本和實驗結果吻合,大球肩部兩個測點處在形成穩定深度之前出現了一個峰值,這是由于大球肩部形成了較為顯著的激波結構,如圖5(b),而隨著顆粒物質補充的逐漸減少,激波會逐漸消失,最后在摩擦阻力的作用下逐漸達到穩定的堆積狀態.在大球兩個側部的兩個測點PS0L 和PS0R,已經基本遠離激波形成區域,不存在明顯的峰值,因此大球側部兩個測點的時程表現為逐漸增高直至穩定堆積狀態,本文模擬的穩定深度略高于實驗結果,但是模擬效果略好于Juez 等[20]的模型.

圖7 4 個測點的顆粒流深度變化時程曲線Fig.7 Depth of granular flow over time at the four probe positions

需要指出的是本文模型基于沿程坐標,能夠較好的表征地形的平均曲率效應,從而準確反映顆粒動力學過程在陡峭地形情形下受到離心和阻力作用,但是障礙物作為相對于曲面地形的局部偏差,在控制方程中以源項體現,在實際情況下障礙物附近的顆粒流動之三維效應會比較明顯,這是造成激波效應不明顯的大球側部測點計算結果相比實驗結果有偏差的主要原因,而三維效應也導致大球肩部測點處的深度峰值計算結果比實驗結果略高.Juez 等[20]的模擬主要基于全局坐標的垂向積分模式,雖然考慮了薄層顆粒流的約化重力效應,但是無法準確表征陡峭地形導致的離心和阻力作用[33],因此模擬結果不如本文模型.

總體來說,本文發展的基于沿程坐標的數值模型相比于基于全局坐標的垂向積分模式在表征復雜地形條件顆粒流與結構物相互作用問題上更為準確一些,數值模型能夠很好地模擬顆粒流與結構物作用中的反射、繞射以及爬升等現象,捕捉相互作用過程中形成的激波結構及其演化等過程.

2 顆粒流與陣列式結構物的相互作用

2.1 算例設置

為了減輕顆粒流災害的破壞力,一種經濟有效的減災措施是在上游地區設置結構物陣列對顆粒流進行調控[6,10].目前,工程界已經采用了在滑坡泥石流雪崩等災害性顆粒流的流經路徑上布置障礙物的方案防治和減輕災害[10,27,34],其設計準則主要依賴于工程經驗和室內小尺度實驗,現階段對該類陣列結構防護效果的數值模擬和定量評估還鮮有報道.本文主要研究一類典型的四面體結構物陣列[10,34]對顆粒流的防護效果,通過設置了不同行數的四面體結構物陣列,模擬分析顆粒流與結構物陣列之間的相互作用過程,探究顆粒流動過程和堆積狀態的改變.計算區域設置為[0,60 m] × [0,40 m],網格數量為400 × 400.選取一個典型斜坡和典型顆粒材料物理參數,見表2.在坐標(6 m,20 m)處設置一個半徑為2m 的半球形的顆粒堆積物模擬初始滑體.

表2 “四面體結構物陣列”系列算例基礎地形生成關鍵參數數值表Table 2 Parameters involved in generating steep curved terrain of “arrays of tetrahedral obstacles” cases

在坡面區域設置不同行數的四面體結構物陣列,每個四面體的高度為4.5 m,底面為面積是1 m2的正三角形,每行結構物由9 個完全相同的四面體組成,每個四面體的間距為4 m,第一行結構物位于x=13 m,不同行之間交錯排列,行間距為3 m.首先模擬了沒有結構物情況下初始顆粒堆積體沿斜坡的流動堆積過程,然后分別模擬1～ 5 行四面體結構物陣列作用下,顆粒流的流態演化和堆積過程.

2.2 流動過程的演化

初始顆粒堆積體在釋放后,由于重力和壓力梯度的作用,會開始變形并沿坡面運動,迅速演變成沿坡面向下的顆粒流.如果沒有結構物作用,在t=4 s時,如圖8(a)所示,大部分顆粒流已經到達平面流通區域.之后,由于沒有了沿斜坡向下的重力分量作用,流通到平面區域的顆粒因為摩擦作用會逐漸減速,而后續顆粒依然會以較高速度持續流入平面區域,所以在波的前緣會堆積大量顆粒,集中于一個半月牙型的狹窄區域,有激波結構存在.

圖8 t=4 s 時顆粒流與不同行數四面體結構物陣列相互作用的流態Fig.8 Flow patterns of granular flow facing obstacle array of tetrahedrons of different rows at t=4.0 s

當存在四面體結構物陣列時,結構物會與顆粒流發生顯著的相互作用而形成新的流動結構與堆積形態.當設置1 行四面體結構物時,如圖8(b)所示,顆粒在流動過程中正面沖擊到四面體后,由于兩個傾斜側面的存在,顆粒流會發生偏轉,分成兩道繼續向下游流動.四面體結構物的背部沒有顆粒物質,形成一個真空區.被分流后的顆粒流還會與被前排其他四面體結構物偏轉的顆粒流發生匯聚,相互沖擊會消耗大量能量.在t=4 s 時,有相當多的顆粒已經到達平面流通區.當設置2 行四面體結構物陣列時,如圖8(c)所示,相互作用過程更為復雜.顆粒沖擊到第一排四面體結構物后發生上文所述的流動行為,然后會繼續被第2 排的四面體結構物重新分流,再次因為偏轉作用被分成兩股向下游流動,被不同四面體分流的顆粒依然可能再次發生匯聚沖擊.由于顆粒流之間發生的多次匯聚沖擊以及被偏轉后流動路徑的變化,在t=4 s 時,大部分顆粒依然位于斜坡區域,僅有少量顆粒物質到達平面區域.隨著四面體結構物行數增加,如圖8(d)～ 圖8(f)所示,顆粒流之間匯聚沖擊的次數也會增加,流動路徑也會變長,但基本的相互作用模式與2 行結構物情況下(圖8(c))并無明顯區別.圖9 給出了5 行四面體障礙物作用時,快速顆粒流形成的典型流態特征,包括激波特征和偏轉特征.

圖9 t=4 s 時顆粒流與5 行四面體結構物陣列相互作用的局部流動信息Fig.9 Local flow information of granular flow facing obstacle array of tetrahedrons of 5 rows at 4.0 s

2.3 結構物誘導快速顆粒流激波結構特征

為深入考察顆粒流和結構物作用時激波的形成和特征,選取了5 行障礙物陣列中一個典型的四面體附近的流動為研究對象,考察3 個截面1～3 處的顆粒流深度變化特征,3 個截面位置如圖10 所示,同時選取3 半球算例(圖3)在中間截面處(y=0.5 m)的顆粒流深度變化特征作為比較對象,結果如圖11所示.圖11 表明高速顆粒流遇到結構物時會減速甚至堆積形成低速區,后續的快速顆粒流遇到前方低速的顆粒流動時來不及光滑過渡,會形成一種類似水躍的間斷結構,在顆粒流領域也被稱為激波[7,9,11-12,18,30],這種激波會耗散一部分顆粒流系統的能量,而障礙物陣列則形成系列的激波結構.這種激波結構和障礙物的陡峭程度有關,如圖11(a)表明半球附近障礙物接近于垂直墻體,形成了十分明顯的激波向后方運動,而圖11(b)～圖11(d)顯示,四面體結構物對顆粒流的降速和阻礙作用相比于接近墻體的半球體要弱,但依然形成了明顯的接近于間斷的激波結構向后方傳播.

圖10 5 排障礙物算例某單個結構物附近顆粒流深度(m)分布,3 條紅色虛線1,2 和3 為3 個考察截面Fig.10 Depth distribution (m) of granular flow around a typical obstacle in the case of five rows of obstacles.The three red dashed lines represent three sections for examination

圖11 障礙物附近典型快速顆粒流激波形成與演化Fig.11 Shock wave formation and evolution for granular flow around an obstacle

圖11 障礙物附近典型快速顆粒流激波形成與演化(續)Fig.11 Shock wave formation and evolution for granular flow around an obstacle (continued)

單一的障礙物可以在障礙物前緣形成弓形激波,耗散能量;在遇到障礙物陣列的情況下會形成系列的弓形激波,耗散更多的能量,如圖8(b)～圖8(f)、圖9(a)所示.因此障礙物陣列對顆粒流的作用包括激波的耗散作用和改變流向的偏轉作用.

2.4 堆積狀態

如圖12 所示,當各算例模擬到15 s 時,流動行為已經基本全部停止,顆粒物質以不同的形態堆積在平面流通區.在沒有結構物存在的情況下,顆粒物質會堆積成一個相對于y方向中心線完全對稱的堆積體,且只有一個堆積峰.

圖12 t =4.0 s 時顆粒流沖擊不同行數四面體結構物陣列后的堆積厚度分布Fig.12 Deposits of granular flow facing obstacle array of tetrahedrons of different rows at t=4.0 s

圖12 t =4.0 s 時顆粒流沖擊不同行數四面體結構物陣列后的堆積厚度分布 (續)Fig.12 Deposits of granular flow facing obstacle array of tetrahedrons of different rows at t=4.0 s (continued)

在有四面體結構物陣列存在的情況下,最后時刻的堆積狀態有了很大的改變.只有一行四面體結構物時,顆粒最終會堆積為一個大堆積體,其中有兩個堆積峰.觀察顆粒流動全過程,可以發現開始堆積時,這兩個堆積峰實際上是分散的兩部分堆積體,但由于偏轉次數以及流動距離等的限制,兩部分堆積體比較接近,隨著顆粒增多逐漸合并為一個堆積體.

當存在兩行四面體結構物時,最后顆粒會堆積為3 個堆積體,中間的堆積體略大.當初始顆粒堆積體向下游流下后被3 行四面體結構物陣列調控時,最終堆積為一大兩小3 個顆粒堆積體.其中中部較大的堆積體存在兩個堆積峰,兩側的堆積體比較小.中部堆積體兩個堆積峰的形成原因與一行結構物時的堆積狀態形成過程類似,依然是由于偏轉次數較少,大量顆粒堆積在中部,使得原本有被分為兩個堆積體趨勢的顆粒合并為一個.

顆粒流與4 行四面體結構物陣列相作用時,相較于3 行陣列多了若干偏轉機會,原來形成中部大堆積體的兩股顆粒流還會被偏轉一次,其一部分會相中部集中,另一部分會被分向兩側,所以最終會形成五個堆積峰,顆粒流質量分布情況也會比3 行結構物陣列的最終情況更加均勻.顆粒流與5 行四面體結構物陣列相作用時,有著更多的偏轉機會,同時,由于分匯次數的增加,顆粒形成的激波結構耗能增加,而且流動路徑會變得更長,所以大量顆粒在距離斜坡區域更近的位置停止運動,在水平流通區的顆粒百分比有所減小.

3 結果與討論

3.1 評價堆積狀態的兩個指標

為了準確描述顆粒流造成的災害影響范圍,判斷防護結構的效果,需要一些定量指標來直觀有效地進行評價.目前已有的一個評估指標是流通效率(runout efficiency,RE),它定義為初始顆粒堆積體流出距離與初始高度的比值,是一個評價滑坡災害影響范圍的重要指標[28,35].各算例在相同的初始顆粒堆積體條件下,RE值越大,說明流出距離越大,災害范圍越廣;RE值越小,說明流出距離越小,結構物防護效果越好.

這個指標顯然并不足以精確表述最后的堆積狀態.為了更準確地描述,從統計學角度出發,參照流通效率的定義,提出了一個新的無量綱量,即偏轉效率(deflection efficiency,DE),來表述顆粒流的偏轉效應,即

其中,h為堆積狀態各單元格內顆粒堆積深度,YC為堆積體的質心坐標,V為顆粒堆積體的體積,r0為初始顆粒堆積體的底面的半徑.DE可以在一定程度上表征防護結構對滑坡災害的偏轉作用,其值越大,說明顆粒流被分散的越廣.表3 列出了所有計算案例的流通效率和偏轉效率.

表3 顆粒流與多行四面體結構物陣列作用流通效率和偏轉效率Table 3 RE and DE of granular flow facing obstacle array of tetrahedrons

根據數據可以發現,由于結構物陣列的作用,RE值減小了,但減小幅度并不大,在僅有1 行四面體結構物時幾乎沒有變化.隨著行數增加,減小幅度約為6%～ 15%.以沒有障礙物的堆積結果為基準,DE的增大幅度更為顯著.在設置1 行結構物后,就可增大25%;設置2 行結構物,DE增大65%;在結構物陣列為3 行時,增大77%;在結構物陣列為4 行時,增大約111%;當陣列行數為5 行時,增大117%.

綜上所述,無量綱量偏轉效率DE指標可以在一定程度上描述陣列結構物對顆粒流的偏轉效應.結合耗散,結構物陣列對顆粒流有十分有效的綜合調控效果.因為本次模擬的是四面體型結構物,它的幾何形狀決定了對顆粒流最大的作用是使其偏轉,主要依靠增加流動路徑和不同股顆粒流之間的匯合沖擊進行耗能,因沖擊結構物導致的能量損失并不多,所以RE值的減小幅度并不如DE的增大幅度顯著.

3.2 指定區域的顆粒堆積深度變化

假設在平面流通區有一處需要受到保護的建筑物或人類聚居地,其坐標為(40 m,20 m),臨近圓弧過渡區域.本文提取了整個模擬過程中,這一點在上游無結構物情況下和上游存在不同行數結構物陣列以作保護的情況下的顆粒堆積深度的數據進行了對比,如圖13 所示.

圖13 (40 m,20 m)處顆粒流堆積深度時間歷程曲線Fig.13 Depth of granular flow over time at (40 m,20 m)

當上游沒有結構物存在時,顆粒流會在短時內迅速通過這一處,然后由于能量耗盡后,顆粒會在松弛堆積時,再次返回這一點位,使得圖像在7 s 左右出現了第二次急速升高的趨勢.當存在一行結構物時,大量顆粒被分流到兩側,首先到達的顆粒離開了中軸,所以這一點位的顆粒通過時間會延后,堆積峰值也會變小.當兩側的顆粒成堆時,會有部分顆粒堆積松散向四周滑落,兩堆顆粒都有滑向中部的顆粒,所以隨著時間推進,顆粒堆積深度會小幅增大.當存在兩行結構物時,正如3.1 中所述,第2 排結構物對顆粒的偏轉反而會使一部分顆粒向中部集中,所以堆積峰值反而會比無結構物時更大.當存在3 行結構物時,顆粒流首次通過時的深度峰值明顯下降.從上一節的內容可知,3 行結構物陣列并不能使顆粒流得到充分的偏轉分散,仍然有大部分顆粒流位于中軸附近.這些通過后的顆粒會慢慢堆積在平面流通區,同時對后續顆粒流起到阻礙作用,所以隨著時間增加,深度曲線會緩緩升高.當存在4 行結構物時,堆積峰值進一步降低,中軸附近的顆粒在還沒到達點位時就已耗盡能量而停止流動.當存在5 行結構物時,顆粒流得到充分偏轉,大量顆粒被分散到兩側.流經這一處的顆粒大大減少且通過時間延后,起到了良好的保護作用.所以,當存在一定形式的結構物陣列后,對于特定地區,可以起到有效的減災效果.

4 結論

針對顆粒流與結構物相互作用問題,本文使用數值模擬方法,通過對典型實驗[32]的模擬計算,驗證了模型具有良好的適用性.在此基礎上研究了顆粒流與不同行數四面體結構物陣列的相互作用.主要結論如下:

(1) 本文建立了基于沿程坐標積分模式建立的數值模型,可以較為準確的表征陡峭地形條件顆粒流與結構物相互作用的動力過程,如繞流、反射、爬升.本數值模型將流動過程、復雜地形、以及障礙物處理融為一體,基本具備了計算和評估陡峭地形條件顆粒流與復雜障礙物系統相互作用的流態和堆積等動力過程.

(2) 結構物對顆粒流的作用主要包括耗散作用和偏轉作用兩種模式.快速顆粒流遇到障礙物并在其前方誘導一個弓形激波,耗散顆粒的能量.同時,障礙物能夠改變顆粒流的流動方向,達到分隔和偏轉顆粒流的作用.

(3) 結構物陣列對顆粒流產生綜合的耗散和偏轉作用,通過多級作用產生系列的弓形激波耗散能量,通過偏轉作用分隔和偏轉顆粒流改變顆粒流的最終堆積形態,對下游地區產生顯著的防護效果.

顆粒流與結構物陣列相互作用的過程十分復雜,對于如何改變結構物陣列的參數以獲得對顆粒流更好的調控效果還有待進一步研究.此外如何將本數值模型與實際勘測的地質數據相結合也有待解決.未來希望能夠形成具有實際意義的大尺度山地災害中防護結構調控效果預測模型,為災害防治決策提供可靠的科學依據.