溝道二維泥石流運動和沖淤數值模型研究

張萬順，趙琰鑫，崔鵬，彭虹，陳雪嬌

(1.武漢大學資源與環境科學學院，430079，武漢;2.中國科學院成都山地災害與環境研究所，640041，成都;3.武漢大學水利水電學院，430072，武漢)

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

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

筆者以水沙混合流模型為基礎，采用混合流沙量動態變化模式，導出泥石流運動控制方程組，建立適于模擬溝道泥石流的二維非恒定數值模型。利用建立的數值模型，模擬云南蔣家溝流域典型泥石流動態演進過程和溝道的沖淤演變，不僅具有一定的理論意義，也可為流域泥石流風險預測和災害治理提供參考。

1 泥石流運動模型

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

1.1 基本方程

連續性方程

動量方程

泥沙對流擴散方程

河床變形方程

床沙級配調整方程

式中: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為泥沙的紊動擴散系數，m2/s;Em為混合層厚度，m;αs為混合流飽和泥沙恢復系數，m2/s;s*i 為混合流對i 組泥沙的挾沙力，kg/m3;ε1為模型參數，當泥石流沖刷到原始河床時ε1=0，否則ε1=1;Pmi為床沙i 組泥沙的比例，%;Pmi，0為原始床沙i 組泥沙的比例，%;g 為重力加速度，m/s2;α 為積分形狀系，α=1;τxb、τyb分別為泥石流運動在x、y 方向床面的切應力分量，N/m2，計算公式為:

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

式中:τB為泥石流賓漢結構應力，N/m2;γ 為泥石流賓漢體的黏滯系數，N·s/m2。

1.2 模型輔助方程

1.2.1 泥石流賓漢結構應力模型 泥石流賓漢結構應力τB是泥石流的網格結構強度和沙粒間的摩擦共同作用影響的結果。τB與泥石流漿體的濃度和流體內細顆粒和粗顆粒的含量有密切的關系。許多研究人員在實驗研究和野外原型觀測的基礎上，總結出了許多較有代表性的經驗公式(表1)。

表1 泥石流賓漢體結構應力公式Tab.1 Bingham yield stress formula of debris flow slurry

本文選用費祥俊［11］提出的泥石流漿體的結構力關系式

1.2.2 挾沙力公式 本研究以張瑞瑾等［14］公式為基礎，推導出泥石流挾沙力公式

2 泥石流運動模型在云南蔣家溝的應用

2.1 研究區域概況

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

2.2 模型輸入條件

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

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

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

2.3 泥石流運動過程的模擬結果

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

圖3 泥石流運動過程模擬結果Fig.3 Simulated debris flow pattern

2.4 泥石流沖淤的模擬結果

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

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

3 小結

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

模型以云南東川蔣家溝為例，研究了泥石流運動隨時間的演進過程，定量分析了泥石流沖淤影響下溝道形態的演變，模型模擬結果和野外實際觀測到的泥石流運動及沖淤特征基本符合，證明該模型能夠較客觀地反映泥石流龍頭隨時間的動態變化過程和泥石流影響下溝道上下游不同區域的沖淤規律，這說明模型建立的機制正確，參數取值也基本合理。模型對于泥石流災害預測和防治具有重要的現實意義。

［1］ 吳積善，田連權，康志成，等.泥石流及其綜合治理［M］.北京: 科學出版社，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］ 王光謙，倪晉仁.泥石流動力學基本方程［J］.科學通報，1994，39(18):1700-1704

［6］ 余斌.二維定常泥流的數值模擬［J］.自然災害學報，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］ 錢寧，王兆印.泥石流運動機理的初步探討［J］.地理學報，1984，39(1):33-43

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

［10］沈壽長，謝慎良.粘性泥石流的結構模式和流變特性［J］.鐵道工程學報，1986(4):26-33

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

［12］王裕宜，詹錢登，鄒仁元，等.泥石流漿體屈服應力綜合表達式的研究［J］.自然災害學報，1999，8(3):103-110

［13］王裕宜，詹錢登，韓文亮，等.粘性泥石流體的應力應變特性和流速參數的確定［J］.中國地質災害與防治學報，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