基于IBM-LES的贯流泵湍流特性

阚阚，张清滢，黄佳程，李昊宇，郑源，陈龙

(1. 河海大学能源与电气学院，江苏 南京，211100； 2. 南通河海大学海洋与近海工程研究院, 江苏 南通，226019)

贯流泵具有结构紧凑、过流能力强等诸多优点，被广泛应用于中国平原地区的防洪排涝、农业灌溉等水利工程和低扬程水力系统.近年来，许多学者通过数值模拟和试验的方法，对其复杂水力特性进行了深入的研究.石丽建等[1]对考虑回流间隙的全贯流泵装置进行了计算流体动力学(CFD)数值模拟，研究发现，回流间隙对全贯流泵的水力性能的负面影响显著.谢荣盛等[2]通过数值模拟研究分析了竖井贯流泵装置和轴伸贯流泵装置内部三维流动的特点，总结了不同流道形式下水力性能差异产生的原因.周亚军等[3]比较了不同进出水流道的水力特性差异，进而提出了竖井式进水流道和直管式出水流道的优化方案.杨帆等[4]研究不同叶轮转速下贯流泵装置出水流道的流场特征，分析了出水流道进口的平均速度环量随流量系数的变化规律，得出了转速对出水流道水力损失的影响.李四海等[5]研究了竖井贯流泵站机组的振动特性，分析了贯流泵装置在不同工况下不同监测点处的压力脉动特性.ZHU等[6]对导叶叶片固定和可调的贯流泵内部流场进行了数值模拟，提出了一种导叶可调的贯流泵设计方案.

现有研究关于贯流泵水力性能、时均流动、压力脉动等宏观外特性以及平均流动的特点与规律较多，对贯流泵内部复杂三维湍流的流动特性和湍流结构的深入研究较少.文中以某贯流泵为例，采用大涡模拟(LES)和高精度CFD求解器，以湍动能的角度分析圆柱坐标系下该贯流泵设计工况的湍流特性和湍动能各向异性组成结构，拟为深入了解复杂水力机械内部的湍流流动提供依据.

1 数值计算方法

计算所用算法描述不可压缩流体流动的控制方程包括Navier-Stokes方程及连续性方程[7]，即

(1)

∇·u=0,

(2)

式中：u为速度；ρ为密度；p为压强；μ为动力黏度；S为剪切张量；g为重力加速度；在大涡模拟中，亚格子应力τsgs使用动态Smagorinsky模型[8]，而在直接数值模拟中此项为0.

固体在流体中与流体的相互作用通过浸入边界法(Immersed Boundary Method，IBM)来表达.相对于流体本身的流动，流体中存在的固体会表现出额外的力，此时可以视为动量方程右边增加了相应的附加项f，即

ρg-∇·τsgs]+f.

(3)

文中所使用的为基于水平集(level-set，LS)法的IBM[9-10].LS函数为一个距离函数.这里定义为每个网格节点至固体壁面的距离.定义位于流体域中的网格点LS值为正值，位于固体域中的网格点LS值为负值.LS函数值为0的等值面即描述为固体壁面.为了使得速度场在流固交界面处满足无滑移边界条件，固体对流体所施加的额外的力需要准确计算.如图1所示为IBM在二维情况下的示意图.图中红色点为IB点，即为施加力所位于的网格点.通过IB点在其沿壁面法向投影方向获得相应固体边界处的投影点，结合流场中IB点相邻流体点的速度值，可以插值得到固体边界处IB点的速度.

图1 浸入边界法二维示意图

图2所示为壁面函数示意图.为保证施加在IB点的力使得IB点切向速度满足壁面速度规律，考虑到文中算例雷诺数较高，结合文中网格分辨率，采用对数律壁面函数法进行插值计算IB点的壁面切向速度[7]

图2 对数律壁面函数示意图

(4)

式中：UIB，UPP分别为IB点和外部流体点切向速度；uτ为壁面摩擦速度；κ为冯卡门系数；yIB和yPP分别为IB点和外部流体点至壁面的法向距离.而IB点法向速度通过流体点速度与壁面速度进行二次插值获得

(5)

式中：Un,IB和Un,PP分别为IB点和外部流体点法向速度.

文中贯流泵模型由三维软件建模后储存为三角形网格(STL)格式[11]，如图3所示.为构建精确的LS函数场，计算网格节点至固体壁面距离即为计算网格节点至三角形网格上所有三角形单元的最短距离.

图3 以STL离散格式网格表达的贯流泵模型

为区分固体点与流体点上LS函数值的符号，文中使用射线法[12]辨识网格点的位置位于固体边界内外.即从每个网格节点发出一条随机射线，计算其穿过固体边界的次数，如为奇数次，则此点为固体点，LS函数值为负；若为偶数次，此点为流体点，LS函数值为正.

2 研究对象及模拟设置

选择某贯流泵泵段作为研究对象，如图4所示.流动区域包括进水管、前置导叶、叶轮、后置导叶、轴及出水管.前置导叶、叶轮、后置导叶的叶片数分别为5，4，7.运行参数为设计流量10 m3/s, 额定转速250 r/min, 最优工况点扬程 2.2 m.

图4 计算域示意图

建立了基于笛卡尔坐标系(x，y，z)下的三维正交网格，x，y，z方向计算域长度分别为6D，1.1D，1.1D，其中进水管直径D=1.7 m.叶轮中心位于(2.85D×0.55D×0.55D)，轮毂比为0.388，叶顶间隙为0.01D.3个方向的网格数分别为1 024，240，240，共计约5 900万个网格节点.考虑到固体域并不参与N-S方程求解，实际有效的流体计算域网格节点约为4 400万.其中x方向在导叶及叶轮段进行了网格加密，y，z方向使用均匀网格.雷诺数以进口速度和圆管直径作为量纲一化后可得，约为7.15×106.量纲一化后的叶轮转速为11.99.

x方向进口采用均匀入流进口[13-14]，出口采用对流出口.由于包括圆管及叶轮室壁面的流道均由IBM生成，因此y，z方向边界条件并不对流动过程产生影响，本算例简单定义为无滑移边界条件.由于文中泵内流动雷诺数较高，采用基于对数律壁模型的大涡模拟方法进行近壁区的流动求解.

文中计算网格采用交错网格，N-S方程采用投影法求解，空间上使用二阶中心差分格式，时间推进上采用二阶龙格库塔格式.泊松方程在科学计算可移植拓展工具包PETSc上求解完成.程序借助数据并行计算平台信息通过界面(MPI)进行并行计算.计算工作完成于美国圣安东尼瀑布实验室高性能计算中心平台，模拟由768核3.06 GHz的英特尔X5675处理器利用MPI并行计算完成，算例总耗时约15 d.外特性的计算结果表明，设计流量工况下模拟的扬程与效率与试验数据吻合良好，误差分别约为5%和2%，验证了求解器针对文中算例数值模拟方法的可行性与计算精度的有效性.

3 结果分析

图5所示为贯流泵轴向速度(U)和展向速度(V)分布云图.由图可见，在均匀来流的条件下，进水管和前置导叶段流动较为均匀，但当旋转叶轮对水流作用后，尾流表现出了复杂的三维湍流流动特性.显著的平均速度梯度和湍动能的产生与耗散是尾流中水力损失的主要组成部分.随着流动中的对流和扩散，水流瞬时速度波动和不均匀度在流动方向下表现出逐渐下降的趋势.

图5 贯流泵轴向速度和展向速度分布云图

考虑到流体机械的轴对称流动特点，模拟结果转换至圆柱坐标系下进行分析.计算共有24个叶轮旋转周期，取后12个叶轮旋转周期作为后处理数据.叶轮叶片为4片,文中选取48个叶片旋转周期进行数据处理.这里选择对靠近轮毂(r=0.25D)、中等跨度(r=0.35D)、靠近流道内壁(r= 0.45D)3个不同径向距离的圆周截面的展开平面进行分析.湍动能作为湍流强度的度量，对其进行分析可以得到贯流泵湍流流场的强度分布及组成.湍动能公式为

(6)

式中：u′x，u′r，u′θ分别为湍动能k在轴向、径向和圆周方向的脉动速度，〈 〉代表统计下的时间平均物理量.

图6所示为靠近轮毂(r=0.25D)圆周截面的湍动能及其分量分布展开图.图6a为湍动能k的总量分布，图6b—6d为k在流向、径向和圆周方向的组成部分，即为3个方向上的雷诺正应力(下同).由图可见，靠近轮毂位置附近，k的组成主要由流向方向的分量〈u′xu′x〉贡献.因为此处流动本身受流向的对流以及轮毂处壁面的剪切影响，都会使流体产生较大的湍流能量，而这部分能量表现为湍流流向方向的剧烈脉动.而k在径向方向与周向方向的雷诺正应力分量〈u′ru′r〉和〈u′θu′θ〉较小，表现为其脉动强度低于湍动能k在轴向的分量.

图6 r=0.25D圆周截面湍动能分布

图7所示为中等跨度(r= 0.35D)圆周截面的湍动能及其分量分布展开图.此处流动代表了水泵叶轮内主流的湍流特性规律.可以发现，k的组成仍主要由流向方向的分量〈u′xu′x〉贡献.最大值区域出现在叶轮叶片后方的尾流中.此处的尾流是在相对运动的叶轮叶片后方形成，说明叶片绕流后的尾流里包含了较大的湍流能量且主要由流向方向的分量贡献.由于在设计工况下运行，流体质点光滑流过固体叶轮与后置导叶叶片表面，叶轮叶片及后置导叶区域径向及圆周方向脉动强度较低，雷诺正应力〈u′ru′r〉、〈u′θu′θ〉也随之较小；当主流流过后置导叶，圆周方向雷诺正应力〈u′θu′θ〉增加，表明后置导叶对主流圆周方向流动的控制同时增加了圆周方向的湍流脉动强度.

图7 r=0.35D圆周截面湍动能分布

图8所示为靠近流道内壁(r=0.45D)圆周截面的湍动能及其分量分布展开图.由图可见，由于考虑了叶顶间隙的影响，此处的流动较为混乱，湍流具有很大的能量，且3个方向上强度相似.流向的湍动能分量〈u′xu′x〉较大，这是因为间隙处流体在逆压梯度作用下产生回流和间隙泄漏涡，回流、泄漏涡与主流区域的速度梯度很大，在此截面表现出很强的湍流运动.而叶顶处流体的平均圆周速度很大，间隙处叶轮叶片无法继续对流体在圆周方向做功，此处流体的圆周方向转动速度无法与主流速度保持一致.同时受间隙涡的影响，湍动能分量〈u′ru′r〉和〈u′θu′θ〉幅值较大.

图8 r=0.45D圆周截面湍动能分布

图9所示为后置导叶尾流中流向x=1D，2D，3D(叶轮在流向方向的中心为起点)不同位置的湍流能及其组成周向平均后沿半径方向的变化规律，横坐标方向为沿圆管中心线至圆管壁面.由图可见，在后置导叶后方的尾流中，湍动能k在3个方向的幅值大小较为平均.湍动能整体及分量沿主流x方向逐渐减弱，表示湍动能能量的逐步耗散与减弱.湍动能k整体沿径向方向有下降趋势.在x=2D，3D的靠近壁面处，会较快上升再下降为0，说明随着流动中速度混合，尾流内湍动能能量耗散，流动逐渐趋向均匀的圆管流动，且壁面处主要由〈u′xu′x〉贡献.这是因为边界层内流体正应力主要由黏性正应力提供，而湍流正应力很小(壁面处为0).沿壁面从径向向圆管中心线看，由于流动雷诺数较高，边界层以外流体速度增加很快.一般在均匀的圆管流动中，雷诺应力最大会出现在流场中平均速度梯度最大的位置，因此在靠近壁面处发展后的湍流与壁面过渡区域湍动能k会出现最大值.靠近壁面处k主要由〈u′xu′x〉贡献，这是因为此处主流方向平均速度梯度最大，其次为圆周方向.

图9 后置导叶尾流不同位置湍流能及其组成的径向变化规律

4 结 论

1) 湍动能在靠近叶轮叶片根部及中等跨度的圆周面主要由流向的分量提供.而由于后置导叶对圆周方向平均速度的控制，主要引起了湍流在后置导叶段的圆周方向正应力高幅值区域.

2) 叶轮叶片叶顶间隙处，由于间隙回流和泄漏旋涡，在间隙涡区域内流体速度梯度大且流动混乱，湍动能3个方向分量强度均很高.

3) 后置导叶尾流中，湍动能在3个方向幅值相当，且随着中心线向管壁的径向和主流方向湍动能强度均随着距离增大而下降.同时尾流随着流向距离增加，流动中速度逐渐混合而使得湍动能特性趋向均匀的圆管流动.