石家庄西南部外洪水动力模型构建与应用

张晓鹏，陈宇飞，杨艳玲，王会波，黄 腾，丁雪松

(1.河北天和咨询有限公司， 河北 石家庄 050021； 2.华北水利水电大学， 河南 郑州 450000；3.河北省水利水电第二勘测设计研究院， 河北 石家庄 050021)

石家庄市地处河北省中南部、华北平原腹地，是河北省省会和政治、经济、科技、文化中心。石家庄市历史上遭遇的洪水主要来源西南部山洪、西北部滹沱河及太平河洪水过境，防洪体系不健全导致洪水漫溢致灾，其中，西南部山洪主要来源于方台沟、台头沟、金河、大宋楼沟、普莲河及洨河。有历史记载的洪灾主要发生在汛期7月下旬至8月上旬，较为严重的有1939年、1954年、1956年、1963年和1996年大洪水，其中1963年8月洪水最为严重，石家庄地区普降暴雨，致灾严重。伴随着石家庄经济社会发展，石家庄城区不断扩大，原抵御外部洪水的行洪河道逐渐被城市建成区所包围，行洪滞洪区域减小，石家庄市防洪抢险任务加剧。因此，在石家庄现有防洪体系的基础上，建立石家庄外洪水动力模型，对石家庄外洪发生、洪水演进进行仿真模拟，对于指导石家庄防洪抢险，避险转移具有重要意义。

洪水模型主要包含有水文学模型和水动力学模型[1]两类。水文学模型类似黑箱或灰箱系统，以系统分析为途径，以水量平衡为依据，模拟水文汇流过程，建立输入与输出的相关关系；水动力模型则是在微观物理定律为基础上，模拟坡面汇流过程。其中水动力模型应用较为广泛，主要用于模拟洪水演进[2-4]、城市内涝[5]、河道泥沙输移[6]和水环境水污染[7]等等，模拟区域包括河道行洪区及漫溢区，城市、蓄滞洪区等等，如：王扬等[8]对韩江南北堤防溃决后对堤外防洪保护区洪水演进进行模拟；李大鸣等[9]对桃林口水库下游山洪沟洪水演进进行模拟，并进行洪灾评估分析；袁水龙等[10]采用MIKE耦合模型模拟淤地坝对小流域暴雨洪水过程的影响；张靖雨等[11]和田福昌等[12]采用洪水内涝模型耦合模拟防洪保护区在堤防溃决后的洪水与内涝叠加过程，进行成灾分析；倪丽丽等[13]对城市内涝，通过精细化处理进行模拟，为城市进行暴雨预警，风险评估提供了经验；戚蓝等[14]针对沿海城市三亚市，建立了雨、洪 、潮多元耦合精细化洪涝分析数值模型，为沿海城市防汛抢险提供了技术支撑；陈平等[15]和杨艳玲等[16]基于MIKE21数值模拟蓄滞洪区洪水演进过程，为防洪风险图制作提供依据。

本文采用MIKE系列软件，构建了石家庄西南部外洪水动力学模型，模拟了新防洪体系下石家庄西南部洪水演进的动态过程，评估了石家庄西南部洪水对于石家庄主城区的影响程度，为防汛预警、风险识别、进行防洪抢险等提供了依据。

1 模型原理

根据石家庄市主城区所在地理位置及历史水灾等具体情况，确定石家庄市主城区洪水分析计算采用水力学法。采用一维、二维，并将二者耦合的水动力学模型对石家庄主城区洪水进行分析计算。

(1) 一维水动力模型。一维水动力计算模型基于垂向积分的物质和动量守恒方程，即利用圣维南方程组来模拟河流或河口的水流状态。方程形式[17]如下：

(1)

(2)

式中：x为计算点空间坐标、t为计算点时间坐标；A为过水断面面积，m2；Q为过流流量，m3/s；h为洪水水位，m；q为侧向流入流量，m3/s；C为谢才系数；R为水力半径；α为动量校正系数；g为重力加速度，m/s。

(2) 二维水动力模型。二维水动力模型的控制方程为二维浅水方程，包括连续性方程和动量方程。模型可以基于笛卡尔坐标系或球面坐标系，数值计算方法可采用基于非结构网格的有限体积法，即不规则三角形网格法，可提高计算速度，同时在拟合复杂地形方面优点突出。方程形式[17-19]如下：

连续性方程：

(3)

x方向动量方程：

(4)

y方向动量方程：

(5)

式中：h为水深，h=d+ζ，其中ζ、d分别为水位和水深，m；p、q分别为x方向和y方向单宽流量，m3/s；g为重力加速度，m/s；C为谢才系数；Ω为科氏力系数；f为风阻力系数；ρω为水密度，kg/m3；V为风速，m/s；Vx、Vy为风速在x、y方向分量，m/s。

2 模型建立

2.1 模拟区域概况

石家庄西南部洪水来源于为洨河流域，包括方台沟、台头沟、大宋楼坡水区、金河、普莲河、洨河和洨河南支，各河沟洪水先后在石家庄市西南防洪堤外侧汇入南泄洪渠、洨河后，穿越京广铁路桥，向东南泄向下游，主要威胁石家庄城区西南部区域，如图1所示。本文模型模拟区域北部以石获南路为界，西以南水北调右堤为界，南部跨越青银高速公路以洨河、北沙河分水岭为界，东至原京港澳高速公路，计算范围289 km2。

2.2 溃口设置

2.2.1 溃口位置确定

溃口位置确定主要考虑以下因素：历史洪水曾经决口的位置、现状堤顶高程较低处、堤防质量薄弱处、河势顶冲处、行政区划分情况、重点保护对象情况、涵闸等。本文综合考虑上述因素确定：南泄洪渠左堤南泄洪渠与洨河汇合口处拟定溃口—东同平溃口。

图1 石家庄市西南部河流水系图

2.2.2 溃决方式和形态确定

(1) 决口方式。堤防决口方式有漫决、冲决、溃决。漫决是指水位超过堤顶，洪水漫溢造成堤防决口；冲决和溃决是指水位尚未到堤顶，由于淘刷、堤防质量、串漏等原因水流冲垮堤防造成决口。本文决口采用溃决方式。西南防洪堤口门溃决时机：此段西南防洪堤断面较为单薄，口门溃决时机为达到50年一遇洪水位。

(2) 决口形态。根据研究区情况，综合考虑实地观测、专家经验及经验公式等方面因素，确定溃口宽度。南防洪堤东同平溃口宽度为200 m，溃口形状等效为矩形。

2.2.3 溃决、漫溢流量计算方法

(1) 堤防溃决流量计算。一维模型中的溃口公式采用基于美国国家气象局开发的溃坝模型(NWS DAMBREAK)来进行计算。溃坝模型(NWS DAMBREAK)[18]由美国工程师弗莱德(Fread)在1988年研制，该模型在国内外应用较为广泛，证明其适应性较强，结构性能较好。模型中，溃口采用堰流公式(Breach Failure)来进行计算。公式形式如下：

(6)

式中：b为溃口底宽，m；h为上游洪水位，m；hb为溃口底高程，m；s为溃口边坡比；cweir为溃口水平堰系数，等于0.546 430；cslope为溃口边坡部分堰系数，等于0.431 856；Cv为入流收缩损失修正系数；ks为淹没修正系数。

(7)

(8)

式中：cB为无量纲数，等于0.740 256；WR为水库宽度由未破坏的坝顶长度来定义；hb,term为决口最终底高程，其值为时间序列文件中的最低高程，m；Qp为上一次迭代中决口的过流量，m3/s；hds为下游水位，m。

本文计算根据一维河道的水位来决定起溃时间，当河道水位超过设计水位时，堤防开始溃决，溃决方式为瞬溃，溃底高程为河滩地高程。

(2) 漫溢流量计算。模型中一二维模型侧向连接、铁路、公路的漫溢流量采用Villemonte Formula堰流公式进行计算：

(9)

式中：q为通过结构物的流量;W为宽度;C为堰流系数;k为堰的指数；Hus为上游的水位；Hds为下游的水位；Hw为堰顶高程。

2.3 模型构建与参数确定

石家庄西南部外洪水动力模型由一维模型和二维模型共同构成。南泄洪渠渠道行洪采用一维模型进行模拟；方台沟、台头沟等河沟由于河道较小，多数河道漫溢出槽，城区较平坦，易串流，因此采用二维模型进行模拟，并通过MIKEFLOOD进行耦合，共同构成石家庄西南部外洪水动力模型。

2.3.1 南泄洪渠一维模型构建

(1) 上、下边界。上边界为方台沟设计洪水过程线，下边界为与二维模型进行耦合的哑元水位边界。

(2) 断面概化。南泄洪渠上游断面摘自《石家庄市区西北部水利防洪生态工程》，桩号3+200—16+220(五支渠汇入口)，断面间距200 m，跌水上下游及工程位置处断面进行加密处理，建立南泄洪渠一维模型。

2.3.2 二维模型构建

(1) 网格剖分。采用不规则三角形网格进行剖分，为了提高计算精度，最大网格剖分面积不大于0.01 km2，见图2。网格剖分后，采用1∶1000实测的地形图建立的DEM高程数据进行网格高程插值，不足部分采用1∶10000地形图进行插补。

(2) 边界处理。模型上边界条件为洨河、洨河南支、普莲河、金河、大宋楼、台头沟、方台沟设计洪水过程线，下边界为青银高速洨河大桥，按公路桥断面推算的洨河水位流量关系进行控制。

图2 模拟区域网格剖分及示意图

(3)地物处理。模型主要阻水地物有货运铁路编组站、铁路、高速公路、内河堤防等，建立二维模型时，除编组站因宽度较大、建基面较高直接扣除模型外，其余阻水地物由于宽度较窄，本文模型均概化为堰，水位低于阻水地物顶部高程时不过水，水位高于顶部高程时按堰流公式计算过流量。对于公路桥涵等过水建筑物，模型按其宽度将道路等线状地物分段处理。

2.3.3 参数确定

(1) 一维模型参数确定。南泄洪渠河槽为浆砌石单式断面，糙率采用0.016。计算时间步长取15 s。因南泄洪渠为环城水系的一部分，有景观蓄水要求，因此初始条件按水深0.5 m考虑。

(2) 二维模型参数确定。在二维模型中，需设置的主要参数有：糙率、计算时间、计算步长、干边界、湿边界、水密度、风场、涡黏系数等。洪水模拟的目的是计算洪水在淹没区的演进过程，提取淹没水深、淹没水位、流场等要素，对于计算结果无明显影响的计算参数采用默认值，包括水密度、风场、涡黏系数等，其余参数设置如下：

① 糙率。在整个计算区域进行糙率分区处理。结合计算区域土地利用图对计算区内糙率进行取值：村庄、建筑用地采用0.100，农田、绿地采用0.045，河流采用0.03。

② 时间步长。时间步长的设定是为保证提高模型运行效率，并保持模型的稳定性，本次最大时间步长设定为15 s，最小时间步长为 0.01 s，模型依据网格剖分质量、洪水演进情况、区域地形等自动调整运行时间步长。

③ 干湿边界。干湿边界为避免模型计算中出现不稳定性和不收敛而设定的参数，本次干水深取值0.005 m，湿水深取值0.1 m，淹没水深取值0.05 m。

2.3.4 模型率定

利用《石家庄市西北部水利防洪生态工程(二期)初步设计报告》成果对一维模型进行率定，水位相差幅度在0.02 m～0.08 m之间，满足精度要求。

3 计算方案

按照现状堤防、河道断面情况分析南泄洪渠与洨河干支流洪水对石家庄城区的淹没影响。计算方案为现状条件下，设置东平同溃口，南泄洪渠及洨河分别发生50年、100年及200年一遇洪水时的淹没情况，共3种工况。

4 模拟分析

洨河流域遭遇50年一遇洪水，溃口断面最高洪水位65.18 m，溃口溃决分洪，洪水沿低洼地段进入二维计算区域。由于洨河流域上游支流河道方台沟、台头沟、大宋楼沟、金河、普莲河、洨河及洨河南支现状河道断面较小，洪水演进时，除方台沟直接入南泄洪渠以外，其余河道洪水均漫溢出槽；洪水演进2 h后，金河洪水在货运编组站西北附近分为南北两股，分别沿货运编组站南、北两侧向东南方向演进；洪水演进4 h后，普莲河洪水在青银高速公路以西沿青银高速向东南扩散，与其他河沟洪水发生串流；洪水演进7 h后，南泄洪渠洪水在部分河段向右岸漫溢，各支流洪水漫溢后逐渐连片；洪水演进10 h后，洨河溃口开始分洪，分洪后洪水向东演进，由于受京广铁路路基影响，向南北两侧持续扩散；与此同时，洨河洪水在铁路桥以东漫过右堤呈扇形向东扩散；洪水演进11 h后，溃口洪水向北穿过南三环桥梁演进至主城区；洪水演进14 h后，淹没面积持续增大，并趋于稳定。洨河流域遭遇100年及200年一遇洪水演进过程与50年一遇洪水相似，仅淹没区域及水深有少许增加，溃口处河道最高洪水位分别为65.45 m、65.62 m。 50年、100年、200年一遇洪水淹没主城区面积0.29 km2、0.34 km2、0.38 km2，水深较大处集中在京广铁路桥附近，最大水深超过2 m；流速较大处多位于河流弯道及溃口位置，流速大于1.5 m/s。模拟主要成果如图3—图5所示(只列出50年一遇洪水情况)。

图3 现状50年一遇水位图

图4 现状50年一遇水深图

图5 现状50年一遇流场图

上述模拟结果表明，石家庄西南部洨河流域发生50年至200年一遇洪水时，南泄洪渠以西区域，河道过流能力有限，防洪能力不足，洪水发生漫溢。南泄洪渠左堤作为保护石家庄主城区的主要屏障，五支渠以上段左岸为西三环兼防洪堤，标准较高，洪水不会通过几座较大的穿堤桥梁进入市区。但沿线部分穿堤涵洞地势较低，洪水时需要及时封堵。五支渠以下段是堤防的薄弱环节，一旦溃堤，洪水将沿京广铁路桥向北，向东进入城区，威胁石家庄主城区安全，建议按照200年一遇防洪标准，尽快实施五支渠口以下段西南防洪堤除险加固工程，确保主城区防洪安全。

5 结 语

本文针对石家庄西南部洪水特点，采用MIKE系列软件，构建了洪水水动力模型，并对石家庄西南部区域在现有防洪体系情况下，在洨河流域遭遇50年、100年及200年一遇洪水时，对洪水演进及淹没情况进行模拟，发现问题并提出建议。

石家庄西南部洪水动力模型的建立可为防汛抢险部门进行洪水风险识别、发布洪水风险预警、进行防洪抢险、避险转移提供决策依据，同时也可为土地规划、国土开发提供依据。