堰塞坝坝体溃决时空演化的数值测算方法

1.本发明属于山地灾害防治、环境治理技术领域，具体涉及一种堰塞坝坝体溃决时空演化的数值测算方法。


背景技术：

2.岩石崩塌、滑坡和泥石流等大规模地质灾害可导致河道局部堵塞及堰塞坝的形成。堰塞坝的形成和破坏是世界范围山区常见的地球物理地貌过程。堰塞坝通常由松散、分选差的土壤组成，这些物质极易发生溃决。溃坝引发的突发洪水对居民、生计和下游当地基础设施构成持续威胁。这些特大洪水是短暂的，通常比稳定河流的流量大几个数量级。洪水并且可以极大地改变河道沿线的地形，是山区地貌演化的主要驱动力。
3.堰塞坝的溃决机制主要是漫顶溃决，即堰塞坝体上游储存的水位缓慢上升至坝顶处并持续上升，导致水流从坝体顶部形成溢流，溢流逐渐冲刷堰塞坝体形成溃口并不断变宽变深的过程，故此溃坝的核心在于溃坝表面溢流与坝体物质的动力学演化。在漫顶溃决过程中，水的势能逐渐转化为动能顺流而下。随着流速的增加，导致土壤侵蚀加剧并使水体携带更多的固体物质导致密度发生变化，而溢流密度体积等性质的变化又会形成正反馈影响坡面侵蚀速率的大小，加剧溃口的空间演变。随着蓄水量的下降。溃口尺寸持续扩大，纵向形成侵蚀面，水位持续降低，最终上游储水会停止流动，形成最终的溃口形状及纵向侵蚀面。
4.堰塞坝坝体溃坝可因为巨大的出流量转化为危害严重的洪水，对下游地区的基础设施及人员构成严重威胁。因此，全面了解溃坝机理对于风险评估和减灾至关重要。
5.现有技术中对堰塞坝溃决关键参数(溢流量，溃口演化)的测算大多依托传统经验或半经验公式，如专利申请号为202011049901.1的文献中提出以坝体中值粒径作为指标，计算堰塞坝坝体的可蚀性，再通过经验公式计算固定形状的溃口演变及下游洪水c峰值流量，虽然该方法基于一定的物理意义，但是无法反映出溃坝洪水及土体各自及相互的力学变化，动力学演化，且堰塞坝溃口形状应是由表面溢流和坝体性质多方面决定的，该文献中采用固定形状溃口缺乏物理依据。
6.在另一专利申请号为201910803602.3的文献中公开了建立了堰塞坝储水量、放大系数、及出流量的经验关系预测泥石流堰塞坝溃决后的冲出量，但是该方法具有较强经验性缺乏物理意义，该方法也无法体现不同形状、性质坝体带来的特殊性，同时溃口形状也采用固定几何形状，无法反映出堰塞坝溃决时固液相之内及之间的力学作用及动力学演变，以及无法考虑溢流洪水一直保持纯水特征即侵蚀进程中溢流与固相物质的相互掺混导致的性质变化


技术实现要素：

7.有鉴于此，本发明的目的之一在于提供一种堰塞坝坝体溃决时空演化的数值测算方法，该方法可以准确的测算堰塞坝溃决过程中各参数的变化。
8.为实现上述目的，本发明的技术方案为：一种堰塞坝坝体溃决时空演化的数值测算方法，包括以下步骤：
9.使用物理参数及几何尺寸参数建立堰塞坝模型；
10.基于堰塞坝溃决溢流为简化、一维、多相、等温的层流水土混合物特性，建立溢流模型；
11.通过深度积分假设对所述溢流模型进行简化，得到多相溢流模型；
12.获取溢流提供的坝体剪切力与坝体的抗剪切强度，判断侵蚀是否发生；
13.基于时间测算坝体所有位置的自然坡脚与临界土体内摩擦角，若当前位置自然坡脚大于临界土体内摩擦角，则当前位置的土体进入溢流。
14.进一步地，所述溢流模型基于ns方程中质量守恒及动量守恒获取。
15.进一步地，所述溢流提供的坝体剪切力通过曼宁方程计算获得。
16.进一步地，所述抗剪切强度通过所述堰塞坝模型的测量参数及剪切强度按纵向逐渐递增方式获得。
17.进一步地，通过以下方式判断侵蚀是否发生：
18.若所述剪切力大于抗剪切强度，则侵蚀发生；
19.反之，则侵蚀不发生。
20.进一步地，所述抗剪切强度通过以下方式得到：
21.τ
c
＝k
r
x+c
b
22.其中，k
r
为无量纲参数，x为坝体某处距溃口的纵向距离，c
b
为溃口处的抗剪强度值。
23.进一步地，所述堰塞坝模型的水土耦合计算使用水土多相模型。
24.与现有技术相比，本发明具有如下优点：
25.1.本发明中，通过建立堰塞坝模型，对模型的数值进行测量计算，可以科学准确获取堰塞坝溃决导致的侵蚀面变化，溃口时空演变，及上游蓄水水位变化，下游洪水流量、流深的变化。
26.2.本发明弥补了现有技术中侵蚀公式不明确，侵蚀缺乏物理意义等缺点，本发明中提供的堰塞坝坝体溃决测算方式具有较强的物理意义，其中参数都可以由真实条件下的堰塞坝溃决模型测算获得，具有较高的可实施性及较好的重复性和科学性。
27.3.本发明中提供的测算方法适用于相关灾害领域工作者使用，提高泥石流及其次生灾害防治，环境保护的科学化水平。
28.4.本发明中考虑堰塞坝抗剪强度在纵向方向上的分布，提出了一种坝体在纵向上抗剪强度分布的模型。
附图说明
29.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
30.图1为本发明中堰塞坝模型抗剪强度纵向分布图；
31.图2为本发明中真实情况下堰塞坝坝体侵蚀随时间变化图；
32.图3为使用堰塞坝坝体溃决时空演化的数值测算方法得到的数值模拟结果，即模拟堰塞坝坝体高度随时间变化图。
具体实施方式
33.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
34.所举实施例是为了更好地对本发明进行说明，但并不是本发明的内容仅局限于所举实施例。所以熟悉本领域的技术人员根据上述发明内容对实施方案进行非本质的改进和调整，仍属于本发明的保护范围。
35.需要注意的是，本实施例中的各代表式的下标如无特别提出解释，则仅仅用于区分，无特殊含义。
36.实施例1
37.本实施例中公开一种堰塞坝坝体溃决时空演化的数值测算方法，具体包括以下步骤：
38.s1：使用物理参数及几何尺寸参数建立堰塞坝模型；
39.本步骤中，需要模拟真实条件下堰塞坝溃决复杂多相耦合运动，首先，使用物理实验参数及坝体几何尺寸参数建立堰塞坝模型，数值实验1∶1还原物理实验中坝体参数，如坝体长宽高等参数；堰塞坝建立后按照坝体坝高设置上游蓄水，蓄水高度h
f
高于坝体h
dam
0.1m，使得上游纯水可以获得一个较大的初速度开始侵蚀堰塞坝(此设置为关键初始条件设置)，若蓄水相较堰塞坝坝体过高，可能发生洪水初速度，洪水深度过大而引起侵蚀率过大，进而使得侵蚀进度大大加快。而若初始动力条件过小，则侵蚀始终无法发生。即是否发生堰塞坝侵蚀，是由坝体可蚀性，坝体形状，上游蓄水量，地形等的综合影响。
40.s2：基于堰塞坝溃决溢流为简化、一维、多相、等温的层流水土混合物特性，建立溢流模型；
41.本步骤中，基于堰塞坝溃决中溢流为一种简化、一维、多相、等温的层流水土混合物建立溢流模型，其密度由混合物中的水和土的含量决定，由于等温假设，水流由ns方程中质量守恒及动量守恒决定，溢流模型为混合物的质量和动量守恒：
[0042][0043][0044]
其中，ρ
f
是溢流密度，ρ
f
＝c
s
ρ
s
+(1
‑
c
s
)ρ
w
，c
s
是体积含沙量，ρ
s
、ρ
w
是泥沙和纯水的密度，在一具体实施例中可以分别为2650kg/m3和1000kg/m3；u
f
是溢流速度矢量；g是重力加速度；t是应力张量，t＝
‑
pi+τ，p是流体压力，τ是剪切应力张量，i为二阶是单位张量；代表向量的散度算子，体现为表向量的散度算子，体现为其中，f＝(fx，fy，fy)为一个向量，fx，fy，fy分别为该向量在x，y，z方向的分量)；为矢量叉乘算子，其中
向量uf、vf为向量的坐标表示。
[0045]
s3：通过深度积分假设对所述溢流模型进行简化，得到多相溢流模型；
[0046]
本步骤中，设定质量守恒和无散场的不可压缩，即并且因为水流其特征长度大大超过了流动深度，故溢流模型可以通过深度积分假设进行简化，使用基于深度平均法简化后的多相溢流模型的质量守恒和动量守恒方程为：
[0047][0048][0049][0050][0051]
其中，h
f
和u
f
分别表示溢流的水深和沿x方向的流速；g
s
是沿河床的重力加速度分量；h
s
是溃决坝体的高度，是溃决坝体高度在时间上的变化率；ρ
f
是溢流密度，其值由溢流中体积含沙量c
s
决定；p为堰塞坝土体的孔隙率；t为时间；x为笛卡尔坐标x方向空间；f
bs
为流体受的底部摩阻力，k
ap
为主动土压力系数，e
b
为坝体侵蚀率，k
d
为侵蚀系数，τ为溢流提供的剪切力，τ
c
为坝体抗剪强度。
[0052]
进一步地，质量平衡方程右侧源项表示由于河床侵蚀引起系统的质量变化率，而动量平衡方程右侧源项表示由于河床侵蚀及溢流运动所引起的系统动量变化。
[0053]
s4：获取溢流提供的坝体剪切力与坝体的抗剪切强度，判断侵蚀是否发生；
[0054]
在溃坝期间，侵蚀可能会发生可能不会发生，其取决于坝体的可蚀性及上游蓄水库容，地形坡脚等因素决定，本实施例中使用一种判断侵蚀是否发送的模型：若坝体表面形成溢流且溢流提供坝体的剪切力大于坝体抗剪强度，即(τ＞τ
c
)，则表示溢流洪水携带土体的能力大于土体抵抗侵蚀的能力，则侵蚀将会发生；若洪水提供剪切力小于土体抗剪强度，则表示洪水无法提供足够大的剪切力携带堰塞坝中的固体物质，则侵蚀不会发生。其中溢流提供的剪切力τ可根据曼宁方程计算：
[0055][0056]
其中，n为曼宁系数，u为溢流流速，h为溢流流深。
[0057]
本实施例中的坝体的抗剪切强度基于真实条件下堰塞坝模型实验的实际测量及剪切强度按纵向逐渐递增的假设，故抗剪切强度通过以下方式得到：
[0058]
τ
c
＝k
r
x+c
b
[0059]
其中，k
r
为无量纲参数，x为坝体某处距溃口的纵向距离，c
b
为溃口处的抗剪强度值。
[0060]
s5：基于时间测算坝体所有位置的自然坡脚与临界土体内摩擦角，若当前位置自然坡脚大于临界土体内摩擦角，则当前位置的土体进入溢流。
[0061]
本发明中建立一个坝体崩塌模型：
[0062][0063]
其中，为当前坝体倾角，为临界倾角。
[0064]
具体地，在堰塞坝溃决数值计算中，每个时间步系统都会自动测算坝体所有位置的自然坡脚，若自然坡脚小于临界土体内摩擦角，则判断土体相对稳定不会发生崩塌效应，若自然坡脚大于土体临界坡脚(临界土体内摩擦角)，则判定土体在本时间步内处于不稳定态，强制将高于临界坡脚内的土体携带进入溢流；本发明中采用的崩塌模型是一个基于土体临界坡脚的一个再平衡方程，该崩塌模型能在一定程度上维持堰塞坝溃决时溃口及侵蚀面演化的稳定性。
[0065]
进一步地，本实施例中的测算方法中使用的水土耦合计算模型为水土多相模型，即表面溢流随堰塞坝侵蚀作用的进程而与侵蚀土体发生掺混导致溢流性质变化及质量变化，溢流的性质及质量变化会直接或间接的进一步影响溢流的表面流速u
f
，流深h
f
，提供的剪切力τ等关键力学参数，进一步间接的影响堰塞坝的侵蚀速率，溃口及侵蚀面的空间演化等。
[0066]
溢流性质变化为：
[0067]
ρ
f
＝c
s
ρ
s
+(1
‑
c
s
)ρ
w
[0068]
质量变化为：
[0069][0070]
本实施例中使用的多相模型能够时刻反映系统中水土动力学变化所带来的实时影响及相应的正反馈变化，使得固体和液体物质能通过相互作用模型(侵蚀，崩塌)耦合在一起。
[0071]
实施例2
[0072]
基于实施例1中的方法，本实施例中进行相关实验操作，具体为模拟真实条件下堰塞坝溃决复杂多相耦合运动，所有参数将取自物理实验所得，用于比较数值实验及物理时间的结果，用于揭示堰塞坝溃决中水土侵蚀物理机理及运动机理。
[0073]
首先，参考步骤s1中建立堰塞坝模型，同时，指定一定体积的蓄水量可以自动形成上游蓄水区域。同时需要设置一个水位高度参数，使之蓄水最高水位比坝高略高，因为高于坝高的部分纯水处于自由态，没有被坝体束缚，在数值计算开始后这部分纯水将在重力作用下向下游流动，以获得一个适当的初速度形成最初的溃口，本实施例中的堰塞坝模型溃决物理学参数取值可参考表1：
[0074]
表格1堰塞坝模型溃决物理学参数取值
[0075][0076][0077]
然后设置侵蚀相关系数，本实施例中基于堰塞坝溃决模型实验结果，即堰塞坝抗剪强度在纵向方向存在分布规律，得到如图1所示的堰塞坝模型抗剪强度纵向分布图，由此可知，可采用三种纵向抗剪强度分布模式，分别是线性分布及多项式分布，坝体的抗剪强度代表了其抵抗溢流剪切力的能力，间接的影响了侵蚀速率的变化。而侵蚀速率的变化会影响坝体溃口及侵蚀面的空间演化。
[0078]
本实施例中采用图1中的线性分布坝体抗剪强度进行堰塞坝模型纵向侵蚀面变化实验，得到如图2所示的堰塞坝坝体侵蚀随时间变化图，同时根据实施例1中的步骤，采取相关模型进行数值计算，得到图3所示的坝体高度随时间变化图，对比两者可知，堰塞坝坝体溃决时空演化的数值测算结果可以很好的和堰塞坝模型真实实验结果相吻合，主要表现在坝体最前端坡脚较高，侵蚀速率随纵向增大明显，坝中位置形成明显的侵蚀坑，其形成位置和侵蚀速率最大处吻合，而纵向后部形成明显的拱起，由坝体抗剪强度的纵向分布和携沙洪水密度不断增大引起。
[0079]
本发明中，依托于一系列松散土壤形成堰塞坝的坝体溃坝模型试验，以探讨堰塞坝溃坝期间形成的漫顶破坏过程，在此基础上，通过研究了溃坝过程中不同侵蚀变量(溃坝流量、溃坝流量、洪水密度和剪应力)与溃坝过程中纵向和横向侵蚀速率的关系。溃坝流量与溃坝过程中的深度及宽度侵蚀密切相关。并计算了坝内侵蚀率的空间分布，侵蚀率沿水流方向先增大后减小。根据分析数据表明，由于侵蚀作用的发生，导致携沙水流的密度上
升，而含沙量的增加会导致洪水在水流方向的冲刷能力降低，导致溃坝过程中侵蚀率降低，在纵向方向形成一定大小的冲刷坑。并且在实验数据中发现，堰塞坝的抗侵蚀能力随纵向方向逐渐增加。从而提出一套以计算堰塞坝溃决水体及土体时空演变的数值测算方法，在原有基础上进一步明确了堰塞坝溃决过程中固相液相复杂的受力特征及相互影响的正反馈机制，提高了对多流态山地灾害过程的动力过程数值模拟研究的科学性和可靠性，提高了灾害防治的针对性、增强预防效果为泥石流减灾提供技术支撑。
[0080]
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。