时间步长对核主泵非定常计算精度的影响机理*

刘 毅,王秀勇,2*,董 峰,杜永峰

(1.兰州理工大学 能源与动力工程学院,甘肃 兰州 730050;2.甘肃省流体机械及系统重点实验室,甘肃 兰州 730050)

0 引 言

近年来,应用计算流体动力学(computational fluid dynamics,CFD)技术对泵进行水力性能预测及优化,已经成为泵类产品前期研发的重要方法[1,2]。CFD的计算精度一直是被关注的重点,同时也是计算流体动力学的核心问题。

在计算域模型表达、网格划分方式、壁面函数应用、湍流模型选取、边界条件设置、流动状态是否定常等方面,国内外学者已经对影响数值计算精度的因素进行了相关研究[3,4]。李晓俊等人[5]研究了网格质量对离心泵数值模拟精度的影响,发现湍流模型的选取与边界层网格质量密切相关。刘宇宁等人[6]对多级离心泵水力性能数值模拟精度的影响因素进行了研究,发现非定常计算的精度高于定常计算。

在研究泵的非定常流动特性时,时间步长是影响其计算精度的一个重要因素,通常要将叶轮每旋转1°～3°所需要的时间设置为时间步长[7-9]。朱荣生等人[10]对核主泵小流量工况下压力脉动特性进行研究,将叶轮每旋转2°作为时间步长进行了非定常计算。王东伟等人[11]在研究离心泵非定常空化流场时,将叶轮每旋转3°作为一个时间步长。但目前并没有研究对非定常计算过程中如何确定合理时间步长这一问题进行系统分析。

众所周知,减小时间步长能够提高非定常计算的精度,而目前的CFD技术还无法完全真实还原泵内的实际流场结构,也就是说,数值计算始终会存在一定的误差[12]。

为了研究时间步长与计算精度之间的关系,首先需要明确泵内与时间密切关联的主要流动现象。由于各种频率下的湍流脉动是影响核主泵内流场分布特征的主要因素,由此可以推断得出,时间步长所对应的响应频率对湍流脉动频率谱的覆盖率,才是影响非定常计算精度的主要因素。因此,有必要对非定常计算精度随时间步长减小时的变化规律展开探讨研究。

笔者以核主泵为研究对象,对全计算域进行六面体结构化网格划分,应用FLUENT软件,选择基于剪切应力输运(shear stress transport,SST)k-ω模型的延迟分离涡模拟(delayed detached-eddy simulation,DDES)进行非定常计算[13],分析时间步长变化对核主泵水力性能的计算精度、瞬时性能曲线的分布特征及旋涡场的影响,并在此基础上,探讨时间步长对核主泵内复杂流动特征的解析效应,从而为确定合理的时间步长提供理论参考。

1 计算模型及数值计算方法

1.1 计算模型

笔者以某型号核主泵原型样机为研究对象,其设计工况点参数为:流量Qd=630 m3/h,扬程H=35 m,转速n=990 r/min,比转速ns=105,叶轮的叶片数为6枚,导叶的叶片数为10枚。

计算域包括进水段、叶轮、导叶、压水室、前腔、口环间隙及后腔等过流区域。

整体几何模型如图1所示。

图1 计算域模型分解图

1.2 网格划分

笔者应用ICEM软件对计算域进行六面体结构化网格划分。计算域网格划分情况如图2所示。

图2 计算域网格划分情况

此时,计算域网格总数量达到1.5×107,完全满足当模型泵的网格数为9.42×106时,即可实现网格无关性的要求[14]。

近壁面第一层网格尺度用无量纲壁面距离y+来表示。

各过流部件近壁面y+值的分布情况,如图3所示。

图3 壁面y+值分布图

其中,y+≤5的网格数量占全计算域近壁面网格总数的65.77%;510的网格数量为16.76%。

由于该泵的三维模型在结构上比较复杂,在应用ICEM划分网格时,若近壁面第一层网格的尺度在现在的基础上进一步减小,则网格质量变得较差,不能满足计算对网格质量的要求。

受软件网格划分功能的限制,暂且还没有更好的办法来进一步降低y+值,但目前已经具有的y+值分布也基本能够满足DDES模型的要求。

1.3 边界条件设置

笔者采用动量方程及连续性方程对整个流场进行求解,湍流模型选择基于SSTk-ω的DDES,进口采用速度进口,出口采用自由出流边界条件,各流体域之间的耦合面设置为interface。

笔者采用有限体积法对控制方程进行离散,用SIMPLEC算法实现压力与速度之间的耦合,壁面采用无滑移壁面边界条件。

进行定常计算时,动静域交界面之间采用多重参考坐标系模型;非定常计算时,以定常计算结果作为初始化流场,动静域交界面采用滑移网格模型。

笔者分别将叶轮每旋转Δα=1°、1.5°、2°和3°所需要的时间设置为时间步长,以Δα=1°为时间步长完成整个计算所需要的时间,比Δα=1.5°所需的时间增加约50%,近似为Δα=2°所需时间的2倍,Δα=3°所需时间的3倍;

每个时间步长设置迭代20次,收敛精度为10-5;在叶轮旋转第5圈时,泵出口处总压的时均值不再发生明显变化,采用叶轮旋转第5圈的数据进行后处理分析。

2 水力性能计算结果分析

2.1 计算结果与试验结果对比

在设计工况点下,4种时间步长的水力性能计算结果与试验结果对比情况,如表1所示。

表1 各时间步长非定常计算结果与试验结果的对比

表1中,各性能参数的计算结果都为时均值,即叶轮旋转一周的过程中,各性能参数在每个时间步长下监测结果的平均值。

(1)由扬程的计算结果与试验结果对比可知:当Δα≤1.5°时,扬程的计算值大于试验值;Δα≥2.0°时,则扬程的计算值小于试验值;并且,当Δα=1.5°时,扬程的计算值最大,之后随着时间步长的增大,扬程的计算值持续减小;Δα≤2°时,扬程在各时间步长下的相对计算误差大小差异不大,其中,Δα=1.5°的计算误差相对较大,为0.40%;Δα=3.0°时,扬程的相对计算误差最大,为1.14%;

(2)由轴功率的计算结果与试验结果对比可知:Δα=1.0°时的计算值略高于试验值,从Δα=1.5°开始计算值低于试验值,此时的相对计算误差最小,为0.1%;之后,随着时间步长的增大,轴功率的计算值持续减小,而相对计算误差则持续增大,在Δα=3°时最大,为1.02%;

(3)由效率的计算结果与试验结果对比可知:Δα=1.5°和Δα=2.0°时的计算值高于试验值,而Δα=1.0°和Δα=3.0°时的计算值则略低于试验值(其中,Δα=1.5°时效率的计算值最大,之后随着时间步长的增大,效率的计算值逐渐减小;Δα=1.5°时效率的相对计算误差最大,为0.48%,其他3个时间步长下,效率的相对计算误差位于0.12%～0.18%之间)。

相对计算误差大小如图4所示。

图4 不同时间步长相对计算误差

由图4综合对比来看:随着时间步长的减小,各性能参数的相对计算误差没有明确的变化规律,也就是单个性能参数的计算精度并不一定随着时间步长的减小而呈现持续提高的趋势;但Δα≤2°时,扬程和轴功率的计算误差明显低于Δα=3.0°,而效率在所有时间步长下的计算误差均较小;

因此,综合各性能参数的计算误差可以判断:在对泵进行非定常计算时,时间步长不大于叶轮每旋转2°所需要的时间,即可获得较高的计算精度;这意味着时间步长的大小存在临界值,当时间步长小于该临界值时,可以在总体上获得较高的计算精度。

扬程、轴功率和效率在各时间步长下的相对计算误差综合分析情况,如表2所示。

表2 各时间步长相对计算误差分析

由表2总体来看:时间步长越小,扬程、轴功率和效率三者计算误差的平均值越小,误差的离散度也越小;即时间步长越小,各性能参数的综合计算精度越高;

根据各性能参数计算误差的平均值和均方差在不同时间步长之间的差异来看:Δα≤2°时,计算误差的平均值及均方差都明显小于Δα=3.0°;其中,Δα=2°时,计算误差的平均值比Δα=3°减小0.43%,均方差减小0.21%;而计算误差的平均值和均方差在Δα=1°、Δα=1.5°和Δα=2°之间的差异较小,平均值最多只相差0.08%,均方差最多只相差0.09%。

由此可知:当Δα≤2°时,虽然时间步长变小会在总体上提高计算的精度,但提高幅度不明显;

上述分析结果进一步表明:Δα=2°是核主泵非定常计算精度发生质变的一个关键临界值;时间步长大于叶轮每旋转2°所需要的时间时,水力性能的计算精度会明显降低;

由上述分析可知:当时间步长小于某个临界值时,随着时间步长的减小,单个性能参数的计算精度不一定会持续提高。虽然所有性能参数的综合计算精度会有所提高,但提高幅度较小。由于完成整个非定常计算所需要的时间与时间步长成反比的关系,当时间步长减为原来的一半,完成整个计算所需要的时间比原来增加几乎一倍。

综合考虑计算精度和计算所消耗时间这两个因素,笔者将叶轮每旋转Δα=2°所需要的时间设置为非定常计算的时间步长,这样既可以满足一般工程应用的计算精度要求,又可以大幅节约计算所消耗的时间。

2.2 水力性能时域分析

笔者在设计工况点进行非定常计算时,采用了4种不同的时间步长,但从计算误差的平均值和均方差在不同时间步长之间的对比分析来看,Δα=1.5°和Δα=2.0°二者之间具有很高的相似度。故在后续的研究内容中,笔者仅选取Δα=1°、Δα=2°和Δα=3°的计算结果进行分析。

不同时间步长下水力性能计算值对应的时域图,如图5所示。

图5 不同时间步长下水力性能时域图T—叶轮旋转360°所用时间;t—叶轮旋转60°所用时间

图5中,从水力性能是否为定常的角度来看:

在叶轮旋转的过程中,泵的扬程、轴功率和效率并非保持为常数状态,而是在一定的区间范围内持续脉动(例如,当Δα=3°时,扬程的脉动幅度为时均值的±1.7%左右,轴功率为±2.8%左右,效率为±1.2%左右;当Δα=2°时,各水力性能的脉动幅度与Δα=3°相接近;而当Δα=1°时,扬程和轴功率的脉动幅度较前二者明显偏小,但效率略微偏大,其中,扬程在时均值附近的脉动幅度为±1.1%左右,轴功率为±2.1%左右,效率为±1.4%左右);

由于水力性能并非定常状态,因而对泵的水力性能进行试验测试时,针对一个工况点进行多次采样时,取平均值是非常重要的。

由各性能参数在不同时间步长之间的对比来看:

在各时间步长下,扬程和轴功率的脉动周期均不明显,而效率的脉动周期相对明显一些,波峰和波谷分别交替出现30次,为叶轮和导叶的干涉周期;扬程和轴功率的脉动频率随时间步长的减小而明显加快,脉动幅度始终在变,同时曲线分布重心上移,也就是时均值随时间步长的减小而增大;效率的脉动频率在各时间步长下总体都保持一致,曲线分布重心也基本相同;只是,当Δα=1°时,效率的瞬时分布曲线出现了更多的微小波动特征。

总体来看:当时间步长减小时,各性能曲线在分布特征上能表现出更多的、小尺度的波动特性,这也正是核主泵内流动特征复杂的细微体现;导叶和压水室内旋涡运动尤其明显,非定常流动特征加强,从而对核主泵的水力性能产生明显的影响,导致水力性能周期性变化的特征受到干涉,脉动周期不再明显,该现象是由泵的内因引起的,并非数值计算方法的问题。

由于时间步长减小后各水力性能的瞬时曲线在分布特征上会出现更多的小幅度的波动特性,这意味着此时的时间步长对核主泵内更细微的非定常流动特征的分辨能力有所提升,能够识别出更高频率下的湍流脉动对流动的影响,因而计算精度会更高一些;这也正是Δα=2°时,各水力性能的综合计算精度明显高于Δα=3°的原因。

由于湍流的脉动频率主要集中在某段区间内(例如自由射流的湍流脉动频率谱主要分布在3 500 Hz以内,主频率在400 Hz左右[15，16]),当时间步长的响应频率(f=6n/Δα,n为核主泵的转速,r/min)大于该区间的上限时,并不会使计算精度得到明显提高;正如Δα=1°(f=5 940 Hz)时,各水力性能的综合计算精度略高于Δα=2°(f=2 970 Hz),但提高幅度却没有Δα=2°相对于Δα=3°(f=1 980Hz)那么明显;因为此时Δα=2°所对应时间步长的响应频率正好位于湍流脉动频率的上限附近,也就是非定常计算的时间步长存在临界值,当时间步长小于该临界值后,对泵内非定常流动特征的响应频率即可满足高精度计算的要求。但随时间步长的继续,减小计算精度的提高幅度却不再那么明显。

3 内流场分析

不同时间步长下动静叶栅内的旋涡分布状态,如图6所示。

图6 动静叶栅内旋涡分布图

不同时间步长下压水室内的旋涡分布状态,如图7所示。

图7 压水室内旋涡分布图

此处,笔者采用Q准则对泵内的涡结构进行提取[17，18],并选取叶轮完成第5个旋转周期末的结果,对流场进行定性分析。

由图6可知:叶轮内的旋涡运动在不同时间步长下的分布情况非常相似,流场结构相对简单,流动的核心区域几乎没有旋涡运动,只在叶片背面及叶轮外缘附近存在两条细长的带状涡(其中,紧靠叶片背面的带状涡较窄,而贴近叶轮流道出口外缘的带状涡相对较宽);受导叶叶片的干涉作用,带状涡在导叶的叶片入口处开始破裂,并流入导叶内部;由于叶轮和导叶的动静干涉具有径向对称的效果,叶轮内的旋涡分布状态也呈现出良好的径向对称性。

导叶内的流场结构比较复杂,整个流道空间几乎都被大量形状各异、分布不均的小尺度旋涡占据;旋涡在导叶流道入口处呈现为带状结构,沿着流动方向,带状涡逐渐发展,破碎后成为各种形状的小尺度旋涡;受叶轮干涉作用时序效应的影响,导叶相邻两流道内的旋涡场结构并不完全相同;虽然叶轮和导叶的动静干涉具有径向对称的效果,但导叶内的旋涡分布状态并没有呈现出良好的对称性(原因在于压水室内的流场结构是非径向对称的,如图7所示,导叶每个流道的出口边界条件均不相同,使各个导叶流道内的旋涡场结构出现差异)。

总体来看:导叶内的旋涡场在不同时间步长下的宏观结构比较相似,但随着时间步长的减小,捕捉到的涡结构尺度更小、数量更多,且分布更为离散;与Δα=3°的涡结构相比,Δα=2°的涡结构明显与Δα=1°的相似度更高,由此也可表明Δα=2°时的时间步长具有足够精度来解析湍流脉动对导叶内流场的影响。

由图7压水室内的旋涡分布可知:由于叶轮沿逆时针方向旋转,压水室内液体的循环方向也为逆时针,内循环的液体对由导叶左侧射入压水室的液体产生较强的冲击作用,从而使压水室左侧的旋涡运动强于右侧;

Δα=1°与Δα=2°时间步长下的旋涡分布十分相似,且都比Δα=3°时的旋涡数量多,强度也略大,这进一步表明,较小的时间步长能够捕捉到更细致的旋涡运动特征,同时也表明,Δα=2°时的时间步长具有足够的精度来解析压水室内的湍流运动情况。

湍流脉动是由各种不同尺度和不同旋向的旋涡相互叠加共同作用的结果,大尺度旋涡主要引起低频脉动,小尺度旋涡主要引起高频脉动,所以时间步长能够反映的频率谱越宽,就能涉及更多小尺度旋涡对流动的影响,数值计算的精度也就越高;

由不同时间步长下,各过流部件内旋涡场分布情况对比可知:相较于Δα=3°,Δα=2°与Δα=1°所对应的时间步长,均能分辨出更小尺度的旋涡运动,这意味着较小的时间步长能够解析更高频率下的湍流脉动情况,因而其计算精度也更高;由于Δα=1°时所捕捉的旋涡场中小尺度旋涡的数量比Δα=2°略多一些,因而其计算精度也比Δα=2°略高一些;同时可以进一步明确,Δα=2°所对应的时间步长是核主泵非定常计算精度明显提高的临界值[19]。

综上所述,在对核主泵展开非定常计算研究时,时间步长确实存在一个临界值,该临界值对应的响应频率能够覆盖核主泵内湍流脉动频率谱的绝大多数频段,从而使计算精度明显提高;

当时间步长小于该临界值时,响应频率虽然会进一步提高,但由于湍流脉动高频段所占的比例极低,所以其计算精度提升的幅度有限。

4 结束语

笔者通过改变非定常计算过程中的时间步长,对比分析了核主泵水力性能在计算精度、水力性能的脉动特性,及旋涡场分布特征等方面的变化异同点,并探讨研究了时间步长对核主泵非定常计算精度的影响机理,确定了时间步长的合理取值范围。

研究结果如下:

(1)随时间步长的减小,单项水力性能的计算精度不一定能提高,但各水力性能的综合计算精度能得以提高;时间步长对应的叶轮旋转角度Δα=2°是非定常计算精度能够明显提高的临界值;当Δα≤2°时,随时间步长的减小,非定常计算精度的提高幅度较小;

(2)在叶轮旋转的过程中,各水力性能的时域图呈现出脉动的特性;扬程和轴功率的脉动周期不规律,而效率则呈现出周期性脉动的特征;随着时间步长的减小,扬程和轴功率的脉动频率加快,效率的脉动频率始终保持不变;当Δα≤2°时,各水力性能的瞬时曲线在几何特征上出现了小尺度波动的形状,表明此时的时间步长能够反映流场中更细微的流动特征;

(3)时间步长减小能够捕捉到更小尺度的旋涡运动,反映出更高频率的湍流脉动对流动的影响;当Δα=2°时,导叶和压水室内的涡量场结构与Δα=1°时相似,小尺度涡的数量均比Δα=3°多,并且分布状态更加离散;

(4)时间步长的响应频率所能覆盖的湍流脉动频率谱的范围,决定了非定常计算的精度,覆盖范围越宽,计算精度也越高;以叶轮每旋转Δα=2°所需要的时间作为时间步长,既可以获得较高的非定常计算精度,又可以节约计算所消耗的时间。

在后续的工作中,笔者将对核主泵内影响湍流脉动频率的主要因素,以及泵内湍流脉动频率区间的确定展开深入研究,为准确定位非定常计算的精度提供理论支撑。