雨水管道末端与河道水动力相互作用模拟研究

摘    要：本文利用开源计算流体力学软件OpenFOAM搭建二维数值模型研究管道-河道耦合水动力的作用机制，通过设置不同的管道上下游边界水位和管道坡度探究不同降雨时期河道水位和管道坡度对雨水管道末端排水能力的影响。研究发现下游河道水位上升不一定降低管道排水能力，在一定强度降雨条件和管道坡度下，河道水位淹没雨水管道末端一定程度时能促进管道提前形成满管流，进而提高管道排水能力，而且坡度小的管道更容易在更低的河水位下形成满管流。研究结果有助于更全面地揭示雨水管道末端与河道的水动力相互作用机制，为三维模拟与实验研究的开展提供依据，也可为城市雨水管道改造及河道水位调控提供指导。

关键词：城市内涝;雨水管道;耦合;数值模拟;OpenFOAM;

城市洪水是最常见和最具破坏性的自然灾害之一[1]。极端暴雨的频发、城市化进程的地面硬化、雨水管道的排水能力不足等都是导致城市内涝的重要因素。

在过去的几十年里，国内外许多研究者通过理论分析、模型试验和数值模拟的方法，对城市内涝的致灾机理及解决方案进行了大量研究。在极端气候方面，郑飞飞等[2]分析了澳大利亚69个雨量计在1966—2012年间的数据，结果表明极端降水在不同时间尺度上发生变化会对洪水风险产生巨大的影响；吕洪等[3]在对郑州进行城市洪涝灾害预测时，利用极端降雨事件的影响构建新的风险指标体系，得到了更全面准确的预测结果；M. Bermúdez等[4]在西班牙西北部沿海河段进行的测试表明，气候变化对河流流量的影响在具有复合洪涝潜力的地区具有非常重要的影响。在地表径流方面，陈洁云等[5]总结了城市雨水径流估算分析模型的发展，证明封闭式分析模型能高效准确地估算城市雨水径流量；Martin Bruwier等[6]系统分析了9个城市特征在洪泛情况下对地表径流的影响，表明洪水严重程度主要受建筑面积的影响；Francisco Peňa等[7]采用 FLO-2D 和MODFLOW-2005提出了一种基于物理、松耦合的建模框架，对降雨过程中地表水和地下水的相互作用导致的城市洪涝进行细化研究。

对于地表径流和管道流耦合的过程，Noh Seong Jin等[8]在城市淹没分析中研究了不同形状雨水口的流量系数；Ricardo Martins等[9]采用计算流体力学(computational fluid dynamics, CFD)和试验相结合的方法对不同排水条件下沟渠的水力特性进行了研究；陈倩等[10]对不同工况雨水口的泄流能力进行了试验研究；还有一些研究人员在城市洪水耦合模型中对进水口进行了不同的模拟分析，并通过试验进行了验证[11,12]。在管道流动方面，Heekyung Park等[13]早在1998年就使用水动力学模型对环状和树状管道系统进行了详细的水动力评估；Pachaly Robson Leo等[14]评估了在雨水管理模型5.1(Storm Water Management Model 5.1,SWMM)中加入Preissmann槽算法后，在复杂、高度动态流入场景中的性能；叶家强等[15]研究在管道建模过程中通过入渗效率来估算管道排水量，其中入渗效率和降雨高潮位都有关系。

在雨水管末端区域，河道水位的顶托效应对其排水性能具有重大影响。Félix L等[16]综述了前人多复合事件驱动因素耦合方法，并提出未来的研究应集中于开发一种紧密耦合的程序，以准确地反映风暴潮和降雨径流之间复杂的物理相互作用；张兆祥等[17]采用水动力模拟，研究淹没出流时水位、坡度、管径对内涝的影响，结果表明，河道水位上升严重制约雨水管道排涝能力；罗鸣等[18]也研究了河道边界水位对管道排水能力的影响，探究不同降雨重现期下出水口管道分别为自由出流、半淹没出流和完全淹没出流3种情形的管道相对排水能力变化。

针对城市内涝的成因，研究人员对降雨、地表径流、地表和管道耦合、管道流动、管道和受纳水体耦合的过程进行了细致的分析。尽管已有不少文献对下游高水位和极端降雨的复合事件造成的城市内涝进行研究，但是大部分文献仅是从宏观上描述高潮位会阻碍管道排水从而导致城市内涝，没有量化高潮位对管道排水能力的影响程度。且大多研究采用SWMM进行模拟，而SWMM在描述管道内水动力过程方面存在严重的不足，无法精准刻画下游顶托情况下管道内的水流复杂情况。

本文主要通过二维模型对末端雨水管道进行模拟，以初步探究不同的降雨强度、水位高度、管道坡度对管道排水能力的影响。采用二维模型的主要原因是其可以高效地捕捉到管道中的明满流交替现象，进而阐释雨水管道末端与河道之间的水动力相互作用趋势和规律，为复杂三维模拟及实验研究的开展提供依据。而且，本文的二维模型模拟结果通过经验公式验证精度可靠。

1 数学模型

1.1 控制方程

OpenFOAM中的interFoam求解器可以求解2个不可压缩相的雷诺平均Navier-Stokes方程，包括连续性方程(公式1)和动量守恒方程(公式2),并通过VOF(volume of fluid)方法捕捉自由表面的运动。在笛卡尔坐标系中，Navier-Stokes控制方程由以下几个部分组成：

式中U表示速度矢量，ρ是流体密度，p*是伪动压，g是重力加速度，X是位移矢量，μeff是有效动力黏度，μ是水或空气的动力黏度，σ是表面张力系数，κ是界面的曲率。此外，为封闭求解雷诺时均N-S方程，本文采用k-omega SST湍流模型。

对于VOF方法，必须求解一个额外的方程来描述相的运动。相函数α被定义为每个单元中水体积的比例，需满足一个经典的平流方程，如下式所示：

本文采用的VOF方法与压缩技术相结合，有助于抑制数值扩散，保留自由表面的锐度。为了实现上述期望，Weller[19]在公式(3)中增加了一个人工压缩项，最终表达式如下：

其中Ur是只在界面上工作的水相和空气相之间的压缩速度。此外，为了实现公式(4)的约束性，显式求解的多维通用限制器算法(MULES)用来确保α的计算结果总是约束在0和1之间。上式中，ρ和μ都是通过使用相位函数对水和空气进行加权计算：

ρ=αρwater+(1-α)ρair (5)

μ=αμwater+(1-α)μair (6)

1.2 数值方法

在求解控制方程(公式1～3)的过程中，采用了有限体积法(finite volume method, FVM)进行数值离散化。具体而言，计算域被离散成一系列的小网格，并使用Rhie和Chow[20]开发的拼合网格方法，将所有的流场信息都存储在各个网格的中心。

采用OpenFOAM开展研究的优势在于控制方程中具有多种不同项的离散化方法和插值方法。基于前人研究基础[21]:本文的瞬态项采用隐式欧拉方法，平流项采用Gauss limitedLinearV1方法，黏性扩散项采用线性校正方法；除上述之外，所有的剩余项都采用了线性内插法；在求解压力和速度耦合问题时，采用PIMPLE算法进行处理。PIMPLE是一种混合了压力关联方程(semi-implicit method for pressure linked equations, SIMPLE)算法的压力隐式算子分裂法(pressure implicit with splitting of operators, PISO)和半隐式方法，主回路继承自PISO。虽然本文中没有使用欠松弛，但其允许方程欠松弛以保证所有方程在每个时间步收敛[22]。

2 模型验证

2.1 数值模型

本文采用的二维数值模型如图1所示。由于数值模拟要求上下游水位边界水位尽量稳定，所以将水塔和水槽的长度设置得尽量长以满足模拟要求，如图1所示，水塔长度为350 m, 高度为1 m, 水槽长度为100 m, 高度为1 m。除此之外，排水管的直径D为0.1 m, 水平方向的长度L为8 m, 坡度i取1%、0.5%和0%,管道下缘距离水槽底部0.25 m, G1,G2和G3是流量监测面。这些参数的选择都是根据工程经验以及实验平台构建方案决定。

图1 二维数值模型示意图 

2.2 网格收敛性

对于基于有限体积法的数值模拟，其计算精度不仅取决于离散化方案和计算算法，还与网格质量和网格大小有关，因此需要对网格的灵敏度进行验证。

为了检查网格的灵敏度，首先比较了管道中4种不同网格大小的数值模拟结果。为了捕捉自由表面，水塔中从排水管上沿至上方0.7 m处为加密区域，网格的纵向尺寸为2×10-3 m; 同理，水槽中从排水管上沿至上方0.6 m处为加密区域，网格纵向尺寸为1.85×10-3 m。由于管道进出口处水流状态复杂，在水塔和水槽中靠近管道的水平方向设置10 m的网格加密段，该段网格水平方向的最小尺寸与管道水平方向最小尺寸相同。由于本数值模拟中主要关注的是管道部分，所以管道中需要有一个足够详细的网格尺寸来捕捉相关的流场特征。如图2所示，在研究中，考虑了4种不同离散度的网格，分别是(a)粗网格、(b)中网格、(c)细网格、(d)最细网格。管道纵向分别划分30、40、50、60份，最小网格尺寸分别为3.3×10-3、2.5×10-3、2×10-3 、1.7×10-3 m。网格的长径比最大(Δx/Δz)为5,最小为1,在靠近管道两端接口处最小。从两端到管道中间各有水平方向1 m的网格长度渐变区。数值模拟采用时间步长为1×10-3 s, 计算至80 s, 取稳定后的流量平均值作为管道排水量。图3比较了不同网格尺寸在测试工况下的管道排水量结果。粗网格和中网格的结果最大偏差接近1.01%,而对于较细和更细网格的结果偏差可以忽略不计。此外，细网格可以节省12%的计算量。因此考虑到计算精度和成本，在接下来的章节中，都采用了管道纵向网格尺寸最小为2×10-3 m的细网格。

图2 4种不同离散度的网格用于数值收敛性研究 

图3 不同尺寸网格的数值收敛研究结果 

2.3 经验结果对比

本文模型在管道自由出流情况下可类比闸孔出流，闸孔出流流量计算经验公式可参照Brater E F等[23]编制的水力学手册，如公式(7)所示：

式中，h是自由出流下闸孔中心点水头，a是闸孔长度，g是重力加速度，C=0.6是排放系数[24]。

为了描述方便，参考前人方法[18],针对排水管出口处的水位情况定义无量纲的淹没度S,S指的是水槽的水位和排水管出口下缘的高度之差与排水管直径的比值，具体如下：

其中，hc是水槽的水位；h0是排水管出口下缘的高度，h0=0.25 m; D是排水管的直径。

选取管道坡度1%,Δh=0.30 m(Δh是水塔和水槽的水位差),管道两端水位同时上涨的一组工况，对比数值结果和经验结果之间的差异。其中由于只能在下游顶托较低时将模型看作自由出流，进而采用经验公式计算流量值，故只计算S=0和S=0.25 2个工况下的经验流量进行对比。如图4所示，经验结果与数值结果很接近，前者略大是因为经验公式对应的是闸孔出流，闸孔后没有管壁存在，而二维模型包含管道管壁对水流的阻尼效应。根据以上分析推断本文中的二维模型可开展雨水管道末端与河道水动力相互作用规律的探析研究，模拟结果满足工程精度要求。

图4 数值和经验结果的比较 

3 排水管道-河道耦合水动力过程

3.1 降雨初期

降雨初期，上下游水位同时上升。在数值模拟中假设水塔和水槽的水位差保持不变以简化问题。降雨初期管道排水量随顶托变化的模拟结果如图5所示。

图5(a～c)展示了当坡度为1%(工程中常见的雨水管道坡度),不同Δh条件下管道排水量随S的增加先变大后趋于平稳。以Δh=0.20 m工况中管道流态的变化分析原因。如图6所示，S=0.25时管道未形成满管流，S=0.5时管道接近满管流，S=0.75时管道已呈满管流。管道中平均流速主要受Δh影响，Δh不变，S从0.25变化到0.75的过程中管道过流面积不断增加，因此排水量增加；当 S>0.75 时，管道中均为满管流，此时过水面积不变，在Δh相同的条件下，管道排水量基本保持不变。

图5(d～f)展示了不同坡度管道在Δh=0.30 m 情况下排水量随S的变化特征。i=0%和i=0.5%的管道在S=0.25后随顶托的增加排水量保持不变，i=1%的管道在S=0.75后排水量保持不变，说明管道坡度越小越容易在较低的河道水位下形成满管流。

图5 降雨初期管道排水量与S的关系 

图6 稳定后排水管入口处的水流状态 

3.2 降雨中期

降雨中期，排水管上游来水基本稳定，下游河道水位持续上涨。在数值模拟中假定水塔水位固定，水槽水位上涨。降雨中期管道排水量随顶托变化的模拟结果如图7所示，增加管道排水量随着下游顶托的呈现先增加后下降的趋势，这与工程中普遍认为的管道排水量随顶托增加而下降不同。

图7 降雨中期管道稳定排水量与S的关系 

图8 降雨末期管道稳定排水量与h水塔的关系 

图7(a～c)展示了1%坡度的管道排水量随S的增大先基本稳定，再突然增大，最后逐渐减小。在S=0到S=0.35之间排水量不变，这是因为随着顶托的增加，水位差减少，导致流速降低；但同时下游水位上升，造成过流面积增加，在这两者耦合影响下，最终流量变化不显著。在S=0.5左右管道形成满管流，此时排水量达到最大值，对于水塔水位为0.55、0.65 m和0.7 m的情况，管道排水量比自由出流时分别增加23%、35%和25%。在此之后，管道都是满管流，过流面积不变，顶托增加导致水位差减小，流速减小，流量减小，水塔水位为0.55、0.65 m和0.7 m时，在S=1.5的情况下排水量比排水量最大时分别减少20%、28%和17%。

水塔水位不变，不同坡度管道排水量随S的变化趋势如图7(d～f)所示。坡度为0的管道在S=0.05左右就已经达到排水量最大值，排水量比自由出流时增加8%,坡度为1%和0.5%的管道分别在S=0.25和S=0.4达到排水量最大值，排水量比自由出流分别增加17%和28%,这也证明了小坡度管道更容易形成满管流。

3.3 降雨末期

降雨末期，排水管上游水位下降，下游河道水位基本稳定。在数值模拟中假定水槽水位固定，水塔水位下降。图8展示了降雨末期管道排水量随水塔水位变化的数值模拟结果。从结果可知，降雨末期管道顶托不变，排水量随着水塔水位的降低而减小，这与普遍的认知一致。

4 结语

本文主要发现在降雨初期，顶托增加的同时上游来水量也在增大，导致管道排水量先增加，直到管道形成满管流后排水量趋于稳定；在降雨中期随着顶托的增加排水量会先达到一个最大值，达到排水量最大值的顶托程度和管道坡度有关；在降雨末期，顶托不变，管道排水量随水塔水位的降低而减小。总的来说，管道坡度越小越容易在较小的顶托程度下形成满管。研究结果可以为城市雨水管道改造及河道水位调控提供指导，对后续的三维模拟与试验研究的开展具有重要指导意义。

文章来源：科技通报