1、研究背景

泥石流是山区特殊自然环境在地质运动过程中伴生的、或与人类不合理工程活动共同作用的、对人类生产活动和生存发展具有严重威胁的一种地质灾害[1]。泥石流具有迸发性、时间短、流速快、危害严重等特征[1]。目前,国内外有很多学者对泥石流灾害进行了研究。舒安平等[2]通过分析泥石流的沟道坡度变化,采用室内物理模拟实验的方法,探究了不同坡度下泥石流的威胁区范围、堆积区泥位等特征变化规律。左正金等[3]基于室内物理试验获取了泥石流相关力学参数,利用PFC2D离散元模拟软件对其运动过程进行数值模拟分析,研究了泥石流灾害启动及形成的过程。宋兵等基于现场调查和泥石流形成机理综合分析,运用RAMMS数值模拟,探究了在不同频率条件下泥石流的运动状况[4]。高云建等以汶马高速的泥石流为例,利用遥感解译和野外调查获取相关数据,从理论上分析了泥石流的危害结果[5]。林斌等利用现有青林沟野外资料,考虑3种影响因素开展室内模型实验,得到了青林沟泥石流的一次性固体冲出量的计算模型[6]。田树峰等以九寨沟至绵阳高速公路7条曾暴发泥石流为例,计算了泥石流运动参数特征并分析了其危害模式,提出了相应的治理措施[7]。

这些前人研究成果多数在研究泥石流灾害的启动机理、物理动力学参数等方面,对于泥石流灾害的研究多局限于某个点的结果讨论分析,缺乏对泥石流连续运动过程的描述。本文选取克枯沟泥石流进行分析,首先对该沟泥石流的动力学参数形成因素、特征性质等进行了分析,然后运用CFX软件对20a一遇频率下克枯沟泥石流堆积区的运动连续过程进行数值模拟,通过生成连续流动时的堆积范围、速度云图等,预测泥石流威胁区范围,为泥石流灾害预警、防治及大桥工程施工提供可靠的安全保障。

2、克枯沟泥石流概况

克枯沟位于拟建汶川至马尔康高速公路K54+900～K55+150段,属四川省阿坝藏族自治州汶川县克枯乡,是杂谷脑河左岸的一级支流,距离汶川县城约3km,沟口坐标为103°33′50.60″E,31°30′42.69″N,国道G317线从泥石流堆积扇上穿过。杂谷脑河系岷江上游的一级支流,发源于鹧鸪山脉,于汶川县汇入岷江。2013年,克枯沟突发一次“7·10”山洪泥石流,全乡1100余亩农田被淹没,房屋损坏200余户,国道受损交通中断,直接经济损失达2000万元以上。

据2018年克枯乡统计资料[8],克枯沟沟口流域内现分布有克枯、大寺、木上3个乡村,人口2000多人,居民区房屋460户,耕田106.67hm2,以及正在修建的汶马高速公路克枯特大桥,如图1所示。经现场调查,若再次爆发大型泥石流,直接经济损失将高于2000万元,修建的克枯特大桥也会受到冲击,且可能堵塞杂谷脑河,从而引发二次灾害。

图1克枯沟威胁对象

克枯沟流域属于崇山峻岭,沟谷横断面形态以“V”型谷为主,流域最高高程3824m,最低高程1358m,相对垂直距离达2466m。克枯沟流域呈树叶状,右岸较左岸高耸。流域内树木较少,基岩出露较多,植被覆盖率较低。克枯沟流域的地形地貌十分有利于泥石流的形成,如图2所示。

流域内出露地层单一,主要为志留系茂县群第二组、第三组及第四系地层,岩性主要为灰绿色、深灰色千枚岩夹砂质灰岩、钙质砂岩,深灰色千枚岩夹灰岩、石英砂岩。汶川县多年平均降雨量为513.4mm,全年84%的降雨量集中在5～10月份,多年平均日降雨量达25.0mm,24h最大降雨量达79.9mm[9]。

在对泥石流的研究分析与防治中,沟床的平均比降是一个关键的影响因素指标,沟床平均比降反映出泥石流沟道的坡度。克枯沟泥石流流经的面积为21.70km2,主沟全长为9.20km,泥石流沟道比降平均值为212.1‰。克枯沟剖面如图3所示。

图2克枯沟全貌

图3克枯沟剖面

根据现场堆积扇沿途调查分析,克枯沟流域的物源主要可分为3类:沟道堆积物、崩滑堆积物和坡面侵蚀物。其中,坡面侵蚀物最为丰富,是形成泥石流的主要物质来源。现场调查数据表明,克枯沟流域共有松散固体物质769.34万m3。虽然全流域松散固体物质十分丰富,但是,这些物源并非全部会参与泥石流活动,且不会一次补给产流,而是逐次补给,多次参与泥石流活动。克枯沟流域可能参与泥石流的物源为96.39万m3,统计分析结果详见表1。

据现场调查,“5·12”汶川地震前,克枯沟近100a未发生过泥石流,受汶川地震和灾害发生前的强降雨影响,于2013年7月10日,克枯沟突发了一次山洪泥石流。根据现场调查的特征分析和查阅当年降雨资料,判定克枯沟于2013年突发的泥石流频率为20a一遇。

3、泥石流数值模拟

3.1 CFX软件及分析过程

本次数值模拟选用的计算机软件ANSYSCFX是一款基于有限体积法的高性能、应用广的软件,现多用于各种各样的流体力学科研领域。应用后处理的结果,对相关流体力学运动进行系统分析[10]。本文CFX分析过程包括如下几个阶段(见图4):

(1)使用Rhinoceros软件建立堆积区地质模型,使用ANSYS中的软件对模型网格进行划分。

(2)前处理(CFX-Pre),对泥石流参数流变模型、边界条件及计算步数等参数进行设置。

(3)求解(CFX-Slover),对控制方程、监视收敛等求解过程进行参数设置。

(4)后处理(CFD-Post),以云图等多种方式显示结果,也可以生成动画或动态图来描述瞬态模拟的连续变化过程,可直观地分析和研究结果。

图4数值模拟流程

3.2 泥石流运动参数及流变模型的确定

3.2.1运动参数的确定

在数值模拟前处理(CFX-Pre)部分,需要确定的泥石流动物理力学参数有泥石流的容重、进口流速和峰值流量[11]。

(1)泥石流容重。

参照DZ/T0220-2006《泥石流灾害防治工程勘查规范》等书籍[12],泥石流100a一遇频率条件下的容重与易发性数量化评分存在着对应关系,根据克枯沟流域实际情况,对其进行评分,得分为101分,查表得克枯沟100a一遇频率下容重为1697kN/m3。

泥石流的容重通常随着其规模的变化而变化,因而泥石流的性质也会随着规模的变化而变化,一般地,泥石流容重随泥石流规模增大而增大。查阅泥石流容重相关资料,由于泥石流容重具有一定的时间概念,设计泥石流容重和泥石流暴发频率之间存在有一定的函数关系[13],即:

式中,γ′c为不同频率泥石流的容重,kN/m3;γc为100a一遇泥石流容重,kN/m3;p′为暴发频率系数;p为泥石流暴发周期,a。

前面已确定了100a一遇泥石流容重值,因此可计算得到不同频率条件下的克枯沟泥石流容重值,具体数值如表2所示。

通过公式计算,得20a一遇的泥石流重度1501kN/m3,即本次数值模拟取值。

(2)泥石流流速。

由于克枯沟泥石流属于稀性泥石流,参照《中国泥石流》和《泥石流防治指南》,根据调查区泥石流特征,本次研究采用以下公式计算流速[9],流速计算数值如表3所示。

式中,Mc为稀性泥石流沟床糙率;a为阻力系数;Hc为泥石流水力半径,m,通常用泥深来代替;Ic为泥石流水力坡度;γc为泥石流流体容重,kN/m3;γw为水的容重,取10kN/m3;γs为泥石流固体物质的容重,kN/m3。

(3)泥石流峰值流量。

目前峰值流量主要使用配比法来计算,现场考察沟道的淤积情况,选择所对应的泥石流堵塞系数Dc[12],从而计算峰值流量(见表4)。

Qc=(1+φ)QpDc(6)

式中,Qc为频率为P的泥石流峰值流量,m3/s;Qp为频率为p的暴雨洪峰流量,m3/s;φ为泥石流泥沙修正系数。

φ=(γc-γw)/(γs-γc)(7)

式中,γc为泥石流的容重,kN/m3;γw为清水的容重,kN/m3;γH为泥石流中固体物质比重,kN/m3;Dc为泥石流堵塞系数。

查阅当年降雨资料,已知克枯沟2013年泥石流暴发规模为20a一遇。所以选择克枯沟泥石流峰值流量数值模拟参数为93.82m3/s。

3.2.2泥石流流变模型的设定

CFX中有多种流变模型供用户选择,本次数值模拟选用经典流变模型—宾汉模型[15],其数学模型为

τ=τB+ηdudyτ=τB+ηdudy(8)

式中,τB为宾汉极限剪应力;η为刚度系数;dudydudy为流体的流速梯度。

3.3 假设条件

本次泥石流连续运动过程数值模拟假设有[14]:①泥石流为流体力学性质完全一致的单一液相流体;②泥石流的坡体、沟道不发生变形。

3.4 建立模型和划分网格

泥石流沟谷的地质模型属于复杂实体模型,无法由ANSYSCFX软件快速生成。因此建模过程中选择使用了AutoCAD、Rhino等软件,遵循从线到面、从面到体的建模方式。

这次数值模拟构建的地质模型是克枯沟堆积区,故只建立500m×750m的堆积区地质模型,模型的高度为15m。地质模型网格划分调用ANSYS中的CAE前处理软件ICEMCFD进行,采用更为适合的四面体网格划分法,检查网格质量后,并单独对堆积区模型的进口、出口以及高度5m以上区域进行加密处理,生成节点数104779个,单元数586421个,划分好的网格模型如图5所示。

3.5 前处理

CFX中有5种边界类型选择,分别为进口式、出口式、墙壁式、开放式、对称式[16]。在克枯沟堆积区模型前处理部分中,进口式边界为高程较高的沟道入口,出口式边界为两侧平缓地段,堆积区上表面为开放式边界,沟道、坡体均设置为墙壁式边界。设置完毕后的模型如图6所示。

图5模型尺寸及网格划分

图6边界设定示意

在前处理流体域设定里,进口边界细节选项设置为质量流量,泥石流质量流量随时间变化的函数关系在Expressions中定义[9]。此次数值模拟的模拟情况为在20a一遇泥石流流量下,克枯沟堆积区泛滥的连续过程,总时间选定为600s。其中初始时刻t=0s流量取洪水流量,t=230s取泥石流峰值流量,泥石流流量随时间变化关系如图7所示。

图7流量随时间变化曲线

出口边界细节选项设置为staticpress,数值为1atm;开放式边界细节选项设置为openingpressureanddirection,数值为1atm;墙壁边界细节选项设置为Roughwall,粗糙度数值为0.2m[10]。

3.6 求解

本次求解为瞬态问题求解,数值模拟全部时长为600s,每一步计算时间为1s,每保存一次计算分析结果的时间设置为5s。在求解处理(CFX-Slover)部分需要设定监视泥石流在进、出口的质量流量曲线变化,通过分析曲线变化,进一步验证数值模拟的可靠性。

求解后得到质量流量与时间步长关系曲线如图8所示。从图8中可以看出,进口边界流量曲线与图7设定的流量随时间变化趋势基本一致,说明进口质量流量与设定一致;出口质量流量曲线在时间上具有滞后特征,说明泥石流从入口流入到从出口流出有较长时间;在t=150s时,出口处曲线开始发生变化,说明此时从出口处有泥石流物质流出;在t=340s时,出口质量流量数值达到顶峰,说明此时出口边界的峰值质量流量出现。进口出口质量流量出现差异,原因是由于流场内部储存了一定量的泥石流。

图8进、出口质量流量曲线

3.7 后处理分析

CFX的后处理内容多样、功能齐全,可以通过云图等多种方式反映整个流场连续运动的细节过程。本次数值模拟结果主要通过变化的流速云图来分析泥石流连续运动的过程。根据图9～13分析泥石流流速随时间的变化。

图9t=40s时速度云图

从时间尺度分析泥石流泛滥过程,泥石流从进口流入,在0～40s时间段,泥石流主要受地形限制,流体沿沟道运动;在t=40s时,泥石流泛滥至威胁区域,威胁对象主要为耕田;在t=70s时,泥石流泛滥至对面坡岸,后迅速往两侧扩散,威胁对象扩展为耕田、国道及居民生活区;在t=150s时,泥石流开始从出口边界流出;在t=340s时,泥石流堆积范围达到最大面积。在70～600s时间段,因两侧地形更为平坦,泥石流不断向两侧运动,导致受灾范围不断扩大,经济损害更为严重;在t=600s时,最终泥石流停止运动。

图10t=70s时速度云图

图11t=150s时速度云图

图12t=340s时速度云图

图13t=600s时速度云图

从时间尺度分析泥石流速度变化:在0～70s时间段,泥石流沿沟道逐渐向下运动,速度不断增大,最大速度出现在沟道入口及中段落差较大处,符合实际;在70～150s时间段,泥石流往两侧扩散,两侧扩散速度明显小于平均速度;在t=150s时,泥石流平均、最大流速均达到最大值,最大流速从10.89m/s扩大至11.74m/s;在150～600s时间段,流速不断减小,流体基本在原先流动范围内流动。从流速云图可看出泥石流速度变化是一个连续的过程,且速度分布也可直观分析。数值模拟泥石流平均速度为4.91m/s,与之前计算数值4.88m/s接近,说明利用数值模结果中的速度值是可行的,从而进一步验证了数值模拟的可靠性。

从力学角度分析:流体的速度变化主要是因为泥石流的重力势能迅速减小,从而导致动能迅速加大。以t=150s的泥石流流速云图(见图11)来分析:在泥石流流动的连续过程中,流速在流经区域沟道入口和中段数值最大,最大数值为11.74m/s,这是由于这两处落差较大,在重力加速度作用下,且泥石流流量不断增加,导致这两处出现最大速度值。由于公路平整,没有高低起伏,泥石流流速明显减缓。

根据现场通过对多位克枯乡当地居民的走访以及在坡表找到泥石流暴发后留下的零星松散堆积物,划定2013年实际堆积范围,将其与数值模拟泥石流流经范围进行对比(见图14)。从平面图可以看出,在坡体沟道位置,数值模拟20a一遇泥石流威胁范围与2013年实际堆积范围基本吻合;而在居民生活区,由于建筑物及耕田等影响,堆积范围比模拟威胁范围略微缩小。这主要是因为数值模拟考虑的是地形高度起伏,而忽略其他因素导致的,但堆积范围在模拟威胁范围内。

图14模拟威胁范围与实际堆积范围对比

4、结论

(1)根据现场调查结论,克枯沟泥石流20a一遇频率下容重为1501kN/m3,属稀性泥石流;据调访,克枯沟地震前近100a来未发生过泥石流,故判其泥石流暴发频率为低频。可将克枯沟泥石流类型综合概括为沟谷+降雨+低频+大型+稀性泥石流。由于受2008年汶川地震影响,流域内岩土体结构被破坏,其泥石流活动可能会进入一个活跃期。2013年7月10日克枯沟已暴发了一次泥石流,从现在土源条件、物源动储量来看,克枯沟仍具备发生大规模泥石流的可能。

(2)数值模拟分析得出:20a一遇频率条件下的克枯沟泥石流在堆积区的平均流速为4.91m/s,最大流速为11.74m/s,最大流速出现在堆积区沟道入口及中段处,这是由于这两处落差较大,重力势能改变引起的;在堆积区公路段由于地势平整,泥石流流速减缓。分析模拟结果表明:泥石流流速与实际值接近,而得到的20a一遇泥石流威胁区范围与2013年堆积范围也基本吻合,所以数值模拟结果可靠度较高。

(3)由于泥石流危害性大,破坏力强,对正在修建的克枯特大桥具有威胁性,因此建议可在克枯沟流域中上游设置雨量监测站,避免雨天在泥石流威胁区进行施工。同时在暴雨等极端天气情况下,应注意观察沟道汇水情况;建议在实施大桥修建工程时应加强桥梁墩台防撞处治,高速公路选择大跨桥梁及必要的排导和防撞工程,以避免工程结束后被泥石流冲击破坏甚至坍塌。

参考文献：

[1]孙明书.云南省永胜县泥石流易发性分区评价[D].昆明:昆明理工大学,2016.

[2]舒安平,张欣,唐川,等.不同坡度条件下非均质泥石流堆积过程与特征[J].水利学报,2013,44(11):1333-1337,1346.

[3]左正金,李才,王硕楠,等.基于PFC2D软件的柿树沟泥石流运动过程数值模拟[J].水电能源科学,2014,32(9):135-138,152.

[4]宋兵,沈军辉,李金洋,等.RAMMS在泥石流运动模拟中的应用:以白沙沟泥石流为例[J].泥沙研究,2018,43(1):32-37.

[5]高云建,陈宁生,田树峰,等.低频泥石流对西南山区交通干线的危害与防治:以汶马高速理县段为例[J].人民长江,2018,49(4):23-28.

[6]林斌,张友谊,罗珂,等.沟道松散物质起动模型试验及冲出量预测:四川省以北川青林沟为例[J].人民长江,2019,50(5):113-118,126.

[7]田树峰,陈宁生,高云建,等.九绵高速平武段泥石流运动参数特征与工程危害[J].人民长江,2018,49(11):64-70.

[8]国家统计局农村社会经济调查司.中国县域统计年鉴·2018(乡镇卷)[M].北京:中国统计出版社,2019.

[9]刘希林,吕学军,苏鹏程.四川汶川茶园沟泥石流灾害特征及危险度评价[J].自然灾害学报,2004(1):66-71.

[10]胡卸文,刁仁辉,梁敬轩,等.基于CFX的江口沟泥石流危险区范围预测模拟[J].岩土力学,2016,37(6):1689-1696..

[11]贺拿,陈宁生,谢万银,等.雅鲁藏布江街需水电站坝址区九锥沟泥石流发育特征[J].中国地质灾害与防治学报,2015,26(3):57-64.

[12]中国地质灾害防治工程协会.泥石流灾害防治工程勘查规范:DZ/T0220-2006[S].北京:中国标准出版社,2006.

[13]沈娜.四川省九龙县石头沟泥石流特征与防治工程措施研究[D].成都:成都理工大学,2008.

[14]胡桂胜,陈宁生,李俊,等.金沙江白鹤滩水电站坝址近场区泥石流运动特征与发展趋势分析[J].水土保持研究,2014,21(2):238-245.

[15]刁仁辉.丹巴县江口沟泥石流发育特征及沟口危险区范围预测研究[D].重庆:西南交通大学,2015.

[16]韩玫.汶川震区“宽缓”与“窄陡”沟道型泥石流致灾机理研究[D].重庆:西南交通大学,2016.

[17]费祥俊,舒安平.泥石流运动机理与灾害防治[M].北京:清华大学出版社,2004.

[18]贾涛,唐川,王纳纳.基于FLO-2D与冲量模型的泥石流危险度分区方法及应用[J].水电能源科学,2015,33(2):152-155,171.

[19]薛强,曾庆铭,徐继维.都江堰市龙池镇煤炭坪泥石流形成条件及特征[J].自然灾害学报,2015,24(2):115-121.

[20]熊斌.安县干沟桥泥石流成因及发展趋势预测研究[D].成都:成都理工大学,2010.

[21]庄建琦,葛永刚,陈兴长.震后都汶公路沿线泥石流沟堵江危险性评估[J]工程地质学报,2012,20(2):195-203.

[22]韩玫,胡涛,王严,等.汶川震区窄陡沟道型泥石流动力学特性及堵江分析:以都汶高速沿线磨子沟为例[J].工程地质学报,2016,24(4):559-568.

[23]胡桂胜,陈宁生,游勇,等.“7.10”连山大桥泥石流运动特征及沟道堵溃分析[J].成都理工大学学报(自然科学版),2015,42(6):643-648.

[24]高云建,陈宁生,赵春瑶,等.干热河谷区泥石流活动特征及对交通干线的危害:以峨汉高速为例[J].人民长江,2018,49(10):54-59.

[25]杨海龙,樊晓一,姜元俊.基于FLO-2D的泥石流工程治理效果分析:以都江堰银洞子沟泥石流为例[J].人民长江,2017,48(17):66-71,95.

[26]焦亮,柳金峰,游勇,等.泥石流沟危险度定量评价研究:以白龙江陇南两河口镇至桔柑乡段为例[J].人民长江,2019,50(6):118-123.

胡佳武,李桥,黄健.基于CFX的克枯沟泥石流动力学特性分析[J].人民长江,2020,51(07):26-32.

基金:四川藏区高速公路重大地质灾害监测预警系统及应用示范研究项目(藏高科技2016-10);国家重点研发计划项目(2016YFC08022003).