泥石流数值模拟方法研究进展

摘要：泥石流作为一种多相混合介质，所包含的物理过程和动力学特征非常复杂。针对泥石流问题的数值模拟方法随着数值计算方法和物理计算模型的发展而发展，基于物理过程的数值模拟方法为探究泥石流复杂物理现象背后的机理提供了一种有效手段。回顾了求解泥石流动力问题的数值模拟方法，从连续介质计算方法、离散介质计算方法和混合介质计算方法3个方面分析了不同数值模拟方法的特点和适用情况，介绍了在泥石流分析中常用的数值模拟软件及其特点，展望了求解泥石流动力学问题的数值模拟方法的发展趋势。结果表明：传统的基于网格的计算方法已有长足发展和较长的应用历史，相对比较成熟，但是在处理大变形、快速运移的自由表面流问题时，存在网格容易畸变等问题；基于粒子的计算方法在处理上述问题时无需网格的划分和维护，易于确定自由表面位置和多相间的界面，但存在边界条件处理困难等问题；混合介质计算方法在较小尺度范围内对固体颗粒物质与液相相互作用机理进行探讨时具有重要作用。

关键词：泥石流；数值模拟；动力学模型；多相；连续介质；离散介质；混合介质；深度积分

中图分类号：P642.23文献标志码：A

0引言

泥石流是一种由土、砂、石等固体颗粒物与水组成的，在重力驱动下沿山坡或沟谷运移的混合流体，具有宽级配、高浓度、直进性、大冲大淤、冲击力大、破坏力强等特点，往往由暴雨、融雪、溃坝、滑坡等引发。成因、物源条件等因素的差异导致不同类型泥石流的动力学特征存在很大差异；沿程侵蚀、堆积的发生导致同一场泥石流在不同阶段也存在较大差异。

泥石流动力学特征受多种因素影响，如固体颗粒物体积分数、粒径组成、黏性物质含量等。根据体积分数可以将狭义泥石流分为3类，即黏性泥石流、稀性泥石流和过渡性泥石流，而广义泥石流还包括泥流和水石流。基于内部应力的特点，可将泥石流分为准静力泥石流和动力泥石流[1]。准静力泥石流是指体积分数大于50%（重度为1 830 kg·m-3）的泥石流，相当于黏性泥石流；动力泥石流是指体积分数在20%～50%之间（重度为1 330～1 830 kg·m-3）的泥石流，相当于稀性泥石流。

侵蚀、堆积的发展和变化伴随着泥石流运移的全过程，进而引起泥石流物质构成比例的变化，并进一步导致内部应力构成上的变化。不同类型的应力在泥石流内部具有一种此消彼长的关系。当颗粒碰撞应力占优时，黏性应力就会减小，如水石流，反之亦然。泥石流动力学数学模型应基于泥石流的物质构成和内部应力构成特点来建立。由于模型的复杂性，数学模型一般需要依靠各种数值方法进行求解。

泥石流动力学问题的数值模拟方法总体上可以分为连续介质、离散介质和混合介质3类。本文分析了采用数值模拟方法求解泥石流动力学问题的特殊性和难点，在此基础上进一步分析了上述3类数值方法在求解泥石流动力学问题方面的适用性及相关研究的最新进展，为系统、快速了解泥石流数值模拟研究提供参考。

1连续介质计算方法

连续介质计算方法以被研究对象的质量和变形连续分布为基本特征，建立描述动力学特征的方程组，并将这些方程基于网格或（和）粒子在计算域上进行离散化，结合初始和边界条件来求解。

1.1基于网格的计算方法

传统的基于网格的计算方法有：有限差分法（FDM）、有限元法（FEM）[2]和有限体积法（FVM）[34]。这类方法将计算域划分为很多连续分布的基本单元（网格或体元），在这些离散化的计算域上，控制方程基于某种数值方法被转化为包含节点未知场变量信息的代数方程组。控制方程有微分和积分两种形式，有限差分法求解微分形式的控制方程，在网格离散的基础上采用某种格式的代数差分代替偏微分控制方程中的偏导数，而有限元法和有限体积法则是基于离散控制体求解积分形式的动力学方程。

基于网格的计算方法在模拟泥石流问题时的主要挑战是：需要处理自由面和锋面移动、侵蚀锋面处物理量的稳定性、激波捕捉、不规则底床和沟道复杂网格的生成和维护等。在泥石流高速运移或岩土体崩滑问题中，急流汇入缓流时会有激波产生。将拉格朗日移动网格格式应用于非守恒形式的动力学方程可以较好地处理自由边界问题，但是无法处理激波。一些改进格式由此诞生，如修正的Godunov算法[5]、PatrovGalerkin法[6]、无振荡算法[7]、近似Riemann求解器[8]、全变差减小格式（TVD）[9]等。上述激波处理方法主要基于有限差分的守恒非振荡形式[10]或者有限体积法对流深守恒的属性[11]。这些算法明显改进了基于网格的计算方法，但这些算法应用于有限体积法时不能保证流深始终为正[12]，为此MangeneyCastelnau基于在粒子微观尺度引入的狄拉克分布函数提出了一种有限体积法的动力格式[12]。这种格式能有效处理不连续问题，再现颗粒崩滑运移中的流动和休止过程，并能保证流深始终为正。任意拉格朗日欧拉法（ALE）使有限差分法和有限元法能够更好地处理运动自由面的位置确定问题，同时该方法与基于固定欧拉网格的处理自由面方法（如流体体积（VOF）法[13]）相比，自由面的离散化精度更高[14]。传统的基于网格的计算方法求解泥石流运移问题的主要困难在于网格的移动和维护以及由网格引起的数值收敛问题。

Savage等研究了有限质量颗粒物质沿粗糙斜面崩滑问题，假设颗粒物质满足库仑摩擦定律，获得了描述崩滑的深度平均方程，即著名的SavageHutter模型[15]。Savage等采用欧拉和拉格朗日两种有限差分形式对模型进行求解，发现拉格朗日型差分格式能更好地处理颗粒物和空气的交界面，求解更简洁有效，结果更稳定可靠[15]。Wieland等采用有限体积法与有限差分法联合的方法求解二维SavageHutter模型[16]。Pitman等在SavageHutter模型中考虑了底床侵蚀的影响，并采用Godunov型有限体积法进行求解，提高了数值稳定性和计算精度[17]。Hungr采用拉格朗日有限差分法求解了滑坡动态模拟（DAN）模型[18]。该模型将运移物质看作等效流体，可以使用不同的流变模型并可近似考虑流变参数受底床物质加入、孔隙压力等因素的影响，近似考虑了运移过程中沟道的侧限作用。Iverson将SavageHutter模型的单相颗粒流模型拓展为二维颗粒流体两相混合模型[19]，又进一步拓展至三维情形，并可以考虑运移过程中孔隙压力的演变[20]，在计算过程中间无需调整模型初始参数便可描述泥石流从起动到堆积的全部动力过程。Denlinger将SavageHutter模型在三维不规则地形上进一步完善，显式地考虑了竖向重力加速度的影响并通过侧应力系数进一步完善了三维空间中库仑应力的表示，采用基于Roe型黎曼求解器的高精度有限体积法求解模型[21]。之后，Iverson等显式地考虑了剪胀和孔隙压力变化之间的相互影响，使模型所包含的方程数增加到5个，分别描述流深、固体体积比、底床孔隙压力和2个速度分量，每个方程包含一个体现颗粒剪胀影响的源项[22]，并采用Godunov型激波捕捉格式的有限体积法进行求解[3]。由于模型假设固液两相具有相同的速度，所以该模型实际是准两相流模型。

相关热词搜索： 泥石流 研究进展 数值 模拟 方法