基于孔隙率和局部时间步长的城市洪水模拟
Urban flood simulation based on porosity and local time step
通讯作者:
收稿日期: 2019-04-11
Received: 2019-04-11
作者简介 About authors
李薇(1984—),女,副教授,从事河流海岸动力学、水沙模拟研究.orcid.org/0000-0002-6016-052X.E-mail:
为提高城市洪水模拟的计算效率,在有限体积法框架下,基于能自动捕捉激波、间断的高精度Harten-Latex-van Leer-Contact (HLLC)近似黎曼算子,并结合局部时间步长技术,建立满足静水平衡特性的各向异性孔隙率浅水模型. 经典城市洪水模拟结果表明,所建立的模型能够精确模拟洪水传播过程的复杂流动现象,并能显著提升计算效率:孔隙率方法降低了建筑物周围网格加密的要求,能使计算效率提升一个数量级;局部时间步长技术让每个网格采用尽可能大的时间步长,减少了循环次数,可进一步提升计算效率约2.0~3.0倍.
关键词:
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:
本文引用格式
李薇, 邹吉玉, 胡鹏.
LI Wei, ZOU Ji-yu, HU Peng.
孔隙率浅水方程通过网格孔隙率反映建筑物和街道的阻水、过水效应,避免了局部加密网格的需要,在城市洪水模拟中受到关注[3-8]. 特别地,能自动捕捉激波和间断、灵活处理干湿变换的高精度数值格式使得孔隙率浅水模型具有更加广阔的应用前景[9-11]. 目前孔隙率浅水模型可分为3类. 第一类是各向同性孔隙率浅水方程[7, 12-13],此类方程假设流动在各个方向性质相同,对于每个网格单元建筑物采用单一的体积孔隙率计算. 第二类是各向异性孔隙率浅水方程[6, 14-15],该类方程在每个网格单元中同时考虑体积孔隙率和边界孔隙率. 体积孔隙率与各向同性孔隙率的定义类似,后者为网格单元边界可流动区域边长与单元边长的比值. Kim等[16]研究表明,各向异性孔隙率方程比各向同性孔隙率方程能更准确反映实际复杂城市街道间的非对称流动问题. 此外,还有部分比较有代表性的其他概念的孔隙率浅水方程. 比如,Guinot[4]基于运动孔隙率和静止孔隙率的概念建立了多重孔隙率浅水方程;Guinot等[17]对单元内和其边界的水流变量、孔隙率分别定义和积分,得到新形式的双积分孔隙率浅水方程. 由于孔隙率参数的引入,孔隙率浅水方程的动量方程需要满足通量、底坡项和建筑壁面静水压力项三者平衡,这使得原本针对传统浅水方程的平衡有限体积法(如表面梯度法[18])难以直接用于孔隙率浅水方程. 针对各向异性孔隙率浅水方程的平衡处理有一些初步进展,Sanders等[14]的研究忽略底坡项,近似认为建筑壁面静水压力项与通量相平衡;Guinot等[17]用各向同性浅水方程中的孔隙率梯度项直接替换各向异性孔隙率方程中的建筑壁面静水压力项,且求解时用单元的边界与体积孔隙率梯度代替各向同性时相邻单元体积孔隙率梯度.
目前,孔隙率方法和局部时间步长技术在浅水数值模拟中独立发展. 本文将这2种高效的数值模拟技术相结合,旨在建立准确高效的浅水计算模型,并将其应用于城市洪水模拟. 具体而言,在有限体积法框架下,基于能自动捕捉激波、间断的高精度HLLC(Harten-Latex-van Leer-Contact)通量计算数值格式,并结合局部时间步长技术,建立满足静水平衡特性的各向异性孔隙率浅水模型. 本文将依次介绍模型的控制方程和数值算法,模型静水平衡特性和城市洪水计算能力和效率的验证,以及网格尺寸、划分方式和局部时间步长技术对孔隙率浅水模型计算的影响.
1. 控制方程
孔隙率浅水方程可写成如下积分形式[15]:
式中:
图 1
图 1 三角网格下孔隙率的概念模型图
Fig.1 Conceptual model diagram of porosity in triangular cell
其中,
连续相位函数的值取决于水面高程与地面高程之间的关系,计算式如下:
式中:
为了更好地捕捉水流的流动方向的差异(即各向异性),将表征孔隙的连续相位函数分解为单元孔隙率
式中:
需要指出的是,令孔隙率浅水方程中所有孔隙率等于1即可得到传统浅水方程,同时额外的源项
2. 数值方法
2.1. 离散格式
式中:
需要指出的是,式(9)中右端建筑物面静水压力的离散必须考虑静水平衡特性,故未直接写出其离散式. 对于该项的处理,考虑到静水时建筑物壁面上的静水压力与单元边界上的静水压力以及底面上的作用力平衡[14],在静水条件下简化式(8)得到如下表达式:
式(10)中的边界静水压力项可离散为
式中:
式中:
式(12)在计算相邻单元之间的界面水深时,通过
式中:
当把式(10)中的项转换到边界计算时,在平衡条件下应同时满足式(14),将式(11)~(13)代入式(10),并通过式(10)~(14)得到建筑壁面的静水压力项离散形式如下:
式(9)中其他项的离散格式参考文献[20],此处不再赘述.
2.2. 局部时间步长
局部时间步长(local time step,LTS)使每个单元以其特定的时间步长更新,与全局时间步长使用同一个最小的时间步长更新相比,可大幅提高计算效率. 目前还没有在孔隙率浅水方程模拟中直接使用局部时间步长的先例,因此各向异性孔隙率结合局部时间步长来提高计算效率值得研究. 孔隙率浅水方程的时间步长计算公式采用Sanders等[14]提出的方法,同时为了避免
式中:
为了计算的稳定,在干湿边界处须作一些修正. 将干单元的
式中:
式中:
通过相邻单元的潜在局部时间步长最终得到各单元的实际局部更新时间步长:
式中:
3. 模型计算
3.1. 参数设置与网格划分
图 2
图 2 理想城市洪水实验布置图
Fig.2 Experimental set-up and dimensions for idealized urban floods
为了对比传统浅水方程和各向异性孔隙率浅水方程的计算精度和效率,以及网格划分对各向异性孔隙率方程计算结果的影响,本研究共设计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. 平衡验证
图 4
图 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类模型的计算水深和流速与实测也吻合较好(见图8、9);特别是,各向异性孔隙率方程与传统浅水方程一样能精细模拟出建筑物和街道间的水深高低起伏,甚至在流速计算上前者优于后者.
图 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的计算结果可看出,若需预测建筑街道间水深和流速的精细波动(见图8、9),孔隙率网格也不能过大,如方案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)
为了定量分析模型的误差,将实验数据和实测结果进行对比,误差计算如下:
其中:
表 1 不同方案的计算值与实测值的误差比较
Tab.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 |
图10显示了各方案在不同的局部时间步长最大更新层级情况下的计算加速效果,这里定义计算时间加速倍数a为方案C_C在
精细网格(如P_RN)的计算结果更加精确,但计算效率没有优势,粗网格(P_CN和P_CG)的计算效率更高,计算精度也令人满意. 值得注意的是,在随机网格中可能出现极小的孔隙率,需要借助工具(如ArcGIS)计算孔隙率以避免数据误差,使用样本点计算孔隙率[17]可能会降低计算稳定性. 总体而言,P_CG是最优选项,但在建筑密集情况下指定网格须作较为复杂的前处理.
4. 结 论
(1)在网格数量与传统浅水方程在同一数量级的情况下,各向异性孔隙率浅水方程可精细反映建筑街道间的水深流速沿程波动,与传统浅水方程具有同等精度;随着网格数量减少,孔隙率浅水方程精度下降,但在保证水深流速趋势和范围正确、损失局部细节精度(如建筑街道间的波动趋于平滑)的情况下,网格数量可减少到原来的1/10.
(2)在精细反映建筑街道间水深流速沿程波动的前提下,局部时间步长技术比各向异性孔隙率方法的计算效率更高,而两者的结合可进一步缩短计算时间;但在大范围城市洪水实时预测预报中,为了保证计算效率往往不得不牺牲计算精度,比如本文各向异性孔隙率浅水方程在平滑了建筑街道间的沿程波动后,计算速度比传统浅水模型(mmax=0)提高了50多倍,具有明显优势,结合局部时间步长,计算速度可进一步提高2.0~3.0倍.
需要指出的是,本文计算效率提高的认识仅基于上述算例,如果采用其他算例,模型的加速效果在定量上可能会有差异. 各向异性孔隙率浅水方程与局部时间步长技术的结合在城市洪水模拟和预报中具有广阔的应用前景,如果能结合GPU并行技术,计算效率还有更大提升空间[29].
参考文献
Applicability of the local inertial approximation of the shallow water equations to flood modeling
[J].DOI:10.1002/wrcr.20366 [本文引用: 1]
Performances and limitations of the diffusive approximation of the 2-D shallow water equations for flood simulation in urban and rural areas
[J].DOI:10.1016/j.apnum.2016.07.003 [本文引用: 1]
A coarse-grid approach to representing building blockage effects in 2D urban flood modelling
[J].
Multiple porosity shallow water models for macroscopic modelling of urban floods
[J].
Macroscopic modelling of urban floods
[J].
Steady-flow experiments in urban areas and anisotropic porosity model
[J].DOI:10.1080/00221686.2016.1238013 [本文引用: 1]
植被群影响下动床洪水数值模拟研究
[J].
Numerical simulation of dam-break flood over mobile and vegetated beds
[J].
城市化地面二维浅水模拟
[J].
2D shallow water simulation for urbanized areas
[J].
The solution of the dam-break problem in the Porous Shallow water Equations
[J].DOI:10.1016/j.advwatres.2018.01.026 [本文引用: 1]
A 1D - 2D shallow water equations solver for discontinuous porosity field based on a generalized riemann problem
[J].
Numerical simulation of dam-break flow and bed change considering the vegetation effects
[J].DOI:10.1016/j.ijsrc.2015.04.004 [本文引用: 1]
Two-dimensional shallow-water model with porosity for urban flood modelling
[J].DOI:10.1080/00221686.2008.9521842 [本文引用: 1]
An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity
[J].
Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling
[J].DOI:10.1016/j.jhydrol.2008.08.009 [本文引用: 10]
Urban flood modeling using shallow water equations with depth-dependent anisotropic porosity
[J].DOI:10.1016/j.jhydrol.2016.08.025 [本文引用: 3]
Urban flood modeling with porous shallow-water equations: a case study of model errors in the presence of anisotropic porosity
[J].DOI:10.1016/j.jhydrol.2015.01.059 [本文引用: 1]
Dual integral porosity shallow water model for urban flood modelling
[J].DOI:10.1016/j.advwatres.2017.02.009 [本文引用: 4]
The surface gradient method for the treatment of source terms in the shallow-water equations
[J].DOI:10.1006/jcph.2000.6670 [本文引用: 1]
Integration of a shallow water model with a local time step
[J].DOI:10.3826/jhr.2008.3243 [本文引用: 1]
Computationally efficient modeling of hydro-sediment-morphodynamic processes using a hybrid local time step/global maximum time step
[J].DOI:10.1016/j.advwatres.2019.03.006 [本文引用: 3]
基于局部分级时间步长方法的水沙耦合数学模拟
[J].
Coupled modeling of sediment-laden flows based on local-time-step approach
[J].
Numerical investigation of a sandbar formation and evolution in a tide-dominated estuary using a hydro-morphodynamic model
[J].DOI:10.1080/21664250.2018.1529263 [本文引用: 1]
Drag, turbulence, and diffusion in flow through emergent vegetation
[J].DOI:10.1029/1998WR900069 [本文引用: 1]
A well-balanced and fully coupled noncapacity model for dam-break flooding
[J].
A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows
[J].DOI:10.1137/S1064827503431090 [本文引用: 1]
A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography
[J].DOI:10.1016/j.cma.2013.01.015 [本文引用: 1]
Dam-break flow through an idealised city
[J].DOI:10.3826/jhr.2008.3164 [本文引用: 1]
Shallow water equations with depth-dependent anisotropic porosity for subgrid-scale topography
[J].
GPU-enhanced finite volume shallow water solver for fast flood simulations
[J].DOI:10.1016/j.envsoft.2014.02.003 [本文引用: 1]
/
〈 |
|
〉 |
