樊赟赟,王思敬,王恩志

（1.東北大學(xué) 資源與土木工程學(xué)院,沈陽 110819;2.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室,北京 100084）

雪崩、滑坡、滾石和泥石流等都是常見的山地災(zāi)害,他們非常危險(xiǎn)且具有很強(qiáng)的破壞性。近年來在人類活動(dòng)和氣候變化等因素的影響下這些山地災(zāi)害發(fā)生的頻率更加頻繁,爆發(fā)規(guī)模也越來越大。為了減小災(zāi)害的危害性除了加強(qiáng)災(zāi)害預(yù)報(bào),還應(yīng)該研究災(zāi)害發(fā)生和發(fā)展的規(guī)律,判斷災(zāi)害可能的運(yùn)動(dòng)路徑及致災(zāi)范圍,并找到經(jīng)濟(jì)而有效的防治方法,這些也將為山區(qū)移民規(guī)劃提供重要的參考和依據(jù)[1-4]。

數(shù)值模擬是研究山地災(zāi)害的一種有效手段,已有許多國內(nèi)外的學(xué)者開展了這方面的研究,并取得了一定的成就[1-17]。雪崩、滑坡、滾石和泥石流等山地災(zāi)害的發(fā)展過程常常表現(xiàn)出一定的顆粒流動(dòng)特征,因而可以將這些山地災(zāi)害的發(fā)展過程以顆粒流動(dòng)來概化[1-4,10-17]。為討論方便論文將上述山地災(zāi)害發(fā)展過程統(tǒng)一視為顆粒流動(dòng),通過研究顆粒流動(dòng)過程的特點(diǎn)來加深對(duì)這些山地災(zāi)害的認(rèn)識(shí),并通過研究障礙物設(shè)置對(duì)顆粒流動(dòng)過程的作用和影響來探討如何選用有效而可靠的災(zāi)害防治措施方案。

設(shè)置障礙物是控制顆粒流動(dòng)過程的有效手段。研究障礙物對(duì)顆粒流動(dòng)過程的影響,一方面可以通過設(shè)置障礙物改變顆粒流動(dòng)的路徑,限制其影響范圍,減小其對(duì)人類居住環(huán)境和基礎(chǔ)設(shè)施的破壞,以達(dá)到防災(zāi)和減災(zāi)的目的,另一方面也可以充分利用已有障礙物,通過優(yōu)化設(shè)計(jì),得到既合理又經(jīng)濟(jì)的防治方案。

近年來國內(nèi)外學(xué)者針對(duì)顆粒流動(dòng)開展了大量的研究,已經(jīng)形成了多種理論和方法。其中Savage和Hutter較早提出了一定質(zhì)量顆粒的流動(dòng)理論具有較大的影響[1-2],理論考慮了由于速度梯度不同而引起的地壓系數(shù)的變化,能夠在正交曲線坐標(biāo)系下模擬一定質(zhì)量顆粒的流動(dòng),因而得到了比較廣泛的關(guān)注和研究[1-4,10-16]。同時(shí)也有學(xué)者比較關(guān)注障礙物對(duì)顆粒流動(dòng)過程的影響,他們通過實(shí)驗(yàn)和數(shù)值模擬討論障礙物的不同形式和高度等對(duì)顆粒流動(dòng)的影響,并討論其中的物理過程[3,8,16-17]。

但上述學(xué)者均未深入討論障礙物的不同設(shè)置對(duì)顆粒流動(dòng)過程的影響。事實(shí)上,在實(shí)際的工程應(yīng)用中障礙物應(yīng)因循具體的地形和地貌布置,其空間設(shè)置將對(duì)顆粒流動(dòng)過程產(chǎn)生較大的影響。該文將應(yīng)用SH理論,采用Roe格式的有限體積離散方法對(duì)流經(jīng)障礙物的顆粒流動(dòng)過程進(jìn)行數(shù)值模擬,并討論障礙物的設(shè)置對(duì)顆粒流動(dòng)過程的影響。

1 模型方程[1-2]

SH理論建立了如圖1所示的正交曲線坐標(biāo)系統(tǒng)[3],其中的虛線平面為曲線參考平面,其水平的坡角為ζ,zb為相對(duì)參考曲面的底高,理論描述了顆粒材料在垂直于參考曲面上的二維運(yùn)動(dòng)。

圖1 正交曲線坐標(biāo)系統(tǒng)

式中的無量綱變量h為顆粒流深度,v=（u,v）為深度平均的平行于坡面的平面速度,而u和v則分別是沿坡向和垂直于坡向的平面速度分量,k=-?ζ/?x為參考平面的曲率,δ為床面摩擦角,K x和K y分別為沿坡向和垂直坡向的地壓力系數(shù)。

方程中的無量綱變量能夠通過映射變換成原有的物理變量,如式（4）所示：

其中g(shù)表示重力加速度,H表示顆粒的特征深度,L表示特征長度,R為溝道沿坡向的特征半徑,SH理論還定義了2個(gè)比值ε=H/L和λ=L/R,他們一般都比較小。

SH理論的一個(gè)關(guān)鍵特點(diǎn)就是引入了CM準(zhǔn)則的地壓力系數(shù),其根據(jù)速度梯度條件的不同可以表示為式（5）-（8）,這些關(guān)系式在x-y平面剪切比x-z和y-z平面剪切小的條件下是正確的。

式中φ是顆粒材料內(nèi)摩擦角,而K yact/pass表示主動(dòng)或被動(dòng)地壓力系數(shù)。

對(duì)于不同類型的山地災(zāi)害,在模化為顆粒流動(dòng)時(shí),顆粒材料主要由不同的內(nèi)摩擦角和床面摩擦角來定義。

2 顆粒流經(jīng)障礙物的過程模擬

2.1 數(shù)值方法及計(jì)算條件

顆粒流動(dòng)過程具有復(fù)雜性,在大梯度物理參數(shù)條件下和顆粒減速堆積或遇到障礙時(shí),流動(dòng)狀態(tài)往往會(huì)發(fā)生變化,因而需要采用穩(wěn)定而有效的數(shù)值方法進(jìn)行計(jì)算模擬。隨著現(xiàn)代數(shù)值方法的發(fā)展,許多優(yōu)秀的數(shù)值方法被應(yīng)用于SH模型方程的計(jì)算中。

論文將以 Roe格式的近似 Riem ann解為基礎(chǔ)[18],采用有限體積方式對(duì)模型方程進(jìn)行數(shù)值離散,在非結(jié)構(gòu)網(wǎng)格上采用 M USCL（Monotonic Upwind Scheme for Conservation Law s）線性重構(gòu)方法對(duì)空間進(jìn)行線性重構(gòu)[19],用 Rung-Kutta方法[20]獲得具有時(shí)空二階精度的格式,對(duì)界面插值進(jìn)行限制以避免出現(xiàn)假振[21],并對(duì)與速度有關(guān)的摩擦阻力項(xiàng)用半隱式的方法離散[22]。

為便于比較分析,選用文獻(xiàn)[3]所用的計(jì)算模型條件,研究障礙物的不同設(shè)置對(duì)顆粒流動(dòng)過程的影響。

如圖2所示,假定在半橢球形容器下有一定質(zhì)量的顆粒材料。顆粒材料在40°的坡面上被突然釋放,經(jīng)過一段光滑的曲面運(yùn)動(dòng)到水平面上。計(jì)算范圍取為沿坡向方向x=[0,30],垂直于坡向方向y=[-10,10]。曲面的坡角可以表示為：

其中橢圓中心坐標(biāo)（x0,y0）=（5,0）,橢圓長軸a=4,橢圓短軸b=2,顆粒材料的最大初始高度hinimax=1。計(jì)算模型中的參數(shù)取為：顆粒內(nèi)摩擦角 φ=35°,床面摩擦角 δ=30°,特征比值 λ=1,ε=1。

假定顆粒粒徑相對(duì)障礙物足夠小,則當(dāng)顆粒連續(xù)的流經(jīng)障礙物后,將形成不同的堆積,而顆粒粒徑的差別則綜合體現(xiàn)在內(nèi)摩擦角和床面摩擦角中。文獻(xiàn)[3]在以上的初始條件下計(jì)算并討論了障礙物的不同高度對(duì)顆粒流動(dòng)過程的影響,研究發(fā)現(xiàn)只有當(dāng)障礙物足夠高時(shí)才能在堆積區(qū)形成無顆粒堆積的“安全區(qū)域”。因而,該文選取能夠形成“安全區(qū)域”的障礙物高度H=5進(jìn)行計(jì)算和討論,障礙物則選擇具有代表性的四面體障礙物。

圖2 計(jì)算模型條件示意圖

2.2 顆粒流經(jīng)朝向上坡障礙物的過程模擬

設(shè)置關(guān)于坡面對(duì)稱的四面體障礙物,障礙物的底面是邊長為4的正三角形,令它的一個(gè)頂點(diǎn)朝向上坡,該頂點(diǎn)的對(duì)邊位于x=17（如圖2中實(shí)線所示）。在計(jì)算中,不考慮障礙物的變形作用,通過改變地形條件實(shí)現(xiàn)四面體障礙物對(duì)顆粒流動(dòng)過程的影響。對(duì)該條件下的顆粒流動(dòng)過程進(jìn)行計(jì)算模擬,得到了在不同時(shí)刻下的顆粒流深度等值線圖,如圖3所示（等值線間距為0.05）。

從圖3結(jié)果中可以看到,顆粒材料在重力作用下開始向x和y方向流動(dòng),且沿坡向方向流動(dòng)較快（t=3）。當(dāng)顆粒到達(dá)障礙物時(shí),由于障礙物的作用,顆粒的流動(dòng)方向發(fā)生了改變,顆粒開始分別向2個(gè)方向流動(dòng)（t=6）,并且在重力作用下繼續(xù)加速流動(dòng)（t=9）。當(dāng)顆粒到達(dá)水平面,前部的流動(dòng)開始減速形成堆積,但尾部的流動(dòng)仍在加速（t=12～18）,直到尾部也到達(dá)堆積區(qū)（t=21）。顆粒材料最終在水平面上堆積形成了兩個(gè)對(duì)稱的顆粒堆（t=24）,由圖中可見,由于障礙物的作用,在對(duì)稱的顆粒堆中間形成了一個(gè)沒有顆粒堆積的“安全區(qū)域”。由此可見,障礙物改變了顆粒流動(dòng)的路徑,并影響了其最終的堆積形態(tài)。在正確的設(shè)置條件下形成了沒有顆粒流動(dòng)和堆積的“安全區(qū)域”,其可以看做是為避免災(zāi)害破壞的實(shí)際可行的應(yīng)用區(qū)域,這對(duì)災(zāi)害防護(hù)的設(shè)計(jì)具有重要的意義。

圖3 顆粒流經(jīng)朝向上坡障礙物的過程

值得一提的是,雖然所用的數(shù)值方法不同,但圖3的模擬結(jié)果與文獻(xiàn)[3]中的相關(guān)結(jié)果相同,這也驗(yàn)證了本文數(shù)值方法的有效性。

2.3 顆粒流經(jīng)朝向下坡障礙物的過程模擬

為了比較研究障礙物的不同角度設(shè)置對(duì)顆粒流動(dòng)過程的影響,將朝向上坡的障礙物沿垂線旋轉(zhuǎn)180°得到了朝向下坡的障礙物（如圖2中虛線所示）。對(duì)顆粒流經(jīng)朝向下坡障礙物的過程進(jìn)行計(jì)算模擬,得到了顆粒流動(dòng)發(fā)展過程中不同時(shí)刻的顆粒流深度等值線圖,如圖4所示（等值線間距為0.05）。

圖4 顆粒流經(jīng)朝向下坡障礙物的過程

從圖4結(jié)果中可以看到,顆粒材料在遇到障礙物前的運(yùn)動(dòng)與圖1中的計(jì)算是一致的（t=3）,當(dāng)顆粒到達(dá)障礙物時(shí),由于障礙物的阻礙作用,前部的顆粒爬上障礙物（t=6）,并在運(yùn)動(dòng)到一定的高度后從障礙物的兩側(cè)流下（t=9）,顆粒被障礙物分開為2部分,這2部分到達(dá)水平面后開始減速（t=12）,在水平面上逐漸形成了兩個(gè)對(duì)稱的顆粒堆（t=12～21）,并在t=24時(shí)顆粒堆積基本完成。與設(shè)置朝向上坡障礙物的結(jié)果相比,由于角度設(shè)置不同,朝向下坡的障礙物截留了部分顆粒,并且由于障礙物的截留作用水平面上的兩個(gè)顆粒堆變小,在2個(gè)顆粒堆中間形成有無顆粒堆積的“安全區(qū)域”,但也因朝向下坡的障礙物缺乏良好的分導(dǎo)作用而較小。由此可見,障礙物的不同角度設(shè)置將對(duì)顆粒的流動(dòng)過程產(chǎn)生重要的影響。

顆粒流經(jīng)朝向下坡障礙物過程中的這些特點(diǎn)已經(jīng)被相關(guān)的試驗(yàn)所證實(shí)[8]。

2.4 障礙物的位置設(shè)置對(duì)顆粒堆積的影響

在實(shí)際的工程中,障礙物往往需要根據(jù)地形地物條件以及防治的需要進(jìn)行布置。這時(shí)障礙物的布置不再像前面討論的2種典型過程那樣具有對(duì)稱性,障礙物的位置設(shè)置將會(huì)對(duì)顆粒流動(dòng)過程及最終的堆積形態(tài)產(chǎn)生影響。因此討論障礙物的位置設(shè)置對(duì)顆粒流動(dòng)和堆積過程的影響,研究其中的變化規(guī)律對(duì)于優(yōu)化防災(zāi)工程設(shè)計(jì)和山區(qū)規(guī)劃設(shè)計(jì)是很有意義的。

圖5 障礙物位置設(shè)置對(duì)顆粒堆積的影響

在相同的計(jì)算條件下將3.2節(jié)中所討論的朝向上坡的障礙物在計(jì)算平面上平移向量V,得到4種不同的障礙物布置。在不同的障礙物布置條件下進(jìn)行計(jì)算,分別得到了各自的顆粒流動(dòng)過程,他們均在t=24時(shí)基本完成了堆積,堆積等值線如圖5所示（等值線間距為0.05）。

從圖5結(jié)果中可以看到,顆粒在流經(jīng)不同位置的障礙物后,均在障礙物的下游形成了2個(gè)顆粒堆,且在2個(gè)顆粒堆的中間形成了一個(gè)沒有顆粒堆積的“安全區(qū)域”。障礙物沿坡向（x方向）的不同位置設(shè)置影響了“安全區(qū)域”的范圍和顆粒堆的形態(tài),這是由顆粒流動(dòng)中不同的能量消耗過程引起的,障礙物沿垂直坡向（y方向）的不同位置設(shè)置則顯著的影響了兩個(gè)顆粒堆的質(zhì)量分配,同時(shí)“安全區(qū)域”具有與障礙物相同的位移趨勢(shì),體現(xiàn)出障礙物對(duì)坡下的保護(hù)作用。

3 結(jié)論

雪崩、滑坡、滾石和泥石流等都是常見的山地災(zāi)害,他們的發(fā)展具有很強(qiáng)的破壞性,對(duì)人類的生命和財(cái)產(chǎn)安全構(gòu)成了嚴(yán)重的威脅。為防治和減小災(zāi)害,就必須加強(qiáng)對(duì)災(zāi)害運(yùn)動(dòng)特征和防護(hù)措施的研究。

這些山地災(zāi)害具有顆粒流動(dòng)的運(yùn)動(dòng)特點(diǎn),因而可以采用顆粒流動(dòng)的相關(guān)理論對(duì)其運(yùn)動(dòng)過程進(jìn)行研究。論文采用了模擬顆粒流動(dòng)的SH理論和Roe的有限體積離散格式對(duì)顆粒流經(jīng)障礙物的過程進(jìn)行了數(shù)值計(jì)算,模擬了在不同障礙物設(shè)置條件下的顆粒流動(dòng)過程,研究了不同的障礙物設(shè)置對(duì)顆粒流動(dòng)和堆積過程的影響。

數(shù)值模擬的結(jié)果表明在顆粒流動(dòng)路徑中設(shè)置障礙物可以改變顆粒流動(dòng)的方向,在障礙物高度足夠的條件下,會(huì)在障礙物的下坡方向形成一個(gè)沒有顆粒堆積的“安全區(qū)域”,其可以看做是為人居環(huán)境和基礎(chǔ)設(shè)施免遭災(zāi)害破壞的實(shí)際可行的應(yīng)用區(qū)域,對(duì)災(zāi)害防護(hù)的設(shè)計(jì)具有重要的意義。但是障礙物的不同設(shè)置會(huì)影響顆粒流動(dòng)和堆積的過程,為了達(dá)到防護(hù)的目的,障礙物的設(shè)置一定要合理。

在實(shí)際的工程設(shè)計(jì)中地形地物條件更加復(fù)雜,結(jié)合數(shù)值模擬的方法可以對(duì)防治障礙物的設(shè)置進(jìn)行優(yōu)化設(shè)計(jì),從而得到經(jīng)濟(jì)而合理的防治措施方案,為山區(qū)災(zāi)害防護(hù)和移民規(guī)劃提供重要的參考和依據(jù)。

[1]SAVAGE S B,HUTTER K.The motion o f a finite mass of granular material down a rough incline[J].Journal of Fluid Mechanics,1989,199：177-215.

[2]SAVAGE S B,HUTTER K.The dynam ics of avalanches of granular materials from initiation to runout.I：Analysis[J].A cta Mechanica,1991,86（1-4）：201-223.

[3]CH IOU M C,WANG Y,HUTTER K.In fluence of obstacles on rapid granu lar flow s[J].A cta Mechanica,2005,175（1-4）：105-122.

[4]IVERSON R M,DENLINGER R P.Flow of variably fluidized granular masses ac ross three-dimensional terrain：1.Coulomb m ixture theory[J].Journal of Geophysical Research,2001,106（B1）：537-552.

[5]王純祥,白世偉,江崎哲郎,等.基于 GIS泥石流二維數(shù)值模擬[J].巖土力學(xué),2007,28（7）：1359-1368.

W ANG CHUN-XIANG,BA ISH I-WEI,ESAKI T,et al.GIS based tw o dimensional numerical simulation of debris flow[J].Rock and SoilMechanics,2007,28（7）：1359-1368.

[6]李珂,唐紅梅,易麗云,等.泥石流溝岸耦合三維數(shù)值仿真[J].重慶建筑大學(xué)學(xué)報(bào),2008,30（1）：68-76.

LI KE,TANG HONG-MEI,YI LI-YUN,et al.Numerical simulation on coup ling problem betw een debris flow and bank of debris flow va lley[J].Journal of Chongqing Jianzhu University,2008,30（1）：68-76.

[7]李同春,李楊楊,章書成,等.泥石流泛濫區(qū)域數(shù)值模擬[J].水利水電科技進(jìn)展,2008,28（6）：1-4.

LI TONG-CHUN,LI YANG-YANG,ZHANG SHUCHENG,etal.Numerical simu lation on inundation area of debris flow[J].Advances in Science and Technology of Water Resources,2008,28（6）：1-4.

[8]馬宗源,張駿,廖紅建.黏性泥石流攔擋工程數(shù)值模擬[J].巖土力學(xué),2007,28（s1）：389-392.

M A ZONG-YUAN,ZHANG JUN,LIAO HONGJIAN.Numerical simulation o f viscous debris flow b lock engineering[J].Rock and SoilM echanics,2007,28（s1）：389-392.

[9]HUNGR O.Numerical modelling of the motion of rapid,flow-like landslides for hazard assessment[J].Ksce Journal o f Civil Engineering,2009,13（4）：281-287.

[10]魯曉兵,王義華,王淑云,等.碎屑流沿坡面運(yùn)動(dòng)的初步分析[J].巖土力學(xué),2004,25（S2）：598-600.

LU XIAO-BING,WANG YI-HUA,WANG SHUYUN,et al.The Primary analysis on the castic gain flow[J].Rock and Soil M echanics,2004,25（S2）：598-600.

[11]YU B,DALBEY K,WEBB A,etal.Numerical issues in com puting inundation areas over natural terrains using Savage-Hutter theory[J].Natural H azards,2009,50（2）：249-267

[12]KALAND C,STRUCKMEIER J.A kinetic scheme for the savage hutter equations[J].M athematicalMethods in the App lied Sciences,2008,31（16）：1922-1945.

[13] BOUCHUT F, FERNANDEZ NIETO E D,MANGENEY A,et al.On new erosion models of Savage-Hutter type for avalanches [J].A cta Mechanica,2008,199（1-4）：181-208.

[14]PUDASA IN IS P,HUTTER K.Rapid shear flow s of d ry granu larmasses down curved and twisted channels[J].Journalof Fluid Mechanics,2003,495：193-128.[15]WANG Y,HUTTER K,PUDASA INI S P.The Savage-Hutter theory：a system of partial differential equations for avalanche flow s o f snow,debris and mud[J].Journal of App lied M athematics and Mechanics,2004,84（8）：507-527.

[16]GRAY JM N T,TA IY C,NOELLE S.Shock w aves,dead zones and particle-free regions in rapid granu lar free-surface flow s[J].Journal o f Fluid Mechanics,2003,491：161-181.

[17]TEUFELSBAUERH,WANG Y,CH IOU M C.Flowobstacle interaction in rapid granu lar avalanches：DEM simu lation and comparison with experiment[J].Granular Matter,2009,11（4）：209-220.

[18]ROE P L.Approximate Riemann so lvers parameter vectors and difference schemes[J].Journal of Computational Physics,1981,43（2）：357-372.

[19]OSHER S.Convergence of generalized MUSCL schemes[J].Siam Journal of Numerical Analysis,1996,22（5）：947-961.

[20]HUBBARD M E.Mu ltidimensional slope limiters for MUSCL type finite volume schemes on unstructured grids[J].Journal of Com putational Physics,1999,155（1）：54-74.

[21]DARW ISH M S,MOUKALLED F.TVD schemes for unstructured grids[J].International Journal of Heat and Mass Transfer,2003,46（4）：599-611.

[22]BRUFAU P,GARCíA-NAVARRO P,VáZQUEZCENDóN M E.Zero mass error using unsteady wetting-drying conditions in shallow flows over dry irregu lar topography[J].International Journal for Numerical Methods in Fluids, 2004,45（10）：1047-1082.

（編輯王秀玲）