非线性空泡模型在轴流泵空化模拟中的评估分析

2020-07-14张天宇苏华山袁建平

流体机械 2020年6期

洪锋，张天宇，苏华山，袁建平，张帆 2，

（1.三峡大学机械与动力学院，湖北宜昌 443002；2.三峡大学水电机械设备设计与维护湖北省重点实验室，湖北宜昌 443002；3.江苏大学国家水泵及系统工程技术研究中心，江苏镇江 212013）

0 引言

空化是轴流泵内部不可避免的一种水动力学现象，空化严重时会破坏过流部件表面和产生振动噪声，严重影响机组的高效安全运行。近年来，基于Navier-Stokes方程和相间质量传输模型框架的空化数值计算方法，由于能描述空化的动量传输和质量传输两个过程，成为空化流动研究手段的主要趋势。Athavale等［1］应用Singhal模型对3种不同比转速船用轴流泵空化流场进行了数值模拟，发现不可凝结气体（non-condensable gas，NCG）能够降低泵的扬程及影响空化的发生和发展，还指出了该空化模型中蒸发与凝结系数限制了其应用。此外，Zwart、Schnerr-Sauer空化模型在模拟单个水翼、螺旋桨与叶片泵叶片表面上的空穴形态均被报道存在不足之处［2-5］。上述几种空化模型的相变率表达式均基于单一球形空泡动力学中常用的Rayleigh-Plesset（R-P）方程的线性形式，该方程作了过多假设，如忽略空泡二阶运动项、液相表面张力等，使得该类空化模型还不能完全且较好地胜任轴流泵叶轮内空化流动问题的研究。

为了改善这种基于线性R-P方程的空化模型模拟空化流动的能力，本文提出了一种非线性空泡模型，联立两相-3种组分的均相流模型，得到了一种非线性空化模型。运用用户自定义语言，把该模型嵌入到ANSYS FLUENT 14.5平台，并对某一型号的轴流泵叶轮内空化流动进行数值计算。通过与传统的Schnerr-Sauer模型的数值计算结果及实验结果作对比，验证非线性空化模型在预测轴流泵空化流动中的可行性与准确性。

1 数值计算模型

1.1 非线性空泡模型

非线性空泡模型基于广义的球形空泡动力学R-P方程建立。液态水及空泡均属于低黏度工质，其运动黏度系数均为 10-6数量级，Ivany等［6］和Knapp等［7］研究单一空泡的增长与溃灭过程发现，水的粘性引起的阻尼效应可以忽略不计。同时，Schnerr等［8］在研究绕三维 Twist-11扭曲水翼空化流场时发现，云空泡脱落、溃灭及随后的冲击脉动动力学特征均是由空泡的惯性效应决定，黏性的阻尼作用十分微小。此外，本文作者在前期工作［9］中还证明了表面张力系数对空泡半径增长的影响明显高于粘性力的影响。因此，非线性空泡模型在线性R-P方程基础上还考虑了空泡半径二阶导数及液相表面张力的影响，线性、非线性R-P方程分别为：

式中pv（t）——空泡内压力；

p∞（t）——无穷远处压力；

t——时间；

ρl——液相密度；

R——空泡半径；

S——表面张力系数。

式（2）中，等号左边为压力驱动项，等号右边依次为二阶导数项、一阶导数项及表面张力项。

为了对比2种不同空泡动力学模型对空泡半径变化的差异，利用MATLAB软件，采用4-5阶RUNGE-KUTTA算法分别对式（1）和（2）进行求解，采用的初始条件如下：空泡初始半径R0=4.0 µm、液相表面张力系数S=0.073 N/m、驱动项pv（t）-p∞（t）=105Pa。计算可以得到该空化核流过低压区时R（t）一阶导数的一个典型解如图1所示，从图中可以看到，线性R-P方程预测的空泡半径增长率为一恒定值，且当空泡半径小于30 µm时，非线性R-P方程预测的空泡半径增长率大于线性R-P方程计算的结果。由此可见，线性R-P方程在空泡初生阶段（R＜30 µm）会过小预测空泡半径的增长速度。

图1 根据线性、非线性R-P方程计算得到的空泡半径增长率

基于两相-3种组分的均相流模型及非线性空泡模型，建立了一种非线性R-P空化模型，该空化模型考虑了不可凝结气体、湍流脉动压力对空化初生的影响。相变率推导详细过程，以及模型系数无关性分析参见文献［9］，相变率基本方程组为：

式中Se——蒸发率；

N0——空化核密度，m-3，取N0=1×1013m-3［8］；

ρg——NCG密度，按照密度形式理想气体状态方程计算；

R0——空泡初始半径，µm，R0=1 µm［10］；

Sc——凝结率；

fg——NCG质量分数，取1.5×10-6［10］；

psat——饱和蒸汽压，Pa，取psat=3 169 Pa；

k——湍动能。

1.2 湍流模型

传统的双方程湍流模型会过度预测闭合空穴尾部的湍流粘度，使得空穴形态预测失真，而在大型复杂三维流动问题中应用直接数值模拟（DNS）或者大涡模拟（LES）时很难找到一个网格无关解，且消耗大量的计算资源。因此，为了保证计算精度及节省计算成本，本文采用一种修正的滤波器湍流模型（MFBM）［11］。该模型通过一个函数fMFBM联立原始滤波器模型（FBM）和密度修正模型（DCM）来共同修正传统双方程湍流模型模拟空化两相流时的不足，修正后的模型兼备FBM和DCM的优点。MFBM的湍流黏度系数为：

式中，n=10、Cμ=0.09、C3=1.0，λ为滤波尺。

2 计算域及网格划分

研究对象为一比转速ns=692的轴流泵，设计流量为Q=312 m3/h、设计扬程H=2.95 m、转速n=1 450 r/min，其过流部件主要包括进出口管路、叶轮和导叶，三维造型过程中保留了叶顶间隙处水体，其厚度为0.5 mm。过流部件如图2所示，整个计算域网格划分在ICEM CFD软件完成，采用结构化网格。

图2 模型轴流泵计算域

为了尽量消除网格数量造成的求解误差，分别对计算域划分了6种不同数量的网格作网格独立性分析，网格节点数分别为N1=1 170 764、N2=1 756 146、N3=2 341 528、N4=2 926 910、N5=3 212 292和N6=3 697 674。在最优工况下（Q=302 m3/h）［12］，计算基于不同网格模型泵的扬程来确定最终网格，计算结果如图3所示，从图中可以看到，当网格节点数大于N5时，扬程的相对变化在1%以内，表明继续加密网格对数值计算结果的影响可以忽略。N5网格计算得到的扬程为3.20 m，对应的试验扬程为3.04 m，二者相对误差为5.26%，计算精度较高。因此，采用网格节点数为N5的网格作为本文最终的计算网格，其过流部件网格总成及叶轮网格示意如图4所示。数值计算时，以非空化流场的计算结果作为空化非定常数值计算的初始值，在非定常数值计算中设置时间步长为Δt=0.11 ms。采用托马（空化）系数来判断空化程度，表达式为：

式中NPSHA——装置空化余量；

H——扬程。

图3 不同数量网格下计算扬程变化

图4 网格示意

3 数值方法验证

该模型轴流泵空化特性测试在江苏大学国家水泵及系统工程技术研究中心完成，按照回转动力泵水力性能验收试验标准GB/T3216-2005进行模型泵性能测试，工作介质为常温清水，试验装置简图如图5所示。

图5 轴流泵空化测试系统示意

图6示出0.8QBEP，1.0QBEP和1.2QBEP3种工况下不同空化模型计算得到的空化性能曲线和临界空化余量（扬程下降3%），与相应试验测量值的对比，其中计算扬程均为时均化后的结果，计算总时间为0.33 s。

图6 空化特性预测值与试验测量值的对比

从图中可以看到，不同工况下，当空化数足够大时，泵内没有发生空化，不同空化模型预测到的扬程十分接近。随着空化数的减小，当空化发生时，各工况下扬程的测量值相比于计算值提前发生断裂，说明数值计算会一定程度上低估泵内空化的初生及发展过程，但非线性模型预测的扬程断裂时的空化数明显大于Schnerr-Sauer模型的计算结果，表明非线性模型会更早预测到空化的发生，造成这种现象的原因有两点：（1）非线性模型中考虑了NCG的影响，该部分气体充当“空化核”的作用，加速了空化的发生；（2）非线性模型中考虑了湍流脉动对空化临界压力的影响。从图6（d）中可以看到，非线性模型预测到的临界空化余量与试验测量值吻合度更高。为了进一步对比两种空化模型捕捉空化特征的能力，截取了空化充分发展阶段（扬程下降10%）不同空化模型预测的叶轮流道空泡形态，邻近比转速（ns=728）相同空化数条件下某一轴流泵的可视化试验结果作为数值计算结果的对比如图7所示。可以发现，不同空化模型的预测结果均捕捉到了尾缘处的小体积空泡团向叶轮流道脱落，此时实验结果中叶顶间隙内附着空穴与延伸至叶轮流道内的空泡混合形成的总空泡团外形上呈三角形，非线性模型预测到的叶顶间隙附近处空泡团的三角形外形明显与试验结果吻合更好，说明非线性模型预测的叶顶间隙内空化发展更加充分。综上说明，非线性空化模型相比于传统的Schnerr-Sauer模型预测轴流泵空化性能及对空化细节捕捉的能力更强。

图7 空化形态与试验结果的对比

4 叶轮内2种主要空化特征

4.1 叶顶间隙空化

轴流泵叶轮内空化主要包含叶顶间隙空化和叶片吸力面空化两大类，张德胜等［13］参考轴流式风机叶轮内流动特点，将轴流泵叶顶间隙处流动进行了归类和描述，认为叶顶间隙空化包括间隙内附着空化、泄漏涡涡带空化、剪切层分离涡空化及尾缘脱离涡空化。鉴于非线性模型对空化捕捉能力优于传统Schnerr-Sauer模型，截取了1.0QBEP，空化数σ=0.59工况下非线性模型计算得到的叶轮内相对速度流线分布，以及叶顶间隙附近处空泡分布分别如图8，9所示。从图中可以看到，非线性空化模型基本上捕捉到了叶顶间隙附近上述4种不同类型的空化。

图8 叶轮内相对速度分布（1.0 QBEP，σ=0.59）

图9 叶顶间隙附近空化分类

图10示出1.0QBEP工况不同空化数下，非线性模型计算得到的叶顶间隙附近处的空泡形态分布，为了使计算结果精确度更高，该图中的所有的空泡分布均采用非定常求解器计算后处理得到，每一工况均运行了2 000迭代步，即每隔t=0.22 s的时间处理一张分图。可以看出，当空化数降低至σ=2.09时，叶顶间隙内初生空化形成；当空化数大幅降低至1.41时，间隙内附着型片状空穴长度出现了较为明显的增长，如图10（b）所示。随着空化数缓慢降低至σ=1.26时，间隙内附着型片状空穴的长度出现了显著的增长；随后，当空化数降低至σ=0.94时，发现附着空穴的长度出现了迅速增长。由此可见，叶顶间隙附着型空化生长速度先极为缓慢，当空化数降低至一定程度时出现快速增长。当空化数低至σ=0.59时，此时叶顶间隙附近的空化由叶顶泄漏涡涡带空化、剪切层分离涡空化以及叶顶翼型脱离涡空化构成。同时还说明，叶顶泄漏涡涡带空化、剪切层分离涡空化及叶顶翼型尾缘脱离涡空化的发生明显滞后于叶顶间隙内的附着型空化。

图10 非线性模型计算得到的不同工况下叶顶间隙附近空穴形态（Δ t=0.22 s，Q=1.0 QBEP）

图11示出1.0QBEP工况不同空化数下，非线性模型计算得到的50%弦长叶顶间隙内的空泡体积分数及速度矢量分布。从图中可以看到，当空化数σ=1.05时，在叶片顶端，由于叶片和泵腔内壁面存在间隙，流体的粘性力不足以克服叶顶翼型工作面和吸力面之间的压差，叶片压力面后的流体收缩被迫加速流入叶顶间隙，当流体通过叶顶间隙区域后以射流的形式逆向流入主流，并与叶片吸力面叶顶附近区域的主流相互作用形成逆时针的泄漏涡，此时因泄漏涡涡心处的压力值高于该温度下的饱和蒸汽压而没有发生泄漏涡涡带空化。当空化数降低至σ=0.91时，叶顶处吸力面发生了空化。随着空化数降低至σ=0.85时，在射流区域内形成了射流剪切层空化，而且此时吸力面上附着空化厚度相比于前一工况出现了急剧增加。当空化数降低至σ=0.78时，泄漏涡涡心处发生空化，泄漏涡涡带空化范围随着空化数的降低而不断扩大，其内部所包含的空泡的体积分数也逐渐增大。在主流的作用下泄漏涡涡带逐渐向吸力面靠近并与叶顶处吸力面空化发生掺混。叶顶泄漏涡是一种高度的三维流动结构，叶顶泄漏涡与涡带空化互相作用。随着空化数的降低，涡带空化呈现出不同的形状，其轴向位置逐渐向吸力面靠近，泄漏涡涡心处的最低压力大小决定了涡带空化区空泡相的体积分。

图11 非线性模型计算得到的不同工况50%弦长处叶顶间隙内流动（Δt=0.22 s，Q=1.0 QBEP）

4.2 叶片吸力面空化

图12示出1.0QBEP工况不同空化数下，非线性模型计算得到的叶片吸力面上空泡体积分数分布。从图中可以看到，当空化数为σ=1.47时，初生空化出现在叶片吸力面进口靠近轮缘处。对比图10（a）所对应的空化数可知，轴流泵叶片吸力面空化明显滞后于叶顶间隙内空化。随着空化数的降低，空化区逐渐向叶片出口及轮毂方向扩展，靠近轮缘处区域内空泡相体积分数急剧增加。当空化数降低至σ=0.85时，吸力面空化已经发展至轮毂，并在2/3倍弦长处形成了一条与叶轮轴向垂直的空泡带状区域，此时轮毂表面区域会因发生空化而附着有空泡团，即毂涡空化。在临界空化工况点时，吸力面上空化带状区分别向叶片进出口方向延伸，并逐渐与轮缘处的空化区掺混一起，靠近轮毂处空化区空泡相体积分数相比前一工况急剧增加。随着空化数继续降低，毂涡空化程度加剧，掺混后的吸力面空穴末端逐渐向叶片出口延伸，如图12（f）所示。当空化数降低至σ=0.63时，吸力面上2/3倍弦长处至叶片出口，在整个叶高范围内完全被高体积分数的空泡覆盖。对比图12（g）和（h）可以发现，当空化数降低至一定程度时，叶片出口靠近轮毂处发生了严重空化，空化区内空泡体积分数极高，在整个吸力面空化区域的后部靠近轮毂的位置，空泡体积分数发生了急剧变化，这种剧烈变化是由于空泡运动至叶片出口附近的高压区发生溃灭所造成。可以推断，当轴流泵叶轮内发生严重空化时叶片出口位置且靠近轮毂处最容易发生空蚀破坏。