多级离心泵水力性能数值模拟精度影响因素研究

刘宇宁，王秀勇，刘志远，白小榜

（1.兰州理工大学 能源与动力工程学院，兰州 730050；2.重庆水泵厂有限责任公司国家级企业技术中心，重庆 400030）

0 引言

多级离心泵在工农业生产中的应用极为广泛，由于多级泵水力模型的开发制造成本较高，因而在设计过程中通常借助于CFD技术对其水力性能进行预测，并对泵内的流场结构进行分析研究，为水力性能的进一步优化提供有效的参考信息［1］。虽然CFD技术在流体机械行业中已被广泛应用，但数值计算精度问题依然是当前研究的重点内容。由于目前关于多级离心泵的数值模拟研究多集中于泵内流动特征的分析，而关于提高水力性能预测精度的研究还很少，因而在既有理论基础上，有必要对影响多级离心泵数值计算精度的因素进行研究，从而为提高水力性能的预测精度提供可靠而又可行的数值模拟方法。

研究结果表明，影响数值计算精度的主要因素有流动状态是否定常、湍流模型的选取以及近壁面网格尺度等［2-5］。在工程应用中，关于流动控制方程的求解通常采用雷诺时均化（RANS）方法，该方法具有较强的计算经济性，对计算机资源要求不高，在离心泵数值模拟中的应用最为广泛［6-7］，但由于时均化过程会导致部分湍流脉动信息丢失，流场保真性较差；大涡模拟（LES）方法在反映湍流脉动特征方面具有RANS无可比拟的优势，常用于捕捉泵内部的涡核分布［8］，计算精度较高，但对计算资源的要求也很高；近年来，结合了LES与RANS各自优点的分离涡模拟（DES）方法在流体机械的数值模拟研究中取得了良好的效果［9-10］，但计算过程中可能出现模化应力不足问题，而延迟分离涡模拟（DDES）方法可以对该问题进行修正，并且在离心泵数值模拟中的应用效果比较明显［11］，有望成为今后离心泵数值模拟所选用湍流模型的主流。

由于目前探讨离心泵数值计算精度影响因素的研究对象主要集中于单级单吸结构，为了系统地研究各因素对多级离心泵数值计算精度的影响，本文以某型号的多级离心泵为研究对象，采用ANSYS FLUENT计算软件，首先对基于RNGk-ε湍流模型的定常和非定常计算的结果进行对比；然后比较分析RNGk-ε和DDES两种湍流模型在多级离心泵数值模拟中的应用效果；最后探讨叶轮近壁面网格尺度的不同对数值计算结果的影响。通过上述研究，为多级离心泵的数值模拟方法提供有价值的参考。

1 数值计算方法

1.1 计算模型

本文以并联结构的节段式多级泵为研究对象，结合试验泵的实际情况，级数取为两级，其设计参数为：流量Q=2 400 m3/h，扬程H=190 m，转速n=1 480r/min，首级和次级叶轮的叶片数均为6枚，导叶的叶片数为8枚。由于并联式结构关于末级压出室中心面两侧对称，为减少数值计算的工作量，计算域在此取其中一侧，包括吸入室、首级叶轮、次级叶轮、导叶、压出室、叶轮前后腔以及口环间隙，整体几何模型如图1所示。

应用ICEM CFD进行网格划分。由于泵体结构复杂，若进行六面体结构化网格划分，则工作量巨大，耗时很长，技术要求极高，难以达到在一般工程应用中进行普遍推广的目的，因而对吸入室、导叶和压出室采用容易实施的四面体非结构网格划分方案。叶轮是影响泵性能好坏的核心部件，叶轮内流场结构的模拟精度是决定泵整体性能计算精度的关键，为了研究叶轮近壁面处的网格尺度对数值计算精度的影响，对叶轮采用两种网格划分方案，一种是进行四面体非结构化网格划分，另一种是进行六面体结构化网格划分并加边界层。叶轮的网格划分情况如图2所示。

图1 几何模型

图2 叶轮网格划分

根据网格无关性验证结果，当网格总数增加到1 100×104后，扬程和效率的计算值趋于稳定，此时吸入室、导叶和压出室的网格尺度设定为12 mm；叶轮采用四面体网格划分时设定叶片表面网格尺度为6 mm，全局网格尺度为8 mm；叶轮采用六面体网格划分时，设定近壁面第一层网格的高度为0.05 mm，网格尺度增长比率为1.2。

1.2 湍流模型

RNGk-ε湍流模型［12］采用重正化群（RNG）的统计技术，它会根据不同强度的旋流来适当地改变湍流黏度，进而对旋流流动进行修正，同时它也兼顾了平均流动中的旋转及旋流情况。

DDES湍流模型［13］是在DES的基础上，将延迟函数引入到方程当中，混合长度尺度被重新构建，从而解决了DES模型中由网格诱导产生的分离流动（Grid Induced Separation，GIS）问题，进而减弱了DES对网格的依赖性。SST-DDES对长度尺度的修正公式为：

式中lDDES——DDES混合长度尺度，m；

lRANS——RANS混合长度尺度，m；

fd——DDES模型的延迟函数；

CDES——DES函数校准常数，CDES=0.61；

Δmax——本地网格单元最大尺度，m；

Cd1，Cd2——经验常数，Cd1=20，Cd2=3；

rd——当地湍流尺度与到壁面距离的比值；

v——运动黏度，m2/s；

vt——旋涡运动黏度，m2/s；

k——卡门通用常数，k=0.412；

d——网格单元到壁面距离，m；

S——应变率张量，s-1；

Ω——涡量张量，s-1。

1.3 计算方法

采用有限体积法对控制方程进行离散。压力和速度的耦合采用SIMPLEC算法。定常计算时，采用动参考系模型；非定常计算时，采用滑移网格模型，并使用前面的定常计算结果作为初始流场。壁面采用无滑移边界条件；入口给定速度边界条件；出口给定自由出流边界条件；静止域与旋转域之间采用Interface边界条件进行数据传递。将叶轮每转过3°圆周角所需要的时间定义为时间步长，其值为3.378 38×10-4s。定常计算时，当泵出口总压趋于稳定时，可视为计算结果收敛；非定常计算时，当泵出口总压呈现周期性波动时，取物理量在最后一个旋转周期内的时均值进行性能计算。

2 计算结果分析

2.1 定常与非定常计算结果对比

图3示出全计算域采用四面体网格时，在全流量工况范围内，RNGk-ε湍流模型分别采用定常与非定常2种计算方法预测得到的扬程、轴功率和效率的计算结果。

图3 RNG k-ε湍流模型定常与非定常计算结果比较

由图3（a）所示的扬程计算结果来看，定常与非定常计算在全流量工况范围内的结果有明显的不同，其中非定常计算的扬程在各个工况点均高于试验值，但扬程计算曲线随流量增加时的变化趋势与试验曲线吻合良好；而定常计算的扬程在（0～0.3）Qd以及（1.1～1.2）Qd工况范围内高于试验值，在（0.3～1.1）Qd工况范围内低于试验值，扬程的计算曲线与试验曲线存在两个交点，导致扬程计算曲线随流量增加时的变化趋势与试验曲线不完全一致，在（0～0.8）Qd工况范围内下降速度较快，在（0.8～1.2）Qd工况范围内下降速度又相对缓慢。从扬程相对计算误差大小的角度来看，在（0.2～0.5）Qd以及（0.9～1.2）Qd工况范围内，定常计算的扬程计算误差小于非定常计算，原因是这两个工况范围处于定常计算曲线与试验曲线的交点位置附近，当偏离这两个工况范围时，定常计算的相对计算误差大于非定常计算；定常计算的扬程相对计算误差最大值在关死工况点，为12.23%，其次是0.8Qd工况点，为9.10%，其他工况点的计算误差位于2.54%～6.61%之间；非定常计算的扬程相对计算误差最大值在1.2Qd工况点，为10.47%，其他工况点的计算误差相差不大，均位于4.13%～7.32%之间。

由图3（b）所示的轴功率计算结果来看，定常与非定常计算的结果差异明显，轴功率的非定常计算值在全流量工况点都高于试验值，但其计算曲线随流量增加时的变化趋势与试验曲线吻合良好，能够准确反映出轴功率与流量之间的特性关系；而轴功率的定常计算值在0.6Qd工况点之后小于试验值，之前大于试验值，并且在关死工况点附近远大于试验值，导致定常计算曲线的变化趋势仅在（0.6～1.2）Qd工况范围内与试验曲线基本一致，但在（0～0.6）Qd工况范围内随流量的增加呈现下降的趋势，与试验曲线的变化趋势相反，不能正确反映轴功率在小流量工况区间内的变化特征。从轴功率相对计算误差大小的角度来看，定常与非定常计算均在关死工况点存在最大计算误差，其中定常计算的相对计算误差为53.36%，非定常计算为22.80%；定常计算的最小误差在0.6Qd工况点，仅有0.09%，该工况点刚好是轴功率的定常计算曲线与试验曲线的交点；非定常计算的误差随着流量的增加逐渐减小，在1.2Qd工况点达到最小值，为5.06%。

由图3（c）所示的效率计算结果来看，定常与非定常计算的效率计算值均在设计工况点之前小于试验值，之后大于试验值；在（0～0.8）Qd工况范围内，效率的非定常计算曲线与试验曲线的吻合度明显高于定常计算；在（0.8～1.2）Qd工况范围内，效率的定常与非定常计算结果比较接近，随着流量的增加效率的计算值也在增大，而试验值则变化平缓，效率的计算曲线与试验曲线在设计工况附近有明显的区别。从效率相对计算误差大小的角度来看，定常与非定常计算的最大误差均在0.2Qd工况点，其中定常计算的相对误差大小为23.72%，非定常计算为7.36%；定常计算误差的次高点在0.3Qd工况点，为14.89%，其他工况点相对计算误差的大小位于1.95%～7.25%之间；非定常计算误差的次高点在1.2Qd工况点，为6.22%，其他工况点相对计算误差的大小位于1.98%～3.97%之间。

总体来看，在全流量工况范围内，多级泵各性能参数的非定常计算曲线与试验曲线的吻合度都远高于定常计算，尤其是效率的非定常计算结果，无论是在计算精度方面还是在计算曲线随流量的变化趋势方面，都具有定常计算不可比拟的明显优势。因此，对多级泵进行性能预测时，采用非定常计算的方法可以获得更精确的计算结果。

2.2 RNG k-ε与DDES湍流模型计算结果对比

图4示出全计算域采用四面体网格时，在全流量工况范围内进行非定常计算的前提下，RNGk-ε与DDES 2种湍流模型关于扬程、轴功率和效率的计算结果。

图4 RNG k-ε与DDES湍流模型非定常计算结果比较

由图4（a）所示的扬程计算结果来看，在全流量工况范围内，DDES湍流模型的扬程计算值明显低于RNGk-ε湍流模型，其计算曲线与试验曲线之间具有极高的吻合度，计算精度远高于RNGk-ε湍流模型。DDES的相对计算误差最大值在1.2Qd工况点，为2.48%，其次是0.4Qd工况点，为2.09%，在其他工况点的相对计算误差值均小于1.30%。

由图4（b）所示的轴功率计算结果来看，在全流量工况范围内，DDES湍流模型的轴功率计算值明显低于RNGk-ε湍流模型，并且也都低于试验值，但二者的轴功率计算曲线随流量增加时的变化趋势仍然保持较高的相似度。从轴功率相对计算误差大小的角度来看，DDES湍流模型在全流量工况范围内的计算误差都明显低于RNGk-ε湍流模型，其相对计算误差的最小值在关死工况点，为0.05%，其次是1.2Qd工况点，为1.95%，在（0.2～1.0）Qd工况范围内的相对计算误差值比较稳定，基本维持在4.22%左右。

由图4（c）所示的效率计算结果来看，在全流量工况范围内，DDES湍流模型的效率计算值明显高于RNGk-ε湍流模型，并且也都高于试验值，但DDES湍流模型的效率计算曲线在设计工况点附近更为平坦一些，可以更准确地反映效率随流量增加时的变化趋势。从效率相对计算误差大小的角度来看，在全流量工况范围内，DDES湍流模型的计算误差值均高于RNGk-ε湍流模型，其计算误差的平均值比RNGk-ε高出1.7%左右。

总体来看，当计算域全局采用四面体网格时，DDES湍流模型在多级泵的扬程和轴功率的计算精度上比RNGk-ε湍流模型具有明显的优势，虽然在预测效率计算曲线随流量增加时的变化趋势上也具有一定的优势，但在计算精度上与RNGk-ε湍流模型相比却处于明显的劣势状态。实际上，DDES湍流模型对近壁面的网格尺度具有较高的要求，而四面体网格划分的方法远达不到其要求，因而下面将进一步探讨近壁面网格尺度的大小对DDES湍流模型计算结果的影响。

2.3 叶轮近壁面网格尺度对计算结果的影响

吸入室、导叶和压出室均采用四面体网格，叶轮分别采用四面体（无边界层）和六面体（有边界层）网格进行数值计算。根据数值计算的统计结果，叶轮近壁面网格的y+值分布情况如图5所示，从图可以看出，网格无边界层时的y+值集中分布在250～1 000的区间内，有边界层时的y+值集中分布在10～30的区间内。虽然叶轮的两种网格划分方案的y+值均不满足DDES湍流模型y+＜5的要求，但仍可以用来定性讨论分析近壁面网格尺度的变化对计算结果的影响，更细致的研究将在后期展开。

图5 不同网格划分方案近壁面y+值分布情况

图6示出采用非定常计算时，DDES湍流模型在叶轮有无边界层网格的情况下关于扬程、轴功率和效率的计算结果。

图6 DDES湍流模型有无边界层计算结果比较

由图6（a）所示的扬程计算结果来看，当在叶轮壁面附近加入边界层使y+值大幅减小时，扬程的计算曲线与试验曲线几乎完全重合，与无边界层的计算结果相比，加入边界层后扬程的计算精度有了明显提高。叶轮网格加入边界层后，扬程的最大相对计算误差在1.2Qd工况点，仅有1.43%，其次是设计工况点，为1.06%，其他工况点的相对计算误差值均在0.50%左右。

由图6（b）所示的轴功率计算结果来看，在叶轮壁面附近加入边界层时，轴功率的计算值比无边界层时有所提高，但在0.2Qd工况点之后仍小于试验值。在考虑了边界层对计算结果的影响后，除了在关死点轴功率的相对计算误差比无边界层时增加了1.17%外，在其他工况点的相对计算误差均明显减小，最低可以减小0.55%，最高减小了2.87%，平均能够减小1.20%。

由图6（c）所示的效率计算结果来看，有边界层时的效率计算曲线与试验曲线的吻合度进一步提高，在全流量工况范围内比无边界层时的计算结果更精确，除了关死工况点的效率在理论上不存在计算误差外，其他工况点有边界层时的相对计算误差均比无边界层时减小了2.50%左右。

总体来看，当减小叶轮近壁面处的网格尺度时，DDES湍流模型关于多级泵各性能参数的计算精度均有明显提高，并且可以推测，当叶轮尤其是所有过流部件近壁面网格尺度的y+值完全满足DDES湍流模型的要求时，其计算精度会有进一步的提高。

2.4 各数值计算方法计算误差综合分析

为了综合比较上述各数值计算方法的优劣，在此对各性能参数在全流量工况范围内相对计算误差大小的平均值及相对计算误差的均方差进行统计，如表1所示。

表1 全流量工况点计算误差分析 （%）

相对计算误差大小的平均值可以反映计算值与试验值之间的偏离程度，而均方差则可以反映计算曲线与试验曲线在曲线形状上的相似度；平均值越小，则计算值越接近试验值；均方差越小，则计算曲线随流量增加时的变化趋势越接近试验曲线。

由表1所示各计算方法相对计算误差的对比来看，在全流量工况点，定常计算方法的计算误差最大，计算曲线与试验曲线在曲线形状的相似度上也最低，而非定常计算方法在上述两个问题上都有明显改善。计算域全局采用四面体网格并且进行非定常计算时，虽然DDES湍流模型关于效率的计算误差比RNGk-ε湍流模型高出了1.72%，但关于扬程、轴功率的计算误差却都比RNGk-ε湍流模型低了4%，并且各性能参数相对计算误差的均方差也都比RNGk-ε湍流模型低很多，也就是DDES湍流模型在预测性能曲线随流量增加时的变化趋势上比RNGk-ε湍流模型要精确很多，所以综合来看，在相同的计算方法下，DDES湍流模型比RNGk-ε湍流模型更适合于预测多级泵的性能参数。当采用DDES湍流模型进行非定常计算时，在叶轮近壁面区域加入边界层网格后，各性能参数的计算误差进一步减小，计算曲线与试验曲线的吻合度进一步提高。

综上所述，对多级泵的性能参数进行预测时，采用DDES湍流模型和非定常计算方法，能够获得较高的计算精度，尤其是对叶轮的近壁面网格进行加密处理时，可以使计算精度进一步提高。在该数值计算方法中，网格划分容易实现，网格数量容易控制，非定常计算所需要的时间也仅为定常计算的4～5倍，但计算精度却比常规的采用RNGk-ε湍流模型和定常计算的方法高出很多，适合于一般工程应用。

3 结论

（1）在全流量工况点，非定常计算方法不仅计算误差比定常计算方法小，并且性能计算曲线的形状也更接近于试验曲线；定常计算时在个别工况点的计算精度极高具有一定的偶然性，与该工况点距离计算曲线和试验曲线交点位置的远近有关。

（2）全计算域采用四面体网格和非定常计算的条件下，DDES湍流模型关于效率的计算精度稍低于RNGk-ε湍流模型，但关于扬程和轴功率的计算精度却远高于RNGk-ε湍流模型，并且各性能参数计算曲线的形状更接近于试验曲线。

（3）当对近壁面网格进行加密处理时，DDES湍流模型的计算精度可以进一步提高。