基于一二维耦合技术的溃坝洪水模拟研究

程 涛，何万成，赵海镜

(1.中国电建集团北京勘测设计研究院有限公司，北京 100024；2.国网新源控股有限公司，北京 100761)

随着社会经济的发展，世界各国对大坝的安全问题都给予了高度的重视。在研究溃坝洪水可能对河道流域下游带来危险的同时，也在采取积极主动的防范措施[1]。如：美国启动的国家大坝安全计划，欧共体则组织了10余个成员国家开展溃坝问题的合作研究，国际大坝委员会则在此基础上进行整合和协调，发布了ICOLD 111号溃坝洪水分析导则。这对于溃坝洪水的分析计算具有明确的规范和指导作用[2]。

巴丹托鲁水电站位于印度尼西亚苏门答腊岛北部的巴丹托鲁河上，苏门答腊岛位于亚欧板块和印度洋板块边缘地区，受地震影响频繁，根据印尼大坝委员会要求，巴丹托鲁水电站需开展溃坝洪水风险研究工作。巴丹托鲁水电站坝址位于山区与海滩交接处，地形条件特殊，不同下垫面条件对模型精度要求亦不同。目前，国内一二维水动力耦合模型在山区[3- 4]以及平原城市河道[5-7]应用较多，基于山地与海滩交接区域的模拟较少。本文根据巴丹托鲁水电站下游地形特点，建立了一二维耦合水动力学模型，针对拟定的溃坝工况进行数值模拟，为后期电站溃坝洪水风险分析和灾害损失评估提供技术支撑。

1 模型构建

溃坝数值模拟主要包括溃口模拟和洪水演进模拟两部分，根据工程特性和地形特点，本次溃口模型选用DamBreak模型，洪水演进模型选用MIKE FLOOD水动力模型，该模型均是在几十年来大量的工程应用经验的基础上持续发展起来的，曾在世界多个国家和地区得到过成功应用[8]。

1.1 模型区域

巴丹托鲁水电站采用混合式开发方式，坝型为重力拱坝，库容1 790万m3。电站厂房位于坝址下游约17 km处。坝址至出海口河道长度约80 km，其中27 km位于山区，53 km位于海滩地区，考虑到居民分布点主要位于下游海滩地区。本次建模范围确定为坝址至河道出海口两岸周边区域，覆盖面积约为1 500 km2。

综合考虑模型适应性和地形特点，山区河道和下游海滩地区分别采用一维和二维模型进行模拟，通过连接工具建立一二维动态耦合水动力学模型。该模型既有一维和二维模型的优点，又避免了采用单一模型时遇到的计算效率、网格精度和准确性方面的问题。模型区域具体分布见图1。

图1 模型范围分布示意

1.2 基本原理

(1)MIKE 11模型。MIKE 11模型控制方程采用一维非恒定流Saint-Venant 方程组，包括垂向积分的物质方程和动量方程，该模型具有稳定性好、计算精度高的特点。

(2)MIKE 21模型。MIKE 21模型控制方程采用平面二维非恒定流N-S方程组，包括水流连续性方程和动量方程。该模型可以计算三角形非结构网格，保证模型边界的光滑性，同时模型自身可以设置简单的水工结构物，通过结构物的组合模拟不同区域的复杂地形。

(3)MIKE FLOOD模型。MIKE FLOOD模型采用了一二维动态耦合技术，它可以将一维模型的末端断面与二维模型里对应一系列地形网格单元连接，在计算过程中，MIKE 11的计算结果作为边界条件加入MIKE 21网格进行二维计算，同时MIKE 21参与连接的网格平均水位也传递回MIKE 11进行一维计算，模拟效果更贴近于实际情况。

1.3 模型建立

一维模型模拟区域采用实测地形数据，根据河道走势和地形变化沿河道选取了50个断面作为模型的输入地形，断面范围覆盖了库区和坝址下游山区河道，长度约17 km。溃口模型DamBreak在MIKE11模型中以水工建筑物的形式添加至坝址断面，可结合断面地形进行溃坝流量计算。

二维模型模拟区域范围较大，根据雷达地形和SRTM数据，依照网格剖分原则和方法，选用非结构三角网格对模拟区域进行剖分和地形差值。模型整体网格数量约为41万个，网格边长范围10～30 m，同时根据地形数据调整网格尺寸，确保整体网格的平滑度，提高计算稳定性。

2 模型参数

2.1 模拟方案

根据国际上发生的混凝土拱坝溃坝案例数据统计，各机构推荐的溃坝特征参数较为接近，混凝土拱坝溃决形式基本为瞬间全溃[9]。

根据巴丹托鲁电站坝型特点和泄流能力分析，发生1 000年一遇洪水情况下，大坝泄流能力依然满足泄洪要求，重力型拱坝遭遇洪水溃坝风险较小，但巴丹托鲁水电站的所在地苏门答腊岛位于亚欧板块和印度洋板块边缘地区，受地震影响频繁，遭遇地震溃坝的风险相对较大。因此，本次溃坝洪水分析计算工况为地震溃坝工况，计算方案特征参数见表1。

表1 模拟方案特征参数取值范围

2.2 边界条件

一维模型上游边界采用多年平均流量106 m3/s，一维模型下游边界给定出口断面的实测水位流量关系。二维模型上下游地区设置为开边界，两侧水流未到达边界设置为陆地边界，二维模型上游流量边界通过耦合工具MIKE FLOOD给定，下游水位边界采用恒定海面高程0 m，模型计算时间步长为1 s。

河道断面糙率根据实测资料推算，取值为0.05～0.07。下游滩地范围较大，参考不同类型下垫面的糙率取值范围，综合实际土地利用情况确定滩地模型糙率，综合糙率取值为0.07。模型经天然河道水位、流量资料率定，模拟精度较高，可用于坝址下游洪水演进计算。

3 模拟成果

3.1 溃口流量过程

图2与图3分别展示了溃口处洪峰流量过程线以及水库水量变化过程线。由于巴丹托鲁水电站库容较小，整个水库水体下泄过程仅持续了约28 min，受溃坝引起的洪水波影响，溃口流量和库水位在溃坝发生后出现短期震荡现象，其中溃口最大洪峰流量为28 300 m3/s，随着库水位逐渐下降，溃口洪峰流量也逐渐减小。

图2 溃口流量过程

图3 水库水量变化过程

3.2 洪水演进过程

由于村庄和农田主要分布在下游滩地地区，因此图4统计了溃坝洪水进行过程中下游地区最大淹没水深，成果显示洪水在上游山区段水深变化较为显著，瞬时局部最大水深接近20 m，随着洪水向滩地行进，洪水大量漫出主河槽，对河道两岸的部分村庄造成了影响，但由于滩地面积较大，淹没水深远小于山区河道，滩地淹没水深在0～3 m范围内。

图4 最大淹没水深分布

为了更直观地展示溃坝洪水沿程的变化过程，从坝址到下游滩地选取了5个典型断面进行了水力要素统计，典型断面水力要素统计成果见表2。

表2 典型断面水力要素统计

由图4及表2可以看出，溃坝发生初期，坝址断面洪峰流量和行进流速最大，随着向下游演进，洪峰流量逐渐坦化，在山区河道行进约48 min即到达位于下游的巴丹托鲁村断面。巴丹托鲁村下游地区地形较为平坦，洪水经过巴丹托鲁村后大量漫滩，由于水库库容较小，下泄水量有限，溃坝水体释放后下游河槽流量急剧减小，同时流速降低，下游海滩地区河道受溃坝洪水影响较小。

3.3 敏感区域淹没分析

巴丹托鲁电站下游地区人烟稀少，根据溃坝洪水淹没水深计算结果，主要影响区域为位于山区和海滩交汇处的巴丹托鲁村，该村沿河分布有大量房屋和公路、桥梁，下游海滩地区沿岸村庄受洪水影响程度相对较轻。巴丹托鲁村断面最大洪峰流量5 880 m3/s，最高洪水位73.2 m，洪峰持续时间20 min，参考实测村子分布高程数据，确定出各设施遭遇溃坝洪水时的影响范围，受影响设施统计结果见表3。

表3 敏感区域淹没设施统计

4 结 语

本文以印度尼西亚巴丹托鲁水电站为例，对坝址下游地形条件复杂的水电站溃坝洪水数值模拟技术进行了探讨，总结出如下几点体会，供同类工程参考。

(1)大区域、复杂地形条件下溃坝洪水数值模拟需同时考虑计算精度与计算效率问题。本文建立的适用于不规则计算域和地形的溃坝洪水一二维耦合模型，引入了最新的动态耦合技术，避免了采用单一模型时遇到的计算效率、网格精度和准确性方面的问题。

(2)溃坝洪水在不同区域衰减速度差异较大，溃坝发生初期，山区河道洪峰流量和行进流速最大，山区河段洪峰坦化幅度相对较小，随着向下游演进，溃坝水体进入海滩区域后，洪水发生漫滩，下游河槽流量急剧减小，同时流速降低，下游海滩地区受溃坝洪水影响较小。

(3)溃坝时的水流运动过程常伴随有泥沙输移等其他物理过程，泥沙冲淤与水流运动会产生相互作用，后期可以水动力模型为基础，建立适用于溃坝模拟的水沙数学模型，实现溃坝洪水多过程耦合模拟。