城市洪涝水文水动力耦合模型构建与评估

2021-06-08黄国如陈文杰喻海军

水科学进展 2021年3期

黄国如,陈文杰,喻海军

(1.华南理工大学土木与交通学院,广东广州 510640；2.华南理工大学亚热带建筑科学国家重点实验室,广东广州 510640；3.中国水利水电科学研究院,北京 100038)

随着全球气候变化加剧以及城市化迅速发展,中国城市洪涝问题十分突出,为减少洪涝灾害对城市区域的影响,亟需建立一套高效稳定的城市洪涝模型,对城市区域的暴雨洪涝过程进行模拟,为城市防洪排涝、抢险救灾等提供决策依据[1- 3]。城市雨洪模拟方法主要可归结为水文学方法、水动力学方法以及水文水动力学方法。水文学方法是最早被采用的方法,其基于水文学原理进行流域产汇流计算,结构简单,效率较高,但仅能得到流域出口处的流量过程,无法给出特定位置的水力特征要素；水动力学方法基于网格单元进行地表产流计算,分别利用二维浅水方程和圣维南方程进行地表及排水管网水流汇流计算,该方法计算精度较高,但计算效率较低；而水文水动力学方法集合了上述2种方法的优势,其以子流域为水文响应单元计算地表产汇流,水流进入地下管网后采用水动力学方法进行计算,溢流到地表的水流则采用二维浅水方程进行计算,该方法具有较好的水文学基础,同时计算效率较高[4- 6]。

针对一、二维模型的耦合方式,国内外已有学者研究提出了一些较为符合实际的耦合理论并构建了耦合模型[7- 9]。对于管网模型与地表二维模型耦合,主要存在垂向水流交换,即水流在检查井、雨水篦处的交换。早期Hsu等[10]将SWMM模型与地表二维模型耦合,模型间水流交换是单向的,当管网能力超载时,水流从检查井处溢出到地表流动；当地下管网有足够的过流能力时,地表水流也不能通过检查井重新回流到管网中。随着研究的不断深入以及技术手段的不断创新,逐渐出现一批能够实现地下地表水流双向交换的耦合模型。Seyoum等[11]将SWMM模型与自主开发的二维地表模型耦合,克服了地表和地下管网水流交换的困难,实现了真正意义上的一、二维水流交换。目前最常见的垂向交换水量计算方法主要为堰流公式和孔口流量公式[12- 14],通过泄流试验比较各种雨水口泄流计算方法的优缺点及适用范围,为垂向交换水量计算方法提供基础[15- 16]。对于一维河网模型与二维模型耦合,其水流交换主要存在水平方向上的正向交换和侧向交换,目前还大多采用堰流公式、互为提供边界法和水量动量平衡法等进行计算[17- 18]。另外,近期国内外学者鼓励在洪涝数值模型的率定和验证中,除了采用通过现场监测等传统途径获取得到管道流量、检查井液位和地表内涝深度等数据集以外,通过电视新闻、微博图片等社交媒体途径获取得到的数据集也可以用于增加与扩展现有的数据集[19- 20]。喻海军等[21]将双时间步法拓展应用于二维浅水模拟中,构建的二维水动力学模型具有良好的精度和可靠性[18]。总的来说,模型的耦合方式与计算方法很多,本研究旨在前人研究的基础上提出一、二维模型耦合方式,为城市暴雨洪涝的“多维和多过程”模拟提供稳健的耦合策略。

本文将一维SWMM模型与二维水动力模型[18,21]进行耦合构建城市洪涝水文水动力耦合模型IHUM(Integrated Hydrology and Hydrodynamics Urban Flood Model),利用传统数据和社交媒体等数据对该模型进行验证,分析该模型在城市流域的适用性,阐述模型的可靠性和精度。

1 城市洪涝水文水动力耦合模型

1.1 地表产汇流模型

在数值模拟中,地表产汇流计算是其中一个重要环节,其计算精度直接影响后续水流在管网汇流以及节点溢流计算的准确性。将计算区域分成若干个子汇水区,根据子汇水区特性单独计算净雨和地表汇流,并假定某个子汇水区产流量均流入某个特定的排水管网节点。将子汇水区分为透水部分、有洼蓄量的不透水部分和无洼蓄量的不透水部分分别计算产流和汇流,利用SWMM模型计算其产汇流过程,具体计算过程参见文献[22]。

1.2 管网或河道一维模型

SWMM模型擅长计算排水管道和渠道等排水设施的水流状态[22],与其他河道水流模型相比,SWMM模型对河道水流模拟存在一些不足之处。在河道模拟模型中,模型假设一个河段内的断面形状沿程变化,而SWMM模型则假设一个河段内的断面形状相同,不同河段之间可以具有不同形状。显然,在断面多变的天然河流模拟中,河道模拟模型的概化方法更符合实际情况,SWMM模型的概化方法将会产生较大误差。但城市河流大多以渠道和整治过的河涌为主,断面形状较为规则,在此情形下,两者概化方法之间存在的差异就会变得很小,意味着SWMM模型在对河段断面概化中所产生的误差可以忽略不计。同时考虑到SWMM模型具有能够处理各种水工建筑物(例如泵站、水闸和堰等)的优势,本文认为采用SWMM模型计算城市河网和地下排水管网的水流过程是较好的解决方法。

1.3 地表二维水动力学模型

二维水动力模型基于二维浅水方程,采用非结构网格中心型的有限体积法,结合Godunov型格式计算界面通量,将分片线性逼近的MUSCL格式应用于变量重构,建立了一个时空均具有二阶精度的二维非恒定流水动力学模型。偏微分方程的数值求解常常受限于显格式小时间步长以及隐格式计算量大等问题,而隐式双时间步法很好地解决了这2个问题,双时间步法自提出后在空气动力学方面得到了广泛应用并取得了良好效果[23- 25]。喻海军[18]将双时间步法拓展应用于二维浅水模拟中,应用双时间步法建立隐式高效的有限体积法数学模型,并采用一系列算例对二维水动力模型处理恒定流、间断流、动边界以及实际地形等方面的能力进行了检验与验证,将隐式双时间步法与显式方法进行对比,结果表明二维模型具有良好的精度和收剑速度,能够有效地处理干湿边界问题、满足和谐性要求和处理复杂水流和实际地形,模型算法的鲁棒性较好,稳定性和效率都较显式格式有了大幅度提升,模型具有处理实际洪水的能力以及良好的计算精度,能够应用于城市洪涝数值模拟计算[21]。

1.4 城市雨洪耦合模型

从水流在城市排水系统中的流动可知,雨水或者洪水在城市内的流动较为复杂,包含了地表、地下排水管网以及河道的水流流动,涉及到地表、地下排水管网和河道三者之间复杂的水流交换机制。

一维水动力学模型具有求解相对简单、计算效率较高以及需要的建模资料相对较少等优点,被广泛地应用于城市地下排水管网和河道水流模拟计算。二维水动力学模型虽然具有结果详细、擅长处理方向不定水流模拟的优势,但模型仍然具有计算过程复杂、计算效率低、对建模资料要求高和不便于概化水工建筑物等不足。而城市区域的水流过程复杂,既有管网、河道、街道等一维属性明显的水流流动,又有地表、街道交叉口等二维属性明显的水流流动,无论采用何种模型对城市区域水流进行模拟都会涉及到模型的不足之处。因此,将一维与二维模型进行耦合,发挥各自的优势对城市水流进行模拟,应用于各自最适合的情形是一个很好、甚至是必然的选择。据此,构建城市洪涝水文水动力耦合模型,其模型结构如图1所示。

图1 城市洪涝水文水动力耦合模型Fig.1 Integrated hydrology and hydrodynamics urban flood model (IHUM)

2 城市雨洪模型耦合方式

城市雨洪模型耦合方式主要涉及到垂向耦合和水平耦合2种方式[17]。

2.1 垂直方向耦合

一、二维模型在垂直方向的耦合连接,即一、二维模型地表地下耦合,主要是针对城市地下排水管网与地表的水流交换问题。在模型模拟中,节点是地下排水管网与地面水流交换的唯一通道,垂向耦合则是计算在节点处发生的水流交换现象。假设节点水头为H1D,与该节点对应的地表网格水位为H2D,根据两者之间的关系将垂直方向的水流交换分为3种情况:①H1D>H2D,此时管网系统中水流通过节点溢出到地表流动；②H1D

由于目前缺乏对垂向连接水流交换机理和计算方法的深入研究,使得垂向连接水流交换的基础理论相对不成熟,计算方法也极为有限,一般都采用堰流公式和孔口流量公式计算垂向连接处的交换水量。因此,本文也采用堰流公式和孔口流量公式计算垂向交换水量。

2.1.1 节点溢流

采用孔口流量公式方法计算溢流量,且考虑该时刻管网的水流状态,具体计算步骤如下:

(1) 假定节点允许溢流,溢流储存面积设置为检查井面积,计算节点水头H1D,通过二维模型模拟结果获取节点对应位置的地表网格水位H2D,利用式(1)计算交换水量。

(1)

式中:Qn,s为计算得到的当前时间步的溢流量,m3/s；co为孔口流量系数,取值范围为[0,1]；Amh为节点蓄水面积,m2；g为重力加速度,m/s2。

(2) 根据SWMM模型的计算结果获取该节点的总蓄水量,限制其在下一时间步的溢流量不能超过该节点的总蓄水量。同时,需设定节点最大允许流量Qem,具体可根据实际情况进行设置。因此,为保证模型稳定性,采用式(2)对模型的溢流量进行限制。

Qn,s=min(Vmh/tn+1,Qem,Qn,s)

(2)

式中:Vmh为节点蓄水量,m3；tn+1为下一步时间步长。

(3) 将溢流量作为SWMM模型节点的外部出流以及二维网格单元的源项,计算下一时间步的结果。

2.1.2 节点回流

采用堰流公式和孔口流量公式计算回流量,具体计算步骤如下:

(1) 通过SWMM模拟结果以及二维模型模拟结果获取节点水头值H1D以及对应位置网格水位值H2D,根据地表与节点的水位差,采用式(3)计算回流量:

(3)

式中:cw为堰流流量系数,取值范围为[0,1]；Qs,n为节点回流量,m3/s；h2D为地表水深,m；w为节点周长或者雨水口宽度,m；Z2D为地面高程,m。

(2) 为保证模型的稳定性,采用式(4)对地表回流量进行限制:

Qs,n=min(Qs,n,Qem,V/tn+1)

(4)

式中:V为与检查井连接的单元格的总水量,m3。

(3) 将回流量作为节点外部入流输入到SWMM模型中,同时也将回流量作为源项代入二维模型中,更新至下一时间步。

2.2 水平方向耦合

针对河道与地面的水流交换问题,采用一、二维模型对水平方向的耦合连接进行处理,即地表一、二维模型耦合,根据耦合位置以及耦合处水流流向等特征的差异,将地表一、二维模型耦合分为正向连接和侧向连接(图2)。正向连接是指河道通过上下游与二维区域连接,耦合位置位于河道两端,水流通过河道两端与二维计算区域进行水流交换,连接处水流方向与河道中水流方向一致,涉及到河道上下游的边界条件问题；侧向连接则是指河道通过两岸与二维区域连接,耦合位置位于河道两岸,水流从河道两岸流向二维区域或者从二维模型计算区域经由两岸流入河道,连接处水流方向与河道水流方向不一致,通常会形成一定的夹角,不涉及河道上下游边界条件问题,因此,2种不同的连接方式通常需要采取不同的连接策略和计算方法。分别采用互为提供边界法和堰流公式法来计算正向连接水流交换和侧向连接水流交换。

图2 地表一、二维模型水平方向耦合示意Fig.2 Horizontal coupling of one and two dimensional surface models

2.2.1 正向连接

正向连接即水流通过河道两端与二维区域进行交换,采用一、二维模型互为对方提供边界条件的方式计算正向连接,该方法关键在于确定连接处2个模型的边界条件,具体思路和步骤如下:

(1) 二维模型以SWMM的下游出流量作为流量边界条件,即:

(5)

式中:Q1D,n为河道与二维区域连接断面的流量,m3/s；M为二维区域与河道连接的单元边数目；lk为k单元边的边长,m；qk为k单元边的单宽流量,m2/s。

根据二维模型处理流量边界的方法,采用曼宁公式将流量分配到连接处的每个单元边上:

(6)

(2) SWMM模型以二维模型的水位作为下边界条件,采用连接处单元的加权平均水位作为SWMM模型的水位边界条件:

(7)

式中:Z1D,n+1为河道下一时间步水位边界条件,m；zk,n为网格水位值,m；L为正向连接边界的总长度,m。

2.2.2 侧向连接

侧向连接即水流通过河道两岸与二维模型计算区域连接处进行交换。河道在一维模型中进行模拟,计算单元为河段；地表在二维模型中进行模拟,计算单元为网格。侧向连接水量交换计算中不涉及边界条件问题,因此,二维模型与河道连接处设置为固边界条件。侧向连接计算具体步骤如下:

(1) 在计算水量交换时,要事先设置好一、二维模型各自的连接单元,通常以二维区域中与河道相连的网格为单位分别进行计算。关键要求出当前时间步与河道相连网格的水位Hc以及网格对应位置河道的水位Hr。网格水位Hc可通过二维模型的计算结果获取,对应位置河道水位则需要在SWMM模型模拟结果中插值获取。SWMM无法提供管道中任意位置的水位结果,为获取与河道相连网格对应位置的河道水位就必须先从SWMM结果中提取河段上下游节点水位,然后采用线性插值计算获取。

(2) 根据不同的地表网格水位Hc以及对应位置的河道水位Hr组合计算侧向连接交换流量。存在以下4种情况:①当Hc和Hr均小于河道堤顶高程Ze时,不发生水流侧向交换,即Q=0；②当Hc>Hr且max(Hc,Hr)>Ze时,发生水流侧向交换,水流方向为从二维区域到河道；③当HcZe时,发生水流侧向交换,水流方向为从河道到二维区域；④当Hc=Hr>Ze时,在实际情况下应根据地表和河道的水流方向来判定水流的方向,但在采用堰流公式计算时,公式中并没有考虑两边水流的初始流速,因此Q=0。采用堰流公式近似计算交换流量Q的方法如下:

(8)

(9)

式中:Hr和Hc分别为堰上、下游水位,分别取河道和二维网格单元的水位值,m；Ze为堰的高程,可根据实际情况设置,m；be为堰的宽度,一般取单元格与河道相连边的边长,m。

(3) 将采用堰流公式计算出来的流量Q作为源项加入到对应网格计算中,同时将流量Q作为管道水量损失(情况③)或者水量增加(情况②)添加到SWMM模型中。

2.3 模型耦合实现方式

耦合模型目的是以水量交换为纽带,通过调用动态链接库文件的方式调用SWMM模型计算引擎进行计算。通过分析SWMM模型代码可以发现,模型提供了9个可被外部程序调用的外部函数,函数可通过dll文件直接调用；而SWMM模型、二维模型和耦合模型均采用C语言编写,因此在不同模型之间的函数调用时无需语言转换设置,所有函数均可在Visual Studio 2010平台使用C语言直接调用。

但仅依靠上述9个外部函数并不能实现上述的一、二维模型耦合计算,还需增加一些外部函数来提取SWMM模型计算过程结果、设置新的时间步以及输入交换的水量。因此,在编译dll文件之前需在SWMM模型中对部分代码进行修改,增加一些外部函数来实现一、二维模型之间数据的通讯,使其更便于一、二维模型耦合。通过调用上述新增的外部函数以及原来的外部函数就可以实现一、二维耦合模型构建,具体参见文献[17]。

3 理论案例验证

为验证所构建的城市洪涝模型的可靠性,采用理论案例和实际例子进行分析评估。在理论案例中,模拟分析了垂向连接和水平连接2种方法的结果,但受篇幅所限,仅给出了垂向连接算例结果。

垂向连接方式如图3所示[17],管网系统由6个节点和6根管道组成,地表区域为1个边长为200 m的方形闭合平原区域,平原区域覆盖节点2、节点3、节点4和节点5。平原区域覆盖的节点均允许溢流或者回流,管网可以通过这4个节点与平原地表进行水流交换,平原高程为0 m,糙率为0.025。管道与节点的连接方式采用管底平接,即管道底端与检查井底部相接。水流通过节点1进入系统,入流量在模拟开始的10 min 内由0逐渐增加至1.0 m3/s,然后保持不变；水流在系统中流动,最后通过节点6离开系统,出口设置为自由出流边界。系统初始为干,即管道和平原的初始水深均为0 m。

图3 垂向连接算例示意Fig.3 Example of vertical connection

图4为管网出水口处的流量过程,从中可以看出,出水口流量在模拟开始时变化较大,从0上涨到0.86 m3/s,在之后很长一段时间内流量的变化速度逐渐减慢,在48 h之后流量上涨速度降为0,即流量已经基本趋于稳定状态,最后稳定在1.0 m3/s。

图4 管网出水口处流量过程线Fig.4 Discharge hydrograph at the outlet of pipe network

采用耦合模型对本算例进行建模计算,模拟时长为48 h。为验证本模型计算结果,采用城市排水商业软件InfoWorks ICM进行对比,InfoWorks ICM为一款将城市排水管网及河道一维模型与地表二维洪涝淹没模型结合在一起的城市排水软件,可以较为真实地模拟地下排水管网系统与地表受纳水体之间的相互作用,可建立与本文耦合模型类似的城市排水模型。

表1为系统处于稳定状态时的管道流量和节点水位,得知:①若节点2无溢流或回流,则管道2与管道3的流量之和应与管道1的流量相等。但表1显示管道2和管道3流量之和小于管道1流量,说明节点2存在溢流,溢流量为0.138 m3/s。②若节点3和节点4无溢流或回流,则管道3流量应与管道5流量相等,管道2流量应与管道4流量相等。但表1显示管道5流量大于管道3流量,说明节点3存在回流,回流量为0.046 m3/s,同理节点4处的回流量亦为0.046 m3/s。③若节点5无溢流或回流,则管道4与管道5的流量之和应与管道6的流量相等。但表1显示管道4和管道5流量之和小于管道6流量,表明节点5存在节点回流,回流量为0.046 m3/s。④纵观整个排水系统,节点2溢流量刚好与节点3、节点4和节点5的回流量之和相等,使整个系统处于稳定状态。

表1 稳定状态时模拟得到的管道流量和节点水位

将管网系统自身属性与模型模拟结果结合分析可知,由于管道2和管道3管径较小,无法及时排走从管道1中流过来的雨水,导致节点2的水位迅速抬升,水流从节点2溢出在地面上流动。而管道4和管道5的管径比管道2和管道3的管径大,管道4和管道5除了能够及时排走管道2和管道3流过来的水量之外,还可以接受地表水量的回流,因此部分水流经节点3和节点4重新汇入管道。而管道6管径与管道1管径一样,排水能力可达到1 m3/s,因此部分水流经节点5重新汇入管道。从系统整体进行分析,在达到平衡状态之后,管网节点溢流量与回流量相等,管网入流量与出流量也相等,管道顶部平原区域中水流处于平衡状态,模型结果总体上符合规律,具有较好的合理性。

本案例研究以InfoWorks ICM的模拟结果为标准,从表1可以看出,本文模型计算结果与InfoWorks ICM计算结果基本吻合,表明本文建立的耦合模型的计算结果是可靠的。2个模型计算得到的管道流量最大相对误差为4.36%,节点水头相差更小,最大不到2.5%,这些差别可能是由于一、二维模型和地表地下耦合算法以及划分和采用的网格不同所造成。总的来说,本模型和InfoWorks ICM软件的计算结果在管道流量和节点水头上的差异均较小,处于可接受的合理范围之内,模型具有较高的可靠性。

4 东濠涌流域验证

4.1 基本概况

东濠涌流域位于广州市越秀区中部,流域属人口密集的老城区,排水系统采用合流制排水管网,流域面积为10.38 km2,受亚热带季风气候影响,流域年均降雨量为1 720.0 mm。东濠涌出口为珠江,为防止珠江河水倒灌,在出口处建设了东濠涌闸门；另外,为了在珠江高水位时仍能进行排水,在东濠涌出口处建设了东濠涌泵站,泵站能力为80.1 m3/s。由于该区域降雨量大以及易受珠江水位顶托影响,该区域常遭受城市洪涝影响。

4.2 基础资料

根据地形、建筑物和排水管道分布,整个研究区域划分成为7 700个子汇水区,一维模型中包含3 928条管道和明渠、3 471个检查井、1个排水泵站和1个出水口。模型计算网格为非结构网格,且流域内的建筑物区域被认为是不积水区域,为保证在模拟时建筑物区域不积水,建筑物在划分网格时被排除在网格区域之外。为获取建筑物轮廓,采用高精度的遥感图进行划分。东濠涌被设置为不积水区域,河涌堤岸被设置为侧向耦合边界,允许河道与地表进行水流交换。一、二维模型的计算步长均设置为0.5 s。

在研究区域上游、中游和下游分别设置3个雨量站以考虑降雨空间分布不均,选择2场实测降雨进行模型率定和验证。第1场降雨(简称降雨1)历时较短,开始于2017年7月15日18:24,结束于18:44,最大降雨强度为84 mm/h；第2场降雨(简称降雨2)发生于2017年9月5日13:00,持续时间为1 h,降雨强度峰值为318 mm/h,远超过第1场降雨峰值。根据多年内涝积水情况,在研究区域上安装了12台电子水尺以监测内涝淹没情况,在管道中设置了2台液位计以监测管道检查井水深,具体位置详见图5。

图5 研究区域示意及内涝淹没分布Fig.5 Study area and distribution of flood inundation

内涝点监控影像、电视新闻和现场监测信息等多种数据源信息被用来作为模型率定与验证数据的获取途径。根据暴雨发生当日的新闻报道,第1场降雨仅造成了极小范围的内涝淹没,而第2场降雨则造成了大面积的内涝淹没。监控影像提供了降雨过程中被监控位置的内涝信息,包括是否产生内涝和内涝淹没的程度。为了充分利用获取到的数据并考虑2个场次降雨内涝的特征,利用降雨1进行模型参数率定,采用检查井水深和内涝淹没位置作为模型参数率定的标准；利用降雨2进行模型参数验证,采用内涝淹没位置和淹没深度作为模型验证的标准。

4.3 模型验证结果与分析

在降雨1中,2个位于主干管道上检查井的水深变化过程被选取作为模型率定的研究对象。采用纳什效率系数(ENS)评价模型率定结果,一般来说,在下垫面特征较为复杂的城市区域,当ENS>0.6时,模型被认为具有良好精度。2个检查井水深变化过程的ENS值分别为0.914和0.602(图6),说明模型模拟精度良好。

图6 检查井水深过程线Fig.6 Water depth hydrograph in the manholes

同时,内涝淹没位置也是模型率定的一个重要标准,表2为模型模拟得到的内涝淹没位置与实地记录的淹没位置对比分析结果。从表2可以看出,在降雨1中,除了站点2以外,其余站点均未发生内涝,且实际发生情况与模拟情况完全吻合；在降雨2中,所有站点均发生了内涝,且除了站点4以外,实际发生情况与模拟情况完全吻合。以上分析表明模型模拟结果与实地监测结果一致,模型参数具有一定的准确性。

表2 内涝淹没位置模拟结果

除此之外,图5绘制了在降雨2情景下研究区域的最大淹没深度。从图中可以发现,模拟结果与新闻报道所反映的情况相符,流域中出现了严重的内涝淹没情况。

表3给出了模拟水深与实测水深的差异,从中可以看出,除S3、S4和S6之外,其余淹没位置水深的相对误差都小于30%。S4在模拟结果中被认为是不发生内涝淹没的位置,这与实际监测结果相反,所以相对误差较大,为100%。除了这些有监测记录的淹没位置之外,模拟结果还存在一些没有监测记录但也发生了内涝淹没的位置,与新闻报道所描写的大范围内涝淹没情况相一致。总的来说,所建立的城市洪涝模型能够反映出研究区域的现状降雨内涝积水情况,并具有较好的可靠性和精度。