CN111046567A - 一种基于Godunov格式的城市排水管网水流数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于Godunov格式的城市排水管网水流数值模拟方法。首先建立各管段和各节点的对应关系，管段采用有限体积单元离散，水力要素存储在单元中心。根据各管段首末端单元水位或压力水头和与其相连的内节点的水位关系来选择内节点的边界条件类型，然后与外边界条件和本时刻各单元水力要素值通过HLL近似Riemann解计算各单元界面处的数值通量，从而更新所有管段单元的水力要素值到下一时刻。同时，根据内节点边界条件计算获得进入或离开内节点的流量，通过内节点处的水量平衡方程计算下一时刻内节点的水位值。该处理方法可以克服传统排水管网模型使用有限差分法带来的局限性，为城市排水管网水流数值模拟提供了一种新的解决方案。

Description

一种基于Godunov格式的城市排水管网水流数值模拟方法

技术领域

本发明涉及市政工程领域，尤其涉及城市防洪领域，具体为一种基于Godunov格式的城市排水管网水流数值模拟方法。

背景技术

城市排水管网是指由一系列排水管段首尾相互连接形成的具有水力联系的统一整体，管段和管段之间通过雨篦子、检查井和雨水井等进行连接，这些连接点定义为节点。在城市排水管网计算中，外边界条件施加的节点一般定义为外节点，其它节点为内节点。

在我国快速城市化的大背景下，随着极端天气频发，城市洪涝问题日益突出。城市洪涝模型作为一种有用的工具，能够为城市洪涝问题的防治工作提供技术支持。在城市洪涝模型中，城市排水管网水流的模拟是必不可少的重要一环。当前的排水管网模型以国际商业化的软件产品为主，其中以美国环保署(EPA)开发的SWMM(StormWaterManagementModel)模型最为流行。国内也有学者从事该方面的研究工作，但还未有成熟的排水管网产品。

当前，在城市排水管网的计算领域使用的主流方法都是始于上世纪70年代之后流行的有限差分法，比较经典的有SWMM模型、MIKEURBAN的MOUSE模型等，应该说这些模型经过多年的发展完善，应用上已相对比较成熟。但是由于数值格式自身的缺陷性，上述这些模型在数值方法上依然存在先天缺陷，如SWMM模型最早为市政规划部门服务，对排水管网水流的计算要求没有那么细致，所以没有对管段进行分段，建立了“LINK-NODE”模式，沿袭到现在该模型依然无法对管段进行分段，对一些比较长的管段，这无疑会带来较大的数值误差；另一方面，当管段埋设的坡度较大或者当管段内出现急、缓流的变化时，由于有限差分法本身的缺陷，为保持计算稳定，当上述情形出现时要采用概化的方法，这不但会损失掉一部分计算精度，同时也会带来数值稳定性问题。

发明内容

本发明的目的在于提供一种基于Godunov有限体积格式的城市排水管网水流数值模拟方法，弥补现有传统方法的缺陷。

以求解Riemann近似解为基础的Godunov格式被称为新一代的水动力模型格式，该格式既能模拟光滑的古典解，又能很好的模拟大梯度的水面流动。该格式在河网水流模拟中有应用，但在城市排水管网水流模拟中尚无应用。城市排水管网与河网有两点突出的不同：1)河网水流全为明渠水流，但排水管网中明渠流和压力流是交替出现的；2)河网的节点范围并不清晰，在河网求解中，一般假定与节点相连接的河段处的水位是相等的，但是排水管网的节点非常明确，为检查井和雨水井等，但节点内的水流形态更为复杂，管段的埋设底高程与节点底高程并不一定一致，增加了节点处理难度。鉴于此，为了突破传统城市排水管网模型应用有限差分法时遇到的困难，将Godunov格式推广应用到城市排水管网水流数值模拟领域具有重要的学术意义和工程应用意义。

本发明是通过以下技术方案实现的：

一种基于Godunov格式的城市排水管网水流数值模拟方法，采用窄缝法模拟明渠流和压力流的变化过程，采用HLL格式计算管段单元界面通量，根据各管段首末端单元水位或压力水头和与其相连接的管网内节点水位的关系确定内节点使用的具体边界条件类型，从而计算流入或流出管网节点的流量，求解管网节点的水量平衡方程后为各管段提供下一时刻的内边界水位。包括如下具体步骤：

1)获取计算范围内城市排水管网数据，包括管网平面几何拓扑数据、各管段的管径尺寸、各管段首末端的管底高程、各管网节点底高程、各管网节点的最大水深和超载水深、各管网节点面积等数据。

2)对各管段和节点进行编号，建立各管段和节点的对应关系，约定水流从节点底高程高的一端流到节点底高程低的一端为正方向，采用有限体积单元离散各排水管段，每条管段的离散单元沿正方向依次编号；水力变量存储在各单元中心点；确定排水管网的入流边界和出流边界，入流边界和出流边界统称为外边界，施加外边界条件的节点称为外节点，其余节点称为内节点。

3)城市排水管网模型初始化：给各管段一维有限体积单元设置糙率和水力条件初始值；给管网内节点设置初始水力条件。

4)确定计算时间步长dt：城市排水管网模型采用显格式的有限体积法求解，该时间步长获取受到CFL(Courant-Friedrichs-Lewy)条件限制；获取当前时刻各管段外边界条件，当前时刻以t时刻表示。

5)采用管道明渠流和压力流统一的方程组描述管道水流运动，具体形式如下：

连续方程：

动量方程：

式中Q为断面流量，A为过水断面面积，g为重力加速度，t为时间变量，s_f为沿程阻力损失，其表达式为s_f＝(n²Q|Q|)/(A²R^4/3)，R为水力半径，n为Manning糙率系数，当管道内为明流时，B为水面宽度，H为水位值；当管道内为压力流时，B＝gA/a²，a为压力波波速，实际计算时，B赋值一很小的窄缝宽度，H为压力水头。

6)根据内节点水位和与其相连的各管段的首(末)端单元的水位值(压力水头值)的关系来判断使用不同的内节点边界条件，从而计算流入或流出内节点的水量，具体过程如下：

a)计算与内节点连接的上游管段流入内节点的流量时，分如下几种情形:

①如果上游管段最后一个单元有水，并且内节点水位值小于该单元底高程值，则此时按自由出流边界条件计算进入内节点的流量；

②如果上游管段最后一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算进入内节点的流量；

③如果上游管段最后一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算进入内节点的流量，此种情况下，水流倒流，进入内节点的流量为负值；

④如果上游管段最后一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算进入内节点的流量，此时进入内节点的流量为0；

b)计算与内节点连接的下游管段流出内节点的流量时，分如下几种情形:

①如果下游管段第一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算流出内节点的流量；

②如果下游管段第一个单元有水，并且内节点水位小于该单元底高程值，则此时按自由出流边界条件计算离开内节点的流量，此种情况下，水流倒流，流出内节点的流量为负值；

③如果下游管段第一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算离开内节点的流量；

④如果下游管段第一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算离开内节点的流量，此时流量为0。

7)通过t时刻各管段单元的水力要素值、外边界条件以及内节点边界条件，采用HLL显格式计算通过单元界面处的数值通量，计算获得各管段t+dt时刻各单元的水力要素值；如果新计算出的某单元的水深值大于管径，则水流变为压力流状态，此时将管段水面宽度设定为窄缝宽度。

8)建立各内节点水量平衡方程，并根据该平衡方程获得t+dt时刻各内节点的水位值；

式中，

表示与内节点k相连接的第i个上游管段流入内节点的流量，

表示与内节点k相连接的第j个下游管段流出内节点的流量，A_k为第k个内节点的面积，

为第k个内节点在t时刻的水位，

为需要求解的第k个内节点在t+dt时刻的水位，m为与第k个内节点相连接的上游管段数，n为与第k个内节点相连接的下游管段数；q_k为第k个内节点与地面的水流交互量，当由内节点流到地面时，q_k为负值，反之为正值。

9)令t＝t+dt，重复步骤4)～8)，直到计算结束。

进一步的，步骤4)中dt的选取受到CFL(Courant-Friedrichs-Lewy)条件限制，具体如式(4)所示：

式中：u为流速；

c为波速，g为重力加速度，A为过流面积，B为水面宽度，当水流变为压力流后，B采用窄缝宽度；Δx为单元空间步长；dt为时间步长。

进一步的，步骤5)当水流为压力流状态时，水面宽度设定为窄缝宽度，理论上窄缝宽度越小，数值解越精确，但过小的窄缝值，计算的稳定性变差，综合考虑数值精度和稳定性的要求，建议窄缝取值范围为B＝0.01m～0.001m。

进一步的，步骤6)中通过内节点边界条件计算流入或流出内节点的流量是通过如下方法完成的：在上游管段的最后一个单元之后或下游管段的第一个单元之前设置一虚拟单元，该虚拟单元的特征跟与其相连的单元特征一致，根据不同的内节点边界条件类型设置虚拟单元的水力变量，步骤6)存在如下几种情形：自由出流边界条件时，虚拟单元的流量和水位等于与其相邻单元的流量和水位；水位边界条件时，虚拟单元的流量等于与其相邻单元的流量，水位等于内节点水位；固壁边界条件时，虚拟单元的流量等于与其相邻单元流量的负值，水位等于与其相邻单元的水位。确定好虚拟单元水力变量后，根据HLL近似Riemann解即可计算通过虚拟单元和与其相邻单元交界面处的质量通量值和动量通量值，该质量通量值即为流入或流出内节点的流量。

附图说明

图1为本发明的一种基于Godunov格式的城市排水管网水流数值模拟方法流程图；

图2为内节点边界条件类型示意图；

图3为城市排水管网算例示意图；

图中C1～C21为各管段的有限体积单元编号；

图4为城市排水管网算例的数值计算结果。

具体实施方式

下面结合附图1和实施例对本发明作进一步说明。

本发明提供的是一种基于Godunov格式的城市排水管网水流数值模拟方法。该方法采用窄缝法将明渠流和压力流的控制方程统一，采用具有良好激波捕捉能力的Godunov有限体积法格式计算界面通量，通过排水管段首末端单元水位或压力水头和与其相连接的内节点的水位关系来判断采用的内节点边界条件，通过求解内节点水量平衡方程来获得下一时刻的内节点水位。该方法能够有效弥补现有传统有限差分法的内在缺陷，提升排水管网的数值计算精度。该方法包括如下具体步骤：

1)获取计算范围内城市排水管网数据，包括管网平面几何拓扑数据、各管段的管径尺寸、各管段首末端的管底高程、各管网节点底高程、各管网节点的最大水深和超载水深、各管网节点面积等数据。

2)对各管段和节点进行编号，建立节点和各管段的对应关系，约定水流从节点底高程高的一端流到节点底高程低的一端为正方向，每条管段的离散单元沿正方向依次编号；采用有限体积单元离散各排水管段，水力变量存储在各单元中心点；确定排水管网的入流边界和出流边界，入流边界和出流边界统称为外边界，施加外边界条件的节点称为外节点，其余节点称为内节点。

3)城市排水管网模型初始化：给各管段一维有限体积单元设置糙率和水力条件初始值；给管网内节点设置初始水力条件。

4)确定计算时间步长dt：城市排水管网模型采用显格式的有限体积法求解，该时间步长获取受到CFL(Courant-Friedrichs-Lewy)条件限制；具体如式(4)所示：

式中：u为流速；

c为波速，g为重力加速度，A为过流面积，B为水面宽度，当水流变为压力流后，B采用窄缝宽度；Δx为单元空间步长；dt为时间步长。

获取当前时刻各管段外边界条件，当前时刻以t时刻表示。

5)采用管道明渠流和压力流统一的方程组描述管道水流运动，具体形式如下：

连续方程：

动量方程：

式中：Q为断面流量，A为过水断面面积，g为重力加速度，t为时间变量，s_f为沿程阻力损失，其表达式为s_f＝(n²Q|Q|)/(A²R^4/3)，R为水力半径，n为Manning糙率系数，当管道内为明流时，B为水面宽度，H为水位值；当管道内为压力流时，B＝gA/a²，a为压力波波速，实际计算时，B赋值一很小的窄缝宽度，理论上窄缝宽度越小，数值解越精确，但过小的窄缝值，计算的稳定性变差，综合考虑数值精度和稳定性的要求，建议窄缝取值范围为B＝0.01m～0.001m，H为压力水头。

6)根据内节点水位和与其相连的各管段的首(末)端单元的水位值(压力水头值)的关系来判断使用不同的内节点边界条件，从而计算流入或流出内节点的水量，具体过程如下：

a)计算与内节点连接的上游管段流入内节点的流量时，分如下几种情形:

①如果上游管段最后一个单元有水，并且内节点水位值小于该单元底高程值，则此时按自由出流边界条件计算进入内节点的流量；

②如果上游管段最后一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算进入内节点的流量；

③如果上游管段最后一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算进入内节点的流量，此种情况下，水流倒流，进入内节点的流量为负值；

④如果上游管段最后一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算进入内节点的流量，此时进入内节点的流量为0。

b)计算与内节点连接的下游管段流出内节点的流量时，分如下几种情形:

①如果下游管段第一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算流出内节点的流量；

②如果下游管段第一个单元有水，并且内节点水位小于该单元底高程值，则此时按自由出流边界条件计算离开内节点的流量，此种情况下，水流倒流，流出内节点的流量为负值；

③如果下游管段第一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算离开内节点的流量；

④如果下游管段第一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算离开内节点的流量，此时流量为0。

通过内节点边界条件计算流入或流出内节点的流量是通过如下方法完成的：在上游管段的最后一个单元之后或下游管段的第一个单元之前设置一虚拟单元，该虚拟单元的特征跟与其相连的单元特征一致，根据不同的内节点边界条件类型设置虚拟单元的水力变量，步骤6)存在如下几种情形：自由出流边界条件时，虚拟单元的流量和水位等于与其相邻单元的流量和水位；水位边界条件时，虚拟单元的流量等于与其相邻单元的流量，水位等于内节点水位；固壁边界条件时，虚拟单元的流量等于与其相邻单元流量的负值，水位等于与其相邻单元的水位。确定好虚拟单元水力变量后，根据HLL近似Riemann解即可计算通过虚拟单元和与其相邻单元交界面处的质量通量值和动量通量值，该质量通量值即为流入或流出内节点的流量。HLL近似Riemann解求解过程可参见文献描述(张大伟,程晓陶,黄金池等.复杂明渠水流运动的高适用性数学模型[J].水利学报,2010,41(4):531-536)。

为便于理解，图2给出了一个内节点边界条件的示例，由图2可以直观的看出，上游管段的最后一个单元有水，并且内节点水位小于该单元的底高程，所以此时上游管段末端为自由出流边界。由于下游管段第一个单元有水，并且内节点水位大于该单元管底高程，所以，此时下游管段首端为水位边界。

7)通过t时刻各管段单元的水力要素值、外边界条件以及内节点边界条件，采用HLL显格式计算通过单元界面处的数值通量，计算获得各管段t+dt时刻各单元的水力要素值；如果新计算出的某单元的水深值大于管径，则水流变为压力流状态，此时将管段水面宽度设定为窄缝宽度。

8)建立各内节点水量平衡方程，并根据该平衡方程获得t+dt时刻各内节点的水位值；

式中，

表示与内节点k相连接的第i个上游管段流入内节点的流量，

表示与内节点k相连接的第j个下游管段流出内节点的流量，A_k为第k个内节点的面积，

为第k个内节点在t时刻的水位，

为需要求解的第k个内节点在t+dt时刻的水位，m为与第k个内节点相连接的上游管段数，n为与第k个内节点相连接的下游管段数；q_k为第k个内节点与地面的水流交互量，当由内节点流到地面时，q_k为负值，反之为正值。

9)令t＝t+dt，重复步骤4)～8)，直到计算结束。

图3为城市排水管网算例示意图，该算例最早由Zhong Ji提出(ZhongZhi.General Hydrodynamic Model for Sewer/Channel Network Systems,Journal ofHydraulic Engineering(ASCE),1998,124(3):307-315.)，用来检验模型处理复杂排水管网的能力。该排水管网由6条管段和6个节点组成，其中管段1管径为0.8m，其余管段管径均为0.5m，管段1、管段2和管段5长度为300m，管段3长度为410m，管段4长度为500m，管段6长度为310m，各管段的划分单元数如图3所示。节点1、2、3为外边界条件所在节点，为外节点，其中节点1为流量外边界，节点2和3为水位外边界，节点4、5、6为内节点，其中内节点5有一入流过程汇入。各管段糙率统一取值为0.0125，窄缝宽取值为0.01m。

初始阶段，各管段中均无水流。该算例中，所有管段都要经历从无水到有水的状态，部分管段会出现明满流交替状态和水流流向往复状态，因此该算例能够有效检验模型处理干湿交替、明满流交替以及流向往复的能力。采用本发明方法与当前最流行的排水管网水流模拟模型SWMM进行比较，SWMM模型采用的版本为当前最新的5.1.013版本。将本发明方法计算结果与采用SWMM模型计算的结果进行比较，结果如图4所示。由该图可以看出，本发明计算结果与SWMM计算结果基本一致。另外值得一提的是，由于本发明采用的方法可以对管段任意分段，所以对管网水流的描述更加细致，对管段3水流到达时间捕捉的更加准确。说明本发明提供的处理排水管网水流模拟的方法是成功的。

上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。

Claims (4)

1.一种基于Godunov格式的城市排水管网水流数值模拟方法，其特征在于：根据管段首末端单元水位或压力水头与内节点的水位关系确定内节点处边界条件，采用Godunov有限体积格式模拟管道内水流运动，采用窄缝法处理管道内明渠流和压力流的变化过程；具体步骤如下：

1)获取计算范围内城市排水管网数据，包括管网平面几何拓扑数据、各管段的管径尺寸、各管段首末端的管底高程、各管网节点底高程、各管网节点的最大水深和超载水深、各管网节点面积数据；

2)对各管段和节点进行编号，建立各管段和节点的对应关系，设定水流从节点底高程高的一端流到节点底高程低的一端为正方向，采用有限体积单元离散各排水管段，每条管段的离散单元沿正方向依次编号；水力变量存储在各单元中心点；确定排水管网的入流边界和出流边界，入流边界和出流边界统称为外边界，施加外边界条件的节点称为外节点，其余节点称为内节点；

3)城市排水管网模型初始化：给各管段一维有限体积单元设置糙率和水力条件初始值；给管网内节点设置初始水力条件；

4)确定计算时间步长dt：城市排水管网模型采用显格式的有限体积法求解，该时间步长获取受到CFL条件限制；获取当前时刻各管段外边界条件，当前时刻以t时刻表示；

5)采用管道明渠流和压力流统一的方程组描述管道水流运动，具体形式如下：

连续方程：

动量方程：

式中Q为断面流量，A为过水断面面积，g为重力加速度，t为时间变量，s_f为沿程阻力损失，其表达式为s_f＝(n²Q|Q|)/(A²R^4/3)，R为水力半径，n为Manning糙率系数；当管道内为明渠流时，B为水面宽度，H为水位；当管道内为压力流时，B＝gA/a²，a为压力波波速，H为压力水头；

6)根据内节点水位和与其相连的各管段的首末端单元的水位或压力水头的关系来判断使用不同的内节点边界条件，从而计算流入或流出内节点的水量，具体过程如下：

a)计算与内节点连接的上游管段流入内节点的流量时，情况如下：

如果上游管段最后一个单元有水，并且内节点水位值小于该单元底高程值，则此时按自由出流边界条件计算进入内节点的流量；

如果上游管段最后一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算进入内节点的流量；

如果上游管段最后一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算进入内节点的流量，此种情况下，水流倒流，进入内节点的流量为负值；

如果上游管段最后一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算进入内节点的流量，此时进入内节点的流量为0；

b)计算与内节点连接的下游管段流出内节点的流量时，情况如下：

如果下游管段第一个单元有水，并且内节点水位大于、等于该单元底高程值，则此时按水位边界条件计算流出内节点的流量；

如果下游管段第一个单元有水，并且内节点水位小于该单元底高程值，则此时按自由出流边界条件计算离开内节点的流量，此种情况下，水流倒流，流出内节点的流量为负值；

如果下游管段第一个单元没有水，并且内节点水位大于该单元底高程值，则此时按照水位边界条件计算离开内节点的流量；

如果下游管段第一个单元没有水，并且内节点水位小于、等于该单元底高程值，则此时按照固壁边界条件计算离开内节点的流量，此时流量为0；

7)通过t时刻各管段单元的水力要素值、外边界条件以及内节点边界条件，采用HLL显格式计算通过单元界面处的数值通量，计算获得各管段t+dt时刻各单元的水力要素值；如果新计算出的某单元的水深值大于管径，则水流变为压力流状态，此时将管段水面宽度设定为窄缝宽度；

8)建立各内节点水量平衡方程，并根据该平衡方程获得t+dt时刻各内节点的水位值；

式中，

表示与内节点k相连接的第i个上游管段流入内节点的流量，

表示与内节点k相连接的第j个下游管段流出内节点的流量，A_k为第k个内节点的面积，

为第k个内节点在t时刻的水位，

为第k个内节点在t+dt时刻的水位，m为与第k个内节点相连接的上游管段数，n为与第k个内节点相连接的下游管段数；q_k为第k个内节点与地面的水流交互量，当由内节点流到地面时，q_k为负值，反之为正值；

9)令t＝t+dt，重复步骤4)～8)，直到计算结束。

2.根据权利要求1所述的一种基于Godunov格式的城市排水管网水流数值模拟方法，其特征在于：步骤4)dt的选取受到CFL条件限制，具体如式(4)所示：

式中：u为流速；

c为波速，g为重力加速度，A为过流面积，B为水面宽度，当水流变为压力流后，B采用窄缝宽度；Δx为单元空间步长；dt为时间步长。

3.根据权利要求1所述的一种基于Godunov格式的城市排水管网水流数值模拟方法，其特征在于：步骤5)当管道内水流为压力流状态时，水面宽度设定为窄缝宽度，窄缝宽度取值范围为B＝0.01m～0.001m。

4.根据权利要求1所述的一种基于Godunov格式的城市排水管网水流数值模拟方法，其特征在于：步骤6)中通过内节点边界条件计算流入或流出内节点的流量是通过如下方法完成的：在上游管段的最后一个单元之后或下游管段的第一个单元之前设置一虚拟单元，该虚拟单元的特征跟与其相连的单元特征一致，根据不同的内节点边界条件类型设置虚拟单元的水力变量，步骤6)存在如下几种情形：自由出流边界条件时，虚拟单元的流量和水位等于与其相邻单元的流量和水位；水位边界条件时，虚拟单元的流量等于与其相邻单元的流量，水位等于内节点水位；固壁边界条件时，虚拟单元的流量等于与其相邻单元流量的负值，水位等于与其相邻单元的水位；确定好虚拟单元水力变量后，根据HLL近似Riemann解计算通过虚拟单元和与其相邻单元交界面处的质量通量值和动量通量值，该质量通量值即为流入或流出内节点的流量。

Images