四川壤塘县茸郎沟泥石流发育特征及危险性

韩 雨 菲， 胡 卸 文， 韩 玫， 熊 冲 冲， 盛 豪

(西南交通大学，四川 成都 610031)

0 引 言

泥石流灾害常见于山区，多是由于集中降雨诱发，具有发生突然、历时短、破坏力大、难以预报等特点[1]。泥石流的危险区范围是指泥石流在流通、堆积等冲淤过程中造成毁灭性灾难的区域。合理可靠的泥石流危险区范围确定，对山区交通选线、城镇村落选址布局以及泥石流综合防治措施的选取等都具有重要的实际意义和科学价值[2]。

目前最为广泛应用的泥石流危险区范围确定主要有物理模拟试验法、数学模型法及数值模拟法。通过查阅国内外文献，表明物理模拟试验及数学模型预测大多以泥石流影响因子为基础进行分析，获得泥石流堆积范围，结果较为单一，而数值模拟法一般不仅能预测泥石流危险范围，也能够体现泥石流运动过程中的动力特征。近年来泥石流流体动力特性及数值模拟方法探索不断深入，胡凯衡[3]根据流团模型原理，假定泥石流摩阻为速度的二次多项式，对泥石流在堆积区的堆积过程和特征进行了数值模拟；韦方强[4]等结合了泥石流流体数值模拟技术和GIS技术，对泥石流的堆积运动过程进行数值模拟，获取了模拟中泥石流堆积范围中每一点的最大动量，以动量作为分区指标建立了模型。马宗源、李珂等人[5-6]应用ANSYS 软件，对粘性泥石流流体进行数值模拟分析，得到了泥石流的流场数据，并初步讨论了泥石流流体对拦挡坝等防治设施的影响。四川阿坝州壤塘县茸郎泥石流沟为则曲河右岸一级支流，茸木达乡于茸郎沟沟口，有乡政府、邮政所、供电站、茸木达中心校、卫生院等基础设施，居民约300人，一旦泥石流爆发则会对当地居民生命财产安全构成严重威胁(见图1)。

图1 茸郎沟沟口乡政府所在地

本文在对茸郎沟泥石流的形成条件及发育特征进行实地野外调查的基础上，运用CFX软件对其运动特征参数及运动过程开展研究，通过对茸郎沟20 a与50 a一遇降雨频率下的泥石流运动过程进行数值模拟，以预测其危险区范围，进而评价该泥石流沟对茸木达乡居民点的危害程度。

1 茸郎沟泥石流发育特征

茸郎沟沟域面积为45.59 km2，主沟纵长14.81 km，整体上可分为清水动力区、形成流通区和堆积区，平均纵坡降51‰，上游清水汇流区及形成流通区呈"V"型沟谷地貌，沟谷较为狭窄，纵坡较陡，水流湍急，且动态变化较大；下游沟段呈“U”型堆积谷地貌，谷宽一般20～50 m左右。茸郎沟支沟较为发育，中游及下游发育9条支沟，流通区两侧岸坡坡度35°～45°，(见图2，图3)。

图2 茸郎沟流域平面图

茸郎沟所处区域高原季风性气候十分显著，干湿季节明显，气温和降雨垂直分带明显。壤塘县多年平均降水量705.6 mm，从最近30 a的降雨资料分析，近年来降雨有逐渐增大的趋势。茸郎沟流域出露地层主要为第四系松散堆积物残坡积层(Q4)及三叠系上统新都桥组(T3xd)板岩-砂岩地层。茸木达地区属构造剥蚀高原地貌，由于局部沟谷段坡面及沟道松散堆积物丰富，在丰水期易形成泥石流灾害。

图3 茸郎沟主沟纵剖面图

经调查走访得知，茸郎沟最近于2015年连续强降雨后发生过一定规模的泥石流灾害，泥石流持续10多分钟，泥石流冲出物质为稀状泥浆并夹杂少量大块石，该次泥石流淹过耕地，漫入当地居民家中并造成少量损毁，直接经济损失达几十万元，所幸未造成人员伤亡。

2 泥石流运动特征参数

泥石流运动特征参数是数值模拟求解设置、分析对比的依据，所需确定参数主要为泥石流重度及峰值流量，可采用雨洪修正法计算获得。

2.1 泥石流重度

现场配浆法与查表法所得泥石流重度基本吻合，重度取值15.30 kN/m3，属于稀性泥石流。

2.2 泥石流流速

茸郎沟泥石流属稀性泥石流，流速计算公式采用西南地区(铁二院陈光曦)公式[7]：

(1)

沟口附近泥石流流速计算结果见表1。

表1 茸郎沟泥石流沟口流速计算结果

(3)泥石流峰值流量

泥石流峰值流量采用雨洪修正法进行计算：

QC=DC(1+φ)QP

(2)

QP=0.278ψAS/tn

(3)

式中QC泥石流断面峰值流量(m3/s)；φ泥沙修正系数；QP暴雨洪峰流量；Dc堵塞系数；ψ洪峰径流系数；S单位历时的暴雨平均强度(mm/h)；A流域面积(km2)；t汇流时间(h)；n暴雨递减指数。

沟口泥石流峰值流量计算结果见表2。

表2 茸郎沟泥石流沟口流量计算结果

3 基于CFX的泥石流运动过程数值模拟

3.1 CFX软件简介

CFX是一种计算流体动力学(Computional Fluid Dynamics)，即CFD软件，通过求解流体控制方程，模拟流体流动时各种相关物理现象，其全隐式耦合多网格线性求解器支持多相流等复杂物理模型，并行计算能力强大，收敛速度快且结果精准。

CFX计算基础为区域离散化数值计算，常用方法有有限差分法、有限元法及有限体积法。CFX求解器采用的为有限体积法，相比于其他两种方法，能够精确满足物理量积分守恒，是目前在流体流动求解中最有效的数值计算方法。

CFX求解过程包括建模、前处理、求解、后处理几个阶段。

3.2 泥石流流场几何设定

泥石流数值模拟的基础是建立泥石流流场几何模型，综合应用Auto CAD，Surfer，Design Modeler等建模工具，通过点-线-面-体的构建方式建立茸郎沟沟口堆积区三维地质模型，模型平面范围为500 m×400 m，流体域高度20 m。在WorkbenchMesh模块中调用ICEM-CFD前处理器进行网格划分，选择生成简单，适合复杂几何图形的四面体(Tetrahedrons)[8]，全模型共生成59 985个节点，293 138个单元(见图4)。

图4 茸郎沟泥石流沟口堆积区模型及网格划分

3.3 前处理

前处理模块中，主要需设定泥石流流体流变模型及边界条件。

3.3.1 泥石流流变模型

流变模型的正确选择对于数值模拟结果的准确性及泥石流流动规律的探索有着重要意义。大量试验结果表明，一般的泥石流流体性质更适合Bingham模型[9]，其数学模型表达如下：

(4)

3.3.2 边界条件设定

CFX中的边界有入口(Inlets)、出口(Outlets)、壁面(Walls)、开放式边界(Opening)、对称面(Symmetry)5种类型。对于茸郎沟沟口模型，设定如下：

(1)入口：堆积区上游，沟道泥石流流入处。在质量和动量(MassandMomentum)模块中设定流体进入方式为质量流量(BulkMassFlowRate)，定义多项式函数Q(t)，表示泥石流流量随时间的变化。设定t=0s时取洪水流量，总模拟时间1/3处为泥石流峰值流量(见图5)。函数通过CFXExpressionLanguage(CEL语言)编译并定义计算中的泥石流流量变化。

(a)P=5%

(b)P=2%图5 茸郎沟泥石流堆积区流量变化过程曲线

(2)出口：模型四侧边界。主要设定出口条件为静态压强，流体性质为亚声速流体(Subsonic)。

(3)开放式边界：模型上部。设置为开放压强与方向，压力1.01×105Pa。

(4)壁面：泥石流流动面。设置为无滑移壁面，粗糙度取0.1m。

设置完成模型(见图6)。

3.4 求解

CFX求解管理器可监视整个求解过程，主要体现在入口与出口质量流量变化，两次模拟进出口质量流量曲线求解结果(见图7)。

图6 模型边界条件设置

(a)P=5%

(b)P=2%图7 质量流量曲线求解结果

由图7可见，进口质量流量曲线与设置的泥石流流量过程曲线吻合。20a一遇泥石流模拟中，出口质量流量曲线分别在t=75s处陡然下降，说明此时泥石流到达了出口边界并流出，在t=320s附近，泥石流峰值流量从入口流入，t=415s时到达出口边界，出口峰值流量较进口峰值流量减小到20.92m3/s，说明建筑体对于泥石流有拦挡与淤积作用。50a一遇泥石流模拟中，泥石流在t=65s处从出口边界流出，较20a一遇泥石流流速较大，与实际情况相符，泥石流进口峰值流量出现在t=450s附近，75s后流出边界，出口峰值流量为34.64m3/s。

3.5 后处理及分析

3.5.1 泥石流运动过程分析

CFX后处理功能多样，可以多种形式直观展示数值模拟结果，如等值线图、矢量图、动画输出等。本文以20a一遇泥石流模拟为例，通过泥石流流速云图描述分析运动过程(见图8)。

(a)t=50 s

(b)t=90 s

(c)t=340 s

(d)t=550 s

(e)t=900 s图8 不同时刻泥石流液面速度云图(a～e)

t=50s时，沿沟道运动的泥石流推进300m，首先淹没了离入口最近的房屋，而后产生分支，大部分仍沿主沟道流动，小部分漫过沟道进入平坦宽阔的堆积区并逐步扩散，由于沟道狭窄，沿沟道运动速度快，最大速度达到8.679m/s，堆积区漫流速度分布在2～4m/s之间。t=90s时，泥石流已经流过出口，在居民区漫流范围逐步扩大，乡政府大院已经完全淹过，下游低矮民房开始被淤埋，此时危险范围沿沟道长度400m，垂直沟道60m。t=340s时，峰值流量已经流入，泥石流整体运动速度增大，平均速度达5.5m/s，最大速度9.130m/s，由于建筑物的拦挡停淤作用，泥石流在居住区内的平均运动速度较小，为1.5m/s。t=550s时，泥石流在堆积区漫流范围达到最大，沿沟道方向长度达430m，垂直沟道最大达280m，此时速度随着流量降低而减小，t=900s时，泛滥范围稍有缩小，入口流量回到洪水流量值，泥石流整体流速达最小值。

3.5.2 泥石流危险区范围分析

取模拟结果中最大泥石流漫流范围作为危险区范围，20a一遇泥石流取t=550s时刻，50a一遇泥石流取t=680s时刻，以及现场访问得到2015年实际受灾范围(见图9)。

图9 泥石流危险区范围

由图9可知，2015年泥石流爆发频率小于20a一遇，受沟道与山体的限制，20a与50a一遇泥石流在垂直沟道方向泛滥范围差别不大，差异主要在于顺沟道方向的最大堆积长度，20a一遇泥石流最大堆积长度为430m，50a一遇泥石流最大堆积长度为470m，泥石流最大泛滥范围分别占居民区总面积的67%、79%。一旦暴发泥石流，茸木达乡政府所辖居民区将面临重大威胁，需采取诸如堆积区上游侧单边或双边防护堤等有效防治措施。

4 结 论

(1)位于壤塘县的茸郎沟沟域面积为45.59km2，主沟纵长14.81km，平均纵坡降51‰，为则曲河右岸一级支流，茸木达乡政府位于茸郎沟沟口，有乡政府、邮政所、供电站、茸木达中心校、卫生院等基础设施，居民约300人，2015年连续强降雨后发生过一定规模的泥石流灾害，对沟口居住区形成部分淤埋，潜在危害性大。

(2)数值模拟结果显示，茸郎沟泥石流20a一遇泥石流通过堆积区平均速度为4.30m/s，最大流速为9.13m/s，50a一遇泥石流通过堆积区平均速度为4.78m/s，最大流速为10.00m/s。由于堆积区平坦开阔，建筑物较多，阻挡淤积作用使得泥石流在居民区内漫流速度较沟道内小得多。

(3)数值模拟所得到的危险区范围显示，20a与50a一遇泥石流均对沟口茸木达乡有较大威胁，最大堆积长度分别为430m、470m，居民区受灾面积分别达67%、79%。说明一旦暴发泥石流，对茸木达乡居民构成严重威胁，需采取有效防治措施。