水流冲击多段滞留气团的三维数值模拟

卢坤铭，周领，刘静

(河海大学水利水电学院，江苏 南京 210098)

实际工程中，输水管道的高程会根据地形条件发生变化，往往会出现气团滞留于管道中来不及排出.水流冲击滞留气团现象，不但会影响系统正常运行，甚至可能引发爆管事故，从而造成极其严重的经济损失，甚至导致人员伤亡[1-3].

目前，针对水流冲击滞留气团瞬变流现象，现有的模型大多为一维数学模型，且主要研究对象为单个气团，对于起伏管道内水流冲击多段气团的瞬变流研究，相关成果很少.MARTIN[4]首次对水流冲击多段滞留气团现象进行了理论研究，基于刚性水体理论，建立了相应的数学模型，但其模型没有考虑水气交界面的运动.刘德有等[5]建立了起伏变特性管道系统中水流冲击单个滞留气团的刚性数学模型.刚性模型具有简单、快捷的优点，但其应用具有一定的局限性，如：在复杂的管网系统中，刚性模型处理起来尤为烦琐和复杂，尤其是出现液柱分离的情况.EPSTEIN[6]将积分法用于水锤基本方程的求解，建立水流冲击两段封闭滞留气团的数学模型，由于没有相应的多段气团的弹性数学模型，该方法仅与刚性数学模型进行了比较，结果显示当气团长度很小时，刚性模型预测的气团压力过高.ZHOU等[7-9]研究了水流冲击滞留气团的瞬变压力及参数变化规律，发现水气耦合作用机理复杂，且可能引起10倍于入口压力的危险峰值压力；同时，提出了“刚性塞”、“虚拟塞”等方法，建立了弹性水体模型，成功避免了特征线法在水气交界面动态追踪时的复杂插值处理，并保持了同样的计算精度；指出在气团含量很小的情况下，刚性水体模型的不适用性.

现有的一维数学模型存在诸多假定，其计算精度往往不能满足实际需求.对于水气耦合两相瞬变流而言，整个过程往往很难通过试验手段准确记录.文中采用三维CFD方法模拟起伏管道内含多段滞留气团的瞬变流，应用Standardk-ε湍流模型，研究两段滞留气团的相互作用规律，并将三维计算结果、一维模型计算结果与试验结果进行对比分析；同时，研究整个过程中水-气作用规律.

1 试验装置

整个系统从上游至下游依次由蓄水池、潜水泵、气罐(压力罐)、电磁流量计、球阀、进气孔口、泄水阀、封闭末端组成.水泵与气罐之间通过不锈钢钢管连接，压力罐至下游由一段1 m长的不锈钢和多段起伏的有机玻璃透明管道组成.从气罐出口至管道末端为水流冲击滞留气团的试验研究管道，其总长为10.97 m，其中，有机玻璃管道内径4 cm，壁厚1 cm，具体如图1所示.

图1 含多个滞留气团试验装置布置图

试验中，将气罐出口和管道入口交界面处定为x= 0，水平管道中心线定为z= 0；图1中给出了管道弯曲段最高点和最低点处的沿线长度x和高程位置z.球阀距上游入口距离为x= 2.236 m；共8个压力传感器，其安装位置分别如下：PT-1#为(x=10.82 m,z=0.632 m)，PT-2#为(x=9.15 m,z=0.075 m)，PT-3#为(x=7.60 m,z=0.82 m)，PT-4#为(x=6.68 m,z=0.24 m)，PT-5#为(x=5.75 m,z=0.32 m)，PT-6#为(x=4.33 m,z=0.96 m)，PT-7#为(x=1.707 m,z=0 m)，PT-8#为(x=1.337 m,z=0 m)；在弯曲管道最顶部和最底部分别安装3个进气孔、3个排水阀，仅用于调节初始状态下气团段数和长度，但在瞬变过程中均处于关闭状态.各压力传感器性能一致，测量范围：0～1 MPa，线性误差和滞后误差范围为0.03%～0.30%(厂家提供).电磁流量计性能参数，公称通径：DN40，流量范围：0～25 m3/h.数据采集系统为美国国家仪器(NI)有限公司产品：PCI-6221，采样率250 kS/s.

试验初始条件：上游球阀完全关闭，球阀上游为有压水体与气罐相通，球阀下游为起伏管段，初始气团的个数和长度通过3个进气孔、4个排水阀进行调节实现的，管道下游末端始终处于封闭状态.水流冲击滞留气团瞬变过程是通过瞬间开启球阀实现的，整个瞬变过程中气团始终滞留于管道内.

表1为3种试验工况的初始状态，各物理参数：xf为起伏管道系统中，阀门开启后t时刻冲击水体的长度，m；x1u为第1 段阻断水体的上游端沿线位置，m；x1d为第1 段阻断水体的下游端沿线位置，m；x2u为第2 段阻断水体的上游端沿线位置，m；x2d为第2 段阻断水体的下游端沿线位置，m；La0,1为第1段气团的长度，m；La0,2为第2段气团的长度，m；La0为气团总长度.试验针对含2个初始滞留气团情况，进行了3种工况试验入口压力均为Pr= 0.16 MPa，初始冲击水体长度均为xf= 3.3 m，具体初始状态值见表1.为保证试验精确性，严格、精确控制试验的初始条件和边界条件；每种工况重复至少3次，直至试验结果基本一致.

表1 3种试验工况的初始状态值

为了便于数值模拟计算比较分析，试验通过高速摄像机的拍摄记录显示，试验中手动球阀的开启时间(从全关至全开)为0.07～0.09 s，通过纯水锤试验得水锤波速约为400 m/s.

2 CFD数值计算

2.1 三维物理模型及网格划分

通过Design Modeler对试验系统进行三维物理建模，导入ICEM划分网格，如图2所示.为了获得更加准确的计算结果，对弯管处、水箱下部、近壁面处进行加密处理，均采用六面体网格.整个管道的网格数为33.2万，最小网格质量为0.350，平均网格质量为0.872.

图2 三维模型及网格示意图

2.2 网格无关性验证

针对工况1，采用3种不同数量网格进行数值模拟计算并与试验结果对比分析验证，得出计算耗时与压力误差.通过改变网格节点数获取不同网格数量，分别为26.2万、33.2万、53.8万，对应的压力峰值与试验结果的误差分别为0.70%,0.45%和0.20%.基于模拟精度与计算时间的考虑，最终选用网格数量为33.2万的网格.

2.3 基本控制方程及湍流模型

基本控制方程包括容积比率方程、动量方程、能量方程以及水体状态方程.

容积比率方程为

(1)

动量方程为

(2)

能量方程为

(3)

式中：ρw为水体密度；αw为水体体积分数；t为时间；v为流体速度；ρ为水-气混合物密度；p为静压；μ为水-气混合物黏度；μw，μa分别为水、气相流体黏度；g为重力加速度；F为外部体积力；E为总能；T为温度；keff为有效传热系数.

模拟中，气体、水体均可压.假定气体遵循理想气体定律，水体状态方程为

(4)

(5)

其中：a为水锤波速；D和E为管径和管壁厚度；E为弹性管道杨氏模量.

湍流模型：通过对3种湍流模型(Standardk-ε，RNGk-ε和Realizablek-ε)进行适用性分析[10]，Standardk-ε模型能更准确地模拟流态变化过程.故文中采用Standardk-ε湍流模型进行数值模拟研究.

2.4 求解控制及边界条件

求解控制：瞬态计算中，采用有限体积法(FVM)对网格进行离散，时间步长为1.0×10-4s，所有参数收敛标准均为1.0×10-6，每一时间步最大迭代步数为100步，经后期观测发现迭代60次以内均能收敛.速度、压力和密度的耦合方式选择PISO方法，压力项采用PRESTO格式，动量项、能量项采用二阶迎风格式，湍动能项、湍流耗散项采用一阶迎风格式，体积分数采用Geo-Reconstruct格式.计算时长5 s.

边界条件：压力罐下部Patch为有压水体，阻断水体及气团均为大气压.阀门区域采用动网格模型，开阀时间为0.08 s.上游管道出口与阀门进口、阀门出口与起伏管道进口设置为2对“interface”边界条件，其他均为默认边界条件.

3 计算结果

3.1 瞬变压力三维CFD计算与试验结果对比

图3为不同气团长度下，PT-1#和PT-4#压力传感器处压力波动计算结果与试验结果.表2给出了PT-1#和PT-4#压力传感器处最大压力(Ha1-max,Ha2-max)的CFD值pCFD和试验值pE.与管道内含一段滞留气团结果相比，管道中含两段滞留气团增加了水力瞬变过程的复杂性.试验结果表明：① 系统的最大压力始终出现在靠近管道末端的滞留气团的活动区域；② 位于管道末端的第2段阻断水体(仅工况3中)，始终保持静止，内部压力处处相等，其与气团压力变化基本一致；③ 在水流冲击过程中，由于两段滞留气团间相互作用，导致气团的压力变化曲线出现峰值波动，其发生位置和剧烈程度与初始两段气团相对体积大小有密切关系.

表2 气团最大压力值的CFD值和试验值

图3和表2表明，2个气团的数值模拟结果与试验数据十分吻合，三维数值模拟可以很好地研究水流冲击多段滞留气团的瞬态压力.此外，水流冲击两段滞留气团压力波动曲线中，两气团的瞬变压力并非同步变化，可能呈现两气团峰值压力交替出现的情况，这与初始气团长度有一定的关系.如图3a所示，第2段滞留气团的压力峰值大于第1段滞留气团，初始的两段滞留气团的压力变化过程中，两者的峰值压力交替出现.如图3b所示，与工况1相比，工况2中初始第1段滞留气团长度增大，其压力变化变得平缓，而第2段滞留气团长度减小，压力变化较为剧烈，2种情况的压力曲线的总体趋势是一致的.图3c显示，与工况1相比，第1段气团长度保持不变，第2段滞留气团缩短91%，导致系统压力峰值增加23%，压力峰值响应时间缩短58%.3种工况峰值压力和初始气团总体含量对比表明，同一入口压力下，随着初始气团总体含量的增加，系统压力峰值相应减小.

图3 2个气团计算结果与试验结果对比

3.2 现有一维模型计算结果对比

图4为3种试验工况下，采用Standardk-ε湍流模型三维模拟得出的系统最大压力与一维模型计算结果、试验观测数据的对比结果.其中，一维模型采用本课题组已发表的弹性水体数学模型[10].

3种工况均表明三维CFD方法能够准确地模拟整个压力波动过程，其压力峰值与试验结果更接近.一维弹性水体模型可以简单模拟整个瞬变过程的压力变化趋势，与试验结果存在较大误差，其压力峰值均大于试验值.一维数值模拟结果表明，工况2在第一个波动周期后存在明显的压力振荡，说明第2段气团受到阻断水体的反复压缩.弹性水体模型假定：① 冲击水体和阻断水体间被气团完全隔开；② 滞留气团占据整个管道截面，水气交界面垂直于管道中心线，水、气互不混掺；③ 管道内滞留气团假定为完善气团，满足气体热力学多变过程方程.试验以及三维数值模拟过程中发现，水体运动至管道弯曲处，可能将气团分成若干部分，能量衰减加快，阻断水体长度发生变化，导致一维数学模型的计算结果与试验有所偏差，无法准确模拟整个压力过程.三维CFD结果明显好于一维模型结果，但与试验结果仍有差别，表现为从第二个周期开始提前于试验结果，可能的原因是很难真实地模拟水-气-管壁的热传递.

图4 一维-CFD-试验结果对比

3.3 动态特性分析

图5为工况1下水-气两相分布图.

图5 不同时刻水-气动态特性

结合试验观察结果，水流冲击管道内含多段滞留气团的瞬变流过程可描述如下：初始时，上游阀门完全关闭，阀门上游段为高压段，与压力罐相连，阀门下游段管道系统中滞留两段气团和阻断水体.当上游阀门瞬间打开之后，管道系统开始充水，靠近上游的第一段滞留气团在冲击水体的作用下开始压缩，其压力逐渐增大；第一段滞留气团压力增大到一定值时，开始推动阻断水体向下游运动.与此类似，第二段滞留气团也受到压缩，由于管道系统末端完全封闭，管道下部阻断水体撞击在管道末端，气团压力迅速增大并反向压缩阻断水体，随着水气不断混掺，最终阻断水体被分隔成两段，滞留气团则变为三段.

数值模拟结果发现，随着气团的压缩-膨胀-变形，水气交界面自由变化，并不垂直于管道中心线，同时，冲击水体和阻断水体互相掺混，当水体运动到管道弯曲处时会产生新的阻断水体，其长度发生变化，将气团分成若干部分.这与试验结果保持一致，并进一步解释三维CFD数值模拟结果优于一维数学模型.

4 结 论

1) 三维CFD模型能够准确地对起伏管道内水流冲击多段滞留气团的瞬变流现象进行可视化分析与处理，且与试验结果吻合较好.

2) 三维CFD模拟结果与试验观察结果显示，系统的最大压力始终出现靠近管道末端的滞留气团的活动区域.在水流冲击过程中，由于两段滞留气团间相互作用，导致气团的压力变化曲线出现峰值波动，其发生位置和剧烈程度与两段气团相对体积大小有密切关系.

3) 一维弹性模型计算结果误差大于CFD结果，这是因为一维模型无法有效地模拟水气交界面变化过程、阻断水体长度的动态变化.

4) 对水流冲击多段滞留气团进行动态模拟并分析误差.整个过程中，水气两相互相掺混，水气交界面自由变化，阻断水体长度时刻发生变化，当水体运动到管道弯曲处时会产生新的阻断水体，将气团分成若干部分.