坝堤瞬间溃决洪水过程研究

Study on Instantaneous Outburst Flood Process of Dambreak

摘要: 坝堤溃决形成的洪水，流量大、传播速度快，其破坏性极大。本文建立了坝堤瞬间溃决与河道洪水波演进耦合的数学模型，并以泉河流域红岩一级水电站为对象，研究分析了红岩一级发生溃坝事故后的溃口洪水过程及其洪水波演进过程，预测溃决后对下游河道造成的洪水影响。结果显示，2号断面流量和红岩一级的出库流量在同一时刻达到峰值，红岩一级出库流量峰值为17,817 m³/s，2号断面的峰值流量为15,281 m³/s；经过28分钟，洪峰到达大峡水库，洪峰流量为10,190 m³/s。该模型可以较好地预测溃坝过程及下游河道洪水波演进过程，模型计算速度较快，满足应急过程中的洪水演进计算快速、高效、实时的要求。

Abstract: The flood formed by dam break has large flow and fast propagation speed, and it is very destructive. In this paper, the mathematical model of the coupling between the instantaneous dam break and the flood wave evolution of the river channel is established. Taking the Hongyan First-Level Hydropower Station in the Quanhe River basin as an example, the flood process of the breach and the flood wave evolution after the dam break accident are studied and analyzed, and the flood control effect on the downstream river channel after the dam break is predicted. The results show that the flow rate of No. 2 section and the flow rate of Hongyan first level reach the peak at the same time, the peak flow rate of Hongyan first level is 17,817 m³/s, and the peak flow rate of No. 2 section is 15,281 m³/s. After 28 minutes, the flood peak reaches the Daxia Reservoir, and the peak flow is 10,190 m³/s. The model can better predict the dam-break process and the downstream river flood wave evolution process, and has a fast calculation speed and meets the requirements of fast, efficient and real-time flood evolution in the emergency process.

1. 引言

近半个世纪以来，坝堤的溃决时有发生，产生的洪水波对坝堤下游生态与环境造成剧烈冲击，甚至造成人类生命和财产安全的巨大损失，影响当地社会经济有效运行 [1]。溃坝洪水波的运动情况十分复杂，影响因素较多，早期多是通过物理模型来进行研究。随着计算技术和数值方法的不断进步，越来越多的研究者对溃坝问题采用了数值计算的方式来模拟，数值计算的优越性也得到了越来越多的体现。

坝堤安全问题逐渐受到广泛的关注，对下游人民的生命财产安全的风险防范意识也不断加强，尽管坝堤失事的概率较小，但其会造成严重的直接或间接损失。要预测和了解坝堤失事后可能造成的影响，必须开展溃坝问题的研究与探讨，因此溃坝水力学研究具有重大的实用价值与理论价值。对于溃坝问题的研究，包含了工程水文学、泥沙冲刷与输运力学、流体力学、边坡不稳定土力学、抢险救灾措施等一系列的相关知识，在机理的认识上还不够准确和完整。现实中的原型观测受到多种条件限制，因此获取原型观测数据相当困难。而受到模型相似性和实际试验条件的约束，物理模型的建立也有较大的难度，因此当前数学模型是探讨溃坝问题的主要手段。溃坝数学模型的研究，能进一步促进溃坝机理的深入了解，能更全面地保证坝堤的安全性与稳定性，预测溃坝洪水的危害程度，为防洪决策、控制和减轻各方面的损失提供可靠的数据支撑。

2. 模型原理

2.1 瞬间溃决

瞬间全溃是溃坝的一种主要形式，也是最简单、最直接的溃坝问题，是研究坝堤逐渐或局部溃决的重要基础，而且一些中小型坝堤通常会在几分钟，甚至更短的时间内溃决、冲刷、坍塌，洪水急泻而下，所以研究瞬间全溃是十分有必要的 [2]。

瞬间溃时水波运动过程如下：溃坝初瞬坝堤上游水位急剧性陡降，形成逆向负波向上游传递，随着大量水体下泄，波形随时间逐渐平缓；坝堤下游水位快速起涨，形成顺流正波，常出现立波(不连续波)，水流湍急，常伴有复杂且强烈的泥沙运动；经过河段的槽蓄以及河道的挠阻，立波逐渐衰退直至完全消散，详见图1。瞬间溃坝的最大流量一般发生在溃决初瞬，其流量过程线大致成下凹形退水曲线，且溃坝洪水波向下游传播的过程是不断坦化的。在探讨溃坝问题时，溃口洪水的计算是研究的关键，当确定完溃口最大流量后，便可进而推求溃口洪水过程线。

Figure 1. Water wave motion process during instant dam collapse

图1. 瞬间溃时水波运动过程

2.1.1. 溃口最大流量

溃口最大流量的计算参考《水力计算手册》 [3] 以及《溃坝洪水数值模拟及其应用》 [4]，取在专业领域经过检验且得到认可的经验公式，以便更准确的反映溃口最大流量 [5]：

黄河水利委员会：

$Q_M=0.296\sqrt{g}{\left(B/b\right)}^{0.4}{\left[\left(11H_0-10h\right)/H_0\right]}^{0.3}b{H}_0^{3/2}$ (1)

式中：Q_M为溃口最大流量(m³/s)；b为溃口平均宽度(m)；H₀为溃坝前上有水深(m)；L为库区长度(m)，一般可采用坝址断面至库区上游端部库面突然缩小处的距离，当L/B > 5时，其影响不在增加，均按L/B = 5计算；B为大坝坝长(m)；h’为溃坝最终残留高度(m)；h为有效水深(m)；k为经验系数(k = 1.4(bh’/BH₀)^1/3)。

根据经验，溃口平均宽度的b值的大小直接影响溃口流量的大小，全溃时等于坝长。当溃坝的蓄水库容V ≥ 100万m³时，按 $b=K_1{V}^{1/4}{B}^{1/7}{H}_0^{1/2}$ 估计(K₁称为坝体的材质系数，对粘土类、粘土心墙或斜墙以及土、石、混凝土等K₁ = 1.19，对均质壤土K₁ = 1.98)，当V < 100万m³时，按 $b=K_2{\left(V{H}_0\right)}^{1/4}$ 估计(坝体施工和管理质量好的K₂取6.6，差的取9.1)。

当b = B，h’ = 0时，为全部溃决；当b < B，h’ = 0时，为横向局部溃决；当b < B，h’ > 0时，为横向与竖向局部溃绝同时存在。

2.1.2. 溃口流量过程线

瞬间溃流量过程线的推求，采用概化典型流量过程线法，基于详算法成果及模型试验资料的整理分析，发现瞬间溃流量过程线与最大流量Q_M、溃坝前下泄流量Q₀ (假定Q₀ = 0，同时考虑上游来水情况)及溃坝可泄库容W (可由溃坝最终残留高度及库容曲线确定)有关。其线形可概化为4次抛物线，也可概化为2.5次抛物线，即溃坝瞬间，流量陡增至Q_M，紧接着流量快速衰减，形成下凹曲线，最终趋近于原下泄流量Q₀。概化典型流量过程线的4次和2.5次抛物线分别如表1、图2所示，其中T为溃坝库容泄空时间，t为任意时刻，目前工程上多用4次抛物线来概化流量过程线。

Table 1. Generalized typical flow process table

表1. 概化典型流量过程线表

Figure 2. Flow process line at dam site

图2. 坝址流量过程线

当Q_M、Q₀及溃坝库容W已知时，就可用试算确定流量过程线，其步骤如下。

1) 根据Q_M及W初步确定泄空时间T。T由下式计算：

$T=K\left(W/Q_M\right)$ (2)

K为系数，对四次抛物线来说，K一般为4~5，对2.5次抛物线K = 3.5。

2) 根据T、Q_M、Q₀由表1或表2初步确定流量过程线。

3) 验证过程线与Q = Q₀直线间的水量是否等于溃坝库容(溃坝底部以上的库容)，如不相等，则需调整初步确定的T值，直到两者相等为止。

2.2. 河道洪水演进

将计算得到的瞬间溃流量过程线输入到河道洪水演进模块。河道洪水演进需要求解一维圣维南方程组，求解一维圣维南方程组方法多样，本研究采用有限体积法和ELM (Euler-Lagrange Method)方式求解一维圣维南方程组，该方法相比传统方法物理意义更明确、求解效率更高、算法稳定性更强的优点 [6] [7]。其求解流程见图3。

2.2.1. 控制方程

一维流体动力学控制方程是圣维南方程 [8]，如下所示：

$\{\begin{array}{l}B\frac{\partial \eta }{\partial t}+\frac{\partial Q}{\partial x}=S\\ \frac{\partial Q}{\partial t}+\frac{\partial \frac{Q^2}{A}}{\partial x}+gA\frac{\partial \eta }{\partial x}+\frac{g{Q}^2}{C^{2}AR}=S\end{array}$ (3)

式中，t为时间，x为沿河道的距离， $\eta$ 为水深，Q为流量，B为水面宽，A为过水断面面积，g为重力加速度，C为谢才系数，S为源项。

2.2.2. 方程离散与求解

在本研究中，我们使用交错网络，水位控制体中心布置在河流断面处，流量控制体布置于水位控制体之间和河道首尾水位控制体外侧 [9]。假定在流量控制体内流体的速度不变，截面积服从梯形分布，流量为流体速度与面积的乘积。为提高稳定性，本研究采用隐式格式求解数值方程 [10] [11] [12] [13]。

Figure 3. Flow chart of river flood routing solution

图3. 河道洪水演进求解流程图

对式(3)进行时间项的隐式差分、压力梯度项的中心差分，并用算子分裂法将动量方程中的对流项用拉格朗日方法的结果代替，可得：

$\{\begin{array}{l}-Q_{i-1}^{n+1}+G_i{\eta }_i^{n+1}+Q_i^{n+1}=G_i{\eta }_i^n+Q_{si}\\ -F_i{\eta }_i^{n+1}+D_i{Q}_i^{n+1}+F_i{\eta }_{i+1}^{n+1}=Q_i^n{}^{\prime }+{Q^{\prime }}_{si}\end{array}$ (4)

其中：

$D_i=1+\frac{2^{\frac{7}{3}}{n}^i{n}^{i+1}g\Delta t{Q}_i^n}{\left(A_i^n+A_{i+1}^n\right){\left(R_i^n+R_{i+1}^n\right)}^{\frac{4}{3}}}$ (5)

$F_i=g\frac{\Delta t}{\Delta {x^{\prime }}_i}\frac{A_i^n+A_{i+1}^n}{2}$ (6)

$G_i=\frac{\Delta x_i{B}_i^n}{\Delta t}$ (7)

将整个河道的方程组联立后求解即可得到各时刻河道的状态。计算完整河道所有控制体的系数后，这些系数即可组成一个稀疏矩阵。本研究采用高斯消元法求增广矩阵的解系数矩阵，即可得到结果。

3. 流域基本情况

红岩一级电站在红岩子附近，河床高程664 m。坝址控制流域面积为340.4 km²，距河口55.4 km。红岩二级电站自红岩一级电站尾水，经9 km的渠道(部分遂洞)引水至狮子口附近建电站，距河口42 km。大峡电站代表坝址初定在小峡，坝址控制流域面积482.7 km²，距河口38 km。龙滩电站自大峡电站尾水，经7.0 km渠道引水至红椿坝附近(白沙河库尾)建站。获得水头46 m，距河口27.8 km。白沙河电站，坝址位于小白沙河河口下游葡萄架附近之峡谷河段，坝址河床高程354.5 m，坝址控制流域面积为812 km²，距河口15 km，泉河流域水系图见图4。

泉河流域红岩一级工程为III等中型工程，大坝混凝土面板堆石坝，坝顶长174.35 m，坝顶宽8 m，最大坝高77.2 m。根据调洪演算成果，1000年一遇洪水位734.56 m，经计算后选定坝顶高程735.20 m，坝顶防浪墙顶高程736.50 m，最大坝底宽220.69 m，上游坝坡1:1.405，下游综合坝坡为1:1.4。水库总库容3670万m³，防洪库容670万m³，兴利库容1420万m³，死库容1580万m³，库容系数5.68%。死水位710 (0.5) m，正常高(蓄)水位730.0 m，设计洪水位为730.84 m (1%)，校核洪水位为734.56 m (0.1%)，最高发电水位729.7 m；总库容3670万m³，防洪库容670万m³，兴利库容1420万m³，死库容1580万m³，库容系数5.68%。

Figure 4. Distribution map of the Quanhe basin

图4. 泉河流域水系图

4. 结果分析

4.1. 溃口最大流量计算

假定红岩一级水电站瞬间溃决，采用上述公式(1)，进行溃口最大流量计算，并采用概化典型流量过程线法推求流量过程线。结合工程上多采用4次抛物线来概化流量过程线，此处假定溃坝前无下泄流量，且考虑上游千年一遇洪水。瞬间溃坝模型的计算参数见表2，溃口流量及洪水过程线见图5：

Table 2. Mathematical model input parameters

表2. 瞬间溃坝模型输入参数

通过模型计算得出的溃口流量过程线及洪水过程线如下：

Figure 5. Flood process line at dam breach site

图5. 溃口洪水过程线

4.2. 河道洪水演进

根据计算得到的溃口洪水过程，将其代入洪水演进模型。由于泉河流域缺乏实测的水位流量数据，因此本研究假定一组下边界条件进行洪水演进。泉河流域坡降较大，且部分河道存在干涸的情况，针对这一干湿边界条件处理的问题，本研究采用加河槽的方式进行处理，对断面加一个细长的底槽，假定初始时刻河槽有水，当流量逐渐增大后，洪水充满水槽，水位上涨，断面由干变湿。

在进行洪水演进之前，先运行一次模型，进行冷处理，上边界为恒定流，下边界为恒定水位，得到河道的初始状态，包括初始水位和初始流量。得到河道的初始状态后，上边界输入洪水流量过程，下边界为水位过程，模型计算可得到洪水演进的模拟过程。

以红岩一级到大峡水电站之间的河段为例，共计19个断面，上边界为溃口计算输出的两天的流量过程，时间间隔为一分钟。假定在前20 h为正常的出库过程，在第20 h发生了瞬间溃坝，流量从20.0 m³/s涨到峰值17,800 m³/s；下边界为恒定的水位过程，恒定水位为565 m；假定河道糙率为0.045，计算时间步长为一分钟，模拟的洪水演进结果如图6所示：

Figure 6. The evolution flood process at dam breach site

图6. 溃口洪水演进过程图

图6展示了红岩一级出库流量过程、大峡的入库流量过程及2号断面的流量过程。2号断面距离红岩一级大坝仅有962 m，且计算中时间步长为分钟，故2号断面流量和红岩一级的出库流量在同一时刻达到峰值，红岩一级出库流量峰值为17,817 m³/s，2号断面的峰值流量为15,281 m³/s；经过28分钟，洪峰到达大峡水库，洪峰流量为10,190 m³/s。

5. 结论

本文建立了坝堤瞬间溃决与河道洪水波演进耦合的数学模型，并以泉河流域红岩一级水电站为对象，研究分析了红岩一级发生溃坝事故后的溃口洪水过程及其洪水波演进过程，预测溃决后对下游河道造成的防洪影响。2号断面流量和红岩一级的出库流量在同一时刻达到峰值，红岩一级出库流量峰值为17,817 m³/s，2号断面的峰值流量为15,281 m³/s；经过28分钟，洪峰到达大峡水库，洪峰流量为10,190 m³/s。计算结果表明该模型可以较好地预测溃坝过程及下游河道洪水波演进过程，模型计算速度较快，满足应急过程的洪水演进计算快速、高效、实时的要求。其结果可用于梯级水库安全评估，为应急调度提供有力的数据支撑，减小溃坝洪水对下游区域造成的灾害影响，提高流域的防洪能力。