基于OpenFOAM的螺旋桨紧急倒车工况数值模拟

郭海鹏, 邹早建, 李广年

(1. 宁波大学 海运学院, 浙江 宁波 315211; 2. 上海交通大学 船舶海洋与建筑工程学院, 上海 200240)

船舶的紧急制动能力与船舶的航行安全性直接相关，而紧急制动能力的优劣很大程度上取决于紧急倒车工况下螺旋桨的水动力特性.相对于目前研究较多的螺旋桨正车前进工况，螺旋桨紧急倒车工况下的绕流场更为复杂，具有明显的非定常及非线性特性，这导致作用在螺旋桨上的水动力载荷同样呈现出复杂的变化规律.因此，从提高船舶航行安全性的角度出发，对螺旋桨紧急倒车工况下的水动力特性及绕流机理的深入研究能够为全面评估船舶的紧急制动能力提供依据.

针对螺旋桨紧急倒车工况，早期开展的相关研究主要基于实验流体动力学(EFD)技术[1-2].近年来，发展迅速的计算流体动力学(CFD)技术成为开展螺旋桨水动力性能研究的有力工具.Chen等[3]基于自主开发求解器对DTMB4381螺旋桨的四象限工况，即紧急正车、正车前进、紧急倒车及倒车后退进行了数值模拟，分析了不同工况下的叶片压力、流线分布及尾流场特征；李理等[4]基于Fluent对AU5-65型螺旋桨的四象限水动力性能进行了数值预报，分析了推力系数的变化规律及相应的流场特征；Verma等[5]和Kumar等[6]基于自主开发求解器对DTMB4381螺旋桨的紧急倒车工况进行了数值模拟，探讨了上游船体对螺旋桨的影响、螺旋桨侧向力与叶片附近流动分离间的联系；王贵彪等[7]基于Fluent对某导管桨的四象限工况进行了数值模拟，分析了不同工况下的推力、转矩及周围流场特征；张文璨等[8]基于Fluent对MARIN B系列螺旋桨的四象限工况水动力性能进行了数值预报，并将数值结果与试验图谱进行了对比验证；Pergande等[9]基于ANSYS CFX对某调距桨的紧急倒车工况进行了数值模拟，探讨了该工况下的流动分离现象；Pontarelli等[10]基于ANSYS REX对DTMB4381和MARIN 7371R螺旋桨的紧急倒车工况进行了数值模拟，对不同进速条件下的流场特征进行了分析；Chen等[11]基于Fluent对DTMB4381螺旋桨紧急正车、紧急倒车工况进行了数值模拟，分析了不同工况下的桨叶压力及流线分布等流场特征；王超等[12-13]基于STAR-CCM+对E1619螺旋桨和SUBOFF潜艇模型的紧急倒车工况进行了数值模拟，分析了桨叶诱导环状涡的产生机理、发展过程以及上游艇体产生的影响；樊翔等[14]基于STAR-CCM+对KCS船模的倒车工况进行了数值模拟，分析了该工况下的螺旋桨载荷变化特性；Wang等[15]基于自主开发的求解器naoe-FOAM-SJTU对KCS船模的紧急停船运动进行了数值模拟，分析了该过程中船体与螺旋桨之间的相互干扰特性.

综上，针对螺旋桨紧急倒车工况的相关数值研究大多借助商业CFD软件实现，这类软件存在价格昂贵、代码不可见等局限性.相较而言，开源CFD平台OpenFOAM以其开源免费、支持高度自定义等优势，得到了越来越多研究人员的关注.本文基于OpenFOAM对5叶螺旋桨DTMB4381的正车前进和紧急倒车工况进行数值模拟，通过与国际公开的模型试验结果进行对比，考察所采用数值方法的有效性，并基于数值模拟结果探讨了螺旋桨紧急倒车工况下局部绕流场特征随进速的变化规律及其与螺旋桨整体水动力特性之间的联系.

1 求解方法

1.1 控制方程

所考虑的螺旋桨绕流为黏性、不可压缩流体流动，其控制方程为雷诺平均的连续性方程和动量方程，即雷诺平均纳维斯托克斯(RANS)方程：

(1)

(2)

式中：ρ为流体密度；U为平均速度；p为平均压力；μeff=ρ(ν+νt)为有效动力黏度，其中ν和νt为运动黏度和湍流涡黏度，后者通过湍流模型确定.

1.2 湍流模型

湍流涡黏度采用Menter等[16]提出的两方程SSTk-ω模型进行求解，其中k和ω为湍流动能和比耗散率.该模型通过剪切应力输运公式将k-ω模型与标准k-ε模型结合，在近壁面采用k-ω模型，而在远场采用标准k-ε模型，充分发挥了两者的优势，是目前船舶与海洋工程领域应用较为广泛的湍流模型.在OpenFOAM中，SSTk-ω模型的输运方程[17]如下：

(3)

(4)

(5)

(6)

式中：S为应变率的不变量；β*、αk、γ、β、αω及αω2为方程常量；F1为结合k-ω模型和k-ε模型的混合函数.变量的具体定义及常数的取值见文献[16-17].

1.3 速度压力耦合求解

控制方程求解过程中的速度与压力的耦合通过OpenFOAM中的PIMPLE算法实现.该方法将Patankar等[18]在1972年提出的压力耦合方程组的半隐式方法(SIMPLE算法)与Issa[19]在1986年提出的PISO算法相结合，可以实现瞬态问题的大尺度时间步求解以提高计算效率.

1.4 物体运动处理

螺旋桨旋转运动模拟采用OpenFOAM中的多重参考系法(MRF)及任意网格界面法(AMI)来实现[20]：MRF方法通过在不同区域内建立不同的参考系与控制方程将螺旋桨动态旋转转化为静态问题进行求解；AMI方法通过在不同区域交界面插值进行通量传递，从而实现螺旋桨动态旋转运动的直接模拟.

2 数值计算

2.1 研究对象及工况

以5叶螺旋桨DTMB4381模型为研究对象，其直径为0.305 m，螺距比为1.2，侧斜角为0°，盘面比为0.725，毂径比为0.2，设计进速系数为0.889，螺旋桨为右旋，其几何模型如图1所示.首先对螺旋桨正车前进工况进行数值模拟，然后进行螺旋桨紧急倒车工况数值模拟，倒车工况的进速方向与正车前进工况相同，但螺旋桨反向旋转.详细模拟工况如表1所示.其中，J为进速系数，UA为进速，n为螺旋桨转速.

图1 DTMB4381螺旋桨的后视图及侧视图Fig.1 Back view and side view of DTMB4381 propeller

表1 数值模拟工况Tab.1 Numerical simulation condition

计算结果包括不同进速系数下的推力系数和转矩系数，定义如下：

(7)

(8)

(9)

式中：D为螺旋桨直径；T为螺旋桨推力；Q为螺旋桨转矩.

2.2 计算域及网格划分

计算域及边界条件设置如图2所示.根据螺旋桨流场特性选取圆柱形计算域，其直径为6D，长度为8D，其中桨盘面前、后分别为3D、5D，同时选取包含螺旋桨的圆柱体作为旋转域，其直径为1.5D，长度为2D，其中桨盘面前、后分别为1D.入口边界设为fixed value边界条件；出口边界设为inlet outlet边界条件，根据速度方向在fixed value与zeroGradient之间进行切换；外侧边界设为no slip边界条件；螺旋桨设为moving wall velocity边界条件；旋转域与静止域的交界面设为cyclic AMI边界条件，可以实现不同区域间的流场信息传递.

图2 计算域及边界条件设置Fig.2 Setup of computational domain and boundary condition

计算域采用OpenFOAM自带的前处理工具进行空间离散，网格划分如图3所示.首先，采用blockMesh工具将计算域离散为均匀的六面体背景网格，然后采用snappyHexMesh工具进行螺旋桨附近区域的网格细化、表面对齐等操作.考虑到螺旋桨绕流场的特性，在螺旋桨附近及下游区域进行体加密，以便捕捉螺旋桨尾流场特征.为了更好地求解螺旋桨近壁面流动，在螺旋桨壁面附近生成4层棱柱形网格，无因次壁面距离y+值控制在30左右.

图3 计算域网格划分Fig.3 Grid generation of computational domain

2.3 迭代计算设置

为了提高数值计算效率及稳定性，在计算过程中首先应用simpleFoam求解器开展定常计算，螺旋桨旋转采用MRF方法模拟，以此获得稳定的初始流场.然后，以获得的初始流场结果为初始值，应用pimpleFoam求解器开展非定常计算，螺旋桨旋转采用AMI方法模拟.在非定常计算过程中，时间步长设为1.0×10-4s，模拟物理时间为5.0 s，对应螺旋桨旋转50圈所需时间.

2.4 网格收敛性分析

图4 桨叶网格划分Fig.4 Grid generation on blade

网格收敛性分析结果见表2，其中KT和KQ的计算值通过对螺旋桨旋转最后5圈的数据进行平均获得.从表中可以看出，KT和KQ的网格收敛率RG为负值且绝对值小于1，因而表现为振荡收敛.同时，网格不确定性无量纲值UG%S1<1，表明当前计算结果受网格尺度的影响较小.此外，不同尺度网格计算对比结果表明，细网格计算结果与中网格计算结果相差很小.因此，权衡数值计算的精度和效率，在后续数值研究中均采用中网格进行计算.

表2 网格收敛性分析结果Tab.2 Results of grid convergence study

3 计算结果与分析

3.1 水动力载荷

螺旋桨正车前进工况下的推力、转矩以及效率的计算结果如图5所示.通过与国际上公开的模型试验结果[22]进行对比发现，本文计算结果与模型试验值吻合良好，推力误差值基本在3%以内，在高进速(J=1.1)工况下误差值最大达到14%；转矩误差整体相对略高，在低进速(J=0.1)工况下误差值最大约为7%；随着进速增大，误差值基本在5%以内.总体而言，数值计算精度与常规的商业CFD计算软件相当[23]，验证了所采用的计算平台及数值方法在预报螺旋桨正车前进工况下的水动力性能方面的有效性.

图5 螺旋桨正车前进工况模拟结果Fig.5 Simulation results of propeller ahead condition

螺旋桨紧急倒车工况下推力、转矩的计算结果如图6所示.通过与国际上公开的模型试验结果[22]进行对比发现，数值计算结果与模型试验结果相比存在一定偏差，且随着进速系数的改变发生显著变化.在低进速(J=-0.5)工况下，推力及转矩的计算值偏高，误差值均在15%左右；在中进速(J=-0.7)工况下，推力及转矩的计算值略微偏低，误差值分别为8%和2%；在高进速(J=-0.9)工况下，推力及转矩的计算值明显偏低，误差值分别达到20%和15%.根据前文网格收敛性分析结果可知，当前螺旋桨紧急倒车工况下的水动力性能预报精度受网格尺度影响较小，而应该与所采用的湍流建模方法难以准确捕捉该工况下存在的大规模流动分离现象有关.总体而言，尽管本文所采用的数值方法对螺旋桨紧急倒车工况下的水动力预报精度相对有限，但数值计算结果仍呈现出与模型试验结果相一致的总体趋势.

图6 螺旋桨紧急倒车工况模拟结果Fig.6 Simulation results of propeller crashback condition

3.2 绕流场特征

3.2.1泄出涡系结构 不同进速系数下数值模拟得到的紧急倒车工况下螺旋桨绕流场中的涡系结构如图7所示.涡系结构采用涡量等值面Q=1 000 表示，同时采用无因次轴向速度Ux/U0着色，其中U0为远场入流速度.结果表明，紧急倒车工况下螺旋桨梢部产生的泄出涡十分明显，并呈现出与正车前进工况截然不同的结构.由于紧急倒车工况下螺旋桨反向旋转，螺旋桨产生的排出流与来流方向相反.在来流影响下，桨叶梢部泄出涡有沿螺旋桨径向发展的趋势，同时，螺旋桨叶背根部也产生了明显的泄出涡.由涡量等值面上的速度分布可知，在靠近叶梢部位，螺旋桨诱导速度的轴向分量与来流方向相反，在叶梢外部区域，来流的影响逐渐占据主导地位，呈现出与来流相一致的流向.

通过对比不同进速下的涡量等值面图可以发现，涡系结构的尺度和位置随入流速度的改变而发生变化.在低进速(J=-0.5)工况下，螺旋桨梢部泄出涡沿周向均匀分布且整体向螺旋桨上游发展.随着进速的增加，螺旋桨叶梢部泄出涡逐步增强并呈现出周向非均匀特性，而且逐渐呈现出向螺旋桨下游发展的趋势.与此同时，叶背根部的泄出涡也愈发显著，叶面根部也有新的泄出涡生成.此外，由于来流速度影响的增强，螺旋桨梢涡上的轴向速度分布也发生了明显变化，来流的影响逐渐在螺旋桨绕流场中占据主导地位.

3.2.2剖面流线分布 不同进速系数下数值模拟得到的紧急倒车工况下螺旋桨附近流场剖面的流线分布如图8所示.此外，图中还给出了剖面上的无因次轴向速度Ux/U0分布.剖面流线分布情况清晰地反映了螺旋桨紧急倒车工况下涡系结构的产生机理及特性：在紧急倒车工况下，螺旋桨反转所诱导产生的流动与来流方向相反，在螺旋桨上游流场产生挤压、剪切，形成沿螺旋桨径向发展的汇聚流，结果与图7中的桨叶梢部泄出涡相对应；同时，由于螺旋桨反转所产生的流动影响，来流在螺旋桨下游发生分离而形成分离流.在汇聚流与分离流的共同作用下，桨盘面后方形成了大尺度环状涡结构.

图7 紧急倒车工况下Q=1 000等值面Fig.7 Iso-surface of Q=1 000 under crashback condition

通过对比不同进速下的流线图可知，由汇聚流和分离流形成的环状涡结构的尺度和位置随入流速度的改变而发生变化.低进速(J=-0.5)工况下螺旋桨诱导速度的影响占据主导地位，汇聚流和分离流的作用范围较广，形成的环状涡结构尺度也较大.随着进速的增加，来流的影响逐渐增大，汇聚流和分离流的作用区域逐渐减小，所形成的大尺度环状涡逐渐收缩，并呈现出向螺旋桨下游方向移动的趋势.

3.2.3叶片压力分布 数值模拟得到的紧急倒车工况下的螺旋桨叶背及叶面压力分布如图9和图10所示.Cp表示无因次化的压力系数，计算公式为

(10)

式中：Δp为相对压力.紧急倒车工况下由于螺旋桨反向旋转，导边与随边也相应改变.螺旋桨叶背以正压为主，表现为压力面，最大正压出现在导边附近区域；相反地，螺旋桨叶面以负压为主，表现为吸力面，最大负压出现在导边及叶梢附近区域.由于压力面与吸力面的压差作用，螺旋桨紧急倒车工况下产生反向的推力.

不同进速下的叶片压力云图对比结果表明，叶片压力分布随入流速度的改变而发生变化.在低进速(J=-0.5)工况下，叶片上的高压区及低压区主要位于导边附近区域，且叶片压力呈现出较好的周向均匀性.随着进速的增大，压力面的高压区及吸力面的低压区均有所增强，特别是压力面的高压区逐渐发展到整个叶背大部分区域.值得注意的是，在压力面靠近导边的叶根部分区域呈现出相对较低的压力，这与高进速工况下的叶背根部位置的泄出涡有关.类似地，吸力面靠近导边的叶根部位形成了新的低压区，同样与叶面根部的泄出涡有关.此外，随着进速的增加，叶片上的压力分布还呈现出一定的周向非均匀性，这与桨叶泄出涡的结构特征相对应.

图9 紧急倒车工况下螺旋桨叶背压力分布Fig.9 Pressure distribution on the back of propeller under crashback condition

图10 紧急倒车工况下螺旋桨叶面压力分布Fig.10 Pressure distribution on the face of propeller under crashback condition

4 结论

基于开源CFD平台OpenFOAM对螺旋桨DTMB4381模型的正车前进及紧急倒车工况进行数值模拟，通过与模型试验数据的对比以及对数值计算结果的分析，可以得到以下结论：

(1) 对于螺旋桨正车前进工况，本文数值模拟结果与模型试验数据吻合良好，除高进速(J=1.1)工况外，推力和转矩的误差均在5%以内；对于螺旋桨紧急倒车工况，数值模拟结果与模型试验数据相差较大，推力和转矩在低进速(J=-0.5)和中进速(J=-0.7)工况下误差在10%左右，在高进速(J=-0.9)工况下误差达到20%.

(2) 紧急倒车工况下，螺旋桨梢部泄出涡在来流影响下呈现沿螺旋桨径向发展的趋势，同时叶片根部也有泄出涡产生.随着来流速度的增大，桨叶泄出涡有所增强并向螺旋桨下游发展，此外还呈现出一定的周向非均匀性.

(3) 紧急倒车工况下，由于螺旋桨的诱导速度与来流相互作用，形成了螺旋桨上游的汇聚流及下游的分离流，两者共同作用在桨盘面后方发展成为大尺度环状涡.该环状涡随着来流速度的增加而逐渐收缩，并呈现出向螺旋桨下游方向移动的趋势.

(4) 紧急倒车工况下，螺旋桨叶背为压力面而叶面为吸力面，压力面的高压区和吸力面的低压区主要分布在导边附近.随着来流速度的增加，桨叶吸力面上的低压区及压力面上的高压区均有所增强，同时压力分布呈现出与桨叶泄出涡相类似的周向非均匀性.

由以上研究结果可知，螺旋桨紧急倒车工况下的流场特征较为复杂，涉及大规模流动分离及回流现象，而本文所采用的RANS方法在该方面求解精度相对有限.因此，需要针对这一问题开展更为深入的研究，将采用具有更高精度的数值方法如大涡模拟(LES)方法、分离涡模拟(DES)方法等进行湍流模拟，从而更准确地揭示螺旋桨紧急倒车工况下的水动力特性及相应的流动机理.