CN110990926B - 一种基于面积修正率的城市地表建筑水动力学仿真方法 - Google Patents

Abstract

本发明提供了基于面积修正率的城市地表建筑水动力学仿真方法，包括：1)基于二阶地形概化精度的非结构网格划分2)基于面积修正率的改进型二维浅水控制方程的确定；3)对步骤2)中控制方程在步骤1划分的网格上进行有限体积法离散，并将2)中提出的改进性控制方程中底坡项转换为线积分进行求解，实现城市地表建筑水动力学仿真。本发明的引入面积修正率的概念，并对引入面积修正率参数后二维浅水方程的改进，基于斜底模式的非结构网格，采用有限体积法离散求解改进后的二维浅水方程，并采用线积分的方式求解处理底坡项，并结合使用面积修正率方法形成的一种基于面积修正率的城市地表建筑水动力学仿真方法。

Description

一种基于面积修正率的城市地表建筑水动力学仿真方法

技术领域

本发明涉及水利工程仿真领域，具体地说，涉及一种基于面积修正率的城市地表建筑水动力学仿真方法。

背景技术

全球气候变化背景下，极端暴雨事件频发，城市暴雨内涝灾害也呈现增多趋势。随着城市化水平的不断提升以及受季风气候影响，我国城市暴雨洪水集中，洪涝灾害问题也非常突出。开展城市暴雨洪水方面的研究，特别是城市洪涝灾害预警预报研究，不仅具有重要的科学价值，对于支撑经济社会可持续发展也具有重要的实用价值。高精度城市暴雨内涝模型是城市洪涝预测预警和灾情评估的重要工具和基础。当模拟分析城市暴雨内涝时，地表大量的建筑物对洪水演进和积水深度的分布具有重要的影响。在大范围城市暴雨内涝模拟应用中，常用的地表建筑物处理方法有固壁边界法和加大糙率法，与固壁边界法类似的还有地形提升法。固壁边界法和地形提升法需要在网格层面来反映地表建筑物的形状和阻水作用，通常会导致网格尺寸非常小，网格数量巨大，构建的模型难以满足实时预报等时效性要求较高的应用。加大糙率法无需在网格层面反映地表建筑物的形状，但糙率取值又成为一个难题，一般不推荐采用。

因而，如何有效和简便地概化地表建筑物，高效地反映城市地表建筑对洪水演进和积水过程的影响是大尺度城市洪水模拟中的一个难题。

发明内容

本发明的目的在于提出一种高效地反映城市地表建筑对洪水演进和积水过程的影响的二维水动力学仿真方法。

为了实现上述目的，本发明提供了基于面积修正率的城市地表建筑水动力学仿真方法，该方法的实现具体包括如下步骤：

1)基于二阶地形概化精度的非结构网格划分

2)基于面积修正率的改进型二维浅水控制方程的确定；

3)对步骤2)中控制方程在步骤1)划分的网格上进行有限体积法离散，并将2)中提出的改进性控制方程中底坡项转换为线积分进行求解，实现城市地表建筑水动力学仿真。

所述步骤1)中二阶地形概化精度的网格是指将高程定义在组成网格单元的节点(即顶点)上，非结构网格包括三角形网格、四边形网格或者同时存在三角形四边形网格的混合网格。

所述步骤2)中的基于面积修正率的改进型二维浅水方程为

其中，t为时间变量；x和y为坐标；h为水深；u和v分别为x和y方向的流速；g为重力加速度；b表示单元格底高程；S_ox和S_oy分别为x和y方向的底坡项分量；S_fx和S_fy分别为x和y方向的摩阻项分量，S_cx和S_cy分别为x和y方向的局部摩阻项分量。S_φx和S_φy分别为x和y方向的修正率比降，U表示守恒变量，E表示x方向的对流能量，G表示y方向的对流能量，S表示源项，S_b表示底坡项，S_f表示摩阻项。

其中，S_cx和S_cy分别为x和y方向的局部摩阻项分量，采用下式近似计算：

类似于底坡项，S_φx和S_φy采用下式计算

所述步骤3)所述仿真方法还包括对底坡项仿真计算，所述底坡项仿真计算的公式为

其中，η为水位，b为单元格底高程。

本发明的有益效果在于：引入面积修正率的概念，并对引入面积修正率参数后二维浅水方程的改进，本发明的计算流程中方法的组合和搭配模式，即基于斜底模式的非结构网格，采用有限体积法离散求解改进后的二维浅水方程，并采用线积分的方式求解处理底坡项，并结合使用本发明提出的面积修正率方法形成的一种基于面积修正率的城市地表建筑水动力学仿真方法。

附图说明

图1为建筑物概化示意图：(a)采用固壁边界法；(b)采用建筑物孔隙率方法；

图2为实施例2中的实验现场布置；

图3为实施例2中的网格示意图(固壁边界法)；

图4为实施例2中的网格示意图(孔隙率法)；

图5为实施例2中的实验入流流量过程；

图6为实施例2中的P10测点处两种方法计算水深与试验值对比；

图7为实施例2中的P9测点处两种方法计算水深与试验值对比；

图8为实施例2中的P4测点处两种方法计算水深与试验值对比；

图9为实施例2中的P3测点处两种方法计算水深与试验值对比。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

本发明的实施提供了基于面积修正率的城市地表建筑水动力学仿真方法，该方法的实现具体包括如下步骤：

第一步基于二阶地形概化精度的网格离散

在构建地表二维洪水模型之前，一般需要将研究区域划分为若干个网格单元进行计算，目前多采用不规则的非结构网格。现有的基于非结构混合网格的二维洪水模型，由于底坡项处理困难，一般采用平底模型，即将高程定义在单元中心，假定整个网格单元是平的。采用平底模型，方法相对简单，但对地形概化上存在着一定的不足，地形概化精度仅为一阶，有些特殊地形难以有效地反映。为了避免这一问题，本发明将高程定义在单元节点上，即每个单元节点都代表一个高程点，地形概化具有二阶精度，并采用改进型的控制方程，成功解决了底坡项求解难题。本发明所提出的方法和步骤均能够应用于三角形网格、四边形网格或三角形与四边形混合的网格。

第二步面积修正率方法的引入与方程的改进

城市地表水流一般都属于浅水，控制方程可以采用二维浅水方程，该方程属于简化的Navier-Stokes方程，一般用来描述具有自由表面的浅水体在重力作用下的流动，可以写成如下所示的矩阵形式：

式中：t为时间；x和y均为笛卡尔坐标；h为水深；u为x方向的流速；v为y方向的流速；g为重力加速度；S_ox和S_oy分别为x和y两个方向的底坡项分量；S_fx和S_fy分别为x和y两个方向的摩阻坡降，U表示守恒变量，E表示x方向的对流能量，G表示y方向的对流能量。

当模拟分析城市暴雨内涝时，地表大量的建筑物对洪水演进和积水深度的分布具有重要的影响。一般来说，最为精确的方法是采用固壁边界法，即在划分网格时，将这些建筑物当作不过水区域，如图1(a)所示，这种方法在小区域模拟中应用较多。当模拟大型区域城市暴雨内涝时，采用固壁边界法时，由于需要较小尺寸的网格来反映地表建筑物的形状，一方面会增加网格数量，另一方面，时间步长会极大的受限于最小尺寸的网格，从而导致模型计算效率低下，特别是在实时预报时，计算效率可能无法满足要求。

为了解决这一问题，本发明在控制方程中引入一个反映地表建筑物的参数对网格面积进行修正，采用下式计算：

式中：φ为建筑物面积修正率；A_c为某一网格内建筑物的面积；A_i为计算网格的面积。

则二维浅水方程可改写为：

式中：S_φx和S_φy分别为x和y方向的修正率比降，S_b表示底坡项，S_f表示摩阻项，S_cx和S_cy分别为x和y方向的局部摩阻项分量，采用下式近似计算：

类似于底坡项，S_φx和S_φy采用下式计算

采用引入面积修正率后的控制方程，在考虑城市地表建筑物时，在网格层面无需特别考虑建筑物的分布和形状，因而网格的尺寸可以划分得更大和更均匀，可以明显地提高计算效率以及减少前期划分网格时的准备时间。

在数值离散方程时采用斜底模式的非结构网格单元以及中心型底坡项计算方法时，如果采用由式(4)和式(5)描述的二维浅水方程，则需要构造动量能量校正项才能保证模型的和谐性，即在静水条件下，模型计算结果能维持流速为零和水深不变的状态。为避免计算修改项，对控制方程进行改时，即在保证方程两边相等的情况下，将对流项中部分组分移动至底坡项中，则引入面积修正率参数后的控制方程可以进一步变为如下形式：

其中，b表示单元格底高程。

第三步有限体积离散与底坡项处理

基于非结构的混合网格，采用单元中心格式(即将守恒变量定义在单元中心)的有限体积法对式(8)进行积分，可以得到：

式中：V为任意控制体。

利用Green公式，将式(10)中的体积分化为沿周边的线积分，则式(3-15)可改写为：

式中：A_i为单元的面积；i为单元编号；n为向外的单元法向量；k是单元边的序号；l为单元边的长度；F_i,k为数值通量。

界面通量采用Godunov型的HLLC格式计算：

式中：S_L、S_M和S_R分别为左波、接触波和右波的波速估计值；F_L和F_R由式(13)计算；

和

由式(14)给出：

式中：φ为建筑物面积修正率；u_N为垂直于单元界面处的流速；u_T,L和u_T,R分别界面左侧和右侧平行于界面的流速；E₁和E₂为向量E_HLL的前两项，E_HLL为采用HLL方程计算的法向通量。

基于非结构网格的二维水动力学模型底坡项的处理一直是计算浅水力学方面的热点问题，为了保证模型的和谐性，一般需要进行特殊处理，常用的处理方法包括水面梯度法、水位方程法、底坡项分解法、底坡项特征分解法等。上述这些方法有各自的特点和适用范围，但均难以应用在斜底模式的非结构混合网格上或是存在着算法复杂等问题，特别是在引入面积修正率这一参数后，问题变得更为复杂。针对这一问题，本发明采用改进形式的控制方程(详见第一步)和将底坡项转换为线积分的求解方式，可以严格保证模型的和谐性。由改进后的二维浅水方程可知，底坡项引起的源项如下：

对上式右端进行合并，可以得到：

对式(16)在网格单元V上进行积分，并转换为线积分，可以得到：

式中：η为水位，b为单元格底高程。

考虑静水情况下，流速为0，单元的对流通量为

将h_i,k＝η_i-b_i,k代入上式，则有

由上式可知，在静水条件下，对流项始终与底坡项相等，不会产生虚假流速，满足和谐性的要求。

实施例2

为了验证本发明提出基于面积修正率方法来考虑城区地表建筑物的精度和可靠度，采用被广泛应用于检验模型在城市密集建筑区域适用性的意大利城区洪水实验中一组数据对本发明提出的方法进行验证和应用。如图2所示，该实验在河道中间设置一系列方形阻水建筑物，并测量了一系列不同入流情况下一些特定测点的水深变化过程。分别采用了两种不同的方法来模拟本算例：固壁边界法和面积修正率法。固壁边界法把阻水建筑当作不过水区域，在划分网格时予以考虑，这种方法相对比较精确，但建筑区域需要将网格划分得很密集来反映过水建筑物的形状，因而网格尺寸相对比较小，网格数量也一般较多，计算效率比较低；面积修正率法，即采用本发明提出的概化方法，在划分网格时不考虑建筑物的影响，但将建筑物的面积在网格面积上予以修正。

采用固壁边界法时，整个研究区域被划分成7381个非结构四边形网格，共15043条网格边，其中建筑分布区域进行网格加密处理(如图3所示)。采用面积修正率法时，整个研究区域被划分成3484个非结构四边形网格，共7089条网格边(如图4所示)，建筑物分布区域不作特殊处理。本次模拟计算模型采用的入流过程如图5所示。采用上述两种不同的概化建模方法，分别模拟计算了60s，并将计算结果与实测值以及Kim等人的计算结果进行对比。选取了4个典型点P3、P4、P9和P10进行验证与对比，对比结果分别如图6至图9所示。

对比分析采用面积修正率法和固壁边界法两种不同建筑物概化方式的计算结果，可以明显看出，两种概化方法的计算结果非常接近，表明本发明提出的地表面积修正率方法结果是可靠的，可以在以较少的网格数量和较大网格尺寸的情况下，获得与固壁边界法较为相近的计算精度。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.基于面积修正率的城市地表建筑水动力学仿真方法，其特征在于，包括：

1)基于二阶地形概化精度的非结构网格划分；

2)基于面积修正率的改进型二维浅水控制方程的确定；

3)对步骤2)中控制方程在步骤1)划分的网格上进行有限体积法离散，并将步骤2)中的改进型二维浅水控制方程中底坡项转换为线积分进行求解，实现城市地表建筑水动力学仿真；所述步骤2)中的基于面积修正率的改进型二维浅水方程为

其中，t为时间变量；x和y为坐标；h为水深；u和v分别为x和y方向的流速；g为重力加速度；b表示单元格底高程；S_ox和S_oy分别为x和y方向的底坡项分量；S_fx和S_fy分别为x和y方向的摩阻项分量，S_cx和S_cy分别为x和y方向的局部摩阻项分量；S_φx和S_φy分别为x和y方向的修正率比降，U表示守恒变量，E表示x方向的对流能量，G表示y方向的对流能量，S表示源项，S_b表示底坡项，S_f表示摩阻项。

2.根据权利要求1中的动力学仿真方法，其特征在于，所述步骤2)中国，S_cx和S_cy分别为x和y方向的局部摩阻项分量，通过如下方程计算：

类似于底坡项，S_φx和S_φy采用下式计算：

3.根据权利要求1中的动力学仿真方法，其特征在于，所述步骤1)中二阶地形概化精度的网格是指将高程定义在组成网格单元的节点上，非结构网格包括三角形网格、四边形网格或者同时存在三角形四边形网格的混合网格。

4.根据权利要求2中的动力学仿真方法，其特征在于，所述仿真方法还包括对底坡项仿真计算。

5.根据权利要求4中的动力学仿真方法，其特征在于，所述底坡项仿真计算的公式为

其中，η为水位，b为单元格底高程。

Images