DOI: 10.11817/j.issn.1672-7207.2015.08.040

碱坪沟泥石流运动模型

张宁^{1, 2}，姚磊华¹，曹中兴³，杨卓²，黄鑫²

(1. 中国地质大学(北京) 工程技术学院，北京，100083；

2. 航天科工惯性技术有限公司，北京，100074；

3.中交水运规划设计院有限公司，北京，100007)

摘要：针对碱坪沟特殊的地理环境和地质条件，通过室内模型实验，研究不同降雨强度下入渗强度变化规律。结合入渗规律，建立适用于碱坪沟泥石流的一维运动模型。该模型在圣维南方程组的基础上，考虑沟谷两侧向中心汇流和入渗因素，并用Preissmann四点隐格式差分方法对方程组进行计算。根据“8.13”碱坪沟泥石流爆发当天的小时降雨量，计算泥石流爆发过程中的流速和流深变化情况，并与监测资料进行对比验证。研究结果表明：该模型较好地反演了“8.13”碱坪沟泥石流运动过程中的参数变化情况。

关键词：碱坪沟泥石流；运动模型；模型槽实验；Preissmann有限差分

中图分类号：P642             文献标志码：A         文章编号：1672-7207(2015)08-3061-06

Movement model of Jianping Gully debris flow

ZHANG Ning^{1, 2}, YAO Leihua¹, CAO Zhongxing³, YANG Zhuo², HUANG Xin²

(1. School of Engineering and Technology, China University of Geosciences (Beijing), Beijing 100083, China;

2. Aerospace Science Inertial Technology Co., Ltd., Beijing 100074, China;

3. China Communications Construction Company Water Transportation Consultants Co.Ltd, Beijing 100007, China)

Abstract: Based on the special geographical and geological conditions of Jianping Gully, a set of laboratory model tank experiments were carried out to study the runoff and infiltration regularity. On the basis of Saint-Venant equations, a movement model of debris flow was built considering the side stream and infiltration. The equations were numerically solved by applying Preissmann four-point linear implicit scheme. According to the actual hourly rainfall happening on August 13th, 2010, the flow depth and speed were calculated and verified by the monitored values. The results show hat the movement process of “8.13 Debris Flow” can be simulated by the model based on parametric inversion.

Key words: Jianping Gully debris flow; movement model; model tank experiment; Preissmann finite difference

碱坪沟位于四川省都江堰市龙池镇龙溪河右岸南岳村4社，流域面积为3.6 km²，河道总长约3 km，河道平均纵坡降为363‰^[1]。主沟河道较为平直，坡度上陡下缓，河道两侧狭窄，上部深切沟谷发育，特别有利于松散固体物源、降水、地下水汇集，这为泥石流的爆发提供了地形条件。沟域内地质构造复杂，断裂极为发育，虹口—映秀断裂的一支从碱坪沟南侧通过。沟域内出露的基岩主要为震旦系安山岩、泥灰岩和玄武岩，经多次构造运动节理裂隙发育，强风化壳厚，可达10余m，在重力和动水压力作用下常形成滑坡、崩塌为泥石流形成提供固体物源，特别是“5.12汶川地震”引发了多处大量的崩塌、滑坡堆积于沟床及两侧坡面，为泥石流暴发储存了丰富物源。固体物质补给集中，数量大，稳定性较差，沟谷时陡时深时狭的地形特点十分有利于泥石流的形成和活动，容易发生泥石流灾害。2010-08-13，该地区形成了短时强降雨，最大雨强为70 mm/h。在暴雨的作用下，碱坪沟爆发了大型泥石流，规模约7万m³，沟口堆积扇淤埋面积约14 600 m²，淤积厚度3~5 m，淤埋河道长度为15~30 m。据现场描述“该泥石流性状犹如混凝土状，多黏稠，其速度根据沟道形态及纵坡降差异较大”，泥石流容重约为2 t/m³，且根据现场抽样实验，黏性土质量分数为12%~15%。综上所述，将“8.13碱坪沟泥石流”定义为大型、沟谷型、暴雨型、黏性泥石流。本文根据碱坪沟特殊的地质、地形条件，结合室内模型槽实验数据及现场资料，对主沟流通区进行数值模拟，力图真实还原“8.13”碱坪沟泥石流启动及运动过程。近些年来，各国专家学者针对泥石流启动、运动、堆积等问题进行了深入研究，取得了大量成果。然而，诸如复杂几何地形、流体的不均匀性、固液两相之间的相互作用、沟谷多方位降雨补给和入渗等因素，大大增加了相关数学模型的复杂程度。目前，人们常用的方法是针对泥石流具体的运动阶段和运动特性，给出相应的模型，并将一些影响不大的项忽略或进行简化处理。较早期的数学模型有Johnson的黏塑模型和Bagnold的膨胀模型^[2]。Johnson^[2]通过模型槽实验研究泥石流浆体的流变特性，认为泥石流是不同于牛顿流体的Bingham流体，增加了宾汉极限剪应力。MacArthur将泥石流作为Bingham流体，运用一维圣维南方程组并添加黏性阻力项，用有限元法对方程组进行了求解。范椿^[3]假设当沟谷长度远大于宽度时，对连续方程、运动方程和Bingham模型的本构方程进行量级比较，推导出一维控制方程。Lin等^[4]进行了颗粒流的膨胀实验，认为流体的剪应力与剪变率的平方有关。Takahashi^[5]提出分散应力的概念，认为在石块型泥石流中，颗粒的碰撞应力远大于黏性力，改进了拜格诺膨胀流模型，建立了惯性粒子流模型，求解了泥石流平均速度和流体深度。1992年，他又研究了侵蚀相的估计，量化了流体运动过程中的侵蚀和堆积作用^[6]。该模型的缺点是没有考虑液相的黏性。Iverson^[7]提出了流体中的黏性孔隙流成分，并推导了二维的固液混合模型，考虑了固体颗粒之间及固体和液体之间的相互作用，并通过实验对其方程组进行了论证。该模型缺点在于需要野外采集诸如流速、流深等参数，这大大增加了模型应用的难度。Iverson等^[8]进一步提出了基于深度均化理论的流-固耦合模型，建立了泥石流仿三维运动模型，模拟了泥石流从启动到堆积的完整过程，并在一定程度上解决了参数取值的问题，但该模型只考虑了泥石流在矩形条件下的运动过程，不适应流通区几何形状发生改变的情况。2004年，他又提出了适应不规则地形条件的计算模型^[9]，但该模型还是没有考虑流体介质的作用，而且计算过程非常复杂，不便于广泛应用。国内一些学者通过区分泥石流固体和液体组分、区域规划、泥位特征值等方法对泥石流进行研究，并建立相应的模型，对泥石流物理性质的研究起到了一定推动作用^[10-14]。上述模型研究的重点不同，并且都有其自身的局限性。碱坪沟泥石流为黏性泥石流，因此，本文作者在推导过程中重点考虑流体中黏性项的影响，忽略颗粒之间的相互碰撞作用。在恒定降雨条件下，考虑入渗及沟谷测流的影响。

1  碱坪沟一维运动模型的建立

1.1  基本假定

1) 泥石流中的固体与液体作为一种混合介质，密度不随时间和空间而改变；2) 碱坪沟泥石流为黏性泥石流，只考虑液相的摩阻坡降；3) 泥石流在地表的运动过程是一维运动，流速与流深在横截面上不发生变化；4) 雨水从沟谷两侧向中心汇集的侧流强度沿纵向均匀分布；5) 入渗为一维垂直入渗且符合室内实验规律，其表达式在实验中获得。

1.2  质量守恒方程

根据Naiver-Stokes方程的守恒形式，当流体密度恒定，其质量守恒方程可以写为

              (1)

式中：h为流体深度，mm；t为时间，s；V为速度向量；M_f为质量流通量。将式(1)简化为一维形式可得

               (2)

式中：u为x方向上的速度，m/s。M_f满足

             (3)

式中：R(t)，I(t)和L(t)分别表示降雨强度、侧流强度和入渗强度。侧流强度为单位时间内沟谷两侧向中心汇流的深度，数值上等于降雨强度减去植物截流强  度f_t：

                (4)

根据Horton经验公式，植物截流强度f_t可以写为

                (5)

式中：为稳定截留强度，mm/h；a为经验系数。

入渗强度为单位时间单位面积上坡面入渗量，根据室内模型槽实验测得。在不同降雨强度条件下(60，90和120 mm/h)进行径流-入渗实验。实验在三峡大学教育部地质灾害重点实验室完成。图1所示为60 mm/h降雨强度下测得的入渗强度曲线。

入渗可分为3个阶段：自由入渗阶段(入渗率为  1 mm/min)、有压入渗阶段(入渗函数为I′(t))和饱和入渗阶段(入渗率为0.22 mm/min)。由图1可知：在50 min前，入渗强度等于降雨强度，没有坡面径流；而在110 min后，入渗趋于稳定，为饱和入渗阶段；在50~110 min内为有压入渗阶段。将有压入渗阶段进行曲线回归分析，如图2所示。

图1  降雨强度60 mm/h时入渗强度曲线

Fig. 1  Infiltration intensity curve at rain fall in tensity of 60 mm/h

图2  有压入渗阶段入渗曲线回归分析

Fig. 2  Regression analyze of infiltration curve in pressured water infiltration period

将50~110 min内测得的20个散点进行曲线拟合，其中生物生长模型的拟合度最高，并得到有压入渗阶段的入渗函数：

        (6)

将60 mm/h雨强条件下的入渗率归纳为分段   函数：

  (7)

同理可得90 mm/h和120 mm/h降雨强度下的入渗函数，将2个入渗函数进行归纳，可得

  (8)

1.3  动量守恒方程

考虑作用在流体微团上的合力，包括两过水断面的压强差、重力在x轴的分量、沿坡面的切应力和流体黏性剪应力，得到运动方程并将其化简为：

      (9)

                (10)

             (11)

其中：μ为流体黏性系数，N·s/m²；θ为流动方向上的坡度，(°)；S_f为流体液相摩阻坡降^[15]；τ为宾汉极限剪应力，N/m²。

最终，适用于碱坪沟泥石流的一维单相流运动方程组为

    (12)

2  模型的计算与验证

2.1  计算方法

运用pressimann四点隐格式差分方法求解一维运动方程组，差分格式如下：

   (13)

     (14)

         (15)

其中，ψ为权重系数，取值为0.6~1.0，项用二阶中心差分来处理，

          (16)

定义同一节点上相邻时段的差值为。根据式(13)~(16)，对运动方程组(12)进行离散并线性化，得到离散方程的线形表达式为

  (17)

其中：A_j，B_j，C_j，D_j，E_j和A_j′，B_j′，C_j′，D_j′和E_j′分别为质量守恒方程和动量守恒方程各项的系数。假设如下线性关系：

      (18)

通过运用双循环迭代法对方程组进行求解，具体求解过程见文献[16-17]。

2.2  计算条件

初始时刻坡面上每个节点的流速和流体深度均为零。

边界条件包括上游和下游2组边界条件。

上游：       

下游：       

模型选取碱坪沟主沟流通区进行计算，主沟流通区长约为1 200 m，平均坡度约为30°，这与室内模型实验的坡度保持一致。流体的物理力学参数均按《都江堰市龙池镇8.13地质灾害应急排查报告》设定，具体参数如下：流体密度ρ=2 t/m³；宾汉剪应力τ=30 N/m²；黏性系数μ=0.3 N·s/m²；权重系数ψ=0.6；△x= 1 m；△t=10 s。

根据模型总长和步长将模型划分为1200个单元，1201个节点，这样共有2400个离散方程和2402个未知量，利用上、下游边界条件正好满足方程组的求解，并运用Delphi-XE软件进行编程计算。

2.3  计算结果

碱坪沟“8.13”泥石流爆发当天，从14:00开始降雨，17:00降雨强度增大至70 mm/h，18:00降雨强度为50 mm/h。设定17:00之前，沟内岩土体已为饱和状态，初始入渗强度调整为饱和入渗。计算时长为2 h，0~1 h的降雨强度为70 mm/h，1~2 h的降雨强度为50 mm/h。碱坪沟主沟流通区已安装2个超声波泥位计，一个位于主沟中部(约x=600 m处)，另一个位于主沟出口处(约x=1 200 m处)。分别计算2泥位计处流体深度变化情况并与监测值进行对比，如图3和图4所示。泥位计算相对误差如表1所示。

图3  主沟中部(x=600 m)泥位计算值与监测值对比

Fig. 3  Comparison of sludge level between calculated values and monitored values in the middle of gully (x=600 m)

图4  主沟出口处(x=1 200 m)泥位计算值与监测值对比

Fig. 4  Comparison of sludge level between calculated values and monitored values at the exit of gully (x=1 200 m)

表1  泥位计算相对误差

Table 1  Calculation relative error of sludge level

由图3和图4可知：2位置处的泥位计算值变化趋势基本符合实际监测结果，即先迅速增大而后趋于平缓。在初始时刻(t=0 min)和终止时刻(t=120 min)时，泥位计算值与监测值差别较大。一方面，由于前期已有持续降雨，因此在计算初始时刻已有少量地表积水；另一方面，泥石流爆发当天的实际降雨过程中，在14:40左右雨量已明显减小，该时刻之后监测值迅速减小。由表1可见：2处的平均相对误差分别为9.62%和7.96%(平均相对误差为不考虑0时刻的相对误差平均值)，计算结果基本符合实际情况。

3  结论

1) 将模型实验获得的渗流强度加入到圣维南方程组中，改变了质量守恒方程中的质量流通量，并且以经验参数形式加入侧流强度，建立一维黏性、沟谷型、暴雨型泥石流运动模型。

2) 将碱坪沟各参数条件代入模型中进行计算，并与监测资料进行对比验证，计算结果与实际情况较符合。该模型可以较好地描述碱坪沟泥石流运动过程中各参数的变化情况，为今后该地区泥石流预测与防治提供一定理论依据。

3) 在该类模型中，需进一步研究以解决泥石流的侵蚀挟带机理、岩土体的非均质性和地形的不规则性等问题。

参考文献：

[1] 褚胜名, 余斌, 李丽, 等. “8.13”碱坪沟泥石流形成机理及特征研究[J]. 中国水土保持, 2011, 11(8): 52-54.

ZHU Shengming, YU Bin, LI Li, et al. Forming mechanism and characteristics of debris flow happened on August 13, 2010 in Jianping gully[J]. Soil and Water Conservation in China, 2011, 11(8): 52-54.

[2] Johnson A M. 地质学中的物理过程[M]. 张之立, 李兴才, 译. 北京: 科学出版社, 1983: 54-63.

Johnson A M. The physical process in geology[M]. ZHANG Zhili, LI Xingcai, Trans. Beijing: Science Press, 1983: 54-63.

[3] 范椿. 一维定常泥石流的数学模型[J]. 力学与实践. 1994, 16(1): 50-52.

FAN Chun. One-dimensional mathematical model of debris flow[J]. Mechanics and Practice, 1994, 16(1): 50-52.

[4] Lin D G, Hsu S Y, Chang K T. Numerical simulation of flow motion and deposition characteristics of granular debris flows[J]. Natural Hazards, 2009, 50(3): 623-650.

[5] Takahashi T. Debris flow mechanics, prediction and countermeasures[M]. London: Taylor & Francis Group, 2007: 126-184.

[6] Takahashi T, Nakagawa H, Harada T, et al. Routing debris flows with particle segregation[J]. Journal of Hydraulic Engineering- ASCE, 1992, 118(11): 1490-1507.

[7] Iverson R M. The physics of debris flow[J]. Reviews of Geophysics, 1997, 35(3): 245-296.

[8] Iverson R M, Denlinger R P. Flow of variably fluidized granular masses, across three-dimensional terrain[J]. Journal of Geophysical Research, 2001, 106(B1): 537-552.

[9] Iverson R M, Logan M, Denlinger R P. Granular avalanches across irregular three-dimensional terrain, 2: Experimental tests[J]. Journal of Geophysical Research, 2004, 109(F1): 549-562.

[10] Zhang S C. A comprehensive approach to the observation and prevention of debris flows in China[J]. Nature Hazards, 1993, 3(7): 1-23.

[11] 李志威, 王兆印, 赵娜, 等. 弯曲河流斜槽裁弯模式与发育过程[J]. 水科学进展, 2013, 24(2): 161-168.

LI Zhiwei, WANG Zhaoyin, ZHAO Na, et al. Development processes of chute cutoffs on meandering rivers[J]. Advances in Water Science, 2013, 24(2): 161-168.

[12] 司伟, 包为民, 瞿思敏. 洪水预报产流误差的动态系统响应曲线修正方法[J]. 水科学进展, 2013, 24(4): 497-503.

SI Wei, BAO Weimin, QU Simin. Runoff error correction in real-time flood forecasting based on dynamic system response curve[J]. Advances in Water Science, 2013, 24(4): 497-503.

[13] 杨寅群, 李怀恩, 杨方社. 基于数学模型的陕西黑河水源区植被过滤带效果评估[J]. 水科学进展, 2013, 24(1): 42-48.

YANG Yanqun, LI Huaien, YANG Fangshe. An assessment of the effectiveness of vegetated filter strips for Heihe River headwaters area using numerical simulation[J]. Advances in Water Science, 2013, 24(1): 42-48.

[14] 周伟, 唐川, 周春花. 汶川震区暴雨泥石流激发雨量特征[J]. 水科学进展, 2012, 23(5): 650-655.

ZHOU Wei, TANG Chuan, ZHOU Chunhua. Critical rainfall characteristics for rainfall-induced debris flows in Wenchuan earthquake affected areas[J]. Advances in Water Science, 2013, 24(1): 42-48.

[15] 王光谦, 邵颂东, 费祥俊. 泥石流模拟I: 模型[J]. 泥沙研究, 1998(3): 7-13.

WANG Guangqian, SHAO Songdong, FEI Xiangjun. The simulation of debris flows: Model[J]. Journal of Sediment Research, 1998(3): 7-13.

[16] 杨国录. 四点时空偏心Preissmann格式的应用问题[J].泥沙研究, 1991(4): 88-98.

YANG Guolu. The application research of Preissmann four-point linear implicit scheme[J]. Journal of Sediment Research, 1991(4): 88-98.

[17] 张国华, 张展羽, 左长清, 等. 坡地自然降雨入渗产流的数值模拟[J]. 水利学报, 2007, 38(6): 668-673.

ZHANG Guohua, ZHANG Zhanyu, ZUO Changqing, et al. Numerical simulation of runoff generation on slope under the condition of rainfall[J]. Journal of Hydraulic Engineering, 2007, 38(6): 668-673.

(编辑  赵俊)

收稿日期：2014-08-06；修回日期：2014-11-10

基金项目(Foundation item)：国家科技重大专项(2010ZX03006-007)；国家重点基础研究发展计划(973计划)项目(2013CB035401)(Project (2010ZX03006-007) supported by the National Science and Technology Major Program of China; Project (2013CB035401) supported by the National Basic Research Program of China)

通信作者：姚磊华，博士，教授，从事地质灾害防治理论与方法研究；E-mail：yaolh@bjtu.edu.cn