当发生艏摇、回转或者处于洋流中时，安装在船舶艉部的螺旋桨和舵将不可避免地在斜流中工作，斜流会对螺旋桨和舵的水动力性能产生影响^[1]，斜流同时也会对螺旋桨的尾涡及尾流场产生影响，螺旋桨的尾涡不再对称地沿螺旋桨轴线两端发展，这会导致舵表面原有的受力情况发生改变。张文照等^[2]发现在大攻角的情况下，攻角引起的螺旋桨进速降低，会导致螺旋桨的推力和扭矩随着攻角的增大而增加。文献[3-5]的研究均表明随着角度的增加，螺旋桨的水动力系数不断增加。常欣等^[6]的研究结果表明斜流下桨叶表面的压力分布是不均匀的。螺旋桨的尾涡会引起流场的改变^[7]，并在水动力性能、空泡性能、激振力和噪声扮演者重要角色^[8]。Muscari等^[9]对螺旋桨尾流场中的脉动压力进行了分析，不同的监测位置上脉动压力最大幅值均出现在叶频频率下。频域分析的结果表明，斜流下的脉动压力在整数倍叶频上的波动幅度更小，在非整数倍叶频上的波动幅度更大^[10]。而当螺旋桨尾流场中存在舵等水动力装置时，螺旋桨的梢涡会发生弯曲、断裂和重连^[11]，此时，螺旋桨尾涡中的脉动压力将传递到水动力装置上，导致水动力装置表面同样发生和叶频高度对应的脉动压力^[12-13]。

螺旋桨尾涡对舵结构产生的脉动压力会导致结构振动，严重时会导致共振和结构破坏，因此有必要对这种脉动压力展开研究，本文在前述研究的基础上，使用STAR-CCM+软件平台对斜流与INSEAN E779A螺旋桨尾涡下NACA0020舵表面的脉动压力进行分析，对舵迎流侧和背流侧点的脉动压力进行比较。

1 大漩涡模拟方法

作为一种重要的流动现象，湍流广泛存在于自然界的各种流动中，但是在数值模拟中很难使用直接求解的方式对其进行模拟，因此使用大漩涡模拟(large eddy simulation, LES)方法进行模拟，三维流动的质量守恒方程和动量守恒方程为：

$ \frac{{\partial \rho }}{{\partial t}} + {\rm{div}} (\rho \mathit{\boldsymbol{u}}) = 0 $

(1)

$ \begin{array}{*{20}{c}} {\frac{{\partial (\rho {u_i})}}{{\partial t}} + {\rm{div}} (\rho {u_i}\mathit{\boldsymbol{u}}) = }\\ { - \frac{{\partial p}}{{\partial {x_i}}} + \frac{{\partial {\tau _{ii}}}}{{\partial {x_i}}} + \frac{{\partial {\tau _{ji}}}}{{\partial {x_j}}} + \frac{{\partial {\tau _{ki}}}}{{\partial {x_k}}} + {F_i}} \end{array} $

(2)

式中：ρ表示密度；u表示速度矢量；t表示时间；p表示压力；F_i表示体积力；u_i表示速度矢量分量；τ_ii、τ_ji、τ_ki等项表示粘性应力τ的分量。粘性应力τ可由流体的变形率来表示，故式(2)可写为:

$ \frac{{\partial (\rho {u_i})}}{{\partial t}} + {\rm{div}} (\rho {u_i}\mathit{\boldsymbol{u}}) = - \frac{{\partial p}}{{\partial {x_i}}} + {\rm{div}} (\mu {\rm{grad}} {u_i}) + {S_i} $

(3)

使用LES进行模拟需要使用滤波函数得到大尺度平均分量φ和小尺度分量φ′，该平均分量φ是由空间平均得到的。滤波后的平均分量φ为：

$ \bar \varphi (t,x) = \int_V \varphi (t,{x^\prime })G(x - {x^\prime },\Delta ){\rm{d}}{x^\prime } $

(4)

式中：V为流域；x′为实际流域中的空间坐标；x为滤波后的大尺度空间上的坐标系；G(x-x′, Δ)为滤波函数，其中Δ=(Δ_xΔ_yΔ_z)^1/3。使用滤波函数处理瞬态下的Navier-Stokes方程，LES的控制方程为：

$ \frac{{\partial (\rho \overline {{u_i}} )}}{{\partial t}} + \frac{{\partial \rho (\overline {{u_j}{u_i}} )}}{{\partial {x_j}}} = - \frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}}} \right) - \frac{{\partial ({\tau _{ij}})}}{{\partial {x_j}}} $

(5)

$ \frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}(\rho \overline {{u_i}} ) = 0 $

(6)

亚格子尺度模型使用动态Smagorinsky亚格子尺度模型为：

$ {\tau _{ij}} = 2{\mu _t}\overline {{S_{ij}}} - \frac{2}{3}({\mu _t}\nabla \cdot \overline {{u_i}} + \rho k)I $

(7)

式中：μ_t为亚格子尺度的湍流粘度；k为亚格子尺度湍动能；I为张量积分恒等式^[14-15]。

2 几何模型与收敛性分析

数值模拟使用INSEAN E779 A螺旋桨和NACA0020剖面舵进行，螺旋桨为直径为227 mm的四桨叶螺旋桨，舵的弦长为180 mm。模拟在设计进速J=0.88下进行，为避免舵的梢涡对结果产生影响，舵设置为无限长，桨与舵的间距设置为0.5D，D表示螺旋桨直径。螺旋浆和舵相对关系如图 1所示。

流场域的总长度设置为18.0D，其中螺旋桨前方为5.0D，螺旋桨后方为13.0D，侧面边界距离螺旋桨的距离设置为2.5D以避免边界对流场带来影响。考虑到针对斜流开展模拟，故其中一个侧面和入口面设置为速度入口，与速度入口相对的2个面设置为压力出口，上下2个表面设置为对称平面，在螺旋桨和舵附近设置加密区，加密区完整包裹螺旋桨尾流的流经区域。如图 2所示。

文献[12, 16]指出，使用LES方法对INSEAN E779 A螺旋桨进行数值模拟，0.004 m尺寸的加密网格和1.0×10^-4 s的时间步能够满足精度的要求，但需本模拟使用的0.003 2 m网格尺寸和1.0×10^-4 s时间步长进行验证。针对J为0.8、0.88、1.0时3个进速开展验证，使用Dubbioso等^[17]的试验数据进行对比，结果如表 1所示，表中的K_T和K_Q分别表示推力系数和转矩系数：

$ {K_T} = \frac{T}{{\rho {n^2}{D^4}}},{K_Q} = \frac{Q}{{\rho {n^2}{D^5}}} $

(8)

式中：T、Q分别为螺旋桨的推力和转矩；ρ为密度；n为螺旋桨的转速。

由文献[16]试验结果可知，在J=0.88处该网格和时间步有足够的精度，因此0.003 2 m的网格尺寸和1.0×10^-4 s的时间步长能够用来进行数值模拟，0.003 2 m的网格尺度下部分网格细节如图 3所示。

3 舵表面脉动压力分析

通过对斜流下的水动力性能分析，表3为将其与直航下的水动力性能进行比较结果。表中的推力系数K_T和转矩系数K_Q通过式(8)得到，舵在X方向的阻力系数C_R和螺旋桨效率η为：

$ {C_R} = \frac{{{R_x}}}{{\frac{1}{2}\rho |{V_A}{|^2}S}} $

(9)

$ \eta = \frac{{{K_T}J}}{{{K_Q}2\pi }} $

(10)

式中：V_A表示进速；S为舵面积；R_x为舵在X方向的阻力。

由表 2可知，相对比直航下的情况，斜流下的推力系数K_T、转矩系数K_Q、舵阻力系数C_R以及螺旋桨效率η均增加，且该结果与Hu等^[18]的研究结果接近。

在数值模拟的过程中，使用Q-Criterion方法对螺旋桨的尾涡进行研究，如图 4所示为斜流角α=15°时的螺旋桨尾涡，其中，Q-Criterion的取值为1 000/S²，当螺旋桨尾涡到达舵导边时，螺旋桨的尾涡会发生弯曲，随后沿舵的两侧传播，在舵的背流侧(Y < 0)与舵边界层发生粘连，螺旋桨的尾涡管同时发生偏移；舵的背流侧的轴向速度大于迎流侧(Y>0)的轴向速度，舵背流侧的螺旋桨梢涡较迎流侧弱。螺旋桨尾涡在舵导边上的弯曲如图 5所示。可见，在舵的背流侧螺旋桨尾涡上的弯曲程度并不相同，在舵的迎流面侧，螺旋桨的尾涡向前弯曲，而在舵的背流侧螺旋桨的尾涡向后漂移。

在涡的作用下，舵表面的压强是不对称的，如图 6所示，背流侧的压强最大值出现在舵的随边处，而压强的最小值出现在舵的导边处，舵的迎流侧上压强最大值出现在导边处，但压强最小值出现在舵的中部，同时可见背流侧的压强更低而迎流侧的压强更高，相比于迎流侧，背流侧的压强分布范围较迎流侧更大。说明舵表面上的压强分布不仅在展长、弦长方向上不均匀，在背流面和迎流面上的分布也是不同。

为了对舵表面的脉动压力进行研究，在舵表面设置了一系列的压力监测点，为了对压力监测点进行的区分，在舵表面上建立坐标系，并使用(I, J, K)来对压力检测点进行区分，其中I表示舵弦长方向，J表示舵厚度方向；K表示舵展长方向；I的取值范围为0~5，其中I=0表示压力检测点位于舵的导边上，其余I值表示压力检测点位于导边后的舵表面上；J的取值为1、0和-1。K的取值为1、0和-1，分别表示位于舵上Z=0.113 5 m、Z=0 m、Z=-0.113 5 m位置的点，其中K=1(Z=0.113 5 m)、K=-1(Z=-0.113 5 m)位于螺旋桨叶梢的后方，压力监测点的设置如图 7所示。

对斜流角α=15°下的舵表面脉动压力进行研究，如图 8所示为压力监测点(0, 0, 1)和(0, 0, -1)上监测到的压力与时间的关系，其中t表示时间，P表示压力。由图可知，(0, 0, 1)点和(0, 0, -1)点上监测到的脉动压力均具有强烈的周期性，其波动周期均约为0.01 s，且存在明显的波峰与波谷。

为了进一步研究这种脉动压力的脉动特征，使用快速傅里叶变换对K=1、K=0、K=-1上的脉动压力进行频域分析，频域分析的结果如图 9所示，其中p表示压力脉动幅值，f表示频率。相比于(0, 0, 1)点和(0, 0, -1)点，(0, 0, 0)点的峰值最高，但是此3点在f=100 Hz的频率上均有较高的峰值，该频率对应螺旋桨的叶频，与图 8中的结论一致，这说明舵前缘点上脉动压力的波动频率与螺旋桨叶频具有非常高的相关性。(0, 0, 1)点在所有整数倍叶频的频率上均具有较高的峰值，而(0, 0, 0)点仅在1倍叶频和2倍叶频上具有较高的峰值，(0, 0, -1)在整个频率范围内具有接近的峰值。

由图 9可知，K=1和K=0位置的点监测到的脉动压力在叶频的整数倍上具有明显的波动峰值，最大波动幅值出现在叶频处，波动幅值随着频率的增加而降低，且同一频率下的波动峰值沿弦长方向逐渐降低。在K=0的位置上，在(4, -1, 0)和(5, -1, 0)处观察到了强烈的非整数倍叶频波动。在K=-1上，这种非整数倍叶频波动更加明显，而叶频频率上的波动则被削弱。

为了进一步研究斜流对舵表面脉动压力的影响，对同一展长位置K=1和K=-1的压力监测点监测到的脉动压力的频域特征进行比较，如图 10所示。首先对K=1的点进行分析，在I=1和I=5处，舵的背流面处脉动压力的波动幅度高于舵的迎流面，但是在其他位置则是迎流面J=1上的脉动压力波动幅度大于背流面J=-1；背流面J=-1上的脉动压力在非整数倍叶频上具有更高的波动幅度，在近舵随边的位置，脉动压力在非整数倍叶频上的脉动幅值具有增加的趋势，表明在该斜流角下，由于毂涡完全分布在背流一侧，舵迎流一侧压力脉动完全由梢涡引起，且斜流会抑制舵迎流面上发生的非整数倍叶频波动。

图 10中同时展示了K=-1位置上舵表面脉动压力在舵厚度方向的分布，在舵中间部位，该位置上，舵的迎流面上的脉动压力在叶频频率下具有更高的的幅值，但是在非整数倍叶频和更高的频率上，舵的背流面上的脉动更加剧烈，同时高频脉动压力幅值更高，结合胡健等^[19]的研究，可以发现在舵背流侧发展的毂涡沿舵表面发生了与梢涡偏移方向相反的偏移。

结合文献[12]中的研究，这种出现在叶频及整数倍叶频下的脉动压力是螺旋桨梢涡涡引起的，这些脉动压力的强弱与螺旋桨梢涡高度相关，透过脉动压力的分布，可以观察到螺旋桨-舵系统重斜流对螺旋桨尾涡产生的影响，图 11同时给出了为Z=0平面上的垂向涡量分布，涡量值的正负表示方向。

由图 11可知，在舵迎流面一侧流场中的垂向涡量增强，在舵和和斜流的共同作用下，舵背流面一侧流场中的螺旋桨尾涡明显减弱，导致舵背流面一侧表面脉动压力的波动幅值在整数倍叶频和非整数倍叶频上接近。

4 结论

1) 斜流和螺旋桨尾涡共同作用下舵迎流面上脉动压力均与螺旋桨的叶频特征更加明显，脉动压力在整数倍叶频上具有明显的波动幅值。

2) 舵的背流面侧的脉动压力与叶频的相关性较迎流面侧弱，背流面上的毂涡也发生了偏移。

本文对斜流和螺旋桨尾涡共同作用下的舵表面脉动压力进行了研究，但是未考虑动态变化的斜流角，下一步的研究应当进行考虑。