水平-竖直管道内水流冲击滞留气团的三维数值模拟研究

卢坤铭

(中国三峡建工(集团)有限公司，四川 成都 610041)

0 引 言

实际输水管道系统的启动填充过程中常会发生水流冲击滞留气团的复杂瞬变流，极易引起异常压力波动，影响系统安全运行甚至导致爆管事故[1-3]。

该瞬变流现象引起了很多学者的广泛关注。Martin[4]首次建立了水流冲击滞留气团现象的刚性水体数学模型，但其模型没有考虑水气交界面的运动。Izquierdo等[5]对Martin[4]的刚性数学模型进行了改进，考虑水气交界面位置动态变化，建立充水排气的刚性数学模型，并指出水体间滞留气团的快速压缩会引起异常压力升高。Liou等[6]针对起伏管道系统中初始时上游阀后完全放空的充水过程，建立了相应的刚性数学模型，但未考虑滞留气团的影响。刘德有等[7]建立了起伏变特性管道系统中水流冲击单个气团的刚性数学模型。刚性模型具有简单、快捷的优点，但其应用具有一定的局限性，如滞留气团含量较小时，可能会得到错误结论。Chaiko等[1]考虑水体和管道的弹性，建立了三种弹性数学模型，并给出了各弹性模型的适用范围。其他学者[8-13]研究了水流冲击滞留气团的瞬变压力及参数变化规律，发现水气耦合作用机理复杂，且可能引起十倍于入口压力的危险峰值压力；同时，提出了“刚性塞”、“虚拟塞”等方法，建立了弹性水体模型。

本文针对水平-竖直管道启动填充过程中含滞留气团瞬变现象，考虑了水体弹性、气体可压缩性、水-气交界面的动态运动，建立了相应的三维CFD数学模型。通过模型计算与实验结果对比，验证了所建模型能够有效地预测快速充水过程中含滞留气团的压力变化；同时，研究该瞬变过程中气团长度、湍流模型、管壁粗糙度等参数动态变化规律。

1 实验研究

针对水平-竖直管道内水流冲击滞留气团的瞬变流现象，设计搭建了如图1所示的实验系统。整个系统主要由潜水泵、气压罐、电磁流量计、球阀、管道、空气阀等组成。整个系统从上游至下游依次由蓄水池、不锈钢多级潜水泵、螺纹式球阀、气罐(压力罐)、电磁流量计、球阀、进气孔口、泄水阀、完全封闭的末端组成。水泵与气罐之间通过不锈钢钢管连接，压力罐至下游由一段1 m长的不锈钢和多段有机玻璃透明管道组成。从气罐出口至管道末端为水流冲击滞留气团的实验研究管道，其总长为8.862 m，其中，有机玻璃管道内径4 cm，壁厚1 cm。本实验中，将水平管道中心线定为Z=0，共安装8个压力传感器(Pressure Transducer，PT)和1个温度传感器，各装置位置参数如表1所示。在最底部分别安装2个排水阀，仅用于调节初始状态下气团长度，但在瞬变过程中均处于关闭状态。各压力传感器性能一致，测量范围：0～0.5 MPa。电磁流量计性能参数，公称通径为DN40，流量范围0～25 m3/h。数据采集系统为美国国家仪(NI有限公司产品PCI-6221，采样率250 kSps/s)。

表1 水平-竖直管道传感器位置 m

水流冲击管道内含滞留气团的动态特性与气团长度、湍流模型、管壁粗糙度、入口压力等参数有关，对于水平-竖直管道与各种参数的相互关系研究很少。表2给出了水流冲击水平-竖直管道内含一段滞留气团的模拟工况参数，分析了系统压力峰值、波动周期与各参数的相互关系规律。初始时所有阀门关闭，球阀快速开启以实现水流冲击气团瞬变过程，球阀开启时间0.1 s。

表2 水平-竖直管道模拟工况参数

2 数值模拟

2.1 控制方程

采用VOF模型追踪水气交界面，三维模型的基本控制方程为

连续性方程

(1)

能量方程

(2)

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

(3)

数值模型中，考虑了水气两相间的热传导、流体与固体壁面间的对流以及固体壁面内的热阻等可能存在的热效应，传导热流由傅里叶定律公式确定，即

q″cond=-λ∇T

(4)

式中，λ为介质(空气、水或壁面)的热传导系数，对于气相和液相之间的热传导用λeff代替λ。

对流传热项由牛顿冷却定律确定，即

q″conv=h(Ts-Tf)

(5)

式中，h为局部对流传热系数；Ts为壁面温度；Tf为流体温度。

壁面传导热阻公式为

(6)

式中，e为壁面厚度；kwall为壁面热传导系数。

2.2 湍流模型

本文对比了Standardk-ε、RNGk-ε这两种湍流模型计算结果。两种模型都适用于湍流，分子之间的黏性可以忽略，区别在于计算湍流黏性的方法不同以及控制湍流扩散的普朗特数不同，Standardk-ε湍流模型适合完全湍流的流动过程模拟，RNGk-ε则针对高雷诺数流动问题。

2.3 数值求解格式及边界条件

该模拟选用三维双精度求解器。采用有限体积法(FVM)进行网格离散，速度、压力和密度的耦合方法选择PISO算法，压力项采用PRESTO格式，动量项、能量项采用二阶迎风格式，体积分数采用Geo-Reconstruct格式进行离散，时间差分采用一阶隐式格式，残差收敛值为10-6，时间步长为0.000 1 s，计算时长3 s。

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

2.4 网格无关性分析

本次数值模拟采用六面体结构化网格，其可以更好地适用于复杂的流动边界。针对工况1，采用3种不同数量网格进行数值模拟计算并与实验结果对比分析验证，得出计算耗时与压力误差。通过改变网格节点数获取不同网格数量，分别为26.2万、33.2万、53.8万，对应的压力峰值与实验结果的误差分别为0.7%、0.45%和0.2%。基于模拟精度与计算时间的考虑，最终选用网格数量为33.2万的网格。

3 参数分析

3.1 气团长度的影响

图2给出了表2中所有工况的数值模拟结果和实验结果对比。由图2可知，模拟结果的压力幅值与波动周期都与实验结果吻合较好，特别是第1峰值结果与实验几乎完全吻合。对于第2压力峰值，随着气团长度的增加，模拟结果与实验结果的误差越来越小。原因可能是因为在数值模拟过程中采用了光滑壁面，同时忽略了法兰、流量计等的局部损失，导致计算的摩阻比实际结果小，而这一损失在流体流过管壁时显然不同。

水流快速冲击气团的过程中，能量耗散主要以热传导、水力机械损失为主。当初始气团长度较长时，可以观察到在水平-竖直管道的弯管处有非常剧烈的水气两相作用，水气交界面呈现不规则变化，这将导致能量主要以热量损失为主，热量耗散占总体能量损失的比重较大，而粘性耗散相对较小。当气团长度较短时，因为气团压缩膨胀路径较短，气-水交界面在竖直管道内趋于稳定，且在整个过程中几乎垂直于管道中心线，这表明对于该工况的能量耗散中以热量损失的比例降低，主要以摩阻损失为主。同一入口压力下，随气团长度减小，整个系统压力峰值增大，且增值随着入口压力的增大而逐渐减小，说明长气团工况对整个系统起到一定的缓冲作用。由图2可知，同一气团长度下，不同入口压力气团压力峰值之差分别为169.23、165.24、149.51 kPa，远大于入口压力增值40 kPa，说明入口压力越大系统压力峰值越大，且峰值增值远大于系统入口压力增值。

图2 不同气团长度情况下，水平-竖直管道内PT1处压力计算结果与实验结果对比

3.2 湍流模型的影响

图3给出了工况1、2、3、6时，不同湍流模型下，数值模拟结果与实验结果的对比结果。

图3 不同湍流模型下，水平-竖直管道PT1处压力计算结果与实验结果对比

由图3可知，4种工况在两种不同湍流模型下的数值模拟结果与实验结果都十分吻合，但是仍与实验结果存在一定的误差，为了能够准确地选择湍流模型进行计算分析，采用一定的数据统计方法对各工况进行分析，相对误差采用下式计算：

(7)

式中，Hexp为实验工况下最高压力；HCFD为数值模拟下最高压力。

由图3a、3b、3c对比可知，在同一较小入口压力(Pr=0.08 MPa)作用下，3种工况在Standardk-ε湍流模型下的第一压力峰值相对误差分别为1.99%、0.05%、4.82%，在RNGk-ε湍流模型下的第1压力峰值相对误差分别为7.64%、4.43%、4.02%。可以明显看出对于较短滞留气团，Standardk-ε湍流模型能够更准确地模拟第1压力峰值。对于气团较长的工况3而言，两种湍流模型的相对误差都大于4%，Standardk-ε湍流模型的计算结果低于实验结果，RNGk-ε湍流模型的结果则高于实验值。可能的原因是由于在低入口压力较长气团工况下，整个管道的紊流发展不充分，系统能量的损失主要以摩阻损失为主，但是两种湍流模型的应用场景更多的是充分发展的湍流。从波动周期分析，Standardk-ε湍流模型下3种工况的模拟结果均与实验吻合较好。由图3a、3c对比可知，Standardk-ε湍流模型能够很好地模拟前两个波动周期，从第3个周期开始，工况1滞后实验值，工况3提前实验结果，这是因为气团较短，水气耦合作用增强，气团运动速度加快，而气团较长时，气团压缩、膨胀过程需要更长时间，导致波动周期滞后。

由图3c、3d对比可知，在同一较长气团(La0=1.058 m)不同压力作用下，2种工况在Standardk-ε湍流模型下的第一压力峰值相对误差分别为4.82%和6.8%，在RNGk-ε湍流模型下的第一压力峰值相对误差分别为4.02%和2.33%。可以明显看出对于较长滞留气团，RNGk-ε湍流模型能够更准确地模拟第一压力峰值，且入口压力越大，Standardk-ε湍流模型的误差就越大，其误差已经超过5%，可能的原因是入口压力越大，随着球阀的渐进式开启，较小开度时在管壁方向会产生射流，开度逐渐加大，射流逐渐充分发展，速度越来越大，使得滞留气团在水平和竖直管道内的速度分布不均，气团在轴向和径向的运动速度不同，产生不对称压缩现象，这非常符合RNGk-ε湍流模型的高应变率和大曲率修正，因此RNGk-ε湍流模型的模拟结果更好。从波动周期分析，2种工况中RNGk-ε湍流模型的第一压力响应时间比Standardk-ε湍流模型分别提前0.032、0.048 5 s，且RNGk-ε湍流模型与实验结果的误差都在1%以内，入口压力越大，波动周期提前时间更长，与实验结果吻合更好。

3.3 管壁粗糙度的影响

图4给出了不同粗糙度下，工况4和工况6中PT1压力计算结果和实验结果对比。由图4a、4c可知，随着管壁粗糙度增加，两种工况第1、2压力峰值与实验结果的差值逐渐增大，但两者还是存在一定差异，分析如下：对于工况4，当管壁粗糙度为0.01、0.1 mm时，第1、2压力峰值的数值模拟结果均高于实验值，且由于管壁粗糙度引起的能量耗散可以忽略不计，当粗糙度为0.5、1 mm时，压力峰值开始出现较大衰减，且随着粗糙度的增加衰减值也逐渐增大，第2压力峰值的衰减程度大于第1压力峰值。图4b为粗糙度对第1、2压力峰值的敏感性分析，随着粗糙度的增加，第2压力峰值对粗糙度的敏感性高于第1压力峰值，其平均衰减程度为10.15%，而第1压力峰值的平均衰减程度为5.92%，两者接近两倍关系，说明对于小气团工况，增加壁面粗糙度高度，主要作用于第2个压力峰值。

图4 不同粗糙度，水平-竖直管道PT1处压力计算结果与实验结果对比

对于工况6，类似于工况4，当管壁粗糙度小于0.5 mm时，第1压力峰值数值模拟结果高于实验结果，平均高8.05 kPa，而第2压力峰值的结果与实验值吻合较好，因此管壁粗糙度引起的能量耗散可以忽略不计。随着粗糙度逐渐增大至0.5 mm和1 mm时，两个峰值均开始出现压力衰减，且随着粗糙度的增加衰减值逐渐增大。由图4d敏感性分析可知，第1压力峰值对粗糙度的敏感性高于第2压力峰值，因此第1压力峰值的衰减值大于第2压力峰值衰减值。第1压力峰值平均衰减10.44%，衰减值为33.72 kPa，第2压力峰值平均衰减9.5%，衰减值为23.41 kPa，说明对于大气团工况，增加壁面粗糙度会显著降低两个峰值压力，且对第1压力峰值影响更大。其原因是因为管壁粗糙度增大壁面剪切应力，使得在能量耗散过程中热传递的比例所占比例不同。

4 结 论

(1)建立水流冲击滞留气团瞬变流的三维精细化模型。充分考虑了水气两相压缩性、水-气-管壁热交换、管壁粗糙度等复杂因素进行三维建模，三维CFD模拟结果与实验结果对比验证了所建三维模型能准确地模拟瞬变压力及水气动态变化，可对水气耦合瞬变流的进行精细深入研究。

(2)研究气团长度、初始压力对系统压力峰值、压力衰减(能量耗散)的影响。系统最大压力随气团长度减小而增大，随入口压力增大而增大，且增值远大于系统入口压力增值。当初始气团较长时，系统能量以热量损失为主；当初始气团较短时，热量损失的比例降低，主要以摩阻损失为主。

(3)湍流模型是三维CFD建模的核心部分之一，影响管路系统的阻力、能量衰减的核心因素。对比发现，当气团较长时，无论入口压力大小均可以选择RNGk-ε湍流模型；当气团长度较短，入口压力较小时，可以选择Standardk-ε湍流模型，当入口压力较大时，则可以选择RNGk-ε湍流模型。

(4)管壁粗糙度对瞬态压力的影响和气团长度有关，总体趋势保持一致，即随着管壁粗糙度增加至0.5 mm时，系统压力峰值逐渐降低，粗糙度增加越大，压力衰减值越大。对于小气团工况，压力衰减主要作用于第二压力峰值；对于大气团工况，粗糙度会增加第一和第二压力峰值的敏感性，且对第一压力峰值的敏感性更高。粗糙度对压力波动周期的影响可以忽略不计。