基于FLOW-3D,的泥石流龙头运动过程模拟研究

罗超鹏，常 鸣，武彬彬，刘沛源，余波

（地质灾害防治与地质环境保护国家重点实验室（成都理工大学），四川 成都 610059）

我国西南地区主要以山地丘陵地貌为主，地质条件极为复杂，雨季雨量较大，尤其在近几年西南地区多次暴发地震之后，山体变得极为破碎，在一些泥石流流域范围内中经常发生崩塌、滑坡等地质灾害，大量的松散物质堆积于沟道之中，为泥石流的形成提供了充足的物源条件[1]。在四川西部强震区一些典型的沟道内存在大量的松散堆积物，致使泥石流暴发频率急剧增加，规模增大，成为震后主要的地质灾害类型[2]。因此，准确分析其全过程运动特征、量化其强度指标、提出科学合理的防治措施对今后泥石流防灾减灾工作的开展具有重大意义。

随着国内外学者近年来对泥石流运动特征展开的深入研究，目前已经有多种方法开展泥石流运动特征分析工作，常用的方法有数值模拟法[3]、统计经验模型法[4]、GIS 与经验模型相结合的预测方法[5]。其中数值模拟在泥石流运动特征研究中的应用逐渐广泛，模拟精度逐渐提升。近年来离散元方法[6−7]、光滑粒子方法[8−9]、格子波尔兹曼方法[10]、连续介质力学方法[11]等数值模拟技术大量应用于泥石流灾害动力学过程研究和危害范围评估。Begueria等[12]采用基于泥石流动力学的二维有限差分流体动力学MassMov2D 模型，借助PCRaster 软件实现在复杂地形中模拟泥石流运动。Christen等[13]提出用连续流体介质方法的RAMMS 模型来调节泥石流的质量及运移速度，实现运动过程模拟。He等[14]利用运动控制方程及临界降雨量值模拟了关山沟泥石流的侵蚀过程。Chang等[15]利用FLO-2D计算模型模拟了都江堰龙溪河流域群发性泥石流运动堆积过程。Nocentini等[16]结合DAN-W 与FLO-2D 模型运用动态代码对同一地区潜在泥石流建模，通过动态模型快速评估研究地区潜在泥石流风险。Ouyang等[17]利用Massflow 模型模拟了甘肃武都泥湾泥石流动力运动过程。Horton等[18]利用voellmy 铲刮模型开展震后泥石流模拟。Trujillo-vela等[19]运用光滑粒子流体动力学（SPH）和离散元方法(DEM)建立了耦合数值模拟分析模型。Liu[20]等利用无限边坡稳定模型、增强的两相泥石流模型等组合，研究小流域内物源启动型泥石流的形成和传播过程。为了逐步优化泥石流评价结果，崔鹏等[21]建立了泥石流风险评估体系、熊江等[22]通过对强震区多期精细化遥感解译，优化了西南山区典型的泥石流危险性评价指标。

综上可见，我国西南山区泥石流数值模拟研究主要基于二维数学模型，在三维数学模型对泥石流的运动堆积过程的模拟方面还存在不足，文中根据泥石流基本特征，开展二经里沟泥石流龙头运动过程数值模拟分析，揭示山区典型泥石流灾变过程，为后期防灾减灾工作提供科学依据。

2020 年10 月3 日，二经里沟发生泥石流灾害，破坏高速公路隔离带护墙，经实际调查分析沟口堆积物约有4 000 m3。在调查中发现沟内物源仍然较丰富，在未来极端降雨条件下仍有泥石流暴发的可能，二经里沟一旦再次暴发，将威胁汶马高速与G317，并且极有可能堵塞河道，形成灾害链，对下游居民生命财产安全造成巨大威胁。因此，为了有效降低泥石流带来的损失，文中以二经里沟为研究对象，基于前期调查工作所取得的泥石流各项成果，采用FLOW-3D 数值计算模型，开展该泥石流的动力学运动特征分析，并提出合理的泥石流防灾减灾对策。

二经里沟位于四川省理县米亚罗镇鹧鸪山来苏河流域，汶马高速西侧，沟口与河道连接（图1）。二经里沟流域面积1.82 km2，沟长2.47 km，最高点海拔4 044 m，沟口海拔2 964 m，高差达1 080 m，沟道平均纵比降为437.2‰，主沟横断面呈“V”字形。二经里沟地处青藏高原东缘与四川盆地西北边缘交错接触带，地形复杂，为高山峡谷区域，山势陡峭。区域主要断裂有茂汶断裂、米亚罗压扭性断层、松岗压扭性断层，复杂的地质构造致使山体岩层整体性较差，为泥石流的形成提供了充足的物源条件，是泥石流形成的关键条件之一[23]。研究区内新构造运动较活跃，属地震强烈和邻区强震波及区，地震动峰值加速度为0.15g，测区基本地震烈度为Ⅶ度。

图1 二经里沟流域概况图Fig.1 The overview map of Erjingli gully watershed

通过无人机航摄以及详细调查，确定了沟道内主要有8 处较大规模的物源分布。2020 年10 月3 日暴发泥石流之后，流域内物源储量显著增加，其中沟道物源方量增加较为明显。在上次泥石流暴发过程中使沟道两侧的崩滑体复活，今后在降雨和重力的作用下不稳定的崩滑体将会持续向沟道运动，极有可能堵塞沟道，放大泥石流的灾害效应。图2（a）为位于沟口以上275 m 处右岸的滑坡，滑动距离为24 m，宽为45 m，约2 052 m3；
图2（b）为位于沟口以上519 m 处右岸的滑坡，滑动距离为26 m，宽62 m，约3 869 m3；
图2（c）位于沟口以上978 m 位置，沟道左侧后缘发育的裂隙，裂缝宽约5～10 cm，延伸约70 m，一旦失稳，将向沟道内输入约3 654 m3物源；
图2（d）为沟口以上1 618 m 处沟道物源，约504 m3。

图2 沟道地形地貌以及典型物源照片Fig.2 The gully landform and typical provenance photos

二经里沟流域属于大陆性高原季风气候，根据气象部门监测数据，该区域年平均降雨量为603.7 mm，日最大降雨量51.2 mm，2020 年米亚罗镇降雨明显较往年偏多，尤其是10 月3 日突发暴雨，平均每小时降雨量达到23.4 mm（图3），流域内裸露区在累计降雨作用下产生坡面崩滑体，导致二经里沟暴发泥石流，流体中含有大量黏土颗粒，泥浆呈黏性，并携带少量树枝及大粒径块石。泥沙高速冲出沟口，破坏力极强，持续约12 min，造成汶马高速被阻断约100 m，且沟口桥梁被横向推移20 cm 以上。泥沙在来苏河河道之中形成大规模堆积，最大堆积深度约6 m，且堵塞河道约150 m。

图3 泥石流暴发期间降雨统计Fig.3 Rainfall statistics during the outbreak of debris flow

2.1 模型概述

二经里沟泥石流具有浓度高、黏性大、边界条件明显的特点，固体颗粒和流体存在着较为强烈的耦合关系。FLOW-3D 数值模拟软件中颗粒流与泥石流流体存在较大的相似性，因此利用颗粒流模型开展二经里沟泥石流龙头运动过程模拟能够更加真实的反映其实际情况。

FLOW-3D 是基于有限差分法对控制方程进行求解的数值计算软件，在模拟过程中将计算区域设置三维网格，每一个网格在计算过程中都有多个相关变量的值在每一个网格上体现，其中网格表面包括三个方向的速度u、v、w，网格中心包括压强p、流体体积分数F、密度ρ、湍流动能q、紊流动能消耗率以及黏度µ 等变量，能够更加真实的还原泥石流实际运过过程。研究利用颗粒流模型、漂移-通量模型、重力和非惯性框架模型以及黏流和湍流模型相互组合，分析泥石流形成、运动、堆积特征，其中控制方程包括连续方程和动量方程，见式（1）—（4）。

Ax、Ay、Az——流体所流经单元格x、y、z三个方向 上的面积比例/%；

Gx、Gy、Gz——流体在网格x、y、z三个方向上的加 速度/（m·s−2）；

fx、fy、fz——流体在网格x、y、z三个方向上的阻力 加速度/（m·s−2）；

VF——流体体积分数；

ρ——流体密度/（g·cm−3）；

p——作用在流体上的压强/Pa。

模拟中的颗粒流模型是将物源考虑为若干个三维球体颗粒以及液体的混合物，分析每一个三维球体在运动过程中的受力情况。当颗粒之间的静力平衡被破坏时，粒子产生运动，泥沙运动特征由内部粒子之间状态变化体现出来。由于实际泥石流流体的非均一性，泥石流在运动过程中颗粒与颗粒、颗粒与液体之间也必将存在相对运动。为了使模型可以允许粒子与粒子之间发生有限位移，并且各粒子之间能够自动辨识新的接触关系，因此在研究中将固体颗粒和液体之间定义为自由表面边界来近似表达黏性泥石流运动特征。由于固体颗粒和周围液体流动速度不同，混合物的两相效应可以用漂移-通量模型近似分析，在漂移通量模型中，可以使用一个选项来定义可能发生的最大固体分数，当达到极限体积分数时，迫使其漂移速度为零。

2.2 泥石流基本参数及三维建模

式中：u、v、w——x、y、z三个方向上的速度/（m·s−1）；

此次在二经里沟堆积区、流通区及形成区各取一组具有代表性的样品，针对每组样品在室内进行颗分试验，绘制颗粒级配曲线（图4）。由级配曲线可以分析得出，堆积区土样颗粒较为均匀，颗粒普遍偏细，平均粒径为1.8 mm；
流通区土体颗粒粒径较堆积区更大，大小不均，级配相对较好，平均粒径为4.3 mm；
形成区土样中大颗粒偏多，土颗粒不均匀，级配良好，颗粒的平均粒径为10.1 mm。

图4 二经里沟各样品级配曲线Fig.4 The grain size distribution of Erjingli gully

通过采取多组有代表性的泥石流堆积物样品加水搅拌成泥石流浆体，测量计算得到泥石流流体的重度以及泥沙体积浓度（式5—6）：

式中：γc——泥石流流体重度/（t·m−3）；

Gc——配置浆体重量/t；

V——配置浆体体积/m3；

CV——泥沙体积浓度；

γs——泥石流中固体颗粒比重/（g·cm−3）；

γw——水的重度/（g·cm−3）。

通过计算得到γc=1.82 t/m3，CV=0.5。

由于缺乏实测数据，粗糙系数n值按照前人给出的天然河道的粗糙系数参考值确定，文中结合二经里沟实际情况取n=0.18[24]。其余参数均采用软件默认值。

通过物理实验及相关分析，在数值分析的颗粒流模型中，泥石流特征控制参数见表1。

表1 二经里沟特征控制参数Table 1 Main parameters of numerical simulation of Erjingli gully

基于无人机航拍以及人员实地调查，对研究区物源规模及分布特征进行了较为精细的记录。但是在调查过程中发现除了少量沟道物源可以根据泥痕测量其厚度以外，大多数地方无法获得物源的空间厚度数据。在本次模拟中参考前人研究计算各物源平均厚度[25]（式7）：

式中：H——物源平均厚度/m；

A——物源面积/m2。

在详细调查中实测了1、2、3、4 处物源的实际体积，其他物源在调查时由于野外测量较为困难，因此模拟中采用经验公式计算确定。最终确定各物源的平均厚度及体积见表2，将计算结果引入Rhino 三维建模软件，构建流域物源三维模型（图5）。

图5 二经里沟地形及物源三维模型示意图Fig.5 Three-dimensional model of topography and provenance of Erjingli gully

表2 二经里沟物源面积、平均厚度及体积Table 2 Provenance area and average thickness of Erjingli gully

2.3 边界条件

在模拟过程中，需要根据实际情况考虑各网格合理的边界条件。一共设置了6 块计算网格，所有的边界条件设置均满足实际要求，包括出流边界O、壁面边界W、连续边界C以及对称边界S等4 种边界类型（表3）。

表3 模拟边界条件设定结果Table 3 Boundary condition setting results of numerical simulation

3.1 泥石流运动特征分析

为了保证模拟结果的完整性，需要设定合适的模拟时间，文中根据余斌[26]提出的黏性泥石流运动平均速度经验公式，见式8，计算出泥石流在沟道中的平均运动速度，进而估算泥石流龙头运动过程的时间。

式中：U——黏性泥石流运动平均速度/（m·s−1）；

g——重力加速度/（m·s−2）；

R——黏性泥石流运动水力半径/m；

S——黏性泥石流运动纵比降/‰；

D50——泥沙颗粒中百分比小于50%的颗粒粒径（中 性粒径）；

D20——泥沙颗粒中百分比小于20%的颗粒粒径。

通过计算得到二经里沟泥石流的平均速度为12.41 m/s，最终确定泥石流龙头全过程历时198.97 s，文中计算时间设置为200 s。

考虑到该沟道较窄，此次模拟工作将计算网格设置为1 m×1 m，通过控制方程计算泥沙在每个网格单元的运动参数。对二经里沟2020 年10 月3 日暴发的泥石流开展精细化三维数值模拟，重现该泥石流暴发过程中泥沙形成、运动、堆积的整个运动过程，图6、图7 为不同模拟时间的泥石流运动情况。

图6 二经里沟泥石流堆积深度模拟结果Fig.6 Simulation results of debris flow accumulation depth in Erjingli gully

图7 二经里沟泥石流流动速度模拟结果Fig.7 Simulation results of flow velocity of debris flow in Erjingli gully

t=0～30 s，沟道两岸物源开始启动，松散土体运动速度开始逐渐增大，并且物源失稳边界开始逐渐扩大，两岸泥沙持续向沟道聚集，最大泥深为7.06 m，运动缓慢；
这一阶段属于物源汇集的过程。

t=30～60 s，由于两岸以及上游物源的持续输入，松散固体物质在沟道内部不断汇集，在上游物质强大的动力作用下，沟道内部大规模泥沙开始运动，不断向下游推进，最大泥深和流速分别为8.28 m、11.94 m/s；
这一阶段属于泥石流运动特征突变阶段。

t=60～100 s，由于泥浆在不断向下游推移的过程中也启动了下游沟道内堆积的松散物质，泥石流规模逐渐增大，此时泥石流前端的泥深和流速分别为16.18 m、14.94 m/s；
这一阶段属于泥石流的持续发展阶段。

t=100～130 s，泥石流运动到沟口，强大的冲击力对高速路以及下游其他基础设施造成严重破坏，泥石流冲出沟口之后运动速度急剧降低，形成的堆积体堵塞河道。随着堆积规模的不断扩大，泥沙开始向河流下游运动，增加河道的堵塞程度；
最终泥沙在河道内最大堆积深度为5.99 m，堆积长度为170.6 m，完全堵塞河道，泥石流基本停止运动，其结果与现场调查相符；
这一阶段属于泥石流的泥沙堆积阶段。

模拟结果显示，泥石流最终将在沟口形成最大深度为5.99 m 的堆积扇，野外实地调查结果为6.0 m，泥深模拟结果与实际情况吻合度较高。通过数值模拟获取的二经里沟泥石流堆积危险范围、野外实际测量堆积危险范围及两者重合的危险范围验证泥石流数值模拟的精确程度式如下：

式中：F——模拟精度；

S0—实测与本次模拟结果的堆积范围重叠区域/m2；

SM——实测堆积面积/m2；

SN——模拟堆积范围/m2，具体计算参数见表4。

表4 模拟结果与野外实测对比验证Table 4 The comparation of simulated results and field measurements

经过计算，此次模拟精度达到88.98%，满足数值模拟的精度要求（图8）。其计算结果和实际现象相似，表明此次数值模拟结果对研究该泥石流灾害的运动特征具有较好的参考价值。

图8 二经里沟模拟结果验证Fig.8 Verifying of the simulated results of Erjingli gully

3.2 泥石流强度分析

泥石流强度分析是今后防灾减灾工作中十分关键的步骤，相关研究提出泥石流的强度大小主要由泥石流的堆积深度和流动速度相关[27]。文中主要考虑对承灾体起主要影响因素的泥石流泥深和流速两个关键参数。用强度指数（IDF）[28]开展泥石流堆积区内强度特征分析，见式10，并得到二经里沟泥石流强度空间分布特征图（图9）。

图9 二经里沟泥石流强度空间分布图Fig.9 Spatial distribution characteristics of debris flow intensity in Erjingli gully

式中：IDF——泥石流强度指数/（m3·s−2）；

t——模拟时间/s；

h——泥石流堆积深度/m；

u、v、w——x、y、z三个方向上的流速/（m·s−1）。

根据计算结果，二经里沟泥石流的强度参数取值区间为0～586.8。文中根据归一化之后的值将泥石流强度划分为4 个强度等级，定义I≥0.5 为极高强度区域；
0.5>I≥0.3 为高强度区域；
0.3>I≥0.1 为中强度区域；
0.1>I≥0 为低强度区域。

将计算结果进行统计分析可知，二经里沟泥石流极高强度区域占609 m2，高强度区域占1 295 m2，中强度区域占3 451 m2，低强度区域占20 034 m2。极高、高强度区域主要分布在下游窄陡沟道以及泥沙堆积厚度较大位置。中、低强度区域主要分布在坡度相对较缓的沟道中上游以及堆积区边缘位置。其中，高、极高强度的位置遭受到泥石流冲淤作用强烈，在后期管道、公路等工程选线中应该尽量避开，并且应该考虑在沟道内修建必要的泥石流消能型防治工程设施，降低泥石流强度。

3.3 拟设拦挡坝减灾效果评价

文中研究考虑了在沟道下游拟设一个拦挡坝，运用相同的计算模型分析泥石流在拦挡坝作用下的运动特征（图10）。

图10 拦挡坝作用下泥石流堆积深度模拟结果Fig.10 Simulation results of debris flow accumulation depth under retaining dam

拦挡坝距离沟口300 m，拟定坝高为12 m，将两侧坝肩嵌入坡体，且在模拟中不考虑坝体破坏。结果表明，拦挡坝对泥石流起到了很好的缓冲作用，泥沙将坝体淤满，坝后最大堆积深度达到了11.78 m，坝体在沟道上游堆积长度43 m，约拦挡1 720 m3泥沙，在拦挡坝的作用下，泥沙冲出沟口的流速为2.55 m/s，形成小规模的堆积，最大堆积深度为2.97 m，形成小规模的堆积，造成约1/3 的河道堵塞，降低了对沟口高速公路产生的影响，减灾效果明显。

文中以2020 年10 月3 日二经里沟暴发的泥石流为研究背景，结合无人机航拍摄影以及现场调查取样工作，对二经里沟泥石流发育特征、降雨特征、物源空间分布特征开展详细调查。利用FLOW-3D 对泥石流的运动全过程开展数值模拟分析，并划分了泥石流强度。文中取得的主要结论如下：

（1）采用FLOW-3D 数值模拟的方法，重现了2020年10 月3 日二经里沟泥石流龙头运动过程，根据泥石流运动特征将二经里沟泥石流暴发全过程划分为4 个阶段：物源汇聚、运动特征突变、持续发展阶段、泥沙堆积阶段。在堆积区范围内将模拟结果与实际调查结果进行对比验证分析，计算出模拟精度为88.98%。

（2）基于泥石流最大堆积深度与流速确定其强度评价指标，建立强度等级划分标准，并绘制强度等级划分图。其中，极高强度区占2.4%，高强度区占5.1%，中强度区占13.6%，低强度区占78.9%。在极高、高强度区范围内房屋、道路、河堤及农田容易受到泥石流的破坏，需要加强监测预警及重点防护。

（3）选择距离沟口约300 m 处设置一个高12 m 拦挡坝，模拟结果显示沟口最大堆积深度减小3.02 m，大幅度降低泥石流对沟口高速公路的威胁程度，该方案能为今后防治工程的修建提供参考。

猜你喜欢 沟口物源泥沙 ——以日本经典文学《金阁寺》与电影《火烧金阁寺》为例">文学经典与改编电影的人物塑造差异探析
——以日本经典文学《金阁寺》与电影《火烧金阁寺》为例戏剧之家(2022年4期)2022-11-12泥沙做的父亲中学生天地(B版)(2022年4期)2022-05-17沟口党员生活·下(2020年12期)2020-01-11《金阁寺》中“斩猫”与“纵火”的内在关联探寻文学教育(2019年19期)2019-11-26车西地区陡坡扇体储层发育主控因素探讨科技视界(2019年1期)2019-03-25成都粘土的分层、成因及物源研究综述安徽农学通报(2018年12期)2018-09-18东营三角洲沙三中物源分析探讨中国科技纵横(2018年7期)2018-05-22超大型油船防泥沙设计船舶与海洋工程(2018年1期)2018-03-31第三纪火山沉积硼矿与火山岩关系研究科技资讯(2015年7期)2015-07-02新疆尼勒克县吉仁台沟口墓地和遗址大众考古(2015年9期)2015-06-26

推荐访问:泥石流 龙头 过程