CN110263411A - 复杂河网水动力计算方法、设备、介质 - Google Patents

Abstract

本发明公开了复杂河网水动力计算方法、设备和介质，本发明方法首先为各汊点设置一个初始的预测水位值，然后针对河道进行分类；针对于缓坡河道和陡峭河道，分别采用有限差分法和有限体积法计算各断面的水动力差数；接着根据各汊点处各河道断面的流量，判断各汊点处各河道断面的流量是否满足流量守恒条件；在满足流量守恒情况下或不满足流量守恒但迭代次数达到最大值情况下，将上述计算的水动力参数作为水动力计算结果；在不满足流量守恒且迭代次数还未达最大值情况下，计算汊点水位校正值，在该值小于一定值的情况下，将将校正后的汊点处水位值作为汊点当前的预测水位值，然后进入下一次迭代计算；本发明方法具有计算速度快以及稳定性强的优点。

Description

复杂河网水动力计算方法、设备、介质

技术领域

本发明属于水动力数值模拟技术领域，特别涉及一种复杂河网水动力计算方法、计算设备和存储介质。

背景技术

河流、湖泊是人类生存和发展的重要资源，为人类提供了重要的生产和生活资源。随着我国社会和经济的快速发展，人口数量的增多，工业化和城市现代化进程的不断加快，社会用水量迅速增加，人们对这些水体的开发利用力度越来越大。河网水动力模型是描述河网水域中水利的数学模型。河网水动力计算是水利、航运和环保等部门经常要进行的工作，可对河网地区的防洪、排涝、灌溉、航运以及水污染等水火灾害防治和水环境保护提供可靠的解决方案。

河口区城市水系具有河网密布、河道纵横交错、闸泵工程数量多等特点，属于全要素复杂河网，即由山区河道、平原河网、湖库、堤防工程、闸泵和下游河口、外海组成的复杂河网系统。这类复杂河网系统对一维河网洪水模拟方法的稳定性要求极高。

Preissmann隐式差分格式与分级解法是最为经典的河网水动力模型求解方法，具有计算效率高的优点，目前仍得到了广泛应用。然而，这类传统方法难以适用于跨临界流和间断流态模拟问题，应用于山区陡峭河道或干河床的洪水模拟时常遇到计算不稳定、水位及流速异常等问题。此外，由于数值迭代计算误差累积效应，这类传统方法对初始水位及流量状态要求较高，即需要为不同量级流量边界计算方案配置相应合适的各断面初始水位及流量，否则容易导致计算崩溃。因此，应用于复杂河网系统洪水实时分析时，这类传统方法的稳定性问题较为突出，难以满足全要素复杂河网洪水模拟需求。

Godunov型有限体积法通过构造Riemann问题计算界面数值通量，可有效处理跨临界流和间断流态模拟问题，具有稳定性强的优点。然而，由于其显式计算格式的时间步长受克朗条件制约，有限体积法的计算效率要低于有限差分格式。

发明内容

本发明的第一目的在于克服现有技术的缺点与不足，提供一种计算速度快以及稳定性强的复杂河网水动力计算方法。

本发明的第二目的在于提供一种计算设备。

本发明的第三目的在于提供一种存储介质。

本发明的第一目的通过下述技术方案实现：一种复杂河网水动力计算方法，步骤如下：

步骤S1、针对河网中各汊点，对各汊点分别预设一个水位值作为其初始的预测水位值；

步骤S2、根据河道水面坡降将河道分类为缓坡河道和陡峭河道；

步骤S3、根据河网中各汊点当前的预测水位值，计算各河道各断面的水动力参数；其中：

针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数；

针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数；

步骤S4、针对于河网中各汊点，获取通过步骤S3计算得到的该汊点处各河道断面的流量，判断该汊点处各河道断面的流量是否满足流量守恒条件；

若是，则将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，则进入步骤S5；

步骤S5、首先判断当前迭代次数是否达到一定值；

若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，首先根据汊点处流入流量与流出流量的差值计算汊点水位校正值，然后判断该水位校正值是否小于一定值Z1；若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为河网水动力计算结果；若否，则根据汊点水位校正值对汊点当前的预测水位值进行校正，将校正后的汊点处水位值作为下一次迭代计算时汊点的预测水位值，然后返回步骤S3。

优选的，每一时段执行一次步骤S1至S5；即在当前时段执行完成步骤S1至S5后，进入下一时段重新执行步骤S1至S5，直到最终时段；

步骤S2中根据河道水面坡降将河道分类为缓坡河道和陡峭河道的具体过程如下：

步骤S21、获取上一时段的河网水动力计算结果，得到上一时段各河道各断面的水动力参数，包括各断面的水位；

步骤S22、针对于各河道，计算其中每相邻两个断面之间的水面坡降：

其中，S_i为河道内第i、i+1个断面之间的水面坡降，Z_i、Z_i+1分别为上一时段河道内第i、i+1个断面的水位；Δx_i为河道内第i、i+1个断面之间的距离；M为河道内断面的总数；

步骤S23、比较步骤S22计算得到的各每相邻两个断面之间的水面坡降，选取出其中一个最大值：

S_max＝max(|S₁|,...,|S_i|,...,|S_M-1|)；

其中，S_max为选取出的最大值；M为河道内断面的总数；

步骤S24、根据步骤S23中选取出的最大值S_max对河道进行分类，当S_max大于等于阈值Z2时，判断河道为陡峭河道，否则判断河道为缓坡河道。

优选的，每一时段执行一次步骤S1至S5；即在当前时段执行完成步骤S1至S5后，进入下一时段重新执行步骤S1至S5，直到最终时段；

当前时段在执行步骤S1中，将上一时段最后计算得到的汊点的预测水位值，作为当前时段各汊点初始的预测水位值。

优选的，步骤S4中通过判断以下公式是否成立来判断汊点处各河道断面的流量是否满足流量守恒条件：

其中，p为汊点编号，j为汊点处的河道编号；Q_p,j为第p个汊点处第j条河道的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量。

优选的，步骤S3中针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数具体过程如下：

步骤S31、对单一缓坡河道的水动力方程进行有限差分离散，得第一式如下：

式中：i为河道内断面编号；为第n+1步对应时刻的河道内第i个断面的水位、流量；a_1i、a_2i、c_1i、c_2i、d_2i、E_1i和E_2i为已知的差分系数；n为计算步数；M为河道内断面的总数；

采用追赶法对河道内各断面流量进行递推，得第二式：

式中：为第n+1步对应时刻河道首断面的水位、流量；为第n+1步对应时刻的河道末断面的水位、流量；E₁、E₂、R₁、R₂、T₁、T₂为已知的递推系数；

步骤S32、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道首断面的水位，根据河道后一个汊点当前的预测水位值得到河道末断面的水位；根据河道首断面的水位和河道末断面的水位，通过步骤S31中的第二式计算得到河道首、末断面的流量；

步骤S33、根据河道首断面的水位、河道末断面的水位、河道首断面的流量以及河道末断面的流量，通过步骤S31中的第一式计算出河道其他断面的水位和流量；同时，根据河道各断面的水位、流量以及断面地形计算出各断面的流速。

优选的，步骤S3中针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数的具体过程如下：

步骤SA、对单一陡峭河段的水动力方程进行有限体积离散，得第一式：

其中，下标i为河段内第i个断面对应的控制体编号；n为计算步数；为第n步对应时刻的第i个断面的状态向量；为系数矩阵；Δt为计算时间步长；Δx_i为第i个断面对应的控制体长度；分别为第n步对应时刻的第i个断面对应的控制体的左、右两侧界面的数值通量；

其中，为源项近似：

式中，为第n步对应时刻的第i个断面的侧向单宽流量；g为重力加速度；为第n步对应时刻的第i个断面的过水断面面积；为第n步对应时刻的第i个断面的水位；为第n步对应时刻的第i个断面的流量；为第n步对应时刻的第i个断面的水力半径；n′为曼宁糙率系数；

采用HLL格式计算界面数值通量，如第二式：

式中：为第n步对应时刻的第i个断面对应的控制体的左侧界面的数值通量，即第i-1与i断面之间的第n步对应时刻的界面数值通量；为第i-1与i断面之间的第n步对应时刻的界面的左、右波速；上标L、R分别代表界面左、右侧；和为基于界面左、右侧重构变量的数值通量，按下式计算：

由下式计算：

式中：上、下标L和R代表界面左、右侧；波速由下式计算：

式中：为第i断面第n步对应时刻的流速；为第i断面第n步对应时刻的特征速度；为第i断面第n步对应时刻的过水断面面积；为第i断面第n步对应时刻的过水断面水面宽度；

步骤SB、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道的上游水位边界条件，根据河道后一个汊点当前的预测水位值得到河道的下游水位边界条件；根据上述上、下游水位边界条件，通过步骤SA中的第二式计算得到河道首、末断面的边界界面处当前时刻即第n步对应时刻的数值通量；

步骤SC、获取上一时刻即第n-1步对应时刻所计算的河道各断面在当前时刻即第n步对应时刻的水位和流量；然后通过步骤SA中的第二式计算得到河道各断面当前时刻即第n步对应时刻各控制体之间界面的数值通量；基于当前时刻即第n步对应时刻各控制体之间界面的数值通量，通过然后通过步骤SA中的第一式计算得到河道各断面在下一时刻即第n+1步对应时刻的水位和流量。

优选的，步骤S5中，根据汊点处流入流量与流出流量的差值计算汊点水位校正值为：

其中p为汊点编号，j为汊点处的河道编号，m为当前迭代次数，即校正步数；Q_p,j,m为第p个汊点处第j条河道在第m次迭代计算时得到的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量；ΔZ_p,m为第m次迭代计算时根据汊点处流入流量与流出流量的差值计算汊点水位校正值；

δ_p,m为第m次迭代计算时第p个汊点的校正系数，具体为：

步骤S5中，根据汊点水位校正值对汊点当前的预测水位值进行校正，得到校正后的汊点处水位值为：

为汊点当前的预测水位值，为校正后的汊点处水位值；其中：

为下一次迭代计算时汊点的预测水位值。

更进一步的，Z1的值为10^-4米。

本发明的第二目的通过下述技术方案实现：一种计算设备，包括处理器以及用于存储处理器可执行程序的存储器，其特征在于，所述处理器执行存储器存储的程序时，实现本发明第一目的所述的复杂河网水动力计算方法。

本发明的第三目的通过下述技术方案实现：一种存储介质，存储有程序，其特征在于，所述程序被处理器执行时，实现本发明第一目的所述的复杂河网水动力计算方法。

本发明相对于现有技术具有如下的优点及效果：

本发明复杂河网水动力计算方法中，首先为各汊点分别设置一个初始的预测水位值，然后针对河道进行分类；针对于缓坡河道和陡峭河道，分别对应采用有限差分法和有限体积法计算两种河道各断面的水动力差数；接着根据各汊点处各河道断面的流量，判断各汊点处各河道断面的流量是否满足流量守恒条件；在满足流量守恒条件的情况下或不满足流量守恒但是当前迭代次数达到最大值的情况下，将上述计算的水动力参数作为河道最终水动力计算结果；在不满足流量守恒且迭代次数还没有达到最大值情况下，根据汊点处流入流量与流出流量的差值计算汊点水位校正值，在该值小于一定值的情况下，将将校正后的汊点处水位值作为汊点当前的预测水位值，然后进入下一次迭代计算。由上述可见，本发明计算方法基于流量守恒的条件对汊点水位进行校正，并且将复杂河网中不同分类的河道即缓坡河道和陡峭河道分别采用有限差分法和有限体积法计算河道各断面的水动力参数，实现了有限差分法与有限体积法的优化组合，显著提高了复杂河网水动力模型的适应性和实际应用价值，具有计算速度快以及稳定性强的优点，解决了传统方法无法充分利用有限差分法在平原河网计算中计算效率高的优势以及有限体积法在山区河道计算中处理复杂流态、干河床时稳定性强的特点等问题。

附图说明

图1是本发明方法流程图。

图2是含缓坡河道及陡峭河道的河网概化示意图。

图3是典型河道的深泓高程示意图。

图4是典型河道的水面线计算结果。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1

本实施例公开了一种复杂河网水动力计算方法，如图1所示，每一时段(每一个时间步长)执行一次下述步骤S1至S5；即在当前时段执行完成步骤S1至S5后，进入下一时段重新执行步骤S1至S5，直到最终时段；其中每一时段，本实施例步骤S1到步骤S5的执行过程如下：

步骤S1、针对河网中各汊点，对各汊点分别预设一个水位值作为其初始的预测水位值；在本实施例中，将上一时段最后计算得到的汊点的预测水位值，作为当前时段各汊点初始的预测水位值；若当前时段为初始时段，则针对各汊点任意设置一个水位值作为各汊点的初始的预测水位值。

步骤S2、根据河道水面坡降将河道分类为缓坡河道和陡峭河道；在本实施例中，河道分类的具体过程如下：

步骤S21、获取上一时段的河网水动力计算结果，得到上一时段各河道各断面的水动力参数，包括各断面的水位；

步骤S22、针对于各河道，计算其中每相邻两个断面之间的水面坡降：

其中，S_i为河道内第i、i+1个断面之间的水面坡降，Z_i、Z_i+1分别为上一时段河道内第i、i+1个断面的水位；Δx_i为河道内第i、i+1个断面之间的距离；

步骤S23、比较步骤S22计算得到的各每相邻两个断面之间的水面坡降，选取出其中一个最大值：

S_max＝max(|S₁|,...,|S_i|,...,|S_M-1|)；

其中，S_max为选取出的最大值；M为河道内的断面数量；

步骤S24、根据步骤S23中选取出的最大值S_max对河道进行分类，当S_max大于等于阈值时，判断河道为陡峭河道，否则判断河道为缓坡河道。在本实施例中上述阈值设置为0.01，即本步骤中，若S_max大于等于0.01，则判断河道为陡峭河道，若S_max小于0.01，则判断河道为缓坡河道。

本实施例中，在执行本实施例上述计算方法之前，针对河道每个断面分别设置初始的水位和流量信息；若当前执行计算方法时处于初始时段，则本实施例上述步骤S2在执行河道分类时，可以将上述设置的河道每个断面的初始水位信息作为上一时刻的河网水动力计算结果，然后执行步骤S21至步骤S24，实现河道分类。

步骤S3、获取各汊点当前的预测水位值，根据河网中各汊点当前的预测水位值，计算各河道各断面的水动力参数；其中：

(1)针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数；具体如下：

步骤S31、对单一缓坡河道的水动力方程进行有限差分离散，得第一式如下：

式中：i为河道内断面编号；为第n+1步对应时刻的河道内第i个断面的水位、流量；a_1i、a_2i、c_1i、c_2i、d_2i、E_1i和E_2i为已知的差分系数；n为计算步数；M为河道内断面的总数；

采用追赶法对河道内各断面流量进行递推，得第二式：

式中：为第n+1步对应时刻的河道首断面的水位、流量； 为第n+1步对应时刻的河道末断面的水位、流量；E₁、E₂、R₁、R₂、T₁、T₂为已知的递推系数；

步骤S32、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道首断面的水位，根据河道后一个汊点当前的预测水位值得到河道末断面的水位；根据河道首断面的水位和河道末断面的水位，通过步骤S31中的第二式计算得到河道首、末断面的流量；在本实施例中，可以将河道前一个汊点当前的预测水位值作为河道首断面的水位，河道后一个汊点当前的预测水位值作为河道末断面的水位。

步骤S33、根据河道首断面的水位、河道末断面的水位、河道首断面的流量以及河道末断面的流量，通过步骤S31中的第一式计算出河道其他断面的水位和流量；同时，根据河道各断面的水位、流量以及断面地形计算出各断面的流速：U＝Q/A。

(2)针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数，具体如下：

步骤SA、对单一陡峭河段的水动力方程进行有限体积离散，得第一式：

其中，下标i为河段内第i个断面对应的控制体编号；n为计算步数；为第n步对应时刻的第i个断面的状态向量；为系数矩阵；Δt为计算时间步长；Δx_i为第i个断面对应的控制体长度；分别为第n步对应时刻的第i个断面对应的控制体的左、右两侧界面的数值通量；在本实施例中，河道各断面的中点处定义为控制体的界面；河道相邻断面的界面之间定义为控制体，即控制体是以河道断面为中心、向左右相邻断面之间中点为边界的“微河段”，另外河道第一个和最后一个断面处也定义为界面，属于河道的边界界面。

其中，为源项近似：

式中，为第n步对应时刻的第i个断面的侧向单宽流量；g为重力加速度；为第n步对应时刻的第i个断面的过水断面面积；为第n步对应时刻的第i个断面的水位；为第n步对应时刻的第i个断面的流量；为第n步对应时刻的第i个断面的水力半径；n′为曼宁糙率系数；

采用HLL格式计算界面数值通量，如第二式：

式中：为第n步对应时刻的第i个断面对应的控制体的左侧界面的数值通量，即第i-1与i断面之间的第n步对应时刻的界面数值通量；为第i-1与i断面之间的第n步对应时刻的界面的左、右波速；上标L、R分别代表界面左、右侧；和为基于界面左、右侧重构变量的数值通量，按下式计算：

由下式计算：

式中：上、下标L和R代表界面左、右侧；波速由下式计算：

式中：为第i断面第n步对应时刻的流速；为第i断面第n步对应时刻的特征速度；为第i断面第n步对应时刻的过水断面面积；为第i断面第n步对应时刻的过水断面水面宽度；

步骤SB、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道的上游水位边界条件，根据河道后一个汊点当前的预测水位值得到河道的下游水位边界条件；根据上述上、下游水位边界条件，通过步骤SA中的第二式计算得到河道首、末断面的边界界面处当前时刻即第n步对应时刻的数值通量；其中河道首、末断面的边界界面即为河道的第1个界面和最后一个界面。

步骤SC、获取上一时刻即第n-1步对应时刻所计算的河道各断面在前时刻即第n步对应时刻的水位和流量；然后通过步骤SA中的第二式计算得到河道各断面当前时刻即第n步对应时刻各控制体之间界面的数值通量；基于当前时刻即第n步对应时刻各控制体之间界面的数值通量，通过然后通过步骤SA中的第一式计算得到河道各断面在下一时刻即第n+1步对应时刻的水位和流量。

在本步骤中，所使用的水动力方程为采用一维圣维南方程作为河道水流控制方程建立一维水动力模型，如下式：

式中：x为里程；t为时间；Z为水位；B为过水断面水面宽度；Q为流量；q为侧向单宽流量，正值表示流入，负值表示流出；A为过水断面面积；g为重力加速度；u为断面平均流速；β为校正系数；R为水力半径；c为谢才系数，c＝R^1/6/n′，n′为曼宁糙率系数；

步骤S4、针对于河网中各汊点，获取通过步骤S3计算得到的该汊点处各河道断面的流量，判断该汊点处各河道断面的流量是否满足流量守恒条件；具体为判定以下公式是否成立：

其中，p为汊点编号，j为汊点处的河道编号；Q_p,j为第p个汊点处第j条河道的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量；

若是，则将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，则进入步骤S5；

步骤S5、首先判断当前迭代次数是否达到一定值Z；在本实施例中Z设置为50，即设置最大迭代次数为50次；

若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，首先根据汊点处流入流量与流出流量的差值计算汊点水位校正值，然后判断该水位校正值是否小于一定值Z1，在本实施例中Z1值设置为10^-4米；若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为河网水动力计算结果；若否，则根据汊点水位校正值对汊点处水位值进行校正，将校正后的汊点处水位值作为下一次迭代计算时汊点的预测水位值，然后返回步骤S3。

本实施例方法上述步骤S5中，根据汊点处流入流量与流出流量的差值计算汊点水位校正值为：

其中p为汊点编号，j为汊点处的河道编号，m为当前迭代次数，即校正步数；Q_p,j,m为第p个汊点处第j条河道在第m次迭代计算时得到的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量；ΔZ_p,m为第m次迭代计算时根据汊点处流入流量与流出流量的差值计算汊点水位校正值；

δ_p,m为第m次迭代计算时第p个汊点的校正系数，具体为：

本实施例上述步骤S5中，根据汊点水位校正值对汊点当前的预测水位值进行校正，得到校正后的汊点处水位值为：

为汊点当前的预测水位值，为校正后的汊点处水位值；其中：

为下一次迭代即第m+1次迭代计算时汊点的预测水位值。

在本实施例中，如图2所示，以福州市江北城区水系为例进行说明。福州江北城区共概化67条河道，总河长约135km；概化的湖泊、调蓄池包括琴亭湖、温泉公园湖、斗门调蓄池等共7个。概化的水闸共36个，其中沿江水闸8个，内河水闸16个，湖泊水闸12个。概化泵站12个，其中沿江泵站6个，内河泵站和湖泊泵站各3个。福州江北城区河网共计1317个断面、127条河段、70个汊点。福州市江北城区水系包括新店片水系、晋安河水系、白马河水系、磨洋河水系和光明港水系。其中，光明港和白马河水系属于典型的平原感潮河网；新店片水系为典型的山区河道，河宽较窄且坡度很陡。

如图3所示为新店溪河段的深泓高程，由于坡降为深泓高程与里程的比值，因此由图3可计算得到新店溪河段的坡降约为2％，将本实施例上述复杂河网水动力计算方法应用于新店片水系；在本实施例中，将新店溪河段出口断面的水位边界设置为4.96m，相应过水断面面积为1.0m²，平均水深0.22m。图4表示为基于本发明方法的不同坡降下所对应的水面线，表1表示为不同坡降下所对应的洪峰流量。

表1

设计频率	P＝0.1％	P＝2％	P＝50％
				洪峰流量(m<sup>3</sup>/s)	353	201	54

根据表1中不同量级流量，新店溪河段出口断面的水位计算值为

5.88m(Q＝54m³/s)～7.39m(Q＝353m³/s)。由图4和表1可得，不同坡降的折线与深泓高程折线方向基本相同，且上游断面水位高于下游断面水位，流量越大的情况下断面的水位越高，因此符合水流运动规律。

针对遇洪情况下的北城区水系72h水动力过程进行分析，表2为不同检测方法所用的计算耗时时长。

表2

方法	有限差分法	有限体积法	本实施例方法
				计算耗时(s)	——	216	109

由表2可知，由于有限差分法无法处理陡坡山区河道，故检测失败；有限体积法的计算耗时为216s；本实施例的复杂河网水动力计算方法计算耗时为109s，检测速度约为有限体积法的2倍。可见，本实施例复杂河网水动力检测方法实现了速度与稳定性的优化组合。

实施例2

本实施例公开了一种计算设备，包括处理器以及用于存储处理器可执行程序的存储器，处理器执行存储器存储的程序时，实现实施例1所述的复杂河网水动力计算方法，具体如下：

步骤S1、针对河网中各汊点，对各汊点分别预设一个水位值作为其初始的预测水位值；

步骤S2、根据河道水面坡降将河道分类为缓坡河道和陡峭河道；

步骤S3、根据河网中各汊点当前的预测水位值，计算各河道各断面的水动力参数；其中：

针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数；

针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数；

步骤S4、针对于河网中各汊点，获取通过步骤S3计算得到的该汊点处各河道断面的流量，判断该汊点处各河道断面的流量是否满足流量守恒条件；

若是，则将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，则进入步骤S5；

步骤S5、首先判断当前迭代次数是否达到一定值；

若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，首先根据汊点处流入流量与流出流量的差值计算汊点水位校正值，然后判断该水位校正值是否小于一定值Z1；若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为河网水动力计算结果；若否，则根据汊点水位校正值对汊点当前的预测水位值进行校正，将校正后的汊点处水位值作为下一次迭代计算时汊点的预测水位值，然后返回步骤S3。

在本实施例中，计算设备可以是台式电脑、笔记本电脑、智能手机、PDA手持终端、平板电脑或其他具有处理器功能的终端设备。

实施例3

本实施例公开了一种存储介质，存储有程序，所述程序被处理器执行时，实现实施例1所述的复杂河网水动力计算方法，具体如下：

步骤S1、针对河网中各汊点，对各汊点分别预设一个水位值作为其初始的预测水位值；

步骤S2、根据河道水面坡降将河道分类为缓坡河道和陡峭河道；

步骤S3、根据河网中各汊点当前的预测水位值，计算各河道各断面的水动力参数；其中：

针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数；

针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数；

步骤S4、针对于河网中各汊点，获取通过步骤S3计算得到的该汊点处各河道断面的流量，判断该汊点处各河道断面的流量是否满足流量守恒条件；

若是，则将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，则进入步骤S5；

步骤S5、首先判断当前迭代次数是否达到一定值；

若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，首先根据汊点处流入流量与流出流量的差值计算汊点水位校正值，然后判断该水位校正值是否小于一定值Z1；若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为河网水动力计算结果；若否，则根据汊点水位校正值对汊点当前的预测水位值进行校正，将校正后的汊点处水位值作为下一次迭代计算时汊点的预测水位值，然后返回步骤S3。

本实施例中的存储介质可以是磁盘、光盘、计算机存储器、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、U盘、移动硬盘等介质。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种复杂河网水动力计算方法，其特征在于，步骤如下：

步骤S1、针对河网中各汊点，对各汊点分别预设一个水位值作为其初始的预测水位值；

步骤S2、根据河道水面坡降将河道分类为缓坡河道和陡峭河道；

步骤S3、根据河网中各汊点当前的预测水位值，计算各河道各断面的水动力参数；其中：

针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数；

针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数；

步骤S4、针对于河网中各汊点，获取通过步骤S3计算得到的该汊点处各河道断面的流量，判断该汊点处各河道断面的流量是否满足流量守恒条件；

若是，则将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，则进入步骤S5；

步骤S5、首先判断当前迭代次数是否达到一定值；

若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为当前河网水动力计算结果；

若否，首先根据汊点处流入流量与流出流量的差值计算汊点水位校正值，然后判断该水位校正值是否小于一定值Z1；若是，则终止迭代，将步骤S3中当前计算得到的各河道各断面的水动力参数作为河网水动力计算结果；若否，则根据汊点水位校正值对汊点当前的预测水位值进行校正，将校正后的汊点处水位值作为下一次迭代计算时汊点的预测水位值，然后返回步骤S3。

2.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，每一时段执行一次步骤S1至S5；即在当前时段执行完成步骤S1至S5后，进入下一时段重新执行步骤S1至S5，直到最终时段；

步骤S2中根据河道水面坡降将河道分类为缓坡河道和陡峭河道的具体过程如下：

步骤S21、获取上一时段的河网水动力计算结果，得到上一时段各河道各断面的水动力参数，包括各断面的水位；

步骤S22、针对于各河道，计算其中每相邻两个断面之间的水面坡降：

其中，S_i为河道内第i、i+1个断面之间的水面坡降，Z_i、Z_i+1分别为上一时段河道内第i、i+1个断面的水位；Δx_i为河道内第i、i+1个断面之间的距离；M为河道内断面的总数；

步骤S23、比较步骤S22计算得到的各每相邻两个断面之间的水面坡降，选取出其中一个最大值：

S_max＝max(|S₁|,...,|S_i|,...,|S_M-1|)；

其中，S_max为选取出的最大值；M为河道内断面的总数；

步骤S24、根据步骤S23中选取出的最大值S_max对河道进行分类，当S_max大于等于阈值Z2时，判断河道为陡峭河道，否则判断河道为缓坡河道。

3.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，每一时段执行一次步骤S1至S5；即在当前时段执行完成步骤S1至S5后，进入下一时段重新执行步骤S1至S5，直到最终时段；

当前时段在执行步骤S1中，将上一时段最后计算得到的汊点的预测水位值，作为当前时段各汊点初始的预测水位值。

4.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，步骤S4中通过判断以下公式是否成立来判断汊点处各河道断面的流量是否满足流量守恒条件：

其中，p为汊点编号，j为汊点处的河道编号；Q_p,j为第p个汊点处第j条河道的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量。

5.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，步骤S3中针对于缓坡河道，采用有限差分法计算缓坡河道各断面的水动力参数具体过程如下：

步骤S31、对单一缓坡河道的水动力方程进行有限差分离散，得第一式如下：

式中：i为河道内断面编号；为第n+1步对应时刻的河道内第i个断面的水位、流量；a_1i、a_2i、c_1i、c_2i、d_2i、E_1i和E_2i为已知的差分系数；n为计算步数；M为河道内断面的总数；

采用追赶法对河道内各断面流量进行递推，得第二式：

式中：为第n+1步对应时刻河道首断面的水位、流量；为第n+1步对应时刻的河道末断面的水位、流量；E₁、E₂、R₁、R₂、T₁、T₂为已知的递推系数；

步骤S32、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道首断面的水位，根据河道后一个汊点当前的预测水位值得到河道末断面的水位；根据河道首断面的水位和河道末断面的水位，通过步骤S31中的第二式计算得到河道首、末断面的流量；

步骤S33、根据河道首断面的水位、河道末断面的水位、河道首断面的流量以及河道末断面的流量，通过步骤S31中的第一式计算出河道其他断面的水位和流量；同时，根据河道各断面的水位、流量以及断面地形计算出各断面的流速。

6.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，步骤S3中针对于陡峭河道，采用有限体积法计算陡峭河道各断面的水动力参数的具体过程如下：

步骤SA、对单一陡峭河段的水动力方程进行有限体积离散，得第一式：

其中，下标i为河段内第i个断面对应的控制体编号；n为计算步数；为第n步对应时刻的第i个断面的状态向量；为系数矩阵；Δt为计算时间步长；Δx_i为第i个断面对应的控制体长度；分别为第n步对应时刻的第i个断面对应的控制体的左、右两侧界面的数值通量；

其中，为源项近似：

式中，为第n步对应时刻的第i个断面的侧向单宽流量；g为重力加速度；为第n步对应时刻的第i个断面的过水断面面积；为第n步对应时刻的第i个断面的水位；为第n步对应时刻的第i个断面的流量；为第n步对应时刻的第i个断面的水力半径；n′为曼宁糙率系数；

采用HLL格式计算界面数值通量，如第二式：

式中：为第n步对应时刻的第i个断面对应的控制体的左侧界面的数值通量，即第i-1与i断面之间的第n步对应时刻的界面数值通量；为第i-1与i断面之间的第n步对应时刻的界面的左、右波速；上标L、R分别代表界面左、右侧；和为基于界面左、右侧重构变量的数值通量，按下式计算：

由下式计算：

式中：上、下标L和R代表界面左、右侧；波速由下式计算：

式中：为第i断面第n步对应时刻的流速；为第i断面第n步对应时刻的特征速度；为第i断面第n步对应时刻的过水断面面积；为第i断面第n步对应时刻的过水断面水面宽度；

步骤SB、获取河道前后两个汊点当前的预测水位值，根据河道前一个汊点当前的预测水位值得到河道的上游水位边界条件，根据河道后一个汊点当前的预测水位值得到河道的下游水位边界条件；根据上述上、下游水位边界条件，通过步骤SA中的第二式计算得到河道首、末断面的边界界面处当前时刻即第n步对应时刻的数值通量；

步骤SC、获取上一时刻即第n-1步对应时刻所计算的河道各断面在当前时刻即第n步对应时刻的水位和流量；然后通过步骤SA中的第二式计算得到河道各断面当前时刻即第n步对应时刻各控制体之间界面的数值通量；基于当前时刻即第n步对应时刻各控制体之间界面的数值通量，通过然后通过步骤SA中的第一式计算得到河道各断面在下一时刻即第n+1步对应时刻的水位和流量。

7.根据权利要求1所述的复杂河网水动力计算方法，其特征在于，步骤S5中，根据汊点处流入流量与流出流量的差值计算汊点水位校正值为：

其中p为汊点编号，j为汊点处的河道编号，m为当前迭代次数，即校正步数；Q_p,j,m为第p个汊点处第j条河道在第m次迭代计算时得到的流量，流入为正，流出为负；k_in、k_out分别为流入、流出汊点的河道数量；ΔZ_p,m为第m次迭代计算时根据汊点处流入流量与流出流量的差值计算汊点水位校正值；

δ_p,m为第m次迭代计算时第p个汊点的校正系数，具体为：

步骤S5中，根据汊点水位校正值对汊点当前的预测水位值进行校正，得到校正后的汊点处水位值为：

为汊点当前的预测水位值，为校正后的汊点处水位值；其中：

为下一次迭代计算时汊点的预测水位值。

8.根据权利要求7所述的复杂河网水动力计算方法，其特征在于，Z1的值为10^-4米。

9.一种计算设备，包括处理器以及用于存储处理器可执行程序的存储器，其特征在于，所述处理器执行存储器存储的程序时，实现权利要求1-8任一项所述的复杂河网水动力计算方法。

10.一种存储介质，存储有程序，其特征在于，所述程序被处理器执行时，实现权利要求1-8任一项所述的复杂河网水动力计算方法。