渝河(隆德段)洪水数值模拟研究

杨东旭, 李金燕

(1.宁夏回族自治区水旱灾害防御中心,宁夏 银川 750001;2.宁夏大学 土木与水利工程学院,宁夏 银川 750021)

洪水分析模拟是指通过洪水计算软件计算特定频率洪水的淹没风险特征信息，确定特定洪水重现期或特定方案下的洪水位、可能淹没范围和淹没水深等；洪灾损失评估是指通过灾情损失评估模型，结合区域洪水分析模拟结果，据此计算特定方案的洪灾损失.

一维(1D)水动力学模型主要是对水流运动过程进行模拟，它可以将复杂河网的水位变化及流量变化过程快速、准确地模拟出来，但模拟的流动路线难以避免“一维线性”特征的出现.而对于水平漫流过程[1]的模拟，主要通过二维(2D)水动力学模型来实现，它是将时间进行有序的间隔设定，对淹没深度进行模拟计算.对于二维水动力学模型来说，它具有计算时间长，对地形数据的分辨率要求较高[2]等缺点.由此可见，一维、二维水动力学模型描述水流问题的侧重点是不同的.一般来说，一维水动力学模型是用来模拟河道洪水演进的流动过程，二维水动力学模型则是模拟解决防洪保护区的洪水流动问题.本文将两种模型进行侧向耦合求解，不仅可以将溃堤、漫堤洪水的淹没特征通过二维平面凸显出来，还可以提高水动力学模型的计算效率和模拟精度[3]，有效减少网格数量，增加模型可靠性[4].对于分洪对河道水位产生的影响，以及河道水位由于某些因素下降对蓄滞洪区或者城市分洪量造成的影响等问题，该耦合模型系统内部均可进行自动计算分析.

本文应用洪水仿真、GIS、数据库等技术，获得防洪保护区内洪水演进的可能途径和洪水所造成的淹没情况，便于防汛部门制定防洪调度决策，及时部署指挥防洪抢险救灾任务，补充和完善洪涝灾情评估，为灾后重建和土地开发利用提供管理依据.

1 数值模型

MIKE11模块可以很好地模拟分析一维河道洪水的演进情况，该模型计算稳定且计算精度高，具有较强的可靠性，适用于河道水文要素的计算，例如河道任意时刻、任意断面的水位观测和流量计算.该模型还可用来模拟水库调度、蓄滞洪区运用，分析溃口等构筑物对河道洪水演进的影响.MIKE21能够较为精准地计算该保护区内的二维洪水演进情况，可以具体模拟出水流水平方向和垂直方向的变化.以下为模型的计算原理.

1.1 一维水动力学模型

一维水动力学模型采用MIKE11模块软件模拟.模拟时需构建一维非恒定流Saint-Venant方程组，即水量方程和动量守恒方程[5]，通过Abbott六点隐式格式离散方程组，按序交替计算每个网格节点处的水位和流量.这种离散格式使得大步长计算时能在非常大的常数下保持稳定，同时节省计算时间.

(1)

(2)

式中：B为平均宽度；Z为水位；t为时间；Q为流量；s为距离；q为旁侧入流单宽流量；A为过水断面面积；g为重力加速度；C为谢才系数；R为水力半径；i为坡度.

1.2 二维水动力模型

二维水动力学模型采用MIKE21模块软件模拟，该模型的控制方程采用改进形式的二维浅水方程[6]

(3)

动量方程为

(4)

(5)

式中：H为水深；Z为水位；q为旁侧入流单宽流量；t为时间；M与N分别为x和y方向的垂向平均单宽流量；u和v分别为垂向平均流速在x与y方向的分量；n为曼宁糙度系数.

1.3 一维二维耦合模型(MIKE FLOOD)

2 渝河(隆德段)洪水模拟

图1 渝河计算区建模范围示意图

2.1 河道断面设置与网格剖分

2.1.1河道断面设置 河道断面是一维水动力学模型的重要基础数据[10]，将已收集的河道断面数据按照MIKE11的断面输入要求整编处理，然后通过内插加密，提高计算精确度，使计算结果最大限度贴近实际值.对于河道形态变化显著的地方以及涉及工程(桥、闸、坝、堰等)的位置，断面需进一步加密处理.经一维水动力学模型计算，渝河河段长度为23.65 km，共设置了265个河道断面，河道布置见图2.

图2 渝河河道断面布置

2.1.2网格剖分 将二维非恒定流方程作为网格剖分时的基本控制方程，对模拟范围进行无结构不规则网格划分，要求最大的不规则网格面积不能大于0.5 km2.网格剖分时需要根据地形、历史洪水淹没范围、防洪工程分布情况等因素划定风险区域，以便于后期分析计算，因此,地形变化较大的区域以及重要地区的网格要适当加密.本文利用网格剖分将渝河区域的二维平面划分成17万个小网格，其中面积最小为20 m2，最大为800 m2，总的计算面积达到43.21 km2，网格剖分示意图见图3.

图3 二维区域网格剖分示意图

2.2 参数确定

2.2.1河道一维水动力模型计算参数 本文将河道粗糙率作为影响渝河河道一维水动力学模型分析精度的计算参数[11]，参考相关设计报告[12—13]，设定渝河河道综合糙率为0.03.将16座过水路面作为河道建筑物处理，概化为堰，鉴于路面具有阻水作用，根据设计参数设定过水尺寸时不可忽略路面阻水因素.综合考虑模型运算效率及模型本身的稳定性等多种因素[14]，将该模型计算迭代步长设定为1 s.根据渝河历史洪灾数据记载，结合渝河的堤防现状，在渝河(三里店水库至小河口段)设置2处溃口，将溃口处河道水位与堤防设计水位持平时设为堤防溃决时刻临界点，即此时堤防瞬间溃决，溃口宽度为60 m.和平村溃口底高程为1 901.7 m，神林村溃口底高程为1 818.3 m.

2.2.2计算区二维水动力模型计算参数 参考相关文献，将防洪保护区综合糙率值设为0.06[15—16]，综合考虑模型运算效率及模型本身的稳定性等多种因素[17]，将该模型计算迭代步长设定为10 s，最小迭代步长设为0.01 s.

2.3 边界条件

根据支流汇入情况，将渝河分为3个河段分别进行洪水风险分析，本文选择的洪水分析量级为20 a一遇、50 a一遇和100 a一遇，在渝河(三里店水库至小河口段)设置2处溃口.根据实际调研结果，渝河洪水溃堤之前，保护区内初始水深为0 m.因此，根据干水深和湿水深理论，可将计算区域的网格设定为干单元.

2.3.2渝河计算区二维水动力学模型边界条件 渝河计算区二维水动力模型外边界具有闭边界和开边界两种.闭边界即为与外界无水流交换的边界，开边界为一维模型和二维模型水流动态交互的边界，河道溃口侧向建筑物耦合的二维区域边界可看作模型开边界.以地面为原点，将高于地面0.5 m的阻水建筑物或过水建筑物作为水动力学模型的内部边界，例如堤防、线状地物、涵洞、桥梁等.一般来说，洪水演进是指可以正常通过涵洞、桥梁等过水建筑物或者具有沿程缺口的线状建筑物，而对于堤防等阻水建筑物来说，洪水的演进是以漫溢形式通过的.河道一维模型将为计算区平面二维模型实时提供漫溢或溃口分流流量，作为固定时间步长内二维模型的入流条件.

3 模拟结果与合理性分析

通过对渝河(隆德段)10 a一遇、20 a一遇、50 a一遇以及100 a一遇的洪水进行模拟，最终建立渝河(隆德段)的洪水预案库.本文将以100 a一遇洪水为代表进行模拟结果分析.

3.1 模拟结果分析

3.1.2溃口分流分析 和平村、神林村溃口(上游临近断面)河道流量变化过程见图4，溃口处河道水位变化过程见图5，溃口分流流量过程见图6.

图4 渝河100 a一遇洪水在和平村、神林村溃口(上游临近断面)河道流量变化过程

图5 渝河100 a一遇洪水在和平村、神林村溃口处河道水位变化过程

本文将河道水位与设计水位相等时设定为河道溃堤临界时刻[15].由图4和图5分析可知，和平村溃口处起溃时刻河道流量为120.3 m3/s，河道洪水水位为1 901.9 m，分流历时为5.75 h；神林村溃口处起溃时刻河道流量为141.7 m3/s，此时河道洪水水位为1 818.7 m，分流历时为5.5 h.由图6分析可知，和平村溃口起溃时刻分流流量为16.7 m3/s，此时分流洪峰流量达到63.1 m3/s，而神林村溃口由于临近上游断面，起溃时刻分流流量是和平村的2.1倍，约为35 m3/s，同时其分流洪峰流量也达到了118.9 m3/s.

图6 渝河100 a一遇洪水在和平村、神林村溃口分流流量过程

3.1.3漫溢分流分析 通过对渝河河道百年一遇洪水进行一维模拟发现，共有18处发生洪水漫溢现象，其中16处为过水路面.漫溢上段漫溢位置6处，下段漫溢位置12处.

3.1.4防洪保护区洪水淹没分析 洪水进入二维平面计算区域的形式有溃流和漫溢两种形式，选取1, 2, 4, 6.5 h共4个时段分别对两个村溃口洪水淹没水深的分布情况(图7)进行比较分析.

由图7可知，洪水演进1 h，和平村溃口已起溃，洪水由溃口向保护区演进.由于过水路面使水位壅高，渝河部分过水路面处出现洪水漫溢现象，计算区内淹没面积1.13 km2.洪水演进2 h，神林村溃口已起溃，洪水由溃口向保护区演进.渝河上段漫溢洪水水深大部分在0.5 m以下，下段河道洪水继续向外漫溢，多处水深在0.5～1 m之间，由于溃堤洪水叠加漫溢洪水，神林村附近水深达到1 m以上，洪水向南北两侧演进过程中受到国道312的阻挡，路前最大积水深1.2 m，计算区淹没面积4.9 km2.洪水演进4 h，漫溢洪水继续向下游演进，国道312阻水效果明显，淹没面积6.16 km2.由于河道来洪减少，河道内水位降低，二维保护区下游洪水通过堤防较低处流向河道.最终，淹没面积为6.22 km2，积水量300.67万m3.受大洪水冲刷影响，沿渝河河道桥梁、护坡护岸及堤防工程、过水路面、引退水工程等工程设施将遭受严重损害.

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

3.2 合理性分析

3.2.1水量平衡分析 本文的水量平衡关系可通过渝河防洪保护区进洪量和区内蓄水量是否相等得以验证.针对渝河防洪保护区设定的6个洪水方案分别进行进洪量和蓄水量的对比分析(见表1)，结果表明渝河不同洪水方案下进洪量等于蓄水量，误差为0 m3，满足水量平衡要求.

表1 渝河防洪保护区进洪量与蓄水量对比表

3.2.2流场分布 1)整体流场分布.通过对洪水淹没的二维平面区域进行模拟，最终将模拟计算出的整体流场分布情况与DEM整体高程进行对比分析(图8和图9).对比图8和9可知,淹没区域流场分布均匀一致，而低洼地带是较大流速主要集中的区域，根据洪水的自然流动原则，即从高到低流动，洪水呈现出的流动态势较为准确.

图8 和平村DEM地形示意图

图9 渝河遭遇100 a一遇洪水在和平村溃堤洪水最大流速示意

2)局部流场分布.通过对局部流场进行模拟分析，将模拟计算出的流场分布情况与线状地物进行比较，对比分析结果见图10和图11(以和平村溃堤为例).结果表明，局部区域的流场分布均匀一致，且线状地物的阻水效果比较明显，洪水呈现出的流动态势较为准确.

图10 二维平面区域流场均匀分布

图11 线状地物阻水效果流场分布

3.2.3同一方案风险信息比较 如图12～图14所示，以和平村溃堤为例，将同一方案的DEM、洪水淹没水深及洪水流速进行比较分析.结果表明，DEM地势较高地区淹没水深较小，相反，地势较低洼地区淹没水深较大，而对于处在地势由高到低的区域的洪水，即DEM过渡地区的洪水，洪水流速偏大.这初步说明该方案计算具有合理性.

图12 和平村溃堤DEM分布

图13 和平村溃堤淹没水深分布

图14 和平村溃堤洪水流速分布

3.2.4不同方案风险信息比较 选定某一固定位置在同一时刻对不同方案的洪水流速、淹没水深等洪水风险信息进行比较分析.以和平村溃堤为例，对渝河50 a和100 a一遇洪水进行风险信息比较.结果表明，同一位置同一时刻100 a一遇的洪水流速略大于50 a一遇洪水流速(图15)；对于淹没水深，100 a一遇洪水略大于50 a一遇洪水，比较结果见图16.通过对不同方案风险信息进行验证分析，可以发现本文模型计算较为合理.

图15 同一位置同一时刻洪水流速比较

图16 同一位置同一时刻淹没水深比较

4 结语

综上，对渝河(隆德段)洪水的耦合模拟，无论是模拟计算过程还是对结果的合理性分析，都对渝河区域洪水治理提供了一定的理论支撑，为水利防汛工作的开展提供依据，有利于防汛应急预案的快速制定及部署.对于水利工程建设与管理相关部门来说，本文的模拟结果可以为渝河防洪治涝工程体系的规划提供基本依据，便于开展有关洪水风险的宣传及相关教育培训等.