首页> 中国专利> 一种含多种结构物渠道的水动力条件并行化数值模拟方法

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

摘要

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

著录项

  • 公开/公告号CN102890732A

    专利类型发明专利

  • 公开/公告日2013-01-23

    原文格式PDF

  • 申请/专利权人 清华大学;

    申请/专利号CN201210333345.X

  • 发明设计人 陈永灿;朱德军;俞茜;刘昭伟;

    申请日2012-09-10

  • 分类号G06F17/50;

  • 代理机构北京纪凯知识产权代理有限公司;

  • 代理人徐宁

  • 地址 100084 北京市海淀区北京市100084信箱82分箱清华大学专利办公室

  • 入库时间 2024-02-19 16:49:45

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2018-08-28

    未缴年费专利权终止 IPC(主分类):G06F17/50 授权公告日:20150422 终止日期:20170910 申请日:20120910

    专利权的终止

  • 2015-04-22

    授权

    授权

  • 2013-03-06

    实质审查的生效 IPC(主分类):G06F17/50 申请日:20120910

    实质审查的生效

  • 2013-01-23

    公开

    公开

说明书

技术领域

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

背景技术

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

发明内容

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

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

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

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

a2j-1,1ΔQj+a2j-1,2ΔAj+a2j-1,3ΔQj+1+a2j-1,4ΔAj+1+R(FjC)=0,

a2j,1ΔQj+a2j,2ΔAj+a2j,3ΔQj+1+a2j,4ΔAj+1+R(FjM)=0,

上式中,a2j-1,1=FjC/Qj,a2j-1,2=FjC/Aj,a2j-1,3=FjC/Qj+1,a2j-1,4=FjC/Aj+1,a2j,1=FjM/Qj,a2j,2=FjM/Aj,a2j,3=FjM/Qj+1,a2j,4=FjM/Aj+1,和分别表示连续和动量方程,表示余量,表示余量,下 标j是断面编号,Δ表示两个连续Newton-Raphson迭代步间的增量;

3)为渠道中每一结构物设置上游边界节点Gi和下游边界节点Hi; 4)给每一结构物上游边界节点的流量QGi赋予初始值

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

6)将每一结构物当前上游边界节点流量QGi,以及与其对应的上下游边界节点 水位ZGi和ZHi,代入已知的结构物流量关系式Qi=g′i(ZGi,ZHi),计算每一结构物的流量 偏差ΔQi

ΔQi=g′i(ZGi,ZHi)-QGi

7)检查每一结构物的流量偏差是否满足ΔQi=0:

如果所有结构物的流量偏差均满足ΔQi=0,进入步骤11);

如果有结构物的流量偏差未满足ΔQi=0,进入步骤8);

8)计算每一结构物的流量相对误差|ΔQi/QGi|;

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

如果所有结构物的流量相对误差|ΔQi/QGi|均小于等于给定的误差阈值,进入步 骤11);

如果有结构物的流量相对误差|ΔQi/QGi|大于给定的误差阈值,进入步骤10);

10)校正每一结构物的上游边界节点流量QGi,返回步骤5);

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

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

QGi=QGi+ΔQGi

ΔQGi与ΔQi成线性关系,其比例系数为

-1/[α(gZGdfGdQG+gZHdfHdQH-1)],

和根据流量关系式Qi=g′i(ZGi,ZHi)确定,且

dfGidQGi1QGiBGi/AGi-gAGiBGi,

dfHidQHi1QHiBHi/AHi+gAHiBHi,

上式中,BGi和BHi为第i个结构物上下游边界节点的水面宽度,AGi和AHi为第i个 结构物上下游边界节点的过水面积,g为重力加速度。

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

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

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

a2j-1,1ΔQj+a2j-1,2ΔAj+a2j-1,3ΔQj+1+a2j-1,4ΔAj+1+R(FjC)=0,

a2j,1ΔQj+a2j,2ΔAj+a2j,3ΔQj+1+a2j,4ΔAj+1+R(FjM)=0,

上式中,a2j-1,1=FjC/Qj,a2j-1,2=FjC/Aj,a2j-1,3=FjC/Qj+1,a2j-1,4=FjC/Aj+1,a2j,1=FjM/Qj,a2j,2=FjM/Aj,a2j,3=FjM/Qj+1,a2j,4=FjM/Aj+1,和分别表示连续和动量方程,表示余量,表示余量,下 标j是断面编号,Δ表示两个连续Newton-Raphson迭代步间的增量;

3)为渠道中每一结构物设置上游边界节点Gi和下游边界节点Hi; 4)给每一结构物上游边界节点的流量QGi赋予初始值

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

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

ΔZi=ZGi-wi(t);

7)检查每一结构物的水位偏差是否满足ΔZi=0:

如果所有结构物的水位偏差均满足ΔZi=0,进入步骤11);

如果有结构物的水位偏差未满足ΔZi=0,进入步骤8);

8)计算每一结构物的水位相对误差|ΔZi/wi(t)|;

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

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

如果有结构物的水位相对误差|ΔZi/wi(t)|大于给定的误差阈值,进入步骤10);

10)校正每一结构物的上游边界节点流量QGi,返回步骤5);

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

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

QGi=QGi+ΔQGi

ΔQGi与ΔZi成线性关系,其比例系数为

-1/[α(dfGdQG)],

dfGidQGi1QGiBGi/AGi-gAGiBGi,

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

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

附图说明

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

具体实施方式

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

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

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

A/t+Q/x-q=0---(1)

Q/t+(Q2/A)/x+gAZ/x+gASf=0---(2)

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

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

a2j-1,1ΔQj+a2j-1,2ΔAj+a2j-1,3ΔQj+1+a2j-1,4ΔAj+1+R(FjC)=0---(3)

a2j,1ΔQj+a2j,2ΔAj+a2j,3ΔQj+1+a2j,4ΔAj+1+R(FjM)=0---(4)

上式中,a2j-1,1=FjC/Qj,a2j-1,2=FjC/Aj,a2j-1,3=FjC/Qj+1,a2j-1,4=FjC/Aj+1,a2j,1=FjM/Qj,a2j,2=FjM/Aj,a2j,3=FjM/Qj+1,a2j,4=FjM/Aj+1.

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

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

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

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

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

1)根据当前上游边界节点流量QG,以及与当前上游边界节点流量QG对应的上 下游边界节点水位ZG和ZH,计算流量偏差ΔQ:

ΔQ=g′(ZG,ZH)-QG(5)

2)检查流量偏差ΔQ是否满足ΔQ=0:

如果未满足,计算流量相对误差|ΔQ/QG|,进入步骤3);

如果满足,进入步骤6);

3)将流量相对误差|ΔQ/QG|与给定的误差阈值进行比较:

如果|ΔQ/QG|大于给定的误差阈值,进入步骤4);

如果|ΔQ/QG|小于等于给定的误差阈值,进入步骤6);

4)对上游边界节点流量QG进行校正:

QG=QG+ΔQG

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

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

α(gZGdfGdQG+gZHdfHdQH-1)ΔQG+R(g-QG)=0---(6)

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

根据定义和式(5)可知:

R(g′-QG)=g′-QG=ΔQ

代入关系式(6),即有ΔQG与ΔQ成线性关系,且系数为:

-1/[α(gZGdfGdQG+gZHdfHdQH-1)]

其中,和可以根据流量关系式Q=g′(ZG,ZH)确定,由特征线理 论有:

dfGdQG1QGBG/AG-gAGBG

dfHdQH1QHBH/AH+gAHBH

其中,BG和BH为结构物上下游边界节点的水面宽度,AG和AH为结构物上下游 边界节点的过水面积,g为重力加速度。

5)将校正后的QG作为结构物上游渠段的下边界条件,QG=QH作为结构物下 游渠段的上边界条件,进行新一轮的Newton-Raphson迭代,求解出新一轮的与校 正后的QG对应的上下游边界节点水位ZG和ZH后,返回步骤1);

6)结束校正。

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

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

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

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

ΔZ=ZG-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)对当前上游边界节点流量QG进行校正:

QG=QG+ΔQG

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

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

αdfGdQGΔQG+R(ZG-w(t))=0---(8)

根据定义和式(7)可知:

R(ZG-w(t))=ZG-w(t)=ΔZ

代入关系式(8),即有ΔQG与ΔZ成线性关系,且系数为:

-1/[α(dfGdQG)]

由特征线理论有:

dfGdQG1QGBG/AG-gAGBG

5)将校正后的QG作为结构物上游渠段的下边界条件,QG=QH作为结构物下 游渠段的上边界条件,进行新一轮的Newton-Raphson迭代,求解出新一轮的与校 正后的QG对应的上下游边界节点水位ZG后,返回步骤1);

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

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

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

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

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号