中国拥有水库大坝近10万座，其中土石坝占比超过93%^[1]。中国的水库大坝大多兴建于20世纪50到70年代，期间建成水库占已建水库的81.8%，且修建的水库大坝多为均质黏性土坝^[2]。受当时经济、技术条件限制等影响，中国病险水库问题突出。据水利部大坝安全管理中心统计^[3]，1954—2017年，中国共有3 539座水库大坝发生溃决，其中85%以上为均质黏性土坝，且漫顶溃决的大坝占总数的50%以上。近年来，国内外研究学者围绕均质黏性土坝的漫顶溃坝问题开展了大量的研究^[4-9]，并提出了一系列的溃坝过程数学模型，但目前的溃坝数学模型尚无法合理描述均质黏性土坝漫顶溃决过程中溃口的发展规律，因此无法准确模拟溃口的流量过程。本文基于均质黏性土坝大尺度漫顶溃决模型试验揭示的溃口形状和冲蚀形式的控制因素，建立了均质黏性土坝漫顶溃坝过程数学模型，为合理预测溃坝过程提供了科学有效的手段。

1 溃坝机理

对于均质黏性土坝漫顶溃决机理的研究，国内外学者开展了大量不同尺度的物理模型试验。其中，最具代表性的是欧盟IMPACT项目资助的溃坝试验（坝高6 m）^[10]和美国农业部Hanson等^[11]开展的溃坝试验（坝高1.5～2.3 m），以及中国南京水利科学研究院开展的世界最高（坝高9.7 m）的实体坝溃决试验^[12]和原型坝高32 m的溃坝离心模型试验^[13]。通过国内外开展的大量模型试验，揭示了均质黏性土坝的漫顶溃决机理。本节主要从漫顶溃坝过程中均质黏性土坝的坝体上下游方向、坝轴线方向及下游坡的溃口发展规律阐述其溃坝机理。

1.1 坝体上下游方向溃口发展

20世纪80年代，Ralston^[14]、Ploey^[15]在对历史上溃坝案例观测资料和模型试验数据的分析研究中发现了均质黏性土坝漫顶溃决过程中坝体上下游方向以“陡坎”冲蚀为主。“陡坎”是指河床面在高程上突降，类似于瀑布状的地貌形态（图1）^[16]。

欧美学者通过溃坝模型试验发现^{[10, 17]}，对于均质黏性土坝的漫顶溃决，下游坡面先出现初始冲坑，然后合并形成近似直立状的“陡坎”，并向上游发展。欧美学者一般都认为“陡坎”的高度与坝高大致相当。南京水利科学研究院开展了坝高9.7 m均质黏性土坝溃坝模型试验^[18]，通过试验可以看出，“陡坎”形成于下游坡中部某一位置，而不是下游坡脚处，究其原因，应为试验坝高较大所致。因此，在溃坝过程数值模拟时应关注“陡坎”的形成位置与发展过程。

1.2 坝轴线方向溃口发展

国内外均质黏性土坝漫顶溃决模型试验显示，坝顶溃口一般呈直立形或倒梯形。南京水利科学研究院^[12]开展的不同黏性（黏粒含量为11.5%～33.0%）的大比尺漫顶溃决试验进一步揭示了坝顶溃口的发展规律。

模型试验显示：对于黏粒含量较低的均质土坝，溃口一般呈倒梯形，由于溃口水流只对水面线以下的坝体进行冲蚀，水面线以上的坝料将滑落进溃口，从而使溃口边坡的坡度基本保持不变（图2）；随着冲蚀深度的增加，溃口边坡发生剪切破坏，且破坏面呈近似平面（图3）。

对于黏粒含量较高的均质土坝，溃口一般呈矩形，如前所述，溃口水流只对水面线以下的坝体进行冲蚀，由于黏性较大，水面线以上的土体呈悬空状态（图4），悬空土体长度逐渐增加而发生倾倒破坏，且破坏面呈近似平面（图5）。

1.3 坝体下游坡溃口发展

目前国内外的研究一般不单独考虑下游坡溃口的发展，而将其与坝顶溃口统一分析，认为溃口呈倒梯形、矩形或正梯形^{[5, 19-20]}。其实溃坝水流翻越坝顶后，受重力作用影响，对下游坡的冲击力明显大于坝顶处，因此下游坡的溃口宽度往往大于坝顶溃口，且下游坡溃口与坝顶溃口也不在同一平面上，不应统一考虑其形状。国内外大尺度模型试验均揭示出共同的现象，即下游坡的溃口宽度明显大于坝轴线的溃口宽度。

综上可知，在漫顶水流的作用下，均质黏性土坝的溃口在3维空间内表现出不同的发展规律。在坝体上下游方向，表现为初始“陡坎”的形成、“陡坎”溯源冲刷现象；在坝轴线方向，溃口的发展表现出表层冲蚀与边坡失稳的现象，且溃口的断面形状随着坝料黏性的增加逐渐由倒梯形变化为矩形；对于下游坡溃口，由于水流在跌落过程中势能转化为动能，冲蚀能力有所增加，下游坡的溃口宽度一般比坝轴线方向的溃口更宽。

2 溃坝数学模型

通过对国内外开展的均质黏性土坝漫顶溃决模型试验的分析研究，发现以往均质黏性土坝漫顶溃决机理的研究存在一定的偏差^[21-22]，因此有必要提出一个可合理反映均质黏性土坝漫顶溃决机理和溃坝过程的数学模型。

本文基于均质黏性土坝漫顶溃决模型试验揭示的机理，建立了一个可模拟其溃决过程的数学模型（图6），模型的具体模块如下。

2.1 溃口流量

溃坝过程是库水位动态变化的过程，且应遵循水量平衡关系：

${A_{\rm{s}}}\frac{{{\rm{d}}{{\textit{z}}_{\rm{s}}}}}{{{\rm{d}}t}} = {Q_{{\rm{in}}}} - {Q_{\rm b}}$

(1)

式中， $ A _{\rm{s}}$ 为水库库面面积， $ {\textit{z}}_{\rm{s}}$ 为水库库水位， $ t $ 为时间， $ Q_{\rm{in}} $ 为入库流量， $ Q _{\rm{b}}$ 为溃口流量。

漫顶溃决时，溃口流量可采用堰流公式计算^[19,23]：

${Q_{\rm{b}}} = 1.7{B_{\rm{b}}}{H^{1.5}} + 1.3\cot \;\alpha {H^{2.5}}$

(2)

式中： $ B _{\rm{b}}$ 为溃口底宽；H为溃口处水深， $ H ={\textit{z}} _{\rm{s}}- {\textit{z}}_{\rm{b}}$ ， $ {\textit{z}} _{\rm{b}}$ 为溃口底部高程； $ \alpha $ 为溃口边坡坡角。

2.2 “陡坎”形成与溯源冲刷

荷兰学者Visser通过研究发现^[24]，漫顶水流流经下游坝坡流速逐渐增大并趋于定值，此位置与漫顶水流深度、下游坝坡坡比等因素有关。将此流速最大的位置定义为“陡坎”形成位置，该位置对应的下游坡长度 $ l_{\rm{n}} $ （即从坝顶沿下游坡向下直至流速最大位置的长度）可表示为：

${l_{\rm{n}}} = \frac{{2.5\left( {Fr_{\rm{n}}^2 - 1} \right){d_{\rm{n}}}}}{{\tan\; \beta }}$

(3)

式中： $ \;\beta $ 为下游坝坡的坡角； $ d_{\rm{n}} $ 为下游坝坡上的水流深度，其中n表示法向方向； $F{r_{\rm{n}}}$ 为弗劳德数。 $ l_{\rm{n}} $ 的计算方法可参阅文献[25]。

初始“陡坎”形成后，在溃口水流的作用下，坝顶和下游坡均遭受冲蚀。假设“陡坎”保持直立状，并以溯源冲刷的形式向上游发展（图7）。“陡坎”向上游的移动速率速度d $ x $ /d $ t $ 可表示为^[26]：

$\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = {C_{\rm{T}}}{q^{1/3}}H_{\rm{e}}^{1/2}$

(4)

式中： $ C_{\rm{T}} $ 为溯源冲刷系数，计算方法可参阅文献[27]； $ q $ 为溃口单宽流量，可由溃口流量除以溃口平均宽度获得； $ H_{\rm{e}} $ 为“陡坎”高度（图7）。

当“陡坎”向上游发展到一定程度后，在上游库水压力和漫顶水流剪应力的共同作用下，溃口处“陡坎”上游坝体发生坍塌，破坏楔体的受力情况如图7所示。当上游楔形体坍塌后，由于溃口水头突然增大，“陡坎”冲蚀过程结束，转化为表层冲蚀。

2.3 坝顶与下游坡溃口发展

对于黏性土的冲蚀，一般采用如下的坝料冲蚀率公式计算^[28]：

$\frac{{{\rm{d}}{{\textit{z}}_{\rm{b}}}}}{{{\rm{d}}t}} = {k_{\rm{d}}}\left( {{\tau _{\rm{b}}} - {\tau _{\rm{c}}}} \right)$

(5)

式中：d $ {\textit{z}} _{\rm{b}}$ /d $ t $ 为溃口底床冲蚀率，即溃口下切速率； $ k_{\rm{d}} $ 为冲蚀系数，计算方法可参阅文献[29]； ${\tau _{\rm{b}}}$ 为水流剪应力； ${\tau _{\rm{c}}}$ 为坝料临界剪应力，可通过希尔兹曲线确定^[30]。 ${\tau _{\rm b}}$ 的计算公式如下：

${\tau _{\rm{b}}} = \frac{{{\rho _{\rm{w}}}g{n^2}Q_{\rm{b}}^2}}{{A_{\rm{w}}^2{R^{\frac{1}{3}}}}}$

(6)

式中， $ {\rho _{\rm{w}}} $ 为水的密度， $ g $ 为重力加速度， $ n $ 为溃口处的糙率系数， $ A_{\rm{w}} $ 为水流面积，R为溃口处的水力半径。

对于坝顶的溃口，无论黏性如何，溃口顶宽与底宽的扩展与底床冲蚀存在式（7）～（8）的关系（图8），只是对于黏性较大的坝体，溃口边坡坡角 $ \;\beta $ =90°（图8（b））。

$\Delta {B_{\rm{t}}} = \frac{{{n_{{\rm{loc}}}}\Delta {{\textit{z}}_{\rm{b}}}}}{{\sin\; \beta }}$

(7)

$\Delta {B_{\rm{b}}} = {n_{{\rm{loc}}}}\Delta {{\textit{z}}_{\rm{b}}}\left( {\frac{1}{{\sin\; \beta }} - \frac{1}{{\tan\; \beta }}} \right)$

(8)

式中：Δ $ B_{\rm{t}} $ 为溃口顶宽增量；Δ $ B _{\rm{b}}$ 为溃口底宽增量； $ n_{\rm{loc}} $ 溃口所在位置（ $ n_{\rm{loc}} $ =1表示溃口位于坝肩，溃口只能朝一个方向发展； $ n_{\rm{loc}} $ =2表示溃口位于坝顶中部，溃口可向两侧发展）； $ \Delta {{\textit{z}}_{\rm{b}}} $ 为溃口深度增量。

对于坝体下游坡的溃口，引入考虑溃口特征的修正系数计算下游坡溃口的顶宽与底宽。模型假设下游坡溃口顶宽与底宽的扩展与底床冲蚀存在以下关系（图9）：

${B_{{\rm{down}}}} = \frac{{{n_{{\rm{loc}}}}{c_{\rm{b}}}{{\textit{z}}_{\rm{b}}}}}{{\sin\; \gamma }}$

(9)

${b_{{\rm{down}}}} = {n_{{\rm{loc}}}}{{\textit{z}}_{\rm{b}}}\left[ {\max \left( {\frac{{{c_{\rm{b}}}}}{{\sin\; \gamma }},\frac{1}{{\tan\; \gamma }}} \right) - \frac{1}{{\tan\; \gamma }}} \right]$

(10)

式中： $ B_{\rm{down}} $ 为下游坡溃口顶宽； $ b_{\rm{down}} $ 为下游坡溃口底宽； $ n_{\rm{loc}} $ 为溃口位置参数（溃口位于坝体中部取2，溃口位于坝肩取1）； $ {\textit{z}}_{\rm{b}}$ 为下游坡溃口深度； $ \gamma $ 为下游坡溃口边坡坡角； $ c _{\rm{b}}$ 为修正系数， ${c_{\rm{b}}} = \min \Bigg[ {1,}\max ( {0,}$ ${1.8\dfrac{{{B_{\rm{b}}}}}{{{b_{{\rm{down}}}}}} - 0.8} ) \Bigg] $ 。

2.4 溃口边坡失稳

随着坝顶溃口深度的不断增加，溃口边坡可能会发生失稳，采用极限平衡法分析边坡稳定性（图10）。当滑动楔形体的驱动力 $ F_{\rm{d}} $ 大于抗滑力 $ F_{\rm{r}} $ 时，边坡失稳。 $ F_{\rm{d}} $ 和 $ F_{\rm{r}} $ 可分别表示为：

${F_{\rm{d}}} = {W_{\rm{s}}}\sin\; \theta = \frac{1}{2}{\gamma _{\rm{b}}}H_{\rm{s}}^2\left( {\frac{1}{{\tan\; \theta }} - \frac{1}{{\tan\; \beta }}} \right)\sin\; \theta $

(11)

$ \begin{aligned}[b] \;\;{F_{\rm{r}}} =& {W_{\rm{s}}}\cos\; \theta \tan\; \varphi + \frac{{C{H_{\rm{s}}}}}{{\sin\; \theta }} = \\ &\frac{1}{2}{\gamma _{\rm{b}}}H_{\rm{s}}^2\left( {\frac{1}{{\tan\; \theta }}\! -\! \frac{1}{{\tan\; \beta }}} \right)\cos\; \theta \tan\; \varphi \!+ \!\frac{{C{H_{\rm{s}}}}}{{\sin\; \theta }} \end{aligned}$

(12)

式中， $ W _{\rm{s}}$ 为破坏体重量， $ H _{\rm{s}}$ 为溃口边坡高度， $ C $ 为土体的黏聚力， $ \varphi $ 为土体的内摩擦角， $ {\gamma _{\rm{b}}} $ 为土体容重， $ \theta $ 为溃口边坡失稳后的坡角。

2.5 不完全溃坝与坝基冲蚀

对于黏粒含量较高的均质土坝，水流漫顶后可能会发生不完全溃坝，即溃口处存在残留坝体；对于黏粒含量较低的均质土坝，水流漫顶后可能会发生坝基冲蚀，即溃口达到坝体底部并对坝基进行冲蚀。模型假设不完全溃坝与坝基冲蚀时溃口底部呈水平层状分布（图11），结合坝体或坝基材料的物理力学性质，预设溃口底部最终高程，并保证溃坝过程中的水量平衡，从而模拟溃坝过程。

2.6 模型计算流程

模型采用按时间步长迭代的数值计算方法模拟溃口发展过程与溃口流量过程之间的耦合关系，具体求解过程如图12所示。该计算流程图给出了各模块之间的相互联系，并且输出每个时间步长溃口流量、溃口尺寸和库水位等计算结果。

3 案例分析

选择国内外3个具有翔实监测资料的大尺度均质黏性土坝漫顶溃决模型试验对本文模型进行验证。3个案例分别为欧盟IMPACT项目大尺度模型试验（Test 1）^[31-32]、美国农业部（USDA）大尺度模型试验（Test 2）^{[11, 33]}，以及南京水利科学研究院（NHRI）大尺度模型试验（Test 2）^{[12, 18]}。模型计算参数，如水库特征信息、入库流量信息、坝体形状信息、筑坝材料物理力学指标等见表1。

表2为3个溃坝案例的模型计算与实测结果，比较的模型参数主要包括溃口峰值流量（ $ Q_{\rm{p}} $ ）、溃口最终平均宽度（ $ B_{\rm{ave}} $ ）及溃口峰值流量出现时间（ $ T _{\rm{p}} $ ）。

由表2的比较可以得出：对于IMPACT#1，溃口峰值流量、溃口最终平均宽度和溃口峰值流量出现时间的相对误差均在±25%以内，其中溃口峰值流量和溃口峰值流量出现时间的误差更控制在±5%以内；对于USDA#2，溃口峰值流量、溃口最终平均宽度和溃口峰值流量出现时间的相对误差均在±25%以内，其中峰值流量出现时间的误差更控制在±5%以内；对于NHRI#2，溃口峰值流量、溃口最终平均宽度和溃口峰值流量出现时间的相对误差均在±25%以内，其中峰值流量的误差更控制在±5%以内。计算结果表明，本文模型能较好地反映溃坝输出参数。

由于3组模型试验有完整的溃坝流量过程监测数据和最终溃口尺寸的数据，故图13对溃口流量过程的计算值与实测值进行了比较，图14仅给出了溃口平均宽度发展过程的计算值。通过对比发现，溃口流量过程线的计算值与实测值基本吻合，溃口平均宽度的发展过程与溃口流量过程趋势一致，表明本文模型可较好地反映3组均质黏性土坝的漫顶溃决过程。

4 结　论

基于大尺度漫顶溃坝模型试验，对均质黏性土坝的漫顶溃决机理进行了研究，揭示了溃口的发展规律。在此基础上建立了均质黏性土坝漫顶溃决过程数学模型，该模型充分考虑了坝料的物理力学特性和漫顶水流特征，选择国内外3个具有实测数据的大尺度溃坝模型试验来验证模型的合理性。通过比较发现，本文模型可较好地反演溃口峰值流量、溃口最终平均宽度与溃口峰值流量出现时间等重要参量，且溃口流量过程线与实测值吻合良好，验证了模型的合理性。