CN110147646B - 一种数值模拟框架下线性挡水构筑物的过流处理方法 - Google Patents

Abstract

本发明公开了一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法。获取计算区域的地形、地貌数据和线性挡水构筑物的几何特征数据后，采用四边形非结构网格离散计算区域，线性挡水构筑物的中心线作为网格剖分的控制线。网格剖分完成后，将代表线性挡水构筑物的网格单元边定义为特殊边，其余的单元边则定义为边界边和普通边。根据线性构筑物的顶部宽度和其上游水深的比值来判断堰流形态，采用对应的堰流公式计算通过特殊边的单宽流量，得出通过该特殊边外法线方向的质量通量和动量通量。该处理方法在兼顾计算效率的同时，能够考虑线性构筑物两侧的质量交换和动量交换过程，弥补了现有方法的不足。

Description

一种数值模拟框架下线性挡水构筑物的过流处理方法

技术领域

本发明属于水利工程技术领域，特别涉及一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法。

背景技术

所谓线性挡水构筑物指的是在计算区域内经常遇到的高速公路路基、铁路路基、挡水堤等构筑物，这些构筑物的顶部通常高于地面，在空间形态上呈线性特征。当洪水发生时，这些构筑物会起到明显的阻水作用，当一侧水位大于顶部高程时，会发生过流。

在工程实践中，通常采用二维浅水模型来模拟大范围的洪水水流运动过程。在洪水二维数值模拟计算中，由于非结构网格良好的边界适应特征，使其成为当前洪水数值模拟的主流网格。在二维非结构数值模拟框架下如何精确处理这些线性阻水构筑物的影响，直接关系到洪水模型计算结果的可靠程度。当前，在二维非结构数值模拟框架下处理这些线性构筑物时有两种常用的方法，一种是采用精细网格将线性构筑物进行剖分离散，这需要很精细的地形数据来支撑，但是当进行大范围的洪水分析时，一方面获取这种精细地形的成本太高，另一方面，虽然这种处理的数值精度较高，但是同时会带来灾难性的计算效率问题，所以在实际应用中较少使用；另一种是非常简化的处理方法，通过宽顶堰流公式将通过线性挡水构筑物顶部的水流直接汇入下游网格，仅考虑了线性挡水构筑物两侧的水量交换，而忽略了两侧的水流动量交换，通常也没有区分线性挡水构筑物顶部不同的堰流特征，这些都会给实际的计算精度带来很大的影响。

发明内容

本发明的目的是提供一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法，该方法区别于现有的处理方法，能够区分线性挡水构筑物顶部过流的堰流形态，不仅考虑挡水构筑物两侧的流量交互，同时能够考虑两侧的动量交互，并给出了线性挡水构筑物两侧数值通量的计算公式。该方法在兼顾计算效率的同时，能够提高现有处理方法的计算精度。

本发明的目的是通过以下方案实现的：

本发明为一种二维非结构洪水数值模拟框架下线性挡水构筑物过流处理方法，充分考虑计算区域内线性挡水构筑物的过流特点，将线性挡水构筑物概化为二维非结构网格单元的特殊边，根据线性挡水构筑物顶部过流特征的不同，采用不同的堰流流量系数计算通过该特殊边的过流流量，之后根据该流量值计算通过该特殊边的流量通量和动量通量，保证通过所有网格单元边的数值通量在存储形式上的统一；其具体的方法步骤如下：

1)获取二维计算区域的全要素DLG数据(含地形数据和地貌类型数据)，比例尺1:10000以上(含1:10000)；从设计图纸或通过现场测绘的手段获取线性挡水构筑物(如高速公路路基、铁路路基、挡水堤等)的几何信息，主要包括中心线坐标、顶部高程Z_w和顶部宽度D_w等信息。

2)计算区域离散：线性挡水构筑物的中心线作为网格剖分时的控制线，采用四边形非结构网格单元离散计算区域，剖分完成后，线性挡水构筑物中心线离散成的网格单元边称为特殊边，计算区域边界线离散成的网格单元边称为边界边，其它边称为普通边。

3)计算区域赋值：采用地形数据对网格形心点进行高程插值，获取网格单元的高程数据；对特殊边进行属性赋值，主要包括线性挡水构筑物的顶部高程和顶部宽度；根据地貌类型数据对网格单元进行糙率赋值；对各单元水力要素变量赋初值。

4)确定时间步长dt：二维非结构洪水数值模拟方法采用显格式进行，计算时间步长选取受CFL(Courant-Friedrichs-Lewy)条件限制。

5)边界边数值通量计算：根据边界边的过流类型(流量边界、水位边界、固壁边界、自由出流边界等)，采用在边界边外侧布设虚单元的方法来计算t时刻通过边界边外法线方向的数值通量。

6)普通边数值通量计算：采用完整二维浅水方程组描述洪水运动过程，采用Roe格式的近似Riemann解来显式计算t时刻通过普通边外法线方向的数值通量。

7)特殊边数值通量计算：当特殊边对应的线性挡水构筑物顶部不过流时，通过该特殊边的质量通量和动量通量均为0；当特殊边对应的线性挡水构筑物顶部过流时，采用堰流公式来计算t时刻通过该特殊边的单宽流量q；根据该特殊边两侧单元的水位判断上下游关系，水位高的一侧定义为上游，水位低的一侧定义为下游，用q除以下游单元t时刻的水深h_d得到t时刻通过该特殊边外法线方向的流速u_d；之后用h_d和u_d获得通过该特殊边t时刻外法线方向的数值通量，如式(1)所示：

其中，F^*为通过特殊边外法线方向的数值通量，n_x，n_y分别为特殊边外法线方向的单位向量在x方向和y方向的分量，g为重力加速度。

8)对每个网格单元4条边的外法线方向数值通量进行求和计算，更新网格单元的水力要素值到t+dt时刻。

9)令t＝t+dt，重复步骤4)-8)，直到计算结束。

进一步的，步骤7)特殊边数值通量计算中，当特殊边对应的线性挡水构筑物顶部过流时，采用通用的堰流公式计算通过该特殊边的单宽流量q：

式中，m为堰流流量系数，g为重力加速度，H₀为线性挡水构筑物顶部有效水头。根据线性挡水构筑物顶部宽度D_w与其上游水深h_u的比值来确定m的取值，D_w/h_u≤0.67时，m取薄壁堰流量系数；当0.67＜D_w/h_u≤2.5，m取实用堰流量系数；D_w/h_u＞2.5时，m取宽顶堰流量系数。

本发明的有益效果：在该发明中，将充分考虑线性挡水构筑物的空间分布特点，将其几何特征概化在非结构网格单元边上，根据线性挡水构筑物顶部宽度与上游水深的比值判断堰流形态，选取不同的堰流流量系数来计算通过挡水构筑物顶部的单宽流量，根据该单宽流量值推导得出通过该特殊边外法线方向的数值通量，不仅考虑了质量交互，同时考虑了动量交互。该处理方法能够很好的兼顾模型整体计算效率和计算精度，弥补现有方法的不足。

附图说明

图1为本发明的一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法流程图。

具体实施方式

下面结合附图1对本发明作进一步的描述。

本发明提供的方法是一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法。将线性挡水构筑物概化为二维非结构单元的特殊边，将其几何特征参数赋值为特殊边的属性，根据构筑物顶部不同的堰流形态选择不同类型的堰流计算公式，如薄壁堰、实用堰和宽顶堰公式。根据堰流公式计算得出的单宽流量推导得出通过该特殊边外法线方向的数值通量，包括质量通量和动量通量，保证通过所有网格单元边的数值通量在存储形式上的统一。该方法的具体过程如下：

1)获取二维计算区域的全要素DLG数据(含地形数据和地貌类型数据)，比例尺1:10000以上(含1:10000)；从设计图纸或通过现场测绘的手段获取线性挡水构筑物(如高速公路路基、铁路路基、挡水堤等)的几何信息，主要包括中心线坐标、顶部高程Z_w和顶部宽度D_w等信息。

2)计算区域离散：线性挡水构筑物的中心线作为网格剖分时的控制线，采用四边形非结构网格单元离散计算区域，剖分完成后，线性挡水构筑物中心线离散成的网格单元边称为特殊边，计算区域边界线离散成的网格单元边称为边界边，其它边称为普通边。

3)计算区域赋值：采用地形数据对网格形心点进行高程插值，获取网格单元的高程数据；对特殊边进行属性赋值，主要包括线性挡水构筑物的顶部高程和顶部宽度；根据地貌类型数据对网格单元进行糙率赋值；对各单元水力要素变量赋初值。

4)确定时间步长dt：二维非结构洪水数值模拟方法采用显格式进行，计算时间步长选取受CFL(Courant-Friedrichs-Lewy)条件限制，限制条件如下：

其中，其中u，v为x，y方向的流速分量，h为水深，g为重力加速度，N_cfl为CFL数，推荐在计算中取为0.8，dt为计算时间步长，L_L,LR为四边形网格单元中心到对应边中点之间的距离。

5)边界边数值通量计算：根据边界边的过流类型(流量边界、水位边界、固壁边界、自由出流边界等)，采用在边界边外侧布设虚单元的方法来计算t时刻通过边界边外法线方向的数值通量，具体的求解方法可参见如下论文(Ying X Y,Jorgeson J and Wang S SY.Modeling dam-break flows using finite volume method on unstructuredgrid.Engineering application of computational fluid mechanics,2009,3(2):184-194.)。

6)普通边数值通量计算：采用完整二维浅水方程组描述洪水运动过程，采用的守恒型方程组如下：

其中：

_h为水深，u，_v分别为x，_y方向的流速，t为时间；和分别为x方向和y方向的单位向量；

分别为x，y方向的坡度，_Zb为地面高程，g为重力加速度；

分别为x，y方向的摩阻项，其中n为Manning糙率系数。

采用Roe格式的近似Riemann解来显式计算t时刻通过普通边外法线方向的数值通量，具体的求解过程可参见如下著作P73-P75页(张大伟.基于Godunov格式的堤坝溃决水流数值模拟,2014,中国水利水电出版社)。

7)特殊边数值通量计算：当特殊边对应的线性挡水构筑物顶部不过流时，通过该特殊边的质量通量和动量通量均为0；当特殊边对应的线性挡水构筑物顶部过流时，采用堰流公式来计算t时刻通过该特殊边的单宽流量q，采用的通用堰流公式如下：

式中，m为堰流流量系数，g为重力加速度，H₀为线性挡水构筑物顶部有效水头。根据线性挡水构筑物顶部宽度D_w与其上游水深h_u的比值来确定m的取值，D_w/h_u≤0.67时，m取薄壁堰流量系数；当0.67＜D_w/h_u≤2.5，m取实用堰流量系数；D_w/h_u＞2.5时，m取宽顶堰流量系数。各类堰流形态下，m的具体取值可参考如下教材P6-P35页(李家星,赵振兴.水力学(下册),2001,河海大学出版社)。

根据该特殊边两侧单元的水位判断上下游关系，水位高的一侧定义为上游，水位低的一侧定义为下游，用q除以下游单元t时刻的水深h_d得到t时刻通过该特殊边外法线方向的流速u_d；之后用h_d和u_d获得通过该特殊边t时刻外法线方向的数值通量，如下式所示：

其中，F^*为通过特殊边外法线方向的数值通量，n_x，n_y分别为特殊边外法线方向的单位向量在x方向和y方向的分量，g为重力加速度。

8)对每个网格单元4条边的外法线方向数值通量进行求和计算，更新网格单元的水力要素值到t+dt时刻。

9)令t＝t+dt，重复步骤4)-8)，直到计算结束。

上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。

Claims (2)

1.一种二维非结构洪水数值模拟框架下线性挡水构筑物过流处理方法，其特征在于：包括以下步骤：

1)获取二维计算区域的全要素DLG数据：所述全要素DLG数据含地形数据和地貌类型数据；从设计图纸或通过现场测绘的手段获取线性挡水构筑物的几何信息，所述几何信息包括中心线坐标、顶部高程Z_w和顶部宽度D_w，所述线性挡水构筑物包括高速公路路基、铁路路基、挡水堤；

2)计算区域离散：线性挡水构筑物的中心线作为网格剖分时的控制线，采用四边形非结构网格单元离散计算区域，剖分完成后，线性挡水构筑物的中心线离散成的网格单元边称为特殊边，计算区域边界线离散成的网格单元边称为边界边，其他边称为普通边；

3)计算区域赋值：采用地形数据对网格形心点进行高程插值，获取网格单元的高程数据；对特殊边进行属性赋值，包括线性挡水构筑物的顶部高程和顶部宽度；根据地貌类型数据对网格单元进行糙率赋值；对各单元水力要素变量赋初值；

4)确定时间步长dt：二维非结构洪水数值模拟方法采用显格式进行，计算时间步长选取受CFL条件限制；

5)边界边数值通量计算：根据边界边的过流类型，采用在边界边外侧布设虚单元的方法来计算t时刻通过边界边外法线方向的数值通量；所述边界边的过流类型包括流量边界、水位边界、固壁边界、自由出流边界；

6)普通边数值通量计算：采用完整二维浅水方程组描述洪水运动过程，采用Roe格式的近似Riemann解来显式计算t时刻通过普通边外法线方向的数值通量；

7)特殊边数值通量计算：当特殊边对应的线性挡水构筑物顶部不过流时，通过该特殊边的质量通量和动量通量均为0；当特殊边对应的线性挡水构筑物顶部过流时，采用堰流公式来计算t时刻通过该特殊边的单宽流量q；根据该特殊边两侧单元的水位判断上下游关系，水位高的一侧定义为上游，水位低的一侧定义为下游，用q除以下游单元t时刻的水深h_d得到t时刻通过该特殊边外法线方向的流速u_d；之后用h_d和u_d获得通过该特殊边t时刻外法线方向的数值通量，如式(1)所示：

式中，F^*为通过特殊边外法线方向的数值通量，n_x，n_y分别为特殊边外法线方向的单位向量在x方向和y方向的分量，g为重力加速度；

8)对每个网格单元4条边的外法线方向数值通量进行求和计算，更新网格单元的水力要素值到t+dt时刻；

9)令t＝t+dt，重复步骤4)-8)，直到计算结束。

2.根据权利要求1所述的线性挡水构筑物过流处理方法，其特征在于：步骤7)特殊边数值通量计算中，当特殊边对应的线性挡水构筑物顶部过流时，采用通用的堰流公式计算t时刻通过该特殊边的单宽流量q：

式中，m为堰流流量系数，g为重力加速度，H₀为线性挡水构筑物顶部有效水头；根据线性挡水构筑物顶部宽度D_w与其上游水深h_u的比值来确定m的取值，D_w/h_u≤0.67时，m取薄壁堰流量系数；当0.67＜D_w/h_u≤2.5，m取实用堰流量系数；D_w/h_u＞2.5时，m取宽顶堰流量系数。