工程科学与技术   2019, Vol. 51 Issue (4): 37-46
堰塞坝水力驱动溃滑过程模拟方法研究
王琳1, 段庆伟2, 孙平2, 张强2, 刘立鹏2, 王乃欣3     
1. 西安理工大学 省部共建西北旱区生态水利国家重点实验室,陕西 西安 710048;
2. 中国水利水电科学研究院 岩土工程研究所,北京 100048;
3. 中水北方勘测设计研究有限责任公司,天津 300222
基金项目: 国家重点研发计划项目(2017YFC1501103);陕西省自然科学基础研究计划(2019JQ–577)
摘要: 为评估堰塞坝安全,需要分析其受水力驱动溃滑导致的溃决变化过程。传统的分析方法或没有运用极限平衡法,或认为溃滑过程中其滑裂面倾斜角度不变,或没有考虑孔隙水压力。实际中,堰塞坝水力驱动溃滑过程是一个需要由极限平衡理论解决的问题。通过采用岩土工程中圆弧形式的边坡稳定分析方法,考虑孔隙水压力影响,运用总应力法和有效应力法模拟堰塞坝不断溃滑的过程(包含溃口横向的不断扩展过程及其垂直下切过程);开发DBS–IWHR的电子表格分析堰塞坝溃滑过程,该电子表格已被耦合到溃决洪水分析电子表格DB–IWHR,用于分析其受水动力驱动溃滑导致的溃决洪水变化过程。结果表明:1)DBS–IWHR提出确定堰塞坝水力驱动溃滑过程的模拟步骤:确定特定圆弧滑面的安全系数FS;在各种可能的滑裂面中,确定与最小安全系数Fm相关的临界滑裂面;确定与Fm=1相关的坡脚失稳的临界深度;连续溃滑过程的模拟。以上过程只需手动4步即可进行溃滑过程模拟,基于VBA编制的程序具有良好的交互性,利于读者和使用者进行矫正和二次开发。2)选取唐家山案例,分析了水动力驱动溃滑过程及溃决洪峰流量过程,其洪峰为6 500 m3/s,溃口宽度为150 m;预测的溃决洪峰流量为7 610 m3/s,溃口宽度为139.6 m。预测值与实测值的误差在允许范围内,验证了改进的水力驱动溃滑过程模拟方法的可靠性。
关键词: 堰塞坝    滑裂面    有效应力法与总应力法    溃滑过程    
Simulation Method of Hydraulically Driven Breaking-sliding Process for Landslide Dam
WANG Lin1, DUAN Qingwei2, SUN Ping2, ZHANG Qiang2, LIU Lipeng2, WANG Naixin3     
1. State Key Lab. of Eco-hydraulics in Northwest Arid Region of China, Xi’an Univ. of Technol., Xi’an 710048, China;
2. Dept. of Geotechnical Eng., China Inst. of Water Resources and Hydropower Research, Beijing 100048, China;
3. BeiFang Investigation, Design and Research Co. Ltd., Tianjin 300222, China
Abstract: Valuating breach flood of landslide dam is usually through analyzing the breach process caused by hydraulically driven breaking-sliding process.The traditional analytical method does not use the limit equilibrium method, or the slope angle of the slip surface is not changed during the collapse process, or the pore water pressure is not considered. In fact, the hydraulically driven breaking-sliding process is a problem that needs to be solved by the limit equilibrium. An approach to determine the critical slip surface by analytical method with circular slip surfaces, the effective and total stress methods dealing with different dam materials, and a procedure to model the stepped failures of the breaking-sliding process due to continuous toe cutting were proposed. A spreadsheet entitled DBS–IWHR was developed to the breaking-sliding process analysis. This spreadsheet was incorporated into another spreadsheet entitled DB–IWHR for the calculation of flood hydrograph caused by hydraulically driven breaking-sliding. Research indicated that DBS–IWHR proposed a simulation step to determine the hydraulically driven breaking-sliding process of the dam, which included finding the factor of safety for a specific circular slip surfaceFS; among a variety of possible slip surfaces, determining the critical one associated with the minimum factor of safety Fm; determining the critical depth of toe cutting associated with Fm = 1; modeling the stepped breaking-sliding process. The simulation of the collapse process could be performed in four steps manually. The spreadsheets are transparent and self-explanatory, easy for practitioners to understand and check for secondary development through the web. The developed approach was tested by back analysis of the Tangjiashan landslide dam breached in 2008 with a flood peak of 6 500 m3/s, breach width of 150 m. The predicted flood peak is 7 610 m3/s and the breach width is 139.6 m. The calculated results of the final breach base level and the peak discharge were in good agreement with the field data. Further, the results were shown to verify the reliability of the improved hydraulically driven breaking-sliding process simulation method.
Key words: landslide dam    slip surface    the effective versus total stress analysis    breaking-sliding process    

堰塞坝多为松散堆积体,具有级配宽、结构不均匀、较易冲蚀的特点,93%的堰塞坝在1年内溃决,且一旦溃决,将对下游居民带来生命和财产损失[1]。1786年6月10日中国大渡河流域康定–沪定区域堰塞坝溃决,洪水流速超过2×105 m3/s,导致超过10万人死亡[2]。2000年山体滑坡堵塞易贡藏布形成易贡堰塞坝,导致约20 亿m3洪水下泄[3]。2008年5月12日汶川地震造成的大型堰塞湖共有256处[4]

为评估堰塞坝安全,需分析其水力驱动溃滑导致的溃决变化过程[5]。一旦溃口底部发生削切,溃口两侧的边坡将发生坍塌,溃口进一步变宽,其受水力因素影响,此为堰塞坝的水力驱动溃滑过程。溃滑过程中,不断发生削切、坍塌的循环过程,直至库水泄完停止。为模拟堰塞坝的溃滑过程,现在大多数的方法均采用直线型的滑楔体模式,如Breach模型[6]、Beed模型[7]、Osman模型[8]等均属于此类。现有的溃决变化过程分析通常结合水力学和岩土工程方法,溃口流量采用宽顶堰公式。

基于对唐家山堰塞坝溃决时现场监测资料的深入分析,作者团队[9]提出改进的溃决变化过程分析方法。改进后的溃决分析方法对输入参数不敏感,可使用简单的电子表格(即DB–IWHR)进行模拟计算。由于篇幅限制,早期主要关注水力学方法的改进。作者将详细介绍溃决过程中改进的水力驱动溃滑过程模拟方法信息。

由于堰塞坝的溃滑过程分析包含了极其繁琐的程序(如寻找临界滑裂面和总应力分析方法),本文拟开发能够模拟溃滑过程的电子计算表格(DBS–IWHR),并实现DBS–IWHR和DB–IWHR的耦合。该方法应用于2008年唐家山堰塞坝溃滑过程的反演分析,为堰塞坝的溃决变化过程分析提供理论支持。

1 堰塞坝溃滑过程 1.1 溃滑机理

堰塞坝的溃滑过程在现场和试验室非常常见[10],并受水力驱动。在唐家山堰塞坝溃决时,溃口宽度从8 m扩大至150 m。图1显示了9.7 m高模型大坝破坏过程中的扩展过程[11]图2捕获了360 s的瞬间,如图1(b)所示,当土体失稳坍塌时,虚线显示了溃决的过程。

图1 不同历时模型大坝的破坏扩展过程 Fig. 1 Enlargement process of a model dam in different elapsed times

图2 图1(b)中溃滑详细变化过程 Fig. 2 Detailed landslide breaking-sliding process shown in Fig. 1(b)

1.2 已有的溃滑过程模拟方法探讨

现有堰塞坝溃滑过程分析方法多采用楔形体破坏模式[1219]。但该方法在细节上有所不同,一般涉及4个方面:1)是否考虑坡脚的垂直削切;2)是否执行临界滑裂面的搜索;3)是否考虑了水力驱动因素;4)是否有一个合理的程序模拟堰塞坝的不断溃滑过程。本文主要介绍代表性的溃滑过程模拟方法—Osman法。

Osman和Thorne[8]提出的方法似乎已经考虑了上述大多数细节,因此得到广泛应用。

首先计算由于水流冲刷作用而导致的溃口展宽距离。假定其临界剪应力为 $ \tau _c$ ,则溃口口门展宽为:

$\Delta w=0.022\;3\left( {\frac{{\tau - {\tau _{\rm c}}}}{\gamma }} \right){{\rm e}^{ - 1.3{\tau _{\rm c}}}}\Delta t$ (1)

式中: $ \Delta w$ $ \Delta t$ 时间内溃口坡脚因水流横向冲刷而变化的距离,m; $ \tau$ 为坡脚处的水流切应力,N/m2 $ \tau _{\rm c}$ 为土体的启动切应力,N/m2 $ \gamma $ 为土体的容重,kN/m3

此方法认为在溃口发生初次崩塌后,其失稳角度恒为 $\;\beta $ ,再采用极限平衡理论计算边坡是否发生后续失稳。假定斜面 $ABCDE$ 相对于平面 $AB$ 的角度为 $i$ 图3)。滑面通过坡脚,在滑面垂直方向沿 $ \Delta {\textit{z}}$ 方向变化,水平方向沿 $ \Delta w$ 变化, $ \Delta {\textit{z}}$ $ \Delta w$ 分别表示溃口纵向冲刷的深度变化和横向冲刷的宽度变化。随着滑面的变化,将形成新的滑面 $FGHCDE$ 。假设直线滑裂面 $HE$ 与水平线的角度为 $\; \beta $

图3 溃滑过程简图 Fig. 3 Sketch of breaking-sliding process

当滑动体 $EDCH$ 处于极限平衡状态时,临界滑动高度 $H'$ 表示为:

$\begin{aligned}[b] &\quad\quad{\left(\frac{H}{{H'}}\right)^2}(\sin\; \beta \cos\; \beta - {\cos^2}\beta \tan\;\phi') -\\ &\frac{H}{{{H^{'}}}}\frac{{2c'}}{{\gamma {H^{'}}}} + (\sin\; \beta \cos\; \beta \tan\;\phi' - {\sin^2}\beta )\frac{1}{{\tan\;i}}=0 \end{aligned} $ (2)

式中, $H$ 为初始溃口高度, $H'$ 为冲蚀后的溃口高度, $c'$ 为土体黏聚力, $ \phi'$ 为土体的摩擦角, $\gamma $ 土体容重。

临界倾角 $\;\beta $ 通过对黏聚力 $c$ 求导并等于0求得:

$ \frac{{{\rm{d}}\beta }}{{{\rm{d}}c}}=0 $ (3)

$\;\beta $ 通过式(4)确定:

$\tan\;2\beta =\frac{{{{\left( {\dfrac{H}{{H'}}} \right)}^2}\tan\;i + \tan\;\phi'}}{{1 - {\dfrac{H}{{H'}}}\tan\;i\tan\;\phi' }}$ (4)

Osman法的缺点也很明显:

1)岩土工程中运用圆弧或更广义的滑面形式描述滑裂面有很长的历史,而此方法认为溃口在溃滑时的滑裂面为直线滑裂面。

2)没有理论支持式(3)的假设,显然黏聚力 $c'$ 是输入的,不是变化的。正如Duncan[20]提到的确定临界滑裂面的常规方法应是采用非线性规划寻找具有最小安全系数的滑裂面。

3)由于溃滑过程受水力驱动,必然会存在水压力的影响,则必须考虑孔隙水压力,但此方法并没有考虑。

Darby等[21]为弥补这些缺陷,考虑了水压力的影响,从理论角度提高了计算精度,但由于其繁琐的步骤,很难应用于实际工程。

2 改进的溃滑过程模拟方法

之前的溃滑过程模拟方法由于缺乏来自岩土工程专业方面的足够支持而导致溃决结果出现偏差。例如,Zhong等[22]在对唐家山溃决过程分析时发现,当底高程由719 m变化到650 m时,溃决流量从6 835.53变化到41 170.75 m3/s。

以上的溃滑过程模拟方法均面临以下问题:1)采用楔形体模式;2)没有运用岩土工程中通常采用的圆弧形式的分析方法;3)没有考虑水力驱动因素;4)没有一个合理的程序模拟堰塞坝的溃滑过程(包含溃口横向的不断扩展过程及其垂直下切过程)。以上方法大部分认为溃口的倾角不变,没有不间断地模拟其不断溃滑过程。

基于此,做如下改进:

1)总应力分析方法和不排水抗剪强度参数。

堰塞坝在溃口水位骤降时,水无法排出,造成了无法通过分析或经验方法获取孔隙水压力数值。美国陆军工程师兵团提出可采用不排水抗剪强度参数的总应力法应用于水位快速骤降的坝体计算。

2)简化Bishop法。

岩土工程界通常采用圆弧滑裂面法进行边坡稳定分析。计算过程中,搜索不同的滑裂面,寻找安全系数 $FS$ 等于最小值 $F_{\rm m}$ 时的滑裂面即为临界滑裂面。

2.1 使用电子表格法模拟边坡稳定性分析

商业程序虽然可以获得快速准确的计算结果,但并不能查看到源代码,且程序提供者对计算过程和结果并不负法律责任;且这些程序的参数输入繁多,无法应用于每个工程师。考虑到堰塞坝的溃决变化过程需要及时预警和决策,故下一代溃滑过程数值方法应采用简单、易于使用的物理模拟方法。电子表格法刚好符合现阶段需要采用快速简单的程序实现岩土过程分析这一要求。作者团队以前在DB–IWHR的工作中提供了一个示例说明如何快速地计算和显示坝体溃决流量过程线[9]。与此类似,与模拟溃滑过程相关的岩土方面的计算工作也可在另一电子表格中执行,不需要简化条件。

使用电子表格法进行边坡稳定性分析可以追溯到Low等[23],他们使用近似形式满足陈祖煜等[24]提出的力和力矩平衡条件,从而允许在电子表格中进行简单的操作。然而滑坡面的条块重量、强度和几何形状的信息仍然需要手动输入和计算。Chen等[25]改进了原来的版本,斜坡的几何信息从具有LISP语言的AutoCad绘图自动获取,再借助Excel提供的VBA编程将信息传输到电子表格中。此版本可通过透明电子表格计算边坡的安全系数,不需要任何手动计算。

2.2 DBS–IWHR的编制 2.2.1 DBS–IWHR电子表格的结构和功能

工作表MAIN为主界面,如图4所示。

图4 程序主页面与流程图 Fig. 4 Worksheets and flow charts

MAIN分为数据输入区、滑面信息区、计算结果区3个区域。A区为数据输入区,用于输入滑坡形状、材料参数及初始的滑裂面形态等信息;B区用于输入滑面重量、水压力等信息;C区用于储存安全系数计算结果及滑面信息。除工作表MAIN,电子表格的最低横杠处一系列单独的工作表提供了计算中涉及的额外细节(图4),分别为SEARCH、FOS、SEARCH 1、SEARCH 2,其他随后说明。

岩土工程界通常选择圆弧滑裂面法开展边坡稳定分析,如简化的Bishop法[26]、美国陆军工程师团法[27]等。计算时需首先确定初始滑裂面的安全系数,本程序开发了“DBS–IWHR”工具栏,实现数据导入、导出,安全系数的求解及搜索滑裂面等工作。

在工作表MAIN中输入基本计算参数并自动计算初始圆弧滑裂面(从CAD导入,或本程序自己导出的数据文件)的各条块参数。需要注意的是,在滑面信息区,用户需要输入圆心坐标,并输入坡趾点作为圆弧滑面的起点,可由此自动求取圆弧半径,进而求得滑面与边界线的交点并作为圆弧滑面的终点,再进行条分法计算。用户可自行输入圆心坐标,也可取默认输入(1/4坝高),对安全系数的求解列于工作表中。

Bishop简化法是考虑力和力矩的条分法:

$\begin{aligned}[b]&\quad\quad FS=\\ &\frac{{\displaystyle\sum\limits_{n=1}^N {[\Delta W(1 \!-\!{r_u})\tan\;\phi '\!+\! c'\Delta x]/[\cos\; \alpha (1 + \tan\;\alpha \tan\;\phi '/FS)]} }}{{\displaystyle\sum\limits_{n=1}^N {\Delta W\sin\; \alpha } }}\end{aligned}$ (5)

式中, $FS$ 为安全系数, $\alpha $ 为条块的倾角, $u$ 为条块的孔隙水压力, $\Delta W$ 为条块的重量, $\Delta x$ 为条块的宽度, $c' $ $ \phi '$ 分别为条块的黏聚力和摩擦角, $r_u$ 为孔压比:

${r_u}=\frac{{u\Delta x}}{{\Delta W}}=\frac{u}{{{\gamma _{\rm avg}}h'}}$ (6)

式中, $h′ $ 为条块的高度, $\gamma _{\rm avg}$ 为条块的单位容重, $u$ 为条块底部的孔隙压力。

图4包含了主流程图的计算步骤:

步骤1 Subroutine INPUT初始化A区域中的输入信息。

步骤2 Subroutine FIND_SLICES运用区域A中的信息定义区域B中圆弧每个土条的 $x$ $y$ $x_c$ $y_c$ $ \Delta x$ $\alpha $ 等信息,并在C区域显示。此过程包括一系列的数学计算,例如计算每个圆弧和滑面连接点处的上部和下部土条的坐标点等数据。输出数据是在区域B中从17行开始的第K列到第Q列之间的单元格范围内反映。

步骤3 Subroutine DATA_UPDATING刷新第S列至第AA列中的单元格,这些单元信息包含式(5)和(6)中定义的 $u$ $h'$ $\gamma _{\rm avg}$ $\Delta W$ $c'$ $r_u$ $\phi'$ 。其中,有的变量是基于先前的计算数据,有的是利用存储在区域A中的可用信息通过编制的VBA编码直接计算的数据。命令“Application.ScreenUpdating”激活并更新这些单元格。

步骤4 Subroutine FIND_FOS执行计算式(5),包含必要中间值的计算及使用VBA动态编程寻找方程中安全系数 $FS$ 的宏。安全系数的求解值存储在第15行、第S列的单元格中,计算细节显示在工作表FOS中,用户可在此工作表中检查求解安全系数的每个细节。

2.2.2 自动搜索临界滑裂面

各种优化方法都可确定最小安全系数,DBS–IWHR使用了最简单的两步网格搜索方法。计算时,首先确定滑裂面的抗滑稳定安全系数,不断重复上述计算步骤,以寻找对应最低安全系数的临界滑裂面。寻找最低安全系数的优化方法,已有很多研究成果[28],本程序选用枚举法求解,并寻找安全系数等于1时的坡趾下切深度。

程序在Search1和Search2工作表中进行了两次搜索,形成了搜索点阵,用于存储 $n$ 个滑面安全系数的信息,其中安全系数最小的即为临界滑裂面。此处 $n$ 可自行选择,程序选择14个滑裂面为基础,实现过程如图5(a)所示。Search1中,在寻找临界滑裂面时,选择圆心的横坐标为自变量,安全系数为因变量作图,最小安全系数对应的滑裂面即为临界滑裂面,如图5(a)所示。经过一次寻找,无法确定其为最小的临界滑裂面,再将上述临界滑裂面信息输入到Search2中。搜索此滑裂面与边界线的交点,连接此滑裂面坡趾点与边界线交点,称为弦线;以此弦线上中垂线的圆心坐标为自变量,寻找其安全系数最小的滑裂面信息,临界滑裂面如图5(b)所示。需要说明的是,若经过两次搜索仍没有寻找到临界滑裂面,可将Search2搜索到的滑裂面信息输入MAIN工作表中的圆弧信息区再次搜索,直至搜索到其最小滑裂面停止。第2次搜索基于第1次搜索,每次仅需很少的工作即可寻找到足够准确的 $F_{\rm m}$

图5 临界滑动面搜索过程 Fig. 5 Process of searching for the critical slip surface

图6中,原始滑坡轮廓由 $ABCEDG$ 表示。当由 $D$ 点表示的坝趾从高程2 261.41 m垂直切割 $\Delta {\textit{z}}$ =5 m,至2 256.41 m时,滑面向具有相同幅度5 m的坡度水平切割到新位置点 $E$ ,轮廓 $ABCEFH$ 的滑面安全系数为0.998。如果发生滑坡,由此得到的表面轮廓由 $ABCJFH$ 表示。此类分析中,圆弧滑面必须通过斜坡坡趾。

图6 堰塞坝的第1次溃滑过程 Fig. 6 First breaking-sliding process for landslide dam

第2次突然扩展是由于连续的溃滑而发生的,可使用与第1次相同的程序模拟计算。再对所剩的堰塞坝继续进行分析,该坝体的一部分由第1次滑坡组成。DBS–IWHR具有自动将从第1次滑动获得的斜面轮廓ABCJFH直接更新到第2次滑动滑面的功能,使得可几乎自动地执行步骤2的横向扩展和下切过程。第2次扩展,新的坡趾深度为0.9 m, $F_{\rm m}$ 为1.010,如图7所示。

图7 堰塞坝的第2次溃滑过程 Fig. 7 Second breaking-sliding process for landslide dam

所选某堰塞坝需计算19步,以模拟从2 256.41到2 210.00 m的溃口剖面不断扩展和下切过程,如图8所示。

图8 堰塞坝的连续19步溃滑过程 Fig. 8 Cross-sections of the 19 breaking-sliding process for landslide dam

2.2.3 总应力法和有效应力法

近代土力学认为影响土体强度的因素是有效应力,而不是总应力。一般通过理论分析及试验确定土体内的孔隙水压力,再利用Mohr–Coulomb强度准则进行边坡稳定分析。但在溃决过程中,水位迅速下降,故应在不排水条件下分析堰塞坝的水动力驱动溃滑过程,类似于分析坝体下游的库水位骤降[29]情况。

DBS–IWHR考虑了溃决分析中经常遇到的两种情况:

1)第1种情况是坝体没有被水渗漏且没有地下水的作用。这种情况涉及到面板坝和堰塞坝。在模拟此类型坝体时选择式(5),与有效应力方法一起使用,式(5)中输入 $r_u$ =0或用户指定的某个其他值。

2)第2种情况是坝体完全饱和,并且在坝体溃决之前坝体内存在浸润线。这种情况涉及由诸如黏土和尾矿渣的不透水材料制成的坝体。

有效应力法[26]和总应力法[30]在本方法中均可用以分析模拟。

由于很难确定孔隙水压力系数 $A$ $B$ [30],故使用有效应力方法模拟溃滑过程较为困难。目前,不排水强度参数的总应力方法已被广泛应用。20世纪80年代,陈祖煜已经就总应力法的概念、强度指标及库水位骤降期的计算步骤进行了介绍[31],并介绍了美国陆军工程师团[27]计算手册中记录良好的计算程序案例。DB–IWHR使用该过程的简化版本,如下所示:

1)快速骤降前确定条块底部的有效应力如下:

${\sigma '_{\rm c}}=\frac{{\Delta W'}}{{\Delta x}}\cos\; \alpha $ (7)

式中, $\Delta W '$ 为条块的有效应力,由浸润线以上和以下土体部分的不饱和部分存在的浮托力确定。

2)通过式(8)确定土体不饱和强度:

$c'={S\!_{\rm u}}={c\!_{\rm cu}} + {\sigma '_{\rm c}}\tan\;{\phi _{{\rm cu}}}$ (8)

式中, $c_{\rm cu}$ $\phi_{\rm cu}$ 为不排水强度参数。

3)基于条件:

$\phi '=0$ (9)

式(5)变为:

$FS=\frac{{\displaystyle\sum\limits_{n=1}^N {{S\!_{\rm{u}}}\Delta x {\rm sec}\;\alpha } }}{{\displaystyle\sum\limits_{n=1}^N {\Delta W\sin\; \alpha } }}$ (10)

在模拟堰塞坝溃滑过程时,认为坝体的浸润线不变。

2.3 DBS–IWHR和DB–IWHR的耦合

DB–IWHR是作者团队自行编制的溃决洪水分析程序。通过溃口的入流量(通过宽顶堰方程计算)与单位时间库容水位之间的平衡建立如下方程:

$Q=CB{(H - {\textit{z}})^{3/2}}=\frac{{\Delta W}}{{\Delta H}}\frac{{\Delta H}}{{\Delta t}} + q$ (11)

式中, $C$ 取1.43~1.69 m1/2/s, $H$ 为库水位, $q$ 为自然入流量。

描述冲蚀的剪应力可使用Manning方程,计算公式为:

$ \tau=\gamma RJ=\frac{\gamma {{n}^{2}}V_{{}}^{2}}{{R}^{1/3}}\approx \frac{\gamma {{n}^{2}}V_{{}}^{2}}{{h}^{1/3}} $ (12)

式中: $ \gamma $ 为水的密度; $n$ 为粗糙度的系数(取0.025 m–1/3·s); $J$ 为坡降;当渠道宽度 $B$ 大于平均流动深度 $h$ 时, $R$ 可以近似为 $h$ 的水力半径。

如果水流没有被尾水淹没,水位通常在坝体的入口处下降。因此,有必要确定流动深度 $h$ 。假设通过溃口渠道的水流稳定,Fread使用Manning方程运用式(13)确定溃口的水深 $h$

$h={\left( {\frac{{nQ}}{{B{J^{0.5}}}}} \right)^{0.6}}$ (13)

将式(11)代入式(13)可得:

$h={\left( {\frac{{nQ}}{{B{J^{0.5}}}}} \right)^{0.6}}={\left( {\frac{{nC{{(H - {\textit{z}})}^{1.5}}}}{{{J^{0.5}}}}} \right)^{0.6}}=\frac{{{n^{1.5}}{C^{1.5}}{{(H - {\textit{z}})}^{0.9}}}}{{{J^{0.3}}}}$ (14)

式中, $H-{\textit{z}}$ 的指数为0.9,可将 $H-{\textit{z}}$ 的指数近似取1,式(14)即简化为:

$h=m(H - {\textit{z}})$ (15)

定义 $m$ 为水位跌落系数,则

$m=\frac{h}{{H - {\textit{z}}}} \approx \frac{{{n^{0.6}}{C^{0.6}}}}{{{J^{0.3}}}}$ (16)

DB–IWHR建议可使用此经验系数 $m$ 值计算溃口的水深 $h$ $m$ 具有明确的物理意义,其值为0.8和0.5,可提高结果的稳定性。

如果 $H$ ${\textit{z}}$ $h$ 都确定了, $m$ 可直接计算得出。流速计算可采用式(17):

$V=\frac{{C{{(H - {\textit{z}})}^{3/2}}}}{h}$ (17)

溃口受水力驱动,底部被不断刷深,两侧边坡发生溃滑,侧面不断扩大,其破坏形式很难统一,因此溃滑过程模拟是一个难以准确实现的过程。以前的分析方法均采用楔形体法,计算过程中采用岩土工程中已经被广泛接受的滑裂面分析方法—简化的Bishop法,如式(5)所示。如图8所示,将溃口不断失稳的数据即每一步的圆心坐标、半径和渠底坐标输入到DB–IWHR中,如此以来,溃决过程分析变得冗长。由于溃决过程发展迅速,张利民指出34%的坝将会在一天内溃决,预警和决策制定的时间十分有限,因此简便并易于应用的程序是十分必要和迫切的。另一方面,使用线性内插近似模拟中间步骤,计算精度并没有太大损失。本文提出如下简化:

只需输入第一个和最后一个滑面的信息,中间断面按连续渐变内插,同时将圆弧滑裂面近似为直线,如图9所示。初始圆弧滑裂面 $MN$ 可近似为直线 $NP$ $NP$ 为弦线 $MN$ 和切线 $NQ$ 的角平分线。

图9 溃滑过程简化 Fig. 9 Simplification of breaking-sliding process

在DB–IWHR计算时,对具有河床高度 ${\textit{z}}$ 和水深度 $h$ 的溃决渠道,基于其初始值 ${\textit{z}}_0$ $h_0$ $B_0$ 及假定高度 ${\textit{z}}_{\rm end}$ 的插值确定高程( ${\textit{z}}+h$ )处的水面 $B$ 的宽度。相应的宽度 $B_{\rm end}$ 由式(18)确定:

${B_{{\rm{end}}}}={B_0} + 2({{\textit{z}}_0} - {{\textit{z}}_{{\rm{end}}}})$ (18)

溃口边坡的角度可通过式(19)确定:

$\beta ={\beta _0} + \frac{{{{\textit{z}}_0} - {\textit{z}}}}{{{{\textit{z}}_0} - {{\textit{z}}_{{\rm{end}}}}}}({\beta _{{\rm{end}}}} - {\beta _0})$ (19)

则水面宽度 $B$ 为:

$B={B_0} + \frac{{{{\textit{z}}_0} - {\textit{z}}}}{{{{\textit{z}}_0} - {{\textit{z}}_{{\rm{end}}}}}}({B_{{\rm{end}}}} - {B_0}) + 2h\tan \left(\beta - \frac{{\text{π}} }{2}\right)$ (20)

在这个步骤中, $h$ 的值未知,但可使用前一步骤中的 $h$ 值,也可满足精度要求。

简化带来了很大的便利性,只有有限的计算精度损失。例如,在唐家山案例DB–IWHR 2014的旧版本和新版本中的计算溃决洪峰分别为7 610.0和7 71.6 m3/s。需要强调的是,DB–IWHR可由用户自行决定使用原有圆弧滑裂面或简化的梯形横截面。

3 唐家山堰塞坝案例 3.1 基本概况

2008年5月12日,四川省北川县城上游4 km的唐家山由于受到汶川8.0级大地震影响,发生大规模高速滑坡。山体冲向左岸掩埋了村庄,高速滑坡形成了唐家山堰塞坝。运用先进的遥感、现场监测等手段推测唐家山堰塞坝长803.40 m,宽611.80 m,高82.65 m,体积为2.437×107 m3。堰塞坝左侧高程为785 m,右侧高程为755 m,坝高约为90~120 m,主要巨石、石块组成,结构比较密实。

文献[9]计算了唐家山的溃口扩展过程到727 m停止。当唐家山堰塞坝发生时,水利部水文局、四川省水文局对整个溃决过程进行了实时监测,根据监测数据,徐文杰[32]采用两种不同的方法(结果分别用Q1、Q2表示)进行了分析,认为溃决流量最终在6 800 m3/s左右,如图10所示。

图10 采用不同方法计算得到的溃决流量过程曲线 Fig. 10 Flood discharge versus time determined by different methods

图10可知:堰塞坝泄洪从6月9日17:30点开始,溃口过流量实测约为93.3 m3/s;到6月10日06:30左右,溃口内的水流速度增大(约2.4 m/s),此时堰塞坝内水位高程为742.1 m;直至6月10日11:30时最大流量约6 800 m3/s;至6月10日20:00左右后溃口内的流速(约2.4 m/s)及流量基本保持不变,表明堰塞坝溃决基本结束,此时堰塞坝内水位高程为720.3 m,整个过程中唐家山堰塞坝的下泄水量约1.53 亿m3。当日下午06:00,高程为740 m时,溃口的宽度开始扩大,宽度为145 m,然后保持不变;河床为720 m高程时,溃口将不再扩大[33]

3.2 溃滑过程模拟方法验证

采用唐家山堰塞坝实例详细说明堰塞坝的水动力驱动溃滑过程。开挖的底高程为740 m,边坡坡率为1.2∶1,土体参数 $c_{\rm cu}$ 为25 kPa, $\phi _{\rm cu}$ 为22°,选择改进的溃滑过程模拟方法。对比商业程序STAB和DBS–IWHR的计算结果。由表12可知,本程序可以得到与商用程序STAB完全一致的结果。

表1 安全系数计算结果对比 Tab. 1 Comparison of safety factor calculation results

表2 搜索临界滑面计算结果 Tab. 2 Calculation results of searching for critical slip surface

选择将溃口的宽度计算到坝趾高程690 m时停止,因为此时大部分湖水在该高程已经泄洪完毕。根据以上方法运用DBS–IWHR程序模拟唐家山溃决过程,如表3图11所示。

表3 DBS–IWHR计算的每一步溃滑过程的横截面详细几何信息 Tab. 3 Cross section geometry information of each breaking-sliding process calculated by DBS–IWHR

图11 唐家山的溃滑过程 Fig. 11 Process of breaking-sliding of Tangjiashan

3.3 溃决洪水分析

将唐家山堰塞坝的实际参数输入到改进的DB–IWHR中,反演其溃决过程,得到其溃决后溃决洪峰流量及溃口宽度(图12)。预测的溃决洪峰流量为7 610 m3/s,溃口宽度为139.6 m(表4),预测值与实测值的误差在允许范围内,验证了改进的水力驱动溃滑过程模拟方法的可靠性。

图12 模拟唐家山堰塞坝溃决过程 Fig. 12 Simulation on the process of Tangjiashan landslide dam breach

表4 计算结果特征数据 Tab. 4 Characteristic data of calculated results

4 结 论

1)堰塞坝的水力驱动溃滑过程需运用力学极限平衡理论分析。在溃滑模拟方法中,考虑库水位骤降时的孔隙水压力,并针对不同坝体材料,运用总应力法和有效应力法模拟其横向的不断扩展过程及其垂直下切过程。运用岩土工程中常用圆弧形式的边坡稳定分析方法搜索临界滑裂面,从而建立新的溃滑过程模拟方法并基于VBA编制了DBS–IWHR程序。

2)DBS–IWHR提出确定堰塞坝水力驱动溃滑过程的模拟步骤:确定特定圆弧滑面的安全系数 $FS$ ;在各种可能的滑裂面中,确定与最小安全系数 $F_{\rm m}$ 相关的临界滑裂面;确定与 $F_{\rm m}$ =1相关的坡脚失稳的临界深度;连续溃滑过程的模拟。只需手动4步即可进行溃滑过程模拟,基于VBA编制的程序具有良好的交互性,利于读者和使用者进行矫正和二次开发。

3)选用唐家山案例,分析了水力驱动溃滑过程及溃决洪峰流量过程。进一步证明了该方法的适用性,其洪峰为6 500 m3/s,预测的溃决洪峰流量为7 610 m3/s,溃口宽度为139.6 m,预测值与实测值的误差在允许范围内,验证了改进的溃滑过程模拟方法的可靠性。

参考文献
[1]
Ermini L,Casagli N. Prediction of the behaviour of landslide dams using a geomorphological dimensionless index[J]. Earth Surface Processes and Landforms, 2003, 28(1): 31-47. DOI:10.1002/(ISSN)1096-9837
[2]
Dai F C,Lee C F,Deng J H,et al. The 1786 earthquake-triggered landslide dam and subsequent dam-break flood on the Dadu River,southwestern China[J]. Geomorphology, 2005, 65(3): 205-221.
[3]
Wang Lin,Chen Zuyu,Wang Naixin,et al. Modeling lateral enlargement in dam breaches using slope stability analysis based on circular slip mode[J]. Engineering Geology, 2016, 209: 70-81. DOI:10.1016/j.enggeo.2016.04.027
[4]
崔鹏,韩用顺,陈晓清. 汶川地震堰塞湖分布规律与风险评估[J]. 四川大学学报(工程科学版), 2009, 41(3): 35-42.
[5]
Wu W M,Altinakar M S,Song C R,et al. Earthen embankment breaching[J]. Journal of Hydraulic Engineering, 2011, 137(12): 1549-1564. DOI:10.1061/(ASCE)HY.1943-7900.0000498
[6]
Brown R J,Rogers D C.BRDAM users’manual[R].Denver:Department of the Interior,1981.
[7]
Singh V P,Scarlatos P D,Collins J G,et al. Breach erosion of earthfill dams (BEED) model[J]. Natural Hazards, 1988, 1(2): 161-180. DOI:10.1007/BF00126613
[8]
Thorne C R,Osman A M. Riverbank stability analysis.Ⅱ:Applications[J]. Journal of Hydraulic Engineering, 1988, 114(2): 151-172. DOI:10.1061/(ASCE)0733-9429(1988)114:2(151)
[9]
Chen Zuyu,Ma Liqiu,Yu Su,et al. Back analysis of the draining process of the Tangjiashan barrier lake[J]. Journal of Hydraulic Engineering, 2015, 141(4): 682-688.
[10]
杨阳,曹叔尤. 堰塞坝漫顶溃决与演变水槽试验指标初探[J]. 四川大学学报(工程科学版), 2015, 47(2): 1-7.
[11]
Zhang Jianyun,Li Yun,Xuan Guoxiang,et al. Overtopping breaching of cohesive homogeneous earth dam with different cohesive strength[J]. Scientia Sinica Technologica, 2009, 52(10): 3024-3029.
[12]
Zhong Qiming,Chen Shengshui,Mei Shiang,et al. Numerical simulation of landslide dam breaching due to overtopping[J]. Landslides, 2018, 15(6): 1183-1192. DOI:10.1007/s10346-017-0935-3
[13]
Morris M W,Kortenhaus A,Visser P J.Modeling breach initiation and growth[R].Wallingford:FLOODsite Consortium,2009.
[14]
Morris M W,Kortenhaus A,Visser P J,et al.Breaching processes:A state of the art review[R].Wallingford:FLOODsite Consortium,2009.
[15]
Mohamed A A A,Samuels P G,Morris M W,et al.Improving the accuracy of prediction of breach formation through embankment dams and flood embankments[C]//Proceedings International Conference on Fluvial Hydraulics (River Flow 2002).Louvain-la-Neuve,2002.
[16]
Wang Guangqian,Liu Fan,Fu Xudong,et al. Simulation of dam breach development for emergency treatment of the Tangjiashan quake lake in China[J]. Science in China Series E:Technological Sciences, 2008, 51(2): 82-94.
[17]
Chang Dongsheng,Zhang Limin. Simulation of the erosion process of landslide dams due to overtopping considering variations in soil erodibility along depth[J]. Natural Hazards and Earth System Sciences, 2010, 10(4): 933-946. DOI:10.5194/nhess-10-933-2010
[18]
Peng Ming,Zhang Limin,Chang Dongsheng,et al. Engineering risk mitigation measures for the landslide dams induced by the 2008 Wenchuan earthquake[J]. Engineering Geology, 2014, 180: 68-84. DOI:10.1016/j.enggeo.2014.03.016
[19]
Shi Zhengming,Guan Shenggong,Peng Ming,et al. Cascading breaching of the Tangjiashan landslide dam and two smaller downstream landslide dams[J]. Engineering Geology, 2015, 193: 445-458. DOI:10.1016/j.enggeo.2015.05.021
[20]
Duncan J M. State of the art:Limit equilibrium and finite-element analysis of slopes[J]. Journal of Geotechnical Engineering, 1997, 122(7): 577-596.
[21]
Darby S E,Thorne C R. Development and testing of riverbank-stability analysis[J]. Journal of Hydraulic Engineering, 1996, 122(8): 443-454. DOI:10.1061/(ASCE)0733-9429(1996)122:8(443)
[22]
Zhong Qiming,Wu Weiming. Discussion of " Back analysis of the draining process of the Tangjiashan Barrier Lake” by Zuyu Chen,Liqiu Ma,Shu Yu,Shujing Chen,Xingbo Zhou,Ping Sun,and Xu Li[J]. Journal of Hydraulic Engineering, 2016, 142(6): 07016001. DOI:10.1061/(ASCE)HY.1943-7900.0001161
[23]
Low B K,Gilbert R B,Wright S G. Slope reliability analysis using generalized method of slices[J]. Journal of Geotechnical & Geoenvironmental Engineering, 1998, 124(4): 43-50.
[24]
Chen Z Y,Morgenstern N R. Extensions to the generalized method of slices for stability analysis[J]. Canadian Geotechnical Journal, 1983, 20(1): 104-119. DOI:10.1139/t83-010
[25]
Chen Lihong,Chen Zuyu,Sun Ping.Slope stability analysis using graphic acquisitions and spreadsheets[C]//Proceedings of 10th International Conference Symposium of Landslide and Engineered Slope.Xi’an,2008:631–637.
[26]
Bishop A W,Morgenstern N. Stability coefficients for earth slopes[J]. Geotechnique, 1960, 10(4): 129-153. DOI:10.1680/geot.1960.10.4.129
[27]
US Army Corps of Engineers.Engineering and design,stability of earth and rock-fill dams:Engineer manual[R].Washington D C:Dept.of the Army,Corps of Engineers,1970.
[28]
Chen Zuyu,Shao Changming. Evaluation of minimum factor of safety in slope stability analysis[J]. Canadian Geotechnical Journal, 1988, 25(4): 735-748. DOI:10.1139/t88-084
[29]
Lowe J Ⅲ,Karafiath L.Stability of earth dams upon drawdown[C]//Proceedings of First Pan-American Conference on Soil Mechanics and Foundation Engineering,Mexico City,1960,537–552.
[30]
Skempton A W. The pore-pressure coefficients A and B[J]. Geotechnique, 1954, 4(4): 143-147. DOI:10.1680/geot.1954.4.4.143
[31]
陈祖煜. 土石坝边坡稳定分析中的总应力法[J]. 水利水电技术, 1984(10): 3-8.
[32]
Xu Wenjie.Study on landslides blocking river events and their disaster effects[D].Beijing:China Institute of Water Resources and Hydropower Research,2010.
徐文杰.滑坡堵江事件及其灾害效应研究[D].北京:中国水利水电科学研究院,2010.
[33]
Chen Zuyu,Ma Liqiu,Yu Su,et al. Closure to " Back analysis of the draining process of Tangjiashan barrier lake” by Zuyu Chen,Liqiu Ma,Shu Yu,Shujing Chen,Xingbo Zhou,Ping Sun,and Xu Li[J]. Journal of Hydraulic Engineering, 2016, 142(6): 07016002. DOI:10.1061/(ASCE)HY.1943-7900.0001162