基于沿程坐標(biāo)積分模式顆粒流與結(jié)構(gòu)物陣列相互作用的數(shù)值模擬1)

楊 肅 張會(huì)琴 余王昕 程鵬達(dá) 劉青泉 王曉亮 ,2)

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

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

** (中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

引言

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

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

在顆粒流與多個(gè)結(jié)構(gòu)物相互作用的研究方面,Juez 等[20]建立了一個(gè)類(lèi)淺水波模型模擬顆粒流通過(guò)3 個(gè)半球結(jié)構(gòu)物的過(guò)程,采用修正重力項(xiàng)表達(dá)陡峭地形效應(yīng),模擬得到的障礙物附近的顆粒流流態(tài)與實(shí)驗(yàn)觀測(cè)基本吻合,但在結(jié)構(gòu)物周?chē)浇纳疃妊莼嬖谳^大的偏差.張睿驍?shù)萚21-23]通過(guò)DEM對(duì)顆粒流和結(jié)構(gòu)物相互作用的影響因素進(jìn)行了研究,分別改變結(jié)構(gòu)物尺寸、攔擋距離主要參數(shù),得到不同形式的攔擋結(jié)構(gòu)對(duì)顆粒流沖擊效應(yīng)的影響.Fei 等[24]在模擬過(guò)程中采用了3 種不同的方式處理結(jié)構(gòu)物,并比較了3 種方法在研究顆粒流與結(jié)構(gòu)物相互作用問(wèn)題上的優(yōu)缺點(diǎn).

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

1 模型與驗(yàn)證

1.1 控制方程

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

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

建立基于沿程坐標(biāo)系的控制方程[25-26]

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

顆粒流是一種特殊的復(fù)雜流體,既表現(xiàn)出流體的性質(zhì),也有一定的固體特性,所以需要構(gòu)建能夠準(zhǔn)確描述這種特性的本構(gòu)關(guān)系.從土力學(xué)角度出發(fā),考慮應(yīng)力的各向異性,kx和ky是由摩爾-庫(kù)倫定律確定的沿坡和跨坡方向的土壓力系數(shù),表征運(yùn)動(dòng)土體水平和垂向壓力的比值,控制顆粒流沿流向和橫向的鋪展速率.對(duì)于主要運(yùn)動(dòng)方向?yàn)檠仄旅嫦蛳碌念w粒流.

方程主要參數(shù)是混合流體的內(nèi)摩擦角 φ 和床面摩擦角 δ .其中,內(nèi)摩擦角主要控制顆粒流體的變形程度,床面摩擦角主要控制其流動(dòng)距離[25,27].

1.2 地形與結(jié)構(gòu)物的參數(shù)化處理

在模擬過(guò)程中,顆粒流流體模型的準(zhǔn)確性和結(jié)構(gòu)物的處理方法都是問(wèn)題的難點(diǎn).在本文的研究中,采用參數(shù)化方案生成基礎(chǔ)地形和結(jié)構(gòu)物陣列.

首先建立了一個(gè)如圖1 所示的沿程坐標(biāo)系.基礎(chǔ)地形由3 部分組成,即藍(lán)色虛線代表的坡面區(qū)域,綠色虛線代表的圓弧過(guò)渡區(qū)域和藍(lán)色實(shí)線代表的水平流通區(qū).其中,坡面區(qū)域與水平面的角度為 θ,圓弧過(guò)渡區(qū)與坡面區(qū)域和水平流通區(qū)的邊界分別為xl和xr.通過(guò)對(duì)這些關(guān)鍵參數(shù),以及x方向和y方向的范圍和網(wǎng)格數(shù)等參數(shù)進(jìn)行簡(jiǎn)單設(shè)置,就可以得到基礎(chǔ)地形.

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

1.3 數(shù)值方法

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

采用基于結(jié)構(gòu)網(wǎng)格的有限體積法的網(wǎng)格劃分如圖2 所示.對(duì)于單元格 (i,j),定義場(chǎng)變量平均值為而方程1 左側(cè)部分的時(shí)間推進(jìn)格式為

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

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

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

采用Strang 分裂方法[31]對(duì)源項(xiàng)進(jìn)行計(jì)算,由于源項(xiàng)不涉及微分運(yùn)算,因此直接采用時(shí)間向前積分,即在對(duì)流算子推進(jìn)的基礎(chǔ)上再進(jìn)行源項(xiàng)直接積分更新動(dòng)量和速度.為保障在單個(gè)時(shí)間步內(nèi)不出現(xiàn)回流現(xiàn)象,單個(gè)單元的速度如果出現(xiàn)回流,則直接在本時(shí)間步的計(jì)算中將其設(shè)置為零.

1.4 模型驗(yàn)證

為了驗(yàn)證本文所建立的基于沿程坐標(biāo)積分模式的數(shù)值模型,模擬Caviedes-Voullième 等[32]開(kāi)展的典型顆粒流與結(jié)構(gòu)物作用實(shí)驗(yàn),即顆粒流與一大二小共計(jì)3 個(gè)半球相互作用的實(shí)驗(yàn).首先將文獻(xiàn)[32]中的全局坐標(biāo)轉(zhuǎn)變?yōu)楸疚牡难爻套鴺?biāo),得到計(jì)算區(qū)域?yàn)閇0,1.88 m] × [0,1 m].初始顆粒堆積體為一個(gè)半徑為104 mm 的半球形堆積體,位于(162.4 mm,500 mm)處.大的半球形結(jié)構(gòu)物的半徑為59 mm,位于(1.083 4 m,0.5 m)處,兩個(gè)小的半球形結(jié)構(gòu)物的半徑為29 mm,分別位于(1.024 8 m,0.41 m)處和(1.024 8 m,0.59 m)處,如圖3 所示.主要地形相關(guān)參數(shù)和物理參數(shù)如表1 所示.

圖3 顆粒初始時(shí)刻及4 個(gè)固定測(cè)點(diǎn)位置示意圖,虛線之間為過(guò)渡區(qū),虛線左側(cè)為斜坡,虛線右側(cè)為流通區(qū)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 “三半球?qū)嶒?yàn)”模擬算例基礎(chǔ)地形生成關(guān)鍵參數(shù)數(shù)值表Table 1 Parameters involved in generating steep curved terrain of “three hemi-spheres experiment”

設(shè)置了5 組不同的網(wǎng)格密度對(duì)“顆粒流沖擊三半球問(wèn)題”進(jìn)行模擬,對(duì)數(shù)值模型進(jìn)行網(wǎng)格無(wú)關(guān)性測(cè)試.提取了障礙物附近的兩個(gè)典型測(cè)點(diǎn)PSL 和PS0L(見(jiàn)圖3),在5 次模擬結(jié)果中的顆粒流深度并進(jìn)行對(duì)比,結(jié)果如圖4 所示.當(dāng)網(wǎng)格密度設(shè)置為50 × 50 和100 × 100 時(shí),固體測(cè)點(diǎn)的時(shí)程曲線變化較大,并沒(méi)有收斂;而網(wǎng)格密度達(dá)到200 × 200 及以上時(shí),計(jì)算結(jié)果達(dá)到收斂.在后續(xù)數(shù)值算例中,為保障收斂性,如無(wú)特殊說(shuō)明,所有的算例中,網(wǎng)格數(shù)量默認(rèn)為400 ×400.

圖4 不同網(wǎng)格數(shù)典型測(cè)點(diǎn)的顆粒流深度時(shí)程曲線Fig.4 Depth of granular flow over time at fixed measuring points under different grid densities

文獻(xiàn)[32] 并沒(méi)有給出床面摩擦角 δ 具體數(shù)值,本文對(duì)其進(jìn)行了多次測(cè)試,發(fā)現(xiàn)當(dāng) δ=34°時(shí),可以更好地描述顆粒堆積與鋪展的狀態(tài),與實(shí)驗(yàn)結(jié)果也更加匹配.而原實(shí)驗(yàn)中,結(jié)構(gòu)物所在區(qū)域的下游地形有一個(gè)逐漸彎曲向上的曲率,但由于那部分區(qū)域幾乎沒(méi)有顆粒流經(jīng)過(guò),也并不影響顆粒流與結(jié)構(gòu)物相互作用的過(guò)程,所以本文的模型將其處理為水平面,并不會(huì)對(duì)最終結(jié)果造成顯著差異.

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

圖5 幾個(gè)典型時(shí)刻顆粒流深度分布文獻(xiàn)結(jié)果[20,32]與本文模擬結(jié)果對(duì)比.(a) 第1 列為實(shí)驗(yàn)結(jié)果,第2 列為文獻(xiàn)中的模擬結(jié)果(圖片使用得到Elsevier 的許可),(b) 本文數(shù)值模型模擬結(jié)果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 幾個(gè)典型時(shí)刻顆粒流深度分布文獻(xiàn)結(jié)果[20,32]與本文模擬結(jié)果對(duì)比.(a) 第1 列為實(shí)驗(yàn)結(jié)果,第2 列為文獻(xiàn)模擬結(jié)果(圖片使用得到Elsevier 的許可),(b) 本文數(shù)值模型模擬結(jié)果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 個(gè)特定測(cè)點(diǎn)的顆粒流深度與實(shí)驗(yàn)結(jié)果和基于全局坐標(biāo)垂向模式的模擬結(jié)果的比較.4 個(gè)測(cè)點(diǎn)分別是PSL,PSR,PS0L 和PS0R,它們分別位于大球的肩部和側(cè)面,如圖3 所示.圖7 中的3 條曲線分別是多次實(shí)驗(yàn)的平均深度時(shí)程演化[32]、Juez等[20]的模擬結(jié)果和本文的模擬結(jié)果.在大球肩部的兩個(gè)測(cè)點(diǎn)PSL 和PSR,本文模擬的峰值深度小于Juez 等[20]的模擬結(jié)果,更加貼近實(shí)驗(yàn)結(jié)果,而最后達(dá)到的穩(wěn)定深度基本和實(shí)驗(yàn)結(jié)果吻合,大球肩部?jī)蓚€(gè)測(cè)點(diǎn)處在形成穩(wěn)定深度之前出現(xiàn)了一個(gè)峰值,這是由于大球肩部形成了較為顯著的激波結(jié)構(gòu),如圖5(b),而隨著顆粒物質(zhì)補(bǔ)充的逐漸減少,激波會(huì)逐漸消失,最后在摩擦阻力的作用下逐漸達(dá)到穩(wěn)定的堆積狀態(tài).在大球兩個(gè)側(cè)部的兩個(gè)測(cè)點(diǎn)PS0L 和PS0R,已經(jīng)基本遠(yuǎn)離激波形成區(qū)域,不存在明顯的峰值,因此大球側(cè)部?jī)蓚€(gè)測(cè)點(diǎn)的時(shí)程表現(xiàn)為逐漸增高直至穩(wěn)定堆積狀態(tài),本文模擬的穩(wěn)定深度略高于實(shí)驗(yàn)結(jié)果,但是模擬效果略好于Juez 等[20]的模型.

圖7 4 個(gè)測(cè)點(diǎn)的顆粒流深度變化時(shí)程曲線Fig.7 Depth of granular flow over time at the four probe positions

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

總體來(lái)說(shuō),本文發(fā)展的基于沿程坐標(biāo)的數(shù)值模型相比于基于全局坐標(biāo)的垂向積分模式在表征復(fù)雜地形條件顆粒流與結(jié)構(gòu)物相互作用問(wèn)題上更為準(zhǔn)確一些,數(shù)值模型能夠很好地模擬顆粒流與結(jié)構(gòu)物作用中的反射、繞射以及爬升等現(xiàn)象,捕捉相互作用過(guò)程中形成的激波結(jié)構(gòu)及其演化等過(guò)程.

2 顆粒流與陣列式結(jié)構(gòu)物的相互作用

2.1 算例設(shè)置

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

表2 “四面體結(jié)構(gòu)物陣列”系列算例基礎(chǔ)地形生成關(guān)鍵參數(shù)數(shù)值表Table 2 Parameters involved in generating steep curved terrain of “arrays of tetrahedral obstacles” cases

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

2.2 流動(dòng)過(guò)程的演化

初始顆粒堆積體在釋放后,由于重力和壓力梯度的作用,會(huì)開(kāi)始變形并沿坡面運(yùn)動(dòng),迅速演變成沿坡面向下的顆粒流.如果沒(méi)有結(jié)構(gòu)物作用,在t=4 s時(shí),如圖8(a)所示,大部分顆粒流已經(jīng)到達(dá)平面流通區(qū)域.之后,由于沒(méi)有了沿斜坡向下的重力分量作用,流通到平面區(qū)域的顆粒因?yàn)槟Σ磷饔脮?huì)逐漸減速,而后續(xù)顆粒依然會(huì)以較高速度持續(xù)流入平面區(qū)域,所以在波的前緣會(huì)堆積大量顆粒,集中于一個(gè)半月牙型的狹窄區(qū)域,有激波結(jié)構(gòu)存在.

圖8 t=4 s 時(shí)顆粒流與不同行數(shù)四面體結(jié)構(gòu)物陣列相互作用的流態(tài)Fig.8 Flow patterns of granular flow facing obstacle array of tetrahedrons of different rows at t=4.0 s

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

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

2.3 結(jié)構(gòu)物誘導(dǎo)快速顆粒流激波結(jié)構(gòu)特征

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

圖10 5 排障礙物算例某單個(gè)結(jié)構(gòu)物附近顆粒流深度(m)分布,3 條紅色虛線1,2 和3 為3 個(gè)考察截面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 障礙物附近典型快速顆粒流激波形成與演化(續(xù))Fig.11 Shock wave formation and evolution for granular flow around an obstacle (continued)

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

2.4 堆積狀態(tài)

如圖12 所示,當(dāng)各算例模擬到15 s 時(shí),流動(dòng)行為已經(jīng)基本全部停止,顆粒物質(zhì)以不同的形態(tài)堆積在平面流通區(qū).在沒(méi)有結(jié)構(gòu)物存在的情況下,顆粒物質(zhì)會(huì)堆積成一個(gè)相對(duì)于y方向中心線完全對(duì)稱(chēng)的堆積體,且只有一個(gè)堆積峰.

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

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

在有四面體結(jié)構(gòu)物陣列存在的情況下,最后時(shí)刻的堆積狀態(tài)有了很大的改變.只有一行四面體結(jié)構(gòu)物時(shí),顆粒最終會(huì)堆積為一個(gè)大堆積體,其中有兩個(gè)堆積峰.觀察顆粒流動(dòng)全過(guò)程,可以發(fā)現(xiàn)開(kāi)始堆積時(shí),這兩個(gè)堆積峰實(shí)際上是分散的兩部分堆積體,但由于偏轉(zhuǎn)次數(shù)以及流動(dòng)距離等的限制,兩部分堆積體比較接近,隨著顆粒增多逐漸合并為一個(gè)堆積體.

當(dāng)存在兩行四面體結(jié)構(gòu)物時(shí),最后顆粒會(huì)堆積為3 個(gè)堆積體,中間的堆積體略大.當(dāng)初始顆粒堆積體向下游流下后被3 行四面體結(jié)構(gòu)物陣列調(diào)控時(shí),最終堆積為一大兩小3 個(gè)顆粒堆積體.其中中部較大的堆積體存在兩個(gè)堆積峰,兩側(cè)的堆積體比較小.中部堆積體兩個(gè)堆積峰的形成原因與一行結(jié)構(gòu)物時(shí)的堆積狀態(tài)形成過(guò)程類(lèi)似,依然是由于偏轉(zhuǎn)次數(shù)較少,大量顆粒堆積在中部,使得原本有被分為兩個(gè)堆積體趨勢(shì)的顆粒合并為一個(gè).

顆粒流與4 行四面體結(jié)構(gòu)物陣列相作用時(shí),相較于3 行陣列多了若干偏轉(zhuǎn)機(jī)會(huì),原來(lái)形成中部大堆積體的兩股顆粒流還會(huì)被偏轉(zhuǎn)一次,其一部分會(huì)相中部集中,另一部分會(huì)被分向兩側(cè),所以最終會(huì)形成五個(gè)堆積峰,顆粒流質(zhì)量分布情況也會(huì)比3 行結(jié)構(gòu)物陣列的最終情況更加均勻.顆粒流與5 行四面體結(jié)構(gòu)物陣列相作用時(shí),有著更多的偏轉(zhuǎn)機(jī)會(huì),同時(shí),由于分匯次數(shù)的增加,顆粒形成的激波結(jié)構(gòu)耗能增加,而且流動(dòng)路徑會(huì)變得更長(zhǎng),所以大量顆粒在距離斜坡區(qū)域更近的位置停止運(yùn)動(dòng),在水平流通區(qū)的顆粒百分比有所減小.

3 結(jié)果與討論

3.1 評(píng)價(jià)堆積狀態(tài)的兩個(gè)指標(biāo)

為了準(zhǔn)確描述顆粒流造成的災(zāi)害影響范圍,判斷防護(hù)結(jié)構(gòu)的效果,需要一些定量指標(biāo)來(lái)直觀有效地進(jìn)行評(píng)價(jià).目前已有的一個(gè)評(píng)估指標(biāo)是流通效率(runout efficiency,RE),它定義為初始顆粒堆積體流出距離與初始高度的比值,是一個(gè)評(píng)價(jià)滑坡災(zāi)害影響范圍的重要指標(biāo)[28,35].各算例在相同的初始顆粒堆積體條件下,RE值越大,說(shuō)明流出距離越大,災(zāi)害范圍越廣;RE值越小,說(shuō)明流出距離越小,結(jié)構(gòu)物防護(hù)效果越好.

這個(gè)指標(biāo)顯然并不足以精確表述最后的堆積狀態(tài).為了更準(zhǔn)確地描述,從統(tǒng)計(jì)學(xué)角度出發(fā),參照流通效率的定義,提出了一個(gè)新的無(wú)量綱量,即偏轉(zhuǎn)效率(deflection efficiency,DE),來(lái)表述顆粒流的偏轉(zhuǎn)效應(yīng),即

其中,h為堆積狀態(tài)各單元格內(nèi)顆粒堆積深度,YC為堆積體的質(zhì)心坐標(biāo),V為顆粒堆積體的體積,r0為初始顆粒堆積體的底面的半徑.DE可以在一定程度上表征防護(hù)結(jié)構(gòu)對(duì)滑坡災(zāi)害的偏轉(zhuǎn)作用,其值越大,說(shuō)明顆粒流被分散的越廣.表3 列出了所有計(jì)算案例的流通效率和偏轉(zhuǎn)效率.

表3 顆粒流與多行四面體結(jié)構(gòu)物陣列作用流通效率和偏轉(zhuǎn)效率Table 3 RE and DE of granular flow facing obstacle array of tetrahedrons

根據(jù)數(shù)據(jù)可以發(fā)現(xiàn),由于結(jié)構(gòu)物陣列的作用,RE值減小了,但減小幅度并不大,在僅有1 行四面體結(jié)構(gòu)物時(shí)幾乎沒(méi)有變化.隨著行數(shù)增加,減小幅度約為6%～ 15%.以沒(méi)有障礙物的堆積結(jié)果為基準(zhǔn),DE的增大幅度更為顯著.在設(shè)置1 行結(jié)構(gòu)物后,就可增大25%;設(shè)置2 行結(jié)構(gòu)物,DE增大65%;在結(jié)構(gòu)物陣列為3 行時(shí),增大77%;在結(jié)構(gòu)物陣列為4 行時(shí),增大約111%;當(dāng)陣列行數(shù)為5 行時(shí),增大117%.

綜上所述,無(wú)量綱量偏轉(zhuǎn)效率DE指標(biāo)可以在一定程度上描述陣列結(jié)構(gòu)物對(duì)顆粒流的偏轉(zhuǎn)效應(yīng).結(jié)合耗散,結(jié)構(gòu)物陣列對(duì)顆粒流有十分有效的綜合調(diào)控效果.因?yàn)楸敬文M的是四面體型結(jié)構(gòu)物,它的幾何形狀決定了對(duì)顆粒流最大的作用是使其偏轉(zhuǎn),主要依靠增加流動(dòng)路徑和不同股顆粒流之間的匯合沖擊進(jìn)行耗能,因沖擊結(jié)構(gòu)物導(dǎo)致的能量損失并不多,所以RE值的減小幅度并不如DE的增大幅度顯著.

3.2 指定區(qū)域的顆粒堆積深度變化

假設(shè)在平面流通區(qū)有一處需要受到保護(hù)的建筑物或人類(lèi)聚居地,其坐標(biāo)為(40 m,20 m),臨近圓弧過(guò)渡區(qū)域.本文提取了整個(gè)模擬過(guò)程中,這一點(diǎn)在上游無(wú)結(jié)構(gòu)物情況下和上游存在不同行數(shù)結(jié)構(gòu)物陣列以作保護(hù)的情況下的顆粒堆積深度的數(shù)據(jù)進(jìn)行了對(duì)比,如圖13 所示.

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

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

4 結(jié)論

針對(duì)顆粒流與結(jié)構(gòu)物相互作用問(wèn)題,本文使用數(shù)值模擬方法,通過(guò)對(duì)典型實(shí)驗(yàn)[32]的模擬計(jì)算,驗(yàn)證了模型具有良好的適用性.在此基礎(chǔ)上研究了顆粒流與不同行數(shù)四面體結(jié)構(gòu)物陣列的相互作用.主要結(jié)論如下:

(1) 本文建立了基于沿程坐標(biāo)積分模式建立的數(shù)值模型,可以較為準(zhǔn)確的表征陡峭地形條件顆粒流與結(jié)構(gòu)物相互作用的動(dòng)力過(guò)程,如繞流、反射、爬升.本數(shù)值模型將流動(dòng)過(guò)程、復(fù)雜地形、以及障礙物處理融為一體,基本具備了計(jì)算和評(píng)估陡峭地形條件顆粒流與復(fù)雜障礙物系統(tǒng)相互作用的流態(tài)和堆積等動(dòng)力過(guò)程.

(2) 結(jié)構(gòu)物對(duì)顆粒流的作用主要包括耗散作用和偏轉(zhuǎn)作用兩種模式.快速顆粒流遇到障礙物并在其前方誘導(dǎo)一個(gè)弓形激波,耗散顆粒的能量.同時(shí),障礙物能夠改變顆粒流的流動(dòng)方向,達(dá)到分隔和偏轉(zhuǎn)顆粒流的作用.

(3) 結(jié)構(gòu)物陣列對(duì)顆粒流產(chǎn)生綜合的耗散和偏轉(zhuǎn)作用,通過(guò)多級(jí)作用產(chǎn)生系列的弓形激波耗散能量,通過(guò)偏轉(zhuǎn)作用分隔和偏轉(zhuǎn)顆粒流改變顆粒流的最終堆積形態(tài),對(duì)下游地區(qū)產(chǎn)生顯著的防護(hù)效果.

顆粒流與結(jié)構(gòu)物陣列相互作用的過(guò)程十分復(fù)雜,對(duì)于如何改變結(jié)構(gòu)物陣列的參數(shù)以獲得對(duì)顆粒流更好的調(diào)控效果還有待進(jìn)一步研究.此外如何將本數(shù)值模型與實(shí)際勘測(cè)的地質(zhì)數(shù)據(jù)相結(jié)合也有待解決.未來(lái)希望能夠形成具有實(shí)際意義的大尺度山地災(zāi)害中防護(hù)結(jié)構(gòu)調(diào)控效果預(yù)測(cè)模型,為災(zāi)害防治決策提供可靠的科學(xué)依據(jù).