CN103218538A - 基于河网关联矩阵的河网一维恒定流计算方法 - Google Patents

Abstract

本发明提出了一种基于河网关联矩阵的河网一维恒定流计算方法，包括根据河网结构绘制相应的河网有向图；根据所述河网有向图构建河网关联矩阵；根据所述河网关联矩阵构造河网一维恒定流综合方程；消除所述河网一维恒定流综合方程的奇异性，并通过牛顿迭代法求解消除奇异性的所述河网一维恒定流综合方程。本发明构造的河网关联矩阵可用于表达任意结构的河网信息，因此适用于单一河道、无环和有环等所有形式的河网一维恒定水流求解；只需要构造一个河网关联矩阵，并通过牛顿迭代法求解，具有较快的收敛速度，可有效减少计算量。

Description

基于河网关联矩阵的河网一维恒定流计算方法

技术领域

本发明涉及一种基于河网关联矩阵的河网一维恒定流计算方法，属于河网水流计算技术领域。

背景技术

自然界中河网结构普遍存在，从河网拓扑结构角度，河网可分为无环河网和有环河网，对这两类河网内一维恒定水流的求解，是河流动力学领域的重要问题。其中，对河网内水力要素（如流量、流速、水位等）的获取，不仅是解决各类水利问题的关键，还是以非耦合方式解决水生生物迁移、污染物扩散等生态和环境问题的先决条件。

在无外力干扰的情况下，自然界中水的流动是由水位差造成的，同一流路内不同两点之间的水位差，使得水在该流路内从高处流向低处，这种水位差实质是水的重力势能差，只要流路通畅，其势能总要转化为动能，水体的重力势能差是水体流动的根本原因。

现有技术对河网一维恒定流的求解，已有利用“连接矩阵”和“迴路矩阵”的求解方法。但此方法只适用于求解有环河网一维恒定水流，并且要构造“连接矩阵”和“迴路矩阵”共两类矩阵反映河网结构，在通过这两类矩阵分别列出节点水流连续方程和迴路方程后，才能通过“迴路法”求解河网内的水力要素。

但是，如前所述，“二矩阵法”仅适用于有环河网水流的求解，对于无环河网并不适用。而且由于“二矩阵法”在求解有环河网水流时，需要根据河网结构构造“连接矩阵”和“迴路矩阵”两个矩阵，导致求解的前期工作量较大。

发明内容

本发明为解决现有河网一维恒定流求解方法存在的不适用于无环河网、前期工作量较大等问题，提出了一种基于河网关联矩阵的河网一维恒定流计算方法。本发明所采用的技术方案如下：

一种基于河网关联矩阵的河网一维恒定流计算方法，包括：

(1)根据河网结构绘制相应的河网有向图；

(2)根据所述河网有向图构建河网关联矩阵；

(3)根据所述河网关联矩阵构造河网一维恒定流综合方程；

(4)消除所述河网一维恒定流综合方程的奇异性，并通过牛顿迭代法求解消除奇异性的所述河网一维恒定流综合方程。

本发明的具体实施方式构造的河网关联矩阵可用于表达任意结构的河网信息，因此适用于单一河道、无环和有环等所有形式的河网一维恒定水流求解；只需要构造一个河网关联矩阵，并通过牛顿迭代法求解，具有较快的收敛速度，既减少了前期工作量，又可有效减少计算量。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面对实施例描述中所使用的附图作简单地介绍，显而易见地，以下描述的附图仅用于对本发明实施例的阐述，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他附图。

图1是本发明提出的基于河网关联矩阵求解河网一维恒定流计算方法的具体实施过程示意图；

图2是本发明的具体实施方式提供的有环河网和无环河网的示意图，其中图2a表示有环河网，图2b表示无环河网；

图3是本发明的具体实施方式提供的有环河网和无环河网的有向图；

图4是本发明的具体实施方式提供的有向图节点上的河网断面示意图；

图5是本发明的具体实施方式提供的河段示意图；

图6是本发明的具体实施方式提供的求解河网一维恒定水流的具体求解流程示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明技术方案进行清楚、完整地描述，显然，以下描述的实施例仅是本发明的部分实施例，并非全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的具体实施方式提供了一种基于河网关联矩阵的河网一维恒定流计算方法，如图1所示，包括:

步骤1，根据河网结构绘制相应的河网有向图。

具体的，有向图是用于表示物件之间的关系的拓扑结构，以G代表有向图，其数学定义可写为G=(V,E)，V和E分别为节点和有向边集合。图3a和图3b分别是与图2a和图2b对应的有环河网和无环有向图，有环有向图中存在环路，无环有向图仅由支汊组成。若以|V|表示有向图G中的节点数，以|E|表示G中的有向边数，则图3a是一个|V|=4、|E|=5的有环有向图，图3b是一个|V|=6、|E|=5的无环有向图。

在以有向图描述河网结构时，有向边与两断面间的河段是“1对1”的关系，节点与河网断面一般是“1对n”的关系(如图4所示)，n是在节点上交汇的河段总数。对于单一河道，其首尾两个节点上的断面分别由首尾两个河段独占，其内每个节点上，一般只对应1个由上下游河段共享的断面。以有向图描述河网的方法如下：

(1)确定河网的节点总数|V|；

(2)确定河网的微河段总数|E|；

(3)以圆代表节点，根据河段的空间分布绘制连结各节点的有向边，并任意指定各有向边的正方向，以从1开始的连续自然数为节点编号，保证每个节点均具有唯一编号，以同样的方法给各条有向边编号，即得反映河网结构的有向图G。

河网有向图中河段方向可任意规定，当河段内流量计算结果为负时，表明该河段内实际水流方向与河段的规定方向相反，这并不影响河段内的真实流向和流量结果。

步骤2，根据所述河网有向图构建河网关联矩阵。

具体的，在绘制出河网结构的有向图后，可根据河网有向图构建河网关联矩阵，此矩阵包含了河网结构的全部信息。以图3a为例，此有环有向图的关联矩阵A如式(1)示：

式(1)中矩阵A有4列5行，表示图3a所示的有环有向图由4个节点和5条有向边组成，关联矩阵A的第i行与有向图的第i条有向边相应，第j列与有向图的第j个节点相应，关联矩阵A的第i行第j列元素a_i,j的确定方法如式(2)所示：

关联矩阵A反映了有向图中节点和有向边的关系，如式(1)中矩阵A的第一行表示有环路河网图3a的1号有向边的起点为节点1、终点为节点2。

步骤3，根据所述河网关联矩阵构造河网一维恒定流综合方程。

具体的，若以A表示任意河网的关联矩阵，则该河网的水流综合方程构造如式(3)所示：

$A^T\mathrm{diag}\left(\frac{1}{N_i\left|Q_i\right|}\right)\mathrm{AZ}=q,(i=\mathrm{0,1},...|E\left|\right)---\left(3\right)$

式(3)及其中各量的含义见以下推导过程：

求解河网水流即求解河网中各河段内流量Q和各节点的水位Z，这要基于水流连续方程和能量方程。

对于河网内的节点j，其水流连续方程如式(4)所示：

$\underset{i=1}{\overset{\left|E\right|}{\mathrm{\Σ}}}{a}_{j,i}{Q}_i=q_j,(j=1,..,|V\left|\right)---\left(4\right)$

式(4)中，i是河段索引，j是节点索引，a_j,i是河网关联矩阵第j行第i列元素值，Q_i是河段i内流量，q_j是节点j上的流量边界条件，即该节点与河网外部的流量交换值。

对于节点j₁和节点j₂之间的河段（如图5所示），其能量方程为：

$Z_{j_1}+{\mathrm{\α}}_{j_1}\frac{V_{j_1}^2}{2g}=Z_{j_2}+{\mathrm{\α}}_{j_2}\frac{V_{j_2}^2}{2g}+h_w---\left(5\right)$

式(4)中，Z为水位，α为动能校正系数，α的值取决于过水断面上流速分布的不均匀程度，流速越不均匀，其值越大。对于一般的均匀流或渐变流，α=1.05～1.1，可近似取1，对于流速极不均匀的情况，α可取2或更大的值。

式(5)中g为重力加速度，h_w是沿程水头损失与局部水头损失之和，可由式(6)计算。

$h_w=\frac{\mathrm{\Δ}{x}_i}{{\stackrel{\‾}{K_i}}^2}{Q}_i^2+\frac{{\mathrm{\ζ}}_i(A_{j_1}^2-A_{j_2}^2)}{2g{A}_{j_1}^2{A}_{j_2}^2}{Q}_i^2=(\frac{\mathrm{\Δ}{x}_i}{{\stackrel{\‾}{K_i}}^2}+\frac{{\mathrm{\ζ}}_i(A_{j_1}^2-A_{j_2}^2)}{2g{A}_{j_1}^2{A}_{j_2}^2})Q_i^2---\left(6\right)$

式(6)中i为河段号，

为河段平均流量模数，Δx为河段长，ζ为局部水头损失系数，A_j为过水面积。

可由式(7)和式(8)计算。

$\stackrel{\‾}{K}=\frac{K_{j_1}+K_{j_2}}{2}---\left(7\right)$

$K_j=C_j{A}_j\sqrt{R_j},(j=j_1,j_2)---\left(8\right)$

式(7)中C_j=R^1/6/n，为第j断面的谢才系数，R为水力半径，n为河段糙率。式(6)中局部水头损失系数ζ的取值规则如下。

将式(6)代入式(5)并整理可得：

$\mathrm{\Δ}{Z}_i=Z_{j_1}-Z_{j_2}=N_i{Q}_i^2---\left(10\right)$

若考虑水流方向，则有：

$\mathrm{\Δ}{Z}_i=Z_{j_1}-Z_{j_2}=S_i{N}_i{\left|Q_i\right|}^2\⇒Q_i=\mathrm{\Δ}{Z}_i/\left(N_i\right|Q_i\left|\right)---\left(11\right)$

其中S_i是符号函数，当Q_i≥0时S_i=1,否则S_i=-1。N_i反映了河段对水流的制约作用，由式(12)计算。

$N_i=\frac{1}{2g}(\frac{{\mathrm{\α}}_{j_2}}{A_{j_2}^2}-\frac{{\mathrm{\α}}_{j_1}}{A_{j_1}^2})+\frac{\mathrm{\Δ}{x}_i}{{\stackrel{\‾}{K_i}}^2}+\frac{{\mathrm{\ζ}}_i(A_{j_1}^2-A_{j_2}^2)}{2g{A}_{j_1}^2{A}_{j_2}^2}---\left(12\right)$

河网恒定流控制方程是式(4)和式(5)，这两式可分别以河网关联矩阵A表达，即：

A^TQ＝q(13)

$\mathrm{AZ}=\mathrm{\ΔZ}=\mathrm{SN}{\left|Q\right|}^2\⇒Q=\mathrm{diag}\left(\frac{1}{N_i\left|Q_i\right|}\right)\mathrm{AZ}---\left(14\right)$

上式中A^T是河网关联矩阵A的转置，向量Q、q和ΔZ的具体形式如式(15)。

$Q=\left[\begin{array}{c}{Q}_1\\ ...\\ Q_i\\ ...\\ Q_{\left|E\right|}\end{array}\right],q=\left[\begin{array}{c}{q}_1\\ ...\\ q_j\\ ...\\ q_{\left|V\right|}\end{array}\right],\mathrm{\ΔZ}=A\left[\begin{array}{c}{Z}_1\\ ...\\ Z_j\\ ...\\ Z_{\left|V\right|}\end{array}\right]=\left[\begin{array}{c}{S}_1{N}_1|Q_1{|}^2\\ ...\\ S_i{N}_i|Q_i{|}^2\\ ..\\ S_{\left|E\right|}{N}_{\left|E\right|}|Q_{\left|E\right|}{|}^2\end{array}\right]---\left(15\right)$

式(14)和式(15)中，向量Q的元素Q_i是河网中河段i的过流量，其值为正表明水流方向与该河段的规定方向相同，为负则反之。向量q中的元素qj表示河网节点j与外部的流量交换值，q_j>0时表明j节点有流量汇入，q_j<0则反之，q_j=0表示j节点与河网外部无流量交换。ΔZ是河网中各河段的水位落差向量。

将式(14)代入式(13)可得：

$A^T\mathrm{diag}\left(\frac{1}{N_i\left|Q_i\right|}\right)\mathrm{AZ}=q,(i=\mathrm{0,1},...|E\left|\right)$

上式即为以河网关联矩阵表达的河网恒定流综合方程(3)，该方程既反映了节点处交汇河段的流量平衡，又反映了各河段自身应满足的能量方程。

步骤4，消除所述河网一维恒定流综合方程的奇异性，并通过牛顿迭代法求解消除奇异性的所述河网一维恒定流综合方程。

具体的，求解河网水流的定解条件，主要有地形条件和水文条件两类。对于一维数学模型，地形条件以断面间距、断面点的起点距和高程对的形式给出。水文条件包括流量和水位两部分，流量边界以式(15)中向量q的形式给出，水位条件一般在河网的出口断面指定。

求解方程(3)需要解决两个问题：一是方程(3)是非线性的；二是河网结构的关联矩阵A是奇异的。对于第一个问题，可选取河网出口断面的已知水位Z作为基准高程值，即规定出口断面水位为0，由此便可删除方程(3)系数矩阵与出口节点相应的列，从而可得式(16):

$A^{*T}\mathrm{diag}\left(\frac{1}{N_i\left|Q_i^k\right|}\right)A^{*}{Z}^{*}=q^{*},(i=\mathrm{0,1},...|E\left|\right)---\left(16\right)$

式(16)中带*号上标的矩阵和向量，均比方程组(3)中的相应量少了1列或1行，经如此处理后，可保证方程组(16)的系数矩阵是正定的。

对于第二个问题，可借助牛顿迭代法。由于方程(16)所求解的未知量是各节点的水位Z，Z又与方程(16)的系数矩阵中的Q相关，由此可知需要构造的迭代式应是Z和Q之间的关系式。根据牛顿迭代法，可由式(10)构造河网中河段i内流量Q_i的函数F(Q_i)，如式(17)示:

$F\left(Q_i\right)=S_i{N}_i{\left|Q_i\right|}^2-\mathrm{\&Delta;}{Z}_i,(i=1...|E\left|\right)---\left(17\right)$

满足方程组(16)的Q和Z一定是满足式(10)的，这就意味着方程组(16)的解应使式(17)定义的函数F取0值。求解之初根据节点处的流量连续方程分配给河网内各河段的流量值多存在误差，因此需要迭代试算。对于第k次和第k+1次迭代，根据Taylor公式有：

$F\left(Q_i^{k+1}\right)=F\left(Q_i^k\right)+\mathrm{\&Delta;}{Q}_i\&times;F^{\&prime;}\left(Q_i^k\right)---\left(18\right)$

迭代应使F

趋于0，故由式(18)得:

$\mathrm{\&Delta;}{Q}_i=-\frac{F\left(Q_i^k\right)}{F^{\&prime;}\left(Q_i^k\right)}---\left(19\right)$

式(17)两端对Q_i求导，可得:

$F^{\&prime;}\left(Q_i\right)=2S_i^2{N}_i\left|Q_i\right|=2N_i\left|Q_i\right|---\left(20\right)$

将式(20)代入式(19)，整理可得:

$Q_i^{k+1}=\frac{1}{2}(Q_i^k+\frac{\mathrm{\&Delta;}{Z}_i}{N_i\left|Q_i^k\right|}),(k=\mathrm{0,1},...)---\left(21\right)$

式(21)即为河网中各河段流量的迭代公式，式中k为迭代次数。

下面通过具体的实施例对本具体实施方式提供的基于河网关联矩阵的河网一维恒定流计算方法作详细说明。

有环路河网水流的求解流程如图6示，以图3a中的有环路河网为例，规定其出口断面4的水位Z=0，则式(16)中的各矩阵和向量应为式(22)所示的形式。

$A^{*T}=\left[\begin{array}{ccccc}1& 1& 0& 0& 0\\ -1& 0& 1& 1& 0\\ 0& -1& -1& 0& 1\end{array}\right],Z^{*}=\left[\begin{array}{c}{Z}_1\\ Z_2\\ Z_3\end{array}\right],q^{*}=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\end{array}\right]---\left(22\right)$

将式(22)中各项代入式(16)，并按照图6的计算流程，便可算得以出口断面水位Z为基准高程的其他各断面水位值和河网内的各河段流量。表1给出了求解图2a中的有环路河网水流的概化算例所采用的定解条件和计算结果，概化算例中，假定各断面具有相同的断面形态。

表1求解图1a中的有环路河网水流所采用的定解条件和计算结果

表1中的底高条件和水位结果是以出口断面的已知水位为基准的。经验证，表1中的水位和流量结果满足式(11)至式(13)和式(16)等河网水流控制方程。河段3的计算流量为负值，说明该河段内的水流方向与河网有向图中有向边的规定方向相反，这与水位结果是一致的。

当给定的边界条件是无环路河网节点的水位时，无环路河网水流的求解流程同上。

本发明实施方式提供的基于河网关联矩阵的河网一维恒定流计算方法与现有技术相比具有以下的优点：

(1)本发明适用范围广，现有“二矩阵法”仅适用于有环河网的一维恒定流求解，本发明适用于单一河道、无环和有环等所有形式的河网一维恒定水流求解。本优点的实现，是由于本发明构造的“河网关联矩阵”可用于表达任意结构的河网信息，这是现有技术所不具备的能力。

(2)本发明工作量小，与现有的“二矩阵法”相比，本发明只需要构造一个河网关联矩阵，并且此矩阵与二矩阵法中的“连接矩阵”是同型矩阵，构造过程工作量相同，但本发明节俭了现状技术构造“迴路矩阵”的工作量。本优点的实现，同样是由于本发明的“河网关联矩阵”信息量更全面，它是河网结构的数学本质。

(3)本发明使用牛顿迭代法求解，与现有技术的“迴路法”相比，收敛速度更快。本优点的实现，是由于在本发明中更合理的应用了Taylor展开式等数学原理，使得求解非线性方程组所需要的流量迭代的理论基础更强，且牛顿迭代法本身更具普适性，均使本发明的方法更先进。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明实施例揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种基于河网关联矩阵的河网一维恒定流计算方法，其特征在于，包括：

(1)根据河网结构绘制相应的河网有向图；

(2)根据所述河网有向图构建河网关联矩阵；

(3)根据所述河网关联矩阵构造河网一维恒定流综合方程；

(4)消除所述河网一维恒定流综合方程的奇异性，并通过牛顿迭代法求解消除奇异性的所述河网一维恒定流综合方程。

2.根据权利要求1所述的方法，其特征在于，所述根据河网结构绘制相应的河网有向图包括：

(1)确定相应河网的节点总数；

(2)确定相应河网的微河段总数；

(3)以圆代表节点，根据河段的空间分布绘制连结各节点的有向边，并任意指定各有向边的正方向，以从1开始的连续自然数为节点编号，保证每个节点均具有唯一编号，以同样的方法给各条有向边编号，制得所述河网有向图。

3.根据权利要求1所述的方法，其特征在于，所述根据所述河网有向图构建河网关联矩阵包括：

将所述河网有向图的第i条有向边与所述河网关联矩阵的第i行相对应，将所述河网有向图的第j个节点与所述河网关联矩阵的第j列相对应，其中i和j都属于自然数，所述河网关联矩阵的第i行第j列元素a_i,j的确定方法为：

4.根据权利要求1所述的方法，其特征在于，所述河网一维恒定流综合方程为：

$A^T\mathrm{diag}\left(\frac{1}{N_i\left|Q_i\right|}\right)\mathrm{AZ}=q,(i=\mathrm{0,1},...|E\left|\right)$

其中，矩阵A表示是河网关联矩阵，A^T表示是河网关联矩阵A的转置，Q_i表示是河网中河段i的过流量，Z表示各节点的水位向量。

5.根据权利要求1所述的方法，其特征在于，所述消除所述河网一维恒定流综合方程的奇异性包括：

选取河网出口断面的已知水位Z作为基准高程值，并规定出口断面水位为0，则将所述河网一维恒定流综合方程变换为：

$A^{*T}\mathrm{diag}\left(\frac{1}{N_i\left|Q_i^k\right|}\right)A^{*}{Z}^{*}=q^{*},(i=\mathrm{0,1},...|E\left|\right).$

6.根据权利要求5所述的方法，其特征在于，所述通过牛顿迭代法求解的迭代公式为：

$Q_i^{k+1}=\frac{1}{2}(Q_i^k+\frac{\mathrm{\&Delta;}{Z}_i}{N_i\left|Q_i^k\right|}),(k=\mathrm{0,1},...)$

其中，k表示迭代次数。

Images