泥石流动力学模型与数值模拟

胡凯衡 崔鹏 李浦

①研究员，②中国科学院院士，③博士研究生，中国科学院山地灾害与地表过程重点实验室，中国科学院水利部成都山地灾害与环境研究所，成都 610041

泥石流动力学模型与数值模拟

胡凯衡①崔鹏②李浦③

①研究员，②中国科学院院士，③博士研究生，中国科学院山地灾害与地表过程重点实验室，中国科学院水利部成都山地灾害与环境研究所，成都 610041

泥石流；动力学；数值模拟；本构关系

泥石流是介于崩塌、滑坡等块体运动和山洪水流之间的一种物理过程，既具有土体的结构性，又具有水体的流动性。从物质组成上看，泥石流是一种由水、岩土体和气体组成的多相介质，具有多种内部结构。其中的固相组分颗粒形状极其不规则，尺度跨越范围较大。浆体与固相颗粒以及颗粒间的相互作用非常复杂。因此，完全考虑各种因素建立全描述的泥石流动力学模型比较困难。从描述组成物质和运动的观点来看，目前泥石流的动力学模型可划分为连续介质、离散介质和混合介质模型。基于不同动力学模型的泥石流运动数值模拟，可广泛应用于流量过程反演、危险范围预测、风险评估、防治工程评估等方面；但泥石流本身启动和产汇流机理涉及多门学科，是需要进一步研究的重要课题。同时，泥石流形成过程和运动过程的数值模拟存在着时空尺度上的差异，如何实现两者的耦合求解仍需深入探讨。

泥石流是由水、土、岩石等多相物质在山区坡地上或沟道内相互作用发展而形成的一种自然灾害，可以视为介于崩塌、滑坡等块体运动与山洪水流之间的一种物理过程。水多土(石)少，容易形成山洪；水少土(石)多，容易形成滑坡碎屑流。当水土达到一定比例，既可以在一定动力作用下受剪流动，又可以由内部的屈服强度抵抗一定的剪切力而停积，形成所谓的泥石流。因此，泥石流体既具有土体的结构性，又具有水体的流动性，但与这两者又有显著的差异。除了多相、碎散、流动和屈服等特征外，泥石流还具有多尺度的特征。组成泥石流的固体组成颗粒形状极其不规则，尺度从黏土、粉沙、细沙、砾石一直到直径几米的漂砾，跨越了数个量级。这些特征给泥石流运动的物理描述和动力学建模带来诸多的困难。现有的泥石流动力学模型可以描述泥石流的部分性质或者现象，但是，迄今为止还没有一个通用的、能够完全反映泥石流物理本质的理论模型。

1 泥石流的物理特征

泥石流体虽然也混有少量气体，但一般情况下将其视为固液两相混合体。它具有一定的屈服强度，能静止在相对平缓的坡度上。其主要特征是：充分饱和、密度大、静切力大、黏度高、固体颗粒组成不均匀等。泥石流物理特征主要指泥石流的固液组成、颗粒级配、内部结构、流变等特征。

1.1 物质组成

泥石流最本质的特征在于它是一种多相介质，主要由水、岩土体和气体组成。在实际沟道运动过程中，泥石流中水与松散岩土组成的混合体受到地形和自身结构的影响会上下翻腾，裹挟附近的空气使之成为泥石流体内部包裹的气体。尤其在沟道跌坎处，高速行进的泥石流往往会与周围空气发生剧烈的相互作用，水石飞溅，表层停留大量的气泡。整个过程水、岩土体和气体三相组成变化大，空间分布很不均匀。

实际上，与土体和水相比，泥石流体中气体含量较少，所以，一般将泥石流体的物质组成简化为水、土和砂石三种成分，其中土砂组分比例、砂石的级配构成因地域、地形条件各异。从物理性质上看，泥石流是一种碎散物质集合体，内部存在无数多个物质界面。泥石流中水的含量和赋存形式在很大程度上决定了泥石流的力学性质，不同性质的泥石流，其各种形式赋存水的含量有明显差异。

1.2 内部结构

由于物质组成和矿物成分不同，在岩土体与水的结合方式、物质的搬运方式的影响下，泥石流体内黏粒、粉粒、沙粒和石块等土体颗粒与含电解质水之间产生各种连接和排列形式，可以形成三种密切相关的结构：网格结构、网粒结构、格架结构。泥石流的物理力学特性不仅依赖于泥石流的物质组成，而且受泥石流结构(即不同物质在泥石流体中的空间排列和组合)的强烈影响。各类泥石流体均有不同程度的结构性。稀性泥石流的结构性体现在细粒浆体上，即由黏粒和物理性黏粒(粉粒)与含有电解质的水构成细粒浆体的网格结构(图1)；泥质泥石流或泥流的结构体现在粗粒浆体上，即沙粒与具有网格结构的细粒浆体结合构成粗粒浆体的网粒结构；泥石流体的结构则是石块与具有网粒结构的粗粒浆体结合构成格架结构[1](图2)。

1.3 物理描述

泥石流体中的水和固体物质并不是截然分明的，是泥浆和大小石块两部分组成的结构性流体。实际上泥石流的物质体系有三部分：①由固相中的细颗粒和液相中的水组成的泥浆体，一般认为细颗粒的粒径上限为2 mm，细颗粒的级配比例影响浆体性质；②粒径大于2 mm而又小于泥石流运动特征尺度的粗颗粒，它充填在泥浆液中对泥石流体的性质有显著影响；③泥石流中的漂砾，粒径大于泥石流运动特征尺度，是泥石流的搬运物。虽然泥石流的运动受其影响，但泥石流体的性质不受它影响。

2 泥石流的动力学模型

泥石流运动过程中，出现许多与洪水、挟沙水流不同的现象，比如阵性来流、物质分选、高速低阻等。Iverson[2]归纳了泥石流所特有的11种现象，认为衡量泥石流动力学模型的优劣取决于描述这些现象的能力。但是，泥石流运动时浆体与粗颗粒、粗颗粒与粗颗粒之间的相互作用非常复杂，完全考虑各种因素建立全描述的泥石流动力学模型较为困难。所以，现有的泥石流动力学模型都是在一定简化和假设的基础上建立的，只能描述泥石流的部分特有现象。

2.1 内部作用力

除了重力、床面阻力等外部作用力之外，泥石流内部还存在浆体黏滞力、孔隙压力、颗粒间摩擦力和碰撞力、紊动上举力等。液相的浆体切应力由三部分组成：宾汉应力、黏性应力、紊动应力。固相的颗粒切应力包括颗粒间静态支撑或低速流动下的准静止摩擦应力、中速流动下由颗粒间的相互滑动或挤压产生的作用力以及高速流动下颗粒间的碰撞应力。这些内力和外力共同作用而维持泥石流的运动。

Iverson等[3]提出了泥石流的库仑混合(Coulomb mixture)模型，考虑了颗粒物质的本构关系和两相之间的相互作用力，尤其是孔隙水压力的存在对于颗粒间相互作用力的影响。图3表示泥石流体中粗、细颗粒的受力情况，泥石流体中细颗粒(颗粒A)受到黏性力的作用而悬浮在流体中，成为液相的一部分；粗颗粒(颗粒B)只有依靠与其他颗粒间的相互作用力才能维持平衡。图4表示斜坡上在重力作用下稳定的颗粒流纵截面，高浓度流动区以颗粒间摩擦为主，低浓度流动区以颗粒间碰撞为主。

图3 泥石流体中粗、细颗粒的受力情况[2]

图4 斜坡上稳定的颗粒流纵截面

这些内部作用力之间的大小和量级随泥石流的物质组成不同而不同。比如：对于粗颗粒浓度较大的泥石流，颗粒间的摩擦力占主导地位；对于粗颗粒浓度较低的泥石流，颗粒间的碰撞力为主。Iverson[2]定义了Savage数、Bagnold数和Darcy数三个无量纲数来衡量不同作用力之间的量级。

2.2 本构关系

泥石流的应力本构关系是指泥石流的剪应力与应变、应变率或法应力之间的函数关系。泥石流是多相介质，内部存在各种力的作用，比如颗粒摩擦力、碰撞力、浆体黏聚力、紊动应力、黏滞力等。不同研究者往往侧重考虑不同类型的作用力，进而提出了各种描述泥石流体本构关系的模型。泥石流体本构模型主要分为单相连续介质模型和两相模型。单相介质本构模型即假定泥石流中的固体颗粒都是均匀分布于流体中，常用的有宾汉体(Bingham)、Bagnold的膨胀体和H-B(Herschel-Bulkley)流体等非牛顿体本构模型。鉴于单相连续介质模型存在的缺陷以及泥石流研究的不断深入，两相模型越来越受到重视。典型的两相本构模型有王光谦等[4]建立的泥石流两流体模型，即液相采用宾汉体模型，固相采用膨胀体模型，Iverson等[3]提出库仑颗粒流模型(表1)。此外，傅旭东等[5]建立的泥石流的动力学模型，通过分析低浓度固液两相流颗粒相速度的分布函数，也对泥石流颗粒相的本构关系进行了阐述。

表1 主要的泥石流本构模型

2.3 动力学模型

从描述组成物质和运动的观点来看，泥石流的动力学模型可划分为连续介质、离散介质和混合介质模型。连续介质模型假设泥石流体在空间连续而无空隙地分布，其宏观物理量如速度、密度等都是空间和时间的连续函数，满足质量守恒、动量和能量守恒定律。虽然实际的泥石流体内部是不连续的，但是，对于主要由细颗粒组成的泥流或者两相速度差很小的泥石流，连续介质在一定程度上可以描述泥石流的运动。离散介质模型是将泥石流简化为由大量的具有一定大小的物质颗粒组成的体系，颗粒之间遵循一定的碰撞或摩擦规律，其宏观物理量是颗粒质量和速度分布等系综平均的结果。混合介质模型用连续介质和离散介质模型分别来描述泥石流体的液相和固相部分的运动。实际中的泥石流都是由两相组成，只是对某些特殊类型的泥石流，比如泥流，可以简化为一相流。一般来说，泥石流液相是由小于一定尺寸细颗粒和水构成的泥浆液体，可视为连续介质；固相为一定尺寸以上的大颗粒，应该用离散介质模型来描述。

3 泥石流的数值模拟

泥石流的数值计算研究最初基于一维均质单相的连续介质模型。由于数值计算具有室内模型实验不可比拟的优势，泥石流运动的数值模拟近年来发展比较快，目前已经逐渐扩展到对多维、非均质、多相模型的数值求解，在流量过程反演、危险范围预测、风险评估、防治工程评估等方面得到广泛应用。

研究泥石流的不同物理过程涉及不同的作用力，以及不同的时间和空间尺度。比如：模拟泥石流翻越障碍物或撞击结构体涉及泥石流与固壁的撞击和摩擦作用，时间尺度为分钟，空间尺度为几十米；模拟泥石流在堆积扇上的扩散需要考虑底面的阻力、泥石流的堆积和侧向运动，时间尺度为小时，空间尺度为上千米。泥石流的运动方程通常为一组非线性的偏微分方程。其源项远比水流运动方程中的复杂。数值算法的选择要根据不同的动力学模型、问题的时间空间尺度、物理特征、收敛性和精度要求等来决定。常用的数值方法主要有有限元法(FEM)、有限体积法(FVM)和有限差分法(FDM)。比如，瑞士的DFEM泥石流计算软件采用的是有限元法，美国的FLO-2D软件采用的是有限差分算法。也有些研究者使用一些可处理离散介质和复杂界面的先进算法，比如格子Boltzmann法(LBM)、SPH(smoothed particle hydrodynamics)法。

连续介质模型中采用不同的泥石流本构关系对数值计算结果有很大的影响。Rickenmann等[6]综合比较了Voellmy、Bingham和Herschel-Bulkley三种常见的本构关系对数值模拟结果的影响。Naef等[7]用DFEM-1D算法比较了Bingham、Voellmy、湍流和Coulomb几种流变阻力关系不同组合对泥石流数值计算的影响，发现含有湍流项的阻力关系结果最好。地形数据的精度以及地形突变点的处理对数值计算结果的影响也是非常大的。DEM高程的1～2 m误差就可能导致泥石流的流向发生很大偏离。在地形的突变点处，比如沟道的跌坎、陡壁、陡坡和堆积扇的洼地处往往会出现流速或流深的奇异值。建筑物等障碍物的存在对计算结果也有直接的影响(图5)。此外，泥石流的运动边界在很多情况下是随时间而变化的。在计算中往往要采用Lagrange的方法来追踪泥石流的自由边界。比如，胡凯衡等[8]采用质点网格法来追踪泥石流在每一时刻的运动范围。等值面函数方法、流体体积法和波前追踪法等运动界面追踪方法，也在泥石流数值计算中有所应用。

图5 舟曲特大泥石流泛滥范围的模拟结果：(a) 不考虑建筑物影响；(b) 考虑建筑物影响

4 发展趋势和展望

现有的泥石流动力学模型和数值模拟涉及的主要是泥石流形成后在沟道和堆积扇运动的过程。浅水波、颗粒流和多相流的动力学模型和计算方法针对这一过程相对成熟，已经发展到可以实际应用的程度。在泥石流运动过程方面需要更深入的研究工作是，考虑泥石流体的非均质性以及动床的快速侵蚀作用。目前，相对来说还处于空白的是泥石流形成的动力学模型。暴雨型泥石流的形成过程，首先是在降雨的作用下坡面土体失稳启动，坡面和细沟产流，然后通过上游的支沟汇流，最终形成一定规模的泥石流。这一过程涉及水文学、土力学和流体力学等学科，是非常复杂的一个物理过程。建立一个能描述泥石流启动、产流和汇流全过程的动力学模型是未来泥石流学科的一个重要研究主题。除此之外，泥石流形成与运动的动力学模型在物理变量、时间空间尺度等方面差异非常大，如何将两者耦合在一起数值求解也是一项非常具有挑战性的课题。

(2014年8月24日收稿)■

参考文献

[1] 吴积善, 田连权, 康志成, 等. 泥石流及其综合治理[M]. 北京: 科学出版社, 1993: 29-32.

[2] IVERSON R M. The physics of debris flows [J]. Reviews of Geophysics, 1997, 35(3): 245-296.

[3] IVERSON R M, DENLINGER R P. Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory [J]. J Geophys Res, 2001, 106(B1): 537-552.

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

[5] 傅旭东, 王光谦. 低浓度固液两相流颗粒相本构关系的动理学分析[J]. 清华大学学报: 自然科学版, 2002, 42(4): 560-563.

[6] RICKENMANN D, LAIGLE D, MCARDELL B W, et al. Comparison of 2D debris-flow simulation models with field events [J]. Computational Geosciences, 2006, 10(2): 241-264.

[7] NAEF D, RICKENMANN D, RUTSCHMANN P, et al. Comparison of flow resistance relations for debris flows using a one-dimensional finite element simulation model [J]. Natural Hazards and Earth System Science, 2006, 6(1): 155-165.

[8] 胡凯衡, 韦方强, 何易平, 等. 流团模型在泥石流危险度分区中的应用[J]. 山地学报, 2004, 21(6): 726-730.

Debris flow dynamic models and numerical computation

HU Kai-heng①, CUI Peng②, LI Pu③
①Professor, ②CAS Member, ③Ph. D. Candidate, Key Laboratory of Mountain Hazards and Earth Surface Process, Chinese Academy of Sciences, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences & Ministry of Water Conservancy, Chengdu 610041, China

Debris flow as one of the most active surface processes is an intermediate form of mass movement between mountain flood and landslide. It is a viscous and highly concentrated fluid-solid mixture, being shear-resistant like soil while flowing like water. As a typical multiphase flow, most of debris flows contain water, soil, rocks and gas, and have special inner structures. The solid grains of debris flows have very irregular shape, and their sizes range several orders of magnitude. The interactions between slurry and grain or grain and grain are complex. Therefore, it is very difficult to develop a dynamic model fully describing the debris-flow process. From the point of view of constituent prototype and motion description, the existing debris-flow dynamic models can be classified into continuum, discrete medium, and hybrid medium models. The numerical simulations based on these dynamic models are widely applied to back-analysis of debris flow discharge, predicting inundated region, risk assessment, and evaluating protection engineering etc. Comparing with the motion process of debris flow, the formation process of debris flow is an interdisciplinary research subject, and needs more attentions. Moreover, how to coupling the two processes with different spatial and temporal scales is still a challenging issue.

debris flow, dynamic model, numerical simulation, constitutive relation

(编辑：温文)

10.3969/j.issn.0253-9608.2014.05.001