坝体溃决过程与溃坝洪水演进耦合数值模拟

, ,

(长江科学院 水力学研究所， 武汉 430010)

1 研究背景

研究土石坝的溃决过程及其坝下游洪水演进，对于我国处置溃坝突发性洪水灾害，提升我国水利安全，具有十分重大的意义[1]。

土石坝的溃决涉及到水力学、土力学以及泥沙动力学等多学科，受坝体材料、坝体型式与尺度、水库库容及上下游水位等多种因素的影响。从20世纪60年代开始，世界各国相继开始了对溃坝问题的研究，逐步形成了如“跌坎冲刷(headcut erosion)”、“切沟侵蚀(gully erosion)”等能较为合理地揭示溃坝过程的理论，为溃坝洪水预测奠定了良好理论基础并提供了丰富的验证资料[2-5]。

由于溃坝洪水的巨大危害性，世界各国均开发了计算模型用于预测与预报，如基于溃口线性发展假设的美国国家气象局DAMBRK模型与SMPDBK模型、基于泥沙输移理论并考虑溃口边坡坍塌的Beed与Breach等模型[4]。近年来溃坝洪水平面二维数值模型的研究取得了长足的进展。Fracacrollo[6]建立了概化的室内溃坝物理模型,开展了试验研究并进行了数值模拟，溃口采用瞬溃假定；崔丹等[7]基于有限体积法与具有总变差减小特性MacCormack格式建立溃坝洪水演进数值模型，模型中溃口的发展采用了线性假定；夏军强等[8]建立了基于无结构三角网格下采用有限体积法求解的二维水动力学模型，溃口采用瞬溃假定；杨志等[9]建立了黑河金盆水库大坝溃口近区二维数值模型和下游地区溃坝洪水演进耦合数学模型，模型中溃口的发展及溃口流量采用公式计算；岳志远等[10]建立二维水沙耦合数学模型，采用修正的上扬通量公式来描述泥沙运动和滑坡堰塞体溃决过程，将溃坝过程与洪水演进耦合模拟，相比前述模型有了较大的进展，但未考虑溃口发展过程中溃口边坡的坍塌。

溃口流量过程是坝下洪水预测的核心参数，其精度不够可能导致预警过度或不足。溃口发展过程是决定溃口流量过程的关键因素，它与溃口水力要素相互耦合，具有很强的时变特性与空间分布的不均性。合理地描述坝体溃决过程中溃口“剪剥”式(或“跌坎”式)发展模式、坝面双螺旋流淘刷引起溃口横向展宽以及边坡失稳坍塌的现象，需要详细的溃口区流场随时间变化的信息，如水深分布、流速分布等。上述计算模型由于其一维性很难体现溃口发展的空间分布不均匀性，而基于水深积分的平面二维数值模型虽不能反映溃口区水流的三维特性，但能提供的水位(水深)场、平面流场及对不规则溃口形态的描述能力等，能较好地满足“跌坎”移动速度、坝体下切速度以及边坡失稳坍塌计算所需的水力参数，从而实现对溃口复杂冲蚀过程的模拟。相比一维数值模型仅采用溃口断面平均流速来计算溃口发展过程有着本质的区别，因此构建坝体溃决过程与溃坝洪水演进耦合的数值模型是溃坝洪水预测的发展趋势。此外，溃口区间断波的捕捉、洪水演进过程中引起的干湿边界的处理、模型的守恒性等亦是溃坝数学模型中需要解决的重点问题。

本文采用双曲线型冲蚀速率表达式描述坝体冲蚀、简化Bishop法搜索临界滑裂面描述溃口边坡坍塌、具有总变差不增特性的MacCormack有限体积法离散控制方程，建立了坝体溃决过程与溃坝水流完全耦合维数值模型，相比采用经验公式计算溃口发展[9]、瞬溃假定[6、8]、线性溃决假定[7]和不考虑溃口边坡坍塌[10]的模型，溃口流量过程的合理性及坝下洪水预测精度均可得到较大的提高。

2 数学模型

2.1 基本方程

(1)

(2)

将控制方程统一写成守恒型式，即

(3)

式中：U为守恒性变量，U=(h,hu,hv)；F和T分别代表对流项与扩散项；S代表源项。

2.2 溃口冲蚀

对于黏性和非黏性土体的冲蚀率与切应力之间关系，国内外学者已进行了大量的研究， Gauche基于水槽试验在2010年提出了针对非黏性土的指数形式的冲蚀速率表达式,即

(4)

(5)

(6)

式中：γ为水体重度(9.8 kN/m3)；J为水力坡度；R′为水力半径，当渠道宽深比足够大时，可近似取水深h;V为水流流速。

2.3 溃口横向展宽

溃决过程中，当侧边坡的重力以及由于渗透造成的渗透力大于土体的摩擦力和黏聚力时，溃口侧边坡将变得不稳定而坍塌。现在大多数的模型均采用直线型的滑楔体模型，如Breach模型、Beed模型、Osman模型均属于此类。本文溃口展宽采用式(7)所示的简化Bishop法搜索临界滑裂面模型。

(7)

式中：ci为第i条土条块的黏聚力(kN/m2)；bi为第i条土条块的宽度(m)；Wi为第i条土条块的重力(kN)；φi为第i条土条块的内摩擦角(°)；θi为第i条土条块与水平面的夹角(°)；rw为水体重度(kN/m3)；hti为浸润线以下土条高度(m)；FB为土坡安全系数。

3 求解方法

3.1 离散格式

由于溃口区水流非恒定性强，水位梯度大，需要引入较大的数值黏性才能保持格式的数值稳定，而数值黏性的增大会导致溃坝波峰值被抹平，降低模拟精度。为了提高计算的稳定性模拟精度，基于MacCormack有限体积法，引入间断波捕捉格式，在其校正步实施通量修正。目前TVD,ENO,AUSM等格式得到了广泛的应用，引入由Yee修正后的Harten TVD格式[12-13]，即：

ΔtS(t+Δt) ；

(8)

(9)

式中：目标p,c分别代表预测步和校正步;F代表界面法向通量；A代表界面法向面积；Δt为时间步长；ΔV为控制体体积;i为网格边数。

式(9)中

(11)

ψ(ak+1/2+γk+1/2)αk+1/2。

(12)

式中:ak+1/2为右特征值；αk+1/2为特征列向量元素，利用函数ψ对αk+1/2进行熵修正,ψ的表达式为

式中δ为一极小正数。式(12)中γk+1/2的表达式为

(14)

式(12)和式(14)中限制函数gk采用minmod限制器，即

(15)

3.2 网格划分与干湿边界

采用非规则局部不连续四边形结构网格划分计算区域[7]，使模型具备对复杂地形的适应能力，同时可减少无效计算网格，提高模型的计算效率。针对溃坝洪水波在传播过程中，引起计算网格剧烈干湿变化，采用了与数值解法相匹配的基于单元属性与基于界面属性的动边界处理技术[7]。

后来，宴姝去英国名校利兹大学念书。英国的博物馆资源丰富，上学的3年间，宴姝最大的印象就是跟着老师跑了各种博物馆、画廊。

4 工程算例

4.1 工程概况

洈水水库位于湖北省西南部松滋洈水镇，与湖南省澧县接壤，坝址处于洈水流域中游洈水镇，距松滋城区约40 km。洈水水库主坝为土坝，最大坝顶长1 640.00 m，最大坝高42.95 m，挡水压力大，一旦发生溃决，溃坝洪水将对下游形成严重威胁。

4.2 成果分析

4.2.1 溃口冲蚀与发展过程

图1为沿溃口中心线冲蚀纵剖面过程图，可以看出，溃口形成1～2 min，水流主要冲刷坝轴线下游约55 m、高程71 m以上的坝面； 3 min以后，坝轴线下游约35 m处附近的坝面冲刷加剧，约5 min冲蚀至65 m高程，随后快速向上、下游发展。

图1 溃口中心线冲蚀过程纵剖面Fig.1 Longitudinal process lines of erosion along breach center line

图2为溃口冲蚀与发展平面过程。由冲蚀发展过程可知，0～4 min坝面溃口逐步展宽，同一时刻从坝顶至坝脚溃口宽度沿程减小；6 min左右，坝顶至1/2坝高范围的溃口发展较快，其宽度相对其上、下游较大；8 min左右，上、下游坝面的溃口宽度迅速增大，至10 min左右形成方向相对的“喇叭”型溃口，且沿溃口中心线，坝体均已溃至65 m的坝脚高程；15～35 min溃口向左右两岸展宽，至39 min溃口基本稳定。

图2 溃口冲蚀与发展过程平面图Fig.2 Plane diagrams of breach erosion and developing process

4.2.2 水面过程与稳定性

溃口区0～20 min瞬时水面见图3(图中灰色为水面)，可以看出，随着溃口的逐渐展宽，溃口下泄流量迅速增大，4～8 min，溃坝洪水快速淹没坝下低洼河槽，随即淹没河漫滩；10 min左右，洪水上溯至两溢洪道出口处。

图3 溃口区瞬时水面图Fig.3 Instantaneous water surface images of breach zone

观察溃坝洪水的演进可以看出，模型较好地模拟了蓄洪区在洪水到达时陆域迅速被淹没的过程，以及洪峰过后地面逐渐出露的过程，表明模型的动边界技术是合理可靠的。

4.2.3 断面流量过程

95.77 m库水位条件下主坝溃决后溃口及坝下游各断面流量过程线见图4。从图4可以看出，主坝溃决后溃口处流量迅速上升并达到39 074 m3/s的峰值流量，在峰值附近短时振荡，随后逐渐减小。坝下不同里程5个断面监测到的洪峰流量沿程依次减小，至坝下22 km处的金鸡山断面，洪峰流量减小至6 166 m3/s。各断面流量过程线较为光滑，除峰值附近短时低频振荡外，未出现高频振荡，表明数学模型的稳定性较强。

图4 坝下游各断面流量过程Fig.4 Process lines of discharge in different sections in the downstream of dam

4.2.4 守恒性分析

95.77 m库水位条件下，库区的初始水量8.347 7×107m3，模拟计算3 h后统计计算区域内的水量为8.339 8×107m3，水量误差约-1‰。

5 结 论

本文建立的耦合数值模型合理地模拟了坝体的溃决过程与坝下洪水演进过程，溃口的冲蚀与发展符合土石坝体溃决的一般规律；模型在溃口附近的急缓流过渡区表现出对间断波良好的捕捉能力，在较大水位梯度区表现了良好的稳定性；采用的动边界处理技术对干湿剧烈变化区域具有较强的模拟能力；模型的水量误差约1‰，守恒性良好。

参考文献：

[1] 刘 宁.21世纪中国水坝安全管理、退役与建设的若干问题[J].中国水利,2004，(23):27-30.

[2] HANSON G J, ROBINSON K M, COOK K R. Prediction of Headcut Migration Using a Deterministic Approach[J].Transactions of ASAE, 2001, 44(3):525-531.

[3] TEMPLE D M, MOORE J S. Headcut Advance Prediction for Earth Spillways[J]. Transactions of the ASAE, 1997, 40(3):557-562.

[4] 朱勇辉,廖鸿志,吴中如. 国外土坝溃坝模拟综述[J].长江科学院院报,2003，20(2):26-29.

[5] 陈生水, 方绪顺, 钟启明, 等. 土石坝漫顶溃坝过程离心模型试验与数值模拟[J]. 岩土工程学报, 2014,(5):922-932

[6] FRACCAROLLO L, TORO E F. Experimental and Numerical Assessment of the Shallow Water Model for Two-dimensional Dam-break Type Problems[J]. Journal of Hydraulic Research, 1995, 33(6):843-864.

[7] 崔 丹,槐文信,姜治兵. 溃坝洪水演进的数值模拟[J]. 华中科技大学学报(自然科学版), 2012，(7):34-37.

[8] 夏军强,王光谦,LIN Bin-liang，等.复杂边界及实际地形上溃坝洪水流动过程模拟[J].水科学进展,2010,21(3):289-298.

[9] 杨 志，冯民权. 溃口近区二维数值模拟与溃坝洪水演进耦合[J]. 水利水运工程学报, 2015,(1):8-19.

[10] 岳志远,曹志先,闫 军.滑坡体溃决洪水数学模型研究[J].水动力学研究与进展,2008,23(5):492-500.

[11] 李相南,陈祖煜. 两种溃坝模型在唐家山堰塞湖溃决模拟中的对比[J].水利水电科技进展, 2017，(2):20-26.

[12] YEE H C. Explicit and Implicit Multidimensional Compact High-resolution Shock-capturing Methods: Formulation[J]. Journal of Computational Physics, 1997,131(1):216-232.

[13] YEE H C, SJÖGREE B. Adaptive Filtering and Limiting in Compact High Order Methods for Multiscale Gas Dynamics and MHD Systems[J]. Computers & Fluids, 2008, 37(5):593-619.