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

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

Info

Publication number
CN109657322A
CN109657322A CN201811522183.8A CN201811522183A CN109657322A CN 109657322 A CN109657322 A CN 109657322A CN 201811522183 A CN201811522183 A CN 201811522183A CN 109657322 A CN109657322 A CN 109657322A
Authority
CN
China
Prior art keywords
solid
liquid
mud
phase
rock flow
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.)
Granted
Application number
CN201811522183.8A
Other languages
English (en)
Other versions
CN109657322B (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

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06Power 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;接着将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体、流体参数,以获得本数值方法预处理阶段所有需要准备的参数。
所述固液多相适用于泥石流的动力学数值模拟方法,其中,所述固液多相泥石流动力过程的数值模拟模型的控制方程基于深度平均积分法的纳维叶-斯托克斯方程,分别伴以适用于泥石流中固相颗粒和液相流体的本构方程能够较好地运用在模拟泥石流动力学过程中;
所述泥石流中固相颗粒及液相浆体的质量方程为:
上式(1)、(2)中,ρsf分别为固体颗粒和浆体密度,为固体颗粒体积浓度,vs,vf分别为固液两相速度,Δm1,Δm2为外界引起两相质量变化率。
同理,泥石流固液两相动量守恒包括:
上式(3)、(4)中,分别为两相的体积浓度,Ts,Tf为固液两相的应力张量,fi为两相间的相间作用力。
所述固液多相适用于泥石流的动力学数值模拟方法,其中,所述固液两相间的相间作用力fi是采用如下模型:
fi=fb+fl+fd+fv; (18)
其中,上式(18)中,fb代表流体静力学中的浮力,fl,fd,fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力;
由于泥石流速度Vm绝大对数小于30m/s,所对应的雷诺数Rep较小;抬升力,虚拟质量力fl,fv<<fd,适合多流态的拖曳力系数及拖曳力模型为:
fd=β(us-uf);(23)
其中,上式(19)-(23)中,β是固液相间动量交换系数,(us-uf)为固液相间相对速度,Cd拖曳力系数,n为孔隙度,D为颗粒粒径,Rep为雷诺数。
所述固液多相适用于泥石流的动力学数值模拟方法,其中,所述固液多相泥石流动力过程的数值模拟模型的控制方程为:
将上述控制方程组即公式(5)-(10)写为向量格式:
其中:
上式(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):
其中,上式(13)-(16)中,μs2分别表示颗粒体在准静止和高速剪切状态下的内摩擦角,d为颗粒平均粒径,I为起始数,为平均剪切率,Cs为颗粒体积浓度。
所述固液多相适用于泥石流的动力学数值模拟方法,其中,所述液相摩阻力Tfx,Tyf示采用牛顿体或宾汉体模型:
其中,上式(17)中,μ为液体粘度,为深度方向平均剪切率,c为粘聚力。
所述固液多相适用于泥石流的动力学数值模拟方法,其中,所述控制方程即公式(5)-(10)采用高精度非交错格式有限差分格式进行求解,变量上标,整数n表示计算所在时间步,整数下标i,j为所建矩形网格x,y方向上空间网格节点号,如果表示为i,j则为非交错格式表示法,如果为为交错格式表示法,Δx,Δy表示x,y方向上空间步长,Δt代表时间步长;向量形式的控制方程即上述公式(11)-(12)的具体离散过程如下:
1)直接将连续性方程即公式(1)和动量方程即公式(2)通过交错格式中心差分法离散即
将离散方程组即公式(24)进行Van-Lee型插值,并将离散方程组在时空重区间进行重积分:
上式(25)中的其为tk时刻离散变量交错格式的单元平均,且将单元平均替代后得到:
对流项即方程(26)的左侧形式采用:
2)将交错格式离散后方程组即公式(24)转化为非交错格式离散方程:
上述公式(28)迭代特定时间步长后,输出结果包括:
特定精度网格下网格各空间节点处得到多相泥石流深的空间分布即
hm(i,j)=hs(i,j)+hf(i,j); (29)
多相泥石流速度的空间分布:
多相泥石流颗粒体积浓度的空间分布即
多相泥石流密度的空间分布即ρm=ρsCsf(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,可变时间步满足:
其中,上式(32)中,ρ(A)为谱半径,Cr代表克拉数;为满足离散方程的稳定性,Cr满足:
所述固液多相适用于泥石流的动力学数值模拟方法,其中,求解上述离散方程即公式(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)),体积浓度空间分布其中,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)及泥石流多相平均速度其中,us,vs为固相x,y方向上的速度分量;uf,vf为液相x,y方向上的分量;vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数;Cs(i,j)为固相颗粒浓度,且和每个矩形网格中的固相液相占比有关系,为固相深度hs(i,j)和液相深度hf(i,j)的一个函数;ρsf为固相液相密度且为定值;
(4)更新当前时间t2(t+Δt)=t1+Δt及更新可变时间步 其中,t2为下个时间层的时间,t1为上个时间层的时间,Δt是时间步,Cr为克拉数方程,min为一个函数且值是取得两个参数中的较小值;ρ(A)为A的谱半径,为每个时间步的最大速度;按如上顺序继续迭代变量至当前时间大于等于计算总时间,迭代结束;
(5)迭代特定时间步长后,输出结果为:
特定精度网格下网格各离散点固液两相体积占比hs,hf及多相泥石流密度ρm(i,j)=ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布 x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度其中,us,vs为固相x,y方向上的速度分量;uf,vf为液相x,y方向上的分量;vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数;Cs(i,j)为固相颗粒浓度,且和每个矩形网格中的固相液相占比有关系,为固相深度hs(i,j)和液相深度hf(i,j)的一个函数(31);ρsf为固相液相密度且为定值。
有益效果:
本发明固液多相适用于泥石流的动力学数值模拟方法能在单相模型的基础上将多相介质较为准确的耦合在一起,提高灾害防治的针对性、增强预防效果为泥石流减灾提供技术支撑,更好的服务泥石流减灾;本发明基于深度平均假设的纳维叶-斯托克斯方程能很好的满足泥石流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)、固体颗粒浓度中值粒径d5O、孔隙率nO、内摩擦角(准静止状态高速剪切状态)、初始含水量θO、流体粘度μ,以获得本数值方法预处理阶段所有需要准备的参数。
其中,上述的泥石流中固相颗粒及液相浆体的质量方程为:
上式(1)、(2)中,ρsf分别为固体颗粒和浆体密度,为固体颗粒体积浓度,vs,vf分别为固液两相速度,Δm1,Δm2为外界引起两相质量变化率。
同理,泥石流固液两相动量守恒包括:
上式(3)、(4)中,分别为两相的体积浓度,Ts,Tf为固液两相的应力张量,fi为两相间的相间作用力。
其中,本发明是将泥石流混合物质量守恒动量守恒联立并通过深度平均假设简化方程,以推导出所述固液多相泥石流动力过程的数值模拟模型的控制方程为:
将上述控制方程组(5)-(10)写为向量格式:
其中:
上式(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):
其中,上式(13)-(16)中,μs2分别表示颗粒体在准静止和高速剪切状态下的内摩擦角,d为颗粒平均粒径,I为起始数,为平均剪切率,Cs为颗粒体积浓度。
如图1所示,对于上述的液相摩阻力Tfx,Tyf,本发明将采用牛顿体或宾汉体模型:
其中,上式(17)中,μ为液体粘度,为深度方向平均剪切率,c为粘聚力;
对于上述的固液两相间的相间作用力fi,本发明将采用如下模型:
fi=fb+fl+fd+fv; (18)
其中,上式(18)中,fb代表流体静力学中的浮力,fl,fd,fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力;由于泥石流速度相对较小(Vm绝大对数小于30m/s),所对应的雷诺数Rep较小;抬升力,虚拟质量力fl,fv<<fd;为避免以往模型考了间断的拖曳力系数Cd导致拖曳力不连续,拖曳力只能考虑部分流态等问题,本发明考虑适合多流态的拖曳力系数及拖曳力模型:
fd=β(us-uf); (23)
其中,上式(19)-(23)中,β是固液相间动量交换系数,(us-uf)为固液相间相对速度,Cd拖曳力系数,n为孔隙度,D为颗粒粒径,Rep为雷诺数。
本发明控制方程的连续方程组即上述公式(5)-(10)采用高精度非交错格式有限差分格式进行求解,变量上标,整数n表示计算所在时间步,整数下标i,j为所建矩形网格x,y方向上空间网格节点号,如果表示为i,j则为非交错格式表示法,如果为为交错格式表示法,Δx,Δy表示x,y方向上空间步长,Δt代表时间步长;向量形式的控制方程即上述公式(11)-(12)的具体离散过程如下:
1)直接将连续性方程即上述公式(1)和动量方程即上述公式(2)通过交错格式中心差分法离散:
将离散后方程组即公式(24)进行Van-Lee型插值,并将离散方程在时空重区间进行重积分:
其中,上式(25)中的为tk时刻离散变量交错格式的单元平均,且将单元平均替代后得到:
对流项(即方程(26)的左侧形式)采用:
2)将交错格式离散后方程组即上述公式(24)转化为非交错格式离散方程:
上述公式(28)迭代特定时间步长后,输出结果包括:
特定精度网格下网格各空间节点处得到多相泥石流深的空间分布即
hm(i,j)=hs(i,j)+hf(i,j); (29)
多相泥石流速度的空间分布即
多相泥石流颗粒体积浓度的空间分布即
多相泥石流密度的空间分布即ρm=ρsCsf(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,可变时间步满足:
其中,上式(32)中,ρ(A)为谱半径,Cr代表克拉数;为满足离散方程的稳定性,Cr满足:
求解上述离散方程即公式(24)-(28)构建的不同时间层变量的迭代关系的具体步骤为:
(1)将需要的参数赋予初始值,包括在初始时间步t=0时间步计算区域内特定精度网格高程图z导入;t=0时间步初始固相颗粒体积浓度Cs,根据Cs确定固液两相在固定精度网格内的体积占比hs=hCs;hf=h(1-Cs);若考虑基底侵蚀沿程规模放大效应,考虑可侵蚀层高程图Ze导入;及通过物理及力学确定固体颗粒参数μs2,IO,d5O导入μ(I)本构关系;液相粘度参数μ导入牛顿流体本构方程,及固液相间作用力参数确定;计算参数导入包括总时间、克拉数、可变精度、结果输出间隔时间。
(2)在初始时间步t=0.01s下求解连续性方程,更新计算区域网格各点固液两相体积占比hs,hf,及多相泥石流密度ρm(i,j)=ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布其中,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)及泥石流多相平均速度
其中,us,vs为固相x,y方向上的速度分量;uf,vf为液相x,y方向上的分量;vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数;Cs(i,j)为固相颗粒浓度,且和每个矩形网格中的固相液相占比有关系,为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31);ρsf为固相液相密度且为定值。
(4)更新当前时间t2(t+Δt)=t1+Δt及更新可变时间步 其中,t2为下个时间层的时间,t1为上个时间层的时间,Δt是时间步,Cr为克拉数方程即公式(34),min为一个函数且值是取得两个参数中的较小值;ρ(A)为A的谱半径,为每个时间步的最大速度;按如上顺序继续迭代变量至当前时间大于等于计算总时间,迭代结束。
(5)迭代特定时间步长后,输出结果为:
特定精度网格下网格各离散点固液两相体积占比hs,hf及多相泥石流密度ρm(i,j)=ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布 x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度其中us,vs为固相x,y方向上的速度分量;uf,vf为液相x,y方向上的分量;vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数;Cs(i,j)为固相颗粒浓度,且和每个矩形网格中的固相液相占比有关系,为固相深度hs(i,j)和液相深度hf(i,j)的一个函数即公式(31);ρsf为固相液相密度且为定值。
泥石流动力学方程是一个二维带有适用于泥石流动力学特征本构关系的浅水波方程组,浅水波方程组由连续性方程即上述公式(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 (9)

1.一种固液多相适用于泥石流的动力学数值模拟方法,其特征在于,采用固液多相泥石流动力过程的数值模拟模型,先通过野外科学考察和物理及力学实验获取泥石流沟道启动、运动、堆积动力学信息;再通过地理信息系统或者高精度3D地形扫描仪确定计算区域内高精度的地形数据及物源数据,通过网格坐标转化表达为(x,y,h)网格数据,x,y,h分别为地形点的x,y坐标和高程值h;接着将计算区域内适量沟道内土体样本通过室内物理及力学实验确定和估算土体、流体参数,以获得本数值方法预处理阶段所有需要准备的参数。
2.如权利要求1所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述固液多相泥石流动力过程的数值模拟模型的控制方程基于深度平均积分法的纳维叶-斯托克斯方程,分别伴以适用于泥石流中固相颗粒和液相流体的本构方程能够较好地运用在模拟泥石流动力学过程中;
所述泥石流中固相颗粒及液相浆体的质量方程为:
上式(1)、(2)中,ρs,ρf分别为固体颗粒和浆体密度,为固体颗粒体积浓度,vs,vf分别为固液两相速度,Δm1,Δm2为外界引起两相质量变化率。
同理,泥石流固液两相动量守恒包括:
上式(3)、(4)中,分别为两相的体积浓度,Ts,Tf为固液两相的应力张量,fi为两相间的相间作用力。
3.如权利要求2所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述固液两相间的相间作用力fi是采用如下模型:
fi=fb+fl+fd+fv; (18)
其中,上式(18)中,fb代表流体静力学中的浮力,fl,fd,fv分别代表流体动力学中的抬升力、托曳力及虚拟质量力;
由于泥石流速度Vm绝大对数小于30m/s,所对应的雷诺数Rep较小;抬升力,虚拟质量力fl,fv<<fd,适合多流态的拖曳力系数及拖曳力模型为:
fd=β(us-uf); (23)
其中,上式(19)-(23)中,β是固液相间动量交换系数,(us-uf)为固液相间相对速度,Cd拖曳力系数,n为孔隙度,D为颗粒粒径,Rep为雷诺数。
4.如权利要求2所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述固液多相泥石流动力过程的数值模拟模型的控制方程为:
将上述控制方程组即公式(5)-(10)写为向量格式:
其中:
上式(5)-(12)中,z为处理后的地形高程,g为重力加速度,hs,hf为固液两相的高度;ks,kf为侧向压力;us,vs,uf,vf分别为固液两相速度矢量在x,y方向上分量;Tsx,Tsy,Tfx,Tfy分别为固液两相摩阻力在x,y方向上的分量;fix,fiy分别为固液力在x,y方向上的分量。
5.如权利要求4所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述固相摩阻力Tsx,Tsy是采用适用于滑坡碎屑流的模型μ(I):
其中,上式(13)-(16)中,μs,μ2分别表示颗粒体在准静止和高速剪切状态下的内摩擦角,d为颗粒平均粒径,I为起始数,为平均剪切率,Cs为颗粒体积浓度。
6.如权利要求4所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述液相摩阻力Tfx,Tyf示采用牛顿体或宾汉体模型:
其中,上式(17)中,μ为液体粘度,为深度方向平均剪切率,c为粘聚力。
7.如权利要求4所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,所述控制方程即公式(5)-(10)采用高精度非交错格式有限差分格式进行求解,变量上标,整数n表示计算所在时间步,整数下标i,j为所建矩形网格x,y方向上空间网格节点号,如果表示为i,j则为非交错格式表示法,如果为 为交错格式表示法,Δx,Δy表示x,y方向上空间步长,Δt代表时间步长;向量形式的控制方程即上述公式(11)-(12)的具体离散过程如下:
1)直接将连续性方程即公式(1)和动量方程即公式(2)通过交错格式中心差分法离散即
将离散方程组即公式(24)进行Van-Lee型插值,并将离散方程组在时空重区间进行重积分:
上式(25)中的其为tk时刻离散变量交错格式的单元平均,且将单元平均替代后得到:
对流项即方程(26)的左侧形式采用:
2)将交错格式离散后方程组即公式(24)转化为非交错格式离散方程:
上述公式(28)迭代特定时间步长后,输出结果包括:
特定精度网格下网格各空间节点处得到多相泥石流深的空间分布即
hm(i,j)=hs(i,j)+hf(i,j); (29)
多相泥石流速度的空间分布:
多相泥石流颗粒体积浓度的空间分布即
多相泥石流密度的空间分布即ρm=ρsCsf(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的函数。
8.如权利要求7所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,将多相控制方程即公式(5)-(10),通过离散方程即公式(24)-(28)离散求解,此离散方程即公式(24)-(28)构建了变量在不同时间层的迭代关系,采用可变时间步长进行变量迭代,初始时间步为0.01s,可变时间步满足:
其中,上式(32)中,ρ(A)为谱半径,Cr代表克拉数;为满足离散方程的稳定性,Cr满足:
9.如权利要求8所述的固液多相适用于泥石流的动力学数值模拟方法,其特征在于,求解上述离散方程即公式(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)),体积浓度空间分布其中,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)及泥石流多相平均速度其中,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及更新可变时间步 其中,t2为下个时间层的时间,t1为上个时间层的时间,Δt是时间步,Cr为克拉数方程,min为一个函数且值是取得两个参数中的较小值;ρ(A)为A的谱半径,为每个时间步的最大速度;按如上顺序继续迭代变量至当前时间大于等于计算总时间,迭代结束;
(5)迭代特定时间步长后,输出结果为:
特定精度网格下网格各离散点固液两相体积占比hs,hf及多相泥石流密度ρm(i,j)=ρsCs(i,j)+ρf(1-Cs(i,j)),体积浓度空间分布x,y方向上速度变化us(i,j),vs(i,j),uf(i,j),vf(i,j)及泥石流多相平均速度其中,us,vs为固相x,y方向上的速度分量;uf,vf为液相x,y方向上的分量;vm代表多相泥石流速度为固相速度vs和液相速度vf的一个函数;Cs(i,j)为固相颗粒浓度,且和每个矩形网格中的固相液相占比有关系,为固相深度hs(i,j)和液相深度hf(i,j)的一个函数(31);ρs,ρf为固相液相密度且为定值。
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 true CN109657322A (zh) 2019-04-19
CN109657322B 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)

Cited By (8)

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

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070219725A1 (en) * 2004-09-10 2007-09-20 Tao Sun Method For Evaluating Sedimentary Basin Properties By Numerical Modeling Of Sedimentation Processes
CN107506566A (zh) * 2017-10-16 2017-12-22 中国科学院、水利部成都山地灾害与环境研究所 一种新型泥石流动力学数值模拟分析方法及系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070219725A1 (en) * 2004-09-10 2007-09-20 Tao Sun Method For Evaluating Sedimentary Basin Properties By Numerical Modeling Of Sedimentation Processes
CN107506566A (zh) * 2017-10-16 2017-12-22 中国科学院、水利部成都山地灾害与环境研究所 一种新型泥石流动力学数值模拟分析方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
乐茂华 等: "基于固液两相流模型的泥石流流速垂向分布研究", 《水利学报》 *

Cited By (13)

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

Also Published As

Publication number Publication date
CN109657322B (zh) 2020-08-21

Similar Documents

Publication Publication Date Title
CN109657322A (zh) 一种固液多相适用于泥石流的动力学数值模拟方法
CN106529198B (zh) 一种泥石流全过程数值模拟及数值计算方法
Beyers et al. Modeling transient snowdrift development around complex three-dimensional structures
Lin et al. Evolution of the large landslide induced by Typhoon Morakot: A case study in the Butangbunasi River, southern Taiwan using the discrete element method
CN106599457B (zh) 一种基于Godunov格式一、二维耦合技术的山洪数值模拟方法
Wells et al. An empirical investigation of gully widening rates in upland concentrated flows
Oñate et al. Advances in the particle finite element method for the analysis of fluid–multibody interaction and bed erosion in free surface flows
CN107506566A (zh) 一种新型泥石流动力学数值模拟分析方法及系统
CN108446502B (zh) 一种利用完整二维浅水方程组获得流域单位线的方法
Zhang Numerical simulation of debris flow runout using Ramms: a case study of Luzhuang Gully in China
Zhang et al. A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid
CN106959261A (zh) 一种预测沉积物颗粒分布与配比的方法
Sawadogo et al. Physical and coupled fully three-dimensional numerical modeling of pressurized bottom outlet flushing processes in reservoirs
Fiedler et al. Hydrologic response of grasslands: Effects of grazing, interactive infiltration, and scale
Sanz-Ramos et al. Numerical modelling of dense snow avalanches with a well-balanced scheme based on the 2D shallow water equations
CN112733412B (zh) 一种水动力作用滑坡运动机制研究的速度等效表征方法
Elfeki et al. Application of the random walk theory for simulation of flood hazards: Jeddah flood 25 November 2009
CN113240803B (zh) 一种降雨诱发边坡地质灾害场景模拟分析方法
Calantoni et al. Simulation of sediment motions using a discrete particle model in the inner surf and swash-zones
Meyrat et al. Voellmy-type mixture rheologies for dilatant, two-layer debris flow models
Garcia‐Sanchez et al. Lattice–gas approach to surface runoff after rain
Cuomo et al. Insights from LS-RAPID modeling of Montaguto earthflow (Italy)
LIANG et al. Simulation of debris flow impacting bridge pier tests based on smooth particle hydromechanics method
Daly et al. Estimating forces on an ice control structure using DEM
Vollmöller A shock-capturing wave-propagation method for dry and saturated granular flows

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

Granted publication date: 20200821

CF01 Termination of patent right due to non-payment of annual fee