一种考虑汊点面积的通用河网水流数值模拟方法

摘要

本发明公开了一种考虑汊点面积的通用河网水流数值模拟方法。获取河网平面几何形状数据和各河段横断面数据，建立各河段和汊点的拓扑关系以及各汊点处的水位面积关系。河网采用有限体积单元离散，水力要素值存储在单元中心。通过建立河网汊点处水量平衡方程，求解汊点处的水位值，该水位值将为与该汊点相连接的河段提供边界条件，汊点提供的水位边界条件与各河段的外边界条件配合后，每条河段的首末端均具备了边界定解条件，再通过Godunov显格式有限体积法来分别模拟各河段水流运动过程。该方法既能处理平原区的缓流流态，也能很好处理山丘区河网的复杂流态和大梯度水面的模拟问题，具有很好的通用性，能够有效弥补现有平原河网水流模拟方法的不足。

权利要求

1.一种考虑汊点面积的通用河网水流数值模拟方法，其特征在于：通过建立河网汊点处水量平衡方程，求解汊点处的水位值，为与该汊点相连接的各条河段提供边界条件，将该边界条件与各河段外边界条件配合，通过Godunov显格式有限体积法来模拟各河段水流运动过程，包括以下具体步骤：

1)获取河网平面几何形状和各河段横断面数据，各横断面间的间距不大于1km；

2)对各河段和汊点进行编号，建立汊点和各河段的对应关系；采用有限体积单元离散各河段，水力变量存储在各单元中心点；根据河网平面几何形状和首末端单元的断面形状建立各汊点处的水位面积关系Z-A；

3)模型初始化：给各河段一维河道单元设置糙率和水力条件初始值；根据各河段首末单元的初始水位获取汊点处初始水位；

4)确定计算时间步长dt：一维河道模型采用显格式的有限体积法求解，该时间步长获取受到CFL条件限制；获取当前时刻各河段外边界条件，当前时刻以t时刻表示；

5)采用完整的一维Saint-Venant方程组描述河道水流运动，采用HLL显格式计算通过单元界面处的数值通量，通过t时刻各单元的水力初始值、外边界条件以及汊点初始水位计算获得各河段t+dt时刻各单元的水力要素值；

6)建立汊点处水量平衡方程，并根据各汊点水量平衡方程获得t+dt时刻各汊点的水位值；采用如公式(3)的各汊点水量平衡方程计算t+dt时刻各汊点的水位值：式中： 为第k个汊点在t时刻的水面面积， 为第k个汊点在t时刻的水位， 为需要求解的第k个汊点在t+dt时刻的水位，n为与第k个汊点相连接的河段个数， 为与汊点k相连接的第i个河段在t时刻流入/流出汊点的流量，流入为正，流出为负；

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

2.根据权利要求1所述的一种考虑汊点面积的通用河网水流数值模拟方法，其特征在于：步骤2)在山丘区大坡度河网计算中，将天然汊点的范围向与其相连的各河段延伸，且不超过各河段首末端单元长度的1/2。

3.根据权利要求1所述的一种考虑汊点面积的通用河网水流数值模拟方法，其特征在于：步骤4)dt的选取受到CFL条件限制，具体如式(1)所示：式中：u为流速，c为波速，Δx为单元空间步长，dt为时间步长。

4.根据权利要求1所述的一种考虑汊点面积的通用河网水流数值模拟方法，其特征在于：步骤5)采用完整的一维Saint-Venant方程组为守恒形式，具体如式(2)所示：其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数。

说明书

一种考虑汊点面积的通用河网水流数值模拟方法

技术领域

[0001] 本发明涉及水利工程领域，尤其涉及防洪减灾领域，具体为一种考虑汊点面积的通用河网水流数值模拟方法。

背景技术

[0002] 河网是指由一系列河段首尾相互连接形成的具有水力联系的统一整体，河段和河段之间的连接点定义为汊点。在水利工程领域，河网的问题比较常见，如自然流域水系、城市的水系和大型灌区的灌排系统等。从上世纪五十年代至今，很多学者对河网水流数值模拟问题进行过研究，该问题一直是水利工程界和学术界关注的难点和热点问题。

[0003] 经过多年的发展，平原河网水流数值模拟方法已经发展的相对比较成熟，在这些方法中，隐式的有限差分方法应用最为广泛。早期的隐式河网求解技术需要将所有河段的断面水力要素同时求解，非常耗时，效率低下。后来在直接求解法的基础上逐渐发展出三级解法、四级解法和汊点分组解法等，将隐式河网的求解效率大幅度提升，河网求解技术得以在工程上大规模使用成为可能。近些年，我国政府在全国范围内开展了山洪灾害防治项目和中小河流治理项目，山洪沟和中小河流的坡度一般都很大，水的流态复杂，缓流和急流常同时出现，将传统的平原河网模型应用到这些区域时，往往会遭到失败。如应用广泛的经典Preissmann的四点格式，不能模拟跨临界流动，用于山洪沟或中小河流地区时，数值解常常发散。

[0004] 近些年，以求解Riemann近似解为基础的Godunov格式在浅水数值模拟领域得到了广泛应用，该格式既能模拟光滑的古典解，又能较好的模拟大梯度的水面流动，可以自动模拟流态过度、捕捉激波。目前，Godunov格式较多用于单河道的水流数值模拟，为了突破传统平原河网数学模型处理复杂流态时遇到的困难，将该格式推广到河网水流模拟领域具有重要的学术意义和工程应用意义，使其既能模拟平原河网水流，又能模拟山区河网水流，具有更广泛的通用性。

发明内容

[0005] 本发明的目的在于提供一种通用的河网水流数值模拟方法，既可以模拟缓流流态，也可以自动模拟各种流态过度，既能用于平原河网区，也能用于大坡度的山丘区河网。

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

[0007] 一种考虑汊点面积的通用河网水流数值模拟方法，通过建立河网汊点处水量平衡方程，求解汊点处的水位值，为与该汊点相连接的各条河段提供边界条件，将该边界条件与各河段外边界条件配合，通过Godunov显格式有限体积法来模拟各河段水流运动过程，从而达到模拟各种流态下河网水流运动过程的目的；包括以下具体步骤：

[0008] 1)获取河网平面几何形状和各河段横断面数据，各横断面间的间距不大于1km；

[0009] 2)对各河段和汊点进行编号，建立汊点和各河段的对应关系；采用有限体积单元离散各河段，水力变量存储在各单元中心点；根据河网平面几何形状和首末端单元的断面形状建立各汊点处的水位面积关系Z-A；

[0010] 3)模型初始化：给各河段一维河道单元设置糙率和水力条件初始值；根据各河段首末单元的初始水位获取汊点处初始水位；

[0011] 4)确定计算时间步长dt：一维河道模型采用显格式的有限体积法求解，该时间步长获取受到CFL(Courant-Friedrichs-Lewy)条件限制；获取当前时刻(t时刻)各河段外边界条件；

[0012] 5)采用完整的一维Saint-Venant方程组描述河道水流运动，采用HLL显格式计算通过单元界面处的数值通量，通过t时刻各单元的水力初始值、外边界条件以及汊点初始水位计算获得各河段t+dt时刻各单元的水力要素值；

[0013] 6)建立汊点处水量平衡方程，根据各汊点水量平衡方程获得t+dt时刻各汊点的水位值；

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

[0015] 进一步的，步骤2)建立的各汊点处的水位面积关系Z-A是根据河网平面几何形状和首末端单元的断面形状计算获取的。在山丘区大坡度河网计算中，为保证数值计算的稳定性，可将天然汊点的范围向与其相连的各河段适当延伸，以不超过各河段首末端单元长度的1/2为宜。

[0016] 进一步的，步骤4)确定计算时间步长dt，由于一维河道模型采用的是基于有限体积法的显式格式，因此，dt的选取受到CFL条件限制，具体如式(1)所示：

[0017]

[0018] 其中u为流速，c为波速，Δx为单元空间步长，dt为时间步长。为保证整体数值计算稳定，Ncfl建议取为不超过1.0的数值。

[0019] 进一步的，步骤5)采用完整的一维Saint-Venant方程组为守恒形式，具体如式(2)所示：

[0020]

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

[0022]

[0023] 式中B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，为便于随后的表述，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manni ng糙率系数。

[0024] 进一步的，步骤6)采用如公式(3)的各汊点水量平衡方程计算t+dt时刻各汊点的水位值：

[0025]

[0026] 式中 为第k个汊点在t时刻的水面面积， 为第k个汊点在t时刻的水位，为需要求解的第k个汊点在t+dt时刻的水位，n为与第k个汊点相连接的河段个数， 为与汊点k相连接的第i个河段在t时刻流入(流出)汊点的流量，流入为正，流出为负。

附图说明

[0027] 图1为本发明的一种考虑汊点面积的通用河网水流数值模拟方法流程图；

[0028] 图2为一维有限体积法示意图，

[0029] i-1、i和i+1为有限体积单元的编号，i-1/2和i+1/2为第i个单元的左右两个界面，Δxi为第i个有限体积单元的长度， 和 为通过单元左右两个界面的数值通量；

[0030] 图3为河网汊点示意图，

[0031] 图中，Q1和Q2为流入汊点的流量值，Q3为流出汊点的流量值；

[0032] 图4单河道大梯度水面潮波算例计算结果；

[0033] 图5为河网算例示意图；

[0034] 图6为河网算例的数值计算结果。

具体实施方式

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

[0036] 本发明提供的是一种考虑汊点面积的通用河网水流数值模拟方法。通过建立汊点的水位面积关系Z-A，由各河段流入(流出)该汊点的水量平衡方程推导下一时刻汊点的水位值，该水位值可作为各河段的边界条件；在单河段的水流数值求解时，采用基于Godunov格式的有限体积法进行，该模型既可以模拟缓流流态，也可以自动模拟各种复杂流态的过度，具有很好的通用性。通过这种方法建立的河网模型能够处理各种河网形状，各种水流流态，有效弥补现有平原河网数值求解方法的不足。该方法包括以下具体步骤：

[0037] 1)获取河网平面几何形状和各河段横断面数据，各横断面间的间距不大于1km；

[0038] 2)对各河段和汊点进行编号，建立汊点和各河段的对应关系；采用有限体积单元离散各河段，水力变量存储在各单元中心点；根据河网平面几何形状和首末端单元的断面形状计算获取各汊点处的水位面积关系Z-A，为保证大坡度山丘区河网水流数值计算的稳定性，可将天然汊点的范围向与其相连的各河段适当延伸，以不超过各河段首末端单元长度的1/2为宜。

[0039] 3)模型初始化：给各河段一维河道单元设置糙率和水力条件初始值；根据各河段首末单元的初始水位获取汊点处初始水位；

[0040] 4)确定计算时间步长dt：一维河道模型采用显格式的有限体积法求解，该时间步长获取受到CFL(Courant-Friedrichs-Lewy)条件限制，具体式(1)所示：

[0041]

[0042] 其中u为流速，c为波速，Δx为单元空间步长，dt为时间步长。为保证整体数值计算稳定，Ncfl建议取为不超过1.0的数值。

[0043] 获取当前时刻(t时刻)各河段外边界条件，可通过线性插值的方法，根据当前时刻值从外边界时间序列中插值出当前时刻的外边界值。上游端的外边界一般是流量边界，下边界是水位边界或水位、流量关系边界。

[0044] 5)采用完整的一维Saint-Venant方程组描述河道水流运动，采用有限体积法来离散一维Saint-Venant方程组，水力变量存储在单元中心，如图2所示。在有限体积法求解中，采用HLL显格式计算通过单元界面处的数值通量，通过t时刻各单元的水力初始值、外边界条件以及汊点初始水位计算获得各河段t+dt时刻各单元的水力要素值；采用的完整一维Saint-Venant方程组为守恒形式，具体如式(2)所示：

[0045]

[0046] 其中

[0047] 式中B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，为便于随后的表述，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，x为空间变量，t为时间变量，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数。

[0048] 一维河道模型的具体求解过程可参见文献描述(张大伟,程晓陶,黄金池等.复杂明渠水流运动的高适用性数学模型[J].水利学报,2010,41(4):531-536)。

[0049] 6)根据各汊点水量平衡方程获得t+dt时刻各汊点的水位值，汊点示意图如图3所示；采用如公式(3)的水量平衡公式计算t+dt时刻汊点的水位值：

[0050]

[0051] 式中 为第k个汊点在t时刻的水面面积， 为第k个汊点在t时刻的水位，为需要求解的第k个汊点在t+dt时刻的水位，n为与第k个汊点相连接的河段个数， 为与汊点k相连接的第i个河段在t时刻流入(流出)汊点的流量，流入为正，流出为负。

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

[0053] 图4为一单河道大梯度水面潮波算例的计算结果，该算例由Bermudez和Vazquez提出(Bermudez A,Vazquez M E.Upwind methods for hyperbolic conservation laws with source terms.Computers&Fluid,1994,23(8):1049-1071.)，主要是用来检验模型对大梯度潮波的计算能力。计算区域为一矩形渠道，长度L为14000m，底高程如下式所示：

[0054]

[0055] 计算时，空间步长取为280m，计算区域共离散成50个计算单元，计算CFL数取0.8，计算的初始水位为60.5m，流速为0m/s。上边界给定水位条件，其表达式为：

[0056]

[0057] 下边界给定流量边界，其流量值始终为0.0。

[0058] 该算例的解析解为：

[0059]

[0060] 图4为7552.13s时潮波的计算结果，由该图可以看出，流速的计算值与解析解均合良好，说明该模型可以很好的用于大梯度水面的潮波计算。

[0061] 图5为河网算例的示意图，该算例是由三条河段组成的一个典型河网算例，最早由Aral et al等提出(Aral M M,Zhang Y,Jin S.Application of relaxation scheme to wave-propagation simulation in open-channel networks[J].Journal of Hydraulic Engineering,1998,124(11):1125-1133)，用以检验模型处理河网水流问题的能力。在该算例中，三条河段长均为5000m，断面均为矩形，其中河段1和河段2的宽度为50m，河段3的宽度为100m，三条河段的底坡均为0.0002，Manning糙率系数均为0.025。各河段计算的初始水深均为1.43m，河段1和河段2的初始流量为50m3/s,河段3的初始流量为100m3/s。河段1和河段2的上游流量边界条件如图3所示，下游边界给定恒定的水位值1.43m。各河段的空间离散步长为100m，Ncfl数取0.8。图6为本文计算结果与Aral et al(1998)采用的松弛格式的计算结果对比。由图6计算结果可以看出，本文的计算结果与文献计算结果总体吻合较好，计算值洪峰165.91m3/s，文献值洪峰163.89m3/s。说明本发明提供的方法在处理河网问题上是成功的。

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