溃堤山洪淹没风险评估水动力耦合模型及应用

田福昌， 张兴源， 苑希民

(天津大学 水利工程仿真与安全国家重点实验室， 天津 300072)

1 研究背景

多年来，我国一直高度重视山洪灾害防范与预警工作，山区溃堤洪水风险分析评估是其中一项重要研究内容。利用水动力数值模型模拟计算溃堤山洪演进过程，可为山洪预警和风险调度提供较为准确的风险信息。山区沟道溃堤洪水淹没过程主要包括沟道与保护区洪水演进以及两者之间洪水的动态交互，通过建立沟道与保护区洪水联合计算的一、二维水动力耦合模型，采用侧向耦合衔接方式实现沟道与保护区之间水量与动量的实时动态传递，能够较为准确地分析评估溃堤山洪淹没风险。

近年来，众多学者针对一、二维水动力耦合模型开展了大量研究工作。陈文龙等[1]建立了基于侧向联解的一维-二维耦合水动力模型，并对郁江中下游河段防洪保护区溃堤及漫堤洪水演进进行了模拟；田福昌[2]建立了考虑漫、溃堤位置处流态变化的河道-泛区二维动态耦合水动力数值模型，并对黄河宁夏段洪水漫滩运动以及青铜峡防洪保护区洪水演进过程进行了模拟；付成威等[3]利用建立的一、二维水动力耦合模型模拟了谷堆圩保护区溃堤洪水的演进过程；苑希民等[4]借鉴全二维气相色谱理论建立了漫溃堤联算的全二维水动力模型，并成功用于黄河宁蒙段河道及两岸灌区的漫溃堤洪水模拟计算；陈俊鸿等[5]建立了基于多种优化方法的一、二维耦合水动力模型，模拟了赣西联圩防洪保护区溃堤洪水演进过程；Dutta等[6]建立了有限差分法的一、二维水动力耦合模型，并对湄公河漫溢洪水运动过程进行了模拟；徐祖信等[7]建立了基于有限元法的一、二维水动力耦合模型，并应用于黄浦江平原感潮河网的流场模拟；Cai Xin等[8]与Xie Zuotao等[9]基于元胞自动机理论建立了洪水演进计算模型并模拟了荆江洞庭湖洪水流动过程；陈智洋等[10]建立了横阳支江一、二维水动力耦合数值模型，计算了横阳支江超标准洪水漫堤、溃堤等工况下的洪水演进情况；李长跃[11]、苑希民等[12]建立了考虑溃堤分流与暴雨等多源洪水条件下的一、二维水动力耦合模型，并将其应用于茨南淝左片防洪保护区多源洪水运动的耦合模拟。

本文建立模拟山洪沟溃堤洪水演进过程的一、二维水动力耦合模型，将计算区内道路、灌渠堤防等概化为宽顶堰并作线性处理，以模拟线状地物附近水流流态及其对洪水演进的影响。采用非结构化网格剖分计算区域，沿线状地物缩小网格尺度并加密剖分计算网格，实现线状地物等特殊边界附近网格地形及水流运动的动态耦联，在此基础上利用干湿边界理论对模型进行优化处理，建立具有真实地形的山区沟道溃堤洪水淹没计算模型，并将其应用于宁夏大武口沟溃堤洪水淹没风险的分析评估。

2 一、二维水动力耦合模型基本原理

2.1 一维水动力模型

一维水动力模型采用的基本方程如下：

(1)

(2)

式中：B为过流横断面宽度，m；qs为沟道旁侧流量，m3/s；QS为沟道总流量，m3/s；A为沟道过水断面面积，m2；g为重力加速度，m/s2；t为时间，s；Z为沟道水位，m；C为谢才系数；s为距离坐标；R为水力半径；i为沟道底坡降。

利用Abbott六点隐式格式[13]离散上述方程组，该离散格式在每一个网格节点按顺序交替计算水位和流量，因此能够在相当大的Courant数下保持计算稳定，可采用较大计算步长以节省计算时间。

2.2 二维水动力模型

二维水动力模型采用的基本方程如下[14]：

连续方程：

(3)

动量方程：

(4)

(5)

式中：H为水深，m；Z为水位，m；qc为连续方程中的源汇项；M与N分别为x和y方向的垂向平均单宽流量，m2/s；u和v分别为垂向平均流速在x与y方向的分量，m/s；g为重力加速度，m/s2；n为曼宁糙率系数。

2.3 一、二维水动力模型耦联方式

采用侧向连接方式实现沟道一维模型与保护区二维模型的耦合衔接。溃堤洪水流态与宽顶堰流较为接近，故溃口流量Qb采用宽顶堰公式[3,6]计算:

(6)

h1=max(Z1,Z2)-Zb，h2=max(min(Z1,Z2)-Zb,0)

式中:Z1、Z2分别为一维、二维模型在耦合界面处的水位，m；Zb为耦合界面的底高程，m；lb为溃口宽度，m。

3 研究实例

3.1 研究区域概况

大武口沟位于石嘴山市北侧，为贺兰山东麓最大的山洪沟。大武口沟流域面积576 km2，沟长50 km，平均比降11.50‰，较大支流有北叉沟、塔塔沟、八号泉沟、大灯沟等。沟道范围内地势北高南低、西高东低，海拔1 350～1 400 m。沟道右岸为山地陡岸，左岸地势较平缓，相对高差较小，两岸植被稀疏，岩石裸露，沟道砂砾等堆积严重。

大武口沟为季节性河流，洪水主要由暴雨产生，洪水特性主要表现为年际变化大、季节性特征明显(多出现在7-8月份)、来势凶猛、暴发频繁、洪峰流量随汇水面积的增大而缓慢增加、洪峰陡涨陡落、峰高量小、峰型尖瘦，洪水历时一般为6～12 h。本文计算区域为大武口沟城区段(自上游石大公路桥至下游平汝铁路桥)及其两岸保护区，总面积为395.19 km2。

图1 大武口沟模型计算区域示意图

3.2 溃堤山洪淹没水动力耦合模型建立

3.2.1 地形概化

(1)沟道断面设置。沟道断面是一维水动力模型的重要基础数据，大武口沟城区段计算范围为石大公路桥至平汝铁路桥(大武口拦洪库入库断面)，沟道长度8.02 km，根据沟道实际宽度及其蜿蜒曲折情况，对沟道断面进行内插加密处理，共设置160个计算断面，断面间距均值为500 m。

(2)网格剖分。非结构化网格具有很强的边界适应能力，能够对任意形状和联通区域进行网格剖分，便于控制网格密度，易于修改和调整，更容易获得高质量网格地形。因此本文采用非结构化三角形网格单元剖分研究区域，并在兴民村扬水渠、平汝铁路、S301、第二农场渠沿线以及区域边界进行网格加密处理，计算区域共划分网格31.20×104个，最小网格面积50 m2，最大网格面积3 500 m2，计算总面积395.19 km2。

3.2.2 参数设定

(1)干湿边界。模型计算过程中，设定干水深Hdry和湿水深Hwet的分界可以消除不稳定性并提高计算效率，利用干湿水深与网格淹没水深(H)进行比较，通过网格淹没水深与干湿水深的关系来判别网格的通量(动量通量和质量通量)计算。若网格H>Hwet，该网格为湿网格，质量通量和动量通量同时参与计算；若网格H

(2)糙率和计算步长。根据《宁夏石嘴山市大武口区大武口沟河流治理工程初步设计报告》、《洪水风险图编制技术细则(试行)》及《水力计算手册(第二版)》，确定大武口沟城区段沟道综合糙率为0.034，堤外平面计算区综合糙率值为0.06。综合考虑模型稳定及运算效率等多种因素[15]，设定大武口沟城区段沟道一维水动力模型计算迭代步长为1 s，计算区二维水动力模型最大计算迭代步长为10 s，最小计算迭代步长为0.01 s。为实现一维模型和二维模型固定时间步长内的动态耦合，耦合模型计算时间步长为1 s。

(3)溃口参数。大武口沟城区段堤防设计防洪标准为50年一遇，考虑计算沟道洪水溃堤分流淹没风险的最大化，按超标准洪水量级选取，故本文确定大武口沟城区段洪水分析标准为100年一遇。考虑大武口沟城区段保护区可能遭遇最大风险，综合河势地形、地质状况、工程状况和历史出险情况，确定溃口位置位于平汝铁路上游1 320 m处，并结合现场调查和防汛专家咨询意见等确定溃口宽度为60 m，堤防瞬间溃决，溃堤时水位达到设计洪水位。

3.2.3 边界条件 大武口沟城区段一维水动力模型上边界条件为石大公路桥断面遭遇100年一遇洪水，如图2所示；下边界为大武口拦洪库入库断面，控制条件为该断面水位-流量关系，如图3所示；采用侧向衔接概化方式计算溃口分流流量过程，以实现固定时间步长内沟道与保护区洪水的实时动态交互及耦合计算。

考虑到计算区域内道路、灌渠渠堤等线状地物对洪水演进过程干扰性较强，影响区域流场及淹没风险分布情况。本文考虑将研究区域内兴民村扬水渠、第二农场渠、平汝铁路和S301等线状地物作为模型内边界并概化为宽顶堰，即当洪水水位未达到线状地物顶部高程时，线状地物起阻水作用，区域不过水；当洪水水位超过线状地物顶部高程时，水流漫溢通过。同时考虑线状地物沿程缺口及桥梁等主要过流形式，将线状地物在桥梁或缺口处断开，断开宽度即为过水区域断面宽度。

3.3 溃堤洪水淹没计算结果与风险分析

3.3.1 溃口分流洪水过程分析 根据所建大武口沟城区段溃堤山洪一、二维水动力耦合模型，模拟计算大武口沟100年一遇洪水演进过程以及溃口分流洪水在保护区内的淹没运动过程，根据洪水计算结果提取平汝铁路桥上游处溃口(上游临近断面)沟道流量过程(如图4所示)和水位过程(如图5所示)，溃口侧向分流流量过程如图6所示。

平汝铁路桥上游堤防溃决时溃口处断面洪水水位等于设计水位，同时考虑大武口沟堤防现状防洪标准较高及抗洪抢险能力等多种因素，设定溃堤分流截止时溃口处沟道水位等于地面高程(溃口分流积水后地面高程)。由图4和图5分析可知，溃口处沟道流量为845 m3/s时堤防开始溃决，此时沟道洪水水位为1 108.03 m，溃堤分流截止时沟道洪水水位为1 107.49 m，溃口分流历时为7 h。由图6分析可知，起溃时刻溃口分流流量为35.74 m3/s，分流洪峰流量为197.10 m3/s，经统计分流总量为179.98×104m3。

3.3.2 溃堤洪水演进及风险分析 洪水由平汝铁路桥上游溃口处侧向分流进入二维平面计算区域，由所建大武口沟城区段溃堤山洪一、二维水动力耦合模型计算得到不同时刻洪水淹没水深分布情况，如图7所示。

图2 石大公路桥断面100年一遇设计洪水过程 图3 大武口拦洪库入库断面水位-流量关系图

图4 溃口上游临近断面沟道流量过程 图5 洪水平汝铁路桥上游处沟道水位变化过程 图6 平汝铁路桥上游溃口侧向分流流量过程

图7 洪水演进淹没水深分布图

由图7分析可知：受区域地形分布影响，溃堤洪水进入计算区后向东南方向演进，受到平汝铁路和兴民村扬水渠阻挡；洪水演进2 h，兴民村扬水渠与平汝铁路之间区域逐渐被淹没，洪水向第二农场渠方向演进，计算区内最大水深3.02 m，淹没面积1.12 km2；洪水演进4 h，洪锋北侧到达S301与平汝铁路之间中心区域，计算区内最大水深3.21 m，淹没面积2.49 km2；洪水演进8 h，受S301阻挡，中心区域大面积积水，计算区内最大水深2.71 m，淹没面积6.03 km2；洪水演进16 h时已基本趋于稳定状态，受兴民村扬水渠、平汝铁路和S301阻水作用，淹没区最大积水面积7.18 km2，最大积水深度为3.32 m(位于平汝铁路路基附近)。考虑洪水冲刷作用影响，大武口沟沿程桥梁、护坡护岸工程、水文测站、计算区内渠堤及道路路基等将遭受严重损害，应提前做好汛前检查与应急抢险防御工作。

4 结 论

(1) 建立一、二维水动力耦合模型模拟山区沟道溃堤洪水演进过程。为了准确反映计算区线状地物对水流运动的影响，采用非结构化网格剖分与局部网格加密处理技术对计算区域进行了地形概化，利用干湿边界理论优化模型，并将淹没区具有挡水作用的线状构筑物概化为宽顶堰，实现了淹没区特殊边界与非结构化网格的无缝耦联，以更好地拟合复杂地形，建立具有真实地形的山区沟道溃口分流洪水淹没计算模型。

(2) 将模型应用于宁夏大武口沟保护区溃堤山洪风险分析评估，准确高效地模拟了溃堤洪水演进过程和风险分布特征，明确了主要淹没区域与防洪重点关注区，计算结果可为大武口沟防汛指挥与应急抢险提供支持与参考。

(3) 所建大武口沟溃堤山洪淹没风险评估水动力耦合模型，可推广应用于沟道工程规划设计与防洪评价、洪灾评估与风险区划等领域。由于未有历史大场次山洪溃堤事件的详细记载资料，本文尚无法实现模型更为精准的率定与验证，后续可根据历史及现场资料搜集情况对其进行补充完善，以最大限度提高模型计算精准度，为防汛指挥与应急决策等提供更为准确的风险信息。