MIKE软件技术在水库溃坝洪水数值模拟中的应用研究

管 静

(河北省承德水文勘测研究中心,河北 承德 067000)

0 引 言

水库大坝发生溃坝的概率虽然较低,但由于我国气候和地势因素,超负荷的洪水与强震同样也会导致溃坝的形成[1-2]。溃坝洪水的流动状况十分复杂,很多因素都会影响到溃坝的洪水。目前,对其模拟主要包含物理模拟、数学模拟以及二者结合的模拟方法[3-4]。其中,数学模拟因为其灵活强、资金消耗较少等优势在国内外得到广泛的应用[5]。高普阳等[6]针对河流下游存在相关障碍物的溃坝流动实际问题,在两相流动模型与有限元算法的基础上,提出了牛顿流体溃坝流动相关的数值模拟方法。Maghsoodi R[7]针对溃坝障碍物对河流影响的相关问题,在流体体积和标准k-1模拟水面的相关方法基础上,对其进行了数值模拟。吕松峰等[8]针对三维地形尾矿库溃坝的相关问题,在颗粒流离散元法与短矩阵离散元软件的基础上,构建了实际云南某尾矿库完整的三维地形,并以此对其下泄演进进行了数值模拟。

在此背景下,本文引入土坝失事洪水演算模型(BREACH)与二维数值工程软件(MIKE21),并将二者组合构建成耦合模型,同时将其运用在河北邢台水库溃坝的洪水数值模拟应用中。其目的是为邢台水库在形成溃坝之前制定合理的防洪避险方案,并通过预测溃坝对下游的影响,为下游制定进一步的防洪计划和转移方案提供理论依据。此外,通过将耦合模型应用在溃坝洪水数值模拟过程中,制定7种不同的方案,既能全方位分析不同情况的影响,又具备一定的创新性。

1 基于MIKE技术的水库溃坝洪水数值模拟分析

1.1 水库溃口流量计算模型研究

在实际的水库溃坝中,不同的大坝类型会导致不同的崩塌形式。根据坝体的溃决历史,可以将其分为两类,即暂时性溃决和渐进溃决[9-10]。通常来讲,混凝土坝体可以按瞬间溃决形式进行计算,而土石坝往往采用逐步溃决模型进行模拟分析。因此,在确定溃坝方式之后,需要选取合适的数学模型来动态分析溃坝的大小和流速,并以此来模拟溃决下游的洪水演变。目前,用于水库溃坝洪水数值模拟的数学模型有很多,研究依据模型适应能力水平,选择BREACH模型。

BREACH模型可以对溃口特性的演变过程以及溃决后的泄洪过程进行预测,同时也可以用于模拟由于漫顶或管涌而造成的溃坝,坝体可以是均匀的,也可以是由两种材料组成的心壁和外壁[11-12]。BREACH模型的基本原理是将水力学、泥沙运动、土力学原理等因素进行综合考虑,首先假定溃口是长方形的,然后在溃口处的崩塌和变形中形成一个梯形,当溃口超过一定深度时,洪水就会停止冲刷溃口处的土体。溃口不再发展,而是趋向一个稳定的形态,由此可以预测堤坝的溃口大小和溃口的流速。

BREACH模型可以细分为漫顶溃决和管涌溃决。对于漫顶溃决来讲,当坡面没有植被覆盖时,由于水流的冲刷,在下游坡面上逐渐形成一个矩形的沟道。此时,河流水流流量计算公式如下:

Pa=3A0(h-hb)

(1)

式中:Pa为溃口河道的实际流量;A0为初始矩形河道的瞬间宽度;h为坝前的实际水位高度;hb为河道底部的实际高度。

因此,水坝溃口的实际发展过程见图1。

图1 水坝溃口的发展示意图

由图1可知,在植被覆盖的情况下,如果坡面上的水流速度超过允许速度,下游的坡面就会受到冲刷。随着植被的减少,河流的冲刷会加剧。随着洪水的冲刷,溃口的深度会越来越大,当达到一定的稳定程度后,溃口就会失去稳定,然后慢慢向两边蔓延,最后形成一个梯形。

对于管涌溃口来说,在渗流过程中,水流作用会使土壤中的细颗粒被冲刷出粗大孔洞,在土壤中形成一个贯穿的渗流通道。随着时间的推移,逐渐发展为暗渠式的溃口,当侵蚀加剧时,就会发生溃坝。此时,溃口流量的计算公式如下:

Wa=M[2λ+(I-Iρ)/1+fL/E]0.5

(2)

式中:Wa为渗透通道的实际流量;M为溃口的横断面积;λ为常数;I-Iρ为溃口的水头;f为摩擦因子;L为渗透管道的实际长度;E为渗透管道的实际直径。

确定溃决的相关模式后,可以选取BREACH模型运算所需参数。研究选取的模型运算参数见图2。

图2 BREACH模型运算所需参数示意图

由图2可知,BREACH模型运算时,所需参数包含坝顶与坝底的高度、坝高的最大值、溢洪道的顶部高度、坝顶的实际长度与宽度、上游与下游坝坡比、坝体材料的粒径和黏聚力、校核洪峰流量和干流的平均坡降、坝体材料的湿密度和摩擦角。

1.2 MIKE21模型原理及构建分析

在河流流量计算的基础上,研究选择MIKE21软件作为后续洪水演进的数值模拟软件。MIKE21是一种分析河流、河口等的二维化仿真模拟工具[13-14],其特点见图3。

图3 MIKE21软件技术的特点示意图

由图3可知,MIKE21软件具备无障碍性,即具备完整的人机交互功能,容易使用;高效性,即拥有快速可靠的模拟引擎;便捷性,即具备相应的支持软件用于数据处理和分析;广泛性,即具备满足各种河口模拟需求的相应模块;兼容性,即用于在第三方程序中处理模拟数据的地理信息系统集成和工具;灵活性,即多种计算网格,模块和许可选项保证用户可以根据实际需求选择模型;公认性,即具备长达25年的相关记录以及在全球广泛应用所证明的技术。

MIKE21的水动力数学模型以二维非稳定流动为基本控制方程,其基本假定是以水为不可压缩水体,并且压力在水深上是按静压力分布的[15]。其中,二维非稳定流动基础控制中的空间离散常用的是有限容积方法,即控制容积方法。在实际的二维模型中,由于模拟区的边界往往是不规则的,因此其把二维计算区分成一些无组织的或不规则的小单元,再对各单元进行水量和动量的均衡,从而得出各计算单元的水位和流量。

二维模型一般采用结构化网格和非结构化网格。其中,非结构化网格可以使网格合理分布,适合于复杂地形的网格划分。研究针对河北省某县的边界条件,将其分为不规则的非结构化网格。网格的大小随着地形的变化而变化,对高程变化小、边界规则的地区划分范围更广;对高程变化大、边界复杂的地区划分的面积会相应减少,也会进行局部加密,以增加计算效率和准确度。

因此,研究利用MIKE21软件对河北省邢台水库的下游城区进行二维模型计算。计算面积21.88km2,而网格剖分的相关部分选择不规则的三角形网络,通过剖分得到25 000多个网孔,平均网孔面积870×10-6km2,最大网孔为5 260×10-6km2,最小为12×10-5km2。MIKE21模型一般将网孔尺寸控制在0.05km2以内,因此研究划分的网孔尺寸满足邢台水库溃坝洪水模拟的精度要求。

在溃坝洪水模拟中,边界条件对计算结果的精确度有较大影响。在二维模型计算中,存在着开放边界和封闭边界两种情况,前者通常为陆边界,后者通常为水边界。在实际的溃坝洪水数值模拟中,依据网格划分水库溃坝洪水演进的数值模型流程见图4。

图4 网格划分方法下溃坝洪水演进模型流程示意图

由图4可知,构建溃坝洪水演进模型流程首先是对水库下游相关区域地形进行必要的处理,研究选择的处理软件为地理信息系统(Geographic Information System,GIS)。其次是利用MIKE21软件内部的相关建模工具,对网格进行相应的划分;之后将GIS提出的相关高程点导入网格文件,进行数字高程插值,并利用MIKE21生成模拟区域的实际地形图;然后对BREACH和MIKE21模型的相关参数进行设置;最后将BREACH的数值计算结果视为边界条件,导入MIKE21模型中进行耦合,以此完成整个模型的实际构建。

2 耦合模型在溃坝洪水数值模拟中的应用

为了验证BREACH与MIKE21的耦合模型在溃坝洪水数值模拟中的应用效果,研究选取河北省邢台水库作为研究对象,将模型应用在该水库某年洪水事故下的溃坝洪水数值模拟中,并将该水库下游水库设定为Q。试验前,依据溃决原因,将溃决模式分为因洪水漫顶溃决的邢台水库和Q水库连溃和Q水库单溃(方案1和方案7)、因地震溃决(方案2和方案3)和管涌溃决(方案4至方案6)的邢台水库单溃。其中,因为溃决水位不同又可以将邢台水库单溃分为5种方案,总计7种方案。因此,利用耦合模型模拟的不同方案下邢台Q上下水库溃口处流量过程曲线见图5。

由图5可知,不同方案下的溃口洪峰流量均大于10 000年一遇的入库洪峰高达2.31×104m3/s的流量。其中,方案1最大的实际下泄流量为18.39×104m3/s;方案2为15.31×104m3/s;方案3为13.30×104m3/s;方案4为8.41×104m3/s;方案5为2.37×104m3/s;方案6为2.99×104m3/s;方案7为2.71×104m3/s。

图5 不同方案下的邢台和Q水库溃口处流量过程曲线示意图

综合来看,不同方案下的溃坝洪水是以单一的波浪形式向下游扩散的,在溃坝后,溃坝的下泄流量很快达到一个高峰,随后逐渐下降。根据各种方案的对比发现,溃坝溃决的水位越高,溃坝的洪峰流量越大;洪水决堤时间越短,洪水的冲刷流量越大。随着决口水位的升高,泄洪时间也在逐渐延长。在邢台水库和Q水库连崩时,Q水库的洪峰流量依赖于邢台水库的溃决流量。邢台水库一旦出现单溃,丰宁上下水库就会同时溃堤。

在此基础上,针对各方案下邢台水库位置至Q水库以及Q水库下游沿线溃坝洪水的实际代表断面进行分析。由于篇幅所限,研究选取沿邢台七里河的5个重要断面地点,分别为两水库之间的两个地点以及Q水库下游沿河3个地点,分别用A、B、C、D、E表示。实际的洪峰流量变化趋势见图6。

由图6可知,方案1决口洪峰向Q水库的洪峰流量为1.46×104m3/s,比邢台水库的泄洪速度下降约20%;洪峰向地点C汇入点的洪水最大流量为1.01×104m3/s,比邢台水库坝址下降约44%。地点E的洪峰流量达到7×104m3/s,比邢台水库的洪水下降约57%。而地震溃决模式的方案2和方案3从邢台水库向Q水库传播时,流量分别降低约18%和17.9%;流至地点E时分别降低约60%和59%。管涌溃决模式的方案4和方案5从邢台水库向Q水库传播时,流量分别为7×104和2×104m3/s,分别降低约16%和4%;最后的方案6和方案7表现出同样的效果。

图6 各方案下溃坝洪水沿线洪峰流量的变化趋势

综合来看,在同一溃决水位条件下,溃坝洪峰流量随长度的增大而逐渐减小。而同一溃决时间段内,溃决时间随溃决水位下降而减小;在同一溃决水位上,洪水消退时间愈长,洪水流向下游的速度愈慢。另外,溃坝洪水向下游演进的实际过程中,不同方案下的沿线主要断面最高水位都有所不同。研究利用耦合模型对其进行数值模拟,见图7。

图7 不同方案下河流沿线最大洪水水位的变化趋势示意图

由图7可知,方案1至方案7在Q水库时的最高水位分别为306、300、297、286、272、270以及270m。到地点C时的最高水位分别为279、274、271、261、251、249以及247m。综合来看,决口水位和决口时间是决定河道最高水位的重要因素,随着决口水位的升高,河道的最高水位也随之升高。在同一决口水位下,决口时间愈短,最高水位愈高。

研究利用耦合模型,在5个重要断面中选取其中两个断面,并对其水位与时间的过程进行分析,结果见图8。

图8 不同方案下两个重要断面的水位特征

由图8可知,方案1的溃坝洪水水位始终最高,其中A地最高为306m,B地最高为290m。而方案7的溃坝洪水水位始终最低,其中A地最高为265m,B地最高为263m。综合来看,各断面洪水的水位随着时间不断向前推移而瞬间上涨,达到最高水位之后,又缓慢降低,直到稳定水位。

3 结 论

为了在水库遭遇溃坝之前制定合理的防洪避险预案,本文将BREACH模型与MIKE21融合成耦合模型,并将其运用在河北邢台水库溃坝的洪水数值模拟试验中。结果表明,不同方案下的溃口洪峰流量均大于高达2.31×104m3/s的流量;方案1决口洪峰向Q水库的洪峰流量为1.46×104m3/s,比邢台水库的泄洪速度下降约20%;7个方案在Q水库时的最高水位为方案1的306m;方案1的溃坝洪水水位始终最高。综合来看,本文提出的耦合模型,在模拟河北省邢台溃坝洪水的数值演变上具备较高的有效性。