張萬(wàn)順，趙琰鑫，崔鵬，彭虹，陳雪嬌

(1.武漢大學(xué)資源與環(huán)境科學(xué)學(xué)院，430079，武漢;2.中國(guó)科學(xué)院成都山地災(zāi)害與環(huán)境研究所，640041，成都;3.武漢大學(xué)水利水電學(xué)院，430072，武漢)

泥石流是山區(qū)特有的一種自然地質(zhì)現(xiàn)象，它是在降水、地形地貌、地質(zhì)構(gòu)造、固體堆積物、植被覆蓋度及人類(lèi)活動(dòng)等多種因素共同作用下，發(fā)生在溝谷或山坡上的一種挾帶大量泥砂、石塊和巨礫等固體物質(zhì)的特殊洪流，具有高密度、高流速、高流量、短歷時(shí)和寬級(jí)配等特征［1］。泥石流高速運(yùn)動(dòng)方式往往給泥石流影響范圍內(nèi)的房屋、橋梁等造成毀滅性的破壞，同時(shí)伴隨著泥石流運(yùn)動(dòng)過(guò)程大量顆粒物質(zhì)的搬運(yùn)和堆積，可以在短時(shí)間內(nèi)對(duì)溝道形態(tài)演變產(chǎn)生巨大的影響，主要表現(xiàn)在在溝道上游強(qiáng)烈下切，導(dǎo)致滑坡和崩塌活動(dòng)，而溝道下游由于不同程度的淤積，造成大片農(nóng)田、鐵路、公路嚴(yán)重受災(zāi)。泥石流的運(yùn)動(dòng)過(guò)程和溝道內(nèi)沖淤過(guò)程的數(shù)值模擬研究對(duì)泥石流災(zāi)害防治具有非常重要的作用。

目前研究中采用的泥石流動(dòng)力學(xué)模型可分為2大類(lèi):一類(lèi)是具有統(tǒng)一內(nèi)在屈服應(yīng)力和黏性的單一相或偽一相流體模型，包括賓漢體模型［2］、膨脹體模型［3］及混合流體模型［4］;另一類(lèi)是兩相流體模型［5］，即把泥石流體中的漿體(由固相中的細(xì)顆粒和液相水組成)和粗顆粒各視為一種流體。近年來(lái)，泥石流數(shù)值模型研究工作已經(jīng)取得了一定進(jìn)展。王光謙等［5］將泥石流漿體用賓漢體模型、粗顆粒用膨脹體模型分別進(jìn)行描述，建立泥石流的兩相流二維模型。余斌［6］利用賓漢體模型，應(yīng)用固液兩相的混合流方程對(duì)泥石流進(jìn)行數(shù)值模擬，計(jì)算了河床底部有凸臺(tái)情況下的泥流分布。張萬(wàn)順等［7］以水沙混合流為基礎(chǔ)建立了適合泥石流模擬的二維非恒定流模型，研究泥石流運(yùn)動(dòng)和主河交匯的相互作用。以上研究對(duì)于溝道泥石流運(yùn)動(dòng)模擬進(jìn)行了積極的探索，但是仍然存在不足，對(duì)陣性非恒定泥石流運(yùn)動(dòng)演進(jìn)過(guò)程研究較少，特別是缺乏對(duì)溝道的沖淤及其與泥石流運(yùn)動(dòng)過(guò)程之間的關(guān)系的描述，不能夠系統(tǒng)地、定量地模擬研究泥石流運(yùn)動(dòng)影響下溝道的沖淤動(dòng)態(tài)過(guò)程。

筆者以水沙混合流模型為基礎(chǔ)，采用混合流沙量動(dòng)態(tài)變化模式，導(dǎo)出泥石流運(yùn)動(dòng)控制方程組，建立適于模擬溝道泥石流的二維非恒定數(shù)值模型。利用建立的數(shù)值模型，模擬云南蔣家溝流域典型泥石流動(dòng)態(tài)演進(jìn)過(guò)程和溝道的沖淤演變，不僅具有一定的理論意義，也可為流域泥石流風(fēng)險(xiǎn)預(yù)測(cè)和災(zāi)害治理提供參考。

1 泥石流運(yùn)動(dòng)模型

從宏觀上，泥石流是水沙的混合流體，在演變過(guò)程中由于水沙的含量不同，表現(xiàn)出塑性蠕動(dòng)流、黏性陣流、陣性連續(xù)流和稀性連續(xù)流等不同運(yùn)動(dòng)形態(tài);從微觀上，泥石流屬于固液兩相流。但若直接采用兩相流模型，則泥石流固液兩相間的相互應(yīng)力作用力難以具體地加以量化，實(shí)際應(yīng)用中還存在很大的困難。筆者采用的泥石流運(yùn)動(dòng)模型以水沙混合流模型為基礎(chǔ)，可以避免考慮兩相間的微觀作用力，同時(shí)也可以利用已有的泥石流的應(yīng)力或阻力的研究成果。

1.1 基本方程

連續(xù)性方程

動(dòng)量方程

泥沙對(duì)流擴(kuò)散方程

河床變形方程

床沙級(jí)配調(diào)整方程

式中:u、v 為x、y 方向的速度分量，m/s;h 為泥石流深，m;zb為溝床高程，m;η 為泥石流自由表面高程，m;ci為第i 組泥沙的含沙量，kg/m3;ωi為第i 組泥沙的沉降速度，m/s;ρs為泥沙密度，kg/m3; ρ 為泥石流混合密度，kg/m3為水密度，kg/m3;εs為泥沙的紊動(dòng)擴(kuò)散系數(shù)，m2/s;Em為混合層厚度，m;αs為混合流飽和泥沙恢復(fù)系數(shù)，m2/s;s*i 為混合流對(duì)i 組泥沙的挾沙力，kg/m3;ε1為模型參數(shù)，當(dāng)泥石流沖刷到原始河床時(shí)ε1=0，否則ε1=1;Pmi為床沙i 組泥沙的比例，%;Pmi，0為原始床沙i 組泥沙的比例，%;g 為重力加速度，m/s2;α 為積分形狀系，α=1;τxb、τyb分別為泥石流運(yùn)動(dòng)在x、y 方向床面的切應(yīng)力分量，N/m2，計(jì)算公式為:

式中:nb為河床糙率; τxx、τxy和τyy為泥石流切應(yīng)力分量，N/m2，計(jì)算公式分別為:

式中:τB為泥石流賓漢結(jié)構(gòu)應(yīng)力，N/m2;γ 為泥石流賓漢體的黏滯系數(shù)，N·s/m2。

1.2 模型輔助方程

1.2.1 泥石流賓漢結(jié)構(gòu)應(yīng)力模型 泥石流賓漢結(jié)構(gòu)應(yīng)力τB是泥石流的網(wǎng)格結(jié)構(gòu)強(qiáng)度和沙粒間的摩擦共同作用影響的結(jié)果。τB與泥石流漿體的濃度和流體內(nèi)細(xì)顆粒和粗顆粒的含量有密切的關(guān)系。許多研究人員在實(shí)驗(yàn)研究和野外原型觀測(cè)的基礎(chǔ)上，總結(jié)出了許多較有代表性的經(jīng)驗(yàn)公式(表1)。

表1 泥石流賓漢體結(jié)構(gòu)應(yīng)力公式Tab.1 Bingham yield stress formula of debris flow slurry

本文選用費(fèi)祥?。?1］提出的泥石流漿體的結(jié)構(gòu)力關(guān)系式

1.2.2 挾沙力公式 本研究以張瑞瑾等［14］公式為基礎(chǔ)，推導(dǎo)出泥石流挾沙力公式

2 泥石流運(yùn)動(dòng)模型在云南蔣家溝的應(yīng)用

2.1 研究區(qū)域概況

蔣家溝流域地處云南東川金沙江支流小江河谷右岸，是小江流域內(nèi)泥石流活動(dòng)頻率最高的一條泥石流溝，流域面積48.6 km2，主溝長(zhǎng)13.9 km。溝內(nèi)巖性軟弱、巖層松散、斷裂縱橫、地形陡峻、植被稀疏，地表形態(tài)破碎，滑坡和崩塌活動(dòng)非常強(qiáng)烈，泥石流發(fā)生頻繁。蔣家溝泥石流觀測(cè)資料較多，系列較長(zhǎng)，為從事相關(guān)理論研究提供了良好的基礎(chǔ)［1］。

2.2 模型輸入條件

模型計(jì)算范圍選取蔣家溝中下游區(qū)域作為數(shù)學(xué)模型計(jì)算溝段，包括多照溝、門(mén)前溝溝口下游至入小江口，全長(zhǎng)約5 km。根據(jù)蔣家溝DEM(10 m×10 m)進(jìn)行網(wǎng)格劃分，整個(gè)計(jì)算范圍劃分成56×35=1 960個(gè)網(wǎng)格。采用中國(guó)科學(xué)院東川泥石流觀測(cè)站1975年5 月的1 次實(shí)測(cè)黏性泥石流運(yùn)動(dòng)過(guò)程作為上游輸入的流量過(guò)程，模擬泥石流在溝道內(nèi)的運(yùn)動(dòng)和沖淤過(guò)程。泥石流過(guò)程的含沙量、流量過(guò)程線如圖1 所示。在蔣家溝泥石流主要補(bǔ)給區(qū)采集土樣并分析顆粒物質(zhì)組成，作為輸入的泥石流泥沙顆粒級(jí)配，泥石流泥沙顆粒級(jí)配曲線如圖2 所示。

圖1 泥石流流量、含沙量過(guò)程線Fig.1 Changing process of discharge and sediment concentration of debris flow

圖2 泥石流泥沙顆粒級(jí)配曲線Fig.2 Grain size distribution of the debris flow sample

2.3 泥石流運(yùn)動(dòng)過(guò)程的模擬結(jié)果

圖3 是模擬的泥石流運(yùn)動(dòng)隨時(shí)間的演進(jìn)過(guò)程。由圖可見(jiàn)，在泥石流開(kāi)始運(yùn)動(dòng)時(shí)(t=0.1 h)速度較小，運(yùn)動(dòng)較緩慢，隨后由于河床比降增大，流速增大，同時(shí)泥石流沿著溝道迅速展寬(t=0.5 h)，當(dāng)泥石流流量達(dá)到峰值時(shí)(t=1.0 h)，泥石流在溝道中的流速也達(dá)到最大，并快速運(yùn)動(dòng)流出蔣家溝谷口，隨后隨著流量下降，泥石流流速也隨之降低。模擬的運(yùn)動(dòng)過(guò)程和野外觀測(cè)的泥石流實(shí)際運(yùn)動(dòng)規(guī)律基本吻合。

圖3 泥石流運(yùn)動(dòng)過(guò)程模擬結(jié)果Fig.3 Simulated debris flow pattern

2.4 泥石流沖淤的模擬結(jié)果

圖4 是模擬的泥石流在泥沙含量最大時(shí)(t=2.0 h)和整個(gè)過(guò)程結(jié)束時(shí)(t=3.6)溝道的沖淤狀況。根據(jù)模擬結(jié)果，淤積和沖刷幅度與泥石流的流量、含沙量有密切關(guān)系。泥石流主要在溝口淤積，而且隨著泥石流泥沙的沉積，溝口的淤積厚度不斷增加。在泥石流運(yùn)動(dòng)過(guò)程中，由于溝床比降增大，泥石流流速增加，挾沙力增加，同時(shí)由于泥沙在溝道上游沉積，泥石流含沙量減小，因此下游溝道受到明顯沖刷。模擬的泥石流淤積厚度一般為0 ～0.5 m，沖刷變幅一般為0 ～0.02 m，說(shuō)明泥石流影響下蔣家溝下游溝道的基本上達(dá)到了沖淤平衡的狀況，泥石流過(guò)程前后地形改變不大。沖淤模擬結(jié)果符合蔣家溝溝床泥石流沖淤演變的基本規(guī)律。

圖4 泥石流沖淤模擬結(jié)果Fig.4 Simulated bed erosion-deposition depth

3 小結(jié)

結(jié)合水沙混合流模型和賓漢體模型理論建立了適于描述泥石流的二維非恒定兩相流數(shù)值模型，避免考慮顆粒間的微觀作用力，并從整體上保證模型具有較高的精度。

模型以云南東川蔣家溝為例，研究了泥石流運(yùn)動(dòng)隨時(shí)間的演進(jìn)過(guò)程，定量分析了泥石流沖淤影響下溝道形態(tài)的演變，模型模擬結(jié)果和野外實(shí)際觀測(cè)到的泥石流運(yùn)動(dòng)及沖淤特征基本符合，證明該模型能夠較客觀地反映泥石流龍頭隨時(shí)間的動(dòng)態(tài)變化過(guò)程和泥石流影響下溝道上下游不同區(qū)域的沖淤規(guī)律，這說(shuō)明模型建立的機(jī)制正確，參數(shù)取值也基本合理。模型對(duì)于泥石流災(zāi)害預(yù)測(cè)和防治具有重要的現(xiàn)實(shí)意義。

［1］ 吳積善，田連權(quán)，康志成，等.泥石流及其綜合治理［M］.北京: 科學(xué)出版社，1993:17-19

［2］ Johnson A M，Rahn P H.Mobilization of debris flows［J］.Zeitschrift fur Geomorphologie，1970，9(Sup):168-186

［3］ Takahashi T.Mechanical characteristics of debrisflow［J］.Journal of the Hydraulics Division，1978(104):1153-1169

［4］ Chen Chenglung.Generalized viscoplastic modeling of debris flow［J］.Journal of Hydraulic Engineering，1988，114(3):237-258

［5］ 王光謙，倪晉仁.泥石流動(dòng)力學(xué)基本方程［J］.科學(xué)通報(bào)，1994，39(18):1700-1704

［6］ 余斌.二維定常泥流的數(shù)值模擬［J］.自然災(zāi)害學(xué)報(bào)，1995，4(4):96-99

［7］ Zhang Wanshun，Cui Peng，Qiao Fei.Numerical simulation of interaction between tributary debris flow and main river［J］.Journal of Sichuan University:Engineering Science Edition，2005(S1):75-80

［8］ 錢(qián)寧，王兆印.泥石流運(yùn)動(dòng)機(jī)理的初步探討［J］.地理學(xué)報(bào)，1984，39(1):33-43

［9］ 費(fèi)祥俊，朱平一.泥石流的粘性及其確定方法［J］.鐵道工程學(xué)報(bào)，1986(2):9-16

［10］沈壽長(zhǎng)，謝慎良.粘性泥石流的結(jié)構(gòu)模式和流變特性［J］.鐵道工程學(xué)報(bào)，1986(4):26-33

［11］費(fèi)祥俊.黃河中下游含沙水流粘度的計(jì)算模型［J］.泥沙研究，1991(2):1-12

［12］王裕宜，詹錢(qián)登，鄒仁元，等.泥石流漿體屈服應(yīng)力綜合表達(dá)式的研究［J］.自然災(zāi)害學(xué)報(bào)，1999，8(3):103-110

［13］王裕宜，詹錢(qián)登，韓文亮，等.粘性泥石流體的應(yīng)力應(yīng)變特性和流速參數(shù)的確定［J］.中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào)，2003，14(1):9-13

［14］Zhang Ruijin，Xie Jianheng.Sedimentation research in China: Systematic Selections［M］.Beijing: China Water and Power Press，1993:57-60