CN109657322B - 一种固液多相适用于泥石流的动力学数值模拟方法 - Google Patents

一种固液多相适用于泥石流的动力学数值模拟方法 Download PDF

Info

Publication number
CN109657322B
CN109657322B CN201811522183.8A CN201811522183A CN109657322B CN 109657322 B CN109657322 B CN 109657322B CN 201811522183 A CN201811522183 A CN 201811522183A CN 109657322 B CN109657322 B CN 109657322B
Authority
CN
China
Prior art keywords
solid
phase
liquid
debris flow
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201811522183.8A
Other languages
English (en)
Other versions
CN109657322A (zh
Inventor
谢云旭
周公旦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Mountain Hazards and Environment IMHE of CAS
Original Assignee
Institute of Mountain Hazards and Environment IMHE of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Mountain Hazards and Environment IMHE of CAS filed Critical Institute of Mountain Hazards and Environment IMHE of CAS
Priority to CN201811522183.8A priority Critical patent/CN109657322B/zh
Publication of CN109657322A publication Critical patent/CN109657322A/zh
Application granted granted Critical
Publication of CN109657322B publication Critical patent/CN109657322B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种固液多相适用于泥石流的动力学数值模拟方法，其采用固液多相泥石流动力过程的数值模拟模型，先通过野外科学考察和物理及力学实验获取泥石流沟道启动、运动、堆积动力学信息；再通过地理信息系统或者高精度3D地形扫描仪确定计算区域内高精度的地形数据及物源数据，通过网格坐标转化表达为(x,y,h)网格数据,x,y,h分别为地形点的x,y坐标和高程值h；接着将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体、流体参数，以获得本数值方法预处理阶段所有需要准备的参数。本发明提高了泥石流动力过程数值模拟研究的科学性和可靠性，提高了灾害防治的针对性、增强预防效果为泥石流减灾提供技术支撑。

Description

一种固液多相适用于泥石流的动力学数值模拟方法
技术领域
本发明涉及泥石流灾害防治、环境治理技术领域，具体涉及一种固液多相适用于泥石流的动力学数值模拟方法。
背景技术
泥石流是一种固相颗粒级配宽、容重变化范围大的典型固液两相流，是一种常见分布于山区的地质灾害，由其特殊固液两相流复杂流态特征，具有规模大、爆发迅速，来势凶猛，冲击力巨大等特点，不仅可以在沟道中迅速流动，且能冲出沟道进入河流水库，造成水面抬升、河流堵塞形成堰塞坝，严重威胁水利水电工程，且堰塞湖溃坝在下游造成溃决洪水灾害，危害下游居民安全，破坏生态环境。且泥石流堆积区域往往地势平坦、水源充足，山区居民常常在这种稳定的滩地上进行农业生产，建设村落，修筑道路，但是泥石流堆积区域受到潜在泥石流灾害的巨大威胁，故泥石流特别是大规模泥石流已成为当前防灾减灾工作中亟待解决的突出问题。
运用深度平均(Depth-averaged)纳维叶－斯托克斯(Navier-Stokes)方程并伴以适当高精度数学离散格式解决此类问题是一大热门问题，然而目前绝大多数都运用单相流体介质方程，虽然其简单快捷易于计算，但是无法反应出固液体积浓度占比在时间、空间上的变化导致的密度改变，流态改变等特点，其定义的应力条件，也只能是简化后的加和形式，无法反应出多相介质受到的不同应力条件、应力耦合变化。同时在模型源项上，目前绝大多数模型都运用以经典力学为基础的滑块受力库伦模型(Coulomb)或者基于经验统计的流变Voellmy模型，其模型缺点明显无法反应出泥石流此种特殊流体因为流态变化而反应出不同的摩阻力性质进而无法很好的模拟出其特殊流体的动力学特征。
发明内容
针对目前现有单相流体介质模型无法准确反映泥石流复杂多相介质的动力学过程，本发明提供了一种能取代单相流体介质方程中简化流体密度、颗粒分布、颗粒浓度占比、受力条件等缺陷，能在单相模型的基础上将多相介质较为准确的耦合在一起，提高泥石流动力过程数值模拟研究的科学性和可靠性，提高灾害防治的针对性、增强预防效果为泥石流减灾提供技术支撑，更好的服务泥石流减灾的固液多相适用于泥石流的动力学数值模拟方法。
本发明的技术方案如下：
上述的固液多相适用于泥石流的动力学数值模拟方法，是采用固液多相泥石流动力过程的数值模拟模型，先通过野外科学考察和物理及力学实验获取泥石流沟道启动、运动、堆积动力学信息；再通过地理信息系统或者高精度3D地形扫描仪确定计算区域内高精度的地形数据及物源数据，通过网格坐标转化表达为(x,y,h)网格数据,x,y,h分别为地形点的x,y坐标和高程值h；接着将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体、流体参数，以获得本数值方法预处理阶段所有需要准备的参数。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述固液多相泥石流动力过程的数值模拟模型的控制方程基于深度平均积分法的纳维叶－斯托克斯方程，分别伴以适用于泥石流中固相颗粒和液相流体的本构方程能够较好地运用在模拟泥石流动力学过程中；
所述泥石流中固相颗粒及液相浆体的质量方程为：
Figure BDA0001903476250000021
Figure BDA0001903476250000022
上式(1)、(2)中，ρs,ρf分别为固体颗粒和浆体密度，
Figure BDA0001903476250000023
为固体颗粒体积浓度，vs,vf分别为固液两相速度，Δm1,Δm2为外界引起两相质量变化率。
同理，泥石流固液两相动量守恒包括：
Figure BDA0001903476250000024
Figure BDA0001903476250000025
上式(3)、(4)中，
Figure BDA0001903476250000026
分别为两相的体积浓度，Ts,Tf为固液两相的应力张量，fi为两相间的相间作用力。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述固液两相间的相间作用力fi是采用如下模型：
fi＝fb+fl+fd+fv； (18)
其中，上式(18)中，fb代表流体静力学中的浮力，fl,fd,fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力；
由于泥石流速度Vm绝大对数小于30m/s，所对应的雷诺数Rep较小；抬升力，虚拟质量力fl,fv＜＜fd，适合多流态的拖曳力系数及拖曳力模型为：
Figure BDA0001903476250000031
Figure BDA0001903476250000032
Figure BDA0001903476250000033
Figure BDA0001903476250000034
fd＝β(us-uf)；(23)
其中，上式(19)-(23)中，β是固液相间动量交换系数，(us-uf)为固液相间相对速度，Cd拖曳力系数，n为孔隙度，D为颗粒粒径，Rep为雷诺数。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述固液多相泥石流动力过程的数值模拟模型的控制方程为：
Figure BDA0001903476250000035
Figure BDA0001903476250000036
Figure BDA0001903476250000037
Figure BDA0001903476250000041
Figure BDA0001903476250000042
Figure BDA0001903476250000043
Figure BDA0001903476250000044
将上述控制方程组即公式(5)-(10)写为向量格式：
Figure BDA0001903476250000045
其中：
Figure BDA0001903476250000046
Figure BDA0001903476250000051
上式(5)-(12)中，z为处理后的地形高程，g为重力加速度，hs,hf为固液两相的高度；ks,kf为侧向压力；us,vs,uf,vf分别为固液两相速度矢量在x,y方向上分量；Tsx,Tsy,Tfx,Tfy分别为固液两相摩阻力在x,y方向上的分量；fix,fiy分别为固液力在x,y方向上的分量。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述固相摩阻力Tsx,Tsy是采用适用于滑坡碎屑流的模型μ(I)：
Figure BDA0001903476250000052
Figure BDA0001903476250000053
Figure BDA0001903476250000054
Figure BDA0001903476250000055
其中，上式(13)-(16)中，μs,μ2分别表示颗粒体在准静止和高速剪切状态下的内摩擦角，d为颗粒平均粒径，I为起始数，
Figure BDA0001903476250000058
为平均剪切率，Cs为颗粒体积浓度。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述液相摩阻力Tfx,Tyf示采用牛顿体或宾汉体模型：
Figure BDA0001903476250000056
其中，上式(17)中，μ为液体粘度，
Figure BDA0001903476250000057
为深度方向平均剪切率，c为粘聚力。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，所述控制方程即公式(5)-(10)采用高精度非交错格式有限差分格式进行求解，变量上标，整数n表示计算所在时间步，整数下标i,j为所建矩形网格x,y方向上空间网格节点号，如果表示为i,j则为非交错格式表示法，如果为
Figure BDA0001903476250000061
为交错格式表示法，Δx,Δy表示x,y方向上空间步长，Δt代表时间步长；向量形式的控制方程即上述公式(11)-(12)的具体离散过程如下：
1)直接将连续性方程即公式(1)和动量方程即公式(2)通过交错格式中心差分法离散即
Figure BDA0001903476250000062
将离散方程组即公式(24)进行Van-Lee型插值,并将离散方程组在时空重区间
Figure BDA0001903476250000063
进行重积分：
Figure BDA0001903476250000064
上式(25)中的
Figure BDA0001903476250000065
其为tk时刻离散变量交错格式的单元平均，且将单元平均替代后得到：
Figure BDA0001903476250000066
Figure BDA0001903476250000071
对流项即方程(26)的左侧形式采用：
Figure BDA0001903476250000072
2)将交错格式离散后方程组即公式(24)转化为非交错格式离散方程：
Figure BDA0001903476250000081
上述公式(28)迭代特定时间步长后，输出结果包括：
特定精度网格下网格各空间节点处得到多相泥石流深的空间分布即
hm(i,j)＝hs(i,j)+hf(i,j)； (29)
多相泥石流速度的空间分布：
Figure BDA0001903476250000082
多相泥石流颗粒体积浓度的空间分布即
Figure BDA0001903476250000083
多相泥石流密度的空间分布即ρm＝ρsCs+ρf(1-Cs)； (32)
上述公式(24)-(32)中，i,j代表离散后网格节点下标，i为x方向下标，j为y方向下标；hm(i,j)代表多相泥石流体在拉格朗日离散法后均匀矩形网格节点上的泥深，是固相泥深hs和液相泥深hf的加和；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数；Cs(i,j)为固相颗粒浓度，和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数；ρm为多相泥石流密度，为固相颗粒密度ρs和液相悬浊液密度ρf的函数。
所述固液多相适用于泥石流的动力学数值模拟方法，其中，将多相控制方程即公式(5)-(10)，通过离散方程即公式(24)-(28)离散求解，此离散方程即公式(24)-(28)构建了变量在不同时间层的迭代关系，采用可变时间步长进行变量迭代，初始时间步为0.01s,可变时间步满足：
Figure BDA0001903476250000091
其中，上式(32)中，ρ(A)为谱半径，Cr代表克拉数；为满足离散方程的稳定性，Cr满足：
Figure BDA0001903476250000092
所述固液多相适用于泥石流的动力学数值模拟方法，其中，求解上述离散方程即公式(24)-(28)构建的不同时间层变量的迭代关系的具体步骤为：
(1)将需要的参数赋予初始值，包括初始时间步t＝0时间步计算区域内特定精度网格高程图z导入；t＝0时间步初始固相颗粒体积浓度Cs，根据Cs确定固液两相在固定精度网格内的体积占比hs＝hCs；hf＝h(1-Cs)；
(2)在初始时间步t＝0.01s下求解连续性方程，更新计算区域网格各点固液两相体积占比hs,hf,及多相泥石流密度ρm(i,j)＝ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布
Figure BDA0001903476250000093
其中，Vols(i,j)和Volf(i,j)为每个网格中固相和液相的体积，因为网格为规则矩形网格，所以可以简化为hs(i,j),hf(i,j)；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数；
(3)求解计算区域内网格各点x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度
Figure BDA0001903476250000094
其中，us,vs为固相x,y方向上的速度分量；uf,vf为液相x,y方向上的分量；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数；ρs,ρf为固相液相密度且为定值；
(4)更新当前时间t2(t+Δt)＝t1+Δt及更新可变时间步
Figure BDA0001903476250000095
Figure BDA0001903476250000096
其中，t2为下个时间层的时间，t1为上个时间层的时间，Δt是时间步，Cr为克拉数方程，min为一个函数且值是取得两个参数中的较小值；ρ(A)为A的谱半径，为每个时间步的最大速度；按如上顺序继续迭代变量至当前时间大于等于计算总时间，迭代结束；
(5)迭代特定时间步长后，输出结果为：
特定精度网格下网格各离散点固液两相体积占比hs,hf及多相泥石流密度ρm(i,j)＝ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布
Figure BDA0001903476250000101
Figure BDA0001903476250000102
x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度
Figure BDA0001903476250000103
其中，us,vs为固相x,y方向上的速度分量；uf,vf为液相x,y方向上的分量；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数(31)；ρs,ρf为固相液相密度且为定值。
有益效果：
本发明固液多相适用于泥石流的动力学数值模拟方法能在单相模型的基础上将多相介质较为准确的耦合在一起，提高灾害防治的针对性、增强预防效果为泥石流减灾提供技术支撑，更好的服务泥石流减灾；本发明基于深度平均假设的纳维叶－斯托克斯方程能很好的满足泥石流H<<L动力学特征，针对目前现有单相流体介质模型无法准确反映泥石流复杂多相介质的动力学过程，能取代单相流体介质方程中简化流体密度、颗粒分布、颗粒浓度占比、受力条件等缺陷，假设多相介质方程中的不同介质充分掺混且单独定义泥石流固相和液相的受力条件；同时引入多种相互作用力的模型，将两相介质通过流体静力学浮力(Buoyancy force)和多种流体动力学力：托曳力、抬升力、虚拟质量力(Drag force,liftforce,virtual mass force)耦合在一起，较为准确的模拟出固相和液相在运动过程中，因为受力条件不同而出现相对速度效应进而出现因浓度占比时空演变下的不同的流态演变和动力学特征，形成基于固液多相的泥石流动力过程的形成-运动-堆积数值模拟方法。
本发明具体优点体现在以下几方面：
(1)通过数值模拟，研究不同地域,不同浓度,不同地形下泥石流动力学特征，将空间地域分布引入到泥石流研究中，为我国特定地区泥石流灾害防治提供技术支持；
(2)本发明具有简单、高效等特点，适合提供给相关灾害领域工作者使用，提高泥石流防治，环境保护的科学化水平；
(3)本发明具体可变高精度，适用于泥石流动力学特征等优点，针对泥石流动力学启动-运动-堆积等过程，能较好地将整个动力学过程完整的模拟出来，提高泥石流动力学过程研究的科学性和可靠性；
(4)本发明在单相控制方程平台的基础上提高在泥石流防治、环境保护等领域的科学性和可靠性。
附图说明
图1为本发明固液多相适用于泥石流的动力学数值模拟方法中液相采用牛顿流体和宾汉流体的应力随应变率变化图；
图2为本发明固液多相适用于泥石流的动力学数值模拟方法的实施例1中香港东涌大屿山北部Yu Tung路泥石流事件总地形图；
图3为本发明固液多相适用于泥石流的动力学数值模拟方法的实施例1中四种不同初始颗粒浓度此次事件沿程流速分布及高程分布图；
图4为本发明固液多相适用于泥石流的动力学数值模拟方法的实施例1中流动距离100米、400米、420米及最终堆积区域处泥石流速度分布数值实验结果俯视图。
图5为本发明固液多相适用于泥石流的动力学数值模拟方法的实施例1中流动到0秒、20秒、50秒、100秒处泥石流固相颗粒浓度分布数值实验结果俯视图。
其中，(Cs＝0.4)，在颗粒浓度(Cs＝0.4下其最终堆积区域为越过Yu Tung公路几米处，和真实情况基本相符。
具体实施方式
本发明固液多相适用于泥石流的动力学数值模拟方法，是采用固液多相泥石流动力过程的数值模拟模型，其控制方程基于深度平均积分法的纳维叶－斯托克斯方程，分别伴以适用于泥石流中固相颗粒和液相流体的本构方程(本构方程是指控制方程即下述公式(3)(4)右边源项中的Ts和Tf,为适用于泥石流的摩阻力本构关系，是固相应力应变关系或者是液相应力应变率的关系。)能够较好地运用在模拟泥石流动力学过程中；首先，通过野外科学考察和物理及力学实验获取泥石流沟道启动、运动、堆积动力学信息；通过地理信息系统或者高精度3D地形扫描仪确定计算区域内高精度的地形数据及物源数据(即区域内每个离散点都有对应的高程值)，通过网格坐标转化表达为(x,y,h)网格数据,x,y,h分为为地形点的x,y坐标和高程值h；将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体及流体参数：颗粒级配(GSD)、固体颗粒浓度
Figure BDA0001903476250000121
中值粒径d5O、孔隙率nO、内摩擦角
Figure BDA0001903476250000122
(准静止状态
Figure BDA0001903476250000123
高速剪切状态
Figure BDA0001903476250000124
)、初始含水量θO、流体粘度μ，以获得本数值方法预处理阶段所有需要准备的参数。
其中，上述的泥石流中固相颗粒及液相浆体的质量方程为：
Figure BDA0001903476250000125
Figure BDA0001903476250000126
上式(1)、(2)中，ρs,ρf分别为固体颗粒和浆体密度，
Figure BDA0001903476250000127
为固体颗粒体积浓度，vs,vf分别为固液两相速度，Δm1,Δm2为外界引起两相质量变化率。
同理，泥石流固液两相动量守恒包括：
Figure BDA0001903476250000128
Figure BDA0001903476250000129
上式(3)、(4)中，
Figure BDA00019034762500001210
分别为两相的体积浓度，Ts,Tf为固液两相的应力张量，fi为两相间的相间作用力。
其中，本发明是将泥石流混合物质量守恒动量守恒联立并通过深度平均假设简化方程，以推导出所述固液多相泥石流动力过程的数值模拟模型的控制方程为：
Figure BDA00019034762500001211
Figure BDA00019034762500001212
Figure BDA0001903476250000131
Figure BDA0001903476250000132
Figure BDA0001903476250000133
Figure BDA0001903476250000134
将上述控制方程组(5)-(10)写为向量格式：
Figure BDA0001903476250000135
其中：
Figure BDA0001903476250000136
Figure BDA0001903476250000141
上式(5)-(12)中，z为处理后的地形高程；g为重力加速度；hs,hf为固液两相的高度；ks,kf为侧向压力；us,vs,uf,vf分别为固液两相速度矢量在x,y方向上分量；Tsx,Tsy,Tfx,Tfy分别为固液两相摩阻力在x,y方向上的分量；fix,fiy分别为固液力在x,y方向上的分量。
对于上述的固相摩阻力Tsx,Tsy,本发明将采用适用于滑坡碎屑流的模型(μ(I)rheology)：
Figure BDA0001903476250000142
Figure BDA0001903476250000143
Figure BDA0001903476250000144
Figure BDA0001903476250000145
其中，上式(13)-(16)中，μs,μ2分别表示颗粒体在准静止和高速剪切状态下的内摩擦角，d为颗粒平均粒径，I为起始数，
Figure BDA0001903476250000148
为平均剪切率，Cs为颗粒体积浓度。
如图1所示，对于上述的液相摩阻力Tfx,Tyf,本发明将采用牛顿体或宾汉体模型：
Figure BDA0001903476250000146
其中，上式(17)中，μ为液体粘度，
Figure BDA0001903476250000147
为深度方向平均剪切率，c为粘聚力；
对于上述的固液两相间的相间作用力fi，本发明将采用如下模型：
fi＝fb+fl+fd+fv； (18)
其中，上式(18)中，fb代表流体静力学中的浮力，fl,fd,fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力；由于泥石流速度相对较小(Vm绝大对数小于30m/s)，所对应的雷诺数Rep较小；抬升力，虚拟质量力fl,fv＜＜fd；为避免以往模型考了间断的拖曳力系数Cd导致拖曳力不连续，拖曳力只能考虑部分流态等问题，本发明考虑适合多流态的拖曳力系数及拖曳力模型：
Figure BDA0001903476250000151
Figure BDA0001903476250000152
Figure BDA0001903476250000153
Figure BDA0001903476250000154
fd＝β(us-uf)； (23)
其中，上式(19)-(23)中，β是固液相间动量交换系数，(us-uf)为固液相间相对速度，Cd拖曳力系数，n为孔隙度，D为颗粒粒径，Rep为雷诺数。
本发明控制方程的连续方程组即上述公式(5)-(10)采用高精度非交错格式有限差分格式进行求解，变量上标，整数n表示计算所在时间步，整数下标i,j为所建矩形网格x,y方向上空间网格节点号，如果表示为i,j则为非交错格式表示法，如果为
Figure BDA0001903476250000155
为交错格式表示法，Δx,Δy表示x,y方向上空间步长，Δt代表时间步长；向量形式的控制方程即上述公式(11)-(12)的具体离散过程如下：
1)直接将连续性方程即上述公式(1)和动量方程即上述公式(2)通过交错格式中心差分法离散:
Figure BDA0001903476250000161
将离散后方程组即公式(24)进行Van-Lee型插值,并将离散方程在时空重区间
Figure BDA0001903476250000162
进行重积分：
Figure BDA0001903476250000163
其中，上式(25)中的
Figure BDA0001903476250000164
为tk时刻离散变量交错格式的单元平均，且将单元平均替代后得到：
Figure BDA0001903476250000165
Figure BDA0001903476250000171
对流项(即方程(26)的左侧形式)采用：
Figure BDA0001903476250000172
2)将交错格式离散后方程组即上述公式(24)转化为非交错格式离散方程：
Figure BDA0001903476250000173
Figure BDA0001903476250000181
上述公式(28)迭代特定时间步长后，输出结果包括：
特定精度网格下网格各空间节点处得到多相泥石流深的空间分布即
hm(i,j)＝hs(i,j)+hf(i,j)； (29)
多相泥石流速度的空间分布即
Figure BDA0001903476250000182
多相泥石流颗粒体积浓度的空间分布即
Figure BDA0001903476250000183
多相泥石流密度的空间分布即ρm＝ρsCs+ρf(1-Cs)； (32)
上述公式(24)-(32)中，i,j代表离散后网格节点下标，i为x方向下标，j为y方向下标；hm(i,j)代表多相泥石流体在拉格朗日离散法后均匀矩形网格节点上的泥深，是固相泥深hs和液相泥深hf的加和；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数即公式(30)；Cs(i,j)为固相颗粒浓度，和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31)；ρm为多相泥石流密度，为固相颗粒密度ρs和液相悬浊液密度ρf的函数即公式(32)。上述公式(11)-(12)中的控制方程向量形式W,F(W),G(W)都参考方程离散方程即上述公式(24)-(28)进行离散化。
其中，本发明将多相控制方程即上述公式(5)-(10)，通过离散方程即上述公式(24)-(28)离散求解，此离散方程即公式(24)-(28)构建了变量在不同时间层的迭代关系；本发明采用可变时间步长进行变量迭代，初始时间步t＝0.01s,可变时间步满足：
Figure BDA0001903476250000184
其中，上式(32)中，ρ(A)为谱半径，Cr代表克拉数；为满足离散方程的稳定性，Cr满足：
Figure BDA0001903476250000185
求解上述离散方程即公式(24)-(28)构建的不同时间层变量的迭代关系的具体步骤为：
(1)将需要的参数赋予初始值，包括在初始时间步t＝0时间步计算区域内特定精度网格高程图z导入；t＝0时间步初始固相颗粒体积浓度Cs，根据Cs确定固液两相在固定精度网格内的体积占比hs＝hCs；hf＝h(1-Cs)；若考虑基底侵蚀沿程规模放大效应，考虑可侵蚀层高程图Ze导入；及通过物理及力学确定固体颗粒参数μs,μ2,IO,d5O导入μ(I)本构关系；液相粘度参数μ导入牛顿流体本构方程，及固液相间作用力参数确定；计算参数导入包括总时间、克拉数、可变精度、结果输出间隔时间。
(2)在初始时间步t＝0.01s下求解连续性方程，更新计算区域网格各点固液两相体积占比hs,hf,及多相泥石流密度ρm(i,j)＝ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布
Figure BDA0001903476250000191
其中，Vols(i,j)和Volf(i,j)为每个网格中固相和液相的体积，因为网格为规则矩形网格，所以可以简化为hs(i,j),hf(i,j)；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31)；
(3)求解计算区域内网格各点x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度
Figure BDA0001903476250000192
其中，us,vs为固相x,y方向上的速度分量；uf,vf为液相x,y方向上的分量；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31)；ρs,ρf为固相液相密度且为定值。
(4)更新当前时间t2(t+Δt)＝t1+Δt及更新可变时间步
Figure BDA0001903476250000193
Figure BDA0001903476250000194
其中，t2为下个时间层的时间，t1为上个时间层的时间，Δt是时间步，Cr为克拉数方程即公式(34)，min为一个函数且值是取得两个参数中的较小值；ρ(A)为A的谱半径，为每个时间步的最大速度；按如上顺序继续迭代变量至当前时间大于等于计算总时间，迭代结束。
(5)迭代特定时间步长后，输出结果为：
特定精度网格下网格各离散点固液两相体积占比hs,hf及多相泥石流密度ρm(i,j)＝ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布
Figure BDA0001903476250000201
Figure BDA0001903476250000202
x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度
Figure BDA0001903476250000203
其中us,vs为固相x,y方向上的速度分量；uf,vf为液相x,y方向上的分量；vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数；Cs(i,j)为固相颗粒浓度，且和每个矩形网格中的固相液相占比有关系，为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31)；ρs,ρf为固相液相密度且为定值。
泥石流动力学方程是一个二维带有适用于泥石流动力学特征本构关系的浅水波方程组，浅水波方程组由连续性方程即上述公式(5)-(6)(又叫质量守恒方程)和空间中x,y方向上动量守恒方程即上述公式(7)-(10)组成，泥石流动力学方程通过变量泥深h的变化来代表空间中z方向z坐标的变化，故可以认为是三维流体运动方程。由浅水波方程本质是非线性复杂双曲型偏微分方程组，本发明采用一种高精度非交错格式中心有限差分法离散求解，此方法和原有方法相比不需要使用繁琐的尼曼求解器求解特征信息，而使用较为简单的单元平均法求得最终的迭代关系方程，在时间和空间都保证二阶精度的同时简便、快捷地计算泥石流动力学过程，可以较为准确的反演泥石流动力学演化过程，通过与真实事件结果对比，在一定程度上进行模型和参数修正，可以很好的揭示泥石流动力学整个过程，为泥石流动力学研究和泥石流放在减灾提供数据和技术支持。
下面结合具体实施例对本发明固液多相适用于泥石流的动力学数值模拟方法作进一步描述：
实施例1
如图2-4，利用本发明的对该事件泥石流动力学模拟，具体步骤如下：Yu Tung滑坡泥石流是香港2008年7月因强降雨条件导致的大屿山北部Yu Tung路沿途山路滑坡所引起的泥石流。这次泥石流是典型的强降雨型沟道泥石流，具有规模大、爆发迅速，来势凶猛，颗粒含量低等特点。此次泥石流方量估算为2,350m2，运动距离为600米，此事件导致Yu Tung公路严重受损，过往车辆受到不同程度损伤，社会影响严重。沟道100及400米处都设有超仰角数据采集及录像系统，将采集到的速度分布信息以及不同颗粒浓度下的数值实验结果做了对比。
表格1:固液多相数值计算参数表
固相颗粒参数 值
颗粒体积浓度C<sub>s</sub>(m) 0.2-0.5
μ<sub>s</sub> 0.22
μ<sub>2</sub> 0.3
I<sub>O</sub> 0.28
中值粒径d<sub>5O</sub>(m) 0.001
液相流体参数 值
粘度(Pa·s) 0.005
将本发明预处理需要的参数导入后，启动泥石流运动计算，可以根据程序附带运动变化监测器监控泥石流体相关动力学信息，并按照预先设定好的时间间隔输出一定格式的结果，通过后处理软件还原整个事件泥石流运动过程，动画中可以继续，暂停，后退，加速播放，输出为特定格式动画等，数据加载过程中自动识别计算区域内网格文件空间坐标系，并且通过一定色泽调整实现整个泥石流运动模拟情景。
本发明在单相的基础上提出的多相适用于泥石流的控制方程即固相考虑适合颗粒流系统的μ(I)流变本构方程，液相通过牛顿或者宾汉体表示；对于相间作用力，考虑适用多种流态下的泥石流体；可变高精度数值算法，在保证精度的同时，避免求解尼曼问题，较为简便地离散方程并求解，精确快速地模拟泥石流动力学过程；在单相方程的基础上更好的服务于泥石流防治，环境保护等领域，具有更高的科学性，时效性，可靠性。

Claims (1)

1.一种固液多相适用于泥石流的动力学数值模拟方法，其特征在于，采用固液多相泥石流动力过程的数值模拟模型，先通过野外科学考察和物理及力学实验获取泥石流沟道启动、运动、堆积动力学信息；
再通过地理信息系统或者高精度3D地形扫描仪确定计算区域内高精度的地形数据及物源数据，通过网格坐标转化表达为(x，y，h)网格数据，x，y，h分别为地形点的x，y坐标和高程值h；
接着将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体、流体参数，以获得本数值方法预处理阶段所有需要准备的参数；
所述固液多相泥石流动力过程的数值模拟模型的控制方程基于深度平均积分法的纳维叶-斯托克斯方程，分别伴以适用于泥石流中固相颗粒和液相流体的本构方程能够较好地运用在模拟泥石流动力学过程中；
所述泥石流中固相颗粒及液相浆体的质量方程为：
Figure FDA0002541837600000011
Figure FDA0002541837600000012
上式(1)、(2)中，ρs，ρf分别为固体颗粒和浆体密度，
Figure FDA0002541837600000013
为固体颗粒体积浓度，vs，vf分别为固液两相速度，Δm1，Δm2为外界引起两相质量变化率；
同理，泥石流固液两相动量守恒包括：
Figure FDA0002541837600000014
Figure FDA0002541837600000015
上式(3)、(4)中，
Figure FDA0002541837600000016
分别为两相的体积浓度，Ts，Tf为固液两相的应力张量，fi为两相间的相间作用力；
所述固液两相间的相间作用力fi是采用如下模型：
fi＝fb+fl+fd+fv； (18)
其中，上式(18)中，fb代表流体静力学中的浮力，fl，fd，fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力；
由于泥石流速度Vm绝大对数小于30m/s，所对应的雷诺数Rep较小；抬升力，虚拟质量力fl，fv＜＜fd，适合多流态的拖曳力系数及拖曳力模型为：
Figure FDA0002541837600000021
Figure FDA0002541837600000022
Figure FDA0002541837600000023
Figure FDA0002541837600000024
fd＝β(us-uf)； (23)
上式(19)-(23)中，vols为固相颗粒的体积，β是固液相间动量交换系数，(us-uf)为固液相间相对速度，Cd为拖曳力系数，n为孔隙度，D为颗粒粒径，Rep为雷诺数；
所述固液多相泥石流动力过程的数值模拟模型的控制方程为：
Figure FDA0002541837600000025
Figure FDA0002541837600000026
Figure FDA0002541837600000027
Figure FDA0002541837600000028
Figure FDA0002541837600000029
Figure FDA00025418376000000210
Figure FDA0002541837600000031
将上述控制方程组即公式(5)-(10)写为向量格式：
Figure FDA0002541837600000032
其中：
Figure FDA0002541837600000033
上式(5)-(12)中，z为处理后的地形高程，g为重力加速度，hs，hf为固液两相的高度；ks，kf为侧向压力；us，vs，uf，vf分别为固液两相速度矢量在x，y方向上分量；Tsx，Tsy，Tfx，Tfy分别为固液两相摩阻力在x，y方向上的分量；fix，fiy分别为固液力在x，y方向上的分量。
CN201811522183.8A 2018-12-13 2018-12-13 一种固液多相适用于泥石流的动力学数值模拟方法 Expired - Fee Related CN109657322B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811522183.8A CN109657322B (zh) 2018-12-13 2018-12-13 一种固液多相适用于泥石流的动力学数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811522183.8A CN109657322B (zh) 2018-12-13 2018-12-13 一种固液多相适用于泥石流的动力学数值模拟方法

Publications (2)

Publication Number Publication Date
CN109657322A CN109657322A (zh) 2019-04-19
CN109657322B true CN109657322B (zh) 2020-08-21

Family

ID=66114316

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811522183.8A Expired - Fee Related CN109657322B (zh) 2018-12-13 2018-12-13 一种固液多相适用于泥石流的动力学数值模拟方法

Country Status (1)

Country Link
CN (1) CN109657322B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112749468B (zh) * 2019-10-29 2022-09-02 中国石油化工股份有限公司 固相悬浮物在油气储集层孔隙中通过能力的数值模拟方法
CN112699618B (zh) * 2020-12-18 2023-01-17 赣江新区澳博颗粒科技研究院有限公司 一种离子型稀土矿原地浸矿过程数值模拟方法
CN112733472B (zh) * 2021-01-11 2021-10-26 中国科学院、水利部成都山地灾害与环境研究所 一种泥石流坡面物源启动量动态计算方法及系统
CN112818574B (zh) * 2021-01-27 2022-10-14 江西理工大学 模拟泥石流起动形成、流动发展和再次淤积的数值方法
CN112784504B (zh) * 2021-01-28 2022-08-30 中国科学院、水利部成都山地灾害与环境研究所 一种强耦合固液多相流数值模拟方法
CN114239352B (zh) * 2021-12-14 2024-04-30 西南交通大学 一种深度积分流体模型和块体系统的流固耦合方法
CN117829031B (zh) * 2024-03-01 2024-06-18 中国科学院、水利部成都山地灾害与环境研究所 考虑径流和泥石流相互作用的动力学模拟方法
CN117973271B (zh) * 2024-04-01 2024-06-14 中国科学院、水利部成都山地灾害与环境研究所 一种考虑深度信息的固-液两相泥石流数值模拟方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107506566A (zh) * 2017-10-16 2017-12-22 中国科学院、水利部成都山地灾害与环境研究所 一种新型泥石流动力学数值模拟分析方法及系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EA010085B1 (ru) * 2004-09-10 2008-06-30 Эксонмобил Апстрим Рисерч Компани Способ оценивания свойств осадочного бассейна путем численного моделирования процессов осадконакопления

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107506566A (zh) * 2017-10-16 2017-12-22 中国科学院、水利部成都山地灾害与环境研究所 一种新型泥石流动力学数值模拟分析方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于固液两相流模型的泥石流流速垂向分布研究;乐茂华 等;《水利学报》;20170228;第48卷(第2期);全文 *

Also Published As

Publication number Publication date
CN109657322A (zh) 2019-04-19

Similar Documents

Publication Publication Date Title
CN109657322B (zh) 一种固液多相适用于泥石流的动力学数值模拟方法
CN107506566B (zh) 一种泥石流动力学数值模拟分析方法及系统
Pastor et al. Application of a SPH depth-integrated model to landslide run-out analysis
CN106529198B (zh) 一种泥石流全过程数值模拟及数值计算方法
Zakeri Submarine debris flow impact on suspended (free-span) pipelines: Normal and longitudinal drag forces
Haddad et al. A SPH depth integrated model for Popocatépetl 2001 lahar (Mexico): Sensitivity analysis and runout simulation
Dulal et al. Numerical computation of free meandering channels with the application of slump blocks on the outer bends
CN109992830B (zh) 基于物质点方法的山体滑坡灾害场景模拟方法
Mao et al. A three-phases model for the simulation of landslide-generated waves using the improved conservative level set method
Wang et al. Tsunami Squares modeling of the 2007 Dayantang landslide generated waves considering the effects in slide/water interactions
Marangoz et al. Two-dimensional modeling of flood wave propagation in residential areas after a dam break with application of diffusive and dynamic wave approaches
Issakhov et al. Numerical simulation of dam break waves on movable beds for various forms of the obstacle by VOF method
Froiio et al. A numerical experiment of backward erosion piping: kinematics and micromechanics
Procter et al. Evaluation of Titan2D modelling forecasts for the 2007 Crater Lake break-out lahar, Mt. Ruapehu, New Zealand
Ceccato et al. The effect of the front inclination on the impact forces transmitted by granular flows to rigid structures
Cleary et al. Combining digital terrain and surface textures with large-scale particle-based computational models to predict dam collapse and landslide events
Turmel et al. Parametric analysis of the mobility of debris from flow slides in sensitive clays
Sanchez et al. Modelling of short runout propagation landslides and debris flows
Abdelrazek et al. Numerical simulation of a small-scale snow avalanche tests using non-Newtonian SPH model
Roy et al. A study on hydrodynamic and morphological behavior of Padma river using Delft3d model
Kafle Advanced dynamic simulations of landslide generated tsunami, submarine mass movement and obstacle interaction
Basirat et al. Eulerian–Eulerian model application to simulate scouring downstream of sluice gate
Tachibana et al. Characterization of transition from Darcy to non-Darcy flow with 3D pore-level simulations
Ma et al. Dynamic Simulation and Analysis of Large‐Scale Debris Flow Field
Rattia et al. Numerical simulation of scour below pipelines using flexible mesh methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200821