基于定常RANS方法的螺旋桨毂涡结构分析

李鹏程，刘登成，董郑庆，郭峰山，夏灏超，陈雷强

(1. 中船重工（上海）节能技术发展有限公司，上海 200011；2. 中国船舶科学研究中心，江苏 无锡 214082)

0 引 言

螺旋桨尾涡系取决于螺旋桨负荷分布形式，而螺旋桨根部负荷对于其推进效率的影响一直是螺旋桨设计研究的重要内容之一。Lerbs[1]最早提出求解最佳径向环量分布的变分法时认为叶根处某一桨叶压力面与相邻桨叶吸力面发生压力中和，造成叶根处升力为0。Kerwin[2]采用镜像涡的办法近似处理桨毂，并认为螺旋桨根部存在环量。孙文愈[3]将面元法计算的桨毂诱导速度加入变分法求解最佳环量中，得到的最佳环量分布在叶根处存在负荷。Brizzolara[4]的计算结果也得到了相同的结论。Hong[5]和Shin[6]采用CFD数值模拟技术分析了螺旋桨环量分布形式，结果均显示桨叶根部存在负荷。Jessup[7]应用LDV流场测量技术分析了该问题，证实了螺旋桨根部的环量是真实存在的，并通过与Wang[8]的计算结果对比，阐述了桨毂对于螺旋桨根部环量存在重要影响。

直接针对螺旋桨毂涡结构进行研究的论文并不多见。Kumar[9]采用大涡模拟（LES）方法进行螺旋桨尾流场不稳定性研究时简要分析了毂涡的流动特性。Felli[10-11]通过高速摄影进行了梢涡和毂涡空泡流动显示观测对螺旋桨尾涡的不稳定性进行了细致分析，认为螺旋桨的梢涡和毂涡的相互干扰是尾涡不稳定的重要原因。此外，回收利用螺旋桨毂涡能量的节能装置研究也是毂涡结构研究的一个重要延伸。Kim[12]，李鑫[13]，高德宝[14]等学者对于毂帽鳍、舵球、消涡轮等不同种类的毂涡节能装置进行过分析研究。

本文基于粘流CFD技术进行了螺旋桨根部涡系的研究分析，并将数值模拟结果与螺旋桨水动力测试及LDV流场测量结果进行了对比，验证了其可靠性。通过将计算结果与空泡流动显示进行对比，描述了螺旋桨毂涡系结构形态及沿轴向的发展变化，揭示了毂涡的形成机理，这对于螺旋桨设计中桨毂影响分析及毂涡节能装置设计具有一定的参考意义。

1 CFD模拟方法

1.1 数值方法

本文螺旋桨水动力性能及尾流场的数值模拟以连续性方程和三维不可压缩雷诺应力（RANS）方程为控制方程，并以SSTk-ω二方程湍流模型确保方程封闭，方程离散采用有限体积法进行。对流项和耗散项的离散分别应用2阶迎风格式和2阶中心差分格式。求解过程基于Simple算法，离散方程采用Gauss-Seidel进行迭代求解。

1.2 计算模型和边界条件

计算模型为5叶螺旋桨，半径R为125 mm，三维模型如图1所示，螺旋桨前方为延伸至计算域入口的与桨毂大端等直径的圆柱。数值模拟计算域为圆柱形区域，入口位于螺旋桨盘面上游4R，出口位于螺旋桨盘面下游16R，计算域半径为6R。本文直角坐标系x沿螺旋桨轴线向下游为正，y，z方向符合右手法则，螺旋桨参考线与y轴正向重合；柱坐标系轴向a向下游为正，径向r向外为正，切向t逆时针为正，参考线位置为0相位。螺旋桨旋转方向为右旋，旋向为负。

图 1 螺旋桨三维模型Fig. 1 Propeller 3D model

计算域入口V0为3 m/s的等值速度进口边界，出口为压力出口边界，螺旋桨、桨毂及桨轴表面为无滑移边界。进速系数J=V0/(nD)=0.55，n为螺旋桨转速。

1.3 计算网格

数值模拟计算网格总数约为18 M，桨叶区域采用非结构化网格，叶根和桨毂壁面网格为0.002 R；桨叶以外区域采用六面体结构化网格，网格增长率为1.05。桨叶根部和桨毂表面的网格形式如图2所示。

图 2 桨叶根部及桨毂表面网格形式Fig. 2 Mesh grids on propeller root and hub

本文所列结果速度以来流速度V0，长度以螺旋桨半径R，压力以ρV02进行无量纲化处理。

2 数值模拟验证

2.1 螺旋桨水动力

该螺旋桨在中国船舶科学研究中心拖曳水池和空泡水筒中分别进行过水动力测试[15]，其结果与本文CFD计算结果的对比见表1，其中螺旋桨推力系数KT、扭矩系数KQ表达于式（1），雷诺数Re见式（2）。

式中：T为螺旋桨推力；Q为螺旋桨扭矩；D为螺旋桨直径；C0.75R为螺旋桨r/R=0.75剖面弦长。

表 1 J=0.55时螺旋桨推力系数和扭矩系数CFD与模型试验结果对比Tab. 1 Comparisons of propeller thrust and torque by CFD simulation and model test at J=0.55

2.2 螺旋桨尾流场

该螺旋桨曾经应用LDV流动测试技术进行过尾流场速度测量[16]，本文CFD数值模拟得到的螺旋桨尾流场与之进行对比。在轴向x/R=0.23和0.54两个盘面位置的周向平均速度沿径向分布的对比结果如图3所示，在这2个轴向位置，径向位置r/R=0.36的流场周向速度分布对比结果如图4所示。从图3可以看出，LDV与CFD得到的周向平均速度差异在大部分区域内吻合良好。从图4可以看出，CFD模拟较准确地描绘了桨叶尾流区速度峰值的位置和大小。图5为螺旋桨毂涡区域内的流场周向速度分布对比结果。螺旋桨毂涡区域湍动能较强[9]，流场脉动较大，这从图5中波动较大的LDV测试结果也可以看出。本文CFD模拟得到的速度场在毂涡区域与LDV测试结果吻合程度相对毂涡外的螺旋桨尾流区为低，但是仍可以得到正确的周向平均切向速度峰值出现的径向位置（见图3（b）），周向速度分布的趋势也与LDV结果相符。

图 3 螺旋桨尾流场CFD模拟与LDV测量周向平均速度结果对比Fig. 3 Comparisons of circumferential averaged propeller wake by CFD simulation and LDV measurement

因此，采用本文的数值模拟结果进行螺旋桨尾涡系分析是可行的。

3 结果及分析

3.1 桨叶环量分布和尾涡强度

螺旋桨的环量分布采用斯托克斯定理，通过对桨后流场中的切向速度进行周向积分得到，其计算方法见式（3），图6为计算结果。

图 4 螺旋桨尾流场CFD模拟与LDV测量周向速度分布结果对比（r/R=0.36）Fig. 4 Comparisons of circumferential averaged propeller wake by CFD simulation and LDV measurement（r/R=0.36）

图 5 毂涡区内CFD模拟与LDV测量周向速度分布结果对比（r/R=0.096，x/R=0.54）Fig. 5 Comparisons of circumferential distribution propeller wake by CFD simulation and LDV measurement in hub vortex region（r/R=0.36）

式中，Z为螺旋桨叶数，Γ为桨叶环量。

Hong[5]通过CFD数值模拟分别采用该方法和绕叶片直接速度积分进行径向环量分布计算，除叶梢区域外，2种方法得到的结果吻合很好。董郑庆[17]和Jessup[18]采用类似方法分析了桨叶径向环量分布。

图 6 螺旋桨径向环量分布（x/R=0.16）Fig. 6 Radial circulation distribution

图 7 螺旋桨内半径尾流场轴向涡量Ωa分布（白色圆圈r/R=0.33）Fig. 7 Inner radius propeller wake axial vorticity distribution(white circle r/R=0.33)

涡识别的方法有很多[19-20]，本文采用传统的涡量方法。图7为螺旋桨后方轴向位置x/R=0.23和0.45处的轴向涡量分布。除在叶根处外，在x/R=0.33附近同样出现了较强的轴向涡量，这一特征在x/R=0.45位置更加明显，通过对图6桨叶环量分布求导发现，环量梯度最大值确实出现在这个半径位置附近，结果如图8所示。

3.2 螺旋桨根部涡系形态

图9为螺旋桨表面压力分布及尾流场中涡线。图10为中纵剖面流场视图，图10（a）为轴向涡量分布和涡线图，图10（b）为压力分布图。图11为空泡流动显示试验照片，图11（a）、图11（c）为螺旋桨相同相位下不同时刻的照片，图11（a）～图11（f）分别为图11（a）～图11（c）的局部放大视图。

图 8 螺旋桨径向环量分布梯度（x/R=0.16）Fig. 8 Radial circulation gradient distribution

图 9 螺旋桨表面压力分布及尾流场涡线图Fig. 9 Pressure distribution on propeller and vortex line in wake

图9 中位于桨毂后方轴心附近的涡线1已经基本接近直线，终止于桨毂小端壁面中心位置；涡线2和涡线3有明显的螺旋形结构，涡线2在靠近桨毂小端壁面呈扩散状并终止于该壁面，而涡线3绕过小端面与叶根尾涡相连结。螺旋桨毂涡在向下游发展的过程中直径不断增大，涡线3的直径一直大于涡线2，且一直在r/R=0.1附近以外区域向下游发展。图10（a）中亦可以看出这一明显特征，即r/R=0.1附近以内区域的涡线均终止于小端面，而绕过桨毂小端面与叶根涡相连结的涡线并未进入该区域。

图 10 中纵剖面（Z=0）涡量和压力分布（白色区域为桨毂）Fig. 10 Middle longitudinal section (Z=0) vorticity and pressure distribution

图 11 相同相位不同时刻流动显示照片Fig. 11 Pictures of hub vortex cavitation

图10 （a）中，靠近桨毂端面位置涡线的扩张的同时，等涡量线的直径有所增大，图11中毂涡空泡也出现了在靠近桨毂壁面位置直径有所增大的现象，这与涡管终止于壁面时的扩张形态非常相似。图10（a）中涡量场在向下游发展至x/R= 0.6附近位置时，叶根部尾涡相连结的涡线向内收缩至y/R=±0.1附近，相同区域毂涡的轴向涡量出现一个明显增大的过程，这与图11（f）中标线1右侧毂涡空泡直径明显变大的变化趋势一致。同时，图10（b）压力云图中也可以看到出现了相对应的变化过程。这说明流动发展至x/R=0.6的轴向位置，桨叶根部尾涡与中心涡出现相互作用，使轴向涡量增大，压力降低，空泡直径增加。图11（f）中，标线1右侧空泡表面出现了不规则的白色螺旋状结构，这一形态与图9涡线结构非常相似，涡线3在涡线1和涡线2的表面缠绕发展，这进一步说明了毂涡的结构形式。在图11（d）、图11（e）中也看到了这一现象，只是螺旋结构的螺距存在一定差异。

3.3 涡通量轴向变化

为进一步研究螺旋桨毂涡结构，对不同轴向断面分别进行轴向涡通量Qa计算，计算方法见式（4），在桨毂区域（x/R＜0.5）径向积分范围为桨毂边界层外至半径R。图12为计算结果。

从图12可以看出，由于涡量耗散，在桨毂端面之前的轴向涡通量呈现不断缓慢下降的趋势，而到达桨毂端面位置时，轴向涡通量突然降至0。由于在r/R=1.0的半径位置的环量已经降至0（见图6），由斯托克斯定理（公式（5））易得，在桨毂区域轴向涡通量大小等于桨毂环量，正负相反，而在桨毂后方（x/R＞0.5）涡通量会降为0。对于根部存在负荷的螺旋桨，在桨毂区域，桨毂环量与尾流场涡通量相互抵消，而在离开桨毂以后，从桨毂端面中心位置产生了一根轴心涡与外围轴向涡通量总和为0，以保证涡通量守恒。进而，针对3.2中已经得到的桨叶根部尾涡仅收缩至r/R=0.1附近区域的结果，计算在x/R=0.51，0.54，0.70断面位置桨毂半径范围内的轴向涡通量，其值与桨毂环量相等时的半径r/R位置分别为0.15，0.13，0.09，这描绘了叶根尾涡收缩进入桨毂区域的轨迹，与图10（a）的涡线轨迹非常接近，进一步验证了毂涡是由轴心涡和缠绕在其表面的叶根尾涡共同构成的这一结论。

图 12 轴向涡通量沿轴向位置的变化Fig. 12 Vortex flux variation along axial position

4 结 语

本文采用CFD技术对螺旋桨尾流场进行数值模拟，并与水动力测试与LDV流场测量结果进行了对比，验证了计算结果的可靠性。通过数值模拟与流动显示结果的对比分析，说明了毂涡系的结构形态，论述了毂涡形成机制。

结果表明，螺旋桨毂涡由从桨毂小端面形成的轴心涡及缠绕在其表面的叶根尾涡共同组成；轴心涡强与桨毂环量相等，与轴心涡以外区域的轴向涡通量抵消，保证了涡通量守恒。

从流动显示结果看，本文的数值模拟正确展示了螺旋桨毂涡系结构，但本文计算仅为定常流动模拟，并未针对后端面的分离流动影响详细分析，这有待将来开展进一步的研究工作。