浙江大学学报(工学版), 2020, 54(3): 614-622 doi: 10.3785/j.issn.1008-973X.2020.03.023

水利工程

基于孔隙率和局部时间步长的城市洪水模拟

李薇,, 邹吉玉, 胡鹏,

Urban flood simulation based on porosity and local time step

LI Wei,, ZOU Ji-yu, HU Peng,

通讯作者: 胡鹏,男,副教授. orcid.org/0000-0001-9214-1318. E-mail: pengphu@zju.edu.cn

收稿日期: 2019-04-11  

Received: 2019-04-11  

作者简介 About authors

李薇(1984—),女,副教授,从事河流海岸动力学、水沙模拟研究.orcid.org/0000-0002-6016-052X.E-mail:lw05@zju.edu.cn , E-mail:lw05@zju.edu.cn

摘要

为提高城市洪水模拟的计算效率,在有限体积法框架下,基于能自动捕捉激波、间断的高精度Harten-Latex-van Leer-Contact (HLLC)近似黎曼算子,并结合局部时间步长技术,建立满足静水平衡特性的各向异性孔隙率浅水模型. 经典城市洪水模拟结果表明,所建立的模型能够精确模拟洪水传播过程的复杂流动现象,并能显著提升计算效率:孔隙率方法降低了建筑物周围网格加密的要求,能使计算效率提升一个数量级;局部时间步长技术让每个网格采用尽可能大的时间步长,减少了循环次数,可进一步提升计算效率约2.0~3.0倍.

关键词: 各向异性孔隙率 ; 浅水方程 ; 城市洪水模拟 ; 局部时间步长 ; 有限体积法

Abstract

A well-balanced shallow water model of anisotropic porosity was developed under the framework of finite volume method, in order to increase the computational efficiency of urban flood simulation. The high-resolution Harten-Latex-van Leer-Contact (HLLC) approximate Riemann solver which could automatically capture shockwaves and discontinuities was used for flux computation, together with the local time step technique for time updating. The application to classic idealized urban floods shows that the represent model can accurately reproduce the complex hydrodynamics of urban floods and increase the computational efficiency markedly: the anisotropic porosity method reduces the requirement of local grid refining around obstructions and increases the computational efficiency by an order of magnitude; the local time step technique allows each grid to use a relatively large time step and reduces the temporal iteration, and further saves the computation cost by about two to three times.

Keywords: anisotropic porosity ; shallow-water equation ; urban flood simulation ; local time step ; finite volume method

PDF (2626KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

李薇, 邹吉玉, 胡鹏. 基于孔隙率和局部时间步长的城市洪水模拟. 浙江大学学报(工学版)[J], 2020, 54(3): 614-622 doi:10.3785/j.issn.1008-973X.2020.03.023

LI Wei, ZOU Ji-yu, HU Peng. Urban flood simulation based on porosity and local time step. Journal of Zhejiang University(Engineering Science)[J], 2020, 54(3): 614-622 doi:10.3785/j.issn.1008-973X.2020.03.023

城市洪涝灾害已成为影响城市安全和经济发展的重大问题之一,建立准确、高效的城市洪水预报数学模型意义重大. 为准确反映建筑群及其间的道路特征,传统洪水模型在建筑物周围局部加密网格,网格规模大、网格尺度小,使得模型计算效率很低. 为缩短计算时间,常将浅水方程简化为动力波模型[1]或扩散波模型[2],但其无法捕捉洪水传播的复杂水动力过程.

孔隙率浅水方程通过网格孔隙率反映建筑物和街道的阻水、过水效应,避免了局部加密网格的需要,在城市洪水模拟中受到关注[3-8]. 特别地,能自动捕捉激波和间断、灵活处理干湿变换的高精度数值格式使得孔隙率浅水模型具有更加广阔的应用前景[9-11]. 目前孔隙率浅水模型可分为3类. 第一类是各向同性孔隙率浅水方程[7, 12-13],此类方程假设流动在各个方向性质相同,对于每个网格单元建筑物采用单一的体积孔隙率计算. 第二类是各向异性孔隙率浅水方程[6, 14-15],该类方程在每个网格单元中同时考虑体积孔隙率和边界孔隙率. 体积孔隙率与各向同性孔隙率的定义类似,后者为网格单元边界可流动区域边长与单元边长的比值. Kim等[16]研究表明,各向异性孔隙率方程比各向同性孔隙率方程能更准确反映实际复杂城市街道间的非对称流动问题. 此外,还有部分比较有代表性的其他概念的孔隙率浅水方程. 比如,Guinot[4]基于运动孔隙率和静止孔隙率的概念建立了多重孔隙率浅水方程;Guinot等[17]对单元内和其边界的水流变量、孔隙率分别定义和积分,得到新形式的双积分孔隙率浅水方程. 由于孔隙率参数的引入,孔隙率浅水方程的动量方程需要满足通量、底坡项和建筑壁面静水压力项三者平衡,这使得原本针对传统浅水方程的平衡有限体积法(如表面梯度法[18])难以直接用于孔隙率浅水方程. 针对各向异性孔隙率浅水方程的平衡处理有一些初步进展,Sanders等[14]的研究忽略底坡项,近似认为建筑壁面静水压力项与通量相平衡;Guinot等[17]用各向同性浅水方程中的孔隙率梯度项直接替换各向异性孔隙率方程中的建筑壁面静水压力项,且求解时用单元的边界与体积孔隙率梯度代替各向同性时相邻单元体积孔隙率梯度.

除采用孔隙率浅水方程外,如何在保证计算稳定的同时尽可能采用较大的时间步长、减少变量更新次数,也是解决城市洪水模型计算效率的关键. 在浅水方程显格式求解时,为保证计算稳定需要用CFL(Courant-Friendrichs-Lewy)条件对时间步长加以限制. 当网格或流态分布不均时,不同网格间的时间步长会出现很大差异. 当采用传统的全局时间步长时,模型会将最小的时间步长作为整个计算区域的时间步长,从而降低整体的计算效率. 局部时间步长概念[19-22]可使每个计算网格根据各自的稳定条件独立更新,大大减少计算区域内的变量更新次数,显著提高非均匀网格和流态情况下的计算效率.

目前,孔隙率方法和局部时间步长技术在浅水数值模拟中独立发展. 本文将这2种高效的数值模拟技术相结合,旨在建立准确高效的浅水计算模型,并将其应用于城市洪水模拟. 具体而言,在有限体积法框架下,基于能自动捕捉激波、间断的高精度HLLC(Harten-Latex-van Leer-Contact)通量计算数值格式,并结合局部时间步长技术,建立满足静水平衡特性的各向异性孔隙率浅水模型. 本文将依次介绍模型的控制方程和数值算法,模型静水平衡特性和城市洪水计算能力和效率的验证,以及网格尺寸、划分方式和局部时间步长技术对孔隙率浅水模型计算的影响.

1. 控制方程

孔隙率浅水方程可写成如下积分形式[15]

$ \frac{\partial }{\partial t}\int_{\varOmega }{i{ U}{\rm d}\varOmega }+\oint_{\partial \varOmega }{i{ E} { n}{\rm d}r}{=}\int_{\varOmega }{i{ s} {\rm d}\varOmega }-\oint_{\partial {{\varOmega }^{{*}}}}{{{{ s}}^{{*}}}{\rm d}{{r}^{{*}}}}. $

$\begin{split} &{ U} = \left[ {\begin{array}{*{20}{c}} h \\ {hu} \\ {hv} \end{array}} \right], \;{ E} { n} = \left[ {\begin{array}{*{20}{c}} {h\left( {u{n_x} + v{n_y}} \right)}\\ {hu\left( {u{n_x} + v{n_y}} \right) + g{h^2}{n_x}/2}\\ {g{h^2}{n_y}/2 + hv\left( {u{n_x} + v{n_y}} \right)} \end{array}} \right],\\ & { s} = \left[ {\begin{array}{*{20}{c}} 0 \\ {gh\left( {{s_{{\rm b},x}} - {s_{{\rm f},x}}} \right)} \\ {gh\left( {{s_{{\rm b},y}} - {s_{{\rm f},y}}} \right)} \end{array}} \right]. \end{split} $

式中: $t$为时间, $i$为表征除建筑物以外孔隙的连续相位函数(详见式(5)), $\varOmega $为单元(如图1所示),U为守恒变量向量, $\int_\varOmega {f{\rm{d }}\varOmega} $表示在单元 $\varOmega $对变量 $f$积分( $f$为任意函数或变量组合),E为通量向量, ${{n}} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {{n_x},}\;{{n_y}} \end{array}}\!\!\!\! \right]$为单元边界的外法向向量(其中 ${n_x}$${n_y}$分别为向量n$x$$y$方向上的组成) $\oint_{\partial \varOmega } {f{\rm{}}{\rm d}r} $为沿着单元边界 $\partial \varOmega $对变量 $f$积分,r为沿逆时针方向环绕单元边界的距离,s为包括摩阻项和坡度项的源汇项向量, ${\varOmega ^*}$代表建筑物, ${{ s}^*}$为建筑物阻力项(详见式(7)), $\oint_{\partial {\varOmega ^*}} {f{\rm{ }}{\rm d}{r^*}} $为沿着建筑物边界 $\partial {\varOmega ^*}$对变量 $f$积分,r*为沿逆时针方向环绕建筑物边界的距离, $h$为水深, $u$$v$$x$$y$方向上的流速, $g$为重力加速度, ${s_{{\rm b},x}}$${s_{{\rm b},y}}$为底部坡度项, ${s_{{\rm f},x}}$${s_{{\rm f},y}}$为阻力坡度,其表达式分别为

图 1

图 1   三角网格下孔隙率的概念模型图

Fig.1   Conceptual model diagram of porosity in triangular cell


$ {s_{{\rm b},x}} = - \partial {Z_{\rm b}}/\partial x,\;\;\; {s_{{\rm b},y}} = - \partial {Z_{\rm b}}/\partial y, $

$ {s_{{\rm f},x}} = {h^{{\rm{} - }{{\rm{4}}/{\rm{3}}}}}{n^2}uV, \;\;\;{s_{{\rm f},y}} = {h^{{\rm{} - }{{\rm{4}}/{\rm{3}}}}}{n^2}vV . $

其中, $n$为曼宁系数, ${Z_{\rm b}}$为底部高程, $V =\sqrt {{u^2} + {v^2}} $.

连续相位函数的值取决于水面高程与地面高程之间的关系,计算式如下:

$ i\left( {x,y} \right) = \left\{ {\begin{array}{*{20}{l}} {1,}&{\eta \left( {x,y} \right) > {Z_{\rm b}}\left( {x,y} \right);}\\ {0,}&{{\text{其他}}.} \end{array}} \right. $

式中: $\eta $为水面高程.

为了更好地捕捉水流的流动方向的差异(即各向异性),将表征孔隙的连续相位函数分解为单元孔隙率 $\phi $(即单元内可流动区域面积与单元总面积的比值)和边界孔隙率 $\psi $(即边界上可流动区域长度与边界长度的比值)[14]

$ \phi = \frac{{\mathop \smallint \nolimits_\varOmega i\left( {\eta - {Z_{\rm b}}} \right){\rm d}\varOmega }}{{\mathop \smallint \nolimits_\varOmega \left( {\eta - {Z_0}} \right){\rm d}\varOmega }},\;\psi = \frac{{\mathop \oint \nolimits_{\partial \varOmega } i\left( {\eta - {Z_{\rm b}}} \right){\rm d}r}}{{{\rm{}}\mathop \oint \nolimits_{\partial \varOmega } \left( {\eta - {Z_0}} \right){\rm d}r}}. $

式中: ${Z_0}$为单元内地形高程最小值,Sanders等[14]将建筑物阻力 ${{ s}^*}$分成2个部分:建筑物面的静水压力,以及建筑在动水下的额外作用力,即

$\begin{split} \oint_{\partial {\varOmega ^{\rm{*}}}} {{{ s}^{\rm{*}}}{\rm d}{r^{\rm{*}}}} {\rm{ = }}& \oint_{\partial {\varOmega ^{\rm{*}}}} {\frac{1}{2}g{{({{\left. h \right|}_{{{\rm{\eta}}_{\rm 0}}}})}^2}{ m}{\rm d}{r^{\rm{*}}} + } \\ & \int_\varOmega {iC_{\rm D}^{\rm b}{ u}V {\rm d}\varOmega } . \end{split} $

式中: $C_{\rm D}^{\rm b}$为建筑拖曳力系数, ${ u} = \left[ 0,\;u,\;v \right]$$h{|_{{\eta _0}}}$为假设水面静止时的水深, ${ m}$为建筑物壁面指向建筑内部的法向向量(见图1). 需要注意的是,式(7)把Nepf [23]动水下的建筑物额外作用力平均到了整个单元. 将式(6)、(7)代入式(1)后,可得各向异性的孔隙浅水方程:

$ \begin{split} &\frac{\partial }{\partial t}\int_{\varOmega }{\phi { U}{\rm d}\varOmega }+\oint_{\partial \varOmega }{\psi { E} { n}{\rm d}r}{=} \\ & \;\;\;\;\;\;\;\;\;\;\;\;\int_{\varOmega }{\phi \left( { s}-C_{{\rm D}}^{{\rm b}}{ u}V \right){\rm d}\varOmega }-\oint_{\partial {{\varOmega }^{*}}}{\frac{1}{2}g{{({{\left. h \right|}_{{{{\rm{\eta}}}_{{\rm 0}}}}})}^{2}}{ m}{\rm d}{{r}^{*}}} {.} \end{split} $

需要指出的是,令孔隙率浅水方程中所有孔隙率等于1即可得到传统浅水方程,同时额外的源项 ${{ s}^*}$(式(7))也无须再计算,这是因为没有建筑, $C_{\rm D}^{\rm b}$$\oint_{\partial {\varOmega ^*}} {0.5g{{({{\left. h \right|}_{{{\rm{\eta}}_{\rm 0}}}})}^2}{ m}{\rm d}{r^*}} $都为0[14]. 本文中传统浅水方程的求解方法同孔隙率浅水方程.

2. 数值方法

2.1. 离散格式

在有限体积法框架下,参照Hu等[20,24]的研究,将式(8)在非结构三角形网格上离散:

$\!\!\begin{split} { U}_i^{n + 1} = & { U}_i^n + \Delta t\left( {{ s} - C_{\rm D}^{\rm b}{ u}V} \right) + \\ & \!\!\frac{{\Delta t}}{{{\phi _i}{A_i}}}\left( { - \mathop \sum \limits_{j{\rm = }1}^3 {\psi _j}{{ E}_j} {{ n}_j}\Delta {r_j} - \oint_{\partial {\varOmega ^*}} {\frac{1}{2}g{{({{\left. h \right|}_{{\eta _0}}})}^2}{ m}{\rm d}{r^*}} } \right). \\ \end{split} $

式中: $\Delta t$为时间步长, ${A_i}$为单元面积, $\Delta {r_j}$为单元边界长度.

需要指出的是,式(9)中右端建筑物面静水压力的离散必须考虑静水平衡特性,故未直接写出其离散式. 对于该项的处理,考虑到静水时建筑物壁面上的静水压力与单元边界上的静水压力以及底面上的作用力平衡[14],在静水条件下简化式(8)得到如下表达式:

$\begin{split} & \oint_{\partial \varOmega } {\frac{1}{2}g\psi {{(h{|_{{{\rm{\eta}}_{\rm 0}}}})}^2}{ n}{\rm d}r} {\rm{}} = - \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\oint_{\partial {\varOmega ^*}} {\frac{1}{2}g{{({{\left. h \right|}_{{{\rm{\eta}}_{\rm 0}}}})}^2}{ m}{\rm d}{r^*}} + \int_\varOmega {\phi gh{{\rm |}_{{{\rm{\eta}}_{\rm 0}}}}{{ s}_{\rm b}}{\rm d}\varOmega }. \\ \end{split}$

式(10)中的边界静水压力项可离散为

$\mathop \oint \limits_{\partial \varOmega } \frac{1}{2}g\psi {({\left. h \right|_{{{\rm{\eta}}_{\rm 0}}}})^2}{ n}{\rm d}r = \sum\nolimits_{j = 1}^3 {\frac{1}{2}\left( {g{\psi _j}{h_j}^2{{ n}_j}\Delta {r_j}} \right)} .$

式中: ${h_j}$为界面中心处水深,为避免负水深,按照文献[25]的方法修正,

${h_j} = \max\; \left( {0,{\eta _{{\rm 0}i}} - {Z_{{\rm b}j}}} \right).$

式中: ${\eta _{0i}}$为静止时单元水面高程,取值为单元中心点高程和水深之和, ${Z_{{\rm b}j}}$为界面中心处高程.

式(12)在计算相邻单元之间的界面水深时,通过 ${\rm{max}}$函数保证了界面水深≥0. 此处界面水深是计算边界静水压力需要用到的临时变量:通过式(12)对边界静水压力的修正不影响动量或质量守恒. 这能确保单元之间的物质或动量输运不过大或过小,从而避免出现负水深. 在传统浅水方程中,底坡源项可以转化为边界量计算,在孔隙率浅水方程中参考文献[26]的Slope flux method方法将式(10)右边的底坡源项转化到单元的3条边进行离散:

$\int_\varOmega {\phi g{{\left. h \right|}_{{{\rm{\eta}}_{\rm 0}}}}{{ s}_{\rm{b}}}{\rm d}\varOmega } = - \mathop \sum \limits_{j = 1}^3 \frac{1}{2}g\phi \left( {{h_i} + {h_j}} \right)\left( {{Z_{{\rm b}j}} - {Z_{{\rm b}i}}} \right){{{n}}_j}\Delta {r_j}.$

式中: ${h_i}$为单元水深, ${Z_{{\rm b}i}}$为单元中心点高程. 由传统浅水方程的静水平衡可知在平衡条件下满足:

$\mathop \sum \limits_{j = 1}^3 \frac{1}{2}g\left( {{h_i} + {h_j}} \right)\left( {{Z_{{\rm b}j}} - {Z_{{\rm b}i}}} \right){{{n}}_j}\Delta {r_j} + \mathop \sum \limits_{j = 1}^3 \frac{1}{2}\left( {g{h_j}^2{{ n}_j}\Delta {r_j}} \right) = 0.$

当把式(10)中的项转换到边界计算时,在平衡条件下应同时满足式(14),将式(11)~(13)代入式(10),并通过式(10)~(14)得到建筑壁面的静水压力项离散形式如下:

$\begin{split} \oint_{\partial {\varOmega ^*}} &{\frac{1}{2}g{{({{\left. h \right|}_{{{\rm{\eta}}_{\rm 0}}}})}^2}{ m}{\rm d}{r^*}} \!\!=\!\!\mathop \sum\nolimits_{j = 1}^3 -\frac{1}{2}\left[ {g\left( {{\psi _j} - 1} \right){h_j}^2{{{n}}_j}\Delta {r_j}} \right] - \\ &\;\;\;\;\;\;\mathop \sum\nolimits_{j = 1}^3 \frac{1}{2}g\left( {\phi - 1} \right)\left( {{h_i} + {h_j}} \right)\left( {{Z_{{\rm b}j}} - {Z_{{\rm b}i}}} \right){{{n}}_j}\Delta {r_j}. \end{split} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$

式(9)中其他项的离散格式参考文献[20],此处不再赘述.

2.2. 局部时间步长

局部时间步长(local time step,LTS)使每个单元以其特定的时间步长更新,与全局时间步长使用同一个最小的时间步长更新相比,可大幅提高计算效率. 目前还没有在孔隙率浅水方程模拟中直接使用局部时间步长的先例,因此各向异性孔隙率结合局部时间步长来提高计算效率值得研究. 孔隙率浅水方程的时间步长计算公式采用Sanders等[14]提出的方法,同时为了避免 ${\phi _i}/\max \; $ $({\psi _j})$出现极小值,将公式修正如下:

$ \Delta {t_i} \leqslant {\rm{CFL}}\left( {\frac{{{\phi _i}{\Omega _i}}}{{{\rm{mean}}\;\left( {{\psi _j}} \right){\rm{max}}\;{{\left( {\partial {\Omega _j}} \right)}_{j { = 1,2,3}}}}}\;\frac{1}{{{c_j}}}} \right). $

式中: $\Delta {t_i}$为时间步长; ${c_j} = \sqrt {gh_j} + {v_j}$vj为单元边界流速;CFL为Courant数,为了避免失稳取一个较小值[15],本文取为0.45.

为了计算的稳定,在干湿边界处须作一些修正. 将干单元的 $\Delta {t_i}$定为极大值,当 ${\phi _i} < 0.001$时,令 ${\phi _i}$=0,将所有孔隙率等于0的单元作干单元处理,如果在干湿边界处水深小于临界值 ${10^{ - 6}}$,则令该单元的时间步长等于计算区域的最大时间步长. 计算更新步骤引用Hu等[20]提出的方法,首先选取最小时间步长 $\Delta {t_0} = \min\; (\Delta {t_i})$作为整个计算区域的基础更新时间步长,再定义各单元的潜在局部时间步长 $\Delta {t_{{ L} - i}}$为最小时间步长的幂函数,即

$\Delta {t_{{ L} - i}} = {2^{{m_{{{L} - }i}}}}\Delta {t_0}.$

式中: $\Delta {t_{{ L} - i}}$为每个单元的潜在局部时间步长(L), ${m_{{{L} - }i}}$为单元的潜在更新层级,

${m_{{ L} - i}} = \min \;\left[ {{\rm int}\; \left( {\frac{{\log\; \left( {\Delta {t_i}/\Delta {t_0}} \right)}}{{\log {\rm{2}}}}} \right),\;{m_{{\rm max}}}} \right].$

式中: ${\rm int} $为取整函数, ${m_{{\rm max}}}$为用户设定的最大更新层级( ${m_{{\rm{max}}}}$=0时代表传统求解方法).

通过相邻单元的潜在局部时间步长最终得到各单元的实际局部更新时间步长:

$\Delta {t_{{{L}^{'}} - i}} = \min\; \left( {{m_{{ L} - i}},{m_{j1}},{m_{j2}},{m_{j3}}} \right)\Delta {t_0}.$

式中: $\Delta {t_{{{L}^{'}} - i}}$为第 $i$个单元的实际局部更新时间步长( ${{L}^{'}} $), ${m_{j1}},{m_{j2}},{m_{j3}}$为三邻边对应单元的潜在更新层级. 对整个计算区域而言,完成一轮更新的时间步长为 ${2^{{m_{\rm{max}}}}}\Delta {t_0}$,其间每次时间递进 $\Delta {t_0}$,即递进 $l$次时间前进 $l\Delta {t_0}$;任意单元 $i$仅在 $l\Delta {t_0}$$\Delta {t_{{ {{L}^{'}} - i}}}$的整数倍时进行更新,其更新时间步长为 $\Delta {t_{{ {{L}^{'}} - i}}}$. 当整个计算区域各单元均更新了 ${2^{{m_{{\rm{max}}}}}}\Delta {t_0}$后再进行下一轮整体更新.

3. 模型计算

3.1. 参数设置与网格划分

本文选取比利时UCL的理想城市溃坝洪水实验[27]来验证所建立的模型. 如图2(根据文献[28]修改)所示,实验水槽长36 m、宽3.6 m,底部水平、横断面为梯形. 闸门上游为水库,初始水深为0.4 m,下游为理想化城市区域,初始水深为0.011 m. 理想城市与溃坝水流方向一致,建筑物尺寸为0.3 m ×0.3 m、街道宽度为0.1 m. 通过瞬间开启闸门形成向下游理想城市传播的溃坝洪水. 除水槽下游考虑为自由出流外,其他为固壁边界. 糙率取为0.01,计算时长为15.5 s,使用的处理器为intel(R)Corel(TM)i7-7700CPU@3.60 GHz.

图 2

图 2   理想城市洪水实验布置图

Fig.2   Experimental set-up and dimensions for idealized urban floods


不同于各向同性孔隙率基于计算域的统计数据定义孔隙率,各向异性孔隙率取决于网格和建筑轮廓相交的部分,因此各向同性孔隙率浅水方程对网格划分不敏感,而各向异性孔隙率浅水方程则对网格较为敏感[17]. Sanders等[14]认为不同的网格划分方式会导致不同的边界孔隙率,因此本文设置不同的网格划分方式以探究其对结果的影响.

为了对比传统浅水方程和各向异性孔隙率浅水方程的计算精度和效率,以及网格划分对各向异性孔隙率方程计算结果的影响,本研究共设计5种计算方案. 方案1(记作C_C,见图3(a))采用传统浅水方程,通过收敛性分析确定满足计算精度的最少网格单元为13 627个,其中街道网格单元边长为0.03 m、上游水库网格边长为0.3 m、其余边壁网格边长为0.15 m. 方案2~5均采用各向异性孔隙率方程,其中方案2(记作P_RN,见图3(b),图右侧色阶标尺表示体积孔隙率大小)各边界尺寸与方案1相同、建筑区域尺寸为0.05 m,网格单元为13 500个,与方案1数量相当,旨在对比网格数量相近时2种方程的计算效果. 方案3(记作P_MN)在方案2基础上粗化建筑物区域网格为0.08 m,共10 259个网格单元,反映孔隙率方程计算精度与精细网格相同时所需的最少网格量. 方案4(记作P_CG,见图3(c))采用gap-conforming网格划分方式[14]将网格点有序地放在建筑物内部,网格边长为0.4 m、共1 545个网格. 方案5(记作P_CN,见图3(d))网格数(1 543个)、网格尺寸与方案4相当但为随机网格,随机网格不指定网格顶点,边界尺寸与P_CG一致,内部网格由SMS软件随机生成,旨在比较网格布置方式的影响.

图 3

图 3   各计算方案在建筑物区域的网格划分图

Fig.3   Diagram of mesh types for each test calculation case in building area


3.2. 平衡验证

为检验模型平衡(well-balance)特性,在x=8.55 m处,地形底坡由平坡变为坡度为1%的斜坡,且四周设为固壁边界,在P_CG网格下开展2个平衡算例. 算例1为溃坝算例,图4展示了测点A1(9.45,2)、A2(11.45,2)、A3(12.45,2)水深随时间的变化. 可以看出,模型能准确反映出溃坝水流由运动到静止的状态. 算例2为初始静水算例,图5显示了模型能一直保持初始静水状态(图右侧色阶标尺表示高程),满足平衡特性.

图 4

图 4   测点水深随时间的变化(平衡算例1)

Fig.4   Water depth evolution over time(balance case 1)


图 5

图 5   不同时刻的水位和地形三维视图(平衡算例2)

Fig.5   3D view of water level and bed topography at different times(balance case 2)


3.3. 计算结果

图6展示了传统浅水方程(方案1:C_C)和各向异性孔隙率方程(方案2:P_RN)的水深计算云图. 可以看出,当网格数量相当时,传统浅水方程和各向异性孔隙率方程的水深等值线基本对称、具有同等精度. 在建筑群前方,水深受建筑群阻挡作用而明显增加;进入建筑间街道后,水深沿程逐渐减小;在建筑群出口处水深达最小值,而后在建筑群后方由于水流的汇聚作用而增加. 具体来看,对于不同代表位置的4个测点1(12.55,2.8)、18(13.3,2.4)、44(14.45,1.6)、55(12.3,2.0)水深,各向异性孔隙率浅水模型与传统浅水模型计算结果吻合良好、均能较好反映实测水深变化过程(见图7). 对于沿街道B-B’(y=2 m),这2类模型的计算水深和流速与实测也吻合较好(见图89);特别是,各向异性孔隙率方程与传统浅水方程一样能精细模拟出建筑物和街道间的水深高低起伏,甚至在流速计算上前者优于后者.

图 6

图 6   方案C_C和P_RN在10 s时的水深等值线对比图

Fig.6   Comparison of water depth contours between case C_C and P_RN at experimental time of 10 s


图 7

图 7   理想化城市洪水实验不同测点处水深随时间的变化

Fig.7   Water depth against time at different gauges for idealized urban flood experiment


图 8

图 8   理想化城市洪水实验不同时刻沿街道y=2 m的水深

Fig.8   Water depth profiles along longitudinal street of y=2 m at different times for idealized urban flood experiment


图 9

图 9   理想化城市洪水实验在不同时刻沿街道y=2 m的流速

Fig.9   Flow velocity profiles along longitudinal street of y=2 m at different times for idealized urban flood experiment


孔隙率浅水方程的计算精度与网格大小密切相关,网格越小计算精度越高、越能细致反映建筑街道间的水流变化情况,但小网格同时限制了计算速度、使洪水预报的实时性降低,因此实际应用中需要在计算精度和效率之间寻找一个平衡.

图7~9的计算结果可看出,若需预测建筑街道间水深和流速的精细波动(见图89),孔隙率网格也不能过大,如方案3(P_MN),此时孔隙率浅水方程的计算时间比传统浅水方程节省23%左右(见图10). 对比方案1(C_C)和方案2(P_RN)可看出,当网格数量接近时,孔隙率浅水方程由于增加了额外的源项计算,可能需要更多的计算时间. 当把网格数量大幅减少时,如方案4(P_CG)的网格数为方案1的1/10,孔隙率浅水方程计算时间仅为传统浅水方程的1.7%;并且孔隙率浅水方程的计算水深和流速变化范围、趋势也与实测吻合良好,仅建筑街道间的精细波动未能反映. 对比方案4和5,粗网格下不同网格布置方式对孔隙率浅水方程计算的水深影响不大,但对流速影响较明显,P_CG方式对流速的计算更准确,这与Sanders等[14]的结论一致.

图 10

图 10   各方案在不同局部时间步长(LTS)最大更新层级情况下的计算时间加速倍数

Fig.10   Acceleration rate of computation time for each approach at different maximal updating levels of local time step (LTS)


为了定量分析模型的误差,将实验数据和实测结果进行对比,误差计算如下:

${{\delta}} = {\left[ {\frac{1}{N}\mathop \sum \limits_{n = 1}^N {{\left( {{f_{{\rm{sim}}}} - {f_{{\rm{exp}}}}} \right)}^2}} \right]^{{\rm{1/2}}}}.$

其中: $\delta $为误差,N为实测数据个数, ${f_{{\rm{sim}}}}$为模拟值, ${f_{{\rm{exp}}}}$为实测值. 计算结果如表1所示,其中 ${\delta_{{\rm{hg}}}}$为测点水深误差, ${\delta_{{\rm{hs}}}}$${\delta_{\rm{V}}}$分别为沿街道B-B’的水深和流速在4、5、6和10 s时刻的误差. 从定量分析看,水深计算的误差较小且分布于0.013 ~0.042 m;虽然流速计算误差稍大,在0.121~0.466 m/s,但仍与实测吻合较好,能合理反映出流速的沿程变化.

表 1   不同方案的计算值与实测值的误差比较

Tab.1  Comparison of errors between calculated and measured values by different solutions

方案 ${ { {\delta} }_{ {\rm{h} }g} }$/m ${ { {\delta} }_{ {\rm{hs} } } }$/m ${ { {\delta} }_{\rm{V} } }$/(m·s–1
测点1 测点18 测点44 测点55 t=4 s t=5 s t=6 s t=10 s t=4 s t=5 s t=6 s t=10 s
C_C 0.042 0.015 0.022 0.032 0.030 0.023 0.022 0.013 0.466 0.253 0.196 0.135
P_CG 0.032 0.019 0.037 0.027 0.034 0.036 0.038 0.024 0.325 0.292 0.317 0.194
P_CN 0.030 0.017 0.028 0.026 0.033 0.041 0.040 0.023 0.297 0.384 0.288 0.263
P_MN 0.038 0.013 0.028 0.023 0.023 0.017 0.021 0.014 0.309 0.212 0.144 0.121
P_RN 0.034 0.013 0.031 0.024 0.023 0.017 0.020 0.014 0.330 0.205 0.174 0.159

新窗口打开| 下载CSV


图10显示了各方案在不同的局部时间步长最大更新层级情况下的计算加速效果,这里定义计算时间加速倍数a为方案C_C在 ${m_{{\rm{max}}}}$=0时的计算时间与其他方案计算时间的比值. 可以看出,孔隙率浅水方程与局部时间步长技术的结合能大大提高模型的计算效率. 仅考虑局部时间步长时会有1~5倍的计算加速,如传统浅水方程使用局部时间步长( ${m_{{\rm{max}}}}$=3)能使计算提速近4倍. 在 ${m_{{\rm{max}}}}$=0时C_C和P_CG的计算时间比值为58,在 ${m_{{\rm{max}}}}$=3时比值为38,说明孔隙率浅水模型具有显著加速效果;同时,不同 ${m_{{\rm{max}}}}$下的计算时间比值差异也表明局部时间步长的时间加速效果在各网格尺寸差异大的情况(方案1:C_C)下更明显. 在保证同等精度、反映水深流速的精细波动条件下,与传统浅水方程相比,孔隙率浅水方程和局部时间步长结合(P_MN,mmax=3)使计算效率提高约5倍. 在保证水深流速趋势和范围正确、损失局部细节的精度(如建筑街道间的波动趋于平滑)但又满足计算效率条件下,局部时间步长与孔隙率浅水方程的结合(P_CG, ${m_{{\rm{max}}}}$=3)能使计算效率提高近140倍.

精细网格(如P_RN)的计算结果更加精确,但计算效率没有优势,粗网格(P_CN和P_CG)的计算效率更高,计算精度也令人满意. 值得注意的是,在随机网格中可能出现极小的孔隙率,需要借助工具(如ArcGIS)计算孔隙率以避免数据误差,使用样本点计算孔隙率[17]可能会降低计算稳定性. 总体而言,P_CG是最优选项,但在建筑密集情况下指定网格须作较为复杂的前处理.

4. 结 论

(1)在网格数量与传统浅水方程在同一数量级的情况下,各向异性孔隙率浅水方程可精细反映建筑街道间的水深流速沿程波动,与传统浅水方程具有同等精度;随着网格数量减少,孔隙率浅水方程精度下降,但在保证水深流速趋势和范围正确、损失局部细节精度(如建筑街道间的波动趋于平滑)的情况下,网格数量可减少到原来的1/10.

(2)在精细反映建筑街道间水深流速沿程波动的前提下,局部时间步长技术比各向异性孔隙率方法的计算效率更高,而两者的结合可进一步缩短计算时间;但在大范围城市洪水实时预测预报中,为了保证计算效率往往不得不牺牲计算精度,比如本文各向异性孔隙率浅水方程在平滑了建筑街道间的沿程波动后,计算速度比传统浅水模型(mmax=0)提高了50多倍,具有明显优势,结合局部时间步长,计算速度可进一步提高2.0~3.0倍.

需要指出的是,本文计算效率提高的认识仅基于上述算例,如果采用其他算例,模型的加速效果在定量上可能会有差异. 各向异性孔隙率浅水方程与局部时间步长技术的结合在城市洪水模拟和预报中具有广阔的应用前景,如果能结合GPU并行技术,计算效率还有更大提升空间[29].

参考文献

DE ALMEIDA G A M, BATES P D

Applicability of the local inertial approximation of the shallow water equations to flood modeling

[J]. Water Resources Research, 2013, 49 (8): 4833- 4844

DOI:10.1002/wrcr.20366      [本文引用: 1]

COSTABILE P, COSTANZO C, MACCHIONE F

Performances and limitations of the diffusive approximation of the 2-D shallow water equations for flood simulation in urban and rural areas

[J]. Applied Numerical Mathematics, 2017, 116: 141- 156

DOI:10.1016/j.apnum.2016.07.003      [本文引用: 1]

CHEN A S, EVANS B, DJORDJEVIC S, et al

A coarse-grid approach to representing building blockage effects in 2D urban flood modelling

[J]. Journal of Hydrology, 2012, 426 (Suppl. C): 1- 16

[本文引用: 1]

GUINOT V

Multiple porosity shallow water models for macroscopic modelling of urban floods

[J]. Advances in Water Resources, 2012, 37 (1): 40- 72

[本文引用: 1]

GUINOT V, DELENNE C

Macroscopic modelling of urban floods

[J]. La Houille Blanche, 2015, (6): 19- 25

VELICKOVIC M, ZECH Y, SOARES-FRAZAO S

Steady-flow experiments in urban areas and anisotropic porosity model

[J]. Journal of Hydraulic Research, 2017, 55 (1): 85- 100

DOI:10.1080/00221686.2016.1238013      [本文引用: 1]

贺治国, 翁浩轩, 高冠, 等

植被群影响下动床洪水数值模拟研究

[J]. 应用基础与工程科学学报, 2017, 25 (1): 17- 27

[本文引用: 1]

HE Zhi-guo, WENG Hao-xuan, GAO Guan, et al

Numerical simulation of dam-break flood over mobile and vegetated beds

[J]. Journal of Basic Science and Engineering, 2017, 25 (1): 17- 27

[本文引用: 1]

周浩澜, 陈洋波

城市化地面二维浅水模拟

[J]. 水科学进展, 2011, 22: 407- 412

[本文引用: 1]

ZHOU Hao-lan, CHEN Yang-bo

2D shallow water simulation for urbanized areas

[J]. Advances in Water Science, 2011, 22: 407- 412

[本文引用: 1]

COZZOLINO L, PEPE V, CIMORELLI L, et al

The solution of the dam-break problem in the Porous Shallow water Equations

[J]. Advances in Water Resources, 2018, 114: 83- 101

DOI:10.1016/j.advwatres.2018.01.026      [本文引用: 1]

FERRARI A, VACONDIO R, DAZZI S, et al

A 1D - 2D shallow water equations solver for discontinuous porosity field based on a generalized riemann problem

[J]. Advances in Water Resources, 2017, 107

HE Z, WU T, WENG H, et al

Numerical simulation of dam-break flow and bed change considering the vegetation effects

[J]. International Journal of Sediment Research, 2017, 32 (1): 105- 120

DOI:10.1016/j.ijsrc.2015.04.004      [本文引用: 1]

SOARES-FRAZAO S, LHOMME J, GUINOT V, et al

Two-dimensional shallow-water model with porosity for urban flood modelling

[J]. Journal of Hydraulic Research, 2008, 46 (1): 45- 64

DOI:10.1080/00221686.2008.9521842      [本文引用: 1]

FINAUD-GUYOT P, DELENNE C, LHOMME J, et al

An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity

[J]. International Journal for Numerical Methods in Fluids, 2010, 62 (12): 1299- 1331

[本文引用: 1]

SANDERS B F, SCHUBERT J E, GALLEGOS H A

Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling

[J]. Journal of Hydrology, 2008, 362 (1-2): 19- 38

DOI:10.1016/j.jhydrol.2008.08.009      [本文引用: 10]

OZGEN I, ZHAO J H, LIANG D F, et al

Urban flood modeling using shallow water equations with depth-dependent anisotropic porosity

[J]. Journal of Hydrology, 2016, 541: 1165- 1184

DOI:10.1016/j.jhydrol.2016.08.025      [本文引用: 3]

KIM B, SANDERS B F, FAMIGLIETTI J S, et al

Urban flood modeling with porous shallow-water equations: a case study of model errors in the presence of anisotropic porosity

[J]. Journal of Hydrology, 2015, 523: 680- 692

DOI:10.1016/j.jhydrol.2015.01.059      [本文引用: 1]

GUINOT V, SANDERS B F, SCHUBERT J E

Dual integral porosity shallow water model for urban flood modelling

[J]. Advances in Water Resources, 2017, 103: 16- 31

DOI:10.1016/j.advwatres.2017.02.009      [本文引用: 4]

ZHOU J G, CAUSON D M, MINGHAM C G, et al

The surface gradient method for the treatment of source terms in the shallow-water equations

[J]. Journal of Computational Physics, 2001, 168 (1): 1- 25

DOI:10.1006/jcph.2000.6670      [本文引用: 1]

SANDERS B F

Integration of a shallow water model with a local time step

[J]. Journal of Hydraulic Research, 2008, 46 (4): 466- 475

DOI:10.3826/jhr.2008.3243      [本文引用: 1]

HU P, LEI Y, HAN J, et al

Computationally efficient modeling of hydro-sediment-morphodynamic processes using a hybrid local time step/global maximum time step

[J]. Advances in Water Resources, 2019, 127: 26- 38

DOI:10.1016/j.advwatres.2019.03.006      [本文引用: 3]

胡鹏, 韩健健, 雷云龙

基于局部分级时间步长方法的水沙耦合数学模拟

[J]. 浙江大学学报: 工学版, 2019, 53 (4): 743- 752

HU Peng, HAN Jian-jian, LEI Yun-long

Coupled modeling of sediment-laden flows based on local-time-step approach

[J]. Journal of Zhejiang University: Engineering Science, 2019, 53 (4): 743- 752

HU P, HAN J, LI W, et al

Numerical investigation of a sandbar formation and evolution in a tide-dominated estuary using a hydro-morphodynamic model

[J]. Coastal Engineering Journal, 2018, 60 (4): 466- 483

DOI:10.1080/21664250.2018.1529263      [本文引用: 1]

NEPF H M

Drag, turbulence, and diffusion in flow through emergent vegetation

[J]. Water Resources Research, 1999, 35 (2): 479- 489

DOI:10.1029/1998WR900069      [本文引用: 1]

YUE Z, LIU H, LI Y, et al

A well-balanced and fully coupled noncapacity model for dam-break flooding

[J]. Mathematical Problems in Engineering, 2015, 1- 13

[本文引用: 1]

AUDUSSE E, BOUCHUT F, BRISTEAU M, et al

A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows

[J]. Journal on Scientific Computing, 2004, 25 (6): 2050- 2065

DOI:10.1137/S1064827503431090      [本文引用: 1]

HOU J M, SIMONS F, MAHGOUB M, et al

A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography

[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 257: 126- 149

DOI:10.1016/j.cma.2013.01.015      [本文引用: 1]

SOARES-FRAZAO S, ZECH Y

Dam-break flow through an idealised city

[J]. Journal of Hydraulic Research, 2008, 46 (5): 648- 658

DOI:10.3826/jhr.2008.3164      [本文引用: 1]

OZGEN I, LIANG D F, HINKELMANN R

Shallow water equations with depth-dependent anisotropic porosity for subgrid-scale topography

[J]. Applied Mathematical Modelling, 2016, 40 (17/18): 7447- 7473

[本文引用: 1]

VACONDIO R, DAL PALU A, MIGNOSA P

GPU-enhanced finite volume shallow water solver for fast flood simulations

[J]. Environmental Modelling and Software, 2014, 57: 60- 75

DOI:10.1016/j.envsoft.2014.02.003      [本文引用: 1]

/