局部逐渐溃坝机理研究及溃口水流模拟

张建国，宁培，殷康，宋振聪

（1.中国电建集团中南勘测设计研究院有限公司，湖南 长沙 410014；2.河南省洛宁抽水蓄能有限公司，河南 洛宁 471700）

1 研究背景

目前就局部逐渐溃坝模型的溃决机理研究，使用的模型主要有两种类型，一种如溃口形状随时间等指数变化的DAMBRK 模型，另一种如OSMAN 的侵蚀崩塌模型，前者主要以经验公式的方法模拟溃决形态，并且DAMBRK 软件的主要目的是为了分析洪水演进问题，溃口形态变化模拟分析方面缺乏足够的理论支持；而后者的侵蚀崩塌理论，它从应力层面出发，考虑溃口侧向和竖向的侵蚀深度，再应用极限平衡理论进行溃口处边坡的稳定分析，从理论层面上说更有说服力，但OSMAN的侵蚀崩塌模型，考虑的外界受力条件比较简化，它仅仅考虑了崩塌土体的重力和粘土的凝聚力，与实际情况相比较缺乏溃口处外界水压力以及土体内部孔隙水压力的影响。

目前就溃坝水流的模拟工作，进行的研究相对较多，主要就二维浅水波方程使用FVM、FDM等方法进行数值模拟，而且已取得了一定的成果，但实际的局部逐渐溃决形式过程中，研究相应时刻复杂溃口形状的情形，使用这些数值计算的模拟方法处理，仍然并不多见，文中考虑CFD 软件Fluent 模拟的便捷性，处理这类问题就比较有优势了。

2 侵蚀崩塌模型的改进分析

土石坝溃决冲刷过程中，溃口处因水流冲刷引起泥沙输运，进而引起溃口的不断侵蚀扩展，最终导致溃口处土体失稳崩塌，循环发展直至水流停止冲刷行为，这整个过程是集水文学、水力学、土力学、挟沙动力学的多学科问题。

采用FLDWAV 模型来计算溃口下泄流量，计算假定是溃口为梯形溃口，考虑行进流速和下游的淹没影响，是目前应用比较广泛的流量计算模型。

式（1）（2）（3）中：Qb是溃口的下泄流量；Qi是初始时刻的溃口下泄流量；cv是行进流速的修正系数；ks是淹没系数；cd是流量系数；bs、hb分别是随溃口发展的某时刻底部宽度、高程；hi是相应时刻上游水位高程；z 是梯形溃口的边坡坡度；Bd是河谷宽度；hbm是溃口发展的最终底高程；R′是淹没影响系数。

选用Meyer-Peter-Muller 公式作为溃口的泥沙输运模型，溃口的泥沙输运量为

式（7）（8）中：△H、△B 分别是指水流经过△t 的侵蚀时间后溃口底部的冲刷深度、侧向的发展宽度；LB是溃口沿水流方向上的过水长度；p为土体的孔隙率；C1为侧向侵蚀系数。

假设初始溃口为梯形，在水流的不断侵蚀冲刷下就会形成侧向展宽和竖向冲深，这样对溃口处边坡的稳定性就会产生不利影响，当侵蚀到一定情况就会发生溃口边坡崩塌，形成新的溃口，但由于水流冲刷作用继续存在，所以这种侵蚀崩塌过程会继续不断出现，直至溃口水流引起的切应力τ小于溃口土体的起动切应力τc，侵蚀停止，这个连续不断的溃坝过程也从而停止。OSMAN针对这个侵蚀崩塌过程的分析，就应用极限平衡理论，在考虑重力W和凝聚力N情况下的坝坡失稳滑动问题，但所考虑的与实际坝坡的受力情况却有相当大的差异，因为溃口处存在一定水位，所以外表面受到水压力P1的影响，同时土壤中存在孔隙水，稳定计算时孔隙水压力P2的影响也就不可忽略了。

在此应用土力学知识，对这个侵蚀崩塌过程进行改进分析。首先确定拉裂区裂缝深度y，由朗肯土压力理论近似得：

式（9）（10）中：y为裂缝拉深深度；c为土体粘性系数；Ka是朗肯理论的主动土压力系数；φ为粘土的内摩擦角。

溃口处坝坡部分浸水，按照瑞典条分法的计算思路，水下土条的重量应按饱和重度计算，同时还要考虑滑动面的孔隙水应力和溃口边坡上的水压力。现以水面下滑动土体为脱离体，其上有边坡上的水压力合力P1、孔隙水应力合力P2，这两个力的合力与在重心位置作用的孔隙水重量和土粒浮力的反作用力Gw(其大小等于水面线下滑动土体的同体积水重)相等，方向相反，所以周界上的水压力对滑动土体的影响就可用水面线以下滑动土体的浮力代替，这实际相当于水下土体重量均按浮重度计算。

滑动土体重计算（参见图1可知滑动土体分为3个区）：

图1 溃口土体稳定分析

（1）区考虑干重的土体重W1：

（2）区考虑干重的土体重W2：

式（11）（12）（13）（14）（15）中：i 为溃口边坡初始角度；β为溃口崩塌的临界角度；H 为坝顶到侵蚀后溃口底部的高度；H1为溃口水面到底部的高度；H2为坝顶到坝坡侵蚀转折点的高度；h1为坝顶到溃口处水面线的高度；γ′为坝体土体的浮重度。

根据OSMAN 研究土体崩塌滑动角的原理，在安全系数fs =1前提下，也就是在保持F抗=F滑的前提下，变换不同的滑动面与水平面的角度，倘若都假定保持临界平衡状态，那么粘滞效果是不同的，凝聚力也是不同的，但是这种假定的保持临界滑动的凝聚力肯定是小于筑坝土体实际的凝聚力，在此以滑动角β为变量，这种临界滑动假定下的凝聚力c为函数求极值，即为真实出现的滑动面。

溃口的稳定计算流程参见图2，水流侵蚀冲刷，计算某时间段△t后溃口侵蚀的变化情况，判断溃口失稳、溃口底部发展最终深度，计算结束的判断条件是溃口切应力τ小于溃口起动切应力τc，在没有达到最终条件下，都要不断的累时迭代重复计算，其中上游水位hi可根据库区水位库容关系以及△t时间段溃口下泄的流量计算获取。

图2 溃口侵蚀崩塌稳定计算流程图

整个侵蚀崩塌过程中的溃口形状变化如图3所示，(a)是初始溃口的形状，经过水流冲刷侵蚀后形成新形状的溃口(b)，溃口(b)达到了临界崩塌状况，发生崩塌，形成(c)，溃口(c)在整个过程中类似如初始溃口，在此之后，继续反复发生如(a)～(c)的过程，当达到最终深度hbm时(通常到坝底基岩处)，形成溃口(d)，此时溃口竖向冲刷就会停止，只发生侧向侵蚀，直到溃口切应力τ小于溃口起动切应力τc，侵蚀崩塌过程结束，形成最终的溃口(e)。这种溃决模式不仅仅有理论依据，而且在对实际的溃决实验观察时，发现也是相当吻合的，具体实例可参见源于香港科技大学研究土石坝溃口发展的文献附图。

图3 侵蚀崩塌过程溃口形状变化图

3 溃坝洪水波的数值模拟

3.1 一维实例验证

考察瞬间全溃问题：等宽矩形断面河道上，坝址位于400 m处，上下游水域都是400 m，溃决前上下游均为静水，水深分别为10 m和2.50 m。

针对该一维实例分别由Fluent模拟和黎曼求解，获得30 s的溃坝洪水水位曲线见图4，带点的曲线为黎曼解，另一条为软件Fluent的模拟效果，从图示可以发现Fluent的模拟效果能比较好的吻合精确解，不足之处就是在捕捉水面线时，激波处有稍许抖动，产生这种原因是黎曼解设置了减小躁动的限制器，而软件Fluent 本身却不具备该功能，但软件模拟结果整体上还是能反应真实水面线情况的，从而验证了软件Fluent模拟溃坝水流的可靠性。

图4 溃坝水流水位图（30 s）

3.2 复杂溃口溃坝实例

溃坝实例主要是应用Fluent模拟的便捷性，模拟该文上述的溃坝发展过程中多边形溃口的水流形态。在图3中，可以看出不同时刻由于侵蚀崩塌而形成的不同形状的溃口，在文中，仅以其过程中某一时刻的溃口形状为例，应用Fluent 进行分析，模拟的为一非对称溃坝问题。

从流域模型图，计算流域的长宽均为200 m，坝体位于中间位置，坝宽为5 m，溃口顶部至左岸为35 m，至右岸为95 m，上下游静水水深分别是30、5 m。溃口形态尺寸图所示的溃口形状是在初始溃口侵蚀一段时间后，考虑到侧向、竖向侵蚀形成的梯形、矩形组合形式的溃口，假定某一时刻的溃口数据为底部宽度40 m，矩形部分高15 m，梯形斜边坡度为1：1，高度为15 m。

用Fluent 模拟获取4.80 s 的溃口处的水流形态，溃口处洪水波向下游纵横双向传播，并有负波向上游传播，坝址的溃口两侧附近形成两个非对称的漩涡，这是由于溃口所处的位置偏向坝体岸坡一侧，并非对称分布，波前以间断波的形式传播，并有壅水现象，这在一维溃坝波中也可以清晰发现，它符合实际的物理现象。

4 总结

文章从溃口泄量、输沙、侵蚀、崩塌以及针对相应溃口形状的水流形态，就溃坝问题进行了系统的模拟分析，在溃口边坡稳定性分析时，首次运用替代法，使溃口边坡受力分析更接近于现实，溃口的侵蚀崩塌机理更趋于合理；对于溃口处的水流形态模拟，引进流体力学软件Fluent进行模拟，操作便捷，并且模拟结果得到相对充足的验证。但目前的研究仍然还有很多不足，如溃口的侵蚀扩展问题仍未突破二维的局限、如何解决软件模拟结果的躁动问题、如何实现溃坝洪水波、输沙、溃口侵蚀崩塌的耦合处理，都将是未来课题研究的趋势。