CN102890732B - 一种含多种结构物渠道的水动力条件并行化数值模拟方法 - Google Patents

Abstract

本发明涉及一种含多种结构物渠道的水动力条件并行化数值模拟方法，它假设结构物处流量，然后根据给定的结构物处的流量与水位关系或者结构物上侧的水位对假设的流量进行校正，并根据渠道结构物两侧流量相等为结构物上下游渠段分别提供下边界条件和上边界条件，从而为每一渠段建立一封闭的方程组，且每一方程组中只含同一渠段的变量，使采用隐式差分法离散而成的系数矩阵仍然保持带状特征，物理上相互联系的渠段在数值上相对独立，再利用带状矩阵求解出渠段变量，通过Newton-Raphson迭代法校正结构物处流量，使流量达到调控目标。根据本发明，原求解单一渠道水动力条件的程序无需大范围改变，就可以实现含多种结构物的渠道的水动力条件的并行求解，达到水量调控的目标。

Description

一种含多种结构物渠道的水动力条件并行化数值模拟方法

技术领域

本发明涉及一种水动力数值计算方法，特别是关于一种求解调水或灌溉渠道水动力条件的含多种结构物渠道的水动力条件并行化数值模拟方法。

背景技术

在河流水利工程中，借助调水渠道或者灌溉渠道能够人为地将水从水资源丰富的区域调往水资源贫乏的区域，从而有效缓解局部地区水资源严重短缺的情况，新中国成立以来的超大型工程——南水北调就是一个典型的应用案例。由于调水渠道和灌溉渠道空间跨度很大，因此在调水过程中常常需要在渠道中设置闸门或电站等多种类型的结构物。求解含有多种类型结构物的渠道的水动力条件则是调水渠道和灌溉渠道实现宏观调度和流量控制的一个必要条件。当前，对于平原渠道水动力条件的分析计算大多采用Saint-Venant双曲方程组建模，通过Preissmann隐式差分法离散成带状矩阵，并利用Newton-Raphson迭代法进行求解。但是这种方法只适于求解单一渠道的流量，对于渠道中存在如闸门或电站等结构物的情况，则需要将渠道在结构物处分隔成两个渠段，按顺序逐级求解，即在求解完第一级渠段得到其下边界条件后，将其作为第二级渠段的上边界条件再进行第二级渠段的计算。又或者采用另一种方法，将渠道的结构物处理成局部损失，然后整体渠道一起求解。这两种方法的求解过程耗时都比较长，效率较低，且当渠道中的结构物发生改变时，需要一定范围地修改相应的计算程序代码，程序的通用性和可移植性差。虽然目前已有一种并行化求解方法可以解决上述问题，参见《一种含梯级水电站河流水动力条件并行化数值模拟方法》（申请号201210172062.1），但是这种方法只适用于求解含有流量只跟上游水位有关的结构物的多级河道的水动力条件。

发明内容

针对上述问题，本发明的目的是提供一种含多种结构物渠道的水动力条件并行化数值模拟方法，基于此方法，无需大范围修改原单一渠道水动力条件求解程序，就能够快速求解含多种结构物的调水渠道或者灌溉渠道的水动力条件。

为实现上述目的，本发明采取以下技术方案：一种含多种结构物渠道的水动力条件并行化数值模拟方法，它包括以下步骤：

1）给定渠道中每一结构物处流量-水位关系式Q_i＝g′_i(Z_Gi，Z_Hi)，下标i表征结构物序号，i≥1，采用Saint-Venant双曲方程组为整条渠道建立数学模型方程组；

2）采用Preissmann隐式差分格式离散步骤1）获得的渠道数学模型方程组，并通过Newton-Raphson迭代法获得对应的线性离散方程组：

$a_{2j-\mathrm{1,1}}{\mathrm{\ΔQ}}_j+a_{2j-\mathrm{1,2}}{\mathrm{\ΔA}}_j+a_{2j-\mathrm{1,3}}{\mathrm{\ΔQ}}_{j+1}+a_{2j-\mathrm{1,4}}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^C\right)=0,$

$a_{2j,1}{\mathrm{\ΔQ}}_j+a_{2j,2}{\mathrm{\ΔA}}_j+a_{2j,3}{\mathrm{\ΔQ}}_{j+1}+a_{2j,4}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^M\right)=0,$

上式中， $a_{2j-\mathrm{1,1}}={\∂F}_j^C/{\∂Q}_j,$ $a_{2j-\mathrm{1,2}}={\∂F}_j^C/{\∂A}_j,$ $a_{2j-\mathrm{1,3}}={\∂F}_j^C/{\∂Q}_{j+1},$ $a_{2j-1,4}={\∂F}_j^C/{\∂A}_{j+1},$ $a_{2j,1}={\∂F}_j^M/{\∂Q}_j,$ $a_{2j,2}={\∂F}_j^M/{\∂A}_j,$ $a_{2j,3}={\∂F}_j^M/{\∂Q}_{j+1},$ $a_{2j,4}={\∂F}_j^M/{\∂A}_{j+1},$ 和分别表示连续和动量方程，表示余量，表示余量，下标j是断面编号，Δ表示两个连续Newton-Raphson迭代步间的增量；

3）为渠道中每一结构物设置上游边界节点G_i和下游边界节点H_i；4）给每一结构物上游边界节点的流量Q_Gi赋予初始值

5）将当前每一级结构物上游边界节点的流量Q_Gi作为其上游渠段的下边界条件，将结构物两侧流量相等Q_Gi＝Q_Hi作为其下游渠段的上边界条件，代入步骤3）获得的线性离散方程组的系数矩阵，对渠道中所有结构物的上下游边界水位进行同步求解，求解出每一级结构物的与上游边界节点流量Q_Gi对应的上下游边界节点水位Z_Gi和Z_Hi；

6）将每一结构物当前上游边界节点流量Q_Gi，以及与其对应的上下游边界节点水位Z_Gi和Z_Hi，代入已知的结构物流量关系式Q_i＝g′_i(Z_Gi，Z_Hi)，计算每一结构物的流量偏差ΔQ_i：

ΔQ_i=g′_i(Z_Gi，Z_Hi)-Q_Gi；

7）检查每一结构物的流量偏差是否满足ΔQ_i＝0：

如果所有结构物的流量偏差均满足ΔQ_i＝0，进入步骤11）；

如果有结构物的流量偏差未满足ΔQ_i＝0，进入步骤8）；

8）计算每一结构物的流量相对误差|ΔQ_i/Q_Gi|；

9）将每一结构物的流量相对误差|ΔQ_i/Q_Gi|分别与给定的误差阈值进行比较：

如果所有结构物的流量相对误差|ΔQ_i/Q_Gi|均小于等于给定的误差阈值，进入步骤11）；

如果有结构物的流量相对误差|ΔQ_i/Q_Gi|大于给定的误差阈值，进入步骤10）；

10）校正每一结构物的上游边界节点流量Q_Gi，返回步骤5）；

11）渠道中所有结构物的上游流量达到调控要求，结束。

上述步骤10）通过下式校正每一结构物上游边界节点流量Q_Gi：

Q_Gi＝Q_Gi+ΔQ_Gi，

ΔQ_Gi与ΔQ_i成线性关系，其比例系数为

$-1/\left[\mathrm{\α}(\frac{{\∂g}^{\′}}{{\∂Z}_G}\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}+\frac{{\∂g}^{\′}}{{\∂Z}_H}\frac{{\mathrm{df}}_H}{{\mathrm{dQ}}_H}-1)\right],$

和根据流量关系式Q_i＝g′_i(Z_Gi，Z_Hi)确定,且

$\frac{{\mathrm{df}}_{\mathrm{Gi}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}\≈\frac{1}{Q_{\mathrm{Gi}}{B}_{\mathrm{Gi}}/A_{\mathrm{Gi}}-\sqrt{{\mathrm{gA}}_{\mathrm{Gi}}{B}_{\mathrm{Gi}}}},$

$\frac{{\mathrm{df}}_{\mathrm{Hi}}}{{\mathrm{dQ}}_{\mathrm{Hi}}}\≈\frac{1}{Q_{\mathrm{Hi}}{B}_{\mathrm{Hi}}/A_{\mathrm{Hi}}+\sqrt{{\mathrm{gA}}_{\mathrm{Hi}}{B}_{\mathrm{Hi}}}},$

上式中，B_Gi和B_Hi为第i个结构物上下游边界节点的水面宽度，A_Gi和A_Hi为第i个结构物上下游边界节点的过水面积，g为重力加速度。

一种含多种结构物渠道的水动力条件并行化数值模拟方法，它包括以下步骤：

1）给定渠道中每一结构物处上侧水位Z_G=w_i(t)，下标i表征结构物序号，i≥1，采用Saint-Venant双曲方程组为整条渠道建立数学模型方程组；

2）采用Preissmann隐式差分格式离散步骤1）获得的渠道数学模型方程组，并通过Newton-Raphson迭代法获得对应的线性离散方程组：

$a_{2j-\mathrm{1,1}}{\mathrm{\ΔQ}}_j+a_{2j-\mathrm{1,2}}{\mathrm{\ΔA}}_j+a_{2j-\mathrm{1,3}}{\mathrm{\ΔQ}}_{j+1}+a_{2j-\mathrm{1,4}}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^C\right)=0,$

$a_{2j,1}{\mathrm{\ΔQ}}_j+a_{2j,2}{\mathrm{\ΔA}}_j+a_{2j,3}{\mathrm{\ΔQ}}_{j+1}+a_{2j,4}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^M\right)=0,$

上式中， $a_{2j-\mathrm{1,1}}={\∂F}_j^C/{\∂Q}_j,$ $a_{2j-\mathrm{1,2}}={\∂F}_j^C/{\∂A}_j,$ $a_{2j-\mathrm{1,3}}={\∂F}_j^C/{\∂Q}_{j+1},$ $a_{2j-1,4}={\∂F}_j^C/{\∂A}_{j+1},$ $a_{2j,1}={\∂F}_j^M/{\∂Q}_j,$ $a_{2j,2}={\∂F}_j^M/{\∂A}_j,$ $a_{2j,3}={\∂F}_j^M/{\∂Q}_{j+1},$ $a_{2j,4}={\∂F}_j^M/{\∂A}_{j+1},$ 和分别表示连续和动量方程，表示余量，表示余量，下标j是断面编号，Δ表示两个连续Newton-Raphson迭代步间的增量；

3）为渠道中每一结构物设置上游边界节点G_i和下游边界节点H_i；4）给每一结构物上游边界节点的流量Q_Gi赋予初始值

5）将当前每一级结构物上游边界节点的流量Q_Gi作为其上游渠段的下边界条件，将结构物两侧流量相等Q_Gi＝Q_Hi作为其下游渠段的上边界条件，代入步骤3）获得的线性离散方程组的系数矩阵，对渠道中所有结构物的上下游边界水位进行同步求解，求解出每一级结构物的与上游边界节点流量Q_Gi对应的上下游边界节点水位Z_Gi；

6）根据当前计算所得的上游边界节点水位Z_Gi和当前已知的上游边界节点水位w_i(t)计算每一结构物的水位偏差ΔZ_i：

ΔZ_i=Z_Gi-w_i(t)；

7）检查每一结构物的水位偏差是否满足ΔZ_i＝0：

如果所有结构物的水位偏差均满足ΔZ_i＝0，进入步骤11）；

如果有结构物的水位偏差未满足ΔZ_i＝0，进入步骤8）；

8）计算每一结构物的水位相对误差|ΔZ_i/w_i(t)|；

9）将每一结构物的水位相对误差|ΔZ_i/w_i(t)|分别与给定的误差阈值进行比较：

如果所有结构物的水位相对误差|ΔZ_i/w_i(t)|均小于等于给定的误差阈值，进入步骤11）；

如果有结构物的水位相对误差|ΔZ_i/w_i(t)|大于给定的误差阈值，进入步骤10）；

10）校正每一结构物的上游边界节点流量Q_Gi，返回步骤5）；

11）结构物的上游流量达到调控要求，校正结束。

上述步骤10）通过下式校正每一结构物上游边界节点流量Q_Gi：

Q_Gi＝Q_Gi+ΔQ_Gi，

ΔQ_Gi与ΔZ_i成线性关系，其比例系数为

$-1/\left[\mathrm{\α}\left(\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}\right)\right],$

且

$\frac{{\mathrm{df}}_{\mathrm{Gi}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}\≈\frac{1}{Q_{\mathrm{Gi}}{B}_{\mathrm{Gi}}/A_{\mathrm{Gi}}-\sqrt{{\mathrm{gA}}_{\mathrm{Gi}}{B}_{\mathrm{Gi}}}},$

上式中，B_Gi为第i个结构物上下游边界节点的水面宽度，A_Gi为第i个结构物上下游边界节点的过水面积，g为重力加速度。

本发明由于采取上述技术方案，具有以下优点：1、本发明对渠道每一结构物两侧的水体流动连接条件进行处理，即通过假设结构物处流量为结构物两侧渠段提供边界条件，并通过迭代校正使结构物处流量满足给定条件，从而为含多种结构物渠道的水动力条件求解提供了一种并行计算方法，该方法与传统的顺序求解渠道水动力条件的方法相比，耗时短，可以很大程度地提高运算效率。2、本发明通过在渠道结构物上游和下游分别补充相同的流量作为边界条件，为每一渠段建立一组封闭的方程组，且每一方程组中只含同一渠道的变量，使采用隐式差分法离散而成的系数矩阵仍然保持带状特征，从而使物理上相互联系的渠段在数值上相对独立，能够运用高效的带状矩阵求解法实现并行计算，提高了整体运算效率。3、本发明基于原传统的求解单一渠道水动力条件方法，根据渠道结构物两侧的水体流动连接条件改变计算节点的边界条件，即可实现对含多种结构物的渠道的水动力条件的并行求解，从而无需大范围地修改原单一渠道水动力条件求解程序，有效地提高了程序的通用性和可移植性。本发明可以广泛用于求解含多种结构物渠道在水量调控中的水动力条件。

附图说明

图1是本发明的原理示意图

具体实施方式

下面结合附图和实施例，对本发明进行详细的描述。

一般认为，渠道结构物处上下游流量相等，根据流量与上下游水位的关系可以分为以下两种类型：一种是流经结构物的流量由结构物上下游的流动情况，如上下游水位共同决定；另一种是流经结构物的流量只跟结构物上游水位有关。本发明提供的方法对上述两种类型的结构物都适用，它通过对两种类型的结构物的两侧水体流动连接条件进行处理，在求解单一渠道水动力条件方法的基础上，实现对含多种结构物的渠道的水动力条件的并行求解。

现有技术中，传统的求解单一渠道水动力条件的方法采用一维Saint-Venant双曲方程组建模，此方程组包括有一连续方程（1）和一动量方程（2）：

$\∂A/\∂t+\∂Q/\∂x-q=0---\left(1\right)$

$\∂Q/\∂t+\∂(Q^2/A)/\∂x+\mathrm{gA}\∂Z/\∂x+\mathrm{gA}{S}_f=0---\left(2\right)$

上式中，A为过水面积；Q为水流流量；Z是水位；t是时间坐标；x是空间坐标；q是单位河长侧向入流量；g是重力加速度；S_f是摩阻坡度，通常用曼宁公式确定。

上述Saint-Venant方程组采用Preissmann隐式差分格式离散化后，可以通过Newton-Raphson迭代法获得以下线性离散方程（3）和（4）：

$a_{2j-\mathrm{1,1}}{\mathrm{\ΔQ}}_j+a_{2j-\mathrm{1,2}}{\mathrm{\ΔA}}_j+a_{2j-\mathrm{1,3}}{\mathrm{\ΔQ}}_{j+1}+a_{2j-\mathrm{1,4}}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^C\right)=0---\left(3\right)$

$a_{2j,1}{\mathrm{\ΔQ}}_j+a_{2j,2}{\mathrm{\ΔA}}_j+a_{2j,3}{\mathrm{\ΔQ}}_{j+1}+a_{2j,4}{\mathrm{\ΔA}}_{j+1}+R\left(F_j^M\right)=0---\left(4\right)$

上式中， $a_{2j-\mathrm{1,1}}={\∂F}_j^C/{\∂Q}_j,$ $a_{2j-\mathrm{1,2}}={\∂F}_j^C/{\∂A}_j,$ $a_{2j-\mathrm{1,3}}={\∂F}_j^C/{\∂Q}_{j+1},$ $a_{2j-1,4}={\∂F}_j^C/{\∂A}_{j+1},$ $a_{2j,1}={\∂F}_j^M/{\∂Q}_j,$ $a_{2j,2}={\∂F}_j^M/{\∂A}_j,$ $a_{2j,3}={\∂F}_j^M/{\∂Q}_{j+1},$ $a_{2j,4}={\∂F}_j^M/{\∂A}_{j+1}.$

其中，和分别表示连续和动量方程，表示余量，表示余量，下标j是断面编号，Δ表示两个连续Newton-Raphson迭代步间的增量。上述方法中，每个计算节点可以是渠道上的一个外边界、一个汊点、一个划分断面或者一个结构物，m个计算节点（m>1）对应有m-1个计算微段和2m个待求变量，每个计算微段均包括两个形如式（3）和式（4）的线性离散方程，从而构成一个含有2（m-1）个方程和2m个待求变量的线性离散方程组。若要求解此线性离散方程组，就必须在计算微段两端补充外边界、汊点或结构物的连接条件，而当计算节点为结构物时，若直接离散则会破坏线性离散方程组系数矩阵的带状特征，而且求解时需要占用大量的计算机内存资源。

为解决上述问题，本发明对结构物两侧的水体流动连接条件进行了处理。如图1所示，字母A～F分别表示河网中渠道之间的节点，Q_D和Q_E渠道DC段和EF段入口处的流量，G、H是渠道中一结构物的上下游边界节点，该结构物可以是一闸门或者一水电站。CG段为结构物上游渠段，HB段为结构物下游渠段。Q_G和Z_G是上游边界节点G的流量和水位，Q_H和Z_H下游边界节点H的流量和水位。

对于第一种类型的结构物，流经结构物的流量由结构物的上下游水位共同决定，有流量关系式Q＝g′(Z_G，Z_H)。且由控制方程可知，上游边界节点G的水位Z_G和下游边界节点H的水位Z_H分别是其流量的函数，有流量-水位函数关系Z_G＝f_G(Q_G)和Z_H＝f_H(Q_H)。又由结构物两侧流量相等可知，当上游边界节点G的流量Q_G为一预先给定的初始值时，下游边界节点H的流量Q_H也应该为

将Q_G和Q_G＝Q_H分别作为CG渠段的下边界条件和HB渠段的上边界条件，采用Preissmann隐式差分格式离散Saint-Venant双曲方程组可以求解出与上游边界节点流量初始值对应的上游边界节点水位和下游边界节点水位由于和是根据设定的流量初始值，采用差分离散计算所得，因此当把它们代入g′(Z_G，Z_H)时，会发现由此确定的流量与预先设定的流量初始值并不相等，不满足流量关系式Q＝g′(Z_G，Z_H)，因此需要重新设定上游边界节点的流量Q_G，也即校正上游边界节点流量。

本发明采用Newton-Raphson迭代法校正上游边界节点流量，其包括以下步骤：

1)根据当前上游边界节点流量Q_G，以及与当前上游边界节点流量Q_G对应的上下游边界节点水位Z_G和Z_H，计算流量偏差ΔQ：

ΔQ=g′(Z_G，Z_H)-Q_G（5）

2）检查流量偏差ΔQ是否满足ΔQ＝0：

如果未满足，计算流量相对误差|ΔQ/Q_G|，进入步骤3）；

如果满足，进入步骤6）；

3）将流量相对误差|ΔQ/Q_G|与给定的误差阈值进行比较：

如果|ΔQ/Q_G|大于给定的误差阈值，进入步骤4）；

如果|ΔQ/Q_G|小于等于给定的误差阈值，进入步骤6）；

4）对上游边界节点流量Q_G进行校正：

Q_G＝Q_G+ΔQ_G

根据特征线理论可知，上式中ΔQ_G与ΔQ成线性关系，其线性关系由以下方法导出：

由牛顿下山法有以下关系式：

$\mathrm{\α}(\frac{{\∂g}^{\′}}{{\∂Z}_G}\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}+\frac{{\∂g}^{\′}}{{\∂Z}_H}\frac{{\mathrm{df}}_H}{{\mathrm{dQ}}_H}-1)\mathrm{\Δ}{Q}_G+R(g^{\′}-Q_G)=0---\left(6\right)$

其中，α为牛顿下山法系数因子，d为微分符号，R（）表示对括号内的结果取余值。

根据定义和式（5）可知：

R(g′-Q_G)＝g′-Q_G＝ΔQ

代入关系式（6），即有ΔQ_G与ΔQ成线性关系，且系数为：

$-1/\left[\mathrm{\α}(\frac{{\∂g}^{\′}}{{\∂Z}_G}\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}+\frac{{\∂g}^{\′}}{{\∂Z}_H}\frac{{\mathrm{df}}_H}{{\mathrm{dQ}}_H}-1)\right]$

其中，和可以根据流量关系式Q＝g′(Z_G，Z_H)确定，由特征线理论有：

$\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}\≈\frac{1}{Q_G{B}_G/A_G-\sqrt{{\mathrm{gA}}_G{B}_G}}$

$\frac{{\mathrm{df}}_H}{{\mathrm{dQ}}_H}\≈\frac{1}{Q_H{B}_H/A_H+\sqrt{{\mathrm{gA}}_H{B}_H}}$

其中，B_G和B_H为结构物上下游边界节点的水面宽度，A_G和A_H为结构物上下游边界节点的过水面积，g为重力加速度。

5）将校正后的Q_G作为结构物上游渠段的下边界条件，Q_G＝Q_H作为结构物下游渠段的上边界条件，进行新一轮的Newton-Raphson迭代，求解出新一轮的与校正后的Q_G对应的上下游边界节点水位Z_G和Z_H后,返回步骤1）；

6）结束校正。

对于第二种类型的结构物，流经结构物的流量只跟结构物的上游水位有关，计算水动力条件的方法与上述第一种类型结构物大体一致，但具体校正过程分以下两种情况：

I）已知结构物处流量-水位函数关系为Q_G＝g″(Z_G)，且有Z_G＝f_G(Q_G)，通过Newton-Rapson迭代法对上游边界节点流量Q_G进行校正，直至流量误差满足给定的误差阈值，该情况校正的具体步骤同第一种类型结构物，此处不再赘叙。

II）已知结构物上游边界节点水位为一事先给的固定值或随时间变化的固定函数，统一表示为w(t)，采用Newton-Raphson迭代法校正上游边界节点流量的过程如下：

1）根据当前计算所得的上游边界节点水位Z_G和当前已知的上游边界节点水位w(t)计算水位偏差ΔZ：

ΔZ=Z_G-w(t)；（7）

2）检查水位偏差ΔZ是否满足ΔZ＝0：

如果未满足，计算水位相对误差|ΔZ/w(t)|，进入步骤3）；

如果满足，进入步骤6）；

3）将水位相对误差|ΔZ/w(t)|与给定的误差阈值进行比较：

如果|ΔZ/w(t)|大于给定的误差阈值，进入步骤4）；

如果|ΔZ/w(t)|小于等于给定的误差阈值，进入步骤6）；

4）对当前上游边界节点流量Q_G进行校正：

Q_G＝Q_G+ΔQ_G

根据特征线理论可知，上式中ΔQ_G与ΔZ成线性关系，其线性关系由以下方法导出：

根据牛顿下山法有以下关系式：

$\mathrm{\α}\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}\mathrm{\Δ}{Q}_G+R(Z_G-w\left(t\right))=0---\left(8\right)$

根据定义和式（7）可知：

R(Z_G-w(t))＝Z_G-w(t)＝ΔZ

代入关系式（8），即有ΔQ_G与ΔZ成线性关系，且系数为：

$-1/\left[\mathrm{\α}\left(\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}\right)\right]$

由特征线理论有：

$\frac{{\mathrm{df}}_G}{{\mathrm{dQ}}_G}\≈\frac{1}{Q_G{B}_G/A_G-\sqrt{{\mathrm{gA}}_G{B}_G}}$

5）将校正后的Q_G作为结构物上游渠段的下边界条件，Q_G＝Q_H作为结构物下游渠段的上边界条件，进行新一轮的Newton-Raphson迭代，求解出新一轮的与校正后的Q_G对应的上下游边界节点水位Z_G后,返回步骤1）；

6）结构物的上游流量达到调控要求，校正结束。

上述方法中，无论是第一种类型或第二种类型的结构物，在求解其上下游边界节点流量的同时，还可以获得与之对应的上游过水面积A_G和下游过水面积A_H，因此流量的校正也可以根据过水面积偏差ΔA进行校正，ΔA＝A_H-A_G。由于步骤完全一致，所以此处不再详细介绍。

上述实施例是针对渠道中只包含一个结构物的情况进行描述，当渠道中包含有两个或两个以上的结构物时，仍然可以利用该方法对多个结构物进行同步求解。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (2)

1.一种含多种结构物渠道的水动力条件并行化数值模拟方法，它包括以下步骤：

1)给定渠道中每一结构物处流量-水位关系式Q_i＝g′_i(Z_Gi,Z_Hi)，下标i表征结构物序号，i≥1，采用Saint-Venant双曲方程组为整条渠道建立数学模型方程组；

2)采用Preissmann隐式差分格式离散步骤1)获得的渠道数学模型方程组，并通过Newton-Raphson迭代法获得对应的线性离散方程组：

$a_{2j-\mathrm{1,1}}\mathrm{\Δ}{Q}_j+a_{2j-\mathrm{1,2}}\mathrm{\Δ}{A}_j+a_{2j-\mathrm{1,3}}\mathrm{\Δ}{Q}_{j+1}+a_{2j-\mathrm{1,4}}\mathrm{\Δ}{A}_{j+1}+R\left(F_j^C\right)=0,$

$a_{2j,1}\mathrm{\Δ}{Q}_j+a_{2j,2}\mathrm{\Δ}{A}_j+a_{2j,3}\mathrm{\Δ}{Q}_{j+1}+a_{2j,4}\mathrm{\Δ}{A}_{j+1}+R\left(F_j^M\right)=0,$

上式中， $a_{2j-\mathrm{1,1}}=\∂F_j^C/\∂Q_j,a_{2j-\mathrm{1,2}}=\∂F_j^C/\∂A_j,a_{2j-\mathrm{1,3}}=\∂F_j^C/\∂Q_{j+1},a_{2j-\mathrm{1,4}}=\∂F_j^C/\∂A_{j+1},$

$a_{2j,1}=\∂F_j^M/\∂Q_j,a_{2j,2}=\∂F_j^M/\∂A_j,a_{2j,3}=\∂F_j^M/\∂Q_{j+1},a_{2j,4}=\∂F_j^M/\∂A_{j+1},$

和分别表示连续和动量方程，表示余量，表示余量，下标j是断面编号，Δ表示两个连续Newton-Raphson迭代步间的增量，A为过水面积；

3)为渠道中每一结构物设置上游边界节点G_i和下游边界节点H_i；

4)给每一结构物上游边界节点的流量Q_Gi赋予初始值

5)将当前每一级结构物上游边界节点的流量Q_Gi作为其上游渠段的下边界条件，将结构物两侧流量相等Q_Gi＝Q_Hi作为其下游渠段的上边界条件，代入步骤2)获得的线性离散方程组的系数矩阵，对渠道中所有结构物的上下游边界水位进行同步求解，求解出每一级结构物的与上游边界节点流量Q_Gi对应的上下游边界节点水位Z_Gi和Z_Hi；

6)将每一结构物当前上游边界节点流量Q_Gi，以及与其对应的上下游边界节点水位Z_Gi和Z_Hi，代入已知的结构物流量-水位关系式Q_i＝g′_i(Z_Gi,Z_Hi)，计算每一结构物的流量偏差ΔQ_i：

ΔQ_i＝g′_i(Z_Gi,Z_Hi)-Q_Gi；

7)检查每一结构物的流量偏差是否满足ΔQ_i＝0：

如果所有结构物的流量偏差均满足ΔQ_i＝0，进入步骤11)；

如果有结构物的流量偏差未满足ΔQ_i＝0，进入步骤8)；

8)计算每一结构物的流量相对误差|ΔQ_i/Q_Gi|；

9)将每一结构物的流量相对误差|ΔQ_i/Q_Gi|分别与给定的误差阈值进行比较：

如果所有结构物的流量相对误差|ΔQ_i/Q_Gi|均小于等于给定的误差阈值，进入步骤11)；

如果有结构物的流量相对误差|ΔQ_i/Q_Gi|大于给定的误差阈值，进入步骤10)；

10)校正每一结构物的上游边界节点流量Q_Gi，返回步骤5)；

所述步骤10)通过下式校正每一结构物上游边界节点流量Q_Gi：

Q_Gi＝Q_Gi+ΔQ_Gi，

ΔQ_Gi与ΔQ_i成线性关系，其比例系数为

$-1/\left[\mathrm{\α}(\frac{\∂g_i^{\′}}{\∂Z_{\mathrm{Gi}}}\frac{d{f}_{\mathrm{Gi}}}{d{Q}_{\mathrm{Gi}}}+\frac{\∂g_i^{\′}}{\∂Z_{\mathrm{Hi}}}\frac{d{f}_{\mathrm{Hi}}}{d{Q}_{\mathrm{Hi}}}-1)\right],$

和根据流量关系式Q_i＝g′_i(Z_Gi,Z_Hi)确定,且

$\frac{d{f}_{\mathrm{Gi}}}{d{Q}_{\mathrm{Gi}}}\≈\frac{1}{Q_{\mathrm{Gi}}{B}_{\mathrm{Gi}}/A_{\mathrm{Gi}}-\sqrt{g{A}_{\mathrm{Gi}}{B}_{\mathrm{Gi}}}},$

$\frac{d{f}_{\mathrm{Hi}}}{d{Q}_{\mathrm{Hi}}}\≈\frac{1}{Q_{\mathrm{Hi}}{B}_{\mathrm{Hi}}/A_{\mathrm{Hi}}+\sqrt{g{A}_{\mathrm{Hi}}{B}_{\mathrm{Hi}}}},$

上式中，B_Gi和B_Hi为第i个结构物上下游边界节点的水面宽度，A_Gi和A_Hi为第i个结构物上下游边界节点的过水面积，g为重力加速度，α为牛顿下山法系数因子；

11)渠道中所有结构物的上游流量达到调控要求，结束。

2.一种如权利要求1所述的含多种结构物渠道的水动力条件并行化数值模拟方法，它包括以下步骤：

1)给定渠道中每一结构物处上侧水位Z_G＝w_i(t)，下标i表征结构物序号，i≥1，采用Saint-Venant双曲方程组为整条渠道建立数学模型方程组；

2)采用Preissmann隐式差分格式离散步骤1)获得的渠道数学模型方程组，并通过Newton-Raphson迭代法获得对应的线性离散方程组：

$a_{2j-\mathrm{1,1}}\mathrm{\Δ}{Q}_j+a_{2j-\mathrm{1,2}}\mathrm{\Δ}{A}_j+a_{2j-\mathrm{1,3}}\mathrm{\Δ}{Q}_{j+1}+a_{2j-\mathrm{1,4}}\mathrm{\Δ}{A}_{j+1}+R\left(F_j^C\right)=0,$

$a_{2j,1}\mathrm{\Δ}{Q}_j+a_{2j,2}\mathrm{\Δ}{A}_j+a_{2j,3}\mathrm{\Δ}{Q}_{j+1}+a_{2j,4}\mathrm{\Δ}{A}_{j+1}+R\left(F_j^M\right)=0,$

上式中， $a_{2j-\mathrm{1,1}}=\∂F_j^C/\∂Q_j,a_{2j-\mathrm{1,2}}=\∂F_j^C/\∂A_j,a_{2j-\mathrm{1,3}}=\∂F_j^C/\∂Q_{j+1},a_{2j-\mathrm{1,4}}=\∂F_j^C/\∂A_{j+1},$

$a_{2j,1}=\∂F_j^M/\∂Q_j,a_{2j,2}=\∂F_j^M/\∂A_j,a_{2j,3}=\∂F_j^M/\∂Q_{j+1},a_{2j,4}=\∂F_j^M/\∂A_{j+1},$

和分别表示连续和动量方程，表示余量，表示余量，下标j是断面编号，Δ表示两个连续Newton-Raphson迭代步间的增量，A为过水面积；

3)为渠道中每一结构物设置上游边界节点G_i和下游边界节点H_i；

4)给每一结构物上游边界节点的流量Q_Gi赋予初始值

5)将当前每一级结构物上游边界节点的流量Q_Gi作为其上游渠段的下边界条件，将结构物两侧流量相等Q_Gi＝Q_Hi作为其下游渠段的上边界条件，代入步骤2)获得的线性离散方程组的系数矩阵，对渠道中所有结构物的上下游边界水位进行同步求解，求解出每一级结构物的与上游边界节点流量Q_Gi对应的上游边界节点水位Z_Gi；

6)根据当前计算所得的上游边界节点水位Z_Gi和当前已知的上游边界节点水位w_i(t)计算每一结构物的水位偏差ΔZ_i：

ΔZ_i＝Z_Gi-w_i(t)；

7)检查每一结构物的水位偏差是否满足ΔZ_i＝0：

如果所有结构物的水位偏差均满足ΔZ_i＝0，进入步骤11)；

如果有结构物的水位偏差未满足ΔZ_i＝0，进入步骤8)；

8)计算每一结构物的水位相对误差|ΔZ_i/w_i(t)|；

9)将每一结构物的水位相对误差|ΔZ_i/w_i(t)|分别与给定的误差阈值进行比较：

如果所有结构物的水位相对误差|ΔZ_i/w_i(t)|均小于等于给定的误差阈值，进入步骤11)；

如果有结构物的水位相对误差|ΔZ_i/w_i(t)|大于给定的误差阈值，进入步骤10)；

10)校正每一结构物的上游边界节点流量Q_Gi，返回步骤5)；

所述步骤10)通过下式校正每一结构物上游边界节点流量Q_Gi：

Q_Gi＝Q_Gi+ΔQ_Gi，

ΔQ_Gi与ΔZ_i成线性关系，其比例系数为

$-1/\left[\mathrm{\α}\left(\frac{d{f}_{\mathrm{Gi}}}{d{Q}_{\mathrm{Gi}}}\right)\right],$

且

$\frac{d{f}_{\mathrm{Gi}}}{d{Q}_{\mathrm{Gi}}}\≈\frac{1}{Q_{\mathrm{Gi}}{B}_{\mathrm{Gi}}/A_{\mathrm{Gi}}-\sqrt{g{A}_{\mathrm{Gi}}{B}_{\mathrm{Gi}}}},$

上式中，B_Gi为第i个结构物上下游边界节点的水面宽度，A_Gi为第i个结构物上下游边界节点的过水面积，g为重力加速度，α为牛顿下山法系数因子；

11)结构物的上游流量达到调控要求，校正结束。