基于修正PANS模型的轴流泵空化特性数值模拟

2023-01-02俞芸芸周大庆于安刘佳佳

排灌机械工程学报 2022年12期

俞芸芸,周大庆,于安,刘佳佳

(1. 南京航空航天大学金城学院，江苏南京 211100；2. 河海大学能源与电气学院，江苏南京 211100; 3. 南水北调东线江苏水源有限责任公司扬州分公司, 江苏扬州 225000)

轴流泵是一种低扬程、大流量且应用广泛的泵型，近年来其受到越来越多的关注.空化是评价泵性能的三大指标之一，与效率、稳定性的研究相比，空化的研究还不是十分的充分，空化问题制约着现代水泵行业的发展，也是国内外学者研究的热点之一[1].轴流泵空化特性的研究一般采用数值模拟和试验验证等方法.随着计算机科学的快速发展，CFD技术被广泛应用于数值模拟中，国内外学者对算法的改进也使得数值模拟的精确度越来越高，数值模拟已经成为研究空化现象的重要手段之一[2].

近年来，用于空化计算的湍流模型得到了进一步修正，其中一种由标准k-ε模型改进而来的局部时均化模型(partially-averaged Navier-Stokes model, PANS模型)被越来越多地应用于空化流动计算中[3-4].GIRIMAJI等[5]认为PANS模型是雷诺平均RANS和直接数值模拟DNS之间的桥梁模型，并提出其控制参数是未分解湍动能比值fk以及未分解耗散率比值fε，并对该模型进行了初步验证.LAKSHMIPATHY等[6]利用PANS模型对圆柱绕流进行了模拟，并采用不同的参数fk和fε值研究了PANS模型的可实施性和准确性.罗大海等[7]采用PANS湍流模型对超声速斜面空腔流动进行了数值模拟，并取定值fε=1，而fk值随流场中的时间和空间变化，结果表明，相比于取定值，采用动态的fk可以更加准确地模拟出三维流动的结构.黄彪等[8]采用基于k-ε模型改进的PANS湍流模型进行水翼的空化流动计算，证明了fk值对空化计算的精度有较大影响.施卫东等[9]采用PANS模型对绕二维Clark-Y型水翼进行了非定常计算，通过对fk取不同的值研究了其对非定常空化的影响.石磊等[10]采用PANS模型对轴流泵叶顶间隙空化进行了数值模拟和试验验证，准确地模拟出了叶顶区的空化特性.HUANG等[11]采用修正的PANS模型对绕台阶流模型进行数值计算，证明了其在流场压力、速度分布和旋涡形态上具有更准确的预测能力.

基于前人研究，文中利用计算软件CFX二次开发技术对标准PANS模型中的fk值(未分解湍动能占总湍动能的比值)进行动态定义，形成修正的PANS模型.利用该湍流模型对轴流泵全流道进行空化数值模拟和试验测试，验证修正PANS模型在轴流泵空化特性模拟上的可靠性.

1 数值计算理论及方法

1.1 基本控制方程

CFD计算方法中涉及的流体动力学控制方程是基于质量守恒定律、动量守恒定律、能量守恒定律形成的3个控制方程，根据常规水力机械运行条件，流动介质为水，且常温常压，不考虑热量交换，一般认为其内部流场为等温不可压缩的三维流动，即不涉及能量交换，只包含质量和动量守恒定律[12]，因此，基本方程包括连续性方程和Navier-Stockes方程.

连续性方程[13]

(1)

Navier-Stockes方程[14]

(2)

式中：ρ为密度；p为流体压力；ui，uj为不同方向流速的瞬时值；μ为动力黏性系数；μt为湍流黏性系数.

1.2 湍流模型修正

文中采用的湍流模型是修正的PANS模型，即基于标准PANS模型，利用CFX的二次开发技术编写CCL语言，对PANS模型中的fk值进行动态定义.

由于轴流泵内湍流流动的高度复杂性，在数值模拟计算中需要对湍流模拟方法进行研究和分析，根据不同的情况选择不同的湍流模型进行数值模拟，以取得与实际流态更为吻合的结果.数值模拟方法分为直接模拟(DNS)和非直接模拟，非直接模拟应用最为广泛的是雷诺时均法(RANS)，雷诺时均法的核心是如何准确合理地求解时均化的雷诺方程，从而取代直接求解瞬时Navier-Stokes方程的方法[15].常见的雷诺时均法包括零方程模型、一方程模型、k-ε模型、k-ω模型等.文中采用的修正PANS模型是基于k-ε模型改变而来，因此对该模型进行详细介绍.

k-ε模型在计算时不直接计算雷诺应力，而是通过湍动黏度(涡黏系数)对湍流应力进行参数表达，即根据Boussinesq涡黏假定得出[16]

(3)

式中：δij的取值方式中，当i=j时，δij=1；当i≠j时，δij=0；k为湍动能，即

(4)

由式(3)可知，k-ε模型的关键在于求解μt.在式(4)关于湍动能k的方程基础上，引入湍动能耗散率ε，形成k-ε模型，ε计算式为

(5)

可得

(6)

式中：Cμ为经验系数，根据LAUNDER的研究，Cμ取值为0.09[17].

k-ε模型比零方程和一方程模型在计算精度和计算成本上有了很大改进，且具有广泛的适用性，但是k-ε模型在推导过程中仍进行了一定的假设和近似，比如它是针对湍流发展充分的高湍流区域建立的，在近壁区等湍流发展不充分的地方可能存在失真，因此需要采用壁面函数法等措施进行修正；此外，k-ε模型在强旋流、壁面弯曲流动时计算存在一定误差.

文中基于k-ε模型，利用局部时均法引入fk和fε对运输方程进行修改，形成PANS模型.fk，fε定义如下

(7)

式中：k，ε分别为总湍动能和总耗散率；下标u代表PANS模型中未分解的物理量.

PANS模型将标准k-ε模型的运输方程改为

(8)

(9)

并针对式(9)中的以下系数进行了修正

(10)

式中：Cε1=1.44；Cε2=1.92；σk=1.0；σε=1.3.

PANS模型是可以从直接模拟模型(DNS)过渡到雷诺时均化模型(RANS)的计算模型，这个过程通过调整模型中的fk和fε进行控制，fk和fε的值取[0,1].根据大量文献，fε在高雷诺数时取值为1更为合适[18].而fk的取值对计算的结果影响较大，当fk=1时，PANS控制方程将复原到RANS模型；当fk=0时，PANS控制方程直接对N-S方程求解，即DNS模型，在fk值由小变大的过程中，PANS模型从DNS平滑过渡到RANS模型.当fk取值较小时，虽然可以使得求解结果更为精确，但对应的网格数量增多，且对计算机能力要求也大大增加，因此在PANS模型计算中fk的取值尤为关键.

目前应用较为广泛的取值方法是将fk值取为[0,1]的一个数或几个数，然后进行计算,从而选出较为适合的值，此方法是在整个计算域取统一的值，需要较大的计算量进行对比且无法得到不同区域内最合适的fk值.

文中利用CFX软件二次开发技术中的CCL语言，对fk值进行了动态定义，即空间上瞬时地根据当地网格条件和湍流长度尺度定义fk的值，形成修正的PANS(modified partially-averaged Navier-Stokes， MPANS)模型，fk公式为[19]

fk=min[1,3(Δ/l)2/3],

(11)

式中：l为当地湍流长度尺度，l=k1.5/ε；Δ为当地网格尺度，Δ=(ΔxΔyΔz)1/3.

由式(11)可知，fk的取值可以根据当地网格和湍流尺度进行动态修改，当湍流尺度较小时，fk=1，此时计算方法改为RANS法；而当湍流尺度较大时，fk<1，此时计算方法趋向于直接计算湍流方程[20]，可更为准确地模拟出该区域的流动特性.

修正的PANS模型解决了标准PANS模型无法精确取到合适的fk值，以及无法在不同区域取不同fk值的问题；与普通雷诺时均化模型相比，修正PANS模型在高湍流流动区域内数值模拟的精度更高；与直接模拟相比，在湍流尺度小的区域进行了时均化处理，减少了计算成本，提高了网格尺寸的适应性.因此，在湍流计算中选取修正的PANS模型更为合适.

2 数值模拟

2.1 计算域及网格划分

文中研究对象为立式轴流泵装置，根据工程要求对全流道进行三维建模，模型由进水流道、叶轮、后置导叶体、出水流道构成，模型额定转速为1 367 r/min，设计流量为0.375 m3/s，叶轮直径为300 mm，轮毂直径为120 mm，叶轮叶片数为4，叶顶间隙为0.3 mm.利用ICEM进行六面体结构化网格划分，为了准确反映边界层内的流动情况，精确模拟靠近叶片表面的空化形态，对网格的Y+值进行了控制，通过调整固体壁面第1层网格的高度来改变Y+值[21]，将其控制在30以内，流体域整体网格数量约为300万.

2.2 求解方法设置

利用CFX软件进行轴流泵全流道空化定常计算，将修正的PANS模型导入湍流控制方程中，计算离散方法采用有限体积法[22]，湍流模型对流项采用高阶精度、中心差分格式，收敛控制设置为最少迭代步数不少于5 000步，残差收敛精度设置为10-5.空化模型采用Zwar-Gerber-Belamri模型[23]，气液两相之间的传输采用基于Rayleigh-Plesset方程[24-25]的计算方法.进口边界条件采用总压进口，以控制在不同汽蚀余量下的空化计算，进口处气相体积分数为0，液相体积分数为1，出口设置为质量流量出口.流体域介质设定为“Water at 25 ℃”，“Water Vapour at 25 ℃”，水在25 ℃时对应的饱和蒸汽压力为3 175 Pa.

2.3 计算结果分析

根据基于修正PANS模型的定常计算结果，对泵的空化性能进行综合分析.为了准确对比不同工况下叶轮内同一位置处的性能，在叶轮内进行分层设置，在叶轮内定义某个径向截面，引入量纲为一的变量Span,

Span=r/R，

(12)

式中：r为该截面距离轮毂的半径长；R为叶轮室距离轮毂的半径长.Span=0时，该截面为轮毂面；Span=1时，该截面为叶轮室壁面.

2.3.1 不同汽蚀余量下的空化性能

选取不同汽蚀余量即NPSH为6.81，5.18和4.52 m时的计算结果，进行空泡体积分数δ、涡量Ω分布、压力p分布的对比，在叶轮内取截面Span=0.95，Span=0.50并进行平铺展开，对比结果见图1，2.

图1 NPSH=6.81，5.18，4.52 m时Span=0.95处空化特性对比

图2 NPSH=6.81, 5.18, 4.52 m时Span=0.50处空化特性对比

根据图1，2的对比可知，从叶轮轮缘至叶轮中部的径向分布上不同汽蚀余量下表现的规律一致，随着汽蚀余量的减小，各截面上空泡体积分数变大，速度涡量变大，压力逐渐减小，而这些变化在叶片的吸力面表现得更为明显，这与空泡的主要发生位置相对应.对比截面Span=0.95和Span=0.50处的结果可知，叶片上的空泡在靠近轮缘处相较于叶片中部更多，叶片轮缘处空化相对较严重，且轮缘处的空泡主要分布在叶片的前端，而叶轮中部的空泡则在叶片的尾部较多；叶轮中部速度涡量值相比于轮缘处较小，且整体分布规则，整体大涡量区域减小，叶片表面压力由中部向轮缘处逐渐增大，叶片中部由于空化较弱，压力相差值减小，压力分布较轮缘处均匀.通过以上的分析可知，叶片表面压力和速度的变化与空化发生的位置和区域大小息息相关，不同汽蚀余量下各变量也表现出不同的规律.

2.3.2 不同流量下的空化性能

在不同流量工况下进行了数值计算，分别选取流量为0.87Q，0.93Q，1.00Q，1.07Q，1.13Q，计算得出对应的临界汽蚀余量分别为4.97，5.16，5.37，5.72 和6.40 m.随着流量的增大，扬程降低，临界汽蚀余量增大，说明泵是否发生空化不仅与进口压力有关，也与流量大小有关，在某一固定的汽蚀余量下，如5.37 m时，泵在设计流量下运行时处于临界空化状态，而在大流量工况下已发生严重空化，小流量工况下则尚未发生空化.

分析流量对空泡的影响，选取相同汽蚀余量，对比分析流量对轴流泵空化性能的影响.当汽蚀余量NPSH=4.87 m时，对不同流量工况下泵叶轮内的空化状态进行分析，由图3各流量下的空泡分布可知，空泡体积分数为0.15.

图3 不同流量下的空泡分布图

由图3可知，在相同汽蚀余量下，不同流量工况表现出不同的空泡分布特点.随着流量的增大，空泡分布的整体量在变大，空化逐渐严重，叶片吸力面的空泡分布由进水边逐渐向出水边转移，而叶片压力面进水边区域的空泡出现了从无到有的变化，当流量达到1.13Q时，叶片吸力面空泡与压力面空泡的产生和溃灭连为一体，使空泡不只是附着于叶片，而是充斥在叶间流道内，此时空化十分严重.

分析相同汽蚀余量(4.87 m)下不同流量工况时的空化对泵内流态的影响，取叶轮和导叶中间截面(Span=0.50)并对其展开，各流量对应的流线分布如图4所示.

图4 NPSH=4.87 m时不同流量下的流线分布

由图4可知，在不同流量工况下，从流线的整体分布来看，叶轮内部的流线变化不大，导叶内流态变化较大.在小流量工况0.87Q下，导叶吸力面出现较大区域的旋涡，随着流量的增加，旋涡逐渐减小直至消失；在大流量工况1.13Q下，导叶吸力面流线较为平稳，但是在导叶的进口处出现扭曲.随着空化逐渐严重，导叶体内的流态逐渐平稳，说明空化未对导叶体内流态产生大的影响，流量的增加使得导叶流道内的水流具有一定的夹持作用，能够减小水流的旋涡.

为了分析不同流量下空化对叶轮内流态稳定性的影响，得到图5所示叶轮内不同截面上平均湍动能值曲线.由图可知，在相同汽蚀余量(4.87 m)下，不同流量工况时叶轮内部的湍动能分布不同，随着流量的增大，湍动能值变大，湍流耗散变严重，这与空化逐渐严重有关，空化的发生导致叶轮内部流动的不稳定性加剧.从湍动能在叶轮的径向分布特点来看，各流量工况下湍动能值具有较为明显的特点，其值整体上随着流量的增加而增大.对于同一流量工况下，湍动能值在Span=0.05即靠近轮毂处最大，沿径向先减小至最小值，而后在Span=0.95即靠近轮缘处增大，说明在轮毂和轮缘处湍流耗散较严重，水流稳定性降低，与空泡的边缘处重合.大流量下空泡体积分数较大，其溃灭对水流湍动能的影响更大，这导致了流量1.13Q下叶轮出口和叶轮内部的湍动能值均比其他工况下大.综上可知，轴流泵在不同流量下发生空化所表现的性能有所不同，空泡的分布、流态的变化以及湍动能值的大小均随着流量的改变呈现出规律性变化.

图5 不同流量下叶轮内各截面平均湍动能值曲线

3 试验验证

在对轴流泵进行空化定常计算并得出轴流泵空化特性以后，对轴流泵模型进行了试验，并对比数值计算和试验结果中的空化特性曲线和空泡分布，验证文中所采取的计算方法的可靠性.

本次空化试验在河海大学水力机械多功能试验台上进行，试验中为了验证数值计算的准确性，对轴流泵的外特性和空化特性均进行了测试，主要涉及的参数有以下几个：扬程、流量、转矩、真空度、压力等.空化试验采用能量法，保持设计流量为0.375 m3/s，通过系统回路内抽真空减小试验系统的空化余量值，从而减小有效汽蚀余量.在抽取不同真空度值时测试轴流泵对应的扬程，并绘制扬程和汽蚀余量的曲线，根据扬程下降2%确定临界汽蚀余量.绘制试验汽蚀余量-扬程曲线，并与数值模拟计算结果进行对比，如图6所示.

图6 试验与数模空化曲线对比

由图6可知，随着汽蚀余量的减小，泵的扬程逐渐减小，空化越来越严重，试验与数值模拟结果曲线趋势相同，但试验测得的扬程较数值计算结果略大，两者误差小于2.9%，这是由试验条件及试验系统引起的误差，且在合理范围内.图中点a为试验所测得的临界汽蚀余量，其值为5.68 m；点A为数值模拟计算所得临界汽蚀余量，其值为5.37 m，试验测得的临界汽蚀余量较数值计算结果偏大，空化性能较数值结果略差.这是由于实际液体中包含某些不可溶气体，形成了空化核，空化核的存在将会降低液体的抗拉强度，试验时若对水体杂质的排除不够彻底，将会增加水体中空化核的含量，而数值模拟时默认的空化核数量较少，这将导致试验空化性能比数值计算结果差.另外，试验时叶轮体各壁面具有一定的粗糙度，当泵内水流流经某些粗糙度较大的地方时，将会引起速度和压力的突变，增加附近区域的湍流强度，导致在粗糙区域内产生局部空化，也在一定程度上降低了泵的空化性能.

试验中除对效率、扬程等数据进行采集外，还在叶轮室壁面上安装了透视窗口，利用闪频仪对叶轮内的空泡进行拍摄，拍摄结果如图7所示.

选取试验和数值计算中空化开始阶段和空化较为严重的2个工况，即图6中的点B，b和点C，c，分别进行空泡分布对比，数值计算结果按照空泡体积分数0.15以上显示，如图7所示.