CN106529198B - 一种泥石流全过程数值模拟及数值计算方法 - Google Patents

Abstract

本发明涉及一种泥石流全过程数值模拟及数值计算方法，其通过研究泥石流形成过程及机理、泥石流沟道物质启动、泥石流运动堆积过程数值计算方法，同时构建泥石流动力学过程和工程调控数值模拟平台，提高泥石流动力过程研究的科学性和可靠性，为泥石流减灾提供技术支撑，更好的服务泥石流减灾；提出的泥石流运动模型，将充分考虑泥石流浆体和大石块的动力特性；首次引入泥石流和动床的物质交换率的概念，考虑泥石流沿程侵蚀与灾害规模变化，形成基于全过程的泥石流形成‑运动‑堆积数值模拟及数值计算方法。本发明能提高泥石流动力过程研究的科学性和可靠性，为泥石流减灾提供技术支撑，更好的服务泥石流减灾。

Description

一种泥石流全过程数值模拟及数值计算方法

技术领域

本发明涉及环境治理技术领域，具体涉及一种泥石流全过程数值模拟及数值计算方法。

背景技术

泥石流是广泛分布于山区的自然灾害，暴发突然，运动速度快，具有巨大的冲击力，所携带的大量泥沙和石块，不仅可以在瞬间淤埋农田、森林、村庄和城镇，冲毁桥墩、房屋、路基，而且会直接进入河流水库，造成堰塞湖或者抬高河床，造成更大的次生灾害，严重威胁水利水电工程，在下游造成溃决洪水灾害，破坏山区的生态环境。泥石流灾害已经成为当前防灾减灾工作中亟待解决的突出问题，尤其是面临着大型泥石流灾害风险的严峻挑战。通过构建泥石流动力学模拟和工程调控平台，主动有效地预防和减轻泥石流灾害。

目前，泥石流频发区水汽特征和复杂地形对山地暴雨的影响，泥石流形成、汇流、放大效应机制，泥石流运动机理等涉及泥石流减灾的关键性问题认识不清楚，泥石流“形成-汇流-堆积”各个过程的物理意义不明确，未能从全过程考虑泥石流动力学机理研究，造成泥石流工程设计参数(流速、流量、冲击力等)计算误差大，防治工程设计主要依赖经验，导致泥石流预报不准、防治工程合理性和安全性差异较大。

以下为现有的三种泥石流机理研究：

一、泥石流“形成-运动-堆积”各个单一过程机理研究

泥石流作为一种由水、泥砂和石块组成的特殊水体，其内部水含量和松散物质之间的相互作用是泥石流形成的核心。当前，泥石流形成方面的研究主要从泥石流发生的外部环境条件进行分析，包括地貌条件、松散颗粒源、水源(降雨)、流域面积及形状、植被覆盖率、坡面或河床坡比以及其他因素等方面的综合研究。泥石流作为固液两相流，其内部作用力包括固相与液相的相互作用，以及固相内部不同颗粒的相互作用，固液两相的作用表现为流体中掺入固体颗粒以后，对流体粘性及其运动的影响；反过来，固体颗粒在高粘性的液体中运动，其流动参数又受到液相的影响。泥石流运动堆积的研究涉及泥石流的形态、流态、流量、流速、冲击、淤积等内容。部分研究从能量角度解释泥石流的运动堆积机理，也有大量学者从泥石流运动堆积模型实验，结合野外调查数据，建立泥石流方量、堆积参数的统计关系模型，尚未发现从泥石流“形成-运动-堆积”全过程的机理研究。

由于天然泥石流形成、运动、堆积各个过程的复杂多变性，使得难以建立通用型泥石流理论模型，从而只能从观测资料、室内试验以及历史资料等出发建立一些地区性经验公式。这些公式适用范围大都局限于某一地区或某一类型的泥石流，地区性较强。目前，泥石流形成、运动、堆积各个过程的物理意义不明确，尚未发现考虑全过程的泥石流机理研究。

二、泥石流“形成-运动-堆积”单一过程数值模拟

泥石流的动力学特征受多种因素影响，如固体颗粒物体积分数、粒径组成、粘性物质含量等，与此同时泥石流运动推移过程中伴随着侵蚀、堆积等物理过程，随机改变泥石流物质构成比例并进一步导致其内部应力构成上的异变。不同类型的应力在泥石流内部具有一种此消彼长的关系。当颗粒碰撞应力占优时,黏性应力就会减小，如水石流，反之亦然。泥石流动力学数学模型的建立是在泥石流的物质构成和内部应力构成特点基础上，由于模型的复杂性，数学模型一般需要依靠各种数值方法进行求解。通过相关学科软件结合相邻学科理论对泥石流的形成、运动和堆积进行模拟计算，旨在通过数字化形式再现泥石流的形成、运动和堆积等过程。

泥石流的数值模拟数学模型通常分为混合流(包括伪一相流和两相混合流)和颗粒流(非常适用于描述含有大量粗大颗粒的泥石流或者水石流运动)两种。目前，关于泥石流数值模拟研究，大多采用假定泥石流为单相流(宾汉体模型、膨胀体、混合流体)，对泥石流的形成运动堆积过程进行反分析，还原已发生灾害事件的方法进行处理，这些方法往往不能反应实际泥石流运动过程中固液相间的作用机理。此外，考虑两相流理论建立的泥石流数值模拟方法缺乏对沟道的冲淤及其与泥石流运动过程之间的关系的描述,不能够准确反映出泥石流的侵蚀和堆积机理,无法系统和定量模拟研究泥石流运动影响下沟道的冲淤动态过程。

三、泥石流数值模拟平台

国外已经初步形成了一批泥石流动力学模拟预测的软件系统产品，部分产品已经进入到商业化阶段，如Flo-2D，PFC2D、PFC3D、GeoDip、Debris-2D、DAN-3D等。这些软件被应用对泥石流进行数值模拟开展了大量研究工作，分析不同降雨工况下堆积体的运动速度场、应力场和位移场，研究崩滑松散碎屑堆积物质启动、形成泥石流的过程，分析降雨作用下松散碎屑物质启动形成泥石流的过程，再现了泥石流运动淤积过程，模拟分析得出了泥石流的威胁范围。

尽管已有多个数值模拟软件得到应用，但这些软件主要针对泥石流“形成-运行-堆积”中的单一个过程进行模拟，且大多数模型是基于恒定流的假设，对实际天然沟道中存在的阵性、流量、流速、含沙量等变化剧烈的非恒定泥石流运动演进过程研究较少,无法较为准确的描述实际泥石流运动状态和泥石流流变性质耦合响应关系的动力学机理。

尽管泥石流研究已经取得了显著进展，由于泥石流灾害形成和运动过程的复杂性和不确定性，目前尚无基于泥石流全过程的数值模拟技术与平台，此外，当前模拟软件对复杂地形进行了简化，采用的泥石流动力学机理相对较简单，未能反映泥石流复杂动力学机理。为此，针对泥石流形成运动堆积过程中的不足，就是本项发明的动因和突破重点。

发明内容

针对上述提到的问题，本发明提出了一种泥石流全过程数值模拟及数值计算平台，其通过研究泥石流形成过程及机理、泥石流沟道物质启动、泥石流运动堆积过程数值计算方法，同时构建泥石流动力学过程数值模拟平台(是基于真实地形对泥石流“形成-运动-堆积”全过程的数值模拟)，提高泥石流动力过程研究的科学性和可靠性，为泥石流减灾提供技术支撑，更好的服务泥石流减灾。

本发明的技术方案如下：

上述的泥石流全过程数值模拟及数值计算方法，首先通过定点观测和野外试验获取源区产流入渗的原始数据、流域降雨数据和不同类型源地的产流数据，研究不同土体与降雨分布条件下坡面产流和坡体破坏的特征；在试验和研究的基础上运用系统动力学方法和随机过程理论研究流域系统特征，进行泥石流源地活动的时空特征分区，分析泥石流流域内坡谷汇流的动力过程，确定泥石流产汇流的特征参数；同时，利用地理信息空间分析技术实现泥石流产汇流的分布式水文模拟和计算，并通过原始观测和模拟实验数据验证，改进理论模型和计算结果；分析在山区洪水作用下沟道水流和松散固体物质的相互作用机制，深入研究泥石流运动的动力学特性，改进泥石流运动方程，揭示泥石流动力和流体性质沿程演化过程，预测给定频率降雨条件下的泥石流规模与危害范围。

所述泥石流全过程数值模拟及数值计算方法，其中：所述泥石流全过程数值模拟及数值计算方法具体包括泥石流形成过程数值模拟、泥石流运动过程数值模拟和泥石流侵蚀与沉积动力过程数值仿真；所述泥石流形成过程数值模拟的具体流程包括形成阶段基础数据准备、构建泥石流形成过程的动力学模型、泥石流形成过程的数值模拟以及泥石流形成阶段输出结果。

所述泥石流全过程数值模拟及数值计算方法，其中，所述形成阶段基础数据准备具体包括：(111)计算域内的降雨历时曲线，即对于每个离散地形点的每一离散时刻均有对应的降雨强度，用坐标表达为(x，y，t，r)；其中x为地形点的x坐标，y为地形点的y坐标，t为时间，r为降雨强度；(112)、计算域内5m精度的地形数据，即对于每个离散地形点均有对应的高程值，用坐标表达为(x，y，z)；其中，x为地形点的x坐标，y为地形点的y坐标，z为地形点的z坐标；(113)计算域内土体的物理及力学参数即包括土体的中值粒径d₅₀、泥沙颗粒比重s、初始含水量θ_i、孔隙率n、水力传导度K、粘聚力c、内摩擦角和土壤基质吸力ψ。

所述泥石流全过程数值模拟及数值计算方法，其中，所述构建泥石流形成过程的动力学模型具体包括：

(121)构建水动力模型

式中，H为水深，P为x方向单宽流量，Q为y方向单宽流量，R为降雨，I为渗流，η为自由水面高程，τ_bx为x方向底部切应力，τ_by为y方向底部切应力，g为重力加速度，ρ为水沙平均密度，初始时刻为清水密度；

其中，τ_bx和τ_by采用如下公式计算：

n为曼宁粗糙系数；

(122)构建泥沙模型

上式(4)中，q_b为单宽推移质体积输沙率，s为泥沙颗粒比重，d₅₀为泥沙中值粒径，θ为Shields数，θ_cr为临界Shields数；

θ的计算式为：

上式(5)中，ρ_s为泥沙颗粒密度，ρ为水沙平均密度，初始时刻为清水密度；

θ_cr的计算式为：

上式(6)中，D_*为泥沙颗粒尺寸参数：

上式(7)中，ν_m为水的运动学粘滞系数；

(123)构建地形演变模型

上式(8)中，Z_b为地面高程的相反数，即有η＝H-Z_b；po_ro为孔隙率，

(124)构建下渗模型

下渗模型采用Green-Ampt下渗模型：

f＝dF/dt＝K(ψ△θ+F)/F (9)

F(t)＝Kt+ψ△θln[(1+F(t))/ψ△θ] (10)

上式(9)和(10)中，F为累积下渗量，f为下渗速度，ψ为土壤基质吸力，△θ＝η-θ_i为湿润锋上下含水量之差，即孔隙率与初始含水量之差；上式(10)为隐函数，其采用不动点迭代或者牛顿迭代法进行计算，然后代入式(9)即求得下渗速度f；

(125)定义变化的水沙平均密度ρ_av及泥沙平均浓度C_av

泥沙与水将进行掺混，考虑一定时间以后，泥沙与水达到充分掺混并具有相同的运动速度，从而定义水沙平均密度ρ_av及泥沙平均浓度C_av如下：

所述泥石流全过程数值模拟及数值计算方法，其中，所述泥石流形成过程的数值模拟的具体流程为:采用基于蛙跳格式的显示有限差分法对上述偏微分方程组进行数值离散；计算网格采用交错网格方式定义，即对于所有标量，如H、η均定义在网格的中心点；与速度相关的矢量，如U、P定义在网格的右边界，而V、Q则定义在网格的上边界；对于综合应力项，如正应力T_xx、T_yy定义在网格中心，而切应力T_xy、T_yx则定义在网格的右上角；其中，变量的下标i,j分别代表X方向和Y方向的空间网格节点号，变量的上标n表示计算所在时间步，△x和△y表示空间步长，△t代表时间步长；具体离散后得到离散方程如下：

(131)水动力方程

连续性方程：

动量方程：

其中：

其中，对流项采用一阶迎风格式进行计算，以x方向为例：

(132)地形演变模型

对于上式(25)中推移质的空间导数项，利用一阶迎风差分在网格中心进行离散如下：

所述泥石流全过程数值模拟及数值计算方法，其中，所述泥石流形成阶段输出结果具体为：5m精度网格下各离散点处水深H，x方向流速U，y方向流速V和泥沙平均浓度C_av。

所述泥石流全过程数值模拟及数值计算方法，其中，所述泥石流运动过程数值模拟具体是动力学方程构建，具体包括：

(211)构建新的泥石流侵蚀率方程

基于泥石流在可侵蚀基底运动所必须遵循的动量间断条件，建立侵蚀率是基底应力之间的关系：

以上公式(32)中，E代表侵蚀速率，代表泥石流平均密度，z_b代表基底地形高度，和分别代表泥石流两个方向上的速度，τ_1b是运动流体总的基底应力，τ_2s是基底的抗力，τ_1b采用如下速度平方相关方程:

上述式(33)中，C_f是一个小于0.1的无量纲参数，ρ_s是固体颗粒材料密度.τ_2s假设遵循库伦摩擦定律，

上述式(34)中，c和φ₂是基底物质的粘聚力和内摩擦角，λ₂是量化基底液化程度的孔隙水压力；

把公式(33)和(34)插入公式(32),得到侵蚀率公式如下：

上述公式(35)中，当上层物质运动速度为0时，基底摩擦应力则满足库伦摩擦应力，提出基底摩擦应力则满足：

上述公式(36)和(37)中，C_f代表紊动系数，s代表固相和液相密度之比，g_z代表重力在z方向上的分量，代表基底物质的内摩擦角，δ代表基底摩擦角，C_f和λ₂为在一定范围能调整的参数，其他参数都是固定的。

(212)构建泥石流运动、侵蚀与沟道演化耦合方程

从Navier-Stokes方程出发，推导出泥石流与基底物质侵蚀所必须满足的质量和动量方程以及间断边界条件，即：

所述泥石流全过程数值模拟及数值计算方法，其中，所述泥石流侵蚀与沉积动力过程数值仿真中侵蚀和沉积效应的泥石流动力学方程如下：

其中：

其中，k代表土压力系数，S_0x和S_fx分别代表重力源项和基底阻力源项在x方向分量，S_0y和S_fy分别代表重力源项和基底阻力源项在y方向分量，ρ_s和ρ_w分别代表固相和液相密度，c代表固相体积分数，p代表基底孔隙率；E为侵蚀效应，D为沉积效应，即

侵蚀方程：

沉积方程：D＝ω(1-C_a)^mC_a (48)

上述的侵蚀方程和沉积方程中涉及的参数是通过实验确定；其中，θ是Shields参数，θ_c是极限Shields参数，d是颗粒直径，是经验系数，ω代表颗粒沉降速度，m代表经验系数，C_a代表固相颗粒在垂向方向上分布规律系数。

有益效果：

因此，在泥石流形成过程及机理、泥石流汇流过程和运动堆积数值模拟方面的定量研究的突破，将为认识泥石流形成过程、定量刻画泥石流运动特征，提高泥石流预报精度，为准确分析泥石流动力响应机制，进而科学布设泥石流防治工程提供理论基础。而这些理论和方法的系统研究成果不仅提升我国泥石流研究的理论水平，而且能更好的服务于国家减灾，保障山区人民生命财产的安全，以及道路、水电和厂矿等重大工程建设的安全。

通过研究泥石流形成过程及机理、泥石流沟道物质启动、泥石流运动堆积过程数值计算方法，同时构建泥石流动力学过程和工程调控数值模拟平台，提高泥石流动力过程研究的科学性和可靠性，为泥石流减灾提供技术支撑，更好的服务泥石流减灾。计算平台提出的泥石流运动模型，将充分考虑泥石流浆体和大石块的动力特性；首次引入泥石流和动床的物质交换率的概念，考虑泥石流沿程侵蚀与灾害规模变化，形成基于全过程的泥石流形成-运动-堆积数值模拟技术，并开发相应的计算平台。与商业计算模拟软件的最大区别在于，这是国际上首个基于泥石流全过程动力学参数确定的模拟系统，可以实现事件之前的计算模拟和分析。

附图说明

图1为本发明泥石流全过程数值模拟及数值计算方法中泥石流侵蚀动力过程计算模拟示意图；

图2为本发明泥石流全过程数值模拟及数值计算方法中地震区某沟泥石流数值模拟结果图；

图3为本发明泥石流全过程数值模拟及数值计算方法中崩塌(碎屑流)运动模拟图；

图4为本发明泥石流全过程数值模拟及数值计算方法中泥石流形成数值模拟结果图；

图5为本发明泥石流全过程数值模拟及数值计算方法中泥石流泥石流运动堆积过程数值模拟技术路线图；

图6为本发明泥石流全过程数值模拟及数值计算方法中泥石流运动数值模拟图；

图7为本发明泥石流全过程数值模拟及数值计算方法中泥石流运动堆积数值模拟结果图。

具体实施方式

本发明泥石流全过程数值模拟及数值计算方法，是采用泥石流形成区野外定点观测、野外及室内模拟实验、理论分析、流域水文模拟、数值模拟、GIS和RS分析相结合的方法。首先通过定点观测和野外试验获取源区产流入渗的原始数据、流域降雨数据和不同类型源地的产流数据，研究不同土体与降雨分布条件下坡面产流和坡体破坏的特征。在这些试验和研究的基础上，运用系统动力学方法和随机过程理论研究流域系统特征，进行泥石流源地活动的时空特征分区，分析泥石流流域内坡谷汇流的动力过程，确定泥石流产汇流的特征参数。同时，利用GIS技术(地理信息空间分析技术)实现泥石流产汇流的分布式水文模拟和计算，并通过原始观测和模拟实验数据验证，改进理论模型和计算结果。进而，分析在山区洪水作用下沟道水流和松散固体物质的相互作用机制，深入研究泥石流运动的动力学特性，改进泥石流运动方程，揭示泥石流动力和流体性质沿程演化过程，预测给定频率降雨条件下的泥石流规模与危害范围，形成基于全过程的泥石流形成-运动-堆积数值模拟技术与计算平台。

本发明泥石流全过程数值模拟及数值计算方法，其具体包括：

S0100、泥石流形成过程数值模拟

S0110、形成阶段基础数据准备

S0111、计算域内的降雨历时曲线，即对于每个离散地形点的每一离散时刻均有对应的降雨强度，用坐标表达为(x，y，t，r)；其中x为地形点的x坐标，y为地形点的y坐标，t为时间，r为降雨强度；

S0112、计算域内5m精度的地形数据，即对于每个离散地形点均有对应的高程值，用坐标表达为(x，y，z)；其中，x为地形点的x坐标，y为地形点的y坐标，z为地形点的z坐标。

S0113、计算域内土体的物理及力学参数，如土体的中值粒径d₅₀，泥沙颗粒比重s，初始含水量θ_i，孔隙率n，水力传导度K，粘聚力c和内摩擦角土壤基质吸力ψ等。

S0120、泥石流形成过程的动力学模型

S0121、水动力模型

式中，H为水深，P为x方向单宽流量，Q为y方向单宽流量，R为降雨，I为渗流，η为自由水面高程，τ_bx为x方向底部切应力，τ_by为y方向底部切应力，g为重力加速度，ρ为水沙平均密度，初始时刻为清水密度。

其中，τ_bx和τ_by采用如下公式计算：

随着程序的计算进行，水沙平均密度ρ会不断改变，从而上述公式可以反映出泥石流形成的不同阶段所对应的底部切应力公式。其中，n为曼宁粗糙系数。

S0122、泥沙模型

其中，q_b为单宽推移质体积输沙率，s为泥沙颗粒比重，d₅₀为泥沙中值粒径，θ为Shields数，θ_cr为临界Shields数。

Shields数θ的计算式为：

其中，ρ_s为泥沙颗粒密度，ρ为水沙平均密度，初始时刻为清水密度。

临界Shields数θ_cr的计算式为：

上式中，D_*为泥沙颗粒尺寸参数：

其中，ν_m为水的运动学粘滞系数。

S0123、地形演变模型

其中，Z_b为地面高程的相反数，即有η＝H-Z_b；poro为孔隙率，

S0124、下渗模型

下渗模型采用Green-Ampt下渗模型：

f＝dF/dt＝K(ψ△θ+F)/F (9)

F(t)＝Kt+ψ△θln[(1+F(t))/ψ△θ] (10)

其中，F为累积下渗量，f为下渗速度，ψ为土壤基质吸力，△θ＝η-θ_i为湿润锋上下含水量之差，即孔隙率与初始含水量之差。

上式(10)式为隐函数，可采用不动点迭代或者牛顿迭代法进行计算，然后代入式(9)即可求得下渗速度f。

S0125、变化的水沙平均密度ρ_av及泥沙平均浓度C_av定义

随着水流的产生与汇集，冲刷也在不断地进行。泥沙与水将进行掺混，考虑一定时间以后，泥沙与水达到充分掺混，二者具有相同的运动速度，从而定义水沙平均密度ρ_av及泥沙平均浓度C_av如下：

S0130、泥石流形成过程的数值模拟

采用基于蛙跳格式的显示有限差分法对上述偏微分方程组进行数值离散。计算网格采用交错网格方式定义。即对于所有标量，如H、η等，均定义在网格的中心点；与速度相关的矢量，如U、P定义在网格的右边界，而V、Q则定义在网格的上边界。对于综合应力项，如正应力T_xx、T_yy定义在网格中心，而切应力T_xy、T_yx则定义在网格的右上角。其中，变量的下标i,j分别代表X方向和Y方向的空间网格节点号，变量的上标n表示计算所在时间步；△x和△y表示空间步长；△t代表时间步长。具体离散后得到离散方程如下：

S0131、水动力方程

连续性方程：

动量方程：

其中：

其中，对流项采用一阶迎风格式进行计算，以x方向为例：

S0132、地形演变模型

对于上式中推移质的空间导数项，利用一阶迎风差分在网格中心进行离散如下：

按顺序求解上述相应的离散后的代数方程组，具体顺序如下：

第一、求解变密度连续性方程组，更新水深H；

第二、求解动量方程组，更新单宽流量P，Q；

第三、求解泥沙模型，更新地面高程Z_b，从而更新水面高程η；

S0140、泥石流形成阶段输出结果

输出结果为5m精度网格下，各离散点处，水深(泥深)H，x方向流速U，y方向流速V，泥沙平均浓度C_av。

S0200、泥石流运动过程数值模拟

S0210、动力学方程构建

S0211、构建了新的泥石流侵蚀率方程

基于泥石流在可侵蚀基底运动所必须遵循的动量间断条件，建立侵蚀率是基底应力之间的关系

以上公式(32)中，E代表侵蚀速率，代表泥石流平均密度，z_b代表基底地形高度，和分别代表泥石流两个方向上的速度，τ_1b是运动流体总的基底应力，τ_2s是基底的抗力，τ_1b采用如下速度平方相关方程：

其中C_f是一个小于0.1的无量纲参数，ρ_s是固体颗粒材料密度.τ_2s假设遵循库伦摩擦定律，

其中c和φ₂是基底物质的粘聚力和内摩擦角；λ₂是量化基底液化程度的孔隙水压力。

把公式(33)和(34)插入公式(32),得到侵蚀率公式如下：

上述公式(35)与Fraccarollo and Capart(2002)类似并且已经在溃坝动力学计算模拟中得到了应用。τ_1b仅仅表示成与速度成平方关系是不合理的，这意味着如果上层物质运动速度为0时，基底摩擦应力也相应的为0，然而这是不真实的。当上层物质运动速度为0时，基底摩擦应力应该满足库伦摩擦应力，提出基底摩擦应力应该满足：

进一步

上述公式(36)和(37)中，C_f代表紊动系数，s代表固相和液相密度之比，g_z代表重力在z方向上的分量，代表基底物质的内摩擦角，δ代表基底摩擦角，C_f和λ₂为在一定范围能调整的参数，其他参数都是固定的。

推出的侵蚀率公式满足边界间断条件，并且修正了原公式只能应用于快速运动流体的限制，具有一定的意义。通过计算与泥石流实验和现场调查结果相比，发现当前的基底侵蚀模型可以很好的描述泥石流侵蚀动力学过程。

S0212、构建了泥石流运动、侵蚀与沟道演化耦合方程

从Navier-Stokes方程出发，推导出了泥石流与基底物质侵蚀所必须满足的质量和动量方程以及间断边界条件。

泥石流动量守恒公式(12)-(16)右端清晰的显示了泥石流与基底物质交换造成的动量交换。当泥石流与沉积物层之间的基底剪切应力采用流体底部基底阻力时，动量交换项必须包含在方程右端，而泥石流与沉积物层之间的基底剪切应力采用沉积层内剪切抗力时，动量交换项隐含在沉积层内剪切抗力内。此前，有关流体与基底物质交换造成的动量交换相互作用关系长期被人们忽视，造成各种不同公式表达式存在，给研究工作带来困扰。而通过严格的深度平均公式推导，能够给出符合物理规律的合理公式表达式，给数值计算仿真准确性带来保障。最终实现了考虑侵蚀和沉积的泥石流动力过程计算程序研发，并通过多种途径验证了计算的可靠性。

S0300、泥石流侵蚀与沉积动力过程数值仿真

考虑侵蚀和沉积的泥石流动力过程计算程序研发，并通过多种途径验证了计算的可靠性。其中考虑侵蚀和沉积效应的泥石流动力学方程如下：

其中：

其中，k代表土压力系数，S_0x和S_fx分别代表重力源项和基底阻力源项在x方向分量，S_0y和S_fy分别代表重力源项和基底阻力源项在y方向分量，ρ_s和ρ_w分别代表固相和液相密度，c代表固相体积分数，p代表基底孔隙率；E为侵蚀效应，D为沉积效应，即：

侵蚀方程：

沉积方程：D＝ω(1-C_a)^mC_a (48)

上述的侵蚀方程和沉积方程中涉及的参数是通过实验确定；其中，θ是Shields参数，θ_c是极限Shields参数，d是颗粒直径，是经验系数，ω代表颗粒沉降速度，m代表经验系数，C_a代表固相颗粒在垂向方向上分布规律系数。

考虑侵蚀和沉积效应的泥石流动力学方程可以认为是一个带有源项的非线性双曲偏微分方程组。当前我们采用有限差分方法TVD-MacCormack求解。此方法的优点是在时间和空间上都同时保证具有二阶精度，并且程序实现起来方便可行。采用开发的泥石流动力过程计算程序，已经开展泥石流运动－堆积过程数值试验研究，并与典型泥石流调查相结合，反演泥石流的动力演化过程，通过修正计算模型及相关参数，再现泥石流的全过程，揭示泥石流的动力学机理与过程。(如图1)

下面结合具体实施例，对本发明作进一步详述：

实施例1

如图2，某泥石流沟位于都江堰龙池镇云华村，属于岷江一级支流龙溪河右岸的次级支沟，该沟是由主沟八一沟，支沟大干沟、小干沟和小湾沟形成的泥石流沟。流域总体自北向南伸展，汇水面积约8.63km²，主沟全长约4.45km，分水岭高程2456m，沟口高程850.7m，落差1605.3m，沟床平均纵坡降376.7‰。受5.12汶川地震的影响，岩体崩解坡积物滑落，产生滑坡和崩塌，进一步增加了沟道流域内的松散堆积物量。在2010年8·13降雨过程中，先后有两次大规模暴发泥石流。

利用本发明的对该沟进行全过程泥石流地理过程模拟，具体步骤如下：

步骤S1000、泥石流形成过程数值模拟

崩塌(碎屑流)运动计算

输入流域内的降雨历时曲线(即对于每个离散地形点的每一离散时刻均有对应的降雨强度)、5m精度的地形数据、域内土体的物理及力学参数(如泥沙颗粒比重，初始含水量，孔隙率等)。数据加载过程中系统自动匹配空间坐标系，并存储各个数据的图例区间、图例值和颜色，为泥石流形成模拟时实时查看泥石流相应的图例变化(图3)。

泥石流形成计算

结合本发明中泥石流形成模型，依据泥沙平均浓度可以判识泥石流(山洪)类型。同时，系统平台产汇流计算之后可以根据产汇流计算数据生成量的变化监控产汇流计算的进度，加载并实时展示产汇流结果，可以实现对泥石流形成过程的监测(图4)。

步骤S2000、泥石流运动堆积过程模拟

根据泥石流形成阶段的计算结果，确定泥石流性质，进一步确定泥石流运行模型(如稀性泥石流、粘性泥石流、洪水模型)。然后输入泥石流运动的相关动力学条件或参数，包括泥石流运动是否考虑密度变化，是否考虑基底侵蚀，设置相关基底侵蚀模型，设置相应的时间步长，数据输出标准等。(如图5)

启动泥石流运动计算之后，可以根据运动计算数据生成量的变化监控运动计算的进度，并按照特定格式存储相应的计算结果，通过加载泥石流运动结果，可以把泥石流运动计算生成的结果加载到系统平台，中途也可以暂停、继续或取消数据加载(图6)。数据加载过程中系统自动匹配空间坐标系，并存储各个数据的图例区间、图例值和颜色，实现展示泥石流运动模拟情景(图7)。

本发明实施例1具有以下优点：

(1)将为中国西南山区暴雨泥石流灾害防治提供技术支撑，为进行生态脆弱区泥石流治理、维护山区河流系统健康提供科技支撑。

(2)全过程的动力学泥石流数值模拟计算平台，可以提供全社会使用，提高我国泥石流防治的规范化和科学化水平，为科学合理设计泥石流减灾工程提供技术平台和具有可操作的通用工具。并通过位于汶川地震极重灾区泥石流减灾示范点，推广减灾技术，培训减灾人员，提升我国泥石流防治水平，在防泥石流灾害防治，保障人民生命财产和重大工程安全方面做出中国科学院具有显示度的贡献。

本发明通过研究泥石流形成过程及机理、泥石流沟道物质启动、泥石流运动堆积过程数值计算方法，同时构建泥石流动力学过程和工程调控数值模拟平台，提高泥石流动力过程研究的科学性和可靠性，为泥石流减灾提供技术支撑，更好的服务泥石流减灾。

Claims (7)

1.一种泥石流全过程数值模拟及数值计算方法，其特征在于：首先通过定点观测和野外试验获取源区产流入渗的原始数据、流域降雨数据和不同类型源地的产流数据，研究不同土体与降雨分布条件下坡面产流和坡体破坏的特征；在试验和研究的基础上运用系统动力学方法和随机过程理论研究流域系统特征，进行泥石流源地活动的时空特征分区，分析泥石流流域内坡谷汇流的动力过程，确定泥石流产汇流的特征参数；同时，利用地理信息空间分析技术实现泥石流产汇流的分布式水文模拟和计算，并通过原始观测和模拟实验数据验证，改进理论模型和计算结果；分析在山区洪水作用下沟道水流和松散固体物质的相互作用机制，深入研究泥石流运动的动力学特性，改进泥石流运动方程，揭示泥石流动力和流体性质沿程演化过程，预测给定频率降雨条件下的泥石流规模与危害范围；

所述泥石流全过程数值模拟及数值计算方法具体包括泥石流形成过程数值模拟、泥石流运动过程数值模拟和泥石流侵蚀与沉积动力过程数值仿真；所述泥石流形成过程数值模拟的具体流程包括形成阶段基础数据准备、构建泥石流形成过程的动力学模型、泥石流形成过程的数值模拟以及泥石流形成阶段输出结果。

2.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述形成阶段基础数据准备具体包括：

(111)计算域内的降雨历时曲线，即对于每个离散地形点的每一离散时刻均有对应的降雨强度，用坐标表达为(x，y，t，r)；其中x为地形点的x坐标，y为地形点的y坐标，t为时间，r为降雨强度；

(112)、计算域内5m精度的地形数据，即对于每个离散地形点均有对应的高程值，用坐标表达为(x，y，z)；其中，x为地形点的x坐标，y为地形点的y坐标，z为地形点的z坐标；

(113)计算域内土体的物理及力学参数即包括土体的中值粒径d₅₀、泥沙颗粒比重s、初始含水量θ_i、孔隙率n、水力传导度K、粘聚力c、内摩擦角和土壤基质吸力ψ。

3.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述构建泥石流形成过程的动力学模型具体包括：

(121)构建水动力模型

式中，H为水深，P为x方向单宽流量，Q为y方向单宽流量，R为降雨，I为渗流，η为自由水面高程，τ_bx为x方向底部切应力，τ_by为y方向底部切应力，g为重力加速度，ρ为水沙平均密度，初始时刻为清水密度；

其中，τ_bx和τ_by采用如下公式计算：

n为曼宁粗糙系数；

(122)构建泥沙模型

上式(4)中，q_b为单宽推移质体积输沙率，s为泥沙颗粒比重，d₅₀为泥沙中值粒径，θ为Shields数，θ_cr为临界Shields数；

θ的计算式为：

上式(5)中，ρ_s为泥沙颗粒密度，ρ为水沙平均密度，初始时刻为清水密度；

θ_cr的计算式为：

上式(6)中，D_*为泥沙颗粒尺寸参数：

上式(7)中，ν_m为水的运动学粘滞系数；

(123)构建地形演变模型

上式(8)中，Z_b为地面高程的相反数，即有η＝H-Z_b；poro为孔隙率，

(124)构建下渗模型

下渗模型采用Green-Ampt下渗模型：

f＝dF/dt＝K(ψΔθ+F)/F (9)

F(t)＝Kt+ψΔθln[(1+F(t))/ψΔθ] (10)

上式(9)和(10)中，F为累积下渗量，f为下渗速度，ψ为土壤基质吸力，Δθ＝η-θ_i为湿润锋上下含水量之差，即孔隙率与初始含水量之差；上式(10)为隐函数，其采用不动点迭代或者牛顿迭代法进行计算，然后代入式(9)即求得下渗速度f；

(125)定义变化的水沙平均密度ρ_av及泥沙平均浓度C_av

泥沙与水将进行掺混，考虑一定时间以后，泥沙与水达到充分掺混并具有相同的运动速度，从而定义水沙平均密度ρ_av及泥沙平均浓度C_av如下：

4.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述泥石流形成过程的数值模拟的具体流程为:采用基于蛙跳格式的显示有限差分法对上述偏微分方程组进行数值离散；计算网格采用交错网格方式定义，即对于所有标量，如H、η均定义在网格的中心点；与速度相关的矢量，如U、P定义在网格的右边界，而V、Q则定义在网格的上边界；对于综合应力项，如正应力T_xx、T_yy定义在网格中心，而切应力T_xy、T_yx则定义在网格的右上角；其中，变量的下标i,j分别代表X方向和Y方向的空间网格节点号，变量的上标n表示计算所在时间步，Δx和Δy表示空间步长，Δt代表时间步长；具体离散后得到离散方程如下：

(131)水动力方程

连续性方程：

动量方程：

其中：

其中，对流项采用一阶迎风格式进行计算，以x方向为例：

(132)地形演变模型

对于上式(25)中推移质的空间导数项，利用一阶迎风差分在网格中心进行离散如下：

5.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述泥石流形成阶段输出结果具体为：5m精度网格下各离散点处水深H，x方向流速U，y方向流速V和泥沙平均浓度C_av。

6.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述泥石流运动过程数值模拟具体是动力学方程构建，具体包括：

(211)构建新的泥石流侵蚀率方程

基于泥石流在可侵蚀基底运动所必须遵循的动量间断条件，建立侵蚀率是基底应力之间的关系：

以上公式(32)中，E代表侵蚀速率，代表泥石流平均密度，z_b代表基底地形高度，和分别代表泥石流两个方向上的速度，τ_1b是运动流体总的基底应力，τ_2s是基底的抗力，τ_1b采用如下速度平方相关方程:

上述式(33)中，C_f是一个小于0.1的无量纲参数，ρ_s是固体颗粒材料密度.τ_2s假设遵循库伦摩擦定律，

上述式(34)中，c和φ₂是基底物质的粘聚力和内摩擦角，λ₂是量化基底液化程度的孔隙水压力；

把公式(33)和(34)插入公式(32),得到侵蚀率公式如下：

上述公式(35)中，当上层物质运动速度为0时，基底摩擦应力则满足库伦摩擦应力，提出基底摩擦应力则满足：

上述公式(36)和(37)中，C_f代表紊动系数，s代表固相和液相密度之比，g_z代表重力在z方向上的分量，代表基底物质的内摩擦角，δ代表基底摩擦角，C_f和λ₂为在一定范围能调整的参数，其他参数都是固定的；

(212)构建泥石流运动、侵蚀与沟道演化耦合方程

从Navier-Stokes方程出发，推导出泥石流与基底物质侵蚀所必须满足的质量和动量方程以及间断边界条件，即：

7.如权利要求1所述的泥石流全过程数值模拟及数值计算方法，其特征在于，所述泥石流侵蚀与沉积动力过程数值仿真中侵蚀和沉积效应的泥石流动力学方程如下：

其中：

其中，k代表土压力系数，S_0x和S_fx分别代表重力源项和基底阻力源项在x方向分量，S_0y和S_fy分别代表重力源项和基底阻力源项在y方向分量，ρ_s和ρ_w分别代表固相和液相密度，c代表固相体积分数，p代表基底孔隙率；E为侵蚀效应，D为沉积效应，即

侵蚀方程：

沉积方程：D＝ω(1-C_a)^mC_a (48)

上述的侵蚀方程和沉积方程中涉及的参数是通过实验确定；其中，θ是S hields参数，θ_c是极限Shields参数，d是颗粒直径，是经验系数，ω代表颗粒沉降速度，m代表经验系数，C_a代表固相颗粒在垂向方向上分布规律系数。