CN109284531A - 一种基于空间拓扑的一二维水动力学耦合方法 - Google Patents

Abstract

本发明提出了一种基于空间拓扑的一二维水动力学耦合方法，包括：构建一维水动力模型，一维水动力模型包括：用于模拟河网的水流和涉水建筑物情况的一维河网模型、用于模拟城市管网的水流情况的一维管网模型；构建二维水动力模型，二维水动力模型用于对待分析的编制范围及控制线内的区域进行网格剖分，根据该区域的地形进行网格插值，并进行网格属性赋值；根据一维河网模型和所述二维水动力模型，执行针对河道与地面的地表一二维模型耦合；根据一维管网模型和所述二维水动力模型，执行针对城市地表地下的一二维模型耦合。本发明提供完善的耦合模式，适应不同洪水耦合计算需求；灵活的一二维耦合概化模式，能很好地解决小河流、街道行洪的概化处理。

Description

一种基于空间拓扑的一二维水动力学耦合方法

技术领域

本发明涉及水动力学分析技术领域，特别涉及一种基于空间拓扑的一二维水动力学耦合方法。

背景技术

作为洪水风险图编制的重要工具，洪水分析软件一直是国外商业软件占据主导地位。我国是一个水利大国，在水利领域的很多方面都取得了举世瞩目的成就，但是我们国内并没有形成一个自己的国产洪水分析软件品牌。山洪和城市洪涝目前仍是对人民生命财产威胁巨大的灾害事件，如果涉及自身的洪水分析方法实现对洪水的可靠分析是当前需要解决的技术问题。

在对洪水进行分析前，如何分别对河网、浅水、城市管网等模型，以及一二维耦合进行处理，提供一种可以适用于不同洪水耦合需求的方式，是当前需要解决的技术问题。

发明内容

本发明的目的旨在至少解决所述技术缺陷之一。

为此，本发明的目的在于提出一种基于空间拓扑的一二维水动力学耦合方法。

为了实现上述目的，本发明的实施例提供一种基于空间拓扑的一二维水动力学耦合方法，包括如下步骤：

步骤S1，构建一维水动力模型，所述一维水动力模型包括：用于模拟河网的水流和涉水建筑物情况的一维河网模型、用于模拟城市管网水流情况的一维管网模型；

步骤S2，构建二维水动力模型，所述二维水动力模型是指用于模拟地表二维浅水流动的二维模型；；

步骤S3，根据所述一维河网模型和所述二维水动力模型，执行针对河道与地面的地表一二维模型耦合；根据所述一维管网模型和所述二维水动力模型，执行针对城市地表地下的一、二维模型耦合。

进一步，在所述步骤S1中，构建所述一维河网模型，包括：

采用圣维南方程作为控制方程，形式如下：

q为旁侧入流，Q、A、B、Z分别为河道断面流量、过水面积、河宽和水位，V_X为旁侧入流流速在水流方向上的分量，一般可以近似为零，K为流量模数，反映河道的实际过流能力，α为动量校正系数，是反映河道断面流速分布均匀性的系数；

设置各类涉水建筑物作为联系要素，模拟涉水建筑物的过流水量；

进一步，在所述步骤S1中，构建所述一维管网模型，包括：

(1)所述一维管网控制方程包括连续方程和动量方程，所述连续方程为： Q为流量；A为过水断面面积；t为时间；x为距离；

所述动量方程为：

H为水头；g为重力加速度；S_f为摩阻坡度。

进一步，在所述步骤S3中，

所述地表一二维耦合可以分为两种连接：

侧向连接：水流从河道两岸流向二维区域或者从二维模型计算区域经由两岸流入河道，包括：

设某时刻河道与二维区域通过侧向交换的方式交换的流量为Q_l，采用堰流公式近似计算交换流量的方法如下：

式中：h_max和h_min分别采用下式计算：

h_max＝max(Z_r,Z_c)-Z_e

h_min＝min(Z_r,Z_c)-Z_e

式中：Z_r和Z_c为堰上下游水位，分别取河道和二维网格单元的水位值；Z_e为堰的高程，一般取河堤岸的高程；b_e为堰的宽度，一般取单元格与河道相连边的边长；

正向连接则：水流通过河道两端与二维计算区域进行水流交换，包括：

第一步：一维河网模型为二维模型提供流量边界，将河道与二维区域相连的那一端的断面流量作为边界条件提供给二维模型，即：

式中：为河道与二维区域连接断面的流量；M为二维区域与河道连接的单元边数目；l_k为单元边的边长；q_k为单元边的单宽流量；

第二步：根据给定的边界条件，二维模型从当前时间步更新至下一时间步；

第三步：二维模型为一维河网模型提供水位边界条件，根据二维模型更新后的单元值，将与河道相连单元的水位作为边界提供给一维模型作为边界条件：

式中：为河道下一时间步的水位边界条件；为更新后的单元水位值；L为垂向连接边界的总长度。

进一步，在所述步骤S3中，通过计算交换水量进行所述城市地表地下的一二维模型耦合，交换水量采用下列公式进行计算：

式中：H_surface为地面水头，H_node为排水管道水头，M为流量系数，H_g为地表高程。

根据本发明实施例的基于空间拓扑的一二维水动力学耦合方法，具有以下有益效果：

1D河网模型既能处理上千条河网+防洪工程调度控制，同时扩展有限体积方法河网计算引擎，能适应山区陡坡河道(网)模拟；

2D模型能够计算大的水面间断，能够捕捉激波；考虑孔隙率，使用大尺度网格考虑房屋等影响；

与SWMM城市排水***管网模型紧密集成，集成强大河网/管网一维显式计算引擎；

强大的网格剖分引擎，采用区域分解法Looping和铺路法Paving，直接生成不规则四边形网格；

完善的耦合模式(侧向、正向、垂向耦合模式)，适应不同洪水耦合计算需求；

灵活的一二维耦合概化模式，能很好地解决小河流、街道行洪的概化处理。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为根据本发明实施例的基于空间拓扑的一二维水动力学耦合方法的流程图；

图2为根据本发明实施例的流量边界条件的示意图；

图3为根据本发明实施例的水平方向水流交换的两种连接方式的示意图；

图4为根据本发明实施例的侧向连接的示意图；

图5为根据本发明实施例的正向连接的示意图。

具体实施方式

下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

如图1所示，本发明实施例的基于空间拓扑的一二维水动力学耦合方法，包括如下步骤：

步骤S1，构建一维水动力模型，一维水动力模型包括：用于模拟河网的水流和涉水建筑物情况的一维河网模型、用于模拟城市管网水流情况的一维管网模型。

在步骤S1中，构建一维河网模型，包括：

采用圣维南方程作为控制方程，形式如下：

q为旁侧入流，Q、A、B、Z分别为河道断面流量、过水面积、河宽和水位，V_X为旁侧入流流速在水流方向上的分量，一般可以近似为零，K为流量模数，反映河道的实际过流能力，α为动量校正系数，是反映河道断面流速分布均匀性的系数。当河道只有一个主槽时，α＝1.0，当河道有若干个主槽和滩地时，在主槽和滩地摩阻比降相等的假定下，可得n为主槽和滩地的分块个数，A_i、K_i为第i分块的过水面积与流量模数，A、K为断面总的过水面积与流量模数；所以α是断面位置及水位的函数，α值也像河道断面资料(河宽、过水面积一样)，可以先整理成α＝α(x，z)作为基本原始资料。对任一由断面i与断面i+1组成的河段，采用四点线性隐式差分格式进行数值离散，得任一河段的差分方程为：

以首节点水位和末节点水位为自由变量，采用三系数追赶法消去中间断面的水位和流量，最后得到首、末断面的流量与首、末节点水位关系的两个方程，即首、末断面流量表示成首、末节点水位的线性关系。

这两个方程形式如下：

其中：Z_(I)为首节点水位，Z_(J)为末节点水位，即首、末断面流量表达为首、末节点水位的线性组合。

依次由后向前把本断面流量表达成本断面水位和末节点水位的线性函数，递推公式如下：

Q_i＝α_i+β_iZ_i+ξ_iZ_(J) (4)

i＝L₂-2,L₂-3,...,L₁

同理从第一河段开始，设法把断面流量表达成本断面水位和首节点水位的线性函数：

Q_i＝θ_i+η_iZ_i+γ_iZ_(I) (5)

i＝L₁+2,L₁+3,...,L₂

因此，由上述递推公式可以得到式。在计算递推式时需要保存六个追赶系数α、β、ζ、θ、η和γ。一旦首、末节点水位求得后，利用式(4)和(5)对同一断面的流量有：

联立求解得：

求得Z_i后，代入到(3)式中即可得Q_i。

过水建筑物的水流模拟河道水流与湖泊内的水流运动，通过堰闸泵等工程筑物设施相连接，在本模型中形象地称其为“联系”；对于联系主要关心其过流流量的大小，根据过水建筑物的类别采用相应的水动力学方法模拟。具体地，设置各类涉水建筑物作为联系要素，模拟涉水建筑物的过流水量。

该特征单元为汇流型单元，主要包括闸、坝、水库、行蓄洪区口门等水工建筑物。该类型单元主要是影响水流的汇流过程，人类通过该类型单元来进行防洪调度、水资源调度。该单元的模拟模型主要是模拟其过水流量过程，下面以典型的宽顶堰为例说明。

宽顶堰上的水流可分为自由出流、淹没出流两种流态，不同流态采用不同的计算公式：

当出流为自由出流时：

当出流为淹没出流时：

式中：B为堰宽，Z_d为堰顶高程，Z_I为堰上节点水位，Z_J为堰下节点水位，H₀＝Z_I-Z_d，h_s＝Z_J-Z_d，m为自由出流系数，一般取0.325—0.385之间。为淹没出流系数，理论最大值为1.0，一般取小于1.0的数。

对自由出流流态，公式离散可得：

Q＝δ_Z1Z_I+β_Z1

对淹没出流流态，公式离散后得：

Q＝δ_Z2(Z_I-Z_J)

式中：δ_Z1、δ_Z2、β_Z1为与Z_I、Z_J有关的系数，一般常采用时段初水位来计算；有时为了提高计算精度，可采用迭代法计算δ_Z1、δ_Z2、β_Z1。

湖泊水流模拟在湖泊内不考虑其水流输运作用，只关心水位的高低，因而在其内只需满足水量平衡方程，采用零维模拟。采用水量平衡方程，零维模拟湖泊水流的水位，其中，水量平衡方程为：

A(z)为节点调蓄面积，∑Q为包括降雨产汇流、河道出入流在内的所有出入节点的流量，Z为水位。

工况控制条件的模拟在流域内部一般有若干工程设施组成的一套防洪控制体系，它们的运行均遵循一定的控制调度原则。这些控制调度使***的运行具有很大的复杂性，如多个工程启用时间、先后顺序均随着水流情况的不同而不同。在某些特殊情况下(如实时调度控制)，预先设定的调度原则难以满足实际要求，需要交互式地实时动态调整；这种复杂动态的控制条件，采用传统的方法无法实现模拟，需专门研制新的控制模拟方式——工程控制运行方式的数值模拟。

采用控制条件模拟水力工程建筑对水流运动的影响：

控制条件模型要素需要作用于堰、闸、泵等水利工程对水流运动进行控制影响。控制条件方式分为：增量控制和Gate控制两种方式。

一、增量控制方式：增量控制将水利工程(堰、闸、泵、口门)的模拟分成三个部分：工程的启用条件、工程开启过程、控制工程。

工程启用条件主要有：水位控制、流量控制、水位流量统计值控制等；

工程开启过程主要有：开启度增量、开启度相对值等；

控制工程：将工程启用条件和开启过程组合，添加到具体的工程形成一个完整的增量控制条件。

二、Gate控制方式：主要适用于平原河网地区，在同一位置闸泵共存，需要联合调度的情况。当满足条件泵站启用时，泵站内外会形成很大的水头差，闸门如果同时也是开启就会倒灌，这种情况下闸门必须关闭。采用Gate控制方式***会自动完成这一过程的模拟。

Gate控制方式采用控制条件表格组成决策树实现，一个控制条件可以同时实现对多个工程的控制。

边界模拟河网区域的边界主要有：区域产水、上游来水以及下游潮位过程，在模型里体现为流量或水位过程，用边界条件来模拟。对于河网水动力模型的边界条件主要有两类：水位边界条件和流量边界条件，由水文站和边界所在的河道断面构成。

一维河网、联系及零维之间的耦合，实际上是各单元交界面上的水量交换问题。反映水流运动的一个重要参数是水位，水位的高低可以直观地反映水流运行的情况，从公式(3)等可见，水位知道后相应的流量等其它水力要素均相应计算出来。从公式中可见河网、联系及零维调蓄的节点的水量平衡方程均相同，因此可将河网节点及零维调蓄单元的节点统称之为水位节点，其相应的水量平衡方程称为节点水位方程，将边界条件代入到相应的节点水位方程中可以得到节点水位线性完备的代数方程组，对节点水位方程采用直接或迭代解法解出所有节点的水位过程，然后回代求解出河道断面水位流量等水力要素。从中可看出建立节点水位方程是耦合模型的关键，节点水位求出后所有面上其它水力要素就很快解出。

此外，为了扩大河网计算引擎的适用范围，在一维河网模型中又引入了基于有限体积法的Godunov格式来处理流态过度的问题，该格式可以很好的处理水面大梯度流动和流态过度的情况。

2、构建一维管网模型，包括：

(1)建立控制方程

一维管网控制方程分为连续方程和动量方程：

连续方程：

式中：Q为流量，m3/s；A为过水断面面积，m2；t为时间，s；x为距离，m。

动量方程：

式中：H为水头，m；g为重力加速度，取9.8m/s2；S_f为摩阻坡度，由曼宁公式求得：

式中：K＝gn²，n为管道的曼宁系数；R为过水断面的水力半径，m；V为流速，绝对值表示摩擦阻力方向与水流方向相反，m/s。

假设v表示平均流速，将代入对流加速度项可得以下方程：

将Q＝Av代入连续方程，方程两边再同时乘以v，移项得方程：

将方程代入动量方程得方程：

忽略S₀项，将上述两个方程联立，依次求解各时段内每个管道的流量和每个节点的水头，有限差分格式如下：

式中：下标1和2分别表示管道或渠道的上下节点；L为管道长度，m。

求得Q_t+Δt：

式中：分别为t时刻的管道末端的加权平均值。

此外，为考虑管道的进出口水头损失，可以从H₂和H₁中减去水头损失。主要未知量为Q_t+Δt、H₂、 H₁、A₂、A₁，变量都与Q、H有关系。因此，还需要有Q和H有关的方程，可以从节点方程得到。

(2)建立节点控制方程，管网和渠道的节点控制方程为：

H为节点水头；Q_t为进出节点的流量；A_sk为节点的自由表面积；

化为有限差分格式为：

求得Δt时段内每个连接段的流量和每个节点的水头。

步骤S2，构建二维水动力模型，二维水动力模型是指用于模拟地表二维浅水流动的二维模型；。

二维洪水模拟模型***采用Godunov算法进行数值计算，其中Riemann问题采用Roe格式的近似 Riemann解进行计算，底坡源项采用特征分级离散，保证模型的守恒性，阻力源项采用隐式离散提高模型的稳定性，采用MUSCL空间重构和预测矫正法使得模型具有时间和空间二阶精度。

洪水演进的计算区域复杂，可能具有各种涉水构筑物，构筑物及其周边的水流不再符合浅水流动，因而无法采用浅水模型进行模拟计算，通常称其为内部边界条件。二维洪水模型对于内部边界条件的处理是其计算难点之一。涉水构筑物的过水能力多进行过大量的研究，通常具有一些成熟的理论或经验公式，本发明结合经验公式和模型的数值解法，通过通量概化计算方法，给出连接涉水构筑物的计算网格边的通量，该通量计算方法既能够保证模型的和谐稳定性，又能够精确计算通过建筑的流量通量和近似计算动量通量。

在步骤S2中，构建二维水动力模型，包括：

建立二维浅水方程，进行空间离散化处理，采用Godunov算法进行数值计算，方程中的底坡源项采用特征分级离散，阻力源项采用隐式离散，采用MUSCL空间重构和预测矫正法使得模型具有时间和空间二阶精度；

采用通量概化计算方法，给出连接涉水构筑物的计算网格边的通量。

水深平均的二维浅水方程可以简写为

式中，h为水深；u为x方向的流速；v为y方向的流速；s_x，s_y为源项，表达式为，

式中，p_a为水面大气压力；z_b为床面底高程；τ_ax，τ_ay为风载的作用力，表达式为，

式中，ρ_a为空气密度；为水面以上10m处的风速；C_Ds为拖曳系数。

c_x，c_y为地转科氏力，其在北半球表达式为，

c_x＝fv

c_y＝-fu

式中，f为科氏系数， 为地球的转动角速度， 为纬度。

τ_bx，τ_by为河底阻力，表达式为，

式中，n为糙率。

写为向量形式为，

U_t+E(U)_x+H(U)_y＝S₀+S

式中，(·)_t，(·)_x，(·)_y分别为对时间，空间平面x，y方向的偏导数，

在对上述的微分方程进行数值离散时，需要确定变量在计算网格中的位置。根据不同的变量位置的定义将计算网格称为Arakawa A-E网格。其中所有变量都定义在单元中心或节点的网格称为Arakawa A。在采用Godunov法离散时，变量定义在单元中心Arakawa A网格，也称作CC(Cell Center)网格，较为常用。

方程可以改写为

式中F＝(E,H)。将上述方程在单元V_i上积分

定义U_i为单元的平均值，存储在单元的中心，即

利用高斯定理把面积分转变为线积分，即

式中，ΔV_i为单元i的面积；为单元的边界；n＝(n_x,n_y)为单元边界的外法线方向，

式中，l_j(i,l)为边j的长度；F_n＝F·n＝En_x+Hn_y为通过第单元i的第j边的数值通量。目前有许多通量F_n的计算方法，也就构成了众多的数值格式。这里将重点介绍采用近似Riemann解计算通量F_n。

二维浅水方程有一个重要的性质，即为旋转不变性。利用这一性质将界面通量计算转换为求解一维 Riemann问题，即

式中，

利用了矩阵T以下性质

T^-1＝T^T

定义是一个局部坐标***，中心在边的中点，为边的外法向方向，为边的方向，即为正交的坐标系。在局部坐标系中的方向速度方向速度分别为，

在此局部坐标系下，二维齐次浅水方程变为一维问题，即为，

式中，为在局部坐标系中的变量。

通量的Jacobian矩阵A为

式中，为波速。

根据矩阵理论，矩阵A可以分解为

A＝RΛR^-1

其中，矩阵R＝[r¹,r²,r³]为

Λ为特征矩阵

其中，λ¹＝un_x+vn_y+cn_y；λ²＝un_x+vn_y；λ³＝un_x+vn_y-cn_y为特征值。特征值反映了特征变量的传播速度和方向，因而可以根据特征值的方向将特征矩阵分解为Λ＝Λ⁺+Λ^-，其中

相应的矩阵A可以分解为

A＝A⁺+A^-

A⁺＝RΛ⁺R^-1

A^-＝RΛ^-R^-1

数值通量为

变量在单元内常数或线性分布的近似下，在单元的边上就构成了一维的Riemann问题，即为，

通过求解该Riemann问题，计算出界面上状态变量值带入方程，就可以得到界面通量

近似Riemann解的计算，首先是将方程线性化

其中，为由和的某种平均值。这里介绍Roe格式近似Riemann解计算数值通量。

Roe格式要求具有所谓的U特性，即：

①相容性，即

②双曲性，和A一样有实数特征值；

③具有

满足U特性的流速和波速为，

将Roe平均的可以得到Roe格式数值通量

式中，

将沿右特征向量方向进行特征分解，

式中，其中

式中Δ()＝()_R-()_L。

得到界面数值通量的计算式

为获得和谐、稳定、守恒的计算结果，底坡源项S₀的离散非常重要，目前有很多的学者对源项的离散方法进行了研究。

底坡和孔隙率源项S₀是水深和底高程的函数，在界面上的间断为，

积分S₀可以得到

将沿着特征方向和进行分解

γ₂＝0

将根据特征矩阵Λ分解为作用于左侧单元和右侧单元的源项

式中，I单位对角矩阵；sign()为符号函数。

离散模型能否保持静水的静止状态，表明了模型是否具有C(conservation)特性。

在静止的水中

h+z_b＝η＝const

u＝v＝0

如果保持水流为静止状态，需要满足的条件之一是界面上的流量通量为0，即

式中，下标C表示通量中的连续方程分量。。

摩擦阻力源项等其它的源项对格式的稳定性也起着重要的作用，为增加格式的稳定性对除底坡源项以外的其它源项进行半隐式离散。用系数来衡量n+1时刻源项的影响系数，1-θ来衡量n时刻的影响系数，源项为

式中，Δt为时间步长，θ＝0时为完全显式，θ＝1时为完全隐式，令

离散为

在河道数值模拟时，通常上游给流量边界条件或单宽流量边界条件，下游给水位边界。如果单元i的第j条边为边界边。设边界上的水位为x,y方向流速为则边界通量为

由特征线理论，沿着正负特征线方向有特征不变量

即，沿特征线不变。R⁺为左特征不变量，R^-为右特征不变量。在边界局部坐标系中，计算区域根据特征线的方向，边界水流分为缓流、流出急流和流入急流三种情况。急流流入情况下，计算区域内信息影响不到边界，在边界上需要给定所有未知变量，即需要同时给定和急流流出情况下，边界外信息够影响不到计算区域内，因而不需要给定任何边界条件，该种情况下 和缓流情况下有

式中，为方向流速。空间一阶精度取单元的值，二阶精度取计算域内的插值。因而至少需要一个给定边界条件才能够计算出边界变量，以下边界条件的讨论都是基于边界缓流情况。

1、单宽流量边界条件

在上游边界给定单宽流量如果边界在计算域的左面有方程成立，将边界条件代入得

迭代求解方程可以获得边上水深h^*＝(c^*)²/g，再根据给定的单宽流量边界条件可以获得外法向流速如果给定的单宽流量为矢量则有否则边界切向流速可以近似的认为等于边界内的切线流速或

2、水位边界条件

在边界上给定水位h^*＝η(t)-z_b。如果边界在计算域的左面，求得外法向流速为

边界内的切线流速

3、流量边界条件

如图2所示，流量Q＝Q(t)分别从单元V₁，V₂，…，V_M的L₁，L₂，…，L_M边流入，流入的角度为θ。各个边界单元的水深和x，y方向的流速分别为{h₁,u₁,v₁}，{h₂,u₂,v₂}，…，{h_M,u_M,v_M}。

根据谢才公式和曼宁公式

式中，U为沿着水流方向流速；J为水力坡度。可见流速近似与水深h的2/3次幂成正比，所以有得到边界边的单宽流量

根据单宽流量边界条件可以计算出边界上的和

4、固壁边界条件

在固壁边界上，可以采用无滑移边界条件和滑移边界条件。无滑移边界条件边界取和滑移边界条件边界取和

在进行非恒定浅水模拟时，采用限制水深方法处理动边界：

限制水深法处理动边界时，把网格分为干、湿和半干三类，

干网格：在n时刻，若网格水深h＜htol1，且相邻单元的水深也小于htol 1，则没有流量和动量通过该公共边，如果所有相邻单元水深都小于htol1，则该单元为干单元；

半干网格：在n时刻，若网格水深htol1＜h＜htol2，且相邻单元水深小于htol2，则在相邻界面上只有流量的通量而没有动量通量；

湿网格：在n时刻，网格水深h＞htol 2；

其中，htol 1和htol2均为预设值。

采用区域分解法或铺路法进行网格剖分，

区域分解法为：在物体的边界上生成边界节点。通过连接两个边界节点，把这个区域剖分成两个子区域。然后在剖分线上增加新的节点。以递归的方式对每个子区域进行剖分，直到所有的子区域不可再分为止，即每个子区域包含六个或四个节点，对于每个包含六节点的子区域采用固定的模板进行闭合处理，生成最后的单元；

在本步骤中，采用区域分解法或铺路法实现对网格的划分。

首先，在待划分区域的边界上生成边界节点，通过连接两个边界节点，将该区域剖分成两个子区域。然后在剖分线上增加新的节点。以递归的方式对每个子区域进行剖分，直到所有的子区域不可再分为止，即每个子区域包含六个或四个节点，对于每个包含六节点的子区域采用固定的模板进行闭合处理，生成最后的单元。

需要注意的是，在区域分解法中，区域的外边界(只有一个)节点必须以逆时针方向给出，内部孔洞的边界(0个或多个)必须以顺时针方向给出。为了生成四边形网格，内外边界节点总数必须是偶数。因此，不仅在生成初始的边界节点时而且在剖分线上生成节点时都要满足这个要求。

在本步骤中，区域分解法主要包括以下步骤：内外边界的合并，边界单元的生成，剖分准则，剖分线上节点的生成，六节点闭合处理，光滑处理等。

1)内外边界的合并

在区域分解法中，要划分网格的区域必须是单连通域，这样区域和子区域都可以用一个连续的边界节点环来表示。对于内部有孔洞的区域，在网格划分之前，必须将多连通域转化成单连通域，即需要对内外边界进行合并。

将切割线上的节点和内部边界节点按一定的顺序***到外部边界节点环中，形成一个新的边界节点环。由于切割线上的节点(包括两端的节点)在新节点环中出现两次，因此只要最初的内外边界节点总数为偶数个，不管切割线上新生成的节点数目是偶数还是奇数，那么最后形成的新节点环中节点的数目肯定为偶数个。这样在切割上生成的节点数目不需要调整，即可满足边界节点总数为偶数的要求。但是为了避免最后生成的网格中出现只有一个单元的边连接内外边界的情况，当在切割线上生成的节点数目为0时，需要将节点数目调整为1。在切割线上生成节点的方法与在剖分线上生成节点的方法完全一致，具体方法见后面的“剖分线上节点的生成”。

2)边界单元的生成

为了提高边界单元的质量，使边界单元尽可能接近正方形，可以采用向内部偏置边界节点的方法在边界上生成一层边界单元。

在区域分解法中，将边界节点向内部偏移预设距离，形成新的边界节点，依次连接外部边界节点和偏置节点生成一层边界单元，使得边界单元的形状接近正方形，其中，所述预设距离为边界节点相邻两边长度的平均值。

具体来说，首先将边界节点向内部偏移一定的距离(等于该点相邻两边长度的平均值)，形成新的边界节点。依次连接外部边界节点和偏置节点生成一层边界单元，形状接近正方形。

根据边界节点的内角α大小，偏置节点的生成可分为以下三种情况：

(1)当120°≤α≤264°时，偏置节点沿边界节点内角平分线方向生成，长度等于该点相邻两边长度的平均值；

(2)当α＜120°时，该边界节点不生成偏置节点，且该点相邻的两边界节点生成的偏置节点应重叠在一起(可分别生成偏置节点，然后对其坐标求平均值)；

(3)当α＞264°时，该边界节点生成3个偏置节点，分别沿内角的1/3、1/2和2/3方向生成，长度等于该点相邻两边长度的平均值。

生成所有的偏置节点后，根据偏置节点的生成情况，按照下面的方法连接边界节点和偏置节点生成边界单元：

(1)对于上述第1种情况生成的偏置节点，连接该偏置节点、边界节点、下一个边界节点及其生成的偏置节点(或第一个偏置节点)，生成一个边界单元；若下一个边界节点属于上述第2种情况则不生成边界单元；

(2)对于上述第2种情况，连接边界节点、前后两个边界节点及其重叠的偏置节点，生成一个边界单元；

(3)对于上述第3种情况，连接边界节点及其生成的三个偏置节点，生成一个边界单元，然后连接边界节点、该点生成的第三个偏置节点、下一个边界节点及其生成的偏置节点(或第一个偏置节点)生成第二个边界单元；若下一个边界节点属于上述第2种情况则不生成第二个边界单元。

由偏移的节点组成的节点环在区域的内部又形成一个新的区域。在进行下一步操作之前，需要判断新边界(即偏移边界)与旧边界是否相交以及新边界自身是否相交。如果出现边界相交的情况，则取消生成的边界单元，不再生成边界单元，下一步仍对原先的边界节点环进行操作；若没有边界相交的情况，则下一步对新的边界节点环进行操作。

本发明提出一种确定最佳剖分线的新方法，在满足网格划分质量的前提下，显著地提高网格划分的效率。该方法确定最佳剖分线的步骤是：

(1)给两个角度临界变量赋初始值，让α＝120°，γ＝60°。

(2)如果节点i，j处的角度α_i、α_j有一个小于α，节点i，j之间就不能形成剖分线。

(3)如果连接节点i，j形成的四个剖分角γ_ij1，γ_ij2，γ_ji1，γ_ji2中有一个小于γ，节点i，j之间就不能形成剖分线。

(4)除了前面(2)、(3)两种情况外，连接节点i，j就可以确定一条剖分线，在所有与节点i相连的剖分线中，选择长度最短的剖分线作为一条候选的剖分线。候选的剖分线确定后，连接候选剖分线的节点就不再与其它节点配对确定剖分线。计算候选剖分线的长度和剖分后两子区域的面积比(≤1)。

(5)从剩余的节点中继续寻找其它可能的候选剖分线，并计算各剖分线的长度和剖分后两子区域的面积比。

(6)若寻找不到候选剖分线，则将临界角度变量α，γ减半，转到(2)处继续寻找。

(7)对于所有的候选剖分线，计算目标函数f＝w₁φ+w₂σ+w₃ι+w₄ω值，选择目标函数值最小且与区域边界不相交的剖分线作为最佳的剖分线。

目标函数f＝w₁φ+w₂σ+w₃ι+w₄ω包含四项φ,σ,ι,ω，为四个量化项(其值范围0～1)，通过加权平均(权重分别为w₁,w₂,w₃,w₄，w₁+w₂+w₃+w₄＝1)的方法得到一个综合的量化值。其中φ项用来衡量剖分线与两角(节点i，j内角)平分线之间的偏差程度，此项值越小，说明剖分线越接近两角平分线的位置。理想的剖分线应与两角的平分线重合。φ项的计算方法如下：

其中，

σ项衡量四个剖分角与π/2之间的偏差程度，此项值越小说明四个剖分角越接近直角，生成的网格越接近结构化网格。此项可以作为网格结构化指数。σ项的计算公式如下：

ι项为剖分线的长度与最长剖分线的比值，一般情况下剖分线的长度越短越好，ι项的计算公式如下：

ι＝l/l_max，

其中l_max为所有候选的剖分线中最长的剖分线。

ω项衡量剖分后两子区域的面积偏差程度，一般情况下，剖分后的两子区域面积大小相差越小越好。ω项的计算公式如下：

其中A为区域的面积，A₁,A₂为剖分两子区域的面积。

确定最佳剖分线后，就可以按照网格密度要求在剖分线上生成新节点，于是原区域一分为二，形成两个新的节点环。对于每个子区域，按照上述方法进行处理，直到所有子区域不可再分为止。

下面对剖分线上节点的生成进行说明。

首先介绍一下网格密度的概念。从直观上看，网格密度大的地方，网格单元尺寸较小；网格密度小的地方，网格单元尺寸较大。从数学上可以定义网格密度为网格单元长度的倒数。

区域(或子区域)边界节点处的网格密度值用该节点相邻两线段长度平均值的倒数来表示。假定剖分线上的网格密度为线性分布。设剖分线两端节点i，j的网格密度值为μ_i和μ_j，剖分线的长度为l_ij，则在剖分线上生成节点的数目另外，剖分后的子区域边界节点总数必须是偶数个，所以为了满足这个要求，必要时还需要对N值进行调整(加1或减1)。当剖分线直接连接最初的内外边界节点且N等于0时，为了避免最后生成的网格中出现只有一个单元的边连接内外边界的情况，需要将N值调整为2。

计算出在剖分线上生成的新节点数目后，下一步就是要确定这些节点在剖分线上的位置。在剖分线生成N个新节点，即将剖分线分割成N+1个线段。剖分线上新节点的位置的确定遵循各线段重量(各线段的平均密度值乘以该线段的长度)相等的原则。设与节点i相邻的第1个节点与节点i的距离为l_i1，则第1个节点的密度值可通过两端节点i，j的网格密度值线性插值得到，即

μ₁＝μ_i+(μ_j-μ_i)×l_i1/l_ij

节点i与第1个节点之间的线段重量为0.5×(μ₁+μ_i)×l_i1，等于剖分线重量的1/(N+1)，即 0.5×(μ₁+μ_i)×l_i1＝0.5×(μ_i+μ_j)×l_ij/(N+1)，

将μ₁＝μ_i+(μ_j-μ_i)×l_i1/l_ij带入上式，得到以l_i1为未知数的一元二次方程，求解该方程，并根据根计算得到第1个节点的网格密度值：

将上述计算得到的网格密度值带入0.5×(μ₁+μ_i)×l_i1＝0.5×(μ_i+μ_j)×l_ij/(N+1)式整理后得到以下等式：

R＝l_i1/l_ij＝(μ_i+μ_j)/(N+1)/(μ₁+μ_i)

设节点i，j的坐标分别为(x_i,y_i)和(x_j,y_j)，则第1个节点的坐标为：

x₁＝x_i+(x_j-x_i)×R，x₁＝x_i+(x_j-x_i)×R，

求出第1个节点的网格密度值和节点坐标之后，将第1个节点看作节点i，N值减1，采用上面的方法求出第2个节点的网格密度值和节点坐标，以此类推，求出剖分线上所有节点的网格密度值和节点坐标。

然后，在剖分线上增加新的节点，以递归的方式对每个子区域进行剖分，直到所有的子区域不可再分为止，即每个子区域包含六个或四个节点，对于每个包含六节点的子区域采用固定的模板进行闭合处理，生成最后的网格单元。

下面对六节点闭合处理过程进行说明。

每当剖分后的子区域包含六个节点，剖分过程结束，程序转入六节点闭合处理子程序，采用固定的模板进行模式匹配，生成最后的单元。进一步，在区域分解法中，在生成最后的网格单元之后，采用Laplace 光滑算法对最后生成的网格进行光滑处理。

(2)采用铺路法生成并优化网格，包括如下步骤：

(1)选择起始点

在待划分区域的多个边界中选择一个边界，并选择网格生成的起始点，其中，取边界上内角最小的节点为起始点；

(2)生成网格

在当前边界节点的基础上生成新的节点，组成新的单元，并更新当前边界；

首先介绍节点夹角的概念，，以节点xi与同层的相邻节点xi-1，xi+1所成的内角α，方向为顺时针。

根据α的大小，可将节点分为：

(1)终止节点：(2)边点，(3)角点：(4)转点：

其中δ取为5°＜δ＜10°。

新节点生成一般是由边界上相邻的三点为基础，设Ni为基点，其夹角为α，从Ni到Ni-1的距离为 d1，Ni到Ni+1的距离为d2，Ni到所生成的新节点的矢量为V。新的网格节点坐标是由上一层网格节点来决定的。

(3)缝合处理边界

在逐层生成网格的过程中，新边界上经常会出现细缝，即边界节点夹角过小的情况，有时甚至会出现负夹角。以上情况将会影响下一层网格生成质量。出现此种情况，***必须自动加以判断，并进行相应的缝合处理。由于初始边界节点为固定节点，无法进行缝合处理，因此缝合处理只用于内部边界节点。

对新生成的边界中相邻边界长短悬殊和节点夹角过小现象，进行边界的缝合处理，以利于网格继续生成。

(a)小角度边界节点的缝合处理

当内部边界节点Ni的夹角α小于π/6时，进行缝合处理，包括：将Ni前一节点Ni-1与其后一节点 Ni+1合并成一个节点Nj，合并后Nj的坐标为：减少细缝处节点的不规则程度，以生成规则的节点。当内部节点Ni的夹角太小时，将在该层网格中形成一条细缝。这种情况的处理是将Ni前一节点Ni-1与其后一节点Ni+1合并成一个节点Nj。合并后Nj的坐标为：

(b)过渡缝合处理

当网格单元的长边与短边的长短比大于2时，进行缝合处理，包括：在长边的中点增加一点，将该点与相邻点连接了，形成一个新的四边形单元。

(4)交叉处理

在新单元生成的过程中，由于边界的不规则，新单元往往会与旧边界发生交叉和重叠现象，中断网格的生成，进行交叉处理。当活动边界向区域内部推进时，往往容易与已生成的网格部分发生交叉或重叠，交叉多产生负内角，因此可在缝合过程中消除掉。然而多数情况下，重叠是由区域初始几何形状造成的。交叉和重叠问题，如果处理不好，不仅影响网格生成质量，而且还会导致网格生成过程无法继续。

(5)边界调整

当边界呈凸起状时，会使后续生成的单元越来越大。通过楔块***法，可改善网格尺寸大小，并控制其增大的趋势。

具体地，网格沿边界生成时，边界的几何形状将影响单元的形成。特别是当边界呈凸起状时，会使后续生成的单元越来越大。通过边界调整，可改善网格尺寸大小，并控制其增大的趋势。采用楔块***的方法来处理此类情况。

若当前边界上有连续3个或以上的节点的内角值大于183，则此边界需进行楔块***处理。楔块***的个数和位置，则根据满足条件的节点个数及其起点和终点的夹角来确定。当节点个数为3时，在中间节点处***一个楔块。当节点个数大于3时，根据起点和终点的夹角的大小，在每π/2角的中点处***一楔块。此算法只适用于内部边界节点。

(6)边界的光滑处理

对所有内部边界节点进行光滑处理，提高边界的光滑程度和后续网格的生成质量。

具体地，光滑处理是网格生成过程中最常用的处理方法。它不仅用于交叉的后处理，而且几乎网格的每次修改后都要进行光滑处理。光滑处理的目的是保证单元尺寸大小，边的垂直度，以及内部边界和全体网格的光滑要求。本节只讨论内部边界的光滑处理。

本算法设定所有不在内部边界上的节点为不动点，只对内部边界节点(不包括初始边界节点)进行光滑处理。内部边界的光滑处理采用修正的等参光滑方法。令Vi为从原点到边界节点Ni的矢量。假设边界节点Ni的周围有n个单元，Vmj，Vmk，Vml为在第m个单元中由原点分别指向Nj，Nk和Nl三个节点的向量，节点必须沿单元以顺时针或逆时针方向排列。由等参光滑后所得的新节点为Ni’，由原点到Ni’的矢量Vi’，可由以下公式得到，

令矢量ΔA代表节点Ni位置的变化，则

Δ_A＝V′_i-V_i

当边界节点周围单元数多于或少于两个时，其位置的变化可由上式求得。

而对于只与两个单元连接的边界节点来说，即n＝2时，其位置的修正既要包括长度光滑要求又要保证角度光滑要求。内部边界节点多为此类节点。首先进行长度光滑处理，如图2，令节点Nj为Ni的相对节点，Vi和Vj分别为从原点到节点Ni和Nj的矢量，lD为由节点生成方法求得的节点Ni的长度，Vij为从 Nj到N’i的矢量，lA为矢量Vij的长度。所以，节点Ni位置变化矢量ΔB为，

然后，进行角度的光滑处理。令节点Nj为Ni的相对节点，Ni-1和Ni+1分别为节点Ni的前后节点， Pi、Pi-1和Pi+1分别为从节点Nj到节点Ni、Ni-1和Ni+1的矢量。矢量PB1平分Pi-1和Pi+1的夹角。令PB2为调整后的节点Ni的新位置矢量，其端点在节点Nj，方向为PB1和Pi的角平分线的方向。接下来进行矢量PB2的长度的计算。首先求得矢量PB2与节点Ni-1和Ni+1连线的交点Q，设点Q到节点Nj的长度为lQ，节点Ni到节点Nj的初始长度为lD，矢量PB2的长度|PB2|由以下公式可得：

在角度光滑过程中，位置矢量的变化量为ΔC，可由如下公式得到

Δ_C＝P_B2-P_i

角度光滑处理有助于保持单元内角的垂直和促进内部边界的光滑。对位于内部边界上周围只有两个单元的节点Ni，其总的位置变化矢量为Δi，可由如下公式得到

尽管以上算法实现起来较复杂，但它对于段的正确形成至关重要，有力地保障了后续网格的生成。

(7)闭合处理

在网格生成的过程中，要进行边界节点总数的判断，当节点总数小于或等于六，则采用闭合处理方法生成最后的单元，并使该边界封闭。

(8)网格整合与光滑处理；

当网格生成完毕后，在局部总存在着一些缺陷，如内角不合格，不规则节点等等，这些缺陷需要通过内部网格的光滑与整合来消除。

网格生成完毕后，其拓扑结构(如节点─单元关系，单元─单元关系，节点─节点关系)也就确定了。这些关系在一定程度上反映着网格的规则程度。例如，在四边形网格中任一内部节点相连的最佳单元个数 NE为4，也就是说，当一个内部节点周围有四个单元，才有可能保证生成的网格无不规则单元。如果某一内部节点周围的单元过多或过少，这个节点周围的单元就有可能非常地不规则。在单元节点生成的公式中，只考虑了节点周围局部的情况，因而在内部多个边界会合时，往往会生成不规则节点。采用整合算法，通过改变网格的拓扑结构来使不规则单元的数量减少到最小。在整合过程中，边界节点为固定点，不做修正。

步骤S3，根据一维河网模型和二维水动力模型，执行针对河道与地面的地表一二维模型耦合；根据一维管网模型和二维水动力模型，执行针对城市地表地下的一二维模型耦合。

如图3所示，根据水流交换方式的不同，地表一二维耦合可以分为两种连接：侧向连接和正向连接。侧向连接是指水流从河道两岸流向二维区域或者从二维模型计算区域经由两岸流入河道，而正向连接则通常是指水流通过河道两端与二维计算区域进行水流交换。两种水流交换方式中，由于一维模型处理的水流和二维模型处理的水流的方向具有明显的差异，侧向连接时二维区域中水流方向在连接处一般与河道中水流形成一定的夹角，且不涉及河道边界条件问题，而正向连接时二维区域中水流方向与河道中水流方向在连接处通常是一致的，涉及到河道上游或者下游的边界条件问题，因而两种不同的连接方式通常需要采取不同的连接策略和计算方法。常用的水平连接计算方法包括堰流公式法，水量平衡法，黎曼问题法，互相提供边界法等。堰流公式法采用堰流公式计算交换流量，具有形式简单的优点，水量平衡法通过流入流出耦合区域的水量平衡关系来迭代求解耦合区域的水位，能够保证耦合区域水量守恒，黎曼问题法通过构建一维黎曼问题来求解数值通量，优势是可以较为方便地计算动量交换，互相提供边界法采用一二维模型互为对方提供边界条件的思想，方法的优点是无需额外的计算，比较容易实现，一般用来计算水平连接中的正向连接。本文采用较为简单，使用最为广泛的堰流公式法来计算侧向连接的水流交换问题，正向连接则采用互为提供边界条件的方式。

在步骤S3中，地表一二维耦合可以分为两种连接：

如图4所示，侧向连接：水流从河道两岸流向二维区域或者从二维模型计算区域经由两岸流入河道，包括：

设某时刻河道与二维区域通过侧向交换的方式交换的流量为Q_l，采用堰流公式近似计算交换流量的方法如下：

式中：h_max和h_min分别采用下式计算：

h_max＝max(Z_r,Z_c)-Z_e

h_min＝min(Z_r,Z_c)-Z_e

式中：Z_r和Z_c为堰上下游水位，分别取河道和二维网格单元的水位值；Z_e为堰的高程，一般取河堤岸的高程；b_e为堰的宽度，一般取单元格与河道相连边的边长；

在计算交换水量时，通常以二维区域中与河道相连的网格为单位分别进行计算。Z_c为与河道相连单元的水位，Z_r为单元格对应位置河道的水位，通过该单元处河道上下游断面水位进行插值得到，根据不同的水位组合，存在如下四种情况：(1)Z_c和Z_r均小于Z_e，则不计算交换水量，即Q_l＝0；(2)Z_c＞Z_r且max(Z_c,Z_r)＞Z_e，则水流从二维单元流向河道；(3)Z_c＜Z_r且max(Z_c,Z_r)＞Z_e，则水流从河道流向二维单元；(4)Z_c＝Z_r＞Z_e，此时依旧有水流交换，但是需要根据二维单元的流速方向来判定水流是流出还是流入。

如图5所示，正向连接则：水流通过河道两端与二维计算区域进行水流交换，包括：

第一步：一维河网模型为二维模型提供流量边界，将河道与二维区域相连的那一端的断面流量作为边界条件提供给二维模型，即：

式中：为河道与二维区域连接断面的流量；M为二维区域与河道连接的单元边数目；l_k为单元边的边长；q_k为单元边的单宽流量；

第二步：根据给定的边界条件，二维模型从当前时间步更新至下一时间步；

第三步：二维模型为一维河网模型提供水位边界条件，根据二维模型更新后的单元值，将与河道相连单元的水位作为边界提供给一维模型作为边界条件：

式中：为河道下一时间步的水位边界条件；为更新后的单元水位值；L为垂向连接边界的总长度。

在步骤S3中，通过计算交换水量进行城市地表地下的一二维模型耦合，交换水量采用下列公式进行计算：

式中：H_surface为地面水头，H_node为排水管道水头，M为流量系数，H_g为地表高程。

城市地表地下一二维模型耦合算法

目前地表洪水模型与地下管网模型的耦合通常做法是计算交换水量，然后代入各自模型中计算，更新到下一步，交换水量采用下列公式进行计算：

式中：H_surface为地面水头，H_node为排水管道水头，M为流量系数，H_g为地表高程。

交换水量计算后，需要作两方面的校核：(1)由于一个网格单元可能对应很多个管网节点，需要以二维网格单元为单位校核拟交换水量是否超过单元现有总水量，即可能出现二维网格单元中总水量不够，无法满足当前与众多管网节点计算的交换流量，出现这种情况，需要按比例减少交换水量；(2)由于交换水量是根据当前步结果显式计算的，未考虑下一时段网格单元以及管道的来水量，可能会出现交换流量过大的情况。若上一时间步二维网格与管网节点之间水流交换方向为网格单元流入管网节点，而上一步计算完成后，该节点出现溢流，说明上一步交换水量过多，需要当前步中将网格单元的水量增加该溢流值，以满足水量平衡。

根据本发明实施例的基于空间拓扑的一二维水动力学耦合方法，具有以下有益效果：

1D河网模型既能处理上千条河网+防洪工程调度控制，同时扩展有限体积方法河网计算引擎，能适应山区陡坡河道(网)模拟；

2D模型能够计算大的水面间断，能够捕捉激波；考虑孔隙率，使用大尺度网格考虑房屋等影响；

与SWMM城市排水***管网模型紧密集成，集成强大河网/管网一维显式计算引擎；

强大的网格剖分引擎，采用区域分解法Looping和铺路法Paving，直接生成不规则四边形网格；

完善的耦合模式(侧向、正向、垂向耦合模式)，适应不同洪水耦合计算需求；

灵活的一二维耦合概化模式，能很好地解决小河流、街道行洪的概化处理。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。本发明的范围由所附权利要求及其等同限定。

Claims (5)

1.一种基于空间拓扑的一二维水动力学耦合方法，其特征在于，包括如下步骤：

步骤S1，构建一维水动力模型，所述一维水动力模型包括：用于模拟河网的水流和涉水建筑物情况的一维河网模型、用于模拟城市管网水流情况的一维管网模型；

步骤S2，构建二维水动力模型，所述二维水动力模型是指用于模拟地表二维浅水流动的二维模型；

步骤S3，根据所述一维河网模型和所述二维水动力模型，执行针对河道与地面的地表一二维模型耦合；根据所述一维管网模型和所述二维水动力模型，执行针对城市地表地下的一、二维模型耦合。

2.如权利要求1所述的基于空间拓扑的一二维水动力学耦合方法，其特征在于，在所述步骤S1中，构建所述一维河网模型，包括：

采用圣维南方程作为控制方程，形式如下：

q为旁侧入流，Q、A、B、Z分别为河道断面流量、过水面积、河宽和水位，V_X为旁侧入流流速在水流方向上的分量，一般可以近似为零，K为流量模数，反映河道的实际过流能力，α为动量校正系数，是反映河道断面流速分布均匀性的系数；

设置各类涉水建筑物作为联系要素，模拟涉水建筑物的过流水量。

3.如权利要求1所述的基于空间拓扑的一二维水动力学耦合方法，其特征在于，在所述步骤S1中，构建所述一维管网模型，包括：

(1)所述一维管网控制方程包括连续方程和动量方程，所述连续方程为：

Q为流量；A为过水断面面积；t为时间；x为距离；

所述动量方程为：

H为水头；g为重力加速度；S_f为摩阻坡度。

4.如权利要求1所述的基于空间拓扑的一二维水动力学耦合方法，其特征在于，在所述步骤S3中，

所述地表一二维耦合可以分为两种连接：

侧向连接：水流从河道两岸流向二维区域或者从二维模型计算区域经由两岸流入河道，包括：

设某时刻河道与二维区域通过侧向交换的方式交换的流量为Q_l，采用堰流公式近似计算交换流量的方法如下：

式中：h_max和h_min分别采用下式计算：

h_max＝max(Z_r,Z_c)-Z_e

h_min＝min(Z_r,Z_c)-Z_e

式中：Z_r和Z_c为堰上下游水位，分别取河道和二维网格单元的水位值；Z_e为堰的高程，一般取河堤岸的高程；b_e为堰的宽度，一般取单元格与河道相连边的边长；

正向连接则：水流通过河道两端与二维计算区域进行水流交换，包括：

第一步：一维河网模型为二维模型提供流量边界，将河道与二维区域相连的那一端的断面流量作为边界条件提供给二维模型，即：

式中：为河道与二维区域连接断面的流量；M为二维区域与河道连接的单元边数目；l_k为单元边的边长；q_k为单元边的单宽流量；

第二步：根据给定的边界条件，二维模型从当前时间步更新至下一时间步；

第三步：二维模型为一维河网模型提供水位边界条件，根据二维模型更新后的单元值，将与河道相连单元的水位作为边界提供给一维模型作为边界条件：

式中：为河道下一时间步的水位边界条件；为更新后的单元水位值；L为垂向连接边界的总长度。

5.如权利要求1所述的基于空间拓扑的一二维水动力学耦合方法，其特征在于，在所述步骤S3中，通过计算交换水量进行所述城市地表地下的一二维模型耦合，交换水量采用下列公式进行计算：

式中：H_surface为地面水头，H_node为排水管道水头，M为流量系数，H_g为地表高程。