管道中滞留气团排放过程的二维数值模拟研究

A 2-D Numerical Simulation Study on the Process of Air Mass Emission in Pipeline

摘要: 在极端暴雨条件下，城市雨洪排水管道中常会发生含滞留气团的瞬变流动。针对有压水流作用下滞留气团运动、排放的动态过程，本文采用VOF (Volume of Fluid Model)模型和k-ε湍流模型进行了二维数值建模和计算，并将模拟结果与已有实验观测结果进行对比分析。结果表明：VOF模型模拟出的滞留气团的运动过程、瞬态压力变化结果均与实验结果吻合。同时，对二维模拟时，主支管连接中可能存在的“等截面比例”和“等管径比例”两种简化方法进行了比较分析，比较结果表明，“等截面比例”方法计算的压力结果更接近实验值。

Abstract: Under severe rainstorm conditions, transient flow of trapped air masses often occurs in urban storm water drainage pipelines. In this paper, VOF (Volume of Fluid Model) model and k-ε turbu-lence model are used to simulate and calculate the dynamic process of the movement and emission of the stranded gas mass under the action of pressure flow. The simulation results are compared with the experimental results. The results show that the motion process and transient pressure variation of the trapped air masses simulated by the VOF model agree well with the experimental results. At the same time, the two simplified methods, such as “equal section ratio” and “equal diameter ratio”, are compared and analyzed in two-dimensional simulation. The comparison results show that the pressure result calculated by the “equal section ratio” method is closer to the experimental value.

1. 引言

暴雨情况下，城市雨洪排水管道将处于满负荷工作状态，其过程中会混入大量的有压空气。滞留气团的存在可能会引起管道堵塞，降低排水系统的总体容积和排水效率。当滞留气团经过竖井的时候，剧烈的气水相互作用甚至会引发竖井处出现爆发性的“间歇泉”现象，冲翻井盖、污水出井口，其对应的异常水锤极有可能导致管道系统的破坏，对城市的道路安全及公共卫生产生较大的影响 [1] 。此外，排水管道系统错综复杂的结构，也将使得城市排水管道中的气液两相流问题变得更为复杂。

随着城市内涝现象的愈发严重，越来越多的学者开始关注城市的排水防涝问题。Jing et al. [2] 就城市间歇泉问题建立了相应的实验系统，从是否有外部压力、管径比值方面，进行了多组实验，证明了外部压力对间歇泉发生的必要性。证实了小竖管直径和体积气团更容易产生间歇泉喷射现象。Muller et al. [3] 采用大尺寸实验装置进行实验，叠加背景流条件，观察到了更细致的喷发情况。刘德有等 [4] - [9] 、Zhou et al. [4] - [9] 、马佳杰等、郭艳惠等 [4] - [9] 针对封闭管道系统中含滞留气团的瞬变流现象进行了大量的实验及数值模拟研究。

长久以来，针对管道系统中涉及滞留气团的复杂气液两相瞬变流，其计算数学模型大多属于一维模型，存在大量的简化。其三维数值模拟，因其漫长的计算时间，相关的研究较少。一些学者 [8] [10] [11] 为了减少计算时间，针对该类型的流体运动问题，在简化一些次要的影响因素后，采用二维数值模拟的方法来还原实验过程。

本文参考Zhou et al. [4] - [9] 封闭式滞留气团的研究成果，对开放式的气团排放过程进行二维数值模拟研究，分析有压气体的释放过程及爆炸原因。并将数值模拟结果与Jing et al. [2] 的实验结果进行对比，验证所提二维模型的准确性。同时也对比分析了在管道二维建模过程中，“等直径比例”和“等截面比例”两种简化方法对计算结果的影响。

2. 数学模型

2.1. 控制方程

本文采用的是VOF模型，这是建立在固定的欧拉网格下的界面跟踪办法，适合使用在几种不相融的流体混合流动的情况。通过求解某一相体积分数的连续性方程来完成各相之间的界面的跟踪。该气液两相流数值模型求解采用k-ε湍流模型方程，气体是可压缩的理想气体，液体为密度998 kg/m³的水。其基本方程如下：

连续性方程：

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0$ (1)

动量方程：

$\frac{\partial \rho u_i}{\partial t}+\frac{\partial \rho u_i{u}_j}{\partial x_j}=\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+\rho g_i$ (2)

容积比率方程：

$\frac{\partial {\alpha }_w}{\partial t}+v\nabla \cdot {\alpha }_w=0$ (3)

式中：ρ是单元平均密度，g_i是i方向上重力体积力。P是静压强，x_j是坐标分量，v和u_i是速度和其分量，τ_ij是应力张量

能量方程：

$\frac{\partial \rho E}{\partial t}+\frac{\partial }{\partial x_j}\left[u_i\left(\rho E+p\right)\right]=\frac{\partial }{\partial x_i}\left(k_{eff}\frac{\partial T}{\partial x_i}\right)$ (4)

式中：E是能量，T是温度，k_eff是有效热传导系数。

对于气水两相流动，水相体积分数的连续性方程为：

$\frac{\partial }{\partial t}\left({\alpha }_w{\rho }_w\right)+\nabla \cdot \left({\alpha }_w{\rho }_{w}U\right)=0$ (5)

式中：α表示体积分数；ρ表示密度；下标w、a分别表示水相、气相。气相的体积分数α_a将受下面方程约束：

${\alpha }_w+{\alpha }_a=1$ (6)

在两相系统中，每个单元的平均密度为：

$\rho =\left({\alpha }_w{\rho }_w\right)+\left(1-{\alpha }_w\right){\rho }_a$ (7)

式中，ρ为气–水混合物密度。

湍流模型采用k-ε模型，其控制方程为：

$\{\begin{array}{l}\frac{\partial \rho k}{\partial t}+\frac{\partial \rho k{u}_i}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b+\rho \epsilon -Y_M\\ \frac{\partial \rho \epsilon }{\partial t}+\frac{\partial \rho \epsilon u_i}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\frac{\rho {\epsilon }^2}{k}\end{array}$ (8)

湍动粘度μ_t可以用k-ε的函数表示：

${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$ (9)

式中：k为湍动能；ε为湍流耗散率；σ_k和σ_ε分别是k和ε的紊流普朗特数，μ是粘性系数，G_k是由平均速度梯度产生的湍流动能，G_b是由浮力产生的湍流动能，Y_M是由可压缩湍流脉动扩散的。方程中的各常数均使用经验值，值取为C_1ε = 1.44，C_2ε = 1.92，C_μ = 0.09，σ_k = 1.0，σ_ε = 1.3。

2.2. 求解方法

采用有限体积法分离式迭代(Segregated Method)的方法求解二维流场，压力速度耦合采用PISO算法。

2.3. 边界条件

1) 进口边界条件设置。进水口断面为压力进口，在实验过程中为系统提供稳定的水压。

2) 出口边界条件设置。出口边界为压力出口，相对压力为零。

3) 管道壁条件设置。壁面为光滑壁面，粗糙度为0。

3. 二维建模

3.1. 算例介绍

为了验证所提二维数值模拟方法在模拟排水管道系统中滞留气团排放过程的有效性，选取Jing et al. [2] 设计搭建的排水管道实验系统，其简化的实验装置示意图如图1所示。整个系统主要由水箱、水平管道、竖管、阀门及压力传感器组成。水平横管总长约6 m，内径D = 0.05 m。竖管高度为1.8 m。上游水箱的尺寸0.51 × 0.63 × 0.69 m。4个PVC直角回转球阀用于将管道分离成不同的部分。管道的下游安装了一个可拆卸的盖子，允许在实验后排水。水平横管尾部和竖管底部各有一个压力传感器，用于监测实验压力数据。

实验开始前，系统水位注至H₀，阀门后为相应体积的空气，气压为一个大气压。阀门打开的同时传感器开始记录实验数据，直至气体从竖管排出。

3.2. 二维建模的简化方法

目前，研究城市排水系统中间歇泉现象的实验方法，主要过程是令恒定体积的空气从水平管道进入到竖直管道并且从竖管口释放。其中，竖管与水平管的直径比(d/D)是一个重要的研究参数，该比值通常小于1。

建立二维模型时，需要将立体的三维实验装置简化到二维平面。在此过程中，对于管网系统，需解决主支管连接处不同的直径比的问题。为此，针对图1所示的管道系统，提出两种简化方案：第一种是按原模型横竖管的直径比直接建立二维模型，称为等直径比例简化；第二种是是按原模型横竖管的截面积比例建立二维模型。各自的简化公式分别为：

等直径比例简化：

$\frac{a}{b}=\frac{d}{D}$ (10)

等截面比例简化：

$\frac{ax}{bx}=\frac{\text{π}{d}^2/4}{\text{π}{D}^2/4}\Rightarrow \frac{a}{b}=\frac{d^2}{D^2}$ (11)

式中：a是二维模型竖管直径，b是二维模型横管直径，d是原实验装置竖管直径，D是实验装置横管直径。

3.3. 二维模型

本文选取Jing et al. [2] 中的两组工况A、B建立二维模型，每种工况根据既定的两种简化方案，分别建立四种计算模型(Ad, As, Bd, Bs)，计算后所得数据与物理实验数据进行对比。四种工况参数见表1。

网格划分采用结构网格技术，在竖管处和横竖管交接处加密，以提高两相流流场计算精度，如图2

Figure 1. Physical experimental apparatus

图1. 物理实验装置

Figure 2. Model grid

图2. 模型网格

Table 1. Four calculation model parameters

表1. 四种计算模型参数

所示。经过网格无关性验证后，得出四个模型的网格数量控制在两万以上模拟效果较好。

4. 计算结果及数据分析

本文旨在：1) 验证二维CFD方法模拟排水管道系统中滞留气团排放过程的可行性；2) 探索滞留气团排放过程中的动态特性；3) 对比两种二维简化方法的模拟效果。首先，分析模拟结果的气液相云图，分析其动态特性。其次，对比二维计算的压力结果与实验结果。最后，对比两种简化方法的模拟效果。

4.1. 滞留气团排放过程的动态分析

对A、B两种工况的滞留气团排放过程进行二维数值模拟，将计算结果与实验数据进行对比分析。

如图3所示，滞留气团的排放可以分为两个阶段，气团在水平管的横向移动(a)和在竖管的纵向移动(b)。A、B两种工况中，气团横向运动到竖管底部一般是8 s左右，在该过程中首先是气团受到有压水体冲击发生的压缩形变过程，在经过初步的气水混合后，开始向竖管方向运动。气体进入竖管后，先是以

Figure 3. Working condition (Bd) simplified gas-liquid phase cloud diagram with constant diameter ratio

图3. 工况(Bd)等直径比例简化气液相云图

段塞流或是弹状流的流态向上运动，气泡逐渐加速，把水体向上推举。段状气泡逐渐融合形成环状流继续向上运动。达到一定高度和速度后，或是水体被推出竖管口，或是气团破裂并且气水高度混合形成雾状流向外喷射。本文从上述的两个阶段对计算数据进行对比。

本文数值模拟所计算的动态过程和结果与Jing et al. [2]已发表论文的实验结果相比较，过程基本一致，也验证了本文数值模拟的正确性。

4.2. 瞬态压力分析

图4、图5、图6、图7为工况(A)、(B)中的在p1、p2两个监测点记录的压力数据。每张图有三组数据进行对比——实验数据、等直径简化模拟数据(d)、等截面简化模拟数据(s)。对比每张压力图的实验数据和等截面模拟数据，结果表明：两组工况数值模拟所得的结果与实验结果基本吻合，压力波动趋势基本一致。尤其在第一峰值，非常接近。气团运动到竖管底部的时间大约为8.4 s，与实验误差在0.8 s以内。压力波动经过第一峰值后与实验数据在峰值上偏小，但周期上完全吻合。8 s后滞留气团进入竖管时造成的剧烈波动比实验值提前一些(1 s)。从上述情况可得结论：二维VOF模型模拟所得的结果与实测结果压力变化过程曲线基本接近，变化趋势相似度高，该模型可用于该类问题仿真计算，可信度高。

4.3. 两种简化方法对比及分析

从图4、图5、图6、图7中可以看出对比后可以发现，在第一峰值左右，三组数据基本吻合，然后按截面比例简化(s)的压力曲线与按直径比例简化(d)模型压力曲线相比，跟实验数据吻合程度更高，虽然在峰值上衰减较快，但在周期上完全吻合。

从工况(A)的模拟可以看出，按截面比例简化方案接近实验数据，更为准确一些。同时在图6、图7中，对比工况(B)后可以发现相似的波形走向，直径比例简化模型的过了第一峰值后，与实验数据开始发生偏移。按截面比例简化方案依然与实验数据较为吻合。

同时从图中可以看到当气体在竖管中运动，即间歇泉过程进入第二阶段时，所有的压力曲线均发生了剧烈波动。按截面比例简化(s)模型在第二阶段的波动与实验相似，而按直径比例简化(d)模型的波动较缓。

从上述的数据图可以直观的看出，按截面比例简化的方案在模拟时更为合理。现实中，城市排水系统和实验装置，大多使用的是圆管。若是采用直径比例简化方案，把圆形管的直径作为二维管道的宽度，加上计算时采用的一个固定的深度，在模拟计算中，整个系统成了由固定深度的矩形管组装的“新的实验装置”。

物理实验竖横管截面比：

$\eta =d^2/D^2$ (12)

新装置的竖横管截面比：

$\mu =\frac{d\ast h}{D\ast h}=d/D$ (13)

由于排水系统中，竖横管的直径比通常小于1，甚至大多情况小于0.5，所以结果 $\mu >\eta$ 。在模拟中，采用直径比例简化会放大竖管的截面，增大流通能力，对模拟结果会产生较大的影响。增大流通能力后，会放大细管截面积，使气体排放更加迅速，总体的排放气柱变短排放时间变短。气体更快更迅速的排出，会时气水相互作用的效果变差，故而，“等直径比例”简化方法的仿真结果较为平缓。

Figure 4. Test data of pressure sensor #1 in working condition A

图4. 工况A的压力传感器#1的检测数据

Figure 5. Test data of pressure sensor #2 in working condition A

图5. 工况(A)的压力传感器#2的检测数据

Figure 6. Test data of pressure sensor #1 in working condition B

图6. 工况(B)的压力传感器#1的检测数据

Figure 7. Test data of pressure sensor #2 in working condition A

图7. 工况(B)的压力传感器#2的检测数据

在气体进入竖管后，压力的剧烈波动是由于气水充分混合后产生的压力震荡，这也发生爆炸式间歇泉的一个标准。在该阶段模拟的数据与实验数据有一定的偏差。这是由于气水的混合度较高，这种复杂的气液两项流对网格划分和计算要就很高。所以在间歇泉反应喷发的时刻，较难进行模拟，有待进一步研究。

5. 结论

采用二维数学模型，对排水系统中的排气问题进行模拟是可行的，其模拟结果与实验结果吻合良好。通过建模过程和模拟结果，可以得出如下结论：

1) 在二维建模时，“等直径比例”与“等截面比例”两种简化方法的对比表明：可以按实际管道的截面比例建立模型，该模型可以有效的控制模型中横竖管的流通能力，令模拟计算结果更接近实际值。按现实管道的直径比例建立的模型，会放大细管截面积，增大细管流通能力，从而影响计算结果。

2) VOF模型可用于模拟滞留气团的运动过程、瞬态压力变化结果等气水动态行为和瞬变压力，模拟结果完全满足研究需要，并且具有一定的参考价值。

基金项目

国家电网公司科技项目资助(合同号：SGBXSJJS1700007)；国家自然科学基金：城市雨洪排水系统中“间歇泉”现象的瞬态特性及发生机理研究。