混流泵在船舶推进、防洪灌溉、城市给排水等领域都有广泛的应用，空化会导致混流泵性能下降、引起泵的振动和噪声，并腐蚀过流部件^[1-2]。空化数值模拟是当前国内外研究的热点之一，很多学者提出了不同的空化模型，如Kunz 空化模型^[3]、Singhal空化模型^[4]、Schnerr-Sauer空化模型^[5]和Zwart-Gerber-Belamri空化模型^[6](简称Zwart空化模型)等，其中Singhal空化模型和Zwart空化模型应用十分广泛^[7]。国内许多学者运用Zwart空化模型对水泵的空化进行了数值模拟，部分学者对于空化性能的预测未得到实验验证^[8-9]，部分学者计算得到的扬程-有效汽蚀余量曲线与试验结果趋势一致，但对临界汽蚀余量的预测误差较大^[10-11]。除空化模型外，空化数值模拟还要选择合适的湍流模型，RNGk-ε湍流模型(RNG湍流模型)和SST湍流模型在空化数值模拟中得到了广泛地应用^[12-13]，但是较少有人对比分析这两种湍流模型对空化数值模拟的影响。

本文运用Zwart空化模型对混流泵的空化性能进行了数值模拟，分析了RNG和SST两种湍流模型对计算结果的影响，并运用湍流粘度修正后的RNG湍流模型进行了定常和非定常空化数值计算，分析了泵内压力脉动的时域特征和频域特征。

1 空化模型

空化模型是描述气相与液相之间相互转化的数学模型，均相流模型认为气相与液相达到了动力平衡和热平衡，把气相与液相作为统一的流体进行研究^[14]。本文采用Zwart空化模型对混流泵的空化性能进行数值模拟，其形式为^[6]

$\dot m = {F_{{\rm{vap}}}}\frac{{3{\alpha _{{\rm{nuc}}}}\left( {1 - {\alpha _v}} \right){\rho _v}}}{{{R_b}}}\sqrt {\frac{2}{3}\frac{{{P_v} - P}}{{{\rho _l}}}} ,P < {P_v}$

(1)

$\dot m = {F_{{\rm{cond}}}}\frac{{3{\alpha _v}{\rho _v}}}{{{R_b}}}\sqrt {\frac{2}{3}\frac{{P - {P_v}}}{{{\rho _l}}}} ,P > {P_v}$

(2)

式中：F_vap为空化系数，取为50；F_cond为空泡凝结系数，取为0.01；α_nuc为汽核体积分数，取为5×10^-4；R_b为空泡半径，取为1×10^-6 m。在空化数值模拟时考虑湍流压力脉动P_turb的影响，对饱和压力P_v进行如下修正:

$$

(3)

式中：k为湍动能，m²/s²；ρ_m为混合密度。

2 时间步长与计算时间

在对水泵进行非定常数值模拟时，时间步长与计算时间存在着因人而异的情况^[15-16]。为了在满足计算要求的前提下，适当地增大时间步长和减少计算步数以提高工作效率，本文根据采样定理，对时间步长和计算步数应满足的基本要求进行了推导和分析。

设时间步长为$\Delta t$，则采样频率f_s为

${f_s} = 1/\Delta t$

(4)

由采样定理可知，在对泵内压力脉动进行采样时，采样频率f_s必须不小于泵内压力脉动最高频率f_c的两倍，即f_s≥2f_c。但是对于混流泵无从知道其最高频率，其内部压力脉动可能含有$0 \sim \infty $的频率成分，但每项研究所关注的最高频率f_H是有限的，因此可以将f_c取为f_H，即

${f_s} \ge 2{f_H}$

(5)

将式(4)代入式(5)得

$\Delta t \le \frac{1}{{2{f_H}}}$

(6)

由式(6)可知，研究所关注的频率f_H越高，则时间步长越小。

频率分辨率$\Delta f$是进行频率分析时的一个重要指标，$\Delta f$越小，对于复杂频率成分就能区分的越清楚。设水泵每分钟的转速为n，叶片数为z₁，导叶数为z₂，则泵轴转动频率(轴频)、叶轮叶片的旋转频率(叶轮叶频)和导叶叶片的相对旋转频率(导叶叶频)分别为：n/60，nz₁/60，nz₂/60。对泵内压力脉动进行分析时，在频域图上至少应区分开频率为零、轴频、叶轮叶频和导叶叶频4个谱峰，则频率分辨率$\Delta f$应满足

$\Delta f < \min \left( {\frac{n}{{60}},\frac{{n{z_1}}}{{60}} - \frac{n}{{60}},\frac{{n{z_2}}}{{60}} - \frac{{n{z_1}}}{{60}}} \right)$

(7)

然而，这只是频率分辨率应满足的最低要求，为了提高分析的精度，频率分辨率应越小越好。

由采样频率f_s和频率分辨率$\Delta f$可得，进行频域分析时所需的最少数据数量N=f_s/$\Delta f$。在进行频率分析时，通常只选取正半轴，故实际所需的最少数据量应为2N。因此，数值模计算时间T应满足:

$T \ge 2N \cdot \Delta t = 2{f_s}/\Delta f \cdot \Delta t = 2/\Delta f$

(8)

式中：计算时间T是指计算结果稳定后的计算时间，计算结果不稳定的时间不包括在内。由计算时间和时间步长，可得数值模拟的计算步数K为

$K = \frac{T}{{\Delta t}} \ge \frac{{4{f_H}}}{{\Delta f}}$

(9)

由式(9)可知，f_H越高$\Delta f$越小，则计算步数越多，所需的计算资源和计算时间也越多。

本研究主要关注水泵空化后中低频成分的变化，故将f_H取为轴频的60倍，即f_H=1 450 Hz，则时间步长$\Delta t$≤1/2 900 s=0.003 45 s。已知水泵的叶片数为4、导叶数为7，由式(7)可知，$\Delta f$<n/60=24.17 Hz，取$\Delta f$=13 Hz，则计算时间T≥0.154 s，计算步数K≥446，即在计算结果稳定后应至少再计算446步。

3 定常计算结果分析 3.1 混流泵的几何及网格

混流泵的几何如图 1所示，叶轮为4片，导叶为7片，泵进口直径为270 mm、出口直径为182 mm，叶轮出口直径为350 mm，设计转速为1 450 r/min。

叶轮和导叶采用Turbo-Grid划分网格，进出口直管用ICEM划分网格，叶轮选用J型拓扑，导叶选用H型拓扑，叶片表面采用O型网格，y⁺保持在50以下，经过网格无关性验证，最终确定计算域总网格数为390万，叶轮、导叶的网格如图 2所示。

3.2 湍流模型对非空化性能计算结果的影响

在研究过程中采用了SST和RNG两种湍流模型，分析了这两种湍流模型对混流泵非空化时水力性能和空化性能的数值模拟结果的影响。非空化时混流泵的功率与流量的关系如图 3所示。

从图 3中可以看出无论是SST湍流模型还是RNG湍流模型计算得出的功率值与实验值都存在一定误差，主要是因为实验测量结果包含了机械摩擦损失。SST湍流模型计算结果和RNG湍流模型计算结果的最大误差分别为-4.47%和-3.98%，但RNG湍流模型的计算结果与实验值更为接近。混流泵的扬程与流量的关系如图 4所示，从图中可以看出RNG湍流模型的计算结果与实验值更为接近，SST湍流模型和RNG湍流模型的最大误差分别为-4.2%和-2.1%。

3.3 湍流模型对空化性能计算结果的影响

空化是一种包含气液相间质量传输的非定常、可压缩的多相湍流流动现象，空化的数值计算涉及空化模型和湍流模型两个方面，选择合适的湍流模型会减小混流泵空化数值模拟的误差。

在计算混流泵的空化性能时，采用总压进口、流量出口，保持混流泵出口流量不变，通过不断降低进口的压力使泵内水体产生汽化。图 5为出口流量为380 kg/s时，混流泵的扬程随进口有效汽蚀余量NPSH的变化曲线。从图 5中可以看出，两种湍流模型计算的扬程-有效汽蚀余量曲线与试验结果趋势都一致，但RNG湍流模型的计算结果与实验值更接近。当NPSH减小到一定程度时，混流泵的扬程急剧下降，主要原因是空化产生的大量气泡堵塞了过流通道，使叶轮的做功能力降低。

临界汽蚀余量为扬程下降3%时的进口有效汽蚀余量，实验时混流泵的临界汽蚀余量为5.09 m，SST湍流模型和RNG湍流模型计算的临界汽蚀余量分别为4.37 m和4.42 m，绝对误差分别为0.72 m和0.67 m。图 6和图 7分别为RNG湍流模型和SST湍流模型的叶片吸力面空泡体积分数分布，从图中可以看出，随着NPSH的减小，空泡覆盖区域逐渐增大，说明泵内空化越来越剧烈；在NPSH相同时，两种湍流模型计算得出的空泡分布差别不大。

3.4 湍流粘度修正对空化计算结果的影响

RNG湍流模型和SST湍流模型都是通过引入湍流粘度，将湍流应力表示成湍流粘度的函数。由于在空化流中，存在着气相和液相两种组分，使得最初仅适用于单相流的湍流模型容易对湍流粘度过度预测，为了提高计算精度，按照Coutier-Delgosha等提出的思路对RNG湍流模型的湍流粘度进行修正^[17]：

${\mu _t} = f\left( {{\rho _m}} \right){C_\mu }\frac{{{k^2}}}{\varepsilon }$

(10)

$f\left( {{\rho _m}} \right) = {\rho _v} + {\left( {\frac{{{\rho _v} - {\rho _m}}}{{{\rho _v} - {\rho _l}}}} \right)^n}\left( {{\rho _l} - {\rho _v}} \right)$

(11)

式中：μ_t为修正后的湍流粘度，ρ_m为混合密度，C_μ为经验系数，k为湍动能，ε为湍动能耗散率，ρ_v为气相密度，ρ_l为液相密度，n取为10。湍流粘度修正前后的扬程-NPSH曲线如图 8所示。

在计算时，通过CEL语言将修正后的湍流粘度公式添加到CFX软件中。从图 8可以看出，湍流粘度修正后计算得出的扬程-NPSH曲线与实验曲线更加接近，湍流粘度修正前计算得出的临界汽蚀余量为4.42 m，修正后计算得出的临界汽蚀余量为4.75 m，绝对误差为0.34 m，比修正前的绝对误差减小了0.33 m。

4 非定常计算结果分析

为了分析空化时泵内的压力脉动的特征，对混流泵进行了非定常数值计算。采用Zwart空化模型和湍流粘度修正后的RNG湍流模型，以定常计算结果为初值，时间步长为0.000 345 s，计算时间为0.3 s，有效数据点数取480个。

分别对有效汽蚀余量NPSH=4.83 m的弱空化状态和NPSH=4.32 m的强空化状态进行了非定常计算，在进行非定常计算时，泵内监测点的分布如图 9所示，R₁、R₂、R₃位于叶轮进口，S₁、S₂、S₃位于导叶进口，P₁、P₂、P₃位于水泵出口。

4.1 压力脉动的时域分析

在NPSH=4.83 m时，各监测点在叶轮旋转一周内压力脉动的时域图如图 10所示，从图中可以看出，叶轮后面的监测点压力大于叶轮进口处监测点的压力，主要是由于叶轮做功，使得泵内压力升高。由于导叶的整流作用，使压能转化为动能，故泵出口处的压力小于叶轮和导叶之间的压力。

在导叶进口处，从轮毂到轮缘的压力逐渐增大，靠近轮缘处的压力最大(如S₁所示)，但在喷口处却是靠近轮缘处的压力最小(如P₁点所示)。由于导叶进口处监测点的压力直接受到叶轮做功的影响，而轮缘处的半径最大、做功能力最强，故压力也最大；而泵出口距离叶轮较远，主要受导叶整流作用的影响，在泵出口断面，由于靠近轮缘处的水流速度大，故压力低。

此外，所有监测点的压力在叶轮旋转一个周期内均呈现出4个小周期，这反应了叶轮旋转对泵内压力的影响，各监测点压力脉动的周期主要跟叶片数(4片)相关。

图 10反映了各监测点的压力值及波动情况，但是各监测点之间压力脉动的相对波动情况难以直观看出。为此，采用归一化压力系数来分析不同位置处压力的相对波动情况，归一化压力系数为：k=P-P-/P-，其中P-为一个周期压力脉动的平均值。选取靠近轮缘的3个监测点(R₁、S₁、P₁)的压力进行分析，3个监测点压力的相对波动情况如图 11所示。

由图 11可以看出，压力波动幅度从大到小依次为R₁、S₁、P₁，即叶轮进口处压力波动最剧烈，泵出口处压力波动最平稳。由于水泵进口前的来流是不旋转的，当不旋转的水流遇到旋转的叶片，就会产生很大的压力波动。泵出口的水流由于受到导叶整流的作用，所以压力波动相对平稳。

4.2 压力脉动的频域分析

为了分析空化时泵内压力脉动的频率情况，对R₁、S₁两点的压力脉动进行了傅里叶变换，分别绘出了两点在弱空化状态(NPSH=4.83 m)和强空化状态(NPSH=4.32 m)时的频域图，如图 12所示。

图 12中Fn表示轴频。从图中可以看出，在两种空化状态时，R₁和S₁两点都以叶轮叶频及其倍频为主。此外，在NPSH=4.32 m时两个监测点的主频幅值均大于NPSH=4.83 m时的主频幅值，主要是因为随着泵进口有效汽蚀余量下降，泵内空化越来越剧烈，空泡大量溃灭产生的压力所导致。

5 结论

运用Zwart空化模型对混流泵的空化性能进行了定常和非定常数值模拟，并对混流泵内部监测点的压力脉动进行了分析，得出以下结论：

1) 对混流泵的非空化和空化性能进行计算时，RNGk-ε湍流模的计算结果都比SST湍流模型的计算结果更接近实验值。通过对湍流粘度进行修正，使数值模拟的计算误差进一步减小。

2) 导叶进口处的平均压力最高，叶轮进口处的平均压力最小；叶轮进口处的压力波动最剧烈，泵出口处的压力波动最平稳；在叶轮出口，压力从轮毂到轮缘逐渐增大；在泵出口，压力从轮毂到轮缘逐渐减小。随着进口有效汽蚀余量降低，泵内压力脉动主频的幅值增加。