《1 前言》

1 前言

蓄滞洪区作为江河防洪体系重要组成部分，是城市防洪的最后屏障，在防洪危机时刻发挥关键作用。近年来随着流域经济社会的发展，滞洪区内洼淀的情况也发生了巨大变化，人口急剧增加，经济发展，财富增多。因此蓄滞洪区的合理调度变得尤为重要和困难，既要有效地拦蓄洪水保证城市的安全，又要最大程度地减少蓄滞洪区洼淀内的洪灾损失。为了实现实时的洪水调度分析，需要建立一个数值仿真系统，对洪水演进过程进行数值模拟，从而选择合理的洪水调度方案，合理、及时地调控滞洪区内各洼淀的使用，正确引导洪水，减少洪灾损失。

洪水演进数值的模拟是掌握洪水演进规律的重要方法之一。对于洪水演进数值模拟，国内外学者都进行了研究，意大利的 Caleffi 和 Valiani [1] 采用二维浅水方程对 Toce 河的洪水演进进行了模拟，我国刘树坤，周孝德，曹志芳等也对洪水演进进行了模拟 [2 ~6] 。滞洪区内铁路、公路、河流、堤防、水闸和各种防洪措施的使用，对适用于河道计算的圣维南方程组、整体求解的有限元或差分方法提出了新的要求。鉴于蓄滞洪区洪水演进具有干涸床面、不连续地面流和河道与洼淀多点水量交换的特征，采用无结构不规则网格的有限体积方法对模拟区域进行剖分计算，真实反映地形地物特点，能够更灵活地进行洪水计算和行洪控制，具有明显的优越性。

笔者综合考虑现有数学模型的优缺点，以一、二维浅水方程为基础进行数学建模，考虑阻水建筑物和动边界等复杂边界条件，采用逆有限体积法（FVM）进行离散和求解。同时本模型作为蓄滞洪区实时洪水调度仿真系统的重要组成部分，结合系统总体框架和调用，建立了完善的数值仿真系统，并在实践中进行了有效的应用和分析。

《2 控制方程组的离散求解》

2 控制方程组的离散求解

《2.1 基本控制方程》

2.1 基本控制方程

1）二维非恒定流基本方程

连续方程：

动量方程：

式中， h 为水深； z 为水位， z ＝z0 ＋h ； z0 为底高程； q 为源汇项； M，N 分别为 x，y 方向上的单宽流量，且 M ＝hu ， N ＝hv ； u，v 分别为 x，y 方向的平均流速； n 为糙率； g 为重力加速度。

2）一维非恒定流基本方程

连续方程：                 

动量方程：               

式中， Q 为截面流量； A 为计算断面的过水面积； Sf 为摩阻坡降。

《2.2 有限体积离散》

2.2 有限体积离散

按照有限体积法思想，采用无结构不规则网格布置方式，单元网格为控制体，在网格中心处计算水深 h ，在网格周边通道计算流量 Q 。即在平衡计算时，沿控制体每一边的法向通量用该边中点处的通量作代表，乘以边长即为通量沿该边的积分。中点的通量用中心格式（取相邻两格子形心处通量的平均）确定。应用水量平衡原理，控制体每一边的法向通量沿环路积分与控制体内的蓄水量平衡，周边通量取相邻控制体形心水位比降所形成的流量来确定。时间层面的计算采用单元、通道交错计算方式。

2.2.1 连续方程的离散

将式（1）改写成矢量形式：

按照有限体积法，将上式在控制体内进行积分：

假定水深 h 和源汇项 q 在同一网格内是均匀的，则式（7）可以写成

由高斯公式，则上式写成

式中， V 为网格周边通道上任意一点的速度矢量； n 为该点的外法线单位方向矢量。

令 Q ＝hV· n ， Q 为网格的各个通道的单宽流量，流入为正，流出为负，则对任一 K 边形网格有

把式（10）代入式（9）得

采用时间中心差分方式，式（11）化成

整理上式，并把时间向后移动一个时间步长。因此，方程（11）离散后得

式中， Ai 为第 i 个网格的单元面积； Lik 为 i 号网格的第 k 号通道的长度； Qik 为 i 号网格的第 k 号通道的单宽流量。

2.2.2 动量方程的离散

控制体周边通量计算依赖河道或滞洪区内各种地形条件如铁路、公路、桥涵、潜堤、堤防、水闸等建筑物情况，对不同类型的通道模型应用不同的理论模式进行概化，如地面型通道、河道型通道、连续堤防、缺口堤防、各类闸门等。通道主要模化为以下几种情况：

1）浅水地面型通道：即通道两侧单元为陆地且水深小于 0.5 m，通道上没有堤防等阻水建筑物。考虑到滞洪区内的地形起伏不大，地面洪水演进主要受到重力和阻力的作用，忽略掉加速度项，利用差分方法离散得到地面型通道的动量离散方程

式中， sign 为符号函数，表示 的正负与（）的正负相同； 为相邻单元形心到通道中心的距离之和； Hj 为通道上的平均水深。

2）深水地面型通道或宽河道型通道：对通道两侧单元为陆地且水深大于 0.5 m，或通道两侧网格均为较宽的河道型网格，通道为过流断面。动量方程沿通道法线方向离散为

式中， 为通道两侧单元的水深； 为通道两侧单元中心沿通道法线的投影速度。

3）窄河道型通道：对于滞洪区内宽度较小的河流，不便于将其划分成独立的单元网格，也不能忽略不计，为了方便计算，模型中将其模拟成具有深度和长宽尺寸的线性单元或称为特殊通道，以反映水流沿河而流以及河道与两侧陆地之间水量交换的现象，如果通道两侧有阻水建筑物，可以将其设为堤防。沿河道动量方程离散为

特殊通道与两侧网格之间的流量计算，采用宽顶堰流公式[7] ，即：

式中， m 为流量系数；σs 为淹没系数。

4）闸门型通道：行洪闸门通道的开启与关闭依赖于洪水调度条件，防潮闸的开启与关闭通过河道水位于海区潮位的差值来调节，闸门过流量由流量和闸上、闸下水位关系曲线确定。

《3 基本约束条件》

3 基本约束条件

《3.1 初始条件》

3.1 初始条件

初始条件分为两部分： a. 有水流流动的计算网格，根据实测资料如水深和河底高程，分析给出水位和流速；b. 没有水流流动的部分，可直接赋水深为零，流速为零。

《3.2 进出口边界条件》

3.2 进出口边界条件

计算域边界分为开边界、闭边界、防潮闸边界和溢流堰边界四类：a. 开边界是指计算域中与外界发生水体交换的那部分边界，边界的选取应满足（）· n ，在开边界处 q ＝q（t）或 z ＝z（t）给定，其中 q（t）和 z（t）为 t 时刻实测流量和水位。在蓄滞洪区一、二维洪水演进模拟中，对进水口门给出流量作为其边界条件，根据一维河网模型的计算结果来确定；对出水口门则可给出流量或水位过程作为边界条件。 b. 而闭边界即为实际存在的陆边界，在计算中可取水力要素法向偏导为 0。 c. 对入海口的出水口用防潮闸闸上、闸下水位差判断闸门开启关闭时机，当防潮闸上水位大于闸下水位时，闸门开启，否则关闭，实时给出流量或水位过程作为防潮闸边界条件。 d. 溢流堰边界是指计算域中与外界发生水体交换处有堰的那部分边界，或入海处的海堤边界，根据堰流公式计算结果来确定。

《3.3 内边界条件》

3.3 内边界条件

在蓄滞洪区内，考虑到铁路、公路、堤防及桥涵等对洪水演进速度影响较大，将他们概化成连续堤或缺口堤，形成有连续堤或缺口堤的特殊通道。对于这些特殊通道或几个蓄滞洪区之间需要通过分洪口门分洪时，其水流流动模拟一般不能按照圣维南方程组进行求解，可采用堰流和闸孔出流公式来描述[7] 。

《3.4 糙率系数》

3.4 糙率系数

糙率 n 是水流和河床相互作用过程中反映河道边界粗糙情况、河道形态等所有影响水流阻力因素的综合参数，其取值的正确与否直接影响到水力计算结果的精度。在蓄洪区内洪水演进模拟中，糙率需要根据区内计算单元的实际地形、地物情况和植被情况等，由实测资料分析并参考对比其他数据[7] 得到。

《4 模型的验证》

4 模型的验证

根据现有的“ 96 · 8 ”洪水实测资料，模拟洪水自大清河北支新盖房分洪道入东淀的全过程。据实际观测 8 月 5 日 17 时洪水开始下泄，洪水进入东淀后，推进速度相当缓慢，8 月 7 日 11 时到达新镇， 8 月 8 日晨 7 时到达苏桥，8 月 12 日 20 时到达廊大公路，8 月 13 日 20 时到达天津市境内的台头镇， 8 月 17 日 6 时到达第六堡。 1996 年 8 月份东淀的来水见图 1。

《图1》

图1 东淀来水示意图

Fig.1 Sketch of flood reaching time measured in Dongdian

模型从 8 月 5 日 17 时开始模拟，计算结果与实测资料对比见表 1。由计算结果可以看出，数值计算得到的时间有所滞后，这可能是由于糙率、渗流以及洼淀内给排水等不确定因素造成的。但从洪水传播速度看，计算结果的误差在可接受范围内，表明模型较好地模拟了东淀内道路等阻水建筑物、桥涵等过流建筑物，可以应用此模型模拟五洼洪水演进及其调度方案。

《表1》

表1 东淀计算与实测的洪水演进过程对比

Table 1 Comparison of flood routing between the calculated with the measured in Dongdian

《5 一、二维洪水演进数值仿真系统的建立》

5 一、二维洪水演进数值仿真系统的建立

一、二维洪水演进数值仿真系统是以蓄滞洪区一、二维洪水演进模型为核心，运用 GIS 技术，按照国家防总制定的蓄滞洪区洪水调度原则，结合洪水演进数值模拟、洪水风险区划、洪灾经济损失估算模型，实现空间和属性信息查询、分析、洪水演进过程数值模拟计算结果可视化、洪灾经济损失统计评估等各项应用。数据流程图见图 2。

《图2》

图2 数据流程图

Fig.2 Data flow diagram

其中一、二维洪水演进模型模块通过 VB 和 Fortran 混合编程来实现，系统设计界面如图 3 所示。主要是由数据输入（包括基本参数、蓄洪区数据、入水口门信息、时段及流量数据和网格属性等）、数据查看、计算和计算结果分析等 4 个子模块构成。

《图3》

图3 系统界面设计

Fig.3 Interface design of the system

《6 系统在实时调度中的应用及分析》

6 系统在实时调度中的应用及分析

天津市南部大清河流域位于海河流域的中部，是海河流域防洪的重点地区，在历次较大洪水总量比重上，都占 30 % ~50 %，是对天津市威胁最大的主要河系。由于大清河水系在天津防洪工作中举足轻重的位置，以及大清河水系存在着方方面面威胁天津市防洪安全的隐患，对大清河滞洪区三洼、五洼联合调度方案选择是否合理，在很大程度上将决定天津市的安全和整体防洪体系的合理性。应用该系统对大清河五洼滞洪区进行了一、二维河道及蓄滞洪区联合洪水演进数值计算，可获得在大清河发生不同频率洪水时，采用不同的调度方案和分洪措施后五洼滞洪区各自的洪水演进过程，实现实时的洪水调度分析，从而选择合理的洪水调度方案，正确引导洪水，减少洪灾损失。文章的洪水模拟分析依据 1994 年国家防办修改后的大清河洪水调度方案，同时结合《海河流域防洪规划》（2004 年 4 月）、《大清河系防洪规划报告》（2004 年 3 月）等相关规划，最终确定调度方案如图 4 所示。

《图4》

图4 大清河滞洪区五洼联合洪水调度图

Fig.4 Flood combined－regulation diagram of five basins in Daqing River basin

《6.1 基本计算条件》

6.1 基本计算条件

按照蓄滞洪区一、二维洪水演进模型网格的布置方式，高程系统采用国家 85 标准，在 1∶100 000 的地形图上，按照无结构网格的布置方式，形成一、二维衔接的河道型网格与洼淀地面型网格。其中，河道共有 148 个单元，405 个河道断面。模型中的河网、闸门、泵站分布示意如图 5。洼淀共有 5 760 个单元。洼淀模型及其上的闸门、泵站与特征位置分布见示意图 5。通过类比各类下垫面情况，确定相应的糙率值，取值在 0.035 ~0.09 之间。分别对 5 年、10 年、20 年、50 年、100 年、200 年 6 种标准洪水按现状独流减河过流2 000 m3 /s和设计独流减河过流 3 600 m3 /s 两种情况进行了模拟计算。来水过程历经 720 h，考虑退水，增加时间长度 380 h，模拟时间总长度为 1 000 h。退水过程采用自然退水，洪水从河口与泄洪道排泄。

《图5》

图5 河网模型、洼淀模型与泵闸及特征位置示意图

Fig.5 Sketch of river model，basin model， pump， brake and characteristic position

《6.2 计算结果分析》

6.2 计算结果分析

由于调度方案较多，计算成果较多，笔者只以 100 年一遇洪水为例，按现状独流减河过流 2 000 m3 /s 和设计独流减河过流 3 600 m3 /s 两种情况分别进行模拟计算，可获得洪水在滞洪区内的演进过程、传播速度及区内任一位置任一时刻的水位、流速、水深等成果。通过统计计算，可以进一步得到各洼淀的水位、流量过程曲线和整个滞洪区最大淹没水深的等值线图，见图 6 至图 8，通过对计算结果进行比较和分析可以实现洪水实时调度。

《图6》

图6 洼淀水位过程曲线

Fig.6 Water height curve of basin

《图7》

图7 洼淀蓄水量过程曲线

Fig.7 Flow curve of basin

《图8》

图8 滞洪区模型计算最大淹没水深分布

Fig.8 The max water －depth distribution of basin model

6.2.1 按独流减河现状 2 000 m3 /s

1）洪水演进过程。从零时刻起算，北支新盖房洪水 33 h 到达刘家铺；约 77 h 洪水到达老堤，同时，大清河在台头西向北扒口；98 h 洪水到达胜芳；由于独流减河现状 2 000 m3 /s 行洪能力的限制，从第 182 h 开始东淀大清河来水在进洪闸前壅高；第 219 h第六埠达到最高水位 6.89 m，独流减河进洪闸上最高水位达到 6.89 m，西河闸上达到最高水位 6.86 m；226 h 独流减河防潮闸上达到最高水位3.44 m，其最大泄量达 2 061 m3 /s；此后，由于上游来水减少，第六埠水位逐渐降低，在 758 h 东淀洪水退至4.0 m3 /s 以下。

南支枣林庄洪水在 25 h 到达王村闸；50 h 赵王新渠通过任庄子口门向东淀分洪；124 h 通过王村分洪闸向文安洼分洪； 148 h 小关开始分洪；第 198 h，东淀第六埠水位达到 6.44 m，启用滩里分洪口门，加大了向文安洼的分洪水量；209 h 东淀、文安洼均已利用且第六埠达到西河右堤超高 2.5 m 水位，达到了通过锅底分洪闸向贾口洼分洪的条件，锅底闸分洪后东淀第六堡水位明显回落，同时由于东淀上游来水逐渐减少，又经过 246 h，东淀水位低于文安洼水位，文安洼洪水经滩里扒口向东淀回流，再经锅底闸向贾口洼分洪，出现了“倒淀”现象，见图 9。

《图9》

图9 滩里及锅底分洪口流量过程

Fig.9 Flow progress of Tanli and Guodi brakes

2）退水过程。在 720 h 后新盖房、枣林庄、小关均停止下泄洪水，依赖现有河道向海区下游排洪。在 768 h，贾口洼梁头达到最高水位5.57 m；此后在独流减河进洪闸连续泄洪和新盖房无后续洪水的条件下，贾口洼水位逐渐回落，不满足利用津浦铁路 25 孔桥分洪的条件。在 1 000 h 东淀洪水已显著回落，但各洼淀仍有蓄水，特别是文安洼及贾口洼，其中文安洼大赵处水位仍高达 6.01 m，贾口洼梁头水位达 5.5 m。

6.2.2 按独流减河设计 3 600 m3 /s

在第 232 h 独流减河进洪闸达到最大流量 3 233 m3 /s，闸上最高水位达到 6.47 m；在 231 h 东淀第六埠达到最高水位 6.49 m；在第 236 h 独流减河防潮闸上达到最高水位 4.38 m，最大泄量 3 255 m3 /s。枣林庄下泄洪水在 124 h 通过王村分洪闸向文安洼分洪；148 h 小关开始分洪；在 227 h 由于东淀第六埠最高水位达到 6.44 m，启用滩里闸及扒口向文安洼分洪；此后由于东淀第六埠及文安洼水位没有达到锅底闸及扒口启用条件，因此未向贾口洼分洪；第 715 h 文安洼大赵处达到最高水位 6.16 m。

6.2.3 对比分析

100 年一遇洪水，在现状条件下需利用东淀、文安洼及贾口洼滞蓄洪水，而设计条件下仅需利用东淀及文安洼就够了。两种情况下计算结果对比情况如表 2 所示。由表 2 可以看出，由于独流减河现状 2 000 m3 /s 行洪能力的限制，相应加大了洼淀分洪运用的风险，增加了贾口洼的使用。现状条件下，从第 182 h 开始东淀大清河来水在进洪闸前壅高，使第六埠在第 219 h 达到最高水位 6.89 m，与设计情况相比增加 0.4 m；同时由于经独流减河入海水量有所减少，大量洪水都利用洼淀滞蓄，相应使东淀及文安洼的最大滞洪量都有所增加，特别是增加了贾口洼的运用，使其滞蓄了 7.78 ×108 m3 水，相应增加了贾口洼的洪水灾害损失。

《表2》

表2 100 年一遇洪水按设计、现状计算结果特征值对比

Table 2 Comparison of flood results between present condition and design condition of one in one hundred years

《7 结语》

7 结语

随着滞洪区经济地位的提高，合理地对蓄滞洪区内各洼淀进行联合调度，有必要建立一套功能强大、使用方便的实时洪水调度仿真系统。为了能够模拟蓄滞洪区内洪水演进过程，获得任一时刻、任一位置的水位、流量等信息，文章考虑多种约束条件，采用有限体积法进行了数学建模，并进行了模型验证，从而实现了一、二维洪水演进数值仿真系统的建立。应用该系统对大清河滞洪区五洼联合调度方案进行了实时模拟，并对不同方案进行了比较和分析，表明该系统可以用于滞洪区内洼淀联合调度，科学地进行防洪调度决策。