一种迎风守恒型河道漫溢出流数值模拟方法

技术领域

本发明涉及水利工程领域，尤其涉及防洪减灾领域，具体为一种迎风守恒型河道漫溢出流数值模拟方法。

背景技术

随着全球气候模式的变化，极端强降雨事件频繁发生，时常出现河道超标准洪水漫溢的现象，洪水漫溢后，会给淹没区的居民带来严重的生命和财产损失。数学模型作为一种有用的工具，可以模拟河道洪水的漫溢过程，为防洪决策部门提供有价值的决策信息。一维数学模型，具有很高的计算效率，在河道水流模拟时可以提供很好的计算结果，但是洪水一旦漫溢后，河道外水流呈现明显的二维运动特征，此时一维模型很难再适用。如果河道内外水流过程模拟统一采用二维数学模型进行，无疑会获得很高的计算精度，但是当前条件下，二维模型的计算效率很难满足工程上实时分析计算的需求。鉴于此，将一维模型和二维模型耦合起来，可以充分利用一维模型的计算效率和二维模型的计算精度。

当前，工程上常用的一、二维模型耦合方式是采用堰流公式的形式来计算河道内外水流的交互，然而该方法无法表达一、二维模型间的动量交换，同时堰流公式中的计算系数选取存在很大的不确定性，这些给洪水漫溢出流的模拟带来一些误差，直接影响漫溢淹没范围和淹没过程的准确程度。另外，当前工程上经常用到的无激波捕捉能力的数值离散格式，如河道模拟时常用到的Preissmann格式，在陡峭地形及复杂流态的条件下，常不能获得正确的计算结果，通常在大江大河的水流模拟中适用效果较好，但是用于中小河流河道漫溢洪水过程模拟时常常失败。

发明内容

本发明的目的在于提供一种高适用性、高精度的河道漫溢出流数值模拟方法，既可以用于大江大河的漫溢出流模拟，也可以用于大坡度山地丘陵区河道漫溢出流的模拟，河道内外的水流交互计算采用具有严格物理数学意义的特征不变量理论完成，改善传统一、二维模型耦合方法的精度。

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

一种迎风守恒型河道漫溢出流数值模拟方法，河道水流采用一维数学模型模拟，河水漫溢出槽后采用二维数学模型模拟，一维数学模型和二维数学模型均采用基于有限体积法的迎风守恒型数值格式离散，水力变量储存在单元中心，河道内外的水流交互通过特征不变量理论进行求解，以保证整套数值离散方法具有很好的精度和数值守恒特性；其具体的方法步骤如下：

1)基础数据获取：获取详细的河道断面数据和河道平面几何形状数据；获取河道外淹没区详细的地形数据。

2)计算区域离散：采用一维有限体积单元离散河道；采用非结构三角形单元离散河道外淹没区域；在漫溢段，需建立一维河道单元和与其相连接的三角形单元的拓扑对应关系。

3)模型初始化：给一维河道单元设置糙率和水力条件初始值；插值网格节点和网格型心点的高程值，给二维网格单元设置糙率和水力条件初始值。

4)获取时间步长dt：由于一、二维数学模型均采用显格式数值离散方法，因此时间步长选取受到CFL条件限制，分别计算一维模型的时间步长dt₁和二维模型的时间步长dt₂，为保持计算整体稳定，取dt₁和dt₂数值较小者为计算所需的时间步长dt，即dt＝min(dt₁，dt₂)。

5)漫溢段水流交互计算：首先判断漫溢段河道内外是否有水流交互，如果有交互，借助特征不变量计算公式，通过一维河道单元的水位值和与其相邻的二维网格单元的水深和流速值计算通过各二维网格单元连接边外法线方向的质量通量和动量通量；如果水流无交互则质量通量和动量通量均为0。

6)一维模型计算：采用完整一维浅水方程组模拟河道水流运动，采用基于有限体积法的迎风守恒型HLL格式计算一维河道单元界面数值通量；步骤5)漫溢段的交互水量在一维模型中作源项处理；计算完成后，一维河道各单元的水力变量更新到t+dt时刻；

7)二维模型计算：采用完整二维浅水方程组模拟河道外水流运动过程，采用基于有限体积法的迎风守恒型Roe格式计算二维单元的界面通量；步骤5)计算的通过每个二维网格单元连接边的质量通量和动量通量作为二维模型的边界条件；计算完成后，二维网格单元的水力变量更新到t+dt时刻。

8)令t＝t+dt，重复步骤4)-7)，直至完成整个计算过程。

进一步的，步骤2)中，在河道漫溢段，一维河道单元称为一维连接单元，单个一维连接单元一般对应多个二维网格单元，二维网格单元与一维连接单元相连接的边称为连接边，每个通过连接边与河道单元直接相连的二维网格单元称为二维连接单元，二维连接单元不能跨两个一维河道单元，即一维河道单元连接边的长度等于与其对应的二维连接单元连接边的长度之和。

进一步的，步骤2)中，在河道漫溢段，需要加密一维河道单元数，同时需要加密二维连接单元数；

进一步的，步骤2)实际应用中漫溢段河道单元长度推荐取不大于200m为宜；

进一步的，步骤2)中二维网格单元连接边长视漫溢处的地形地貌而定，如果河道无堤防，该处的二维网格单元连接边长推荐取小于河道单元连接边长的1/3为宜，如果有河道堤防，二维网格单元的连接边不宜大于堤防顶宽的一半。

进一步的，步骤5)中，当河道内外有水流交互时，通过具有严格物理数学基础的特征不变量公式来计算通过二维网格单元连接边的质量通量和动量通量，特征不变量的计算公式如式(1)：

u_Ln+2c_L＝u_*+2c_*>

式中：u_*,c_*为二维网格单元连接边上的外法线方向流速和波速，u_Ln,c_L为二维网格单元中心水力要素值在连接边外法线方向的流速和波速。

上式中，u_Ln和c_L是t时刻的已知量，c_*可以通过一维河道单元的水位值求得，根据式(1)可直接求出二维连接边上的外法线方向流速u_*。获取u_*后，结合一维河道单元的水位值，可以根据二维浅水方程组直接得出通过二维连接边外法线方向的质量通量和动量通量。

进一步的，步骤6)中，采用完整的一维浅水方程组为守恒形式，具体如式(2)所示：

其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：

式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f₁和f₂分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n²Q|Q|)/(A²R^4/3)，R为水力半径，n为Manning糙率系数，q_l为质量源项。

进一步的，当一维河道单元内外有水流交互时，q_l的求解方法如式(3)所示：

式中：q_j为与该一维河道单元相连接的第j个二维单元连接边外法线方向的质量通量，s_j为第j个二维单元连接边的长度，k为该一维河道单元对应的二维连接单元的个数。

进一步的，步骤1)中获取河道断面数据和河道平面几何形状数据，河道断面间距小于1km；获取河道外淹没区地形数据，比例尺精度大于1:10000。

Godunov格式的迎风守恒型数值计算格式可以很好的适应各种复杂地形，精确的捕捉自由水面，可以完美的取代传统的浅水数值模拟方法。将这种迎风守恒型格式用于河道漫溢出流过程的模拟，具有更广泛的适用性，一方面可适应坡度平缓的大江大河洪水过程模拟，同时可用于坡度陡峭的中小河流洪水过程的模拟。为克服传统的采用堰流公式耦合一、二维模型的缺陷，采用具有严格物理数学基础的特征不变量理论来处理一、二维模型间的水流交互问题。本方法可为河道漫溢出流计算提供一种新的高适应性和高精度的解决方案。

本发明的有益效果：

在进行河道漫溢出流分析计算时，采用迎风守恒型的数值格式使得模型具有更广泛的适用性，不仅能够用于坡度平缓的大江大河漫溢出流过程的模拟，也能用于山地丘陵区大坡度中小河流漫溢过程的模拟；另外，采用具有严格物理数学基础的特征不变量理论处理河道内外的水流交互，可以克服采用堰流公式的传统耦合方式的缺陷，提高模拟精度。综上，本方法能够扩大现有河道漫溢出流数学模型的适用范围、提升数值精度。

附图说明

图1为本发明的一种迎风守恒型河道漫溢出流数值计算方法流程图；

图2为一、二维数学模型耦合示意图，i为一维连接单元，k₁、k₂、k₃为二维连接单元，s₁、s₂、s₃为二维连接单元的连接边。

具体实施方式

下面结合附图1和附图2对本发明作进一步说明。

该发明提供的是一种迎风守恒型河道漫溢出流数值计算方法。一方面，迎风守恒型的数值计算方法具有很广泛的适用性，不仅适用于坡度平缓的大江大河洪水漫溢过程模拟，也适用于大坡度山地丘陵区的中小河流漫溢洪水过程模拟；另一方面，为保证模型的整体计算精度和守恒特性，河道内外的水流交互过程计算采用具有严格物理数学基础的特征不变量理论进行。该方法具体过程如下：

1)基础数据获取：获取详细的河道断面数据和河道平面几何形状数据，河道断面间距需小于1km；获取河道外淹没区详细的地形数据，比例尺精度需大于1:10000。

2)计算区域离散：采用一维有限体积单元离散河道；采用非结构三角形单元离散河道外淹没区域。在漫溢段，需建立一维河道单元和与其相连接的三角形单元的拓扑对应关系。漫溢段的一维河道单元称为一维连接单元，单个一维连接单元一般对应多个二维网格单元，二维网格单元与一维连接单元相连接的边称为连接边，每个通过连接边与河道单元直接相连的二维网格单元称为二维连接单元，二维连接单元不能跨两个一维河道单元二维网格单元，即一维河道单元连接边的长度等于与其对应的二维连接单元连接边的长度之和。在附图2中，一维河道单元i处于漫溢段为一维连接单元，与其相连接的二维连接单元为k₁、k₂和k₃，对应的连接边分别为s₁、s₂和s₃。为保证计算精度，在河道漫溢段，需要加密河道单元数，同时需要加密二维连接单元数；实际应用中漫溢段河道单元长度推荐不大于200m左右为宜，二维网格单元连接边长视漫溢处的地形地貌而定，如果河道无堤防，该处的二维网格单元连接边推荐取小于河道单元边长的1/3为宜，如果有河道堤防，三角形单元的连接边不宜大于堤防顶宽的一半。

3)模型初始化：给一维河道单元设置糙率和水力条件初始值；插值网格节点和网格型心点的高程值，给二维网格单元设置糙率和水力条件初始值。

4)获取时间步长dt：由于一、二维数学模型均采用显格式数值离散方法，因此时间步长选取受到CFL条件限制，分别计算一维模型的时间步长dt₁和二维模型的时间步长dt₂，为保持计算整体稳定，取dt₁和dt₂数值较小者为计算所需的时间步长dt，即dt＝min(dt₁，dt₂)。

5)漫溢段水流交互计算：首先判断漫溢段河道内外是否有水流交互，如果无交互则质量通量和动量通量均为0；如果有交互，借助具有严格物理数学基础的特征不变量计算公式，通过一维单元的水位值和与其相邻网格单元的水深和流速值计算通过各网格单元连接边外法线方向的质量通量和动量通量，计算公式如式(1)：

u_Ln+2c_L＝u_*+2c_*>

u_*,c_*为二维网格单元连接边上的外法线方向流速和波速，u_Ln,c_L为二维网格单元中心水力要素值在连接边外法线方向的流速和波速。

上式中，u_Ln和c_L是t时刻的已知量，c_*可以通过一维河道单元的水位值求得，根据式1)可直接求出二维连接边上的外法线方向流速u_*。获取u_*后，结合一维河道单元的水位值可以直接得出通过二维连接边外法线方向的质量通量和动量通量。

6)一维模型计算：采用完整一维浅水方程组模拟河道水流运动，完整的一维浅水方程组为守恒形式，具体如式(2)所示：

其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：

式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f₁和f₂分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n²Q|Q|)/(A²R^4/3)，R为水力半径，n为Manning糙率系数，q_l为质量源项。当一维河道单元内外有水流交互时，q_l的求解方法如式(3)所示：

式中：q_j为与该一维河道单元相连接的第j个二维单元连接边外法线方向的质量通量，s_j为第j个二维单元连接边的长度，k为该一维河道单元对应的二维连接单元的个数。

采用基于有限体积法的迎风守恒型HLL格式计算一维河道单元界面数值通量，采用TVD—MUSCL格式进行变量的空间数值重构，采用Hancock预测校正格式将时间精度提高到二阶，具体数值离散过程可参考如下文献描述(张大伟，基于Godunov格式的堤坝溃决水流数值模拟，2014，中国水利水电出版社)。

步骤5)漫溢段的交互水量在一维模型中作源项处理；计算完成后，一维河道各单元的水力变量更新到t+dt时刻。

7)二维模型计算：采用完整二维浅水方程组模拟河道外水流运动过程，采用的完整二维浅水方程的守恒形式如式(4)：

式中：h为水深，u，v分别为x，y方向的流速，t为时间，分别x，y方向的坡度，Z_b为地面高程，g为重力加速度，分别为x，y方向的摩阻项，其中n为Manning糙率系数。

采用基于有限体积法的二维迎风守恒型Roe格式离散二维浅水方程组，采用三角形三个顶点高程数据直接积分计算底坡项，采用半隐的数值离散格式处理摩阻项，具体数值离散过程可参考如下文献描述(张大伟,权锦,马建明,向立云.基于Godunov格式的流域地表径流二维数值模拟.水利学报,2018,49(7):787-794.)。

步骤5)计算的通过每个二维网格单元连接边的质量通量和动量通量作为二维模型的边界条件；计算完成后，二维网格单元的水力变量更新到t+dt时刻。

8)令t＝t+dt，重复步骤4)-7)，直至完成整个计算过程。

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