CN109255185A - 一种基于城市地表地下管网的一二维水动力学耦合分析方法 - Google Patents

Abstract

本发明提出了一种基于城市地表地下管网的一二维水动力学耦合分析方法，包括：城市管网模型的建立，城市管网模型又包括产流计算模型、一维管网模型和一二维耦合模型。输入地表水数据和地下水数据，根据所述地表水数据和地下水数据建立地表地下连接关系；生成并优化网格单元；分别计算地表径流和排水管网水流；计算地表径流和排水管网水流的交换水量；对计算得到的交换水量进行校核，以实现地表二维地下管网水动力学耦合。本发明实现对城市暴雨内涝进行地表地下水流耦合分析计算。

Description

一种基于城市地表地下管网的一二维水动力学耦合分析方法

技术领域

本发明涉及水动力学分析技术领域，特别涉及一种基于城市地表地下管网的一二维水动力学耦合分析方法。

背景技术

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

虽然，现有技术中公开了基于一维非恒定数值模型的洪水分析方法等方式，但尚未有提出一种基于城市地表地下管网的一二维水动力学耦合分析方法。

发明内容

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

为此，本发明的目的在于提出一种基于城市地表地下管网的一二维水动力学耦合分析方法。

为了实现上述目的，本发明的实施例提供一种基于城市地表地下管网的一二维水动 力学耦合分析方法，其中，建立城市地表地下管网模型，城市地表地下管网模型包括产流计算模型、一维管网模型和一二维耦合模型，

建立产流计算模型的方法为：将计算区域划分成若干个子汇水区，每一个子汇水区 概化为一个非线性蓄水池，其入流项有降水和来自上游子流域的流出量，只有当蓄水池水 深超过最大洼地蓄水量时地表出流才会发生，其大小通过曼宁公式计算得出：

式中：Q为出流量，m³/s；W为子流域宽度，m；S为坡度；n为曼宁粗率系数；

建立一维管网模型的方法为：计算管渠的汇流，通过控制方程计算连接管渠的连续 性和动量平衡性，通过节点控制方程计算节点的出水量，

建立一二维耦合模型，包括以下步骤：

步骤S1，输入地表水数据和地下水数据，根据所述地表水数据和地下水数据建立地表 地下连接关系；

步骤S2，生成并优化网格单元；

步骤S3，分别计算地表径流和排水管网水流；

步骤S4，计算地表径流和排水管网水流的交换水量；

步骤S5，对计算得到的交换水量进行校核，以实现地表二维地下管网水动力学耦合。

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述产 流计算模型中的每个子流域产流包括三种不同类型的地表产流：

无洼蓄不透水地表上的降雨损失主要为蒸发，产流量表示为：

R₁＝P-E

式中：R₁为无洼不透水地表的产流量；P为降雨量；E为蒸发量，

有洼蓄不透水地表上的降雨损失主要为填洼和蒸发，产流量表示为：

R₂＝P-D-E

式中：R₂为有洼不透水地表的产流量；D为洼蓄量，

透水地表的降雨损失主要包括洼蓄和下渗，产流量表示为：

R₃＝(i-f-E)×Δt

式中：R₃为透水地表的产流量；i为降雨强度；f为入渗强度；Δt为时间间隔。

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述建 立一维管网模型中，控制方程又包括连续方程和动量方程，

连续方程：

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

动量方程：

式中：H为水深；g为重力加速度；S_f为摩阻坡度，

节点控制方程：

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

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述步 骤S5中，对所述交换水量进行校核，包括：

(1)由于一个网格单元对应多个管网节点，以二维网格单元为单位校核拟交换水量是 否超过单元现有总水量，当出现二维网格单元中总水量不够，无法满足当前与众多管网节 点计算的交换流量时，需要按比例减少交换水量；

(2)当上一时间步二维网格与管网节点之间水流交换方向为网格单元流入管网节点， 而上一步计算完成后，该节点出现溢流，说明上一步交换水量过多，需要当前步中将网格 单元的水量增加该溢流值，以满足水量平衡。

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述步 骤S2中，采用区域分解法生成并优化网格，包括如下步骤：

在待划分区域的边界上生成边界节点，通过连接两个边界节点，将该区域剖分成两个 子区域；

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

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述步 骤S2中，采用铺路法生成并优化网格，包括如下步骤：

(1)选择起始点

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

(2)生成网格

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

(3)缝合处理边界

对新生成的边界中相邻边界长短悬殊和节点夹角过小现象，进行边界的缝合处理；

(4)交叉处理

在新单元生成的过程中，由于边界的不规则，新单元往往会与旧边界发生交叉和重叠 现象，中断网格的生成，进行交叉处理；

(5)边界调整

当边界呈凸起状时，通过楔块插入法，改善网格尺寸大小，对边界进行调整；

(6)边界的光滑处理

对所有内部边界节点进行光滑处理；

(7)闭合处理

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

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

当网格生成完毕后，对内部网格进行光滑和整合处理，以消除存在的缺陷。

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述铺 路法中，所述缝合处理边界，包括：

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

当内部边界节点Ni的夹角α小于π/6时，进行缝合处理，包括：将Ni前一节点Ni-1与其后一节点Ni+1合并成一个节点Nj，合并后Nj的坐标为：

(2)过渡缝合处理

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

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述铺 路法中，所述交叉处理，包括：

当判断出现网格交叉重叠情况时，采用将重叠的部分连接的方法来修正网格，连接将 导致形成两个新的内部边界，为确保连接有效同时避免生成不规则节点，进行了以下处理：

(1)偶数限制

保证连接后的边界节点总数均为偶数

(2)连续缝合

交叉边界连接完毕后，进行缝合处理，保证生成的新边界的质量，以利于后续网格的 生成。

本发明一种基于城市地表地下管网的一二维水动力学耦合分析方法，其中，所述铺 路法中，所述闭合处理，包括如下步骤：

当边界包括四个节点时，则将其作为一个单元处理；当边界包括六个节点时，则根据 其中包含的终止点个数进行闭合处理。

根据本发明实施例的基于城市地表地下管网的一二维水动力学耦合分析方法，对产 流计算模型、一维管网模型和一二维耦合模型进行单独建立，对每个模型进行单独的数据 处理，从而完成对城市地表地下管网的一二维水动力学耦合的分析。

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

附图说明

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

图1为根据本发明实施例的地表产流计算模型概念图；

图2为根据本发明实施例的基于城市地表地下管网的一二维耦合模型建立方法的流 程图；

图3为根据本发明实施例的网络生成过程的示意图；

图4为根据本发明实施例的内外边界合并的示意图；

图5为根据本发明实施例的内部含有四个孔洞的区域内外边界合并的示意图；

图6为根据本发明一个实施例的边界单元生成的示意图；

图7为根据本发明另一个实施例的边界单元生成的示意图；

图8为根据本发明实施例的网格剖分线的示意图；

图9为根据本发明实施例的六节点闭合处理模版的示意图；

图10(a)和图10(b)为根据本发明实施例的光滑处理的示意图；

图11为根据本发明实施例的铺路法的实施流程图；

图12为根据本发明实施例的节点夹角的示意图；

图13为根据本发明一个实施例的网格生成算法的示意图；

图14为根据本发明另一个实施例的网格生成算法的示意图；

图15为根据本发明再一个实施例的网格生成算法的示意图；

图16为根据本发明实施例的小角度缝合的示意图；

图17为根据本发明实施例的过渡缝合处理的示意图；

图18为根据本发明实施例的交叉判断的示意图；

图19(a)和图19(b)为根据本发明实施例的交叉处理的示意图；

图20为根据本发明实施例的楔块插入的示意图；

图21为根据本发明实施例的长度光滑处理的示意图；

图22为根据本发明实施例的角度光滑处理的示意图；

图23为根据本发明实施例的节点消去的示意图；

图24(a)和图24(b)为根据本发明实施例的网格整合的示意图；

图25为根据本发明实施例的对角线替换的示意图；

图26为根据本发明实施例的边消去的示意图；

图27为根据本发明实施例的网格自动生成的流程图；

图28为根据本发明实施例的约束线处理的示意图；

图29为根据本发明实施例的约束线与外边界的合并示意图；

图30为根据本发明实施例的多区域边界的确定示意图；

图31为根据本发明实施例的计算区域离散化的示意图；

图32为根据本发明实施例的栅格单元的示意图。

具体实施方式

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

本发明实施例的一种基于城市地表地下管网的一二维水动力学耦合分析方法，建立 城市地表地下管网模型，城市地表地下管网模型包括产流计算模型、一维管网模型和一二维耦合模型。

产流计算模型的建立：将计算区域划分成若干个子汇水区或者排水区，将每个子汇水 区或者排水区概化成如图1所示的概念模型。图1中d为蓄水池水深，d_p为蓄水池最大洼 蓄深，Q为地表出流。

每一个子流域概化为一个非线性蓄水池，其入流项有降水和来自上游子流域的流出量； 流出项包括蒸散发、下渗和出流量。蓄水池的容量为最大洼地蓄水量，蓄水池中的水深由 子流域的水量平衡计算得出，并且随着时间不断更新。只有当蓄水池水深超过最大洼地蓄 水量时地表出流才会发生，其大小通过曼宁公式计算得出：

式中：Q为出流量，m³/s；W为子流域宽度，m；S为坡度；n为曼宁粗率系数。

地表产流是指降雨经过损失变成净雨的过程。每个子流域产流由3部分组成：透水面 积上的产流不仅要扣除洼蓄量，还要扣除下渗和蒸散发引起的初损；有洼蓄不透水面积上 的产流等于其降雨量减去蒸散发和洼蓄量；无洼蓄不透水面积上的产流等于其降雨量减去 蒸发损失。三种类型地表单独进行产流计算，子流域出流量等于三个部分出流量之和。

无洼蓄不透水地表上的降雨损失主要为蒸发，产流量表示为：

R₁＝P-E

式中：R₁为无洼不透水地表的产流量，mm；P为降雨量，mm；E为蒸发量，mm。

有洼蓄不透水地表上的降雨损失主要为填洼和蒸发，产流量表示为：

R₂＝P-D-E

式中：R₂为有洼不透水地表的产流量，mm；D为洼蓄量，mm。

透水地表的降雨损失主要包括洼蓄和下渗，产流量表示为：

R₃＝(i-f-E)×Δt

式中：R₃为透水地表的产流量，mm；i为降雨强度，mm/h；f为入渗强度，mm/h；Δt为时间间隔，h。

一维管网模型的建立：一维管网模型提供三种方法用于管渠的汇流计算，即恒定流法、 运动波法和动力波法。恒定流法假定在每一个计算时段流动都是恒定、均匀的，是最简单 的汇流计算方法。运动波法可以模拟管渠中水流的空间和时间变化，但是仍然不能考虑回 水、入口及出口损失、逆流和有压流动。动力波法按照求解完整的圣维南方程组来进行汇 流计算，是最准确同时也是最复杂的方法。模型建立时，对于连接管渠写出连续性和动量 平衡方程，对于节点写出水量平衡方程。动力波法可以模拟管渠的蓄变、回水、逆流和有 压流动等复杂流态。

1)控制方程

控制方程分为连续方程和动量方程，

连续方程：

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

动量方程：

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

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

假设v表示平均流速，将代入方程(4.3-7)中的对流加速度项可得方程：

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

将方程上述方程代入动量方程中，得方程：

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

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

可求得Q_t+Δt：

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

此外，为考虑管道的进出口水头损失，可以从H₂和H₁中减去水头损失。因此，还需要有Q和H有关的方程，可以从节点方程得到。

2)节点控制方程

管网和渠道的节点控制方程为：

式中：H为节点水头，m；Q_t为进出节点的流量，m³/s；A_sk为节点的自由表面积， m²。

化为有限差分格式为：

上述两方程联立，可依次求得Δt时段内每个连接段的流量和每个节点的水头。

一二维耦合模型的建立：如图2所示，包括如下步骤：

步骤S1，输入地表水数据和地下水数据，根据所述地表水数据和地下水数据建立地表 地下连接关系。

步骤S2，生成并优化网格单元。

下面首先对网格剖分技术进行说明。

对于连续的物理系统的数学描述，通常是用偏微分方程来完成的。为了在计算机上实 现对这些物理系统的行为或状态的模拟，连续的方程必须离散化，在方程的求解域上(时 间和空间)仅仅需要有限个点，通过计算这些点上的未知变量既而得到整个区域上的物理 量的分布。有限差分，有限体积和有限元等数值方法都是通过这种方法来实现的。这些数 值方法的非常重要的一个部分就是实现对求解区域的网格剖分。

到目前为止，结构化网格技术发展得相对比较成熟，而非结构化网格技术由于起步较 晚，实现比较困难等方面的原因，现在正在处于逐渐走向成熟的阶段。有限元网格生成方 法按其所生成的单元类型可分为结构化(Construction)网格生成方法和非结构化(Non-Construction)网格生成方法。

从严格意义上讲，结构化网格是指网格区域内所有的内部点都具有相同的毗邻单元。 结构化网格有很多优点：(1)它可以很容易地实现区域的边界拟合，适于流体和表面应力 集中等方面的计算；(2)网格生成的速度快；(3)网格生成的质量好；(4)数据结构简单；(5)对曲面或空间的拟合大多数采用参数化或样条插值的方法得到，区域光滑，与实际的模型更容易接近。它的最典型的缺点是适用的范围比较窄。尤其随着近几年的计算机和数值方法的快速发展，人们对求解区域的复杂性的要求越来越高，在这种情况下，结构化网格生成技术就显得力不从心了。

同结构化网格的定义相对应，非结构化网格是指网格区域内的内部点不具有相同的毗 邻单元，即与网格剖分区域内的不同内点相连的网格数目不同。从定义上可以看出，结构 化网格和非结构化网格有相互重叠的部分，即非结构化网格中可能会包含结构化网格的部 分。

由于非结构化网格的生成技术比较复杂，随着人们对求解区域复杂性的不断提高，对 非结构化网格生成技术的要求越来越高。非结构化网格生成技术中只有平面三角形的自动 生成技术比较成熟，平面四边形网格的生成技术正在走向成熟。而空间任意曲面的三角形、 四边形网格的生成，三维任意几何形状实体的四面体网格和六面体网格的生成技术还远远 没有达到成熟，需要解决的问题还非常多。主要的困难是从二维到三维以后，待剖分网格 的空间区域非常复杂，除四面体单元以外，很难生成同一种类型的网格，需要各种网格形 式之间的过渡，如金字塔形，五面体形等等。

非结构化网格技术的分类，可以根据应用的领域分为应用于差分法的网格生成技术 (grid generation technology)和应用于有限元方法中的网格生成技术(meshgeneration technology)，应用于差分计算领域的网格除了要满足区域的几何形状要求以外，还要满足 某些特殊的性质(如垂直正交，与流线平行正交等)，因而从技术实现上来说就更困难一些。 基于有限元方法的网格生成技术相对非常自由，对生成的网格只要满足一些形状上的要求 就可以了。

在本步骤中，采用区域分解法或铺路法实现对网格的划分。下面分别对这两种方法的 具体网格剖分过程进行说明。

(1)参考图3，采用区域分解法生成并优化网格，包括如下步骤：

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

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

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

1)内外边界的合并

在区域分解法中，要划分网格的区域必须是单连通域，这样区域和子区域都可以用一 个连续的边界节点环来表示。对于内部有孔洞的区域，在网格划分之前，必须将多连通域 转化成单连通域，即需要对内外边界进行合并。图4所示的区域中，内部含有一个孔洞， 区域的边界不能用一个边界节点环来表示。图4(a)表示最初生成的内外边界节点，图中的 箭头表示内外边界节点的方向，即外边界节点是以逆时针顺序给出，内边界节点是以顺时 针顺序给出。从内外边界节点中分别选出一个节点来，通过连接这两个节点形成一条切割 线，如图4(b)所示，假想用这条切割线将区域切开一条宽度为0的缝隙，如图4(c)所示， 这样就可以将内部含有孔洞的多连通区域转化成只有一个外部边界的单连通区域。最后在 切割线上，按照网格密度要求生成新的节点，如图4(d)所示。

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

选择切割线的位置遵循以下两条原则：

(1)连接切割线的内部边界节点处的内角应大于或等于180度；

(2)切割线是所有连接内外边界节点的线段中长度最短的一条。

当要划分的区域包含多个孔洞，即有多个内部边界时，内外边界合并的方法是：对每 个内部边界选择与外部边界的切割线，在所有的内部边界中，选择切割线最短的内部边界 与外部边界进行合并，形成新的外部边界，以此类推，对剩余的内部边界进行合并，直到所有的内部边界都与外部边界进行合并。

图5为内部含有四个孔洞的区域内外边界合并的示意图。在这四个孔洞与初始外边界 的切割线中，第1个孔洞的切割线最短，因此首先第1个孔洞与初始外边界进行合并，形 成新的外边界，如图5(a)所示。在剩余的3个孔洞中，选择第4个孔洞与新的外边界进行合并，如图5(b)所示。依次类推对剩余的2个孔洞进行合并，如图5(c)、图5(d)所示， 最后得到只有一个外部边界的节点环。

2)边界单元的生成

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

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

具体来说，首先将边界节点向内部偏移一定的距离(等于该点相邻两边长度的平均值)， 形成新的边界节点，如图6(a)和图7(a)所示。依次连接外部边界节点和偏置节点生成一层 边界单元，如图6(b)和图7(b)所示，从图中可以看出边界单元形状接近正方形。

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

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

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

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

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

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

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

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

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

下面对网格剖分的剖分准则进行说明。

从区域边界节点环中选择两个节点，通过连接这两个节点形成一条剖分线将区域一分 为二。如何选择剖分线不仅影响到最后生成的网格质量，也影响到网格生成的效率。目前 大部分文献都是以剖分角度、剖分线的长度以及剖分后两子区域的面积比作为选择剖分线 的准则。如图7所示，连接节点i，j生成一条剖分线ij，将区域分成两个子区域，同时在节点i，j处形成四个剖分角γ_ij1，γ_ij2和γ_ji1，γ_ji2。原则上四个剖分角越接近直角，剖分线的长度越短，两子区域面积比越接近1，那么这样的剖分线就是最佳的剖分线。通常的 做法是将上述三个标准进行量化，采用加权平均的方法得到一个量化值，对所有可能的剖 分线比较其量化值，选择量化值最小的剖分线作为最佳的剖分线。

但实践中发现，上述方法存在如下的缺点：(1)三个标准量化权值的最优值很难确定， 文献中给出的权值各不相同，所给出的经验值，并不一定适用所有的情况；(2)计算时间长， 对于有n个节点的节点环，可能的剖分线条数为n(n-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₂为剖分两子区域的面积。

目前的应用，很难给出适用于各种情况的最优权值。不过在上述四个指标中，前两个 指标与最后生成的网格质量密切相关，因此权重最大。只要保证前两个指标的权重最大， 一般情况下都能生成质量较高的网格。增大w₄的数值可减少网格划分的时间，但网格质量 会下降，减少w₄的数值可提高网格的质量，但划分时间明显增加。参考的w₁,w₂,w₃,w₄值分 别为0.4，0.4，0.1和0.1。采用上述寻找最佳剖分线的方法，即使四个权值取任意值， 也能生成较为满意的网格。

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

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

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

区域(或子区域)边界节点处的网格密度值用该节点相邻两线段长度平均值的倒数来表 示。假定剖分线上的网格密度为线性分布。设剖分线两端节点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个节点的网格密度值和节点坐标，以此类推，求出剖分线上所有节点的网格密度值和节点坐标。

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

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

每当剖分后的子区域包含六个节点，剖分过程结束，程序转入六节点闭合处理子程序， 采用固定的模板进行模式匹配，生成最后的单元。

根据六节点区域几何形状可大致分为矩形(四条直边)、三角形(三条直边)、半圆形(两 条直边)和圆形(五条或六条直边)四类，如图9示。每一类又根据几何形状的每条边包 含网格单元边的数目分成三种情况(圆形类实际为两种情况)。例如矩形图形2-1-2-1表示 四条边包含网格单元边的数目分别为2、1、2和1，通过连接相对边上的两个节点，将区域分解成两个四边形单元；矩形图形3-1-1-1表示四条边包含网格单元边的数目分别为3、1、1和1，在矩形的内部增加两个新节点，将区域分解成四个四边形单元；矩形图形2-2-1-1表示四条边包含网格单元边的数目分别为2、2、1和1，在矩形的内部增加一个新节点， 将区域分解成三个四边形单元。图中每种模板编号的数字之和等于6(即六节点区域边的 条数)，六节点区域可能的模板总数为11。

从图9中可以看出，六点闭合处理生成的四边形单元质量较差，但是这些单元一般在 区域的内部，可以通过后面的光滑处理来提高单元的质量。

进一步，在区域分解法中，在生成最后的网格单元之后，采用Laplace光滑算法对最 后生成的网格进行光滑处理。

由于直线剖分和闭合处理等原因，区域分解法生成的最后单元质量可能较差。为了改 进单元形状，可以用Laplace光滑算法对最后生成的网格进行光滑处理。

如图10(a)和(b)所示，Laplace光滑算法的原理。对于所有的内部网格节点，其 坐标用其相邻节点坐标的平均值来代替。例如图10(a)中的节点P，其相邻的节点有四个 分别是P1、P2、P3和P4，用这四个相邻节点坐标的平均值替换掉P点的坐标值，得到P 点的新位置如图10(b)所示。Laplace光滑算法通常采用迭代的方式对所有的内部网格节 点进行调整，迭代次数通常为4次左右迭代过程就可以收敛。

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

图11示出了生成网格的过程该方法将网格的生成过程分成几个基本步骤：

(1)选择起始点

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

(2)生成网格

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

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

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

(1)终止节点：

(2)边点，

(3)角点：

(4)转点：

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

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

根据对网格节点的不同定义，以不同类型节点为基点来生成新节点的算法是不同的：

以终止点和边点为基点的算法：

如图13所示，此时将由Ni-1，Ni，Ni+1三点生成一个新节点Nj，其中：

V与Ni-1Ni的夹角为α/2。令d1为节点Ni-1和Ni的距离，d2为节点Ni+1和Ni的 距离，以下节点生成算法中d1、d2的意义与此相同。

以角点为基点的算法：

如图14所示，此时将由Ni-1，Ni，Ni+1三点生成三个新节点Nj，Nk，Nl，其中：

|V_l|＝|V_j|

Vj，Vk，Vl与Ni-1Ni的夹角分别为α/3，α/2，2α/3。

以转点为基点的算法：

如图15所示，此时将由由Ni-1，Ni，Ni+1三点生成五个新节点Nj，Nk，Nl，Nm，Nn 其中：

|V_l|＝|V_j|

|V_m|＝|V_k|

|V_n|＝|V_j|

Vj，Vk，Vl，Vm，Vn与Ni-1Ni的夹角分别为α/4，3α/8，α/2，5α/8，3α/4。

(3)缝合处理边界

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

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

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

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

(b)过渡缝合处理

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

通常，对单元尺寸过渡较大的边界来说，在过渡处容易导致单元的长边与短边相差太 大，以此为边界继续往里生成新的网格时，会造成网格的畸变和交叉。所以，当长边与短 边相差到一定程度时，规定在长短比大于2时，需要进行缝合处理。如图所示，过渡缝合的方法是：在长边的中点增加一点d，连接cd及da并消除db线段。这样相当于增加一个 新的四边形单元abcd。

(4)交叉处理

在新单元生成的过程中，由于边界的不规则，新单元往往会与旧边界发生交叉和重叠 现象，中断网格的生成，进行交叉处理。

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

1)交叉的判断

如图18所示，两线段AB，CD其端点分别为A(x1，y1)、B(x2，y2)和C(x3，y3)、D(x4，y4)，假设它们的交点为P(x，y)，由两点式直线方程可得两线段方程分别为：

整理得

(y₂-y₁)x-(x₂-x₁)y＝x₁y₂-x₂y₁

(y₃-y₄)x-(x₄-x₃)y＝x₃y₄-x₄y₃

行列式

根据克莱姆法则，当D≠0时，方程有解，即AB和CD相交。为了分析交点可能出现的情况，我们用参数μ和ω分别建立线段AB和CD的方程如下：

AB段：A+μ(B-A)

CD段：C+ω(D-C)

如果有交叉发生，则可得到以下两个方程：

x1+μ(x2-x1)＝x3+ω(x4-x3)

y1+μ(y2-y1)＝y3+ω(y4-y3)

解方程，当μ和ω的值均在区间[0，1]之间时，两线段相交。

2)交叉的处理

当由边界向域内逐层生成网格时，由于边界的不规则，层与层之间多会发生大范围的 交叉和重叠现象，这种现象不仅难于处理，而且层与层的交汇处还易产生不规则单元。重 新设计了网格生成算法，不再仅基于一个边界按层渐进，而是采用多个动态边界，在每一 个边界上，一段一段，甚至一个一个单元的生成新网格系统。本算法在新单元生成的过程 中，随时判断有无交叉的发生，当发生交叉时，网格生成立刻中断，进行交叉处理，而不 是等层生成完毕时再处理，这样就避免了大面积的网格交叉和重叠现象，只是局部个别网 格的交叉和重叠，简化了交叉重叠问题的处理，也使算法更稳定可靠。

如图19所示，当发生交叉时，采用将重叠的部分连接的方法来修正网格。连接将导致 形成两个新的内部边界。在许多情况下，连接交叉的两边是很好的方法。为确保连接有效 同时避免生成不规则节点，进行了以下处理：

(i)偶数限制。

必须保证连接后的边界节点总数均为偶数。当连接完毕后，必须验证两个新边界的节 点数，若为偶数，则边界不变；若为奇数，则采用调整对应连接边的办法解决。

(ii)连续缝合

交叉边界连接完毕后，有时连接交叉的节点处会出现小角度，若不处理，就会产生不 规则单元。因此，本算法在交叉处理完毕后，还要进行缝合处理，保证生成的新边界的质量，以利于后续网格的生成。

(5)边界调整

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

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

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

如图20所示，楔块插入过程中，生成两个新节点和一个新单元。在节点Ni处插入楔块，首先将节点Ni位置调整到线段NiNi-1的1/3处，新生成节点Ni’位于原线段NiNi+1 的1/3处，同时节点Ni’将替代节点Ni在右边单元中的位置。这样，就形成了一个裂缝。 第二个新节点Nk位置，可将节点Ni，Nj，Ni’将带入光滑公式求得。从而得到如图19(b) 所示的新单元。

(6)边界的光滑处理

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

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

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

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

Δ_A＝V′_i-V_i

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

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

如图21所示，这种调整可简单理解为只是对向量Vij的长度进行调整，而不调整其方 向。长度调整可保证单元的尺寸，但不能保证单元内角的垂直度。

然后，进行角度的光滑处理。如图22所示，令节点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)闭合处理

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

具体地，当边界包括四个节点时，则将其作为一个单元处理；当边界包括六个节点时， 则可根据其中包含的终止点个数n来进行闭合处理。闭合处理后，有些四边形单元较不规 则，这是由于闭合处理多发生在形心所致，可通过网格质量优化得到解决。

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

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

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

1)节点消去

如图23，如果节点A的周围有两个单元，即NEA＝2，则节点A将被消除，相应的两个单元合并成了一个单元。节点消去将对整个内部节点进行循环，直到消除所有的此类节点。

2)单元消去

如图24图(a)中，由节点Ni、Nj、Nk和Nl组成的单元Ei中，只有Ni为规则节点。 单元对角线上的两节点Nj和Nl周围各有三个单元，Nk周围有五个单元。通过将Nj和Nl 合并，删除这个单元。其结果如图23)所示，原来的规则节点Ni变为不规则节点，而原来 的三个不规则节点Nj、Nk和Nl则变为规则节点。因此，网格的整体质量提高了。单元消 去也将对整个内部单元进行循环，直到消除所有的此类单元。

3)对角线替换

对角线替换算法将对所有由内部节点连接的单元边进行检验。如图25所示，由节点A， F，D，B，E，C组成的六边形构成了两个相邻单元Ei和Ej，如果AB边的单元连接数目之和满足下式：

NE_A+NE_B≥9

那么，当下式成立时，六边形对角线AB将被CD替换，如图25b示，

(NE_A+NE_B)-(NE_C+NE_D)≥(NE_A+NE_B)-(NE_E+NE_F)≥3

如果，下式成立，六边形对角线AB将被EF替换，如图25c示，

(NE_A+NE_B)-(NE_E+NE_F)≥(NE_A+NE_B)-(NE_C+NE_D)≥3

如果上述两式都不满足，则不进行替换。

4)边消去

边消去算法将对每一个由两内部节点连接的边进行判别。如图26a所示，EF边的两端 点的单元连接个数都为3，即NEE＝3和NEF＝3，且包含EF边的两个单元Ei和Ej中都不含有初始边界点，则消去EF边，同时也消去了单元Ei和Ej。

当A，B，C，D各点的单元连接个数满足以下条件，则连接AD边形成两个新单元，图25b所示，

NE_A+NE_D≤NE_B+NE_C

否则，连接BC边形成两个新单元，如图26c所示。

图27为根据本发明实施例的网格自动生成的流程图。

步骤S3，分别计算地表径流和排水管网水流。

步骤S4，计算地表径流和排水管网水流的交换水量，包括：

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

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

步骤S5，对计算得到的交换水量进行校核，以实现地表二维地下管网水动力学耦合。

具体地，(1)由于一个网格单元对应多个管网节点，以二维网格单元为单位校核拟交 换水量是否超过单元现有总水量，当出现二维网格单元中总水量不够，无法满足当前与众 多管网节点计算的交换流量时，需要按比例减少交换水量；

(2)由于交换水量是根据当前步结果显式计算的，未考虑下一时段网格单元以及管道 的来水量，可能会出现交换流量过大的情况。若上一时间步二维网格与管网节点之间水流 交换方向为网格单元流入管网节点，而上一步计算完成后，该节点出现溢流，说明上一步 交换水量过多，需要当前步中将网格单元的水量增加该溢流值，以满足水量平衡。

对于区域内部有特征约束线的情况，如图28(a)所示，正方形的区域内部有一条约束线， 首先在边界和约束线上生成网格节点，约束线上的节点为最终的网格节点，其位置固定在 约束线上。将约束线看作是面积为零的孔洞，如图28(b)所示，约束线上节点形成封闭的 节点环，除了端节点外，其余的节点在节点环中出现两次。将约束线看作是内部孔洞后， 就可以采用上面的方法对内部边界与外部边界进行合并，形成一个单连通区域，如图28(c) 所示。

但是，由于约束线上的节点形成的节点环没有方向性，即无法判断其节点顺序是逆时 针还是顺时针，当切割线的一个端点位于约束线的中间时，必须保证最后合并的边界不能 出现交叉的情况。例如，对于图29(a)所示的约束线，如果按照图29(b)所示的节点顺序 合并内外边界，则最后合并的边界会出现交叉的情况，而按照29(c)所示的节点顺序合并 内外边界，则为正确的边界合并。另外，当约束线的一个端点位于区域的边界上时，仍可将其作为孔洞处理，只不过该约束线与边界合并时，切割线的长度为零，重合处的节点采用相同的编号。

对所有的内部特征处理完毕后，要划分的区域就转换成单连通区域，可以采用区域分 解法对其进四边形网格划分。

当区域内部多条约束相交形成一个或多个封闭的区域时，或者约束线连接区域的边界 将区域分隔成多个子区域时，则不能直接采用上面介绍的约束线处理方法，必须先将要划 分网格的区域分解成独立的子区域，然后在每一个区域上生成四边形网格。例如，图30所 示的正方形区域内有5条约束线，这些约束线相交在区域的内部形成2个封闭的子区域， 去掉这2个子区域后剩下的带有孔洞的区域为第3区域。这样约束线和边界将正方形区域 分隔成3个独立的子区域，对于区域1和2可以直接生成四边形网格，对于区域3，孔洞 上的约束线作为区域3的内部边界，其余的约束线作为区域3里的约束线，这样就可以按 照上面介绍的约束线处理方法，对区域3进行网格划分。

对于上面的约束线分隔区域的情况，可以按照下面的步骤自动确定各个子区域的边界：

(1)计算约束线与边界以及约束线与约束线间的交点，用交点将边界和约束线分隔成 小的线段。例如，图中的约束线被交点分隔成17条线段；

(2)对于区域的外边界和内部孔洞用有向线段来表示，外部边界上的线段以逆时针顺 序给出，内部边界上线段以顺时针顺序给出，将这些有向线段加入到集合S中；

(3)将约束线段加入到集合L中，对约束线段的端点进行计数，对于端点计数为1的约束线段从集合L中删除，这些约束线段不形成子区域的边界，而是作为子区域内的约束线；

(4)集合L中的约束线段在形成子区域的边界时将被使用2次，所以先将集合L中的线段加入到集合S中，然后将集合L中的各线段改变方向后再次加入到集合S中；

(5)从集合S中任取一条线段加入到集合M中，并从集合S中删除这条线段，然后从集合S中寻找与这条线段首尾相连的下一条线段，如果有多条线段与这条线段首尾相连，则选择与这条线段夹角最小的那条线段加入到集合M中，并从集合S中删除这条线段，重 复这个过程直到集合M中的线段形成一个有向封闭的环。

(6)若集合S不为空，则重复步骤(5)寻找所有的有向封闭环。

(7)计算所有有向封闭环的面积，对于面积为正的封闭环作为每一个子区域的外边界， 对于面积为负的封闭环则为子区域的内部边界，并确定该封闭环为哪个子区域的内部边界， 这样每一个子区域的边界就确定了。

网格生成程序应能根据用户的要求，对某些区域进行适当的网格加密，有效地控制网 格密度以及最后生成的单元数目。从直观上看，网格密度大的地方，网格单元尺寸较小； 网格密度小的地方，网格单元尺寸较大。从数学上可以定义网格密度为网格单元长度的倒 数。在要划分的区域里，网格密度分布应该是连续的并且过渡光滑，这样根据网格密度生 成的网格单元才会过渡均匀。二维网格密度是坐标x，y的函数，选择Laplace方程来表示 要划分区域内的网格密度分布函数，即

(在Ω内)

式中ρ----网格密度分布函数，Ω----计算区域

边界条件为：

(在Γ上)

式中Γ----计算区域的边界，n----边界Γ的法线方向

选择Laplace方程作为密度函数有两个原因：一是用该方程表示的密度分布过渡光滑； 二是可以按照下面介绍的方法，应用有限元技术求得网格密度的数值解。

用有限元方法求解Laplace方程的数值解，首先需要将计算区域离散化，然后对单元 进行分析，得到单元刚度方程，然后组装成整体刚度方程，最后施加边界条件，求解整体刚度方程即可得到方程的数值解。

为了求解区域上的网格密度，第一步要将区域离散成计算单元。在包含区域的最小矩 形里，均匀划分栅格单元，如图31(a)所示，然后去掉不在区域的栅格单元，即可将区域 离散成矩形的计算单元，如图31(b)所示，栅格单元节点上的密度值为待求的未知量。

第二步对栅格单元进行分析。任取一栅格单元如图32所示。栅格单元内任意一点的密 度值都可以用单元四个节点的密度值经插值而得到。栅格单元四个节点的形函数分别是：

q₁＝(1-ξ)(1-η)/4

q₂＝(1+ξ)(1-η)/4

q₃＝(1+ξ)(1+η)/4

q₄＝(1-ξ)(1+η)/4

则栅格单元内任意一点的密度值为：

式中ξ,η----栅格单元的局部坐标

ρ_i----栅格单元的第i个节点处的密度值

经整理后得到单元刚度方程如下：

Kρ＝0

式中K----单元刚度矩阵

ρ----单元节点处的密度值，其中

K＝∫MM^TdΩ

ρ＝{ρ₁,ρ₂,ρ₃,ρ₄}

将计算区域里所有单元刚度方程进行组装可得到整体刚度方程，对整体刚度方程施加 密度边界条件，然后求解整体刚度方程，就可得到各栅格节点上的密度值。当所有栅格节 点上的密度值确定后，区域内任意一点的网格密度值就可以通过该点所在的栅格单元的四 个节点上的密度值经插值得到。

在一些应用场合中，需要对物体某些区域的网格适当加密，这时采用密度窗口的方式 设定网格密度就显得非常方便灵活。由用户在这些区域指定多边形的密度窗口，给每个窗 口设定一个相对密度值，如果密度窗口与边界相交，计算交点的位置，交点处的密度值等 于窗口的密度值，施加密度边界条件，求解刚度方程，得到各栅格节点上的密度值，最后 将密度窗口内栅格节点上的密度值用该窗口的密度值替换掉。这样根据密度窗口就可以得 到区域各栅格节点上的密度值了。需要说明的是，由于窗口的密度值采用的是相对值，因 此至少需要定义两个密度窗口才能体现出相对值。当密度窗口与边界没有交点时，边界上 的相对密度值假定为1。

网格密度可以采用绝对值，也可以采用相对值。对于相对密度值需要根据划分的单元 数目进行调整，调整系数k可以按下面的公式确定：

式中ρ*---相对密度函数

N---为要划分单元的数目

A----为由栅格单元组成的区域。

上式中的积分项可利用栅格节点上的相对密度值，在各栅格单元上计算其数值解，这 样就可以得到调整系数k值，用k值乘以栅格节点上的相对密度值，就得到了各栅格节点 上的网格密度值。

根据本发明实施例的基于城市地表地下管网的一二维水动力学耦合分析方法，对产流 计算模型、一维管网模型和一二维耦合模型进行单独建立，对每个模型进行单独的数据处 理，从而完成对城市地表地下管网的一二维水动力学耦合的分析。

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

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

Claims (9)

1.一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，建立城市地表地下管网模型，城市地表地下管网模型包括产流计算模型、一维管网模型和一二维耦合模型，

建立产流计算模型的方法为：将计算区域划分成若干个子汇水区，每一个子汇水区概化为一个非线性蓄水池，其入流项有降水和来自上游子流域的流出量，只有当蓄水池水深超过最大洼地蓄水量时地表出流才会发生，其大小通过曼宁公式计算得出：

式中：Q为出流量，m³/s；W为子流域宽度，m；S为坡度；n为曼宁粗率系数；

建立一维管网模型的方法为：计算管渠的汇流，通过控制方程计算连接管渠的连续性和动量平衡性，通过节点控制方程计算节点的出水量，

建立一二维耦合模型，包括以下步骤：

步骤S1，输入地表水数据和地下水数据，根据所述地表水数据和地下水数据建立地表地下连接关系；

步骤S2，生成并优化网格单元；

步骤S3，分别计算地表径流和排水管网水流；

步骤S4，计算地表径流和排水管网水流的交换水量；

步骤S5，对计算得到的交换水量进行校核，以实现地表二维地下管网水动力学耦合。

2.如权利要求1所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在所述产流计算模型中的每个子流域产流包括三种不同类型的地表产流：

无洼蓄不透水地表上的降雨损失主要为蒸发，产流量表示为：

R₁＝P-E

式中：R₁为无洼不透水地表的产流量；P为降雨量；E为蒸发量，

有洼蓄不透水地表上的降雨损失主要为填洼和蒸发，产流量表示为：

R₂＝P-D-E

式中：R₂为有洼不透水地表的产流量；D为洼蓄量，

透水地表的降雨损失主要包括洼蓄和下渗，产流量表示为：

R₃＝(i-f-E)×Δt

式中：R₃为透水地表的产流量；i为降雨强度；f为入渗强度；Δt为时间间隔。

3.如权利要求1所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在所述建立一维管网模型中，控制方程又包括连续方程和动量方程，连续方程：

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

动量方程：

式中：H为水深；g为重力加速度；S_f为摩阻坡度，

节点控制方程：

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

4.如权利要求1所述的基于空间拓扑的地表二维地下管网水动力学耦合方法，其特征在于，在所述步骤S5中，对所述交换水量进行校核，包括：

(1)由于一个网格单元对应多个管网节点，以二维网格单元为单位校核拟交换水量是否超过单元现有总水量，当出现二维网格单元中总水量不够，无法满足当前与众多管网节点计算的交换流量时，需要按比例减少交换水量；

(2)当上一时间步二维网格与管网节点之间水流交换方向为网格单元流入管网节点，而上一步计算完成后，该节点出现溢流，说明上一步交换水量过多，需要当前步中将网格单元的水量增加该溢流值，以满足水量平衡。

5.如权利要求1所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在所述步骤S2中，采用区域分解法生成并优化网格，包括如下步骤：

在待划分区域的边界上生成边界节点，通过连接两个边界节点，将该区域剖分成两个子区域；

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

6.如权利要求1所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在所述步骤S2中，采用铺路法生成并优化网格，包括如下步骤：

(1)选择起始点

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

(2)生成网格

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

(3)缝合处理边界

对新生成的边界中相邻边界长短悬殊和节点夹角过小现象，进行边界的缝合处理；

(4)交叉处理

在新单元生成的过程中，由于边界的不规则，新单元往往会与旧边界发生交叉和重叠现象，中断网格的生成，进行交叉处理；

(5)边界调整

当边界呈凸起状时，通过楔块插入法，改善网格尺寸大小，对边界进行调整；

(6)边界的光滑处理

对所有内部边界节点进行光滑处理；

(7)闭合处理

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

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

当网格生成完毕后，对内部网格进行光滑和整合处理，以消除存在的缺陷。

7.如权利要求6所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在铺路法中，所述缝合处理边界，包括：

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

当内部边界节点Ni的夹角α小于π/6时，进行缝合处理，包括：将Ni前一节点Ni-1与其后一节点Ni+1合并成一个节点Nj，合并后Nj的坐标为：

(2)过渡缝合处理

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

8.如权利要求6所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在铺路法中，所述交叉处理，包括：

当判断出现网格交叉重叠情况时，采用将重叠的部分连接的方法来修正网格，连接将导致形成两个新的内部边界，为确保连接有效同时避免生成不规则节点，进行了以下处理：

(1)偶数限制

保证连接后的边界节点总数均为偶数

(2)连续缝合

交叉边界连接完毕后，进行缝合处理，保证生成的新边界的质量，以利于后续网格的生成。

9.如权利要求6所述的一种基于城市地表地下管网的一二维水动力学耦合分析方法，其特征在于，在铺路法中，所述闭合处理，包括如下步骤：

当边界包括四个节点时，则将其作为一个单元处理；当边界包括六个节点时，则根据其中包含的终止点个数进行闭合处理。