CN114091163B - 一种基于有限体积法的河道构筑物过流数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于有限体积法的河道构筑物过流数值模拟方法。首先获取河道的平面几何数据、河道断面数据和构筑物几何尺寸，河道采用一维有限体积单元离散，水力要素值存储在单元中心，构筑物所在位置为单元界面。根据构筑物过流公式计算每一时间步长内通过构筑物界面的水量，然后采用源项处理的方法在与该界面相邻的上下游两个单元中分别扣除和增加这部分水量，为保证计算的连续性，采用无反射边界条件对构筑物界面处的水力要素值进行重构，进而计算构筑物界面处的数值通量。本发明可以准确反映各构筑物的过流特点，同时保证一维河道水流模型整体计算的精度和稳定性，为一维有限体积法数值框架下的构筑物过流处理提供了一种新的解决途径。

Description

一种基于有限体积法的河道构筑物过流数值模拟方法

技术领域

本发明涉及水利工程领域，尤其涉及防洪风险分析领域，具体为一种基于有限体积法的河道构筑物过流数值模拟方法。

背景技术

所谓河道构筑物指的是河道内存在的挡水堰、拦河闸、桥涵等对河道水流运动存在明显影响的构筑物。在工程实践中，通常采用一维数学模型来模拟河道洪水水流运动过程，而构筑物处的水流流态极度复杂，一维模型所采用的控制方程组已不能描述这种复杂水流流态，构筑物也被称为计算区域内的奇异点，需要专门特殊处理。工程实践中，河道中的构筑物是非常常见的，因此，在一维模型中能否合理模拟河道中的构筑物过流过程，则决定了模型是否能够大范围的推广应用。

在现有的主流洪水一维数值模拟计算中，由于隐式有限差分法发展较早，起源于上世纪六十年代，是现在河道洪水模拟的主流格式，常见的有限差分格式有Preissmann格式、Abbott-Lonescu格式等，在这类模型中处理构筑物的时候，是以添加虚拟断面和求解大型稀疏矩阵的方式完成的。但是这类基于隐式有限差分的数值格式不具备激波捕捉和处理流态过度的能力，在一些河道地形陡峭，坡度变化较大的河道进行应用时，常常造成数值解的发散，甚至模拟失败的问题。

近些年，以求解Riemann近似解为基础的Godunov有限体积法在洪水数值模拟领域得到了广泛应用，这类方法既可以很好的模拟坡度平缓的河道水流，同时能够自适应模拟流态的变化，在坡度陡峭的山丘区河道也能很好的适应，因此，该类方法逐渐发展成为模拟河道水流的新型计算格式，如何在新的有限体积法的数值框架下处理好河道构筑物的过流问题也是当前该领域的一个热点和难点问题。

发明内容

本发明的目的在于提供一种基于有限体积法的河道构筑物过流数值模拟方法，可以在一维显式有限体积法的数值框架下准确的反应各构筑物的过流特点，同时保证一维河道水流模型整体计算的精度和稳定性。

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

一种基于有限体积法的河道构筑物过流数值模拟方法，采用一维有限体积单元离散河道，将河道构筑物所在的位置作为河道计算单元的构筑物界面，根据构筑物过流公式(如堰流公式、闸口出流公式等)计算每一时间步长内通过该界面的水量，采用源项处理的方法在与该界面相邻的上下游两个单元中分别扣除和增加这部分水量，为保证计算的连续性，对构筑物界面处的水力要素值进行数值重构。包括以下具体步骤：

1)获取河道和构筑物基础数据：获取河道平面几何形状和控制断面形状数据，相邻断面间距不大于1km，获取河道内构筑物的空间位置信息以及几何尺寸信息，为保证计算精度，在构筑物上下游各布设两根控制断面，断面间距不大于100m。

2)对河道进行有限体积离散：采用一维有限体积单元离散河道，断面位置为单元中心位置，两断面中点位置为单元界面位置，如该位置无构筑物，该界面称为常规界面，如该位置有构筑物，则该界面称为构筑物界面，河道水力要素值存储在单元中心。

3)一维河道模型初始化：计算条件初始化，给河道各单元赋值初始水力要素值，即初始水位和初始流量值，设置构筑物初始运行状态，获取河道上下游边界的初始值。

4)获取t时刻外边界条件，根据CFL(Courant-Friedrichs-Lewy)条件，获取计算时间步长dt。

5)求解各常规界面的数值通量、构筑物的过流量和构筑物界面的数值通量：采用考虑源项的一维守恒型浅水方程组描述河道水流运动，采用基于HLL近似Riemann解的有限体积法计算t时刻各单元常规界面处的数值通量，根据构筑物过流特点(过流特点如：有的构筑物出流仅受上游水位控制、有的受上下游水位同时控制，通过构筑物的水流有的是明渠流，有的是压力流等)，采用构筑物上下游单元水力要素值计算t时刻通过该构筑物的流量值，该流量值在上下游单元中作为连续方程的源项处理，为保证计算的连续性，该构筑物界面处的水力要素值采用无反射边界条件进行重构，进而计算通过构筑物界面处的数值通量。

6)河道求解，获得t+dt时刻各单元水力要素值：通过t时刻各单元界面处的数值通量值和构筑物过流的源项值更新t+dt时刻各单元中心的水力要素值，河道两端的边界条件同步更新到t+dt时刻。

7)令t＝t+dt，重复步骤4)～6)，直到计算结束。

进一步的，步骤4)dt的选取受到CFL(Courant-Friedrichs-Lewy)条件限制，具体如式(1)所示：

式中：N_cfl为CFL数，u为断面平均流速，c为波速，Δx为有限体积单元空间步长，dt为时间步长。

进一步的，步骤5)采用的控制方程组为考虑构筑物过流源项的一维守恒型浅水方程组，具体如式(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为河道单位长度内的源项值。

进一步的，步骤5)根据构筑物过流特点，采用构筑物上下游单元水力要素值计算t时刻通过该构筑物的过流量，以河道中的挡水堰构筑物为例，其流量值计算公式如式(3)所示：

式中，h_up＝Z_up-Z_weir；h_down＝Z_down-Z_weir；q_w为通过构筑物的流量；Z_up，Z_down分别为构筑物上下游的计算单元水位；Z_weir为堰顶高程，l为堰宽。

进一步的，步骤5)中为为保证所有界面通量计算的连续性，采用无反射边界对构筑物界面处的水力要素值进行重构，具体如下：当分别计算构筑物上游单元和下游单元通过此界面处的数值通量时，界面处的水力要素值分别重构为：

Z^*＝Z_up，Q^*＝0 (4)

Z^*＝Z_down，Q^*＝0 (5)

式中，Z^*,Q^*为单元界面处的水力要素值，Z^*,Q^*分别为单元界面处重构后的水位和流量值。Z_up为构筑物上游单元水位，Z_down为构筑物下游单元水位。

本发明的优点和有益效果是：

本发明的方法可以准确反映各构筑物的过流特点，同时保证一维河道水流模型整体计算的精度和稳定性，为一维有限体积法数值框架下的构筑物过流处理提供了一种新的解决途径。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1为本发明的一种基于有限体积法的河道构筑物过流数值模拟方法流程图；

图2为带构筑物的河道有限体积单元离散示意图；

图中，i-1、i和i+1为有限体积单元的编号，i-1/2和i+1/2为第i个单元的左右两个界面，Δx_i为第i个有限体积单元的长度，

和

为通过单元左右两个界面的数值通量；

图3为带挡水堰河道算例的数值计算结果。

具体实施方式

实施例一：

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

本发明提供的是一种基于有限体积法的河道构筑物过流数值模拟方法。该方法采用一维有限体积单元离散河道，将河道构筑物所在的位置作为河道计算单元的构筑物界面，根据构筑物过流公式计算每一时间步长内通过该界面的水量，采用源项处理的方法在与该界面相邻的上下游两个单元中分别扣除和增加这部分水量，为保证计算的连续性，对构筑物界面处的水力要素值进行数值重构。该方法可以在一维显式有限体积法的数值框架下准确的反应各构筑物的过流特点，同时保证一维河道水流模型整体计算的精度和稳定性。该方法包括如下具体步骤：

1)获取河道平面几何形状和控制断面形状数据，相邻断面间距不大于1km，获取河道内构筑物的空间位置信息以及几何尺寸信息，为保证计算精度，在构筑物上下游各布设两根控制断面，断面间距不大于100m。

2)采用一维有限体积单元离散河道，断面位置为单元中心位置，两断面中点位置为单元界面位置，如该位置无构筑物，该界面称为常规界面，如该位置有构筑物，则该界面称为构筑物界面，河道水力要素值存储在单元中心，如图2所示。

3)计算条件初始化，给河道各单元赋值初始水力要素值，即初始水位和初始流量值，设置构筑物初始运行状态，获取河道上下游边界的初始值。

4)根据CFL(Courant-Friedrichs-Lewy)条件，获取计算时间步长dt，CFL的具体条件如式(1)所示：

式中：u为断面平均流速，c为波速，Δx为有限体积单元空间步长，dt为时间步长。为保证整体数值计算稳定，N_cfl建议取为不超过1.0的数值。

5)采用考虑源项的一维守恒型浅水方程组描述河道水流运动，具体如式(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为河道单位长度内的源项值。

采用基于HLL近似Riemann解的有限体积法计算t时刻各单元常规界面处的数值通量，详细的HLL近似Riemann解求解过程可参见文献描述(张大伟,程晓陶,黄金池等.复杂明渠水流运动的高适用性数学模型[J].水利学报,2010,41(4):531-536；张大伟，基于Godunov格式的堤坝溃决水流数值模拟[M].中国水利水电出版社,北京,2014,12)。

根据构筑物过流特点，采用构筑物上下游单元水力要素值计算t时刻通过该构筑物的过流量，以河道中的挡水堰构筑物为例，其流量值计算公式如式(3)

所示：

式中，h_up＝Z_up-Z_weir；h_down＝Z_down-Z_weir；q_w为通过构筑物的流量；Z_up，Z_down分别为构筑物上下游的计算单元水位；Z_weir为堰顶高程，l为堰宽。

对于不同的构筑物都可以通过其上下游相邻单元的水力要素值获得当前时刻的过流量，需要注意的是有的构筑物有确定的过流公式可以直接采用，有的构筑物设计复杂，没有确定的过流公式可以采用，这种情况则需要首先通过物理试验等手段获取构筑物的过流曲线公式，再通过构筑物上下游的水力要素值计算获得相应的过流量。对确定的构筑物来讲，t时刻从构筑物上游单元流出的水量必然会进入到下游相邻的单元，在控制方程中采用添加源项的方式来考虑构筑物的过流问题，可以从理论上保证数学模型计算整体的质量守恒特性，同时该方法也能很好的保证计算模型的稳定性。

为保证所有界面通量计算的连续性，采用无反射边界对构筑物界面处的水力要素值进行重构，具体如下：当分别计算构筑物上游单元和下游单元通过此界面处的数值通量时，界面处的水力要素值分别重构为：

Z^*＝Z_up，Q^*＝0 (4)

Z^*＝Z_down，Q^*＝0 (5)

式中，Z^*,Q^*为单元界面处的水力要素值，Z_up为构筑物上游单元水位，Z_down为构筑物下游单元水位。

根据构筑物界面处重构出的水力要素值，进而计算通过构筑物界面处的数值通量。

6)通过t时刻各单元界面处的数值通量值和构筑物过流的源项值更新t+dt时刻各单元中心的水力要素值，河道两端的边界条件同步更新到t+dt时刻；

7)令t＝t+dt，重复步骤4)～6)，直到计算结束。

图3为一带挡水堰的河道恒定流算例计算结果，该河道长1000m，宽50m，平底矩形，上游入流Q＝100m³/s，下游恒定水位1.5m，在河道正中间500m处有一挡水堰，该挡水堰宽50m，高1.0m。河道平均离散成20个单元，每个单元长50m，河道糙率设定为0.025。由图3的计算结果可以看出，未考虑挡水堰和考虑挡水堰后，各单元计算出的流量值与理论值均完全符合，稳定在100m³/s，这也说明采用本发明处理构筑物的方法可以很好的保证数学模型的质量守恒特性；从水面线看，不考虑挡水堰时，上游端水位稳定在1.957m，考虑挡水堰的影响后，挡水堰上游有明显的雍水效果，此时上游水位稳定在1.998m，并且在水流通过挡水堰后，可以看到水面有一个明显的跌水效应，随着流路变长，挡水堰下游水位逐渐恢复到原始的状态，计算结果反映出的这些特点与真实挡水堰的过流效果是一致的。综上说明本发明提出的基于有限体积法的河道构筑物过流数值模拟方法是成功的。

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

Claims (6)

1.一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：将河道构筑物所在的位置作为河道计算单元的界面，根据构筑物过流公式计算每一时间步长内通过该界面的水量，采用源项处理的方法在与该界面相邻的上下游两个单元中分别扣除和增加这部分水量；具体步骤如下：

1)获取河道和构筑物基础数据：获取河道平面几何形状和控制断面形状数据，相邻断面间距不大于1km；获取河道构筑物的空间位置信息以及几何尺寸信息，在构筑物上下游各布设两根控制断面，且断面间距不大于100m；

2)对河道进行有限体积离散：采用一维有限体积单元离散河道，断面位置为单元中心位置，两断面中点位置为单元界面位置，如该位置无构筑物，该界面称为常规界面，如该位置有构筑物，则该界面称为构筑物界面，河道水力要素值存储在单元中心；

3)一维河道模型初始化：计算条件初始化，给河道各单元赋值初始水力要素值，即初始水位和初始流量值，设置构筑物初始运行状态，获取河道上下游边界条件的初始值；

4)获取t时刻外边界条件，根据CFL条件，获取计算时间步长dt；

5)求解各常规界面的数值通量、构筑物的过流量和构筑物界面的数值通量：采用考虑源项的一维守恒型浅水方程组描述河道水流运动，采用基于HLL近似Riemann解的有限体积法计算t时刻各单元常规界面处的数值通量；采用构筑物上下游单元水力要素值计算t时刻通过该构筑物的流量值，该流量值在上下游单元中作为连续方程的源项处理；构筑物界面处的水力要素值采用无反射边界条件进行重构，进而计算通过该构筑物界面处的数值通量；

6)河道求解，获得t+dt时刻各单元水力要素值：通过t时刻各单元界面处的数值通量值和构筑物过流的源项值更新t+dt时刻各单元中心的水力要素值，河道两端的边界条件同步更新到t+dt时刻；

7)令t＝t+dt，重复步骤4)～6)，直到计算结束。

2.根据权利要求1所述的一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：所述河道构筑物包括河道内存在的挡水堰、拦河闸、桥涵。

3.根据权利要求1所述的一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：步骤4)dt的选取受到CFL条件限制，具体如式(1)所示：

式中：N_cfl为CFL数，u为断面平均流速，c为波速，Δx为有限体积单元空间步长，dt为时间步长。

4.根据权利要求1所述的一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：步骤5)采用控制方程组为考虑构筑物过流源项的一维守恒型浅水方程组，具体如式(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为河道单位长度内的源项值。

5.根据权利要求1所述的一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：步骤5)根据构筑物过流特点，采用构筑物上下游单元水力要素值计算t时刻通过该构筑物的过流量，当河道构筑物为挡水堰时，其流量值计算公式如式(3)所示：

式中，h_up＝Z_up-Z_weir；h_down＝Z_down-Z_weir；q_w为通过构筑物的流量；Z_up，Z_down分别为构筑物上下游的计算单元水位；Z_weir为堰顶高程，l为堰宽。

6.根据权利要求5所述的一种基于有限体积法的河道构筑物过流数值模拟方法，其特征在于：步骤5)中采用无反射边界条件对构筑物界面处的水力要素值进行重构，具体如下：当分别计算构筑物上游单元和下游单元通过此界面处的数值通量时，界面处的水力要素值分别重构为：

Z^*＝Z_up，Q^*＝0 (4)

Z^*＝Z_down，Q^*＝0 (5)

式中，Z^*,Q^*分别为单元界面处重构后的水位和流量值。

Images