CN112784502B - 一种水文水力学动态双向耦合的洪水预测方法 - Google Patents

Abstract

本发明公开一种水文水力学动态双向耦合的洪水预测方法，包括：A1，建立模型，模型包括动态边界，动态边界的下游为淹没区，上游为非淹没区，形成位于淹没区的水力单元和位于非淹没区的水文单元；A2，对水力单元用水力计算模型计算特征波数据，对水文单元用水文计算模型计算流场数据，根据特征波数据确定淹没区与非淹没区之间的质量和动量传递；A3，根据所述淹没区与非淹没区之间的质量和动量传递更新水深数据，并将模型中水深低于水深阈值的作为水文单元，水深高于水深阈值的作为水力单元，返回A2。本发明在计算中不固定淹没区与非淹没区之间的边界位置，边界推移形成新的动态边界，实现了两种区域的交替转换，模拟更接近真实情境。

Description

一种水文水力学动态双向耦合的洪水预测方法

技术领域

本发明涉及一种洪水模拟方法，尤其是一种水文水力学动态双向耦合的洪水预测方法，可以进行高效准确的洪水淹没模拟。

背景技术

随着城市化进程加快，城镇范围向周边地区快速扩展。为了便于获得水资源，新规划城镇选址往往选在靠近河道的平坦洼地带，城镇不透水面积的增加，导致地面积水量增加，加之城市排水体系设计时对下游河道高水位情况考虑欠缺，每次发生较强的降雨就会出现严重内涝。因此，科学解决洪水问题和内涝问题究具有很强的社会效益和科学意义。

洪涝预报是根据洪涝的形成和运动规律，基于已知的水文、土地、气象、水工建筑等条件，利用数学物理模型预测未来洪水的发生、变化、演进的综合性科学。洪涝预报模型主要分为水文预报模型和水力预报模型，两者的研究方向和计算模拟过程有各自的侧重点。水文预报模型以降雨作为输入驱动，根据流域内土壤、土地利用及河道地形等，分析计算不同频率降雨的产汇流情况，将径流沿河道演算到流域出口，输出结果一般是流域出口的流量变化过程。模型建立需要考虑水文循环的各个环节，对表征各环节的数学方程进行合理选择，计算所需参数种类繁多，包括水文循环各环节的土壤、气象、植被等多个方面，需要大量数据进行率定验证后才能进行实际应用。模型运行速度较快。而水力计算模型侧重于模拟洪水在河道及河漫滩的演进过程，输出结果通常是洪涝淹没范围的变化过程，以及淹没区域各处的水深和流速信息。模型建立的重点是利用数值计算方法对数学方程各项进行严密求解，计算格式通过严密的数学推导验证后模型方可使用。模型需要求解大量的数学物理方程，运行速度较慢。

流域范围较大，根据有水和无水特征，流域可以分为淹没区和非淹没区，如河道、水库、湖泊为淹没区，其余部分为非淹没区。求解淹没区的水流运动常采用水动力类模型，采用水文计算模型往往不能反映流动的动量传递，误差较大；对于范围较大的产汇流非淹没区，常采用水文计算模型确定产汇流水量，如果用水动力模型计算非淹没区的流动要素，往往耗费计算时间，数值计算容易不稳定。同时由于非淹没区的地理信息数据的不确定性和土地利用信息的不确定性，使用高精度的水动力模型求解，获得预期的高精度的结果也会受到各种质疑。

对于流域洪水问题单纯依赖水文计算模型或水动力模型，获得的结果应用范围有限。将水文-水动力模型进行有机耦合，是解决流域洪水及流域水环境问题的有效方法。水文- 水动力耦合模型中应用最成熟、广泛的耦合方法是“串联耦合”，也称“松散耦合”或“单向耦合”。串联式耦合模型虽然应用方便，但是对于真实的物理过程描述仍然不够准确。主要是将水文产汇流水量集中加载到河道上游，作为水动力模型的入流边界条件。除上游地点流量数值不真实外，流量边界条件位置的选取具有人为性。水力学边界条件输入的流量数值误差及边界位置点的误差，会给洪水预测结果带来误差，进而影响模拟结果的准确性。

单向耦合方法的误差来源之一是水文单元与水力单元的位置的固定不变，即水文单元的出流控制断面，也是水力单元的边界条件位置。在山谷、丘陵等地形起伏明显，河道汇流出口明确的情况下，单向耦合方式是可行的。但是在河漫滩附近、平原等地形，汇流进入洪水淹没区的位置并不确定，并且随着洪水涨落，淹没区与汇流区会交替变化。此时应用单向耦合方式，预先设定控制断面(汇流出口即淹没区边界)具有很大的人为性。控制点距离淹没区过远，那么淹没区的计算范围太大，降低计算效率，且进入淹没区的汇流量比实际情况较小，导致模拟得到的淹没范围变小；控制点距离淹没区过近，可能会出现洪水传播范围比模拟范围更大，计算过程中洪水淹没范围突破边界位置，导致模型计算失真。

发明内容

本发明根据水文计算模型和水力计算模型计算动态边界处的特征波传播情况，利用特征波分析判断淹没区与非淹没区边界处的质量(相当于流量)、动量(相当于流速)交换，遵循质量、动量守恒理论。

水文水力学动态双向耦合洪水预测模型(DBCM)首先将流域在空间上分区，分成洪水淹没区和非淹没区，并将两个区域动态边界位置变化范围定义为过渡区。不同区域采用不同的数学模型，不同区域采用不同的网格尺寸和不同的时间步长。在同一时间步内，两个区域通过动态边界连接并耦合，根据水深、特征波方向及水面衔接关系确定动态边界移动及界面水量、动量交换。本发明提出不同区域不同网格尺度的耦合方法，包括不同网格尺度的区域连接方法，不同区域不同时间步长的协调方法，动边界位置及界面计算方法。

本发明提出的DBCM数值稳定性好，计算效率高，适合应用于实际流域尺度的模型。本发明提到的双向耦合洪水预测模型包括非淹没区、淹没区和动态边界。

非淹没区：产流：降雨、蒸发、渗出、土壤储存水量等水量作为质量守恒方程的源项；汇流：利用2D扩散波方程(水文计算模型)，进行汇流计算。通过采用较大尺寸的粗网格和较大的时间步长，达到提高计算效率的目的。扩散波方程如下：

其中，h为水深，Ax和Ay为各自方向的流体面积，Sx和Sy为各自方向的坡度，n 为糙率，Rx和Ry为流动系数，Q_x和Q_y分别为各自方向的流量。

Q_m为降雨、蒸发、渗出、土壤储存水量等水量作为质量守恒方程的源项

q为单宽流量，即qx为QX/dx，qy为QY/d′y。

淹没区：产流：降雨、蒸发、渗出、土壤储存水量等水文要素产生的水量变化作为质量守恒方程的源项。洪水运动：2D运动波方程(水力计算模型)。淹没区受库朗稳定条件限制，采用较小的时间步长；局部区域的网格尺寸选取也较小，以反映以米为单位计量的堤防宽度尺寸。运动波方程如下：

U＝[h，hu，hv]^T

F＝[hu，hu²+gh²/2，huv]^T

G＝[hv，huv，hv²+gh²/2]^T

其中，h为水深，g为重力加速度，u为x方向的速度，v为y方向的速度，z为水位， C为谢才系数，S_X和S_Y为x和y方向的底坡源项和摩擦阻力项；

U为守恒物理量向量；

F为x方向对流项；

G为y方向对流项。

动态边界：以动态边界界面为边界，边界上游为水文产汇流(非淹没区)网格，下游为淹没区网格，根据界面上游、下游的流速、水深、特征波方向、水面线衔接理论确定动边界移位置变化和界面水量、动量交换。流域各处根据水深分布使用不同的计算模型，水深较浅的汇流区使用水文计算模型计算地表产汇流，水深较大的淹没区使用水力计算模型强调洪水的运动以及淹没情况。随着洪水涨落变化，两种区域的范围也处于动态变化中，同一位置在不同时刻适用的计算模型也需要动态调整。汇流进入洪水淹没区的边界位置也在不断变化。

本发明的技术方案如下：

一种水文水力学动态双向耦合的洪水预测方法，包括以下步骤：

A1，建立双向耦合洪水预测模型，所述双向耦合洪水预测模型包括动态边界s，动态边界s的下游为淹没区w，上游为非淹没区d，对双向耦合洪水预测模型划分网格，形成位于淹没区w的水力单元和位于非淹没区d的水文单元；

A2，对于水力单元采用水力计算模型计算特征波数据，对于水文单元采用水文计算模型计算流场数据，并根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递；

A3，根据所述淹没区w与非淹没区d之间的质量和动量传递更新水深数据，并将双向耦合洪水预测模型中水深低于水深阈值的作为水文单元，水深高于水深阈值的作为水力单元，返回A2。

可选地，步骤A2中，根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递所依据的传播条件包括：

第一判定条件：判断动态边界处特征波传播方向，如果动态边界处的特征波传播方向全部指向水力单元，则判定水文单元向水力单元传递质量，不传递动量；

第二判定条件：如果动态边界处的特征波传播方向不是全部指向水力单元，则判定水力单元向水文单元传递质量和动量。

可选地，所述水文计算模型采用扩散波方程进行水文计算，所述水力计算模型采用运动波方程进行水力计算，所述特征波数据包括流速、水深、特征波方向。

可选地，还设置过渡区ds表示动态边界s的位置变化范围，所述过渡区ds采用与淹没区w相同的网格划分形式。

可选地，还设置有网格变化线g，所述网格变化线g位于非淹没区d与过渡区ds之间，用于分隔不同的网格划分。

可选地，步骤A2中，所述对于水力单元采用水力计算模型计算特征波数据，对于水文单元采用水文计算模型计算流场数据，并根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递的过程包括：

对整个双向耦合洪水预测模型划分粗网格，并在淹没区w和过渡区ds进一步划分细网格，其中，h表示水深，u表示流速，

已知t_n＝nΔt时刻，淹没区w的水深和流速

非淹没区d的水深和流速

过渡区ds的水深和流速

动态边界s的位置sⁿ，动态边界s处的水深和流速

淹没区w的运动波方程时间步长为Δt，非淹没区d的扩散波方程的时间步长为 mΔt，且m＞1，根据以下步骤确定淹没区w与非淹没区d之间的质量和动量传递：

S1，在整个双向耦合洪水预测模型的粗网格上进行水文计算，获得淹没区w、非淹没区d的粗网格单元节点在(n+m)Δt时刻的水深和流速变量

其中，对于利用水文计算获得的淹没区的变量标记*；

S2，将粗网格的水深和流速变量在网格变化线g上按距离线性插值，获得在网格变化线g处的

作为边界条件，求解扩散波方程，获得过渡区ds以及淹没区w的细网格水深和流速变量

S3，在过渡区ds的细网格上，将(n+m)Δt时刻的过渡区ds的

在时间方向上线性插值到(n+k)Δt时刻的值：

S4，利用淹没区w的(n+k-1)Δt时刻的

过渡区ds的(n+k)Δt时刻的

确定移动后的(n+k)Δt时刻的动态边界信息s^n+k，进而确定(n+k)Δt时刻动态边界s上的

作为求解运动波方程的边界条件；

S5，利用(n+k)Δt时刻动态边界信息s^n+k以及动态边界s上的

仅在淹没区w求解运动波方程，求出淹没区w的(n+k)Δt时刻水深和流速

S6，利用淹没区w的水深

替换S4中的过渡区的水文计算获得的

返回S2，直到循环m次获得(n+m)Δt时刻的淹没区w的水深和流速

动态边界位置s^n+m及动态边界处的水深和流速

则返回S1，用淹没区w的

替代S1中利用水文计算获得的淹没区的水深

并继续执行。

可选地，在所述动态边界s处，根据两个相邻的网格的动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式。

可选地，在所述动态边界s处，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式包括以下一种情况：

对于水文单元k，i，水力单元j，i在k与j之间，k，i的水位坡度小于i，j的水位坡度，水文单元i的流速指向水力单元j，单元k的水深变化与单元i无关，单元k的流速变化由k在动态边界同一侧的水文单元分析获得；

首先统一用水力计算模型计算，获得各单元的流速和水深数据，单元i的速度采用当前时刻的流速，单元j的速度则采用前一时刻的流速，根据地表高程获得在动态边界处单元i和单元j的有效水深，根据有效水深以及流速进行特征波计算，如果单元i与单元j 符合所述传播条件的第一判定条件，则单元j的动量不传递到单元i，单元i仍采用由扩散波方程计算的流速，并且只有出流量，没有入流量，而单元j的水深变化除了自身根据水力计算模型计算得到的变化，还加上单元i传递的质量，但单元i，j的动态边界处没有动量传递。

可选地，在所述动态边界s上，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式包括以下一种情况：

水文单元k，i和水力单元j，并且i在k与j之间，k，i的水位坡度大于i，j的水位坡度，单元i的流速指向单元k，单元k的水深变化来自于单元i的传递的质量；

首先统一用水力计算模型计算，获得各单元的流速和水深数据，单元i的速度采用当前时刻的流速，单元j的速度则采用前一时刻的流速，对单元i，j在动态边界处进行特征波分析，获得单元i的有效水深和单元j的有效水深，根据有效水深以及流速进行特征波计算，如果单元i与单元j符合所述传播条件的第二判定条件，说明单元j的动量传递到单元i，则单元i纳入水力计算模型计算范围，单元i的水深变化量扣除当前流速向单元k输出的质量，以及与单元j通过求解水力计算方程得到的质量，而单元i的流速则根据当前速度与单元j通过求解水力计算方程得到流速增量并更新。

本发明的水文水力学动态双向耦合的洪水预测方法在计算中不固定淹没区与非淹没区之间的边界位置，汇流进入淹没区的地点是动态变化的，洪水传播也可以淹没原有的汇流点，边界推移形成新的动态边界，实现了两种区域的交替转换，模拟更接近真实情境。

本发明提出的动态双向耦合的洪水预测方法以通量求解分析为理论基础，根据水力计算模型计算单元动态边界处的特征波传播判断这两种模型的单元之间是否产生相互影响。如果动态边界处的特征波传播全部指向水力单元，那么水文单元向水力单元传递质量，不传递动量，两种单元应用各自原有的模型求解，动态边界保持不变或是向水力单元移动。如果动态边界处特征波传播不能全部指向水力单元，那么水力计算模型占据主导地位，所有单元以水力计算模型的计算方法求解，动态边界向水文单元移动。这种方法更加符合物理实际，动态边界实现了两种区域的交替转换，模拟更接近真实情境。利用特征波分析判断淹没区与非淹没区边界处的质量、动量交换，遵循质量、动量守恒理论。对比单向耦合方法要先后计算水文计算模型和水力计算模型，双向耦合方法中两种模型同时计算，模型计算时间可变，非淹没区范围较大时，主要采用水文计算模型计算水流运动，反之淹没区范围较大时，计算域内以水力计算模型为主，提高了整体计算效率。

附图说明

通过结合下面附图对其实施例进行描述，本发明的上述特征和技术优点将会变得更加清楚和容易理解。

图1是表示本发明实施例的水文水力学动态双向耦合的洪水预测方法的流程示意图；

图2是表示倾斜V形集水区的示意图；

图3是表示DBCM的流量过程线与其他现有模型的比较示意图；

图4是表示本发明实施例的流域分区示意图；

图5是表示本发明实施例的第一种情况对应的单元示意图；

图6是表示本发明实施例的第二种情况对应的单元示意图。

具体实施方式

下面将参考附图来描述本发明所述的实施例。本领域的普通技术人员可以认识到，在不偏离本发明的精神和范围的情况下，可以用各种不同的方式或其组合对所描述的实施例进行修正。因此，附图和描述在本质上是说明性的，而不是用于限制权利要求的保护范围。此外，在本说明书中，附图未按比例画出，并且相同的附图标记表示相同的部分。

首先，先通过模拟倾斜V形集水区上的二维地表流(Di Baldassarre等人，1996；Panday 和Huyakorn，2004)来验证一下水文和水动力模型的计算域是否可以动态切换，并比较 DBCM、UCM(单向耦合模型，例如MikeSHE/Mikell)和BCM(双向耦合模型)之间的差异。如图2所示，描绘了示例的地形。

计算区域呈对称的V形，中心区域有一对对称的山坡形成一条通道。河床坡度沿翼展方向S_0x为±0.05，与河道平行的流向S_0y，为0.02。山坡上的曼宁系数n为0.015，主河道的曼宁系数n为0.15。总模拟时间为180min，恒定降雨强度为10.8mmh-1，持续90min。图2中显示了V形集水区的详细尺寸和相关信息。

考虑到水力计算模型比水文计算模型更能描述地表径流，因此采用了相同降雨条件下的HM2D(HydroMPM2D，一款功能强大的水动力模型软件，属于双向耦合模型)和DBCM。当水深小于0.005m时，网格采用水文计算模型计算，当水深大于0.005m时，网格适用于水力计算模型。结果与由(Di Baldassarre等人，1996；Panday和Huyakorn，2004)开发的数值模型进行了比较。

如图3所示，将HM2D和DBCM得到的流量过程线与其他现有模型进行了比较。流量过程线与洪峰流量有很好的一致性。用HM2D和DBCM模拟的初始阶段与其他预测相一致。比较HM2D和DBCM的水文曲线可以看出，它们的上升段和洪峰流量非常吻合。因此，这两种方法都被用来模拟地面流动。HM2D与DBCM的差异在退行肢体逐渐显现。 HM2D在整个计算过程中使用水力计算模型(运动波方程)模拟水运动，而当上游水深低于阈值时，DBCM将从水力计算模型切换到水文计算模型(扩散波方程)。由于水文计算模型中没有时间偏导数项，目前的流速是当前水位梯度的函数，不等于前一时刻的速度加上通量项。因此，当DBCM从水力计算模型切换到水文计算模型时，速度计算方法也随之改变，因此，HM2D和DBCM之间的流量差就会出现。因此，在DBCM中，出口流量稍大，但随后稍小，以确保总体质量是守恒的。

参阅图2至6所示。流域洪水预测动态双向耦合模型将流域分成非淹没区和淹没区，两个区域之间以动态边界s连接。为了便于计算动态边界位置变化，在初始淹没区外围设置过渡区，在过渡区设置较细网格。过渡区是与其它两个区域重叠的区域，过渡区设置与淹没区一样的细网格，便于动态边界移动计算。

非淹没区为水文产汇流区，空间范围一般占流域总面积的80％以上，采用尺寸较大的粗网格。在非淹没区采用扩散波方程组进行产汇流计算，其连续方程右端源项为降雨扣除入渗、蒸发、土壤储存等因素的水量。因水深与时间关联，可采用显式求解；流速计算公式与时间无关，流速与水头坡降和河床糙率有关，可选用稳定性比较好的方法求解，因此非淹没区可以使用较大的时间步长。网格尺寸也相对较大，因此减少了网格数量。在非淹没区拟采用较大的时间步长、较大的网格尺寸，从而大幅提高计算效率。

三个区域的网格设置如图4所示：

1)非淹没区(标记d)，设置粗网格(粗网格可以覆盖淹没区)进行水文计算；

2)淹没区(标记w)，设置细网格进行水力计算；

3)过渡区(标记ds)，设置与淹没区同样的细网格用于动态边界位置计算。

两条内边界线：

1)动态边界(标记s)：位于过渡区，区分淹没区和非淹没区；

2)网格变化线(标记g)：位于非淹没区，用于水文粗网格与细网格变量过渡。

先将整个双向耦合洪水预测模型划分粗网格，并在淹没区w和过渡区ds进一步划分细网格，网格变化线g位于非淹没区与过渡区之间，其中，h表示水深，u表示流速，

已知t_n＝nΔt时刻，淹没区w的水深和流速

非淹没区d的水深和流速

过渡区ds的水深和流速

动态边界s的位置sⁿ，动态边界s处的水深和流速

淹没区w的运动波方程时间步长为Δt，非淹没区d的扩散波方程的时间步长为mΔt(m＞1)，根据以下步骤确定淹没区w与非淹没区d之间的质量和动量传递：

S1，在全双向耦合洪水预测模型的粗网格上进行水文计算，获得淹没区w、非淹没区 d的粗网格单元节点在(n+m)Δt时刻的水深和流速变量

其中，对于利用水文计算获得的淹没区的变量标记*；

S2，将粗网格的水深和流速变量在网格变化线g上插值，获得在网格变化线g处的

作为边界条件，求解扩散波方程，获得过渡区ds以及淹没区w的细网格水深和流速变量

S3，在过渡区ds的细网格上，将(n+m)Δt时刻的过渡区ds的

在时间方向上线性插值到(n+k)Δt时刻的值：

S4，利用淹没区w的(n+k-1)Δt时刻的

过渡区ds的(n+k)Δt时刻的

确定移动后的(n+k)Δt时刻的动态边界信息s^n+k，进而确定(n+k)Δt时刻动态边界s上的

作为求解运动波方程的边界条件；

S5，利用(n+k)Δt时刻动态边界信息s^n+k，仅在淹没区w求解运动波方程，求出淹没区w的(n+k)Δt时刻水深和流速

S6，利用淹没区w的水深

替换S4中的过渡区的水文计算获得的

返回S2，直到循环m次获得(n+m)Δt时刻的淹没区w的水深和流速

非淹没区d的水深和流速

动态边界位置s^n+m及动态边界处的水深和流速

则返回S1，用淹没区w的

替代S1中利用水文计算获得的淹没区的水深

并继续执行。

进一步地，步骤S4中，在所述动态边界s上，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的流量传输形式。根据黎曼问题的定义，所述动态边界处是指两个网格的相邻边即为动态边界处。

动态边界处的流量计算是由两种单元的流动状态所决定的，当两种单元的流动方向一致时，动态边界动态边界处的流量可能等于上游的坡面汇流量，完全由水文单元的汇流计算所决定，与单向耦合方法的计算一致。而当两种单元流动方向相反时，受下游洪水顶托作用，动态边界处动态边界的过流量可能由上游汇流和下游洪水流动共同决定。此时如果仍然将上游坡面汇流量直接赋予边界过流量，必然产生误差。因此，首先必须确定动态边界两侧两种单元的流动状态，进而确定动态边界处的流量。

根据Godunov理论(Godunov，1959)，采用有限体积法求解对流通量可以看作是一个局部Riemann问题。每个网格之间的不连续特征速度(特征波)代表局部流体变量在时间和空间上的传播。如图5、图6所示，单元k，i为水文单元，位于动态边界s的左侧，j 为水力单元，位于动态边界s的右侧，S_L代表左侧的单元产生的特征速度，S_R代表右侧的单元产生的特征速度。

其中，h_R代表左侧的有效水深；

h_L代表右侧的有效水深；

U_L表示左侧的单元的前一时刻的速度，

U_R表示右侧的单元的当前时刻的速度，

g为重力加速度。

其中，有效水深是区别于常规的普通水深的概念的，普通水深为水面高度减去底坡高度即为水深。而有效水深则为水面高度减去网格界面处两侧的底坡高度中最大的那个。

步骤A2中，根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递所依据的传播条件包括：

第一判定条件：判断动态边界处特征波传播方向，如果动态边界处的特征波传播方向全部指向水力单元，则判定水文单元向水力单元传递质量，不传递动量；

第二判定条件：如果动态边界处的特征波传播方向不是全部指向水力单元，则判定水力单元向水文单元传递质量和动量。

用特征速度来描述上述的传播条件是：第一判断条件：当特征速度S_L，S_R均为正时，流量完全取决于左侧的非淹没区的流动状态，当特征速度均为负时，流量也完全取决于左侧的非淹没区的流动状态。第二判断条件：当特征速度一个为负值一个为正值时，必须同时考虑两个网格的当前流动状态。将此方法应用于DBCM，可以指定边界处的计算格式。众所周知，水文计算模型只传递水质量，而水力计算模型同时传递水质量和动量。

本发明建立的动态双向耦合方法在两种模型的动态边界上，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，判断两种模型之间的流量传输。下面对动态边界的不同情况详细说明。

第一种情况：水文、水力单元各自独立计算，对应的实际情景如坡面向河道汇流，动态边界s只有水文计算模型计算的流量通过。如图5所示，其中单元k，i水深小，地表高程大，应用水文计算模型，单元j水深大，地表高程小，应用水力计算模型。首先对单元 k，i，j不作区分，统一应用扩散波方程的坡度分析，可知k，i的水位坡度小于i，j的水位坡度。根据扩散波方程的计算流速指向坡度最大的方向，水文单元i的流速指向水力单元j，因此在该步计算中，单元k的水深变化与单元i无关，流速变化则由k左侧的其他单元再行分析。

将单元i与单元j构成局部黎曼问题应用特征波分析。单元i的速度采用扩散波方程得到的当前时刻的流速，单元j的速度则采用前一时刻的流速。根据地表高程可获得在动态边界处单元i的有效水深和单元j的有效水深。虽然单元j的实际水深较大，但在动态边界处有效水深反而较小。根据有效水深和流速进行特征波计算(即上述的特征速度公式)，如果单元i与单元j符合的第一判断条件，即单元i的特征速度S_L与单元j的特征速度S_R同时为正或同时为负。说明单元j的动量不能传递到单元i。那么单元i仍采用由扩散波方程计算的流速，并且只有纯粹的出流量，没有入流量。而单元j的水深变化除了自身根据水力计算模型计算得到的变化，还加上单元i传递的水量。水文计算模型的动量方程没有对流项说明其动量较小，因此即使传递了质量但动量较小可以忽略，即单元i，j 间断处没有动量传递，两者流速没有产生相互影响。单元j的流速变化决定于本身与右侧其他单元是否产生通量传递，与单元i无关。

第二种情况：水文、水力单元统一计算，对应的实际情景如洪水传播，淹没区扩大，动态边界按照水力计算模型的方法计算通量。如图6所示，依旧是单元k，i水深小，应用水文计算模型，单元j水深大，地表高程小，应用水力计算模型。首先对单元k，i，j不作区分，统一应用扩散波方程的坡度分析，，可知单元i的流速指向单元k，单元k的水深变化来自于单元i的入流量。根据地表高程可获得单元i的有效水深和单元j的有效水深。单元i的速度采用当前时刻的流速，单元j的速度则采用前一时刻的流速，根据有效水深和流速进行特征波计算，如果单元i与单元j符合传播条件的第二判断条件，即单元i的特征速度S_L与单元j的特征速度SR一个为负值一个为正值时。说明单元j的动量可以传递到单元i。单元i纳入水力计算模型计算范围。单元i的水深变化量需要扣除当前流速向单元k输出的流量，以及与单元j通过求解运动波方程得到的流量。而单元i 的流速则以当前速度与单元j通过求解运动波方程得到流速增量并更新。

以上两种情况解释了在动态边界处的流量和流速计算，计算完成后根据设定水深阈值，判断更新后的两种计算模型区域，也即是动态边界的移动情况。

以上所述仅为本发明的优选实施例，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种水文水力学动态双向耦合的洪水预测方法，其特征在于，包括以下步骤：

A1，建立双向耦合洪水预测模型，所述双向耦合洪水预测模型包括动态边界s，动态边界s的下游为淹没区w，上游为非淹没区d，对双向耦合洪水预测模型划分网格，形成位于淹没区w的水力单元和位于非淹没区d的水文单元；

A2，对于水力单元采用水力计算模型计算特征波数据，对于水文单元采用水文计算模型计算流场数据，并根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递；

A3，根据所述淹没区w与非淹没区d之间的质量和动量传递更新水深数据，并将双向耦合洪水预测模型中水深低于水深阈值的作为水文单元，水深高于水深阈值的作为水力单元，返回A2，

步骤A2中，根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递所依据的传播条件包括：

第一判定条件：判断动态边界处特征波传播方向，如果动态边界处的特征波传播方向全部指向水力单元，则判定水文单元向水力单元传递质量，不传递动量；

第二判定条件：如果动态边界处的特征波传播方向不是全部指向水力单元，则判定水力单元向水文单元传递质量和动量，

在所述动态边界s处，根据两个相邻的网格的动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式，

在所述动态边界s处，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式包括以下一种情况：

对于水文单元k，i，水力单元j，i在k与j之间，k，i的水位坡度小于i，j的水位坡度，水文单元i的流速指向水力单元j，单元k的水深变化与单元i无关，单元k的流速变化由k在动态边界同一侧的水文单元分析获得；

首先统一用水力计算模型计算，获得各单元的流速和水深数据，单元i的速度采用当前时刻的流速，单元j的速度则采用前一时刻的流速，根据地表高程获得在动态边界处单元i和单元j的有效水深，根据有效水深以及流速进行特征波计算，如果单元i与单元j符合所述传播条件的第一判定条件，则单元j的动量不传递到单元i，单元i仍采用由扩散波方程计算的流速，并且只有出流量，没有入流量，而单元j的水深变化除了自身根据水力计算模型计算得到的变化，还加上单元i传递的质量，但单元i，j的动态边界处没有动量传递，

在所述动态边界s上，根据动态边界处的水位和地表高程，建立关于流速和有效水深的特征波分析，从而判断水文单元和水力单元之间的质量和动量传输形式包括以下一种情况：

水文单元k，i和水力单元j，并且i在k与j之间，k，i的水位坡度大于i，j的水位坡度，单元i的流速指向单元k，单元k的水深变化来自于单元i的传递的质量；

首先统一用水力计算模型计算，获得各单元的流速和水深数据，单元i的速度采用当前时刻的流速，单元j的速度则采用前一时刻的流速，对单元i，j在动态边界处进行特征波分析，获得单元i的有效水深和单元j的有效水深，根据有效水深以及流速进行特征波计算，如果单元i与单元j符合所述传播条件的第二判定条件，说明单元j的动量传递到单元i，则单元i纳入水力计算模型计算范围，单元i的水深变化量扣除当前流速向单元k输出的质量，以及与单元j通过求解水力计算方程得到的质量，而单元i的流速则根据当前速度与单元j通过求解水力计算方程得到流速增量并更新。

2.根据权利要求1所述的水文水力学动态双向耦合的洪水预测方法，其特征在于，

所述水文计算模型采用扩散波方程进行水文计算，所述水力计算模型采用运动波方程进行水力计算，

所述特征波数据包括流速、水深、特征波方向。

3.根据权利要求2所述的水文水力学动态双向耦合的洪水预测方法，其特征在于，

还设置过渡区ds表示动态边界s的位置变化范围，所述过渡区ds采用与淹没区w相同的网格划分形式。

4.根据权利要求3所述的水文水力学动态双向耦合的洪水预测方法，其特征在于，还设置有网格变化线g，所述网格变化线g位于非淹没区d与过渡区ds之间，用于分隔不同的网格划分。

5.根据权利要求4所述的水文水力学动态双向耦合的洪水预测方法，其特征在于，步骤A2中，所述对于水力单元采用水力计算模型计算特征波数据，对于水文单元采用水文计算模型计算流场数据，并根据所述特征波数据确定淹没区w与非淹没区d之间的质量和动量传递的过程包括：

对整个双向耦合洪水预测模型划分粗网格，并在淹没区w和过渡区ds进一步划分细网格，其中，h表示水深，u表示流速，

已知t_n＝nΔt时刻，淹没区w的水深和流速

非淹没区d的水深和流速

过渡区ds的水深和流速

动态边界s的位置sⁿ，动态边界s处的水深和流速

淹没区w的运动波方程时间步长为Δt，非淹没区d的扩散波方程的时间步长为mΔt，且m＞1，根据以下步骤确定淹没区w与非淹没区d之间的质量和动量传递：

S1，在整个双向耦合洪水预测模型的粗网格上进行水文计算，获得淹没区w、非淹没区d的粗网格单元节点在(n+m)Δt时刻的水深和流速变量

其中，对于利用水文计算获得的淹没区的变量标记*；

S2，将粗网格的水深和流速变量在网格变化线g上按距离线性插值，获得在网格变化线g处的

作为边界条件，求解扩散波方程，获得过渡区ds以及淹没区w的细网格水深和流速变量

S3，在过渡区ds的细网格上，将(n+m)Δt时刻的过渡区ds的

在时间方向上线性插值到(n+k)Δt时刻的值：

S4，利用淹没区w的(n+k-1)Δt时刻的

过渡区ds的(n+k)Δt时刻的

确定移动后的(n+k)Δt时刻的动态边界信息s^n+k，进而确定(n+k)Δt时刻动态边界s上的

作为求解运动波方程的边界条件；

S5，利用(n+k)Δt时刻动态边界信息s^n+k以及动态边界s上的

仅在淹没区w求解运动波方程，求出淹没区w的(n+k)Δt时刻水深和流速

S6，利用淹没区w的水深

替换S4中的过渡区的水文计算获得的

返回S2，直到循环m次获得(n+m)Δt时刻的淹没区w的水深和流速

动态边界位置s^n+m及动态边界处的水深和流速

则返回S1，用淹没区w的

替代S1中利用水文计算获得的淹没区的水深

并继续执行。

Images