敏 政, 朱月龙, 黎义斌, 张 梅
(1.兰州理工大学能源与动力工程学院, 甘肃 兰州 730050; 2.甘肃省流体机械及系统重点实验室, 甘肃 兰州 730050)
目前,常用的湍流数值模拟方法有直接数值模拟(DNS)、Reynolds平均法(RANS)及大涡模拟(LES)。直接数值模拟计算资源消耗大,计算周期长,不容易实现。RANS方法在计算附体和小分离流动时可获得令人满意的结果,然而不能准确模拟大范围的分离流动。RANS方法的缺陷促使人们发展了可直接求解大尺度运动大涡模拟方法。LES方法可精确模拟分离流动以及与几何相关的大尺度非定常运动,所花费的计算资源仅需DNS方法的很小部分,但当模拟边界层流动时,LES所需的计算资源与DNS几乎相当。另外,LES方法的近壁模式尚不成熟,不能完全分辨出高Reynolds边界层的近壁流动结构,所描述边界层的增长和分离不准确[1-2]。
1997年, Spalart等[3]提出分离涡方法(Detached Eddy Simulation, DES),基于局部网格间距和湍流长度尺寸,在RANS和LES模型间转换。DES湍流模型在边界层运行RANS,大分离流动区域转换成LES。2006年,Spalart等[4]又提出了延迟分离涡方法(Delayed DES, DDES),引入一个延迟函数,重新构造了DDES的长度尺度,很大程度避免了MSD(Modeled stress Depletion)问题。
刘若阳等[5]采用DDES研究叶栅分离旋涡的非定常流动特性,发现攻角和叶栅几何参数的改变对于流场旋涡运动有很大影响,叶背分离点的相对位置则是影响旋涡非定常运动形式和漩涡脱落频率的主要内在因素。林敦等[6]通过DDES方法分别研究了绝热壁面条件和等温壁面条件下跨音高压透平导叶LS89中的基本流动现象,证明了DDES方法计算可靠性以及对于非定常精细结构的出色分辨力。在大迎角静态翼型大分离流动模拟中,DDES方法捕获了非定常RANS计算未能获得的机翼背风面的涡脱落现象[7]。在二维扩压流场中存在尾缘脱落涡和吸力面脱落涡两种脱落涡,它们的脱落频率随着冲角和马赫数的变化而变化[8]。压力脉动主要由叶片和蜗壳的动静相干作用,蜗壳内的压力脉动比较明显[9]。进行全局和局部流场诊断,通过和用压力、速度等传统分析方法的对比,结合试验数据,验证了涡动力学分析方法的实用性[10-12]。
本文借助FLUENT软件应用延迟分离涡(DDES)对离心泵进行非定常数值模拟,研究离心泵流道内非定常旋涡结构及涡结构分离、脱落等演化过程,以及不同工况下的旋涡分布情况。
单级离心泵基本参数:流量Q=160 m3/h,扬程H=45 m,额定转速n=2 950 r/min,叶轮出口直径D2=202 mm,叶轮出口宽度b2=20 mm,叶轮出口角β2=25°,叶片包角Φ=130°,叶片数Z=6。
数值计算基于延迟分离涡模拟方法(DDES),引入一个延迟函数,重新构造了DDES的长度尺度,同时考虑了网格尺度和涡黏场,避免了DES的模化应力损耗问题。DDES模型有基于S-A模型的DDES-SA方法和基于SST模型的DDES-SST方法。本文选用DDES-SST方法,其控制方程为:
DDES的长度尺寸为
lDDES=lRANS-fdmax(0,lRANS-lLES)
fd为延迟函数,其表达式为
fd=1-tanh[(8rd)3]
式中:vt、v分别表示涡和分子黏性;S表示应变率张量的大小;Ω表示涡量张量的大小;卡门常数κ=0.41;rd表示湍流尺度与壁面距离的比值。rd≪1时,fd=1,该区域运行LES算法,fd=0时,该区域运行RANS方法。
计算域由进口段、叶轮、蜗壳、出水段组成,使用Pro/E三维建模。网格划分采用ICEM 软件,为提高计算精度,整体采用六面体结构化网格,整体质量在0.5以上,角度18°以上。对叶片壁面网格进行加密,首层网格厚度为0.05 mm,共10层,比率为1.1,通过网格无关性验证,确定计算域网格总数为400万。叶轮、蜗壳及局部加密网格如图1所示。
进口段设为速度进口,出口段设为自由出流。采用无滑移壁面边界条件,近壁面采用标准壁面函数。延迟分离涡(DDES)模拟计算时采用的是全隐式时间格式,压力速度耦合采用SIMPLEC算法,二阶空间离散格式。
图1 叶轮和蜗壳网格
性能试验台如图2所示。
图2 离心泵性能试验台
将数值模拟计算得出的性能曲线与试验得出的性能曲线进行对比,如图3所示。两曲线在趋势上一致。额定工况点时,两者效率相差2.8%,扬程相差4.1%。这是因为数值模拟未考虑泄漏损失和机械损失,故模拟计算的结果高于试验值,误差在允许的范围内,因此该数学模型能够比较准确地预测泵的外特性,能为非定常计算提供保证。
图3 泵外特性曲线
涡量是描述流场的旋转运动特性的一个基本物理量,借助post后处理软件做出离心泵内部截面涡量分布图,Vorticity 取值范围为500~1500 s-1,有助于研究旋涡结构及涡结构分离、脱落等演化过程。
图4表示0.6Q工况下,RNGk-ε、DDES两种湍流模型模拟的中间截面涡量分布图。t0表示某一时刻,T表示叶轮转动一周的时间。t0时刻,旋涡主要集中在叶片背面和蜗壳进口附近,叶片背面附近的旋涡主要呈长条状,蜗壳进口附近的旋涡主要呈块状、团状。比较图4的t0时刻,RNGk-ε和DDES两种湍流模型模拟结果基本一致;t0+T时刻(即叶轮转动一周后),RNGk-ε和DDES两种湍流模型模拟结果差异很大。图4(d)叶轮流道内的旋涡从之前单一长条状变为许多小段,形成不对称的旋涡分裂,蜗壳进口附近的旋涡随着流体逆时针运动,之前块状的旋涡破裂成许多小碎块。图4(c)叶轮流道内的旋涡较之前无大的变化,基本一致,蜗壳进口附近的旋涡也随流体逆时针运动,但旋涡主要呈块状、团状。
通过对比RNGk-ε、DDES两种湍流模型模拟的中间截面涡量分布,发现都能清晰地模拟旋涡,但细节方面有很大差距。DDES湍流模型能很好地模拟旋涡分裂过程,清晰地看到由大旋涡分裂成的小旋涡,而从RNGk-ε湍流模型只能看到块状、团状的整体大旋涡,分裂后的小旋涡无法观测到。故DDES湍流模型较RNGk-ε湍流模型能更好地模拟旋涡演化过程。
图4 不同湍流模型的涡量分布
图5表示额定工况点中间截面(Y=0 mm)的涡量演化过程。叶轮叶片数Z=6,故取1/6T、2/6T、3/6T、4/6T、5/6T、T这6个时间段为研究对象。t0+1/6T时刻,旋涡主要集中在叶片背面和蜗壳进口附近,叶轮流道内的旋涡随着流体运动流入蜗壳,叶片背面附近的旋涡主要呈长条状,蜗壳进口附近的旋涡主要呈块状、团状;t0+2/6T、t0+3/6T时刻,叶轮流道内的旋涡变化不大,但蜗壳内的旋涡发生明显变化,蜗壳内的旋涡随着流体逆时针运动,之前块状的大旋涡破裂为许多连续的小旋涡,涡量同时也变小。当叶片经过隔舌位置后,其尾涡由于附带圆周速度,撞击隔舌,尾涡一分为二,一部分汇入蜗室,与蜗室内其他旋涡相融。同时,随着流体逆时针运动,继续分裂成许多小旋涡,涡量逐渐耗散,另一部分流向蜗壳出口,沿着壁面继续运动衰减;t0+4/6T、t0+5/6T时刻,叶轮流道内,叶片背面产生旋涡的区域逐渐减少;t0+T时刻,叶轮内只剩叶片出口尾涡,蜗壳内的涡量进一步减小,旋涡逐渐达到稳定分布。
图5 涡量的演化过程
图6表示额定工况点,不同截面的涡量分布。观察可得,不同截面的涡量分布有区别。对称截面的旋涡分布相似,但涡量值有区别;比较Y=3、6、9 mm截面的涡量云图,Y=3、6 mm截面的旋涡分布大致相似,旋涡从叶轮随流体进入蜗壳,蜗壳内的旋涡呈长条状,相切于叶片;而Y=9 mm截面内的旋涡主要集中在叶轮流道的出口边附近,蜗壳内的旋涡不明显。由此可得,蜗壳内的旋涡主要集中于Y=-6~6 mm区域之间。
图6 不同截面的涡量分布
图7表示运行稳定时,额定工况点同一时刻不同子午截面的旋涡分布图。从图中可以清晰地看到蜗壳内涡量的分布情况,蜗壳内的旋涡主要集中在叶轮和蜗壳交界面附近,蜗壳进口壁面有明显的附着涡。 b-b图右半部分中,叶轮叶片出口的尾涡进入蜗壳后,沿着径向继续运动,直至碰到蜗壳壁面,然后运动方向改变90°,沿着壁面继续运动,形成两个方向相反的大循环运动,截面越小更容易观察到,运动的过程中涡量逐渐减小。
图7 子午面涡量分布
图8表示0.6Q、1.0Q、1.4Q工况,中间截面的涡量分布图。一般非定常计算需要计算5~8个叶轮旋转周期方可得到可靠的解。本文选用第6个叶轮旋转周期。观察可得,流量不同,流道内涡量分布也不同,呈现一定规律,随着流量的增加,旋涡的分布范围逐渐减少,同时涡量也逐渐减小。对比3种工况,发现流道内的旋涡主要集中在叶轮叶片背面及蜗壳进口附近,隔舌处都有旋涡撞击,体现了旋涡分布的相似性。0.6Q工况时,叶轮和蜗壳流道内充满旋涡,2/3叶片长度的区域产生旋涡,形成类似卡门涡街的现象,涡量大,内部流动性差,随着流量的增加,叶片背面产生旋涡的区域逐渐减少,蜗壳内的旋涡相应减少。1.4Q工况时,流道内只有叶片出口切线状的尾涡。因为大尺度的旋涡不断地从主流获得能量,通过旋涡间的相互作用,能量逐渐向小尺度的旋涡传递,由于流体黏性的作用,小尺度的旋涡不断消失机械能耗散为流体的热能。涡量分布也反映了能量的耗散情况,所以性能曲线上0.6Q工况的泵效率远低于1.0Q、1.4Q工况时的泵效率。
图8 不同工况的涡量分布
图9表示蜗壳第三断面监测点位置。图10表示监测点的压力脉动特性。从时域图可以看出,蜗壳内部的压力呈周期性变化,在一个旋转周期内出现6个波峰和6个波谷,且这3个监测点的压力脉动频率与叶频保持一致,表明蜗壳内压力脉动的变化规律受叶轮旋转和叶片数影响。蜗壳内的压力脉动幅度在进口处最大,随着直径的增大,逐渐减弱。从频域图也可看出,蜗壳进口处压力脉动的剧烈程度最大。这说明叶轮和蜗壳间的动静干涉作用是诱发压力脉动的主要因素。
图9 蜗壳第三断面监测点
图11表示蜗壳第一断面、第八断面及隔舌监测点的位置。图12表示监测点的压力脉动特性。
图10 蜗壳第三断面的压力脉动特性
图11 蜗壳第一断面、第八断面及隔舌监测点
图12 蜗壳第一断面、第八断面及隔舌的压力脉动特性
隔舌(4点)处的压力系数Cp值最大,大约0.8,压力脉动程度剧烈;而5点、6点的压力系数Cp值,均值在0左右,压力脉动程度远小于隔舌点。这是因为叶轮转动经过隔舌位置时,从其流道内流出的流体由于附带圆周速度,故流体会撞击隔舌,同时蜗壳内的流体也会经过隔舌,隔舌处流动复杂,所以其压力脉动幅度剧烈。
1)在涡分裂过程中,相比于RNGk-ε湍流模型,DDES湍流模型能清晰地看到由大漩涡分裂而成的小漩涡,能很好地模拟精细流场结构,利于复杂流动的研究。
2)从涡量云图发现,流道内的旋涡主要集中在叶轮叶片背面、蜗壳进口及隔舌附近。在旋涡随流体运动的过程中,伴随着旋涡的初生、相融、迁移、破裂、耗散消亡等现象。流量对流道内旋涡的分布有明显影响,随着流量的增加,旋涡逐渐减少,涡量也逐渐减小。对于该泵,蜗壳内的旋涡主要集中于Y=-6~6 mm区域之间。
3)蜗壳内的压力脉动在进口处最大,随着直径的增大,逐渐减弱;隔舌位置的压力脉动远远大于蜗壳的其他位置。