四川绵竹走马岭特大泥石流动力特性的数值分析

阴 可，苟印祥，2，李承泽，吴汉辉，3

(1.重庆大学土木工程学院，重庆 400045;2.四川省交通运输厅交通勘察设计研究院，四川成都 610017;3.重庆蜀通岩土工程有限公司，重庆 401147)

0 引言

我国是一个地质灾害多发国家，其中泥石流这种地质灾害对人民生命和财产造成的损失是十分惨痛的，关于泥石流的防治和研究也是长期的、持续的。目前，泥石流的防治工程研究主要有:泥石流的评估及预测预报［1-3］，泥石流流体模型研究［4-5］，泥石流的数值模拟分析［6-10］等，其中泥石流的数值模拟分析随着计算机性能的大幅提高以及泥石流流体模型研究的进一步发展而逐步开始运用，三维泥石流的数值模拟不仅能直观反映泥石流流动的整个过程，而且能计算出流动过程中的泥石流流场变化情况，并且可以与泥石流拦砂坝等防治工程设施进行流固耦合分析，为泥石流相关研究提供一些参考。

四川省绵竹市走马岭特大泥石流位于距绵竹市区约35km的清平乡银杏沟旅游区内盐井村走马岭，为绵远河左岸支流。该沟整个流域面积为5.7km2，最高点海拔高程2044m，最低点位于走马岭沟汇入绵远河口，高程904m，相对高差1140m。主沟长度为3.5km，从沟口至沟源沟床坡降为85‰～361‰，沟源处坡降局部达620‰，平均沟床坡降为145‰，沟坡坡度一般为30°～45°。流域基岩地层依次为:寒武系下统清平组长石云母石英砂岩及页岩、震旦系上统邱家河组白云岩和纯灰岩、砂页岩夹泥质灰岩、石英砂岩等。沟内覆盖层为:坡积—洪积层卵石土，厚约4～15m，泥石流堆积层，沟槽内残留的堆积层厚度为2～3m，沟口的堆积层厚度为3～12m。走马岭流域属四川盆地西北部的龙门山推覆构造带前缘，处于2008年“5.12”地震地震带内。

2009年9月24日，清平乡普降暴雨，引发“9.24”泥石流灾害，2010年8月13日再降暴雨，最大每小时降雨量为70mm(略高于50a一遇)，日降雨量为227.5mm，走马岭再次爆发泥石流，一次性冲出量达8.3×105m3。泥石流将沟谷堵塞填高，堆积体厚度达8～15m，局部甚至达20m，河床平均淤积抬高7～10m。灾害造成多幢房屋被掩埋，公路被冲毁，绵远河改道，直接经济损失3000×104元以上。

走马岭特大泥石流治理工程拟采取“拦蓄+固源+排导+停淤+监测预警+搬迁避让”为主的工程措施。考虑在走马岭主沟上共修建7座拦砂坝(有效坝高8～12m，上游面坡坡度1∶0.70，下游面坡坡度1∶0.20，坝顶轴线长度 50～71m，1#～4#坝为天然地基，5#～7#坝采用桩基础，坝身均采用C35块石混凝土，密度为2400kg/m3)及其他防治措施，包括谷坊坝、两岸防冲墙、泥石流停淤场和泥石流排导槽等。为了给该方案及设计提供一定的参考和优化建议，特进行泥石流动力特性的数值模拟计算和分析。

1 泥石流动力特性的三维数值分析方法与思路

1.1 泥石流数值模拟的理论基础

泥石流三维数值分析的理论基础是计算流体动力学(Computational Fluid Dynamics，简称 CFD)，它可以看作是在流体三大基本方程(质量守恒方程、动量守恒方程、能量守恒方程)的控制下的流体数值模拟，通过模拟计算得到复杂流体问题的流场中各点处的速度、压力、浓度等基本物理量的分布，以及这些物理量随时间变化的规律等。本文采用基于有限体积法(又称控制体积法)的软件CFX进行流场计算和基于有限单元法的ANSYS软件进行流固耦合计算。这里有限体积法的基本思路是:将计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体积，将待解的控制方程对每一个控制体积积分，从而得出一组离散方程，其关键是在导出离散方程过程中，需要对界面上的被求函数本身及其导数的分布作某种形式的假定［11］。

1.2 泥石流浆体的流变模型

泥石流流体是非牛顿流体，近二十年来，泥石流流体模型研究有了较大的发展，常见的有:宾汉体模型、膨胀体模型、粘塑体模型、颗粒流模型、两相流模型［4］等几种。根据走马岭泥石流的特征和数值分析软件的需要采用宾汉体模型，即简化认为泥石流流体的流变特性符合以下应力应变关系:

1.3 泥石流数值模拟的流程

(1)建立泥石流流域几何模型;(2)确定初始条件及边界条件，建立计算域;(3)划分网格，确定求解控制参数;(4)有限体积法求解泥石流流场，分析结果;(5)将计算获得的流场压力等数据导入ANSYS软件，加载到拦砂坝上;(6)有限元法求解拦砂坝的应力、应变及安全系数等;(7)整理计算结果，总结分析成果，提出相关意见和建议等。

2 泥石流数值分析结果

2.1 泥石流流场数值分析

走马岭泥石流数值分析采用的泥石流参数为:重度γ=1800kg/m3，进口泥石流流量依次为，主沟Q0=74.62 m3/s，①号支沟 Q1=725.83 m3/s，②号支沟Q2=40.90 m3/s，④号支沟 Q3=56.71m3/s，⑤号支沟 Q5=124.09 m3/s，⑥号支沟 Q6=64.35 m3/s，⑦号支沟Q7=57.00 m3/s。模拟湍流类型为k-ε模型，泥石流流体的流变特性参数采用CFX软件自带的非牛顿流体参数设置对应的选项。泥石流模拟计算的典型流动状态见图1，对河床压力的等值云图见图2，流域的流速矢量图见图3。

模拟计算获得的走马岭泥石流主沟断面流场数据见表1。

图1 泥石流239.99s和1634.15s时刻的流动状态模拟图Fig.1 Debris flow state simulation at 239.99s and 1634.15s

图2 泥石流对河床压力云图(1634.15s)Fig.2 Debris flow pressure cloud on river bed(1634.15s)

图3 泥石流流域流速矢量图(1634.15s)Fig.3 Debris flow velocity vector diagram(1634.15s)

表1 走马岭泥石流主沟断面流场数据表Table 1 Debris flow Cross-section field data at the master valley of Zoumaling

模拟计算获得的走马岭泥石流5#、6#、7#支沟断面流场数据见表2(各个支沟的横断面从各支沟上游到该支沟与主沟汇流处依次按200～300m间距布置)。

从以上走马岭泥石流数值计算结果分析表明:泥石流的流速从上游到下游呈逐渐增大的趋势，同时还受河床坡降和河面宽度的影响，河床坡降越大或河床宽度越小的泥石流流速越大。泥石流的泥深随流速的增加而逐渐减小，同时也受河床坡降和河床宽度影响，河床坡降越大或河床越宽泥深则越小。

表2 走马岭泥石流5#、6#、7#支沟断面流场数据表Table 2 The slaves valleys of Zoumaling debris cross-section flow field data

2.2 泥石流流场与拦砂坝的流固耦合计算分析

采用CFX软件模拟计算了走马岭泥石流的流场数据，将获得的泥石流流场数据导入ANSYS软件中进行流固耦合计算，即可得出按设计方案考虑设置的泥石流拦砂坝上作用的流体压力和流体的运动特征，从而计算得到拦砂坝随泥石流流动产生的应力、应变等结果。

选取1#泥石流拦砂坝作为代表，满库时泥石流流体压力场分布见图4。

图4 泥石流拦砂坝满库时流体压力形态图(1#拦砂坝)Fig.4 Debris’s pressure and shape when the debris flow through the dam(1#dam)

可获得各拦砂坝最大受力状态见图5。

图5 1#～7#拦砂坝最大受力状态统计图Fig.5 The maximum stress state of 1#～7#dams

泥石流流场数值模拟显示，泥石流在拦砂坝前的泥深随着时间的增加逐渐加大，满库后漫过拦砂坝向下游流出，将泥石流瞬态模拟的各个时间点对拦砂坝的压力和泥石流泥深数据导入ANSYS软件有限元结构计算模块，安全系数计算中，拦砂坝素混凝土采用材料的极限强度，通过计算得到拦砂坝安全系数与泥深关系见图6。

图6 1#～7#拦砂坝安全系数与泥深关系图Fig.6 Safety factors and debris deep data’s relationship of 1#～7#dams

拦砂坝中1#～4#拦砂坝是采用的天然地基扩展基础，基础前部有榫齿状构造，5#～7#拦砂坝主要采用的是桩基础，基础底部平整规则。从图4～图6可知，1#～4#拦砂坝受力状态以及安全系数随泥深、坝高的变化规律具有一致性。5#～7#拦砂坝受力状态以及安全系数随泥深、坝高的变化规律也具有一致性，1#～4#拦砂坝的安全系数总体上小于5#～7#拦砂坝的安全系数。

3 结论与建议

目前，还没有关于泥石流动力特性数值计算完备的理论，以及专门针对泥石流数值分析的成熟软件，因此泥石流的数值分析结果尚不能完全准确地定量反映泥石流流场特性，但是利用CFD软件计算的结果可以反映泥石流流动的一般规律，对泥石流的学术研究和防治工作具有一定的指导意义和参考价值。

本文走马岭泥石流动力特性的数值计算和分析结果表明:

(1)泥石流对河床的最大壁面剪切应力为8.32×103Pa，位于⑤号支沟中上游;5#～7#拦砂坝下游300m段泥石流对河床壁面剪切应力也较大，达4.1×103～8.0×103Pa。

(2)泥石流的最大流速是13.57m/s，位于⑤号支沟中游;在5#～7#拦砂坝下游300m段，泥石流流速也相对较大，达到5～8m/s。

(3)泥石流对河床的最大压力为1.49×105Pa，位于②号支沟与主沟交汇处;其次在⑥号支沟上游，5#～7#拦砂坝下游300m段泥石流对河床的压力也较大，达2.2×104～1.0×105Pa。

(4)随着拦砂坝坝体高度的增加，在泥石流的冲击作用下，拦砂坝的安全系数减小;当拦砂坝过高时，坝体底部局部区域混凝土的应力在自重作用下就接近其材料的极限强度，坝体安全性较低。

(5)采用桩基础的5#～7#拦砂坝安全系数高于采用天然地基或扩展基础的1#～4#拦砂坝;随着泥石流泥深的增加，坝体的安全系数逐步降低。当泥石流满库过流时，各拦砂坝的安全系数最低。

根据上述计算分析结果，建议走马岭泥石流治理工程应注意以下几点:

(1)2#、4#、5#、6#支沟两岸的河流岸坡应加强防冲刷及固岸措施。2#～7#支沟之间应当加强河床的清理疏通工作，这有利于形成相对稳定的流动及降低泥石流下蚀与侧蚀作用。

(2)泥石流出口段有较大面积的堆积区，该区内泥石流的速度逐渐降低，且泥石流流动呈面状铺开，可以开发出来做泥石流的停淤场，但需要将场地进行整理，加大停淤区的地形横向坡度，以利于泥石流物质进一步横向流动以及停淤。

(3)5#～7#拦砂坝坝体结构极限强度安全系数在4.5(坝体素混凝土采用极限强度)以上，也就是说，此时坝体部分的材料强度满足安全性要求，并且有少量的富余。1#～4#拦砂坝坝体结构极限强度安全系数为2～3(坝体素混凝土采用极限强度)，受力位置在坝底转角处，因此建议适当加强坝体强度，或采取降低坝高，分多级建坝，优化坝体结构等措施以增强坝体安全储备，保证拦砂坝的长期可靠性。

［1］韦方强，胡凯衡，J L Lopez.泥石流危险性分区及其在泥石流减灾中的应用［J］.中国地质灾害与防治学报，2007，18(1):23-27.WEI Fangqiang，HU Kaiheng，J L Lopez.Debris flow risk zoning and its application in disaster mitigation［J］.The Chinese Journal of Geological Hazard and Control，2007，18(1):23-27.

［2］唐川，朱大奎.基于GIS技术的泥石流风险评价研究［J］. 地理科学，2002，22(9):300-304.TANG Chuan，ZHU Daku.Assessment of debris flow risk ofYunnan province by using GIS ［J］.Scientia Geographica Sinica，2002，22(9):300-304.

［3］刘涌江，胡厚田，白志勇.泥石流危险度评价的神经网络法［J］.地质与勘探，2001，37(2):84-87.LIU Yongjiang，HU Houtian，BAI Zhiyon.Artificial neural network method for evaluating the dangerous degree of bebris flows［J］.Geology and Prospecting，2001，37(2):84-87.

［4］王光谦，倪晋仁，张军，等.泥石流的颗粒流模型［J］. 山地研究，1992，10(1):1-10.WANG Guangqian，NI Jinren，ZHANG Jun，et al.Grainular flow model for debris flow［J］.Journal of Mountain Research，1992，10(1):1-10.

［5］倪晋仁，王光谦.泥石流的结构两相流模型:I.理论［J］. 地理学报，1998，53(1):66-76.NI Jinren，WANG Guangqian.Conceptual two-phase flow model of debris flow:Ⅰ Theory［J］.Acta Geographica Sinica，1998，53(1):66-76.

［6］马宗源，廖红建，张骏.Bingham型黏性泥石流流体的三维数值模拟［J］.西安交通大学学报，2008，42(9):1146-1150.MA Zongyuan， LIAO Hongjian， ZHANG Jun.Three dimensional numerical simulation of Bingham viscous debris flow fluid［J］.JournalofXi’an Jiaotong University，2008，42(9):1146-1150.

［7］何菊明，彭云雄.液-固两相流的运动描述与数值模拟［J］.湖北师范学院学报(自然科学版)，2006，26(3):26-29.HE Juming，PENG Yunxiong.Motion description and numerical simulation of fluid-solid two-phase flows［J］.Journal of Hubei Normal University(Natural Science)，2006，26(3):26-29.

［8］陈春光，姚令侃.泥石流与主河交汇区三维数值模拟［J］.重庆交通学院学报，2006，25(2):61-65.CHEN Chunguang，YAO Lingkan.Three dimensional numerical simulation of confluence of debris flow and main river［J］.Journal of Chongqing Jiaotong University，2006，25(2):61-65.

［9］刘学，王兴奎，王光谦.基于GIS的泥石流过程模拟三维可视化［J］. 水科学进展，1999，10(4):388-392.LIU Xue， WANG Xingkui， WANG Guangqian.GISBased 3-D visualization of debris process simulation［J］.Advances in Water Science，1999，10(4):388-392.

［10］王纯祥，白世伟，江崎哲郎，等.泥石流的二维数学模型［J］.岩土力学，2007，28(6):1237-1241.WANG Chunxiang，BAI Shiwei，ESAKI Tetsuro，et al.Two-dimensional mathematical model of debris flow［J］.Rock and Soil Mechanics，2007，28(6):1237-1241.

［11］王福军.计算流体动力学分析 CFD软件原理与应用［M］.北京:清华大学出版社，2004:272.WANG Fujun.Analysis of computational fuild dynamics analysis［M］.Beijing Tsinghua University press，2004:272.