本发明公开了一种使用汊点影响区水面梯度的河网水流数值模拟方法。在本方法中，采用的是基于Godunov有限体积格式的显式数值求解方法，水力要素变量存储在单元中心。考虑在河网汊点影响区域内，水面梯度项是水流运动的最主要的驱动力，其它水流驱动项在该区域内可以忽略，所以在该区域内，利用汊点与各连接河段的水面梯度计算各河段流入或流出汊点的流量值，根据水量守恒原理建立汊点方程，通过数值迭代方法求解该汊点方程获得下一时刻的汊点水位值，该水位值将作为与该汊点相连接河段的水位边界条件，进而实现河网模型的求解。该方法提出的河网汊点处理方法具有清晰的物理机制，计算简洁高效，为显格式的有限体积河网求解方法提供了一种新的解决途径。
1.一种使用汊点影响区水面梯度的河网水流数值模拟方法，其特征在于：将汊点与各连接河段首/末端单元中心所连成的区域定义为汊点影响区域，在汊点影响区内利用汊点与各连接河段的水面梯度计算各河段流入或流出汊点的流量值，根据这些流量值建立汊点水量平衡方程，通过数值迭代方法求解该水量平衡方程获得下一时刻的汊点水位值，进而实现河网模型的求解，包括以下具体步骤：
1)测量或收集河网数据，包括河网平面几何形状以及各组成河段横断面信息；
2)对河网汊点及各河段进行编号，建立各汊点与各河段的拓扑关系，采用有限体积单元对各河段进行离散，对各河段离散后的有限体积单元依次进行编号，默认各河段水流的方向为从编号小的单元向编号大的单元流动，水力要素变量存储在有限体积单元的中心；
3)给各河段有限体积单元设置单元糙率，给各河段的有限体积单元设置初始水位值和初始流量值，给各汊点设置初始水位值；
4)当前时刻以T时刻表示，采用线性插值方法获取当前时刻各河段的外边界条件值；
5)根据CFL条件计算满足当前时刻计算稳定的时间步长DT；
6)采用完整的一维浅水方程组描述各河段水流运动，利用T时刻各有限体积单元中心的水力要素值、汊点水位值和外边界条件的值，采用HLL数值格式计算通过各河段有限体积单元界面处的数值通量，进而求得各河段有限体积单元在T+DT时刻的水力要素值；
7)求解T+DT时刻各汊点处的水位值：首先假定T+DT时刻各汊点处的水位值与T时刻汊点的水位值相等，即ZNODE(I，T+DT)＝ZNODE(I，T)，其中I为汊点的编号，ZNODE为汊点水位，计算该汊点与相邻各河段首/末单元的水面梯度，计算公式如下：SNODE(I,J)＝(ZLINK(I,J)-ZNODE(I,T+DT))/LENGTH(I,J) (1)其中SNODE(I,J)表示汊点I与与其相连接的第J条河段在汊点影响区的水面梯度，ZLINK(I,J)为与汊点I相连接的第J条河段在该汊点影响区内的计算单元在T+DT时刻的水位值，LENGTH(I,J)为汊点I与与其相连接的第J条河段的首/末端单元中心的距离；
根据计算出的汊点与各连接河段在汊点影响区的水面梯度值，计算各连接河段流入或流出汊点的流量，公式如下：QNODE(I,J)＝n(I,J)-1*R(I,J)2/3*SNODE(I,J)1/2*ANODE(I,J) (2)其中，QNODE(I,J)为与汊点I相连接的第J条河段流入或流出汊点I的流量值，如果该值为正，则说明是流入汊点，反之，则是流出汊点；n(I,J)为与汊点I相连接的第J条河段在汊点影响区的糙率值；R(I,J)为汊点I与与其相连的第J条河段的连接断面处的水力半径值；
ANODE(I,J)为汊点I与与其相连的第J条河段的连接断面处的过水面积；
处理河网汊点时，忽略汊点自身的水量变化，认为在DT时段内，流入汊点的流量和与流出汊点的流量和相等，如下式所示：其中，K为与汊点I相连接的河段总数；
将(1)式和(2)式带入(3)式整理后，得到汊点方程，如下式所示：
对式(4)采用Newton-Raphson迭代法进行求解，获得满足汊点水量守恒条件的T+DT时刻的汊点水位值ZNODE(I,T+DT),该值作为下一时刻求解时与汊点I相连接的各河段的水位边界条件；
8)令T＝T+DT，重复步骤4)～8)，直到计算结束。
2.根据权利要求1所述的一种使用汊点影响区水面梯度的河网水流数值模拟方法，其特征在于：步骤5)DT的选取受到CFL条件限制，具体如式(5)所示：式中：u为流速，c为波速，Δx为单元空间步长，DT为时间步长。
3.根据权利要求1所述的一种使用汊点影响区水面梯度的河网水流数值模拟方法，其特征在于：步骤6)采用完整的一维浅水方程组为守恒形式，具体如式(6)所示：其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数。
4.根据权利要求1所述的一种使用汊点影响区水面梯度的河网水流数值模拟方法，其特征在于：步骤7)中R(I,J)与ANODE(I,J)求解时，使用与汊点I相连接的第J条河段的首/末单元中心处的断面来作为汊点I与第J条河段的连接断面。
一种使用汊点影响区水面梯度的河网水流数值模拟方法
技术领域
[0001] 本发明涉及水利工程领域，尤其涉及防洪减灾领域，具体为一种使用汊点影响区水面梯度的河网水流数值模拟方法。
背景技术
[0002] 河网指的是一系列河段首尾相互连接形成的具有水力联系的统一整体，河段和河段之间的连接点定义为汊点。我国疆域辽阔，水系纵横交织，河网的问题非常常见。在进行流域洪水调度、农田水利规划、城市水系综合整治等方面，都离不开河网水流数值模拟技术的支持。学术界对河网水流数值模拟的研究起步较早，可追溯到上世纪五十年代，但是，随着水利行业的不断发展，现有的河网模拟技术已不能完全满足水利工程领域的需求。
[0003] 目前，最主流的河网数值模拟方法是基于隐式有限差分的方法，比较代表性的数值计算格式有Preissmann的四点中心隐格式和Abbott六点中心隐格式，各大洪水分析的商业软件也大都采用了这类方法。在河网早期的计算技术中，以直接求解大型稀疏矩阵为主，但是当遇到较大规模的河网时，计算速度将变得非常慢，计算效率低下，为克服这种情况，后来逐渐发展出河网的各类分级解法，其本质都是降低稀疏矩阵的求解规模，加快计算的速度。目前工程实践中采用的河网数值模拟方法大多是这类基于分级思想的隐式河网求解技术。这类方法在应用到水流变化比较平稳的平原地区时，效果较好，但是当遇到水面梯度较大的情形，由于这类数值格式本身不具备激波捕捉的功能，因此，常常会遇到数值结果发散的情形。我们国家当前正开展的山洪灾害防治项目和中小河流的治理项目中，常需用到河网水流数值模拟技术，这些地区的沟道坡度较大，水流流态变化比较剧烈，传统的河网技术已渐渐不再适用了，亟需发展新的适用性更好的河网水流模拟技术。
[0004] 近些年，以求解Riemann近似解为基础的Godunov有限体积格式得到了极大的发展，经典的数值格式有HLL格式，Roe格式等。这类数值格式本身具有激波捕捉能力，能够自动模拟急流和缓流的变化，适应大梯度的水面求解。当前，单河道的Godunov格式数值计算方法发展的已经相对较为成熟，如何将成熟的单河道Godunov格式推广应用到河网求解中是当前的一个研究热点问题。为了保证该类格式的激波捕捉特性，这种类型的数值格式基本上都是显格式的，所以原有的发展的较为成熟的汊点处理技术不再适用。
发明内容
[0005] 本发明的目的在于提供一种适用的汊点求解处理方法，将Godunov数值求解技术应用到河网水流数值模拟中，克服传统的河网求解数值格式不具备激波捕捉功能的缺陷，扩大河网水流数值模拟的适用范围。
[0006] 本发明是通过以下技术方案实现的：
[0007] 一种使用汊点影响区水面梯度的河网水流数值模拟方法，首先将河网汊点附近的区域定义为汊点影响区域，在该区域内，水面梯度项是水流运动的最主要的驱动力，其它驱动项在该区域内可以忽略，所以可以利用汊点与各连接河段的水面梯度计算各河段流入或流出汊点的流量值，根据这些流量值建立汊点质量平衡方程，通过数值迭代方法求解该质量平衡方程获得下一时刻的汊点水位值，进而实现河网模型的求解，包括以下具体步骤：
[0008] 1)测量或收集河网数据，包括河网平面几何形状以及各组成河段的详细横断面信息，一般讲，各横断面间距不宜大于1km，河道断面形状变化剧烈的地方需要适当增加断面数量；
[0009] 2)对河网汊点及各河段进行编号，建立各汊点与各河段的拓扑关系，采用有限体积单元对各河段进行离散，对各河段离散后的有限体积单元依次进行编号，默认各河段水流的方向为从编号小的单元向编号大的单元流动，水力要素变量存储在有限体积单元的中心；
[0010] 3)给各河段有限体积单元设置单元糙率，给各河段的有限体积单元设置初始水位值和初始流量值，给各汊点设置初始水位值；
[0011] 4)当前时刻以T时刻表示，采用线性插值方法获取当前时刻各河段的外边界条件值；
[0012] 5)根据CFL条件计算满足当前时刻计算稳定的时间步长DT；
[0013] 6)采用完整的一维浅水方程组描述各河段水流运动，利用T时刻各有限体积单元中心的水力要素值、汊点水位值和外边界条件的值，采用HLL数值格式计算通过各河段有限体积单元界面处的数值通量，进而求得各河段有限体积单元在T+DT的水力要素值；
[0014] 7)求解T+DT时刻各汊点处的水位值：首先假定T+DT时刻各汊点处的水位值与T时刻汊点的水位值相等，即ZNODE(I，T+DT)＝ZNODE(I，T)，其中I为汊点的编号，ZNODE为汊点水位，计算该汊点与与其相邻各河段首末单元的水位梯度，计算公式如下：
[0015] SNODE(I,J)＝(ZLINK(I,J)-ZNODE(I,T+DT))/LENGTH(I,J) (1)
[0016] 其中SNODE(I,J)表示汊点I与与其相连接的第J条河段在汊点影响区的水面梯度，ZLINK(I,J)为与汊点I相连接的第J条河段在该汊点影响区内的计算单元在T+DT时刻的水位值，LENGTH(I,J)为汊点I与与其相连接的第J条河段的首(末)端单元中心的距离。
[0017] 根据计算出的汊点与各连接河段在汊点影响区的水面梯度值，计算各连接河段流入或流出汊点的流量，公式如下：
[0018] QNODE(I,J)＝n(I,J)-1*R(I,J)2/3*SNODE(I,J)1/2*ANODE(I,J) (2)[0019] 其中，QNODE(I,J)为与汊点I相连接的第J条河段流入或流出汊点I的流量值，如果该值为正，则说明是流入汊点，反之，则是流出汊点；n(I,J)为与汊点I相连接的第J条河段在汊点影响区的糙率值；R(I,J)为汊点I与与其相连的第J条河段的连接断面处的水力半径值；ANODE(I,J)为汊点I与与其相连的第J条河段的连接断面处的过水面积。
[0020] 不考虑汊点自身的水量变化，根据汊点质量守恒原理，可认为在DT时段内，流入汊点的流量和与流出汊点的流量和相等，如下式所示：
[0021]
[0022] 其中，K为与汊点I相连接的河段总数。
[0023] 将(1)式和(2)式带入(3)式整理后，可得到以汊点水位ZNODE(I,T+DT)为自变量的一个非线性函数，该函数即为汊点方程，如下式所示：
[0024] f(ZNODE(I,T+DT))＝0 (4)
[0025] 对式(4)采用Newton-Raphson迭代法进行求解，获得满足汊点质量守恒条件的T+DT时刻的汊点水位值ZNODE(I,T+DT),该值作为下一时刻求解时与汊点I相连接的各河段的水位边界条件。
[0026] 8)令T＝T+DT，重复步骤4)～8)，直到计算结束。
[0027] 进一步的，步骤5)DT的选取受到CFL(Courant-Friedrichs-Lewy)条件限制，具体如式(5)所示：
[0028]
[0029] 式中：u为流速，c为波速，Δx为单元空间步长，DT为时间步长。
[0030] 进一步的，步骤6)采用完整的一维浅水方程组为守恒形式，具体如式(6)所示：
[0031]
[0032] 其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：
[0033]
[0034] 式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数。
[0035] 进一步的，步骤7)中在计算公式(2)里R(I,J)与ANODE(I,J)时，可使用与汊点I相连接的第J条河段的首(末)单元中心处的断面来作为汊点I与第J条河段的连接断面。
附图说明
[0036] 图1为本发明的一种使用汊点影响区水面梯度的河网水流数值模拟方法流程图；
[0037] 图2为河网汊点影响区示意图；
[0038] 图3为河网算例示意图；
[0039] 图4为河网算例的数值计算结果。
具体实施方式
[0040] 下面结合附图1和实施例对本发明作进一步说明。
[0041] 本发明提供的是一种使用汊点影响区水面梯度的河网水流数值模拟方法。根据汊点影响区内的水流主要是由水面梯度项驱动的特点，利用汊点与各连接河段的水面梯度计算各河段流入或流出汊点的流量值，根据这些流量值建立汊点质量平衡方程，通过数值迭代方法求解该质量平衡方程获得下一时刻的汊点水位值，进而实现河网模型的求解。该方法为将Godunov显格式的有限体积法推广到河网水流模拟中提供了一种可行的方法，与现有的显格式汊点处理方法相比，该方法物理意义清晰，数值处理方法简单，具有广阔的推广应用空间。该方法包括以下具体步骤：
[0042] 1)测量或收集河网数据，包括河网平面几何形状以及各组成河段的详细横断面信息，一般讲，各横断面间距不宜大于1km，河道断面形状变化剧烈的地方需要适当增加断面数量；
[0043] 2)对河网汊点及各河段进行编号，建立各汊点与各河段的拓扑关系，采用有限体积单元对各河段进行离散，对各河段离散后的有限体积单元依次进行编号，默认各河段水流的方向为从编号小的单元向编号大的单元流动，水力要素变量存储在有限体积单元的中心；
[0044] 3)给各河段有限体积单元设置单元糙率，给各河段的有限体积单元设置初始水位值和初始流量值，给各汊点设置初始水位值；
[0045] 4)当前时刻以T时刻表示，采用线性插值方法获取当前时刻各河段的外边界条件值；
[0046] 5)根据CFL条件计算满足当前时刻计算稳定的时间步长DT，DT的选取受到CFL(Courant-Friedrichs-Lewy)条件限制，具体如式(5)所示：
[0047]
[0048] 式中：u为流速，c为波速，Δx为单元空间步长，DT为时间步长。
[0049] 6)采用完整的一维浅水方程组描述各河段水流运动，采用完整的一维浅水方程组为守恒形式，具体如式(6)所示：
[0050]
[0051] 其中，x为空间变量，t为时间变量，D、U、F、S为方程组中各变量的向量表述，具体如下：
[0052]
[0053] 式中：B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数。
[0054] 利用T时刻各有限体积单元中心的水力要素值、汊点水位值和外边界条件的值，采用HLL数值格式计算通过各河段有限体积单元界面处的数值通量，进而求得各河段有限体积单元在T+DT的水力要素值。关于一维HLL格式有限体积模型的具体求解过程可参见文献描述(张大伟,程晓陶,黄金池等.复杂明渠水流运动的高适用性数学模型[J].水利学报,2010,41(4):531-536)。
[0055] 7)求解T+DT时刻各汊点处的水位值：首先假定T+DT时刻各汊点处的水位值与T时刻汊点的水位值相等，即ZNODE(I，T+DT)＝ZNODE(I，T)，其中I为汊点的编号，ZNODE为汊点水位，计算该汊点与与其相邻各河段首末单元的水位梯度，计算公式如下：
[0056] SNODE(I,J)＝(ZLINK(I,J)-ZNODE(I,T+DT))/LENGTH(I,J) (1)
[0057] 其中SNODE(I,J)表示汊点I与与其相连接的第J条河段在汊点影响区的水面梯度，ZLINK(I,J)为与汊点I相连接的第J条河段在该汊点影响区内的计算单元在T+DT时刻的水位值，LENGTH(I,J)为汊点I与与其相连接的第J条河段的首(末)端单元中心的距离。
[0058] 根据计算出的汊点与各连接河段在汊点影响区的水面梯度值，计算各连接河段流入或流出汊点的流量，公式如下：
[0059] QNODE(I,J)＝n(I,J)-1*R(I,J)2/3*SNODE(I,J)1/2*ANODE(I,J) (2)[0060] 其中，QNODE(I,J)为与汊点I相连接的第J条河段流入或流出汊点I的流量值，如果该值为正，则说明是流入汊点，反之，则是流出汊点；n(I,J)为与汊点I相连接的第J条河段在汊点影响区的糙率值；R(I,J)为汊点I与与其相连的第J条河段的连接断面处的水力半径值；ANODE(I,J)为汊点I与与其相连的第J条河段的连接断面处的过水面积；该连接断面可使用与汊点I相连接的第J条河段的首(末)单元中心处的断面代替。
[0061] 不考虑汊点自身的水量变化，根据汊点质量守恒原理，可认为在DT时段内，流入汊点的流量和与流出汊点的流量和相等，如下式所示：
[0062]
[0063] 其中，K为与汊点I相连接的河段总数。
[0064] 将(1)式和(2)式带入(3)式整理后，可得到以汊点水位ZNODE(I,T+DT)为自变量的一个非线性函数，该函数即为汊点方程，如下式所示：
[0065] f(ZNODE(I,T+DT))＝0 (4)
[0066] 对式(4)采用Newton-Raphson迭代法进行求解，获得满足汊点质量守恒条件的T+DT时刻的汊点水位值ZNODE(I,T+DT),该值作为下一时刻求解时与汊点I相连接的各河段的水位边界条件。关于Newton-Raphson迭代法的详细求解方法，可参见文献描述(张大伟.基于Godunov格式的堤坝溃决水流数值模拟[M].北京:中国水利水电出版社,2014:130-134.)[0067] 8)令T＝T+DT，重复步骤4)～8)，直到计算结束。
[0068] 图3为河网算例的示意图，该算例是由三条河段组成的一个典型河网算例，最早由Aral et al等提出(Aral M M,ZhangY,Jin S.Application ofrelaxation scheme to wave-propagation simulation in open-channel networks[J].Journal of Hydraulic Engineering,1998,124(11):1125-1133)，用以检验模型处理河网水流问题的能力。后来，张大伟等也曾使用该算例来检验所开发的数学模型结果的合理性(张大伟等,一种考虑汊点面积的通用河网水流数值模拟方法:中国,ZL201910122266.6[P],2019-9-27；张大伟等,应用Godunov格式模拟复杂河网明渠水流运动[J],应用基础与工程学报,2015,23(6):1088-1096.)。在该算例中，三条河段长均为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，CFL数取0.8。图4为本文计算结果与Aral et al(1998)采用的松弛格式的计算结果对比。由图4计算结果可以看出，采用本发明提出的河网水流数值模拟方法的计算结果与原文献的计算结果基本一致，过程吻合较好，计算值洪峰
162.83m3/s，文献值洪峰163.89m3/s。说明本发明提供的方法在处理河网问题上是成功的。
[0069] 上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。
