噪声特性是船舶螺旋桨设计中必须考虑的问题. 噪声的数值预报技术目前已经比较成熟. 最早可以追溯到基于连续性方程和动量方程推导得到的Lighthill方程^[1]，它适用于湍流声源在自由域的传播问题. Curle^[2]基于Lighthill方程进行发展，得到适用于固壁面辐射噪声的声类比方程. Williams等^[3]基于广义格林函数重新推导波动方程，综合考虑壁面运动和流场湍流脉动的贡献，方程右端物理意义清晰，使用最广泛.

混合方法是结合计算流体力学（computational fluid dynamic, CFD）技术和声学边界元理论计算噪声声压级的方法，求解得到计算域以及固壁面上的压力速度分布后，结合声类比理论对监测点的声压进行推算. 混合方法是目前进行声学模拟最广泛使用的方法. Gorji等^[4-7]利用混合法计算螺旋桨非空化噪声声压级，分析低频线谱噪声的特征，研究声源指向性. Sun等^[8-10]采用混合法对螺旋桨推进装置进行噪声声压级数值验证，开展优化降噪分析.

船体辐射噪声主要来自于以下几个部分^[11]：机械振动、螺旋桨的旋转和流动中形成的湍流以及螺旋桨上剧烈的空化作用. 其中由螺旋桨旋转和湍流作用产生的噪声统称为非空化噪声，它可以按照作用形式分为螺旋桨的排水噪声，叶片两侧压差产生的噪声以及叶片附近湍流产生的噪声. 它们分别对应单极子、偶极子和四极子噪声. 对于船舶螺旋桨，湍流产生的四极子噪声声压级相对于单极子和偶极子噪声声压级的贡献较小，通常可以忽略. 综上所述，螺旋桨在非空化状态下，噪声主要来源于偶极子. 这和螺旋桨的表面脉动力有关联，因此叶片附近流动的模拟极大地影响噪声计算的准确性.

Suboff是美国国防高级研究计划署（DARPA）设计的潜艇标模，Liu等^[12]给出实验测得的suboff的阻力和艇表面的压力分布等特性. 相关噪声的实验，特别是在艇与螺旋桨组合体的噪声测量方面，公开发表的文献较少. 已有研究大多集中在简单构型体或者单桨的敞水实验. Felice等^[13-14]对标模桨开展实验和数值方面的研究. 曾赛等^[15-16]对转桨和单桨的噪声进行实验和数值方面的研究.

本文以DARPA suboff模型和DTMB 4119型螺旋桨作为研究对象，借助STAR-CCM+12.06平台计算水动力学性能和非空化噪声声压级. 比较不同湍流模型和网格划分方案所得模拟结果的差异，通过宏观水动力学性能及局部流场参数，评估不同湍流模型的准确性. 采用FW-H方程进行声学计算，分析噪声的频谱特征及指向性特征，研究艇尾迹伴流场对噪声特性的影响规律.

DTMB 4119桨叶模型的直径为0.304 8 m，参照文献[13,17,18]的计算域尺寸，圆柱形计算域直径D₁=1.6 m，旋转区域D₂=0.32 m，长度为0.30 m，在旋转区域指定运动，转速为10 r/s. 距离入口的长度L₁=0.75 m，距离出口的长度L₂=2 m，单桨模型速度入口流速为2.539 m/s，出口为压力出口条件，侧边边界设置为对称平面.

Suboff潜艇模型全长L=4.356 m，艇身最大直径记为D_max=0.508 m，参照文献^[19]的计算域布置，计算域布置与图1（b）一致，将suboff模型放置在圆柱形流场中央，直径D=10D_max=5.08 m，艇上游长度L₃取3.5 m，接近一倍艇长；流场下游长度L₄取8.712 m，为两倍艇长. 入口来流速度U取3.051 m/s，出口为压力出口，侧面边界设置为对称平面. 艇加桨的模型计算域设置与单艇一致，如图1（b）所示. 为了更好地捕捉流场结构细节，在桨叶周围及尾部区域都进行加密处理，图2（a）~（c）分别给出单桨、单艇模型及艇桨组合模型计算域纵截面的网格分布示意图.

为了进行网格无关性验证，划分了4套网格，分别为60万网格单元数的粗网格方案A、170万网格单元数的中等网格方案B、270万网格单元数的中等网格C和450万网格单元数的细网格D. 在进速系数J=0.833，n=10 r/s的工况下，对单桨模型进行计算，时均推力系数k_t和扭矩系数k_q如表1所示. 当网格加密到270万的时候，和实验结果较接近，且计算量在可接受范围内。在随后的模拟中均采用C套网格划分方案.

为了研究声压频谱特性在空间的分布情况，在轴向平面和径向平面内分别布置36个测点，如图3所示，它们距离桨盘中心1 m，以10°为间隔. 在流场计算稳定后，调用FW-H声类比方程计算声压.

式（1）为FW-H方程的微分形式，方程右端分别对应于四极子、单极子和偶极子声源的贡献，其中 $\overline {{T_{ij}}} $为广义形式的Lighthill张量. $f$表示控制面，可以取包裹声源区域的内部界面，也可以直接取固壁面；对于马赫数较小（Ma<0.3）的情形，湍流引起的噪声声压级较小，一般直接选用固壁面.

(1) $\begin{split} &\left( {\frac{{{\partial ^2}}}{{\partial {t^2}}} - {c^2}\frac{{{\partial ^2}}}{{\partial x_i^2}}} \right)\overline {(\rho - {\rho _0})} = \frac{{{\partial ^2}\overline {{T_{ij}}} }}{{\partial {x_i}\partial {x_j}}} + \\ &\frac{\partial }{{\partial t}}\left( {\left[ {{\rho _0}{v_n} + \rho ({u_n} - {v_n})} \right]\delta (f)} \right) - \\ &\frac{\partial }{{\partial {x_i}}}\left\{ {\left[ {{P_{ij}}{n_i} + \rho {u_i}({u_n} - {v_n})\delta (f)} \right]} \right\}. \end{split} $

采用Dunn Farassat Padula Formulation 1A方法，对流场声压进行计算. 测点的声压可以表示为单、偶、四极子声源相叠加：

(2) ${p^\prime}(x,t) = p_{\rm{T}}^\prime(x,t) + p_{\rm{L}}^\prime(x,t) + p_{\rm{Q}}^\prime(x,t).$

式中： $p_{\rm{T}}^\prime$表示单极子项，也称为厚度噪声，

(3) $p_{\rm{T}}^\prime (x,t) = \frac{1}{{4{\text{π}} }}\left( {\left( {\frac{\partial }{{\partial t}}} \right)\int_S {{{\left[ {\frac{Q}{{(r(1 - {M_r}))}}} \right]}_{{\rm{ret}}}}} {\rm{d}}S} \right),$

其中 $Q = {\rho _0}{U_i}{n_i}$， ${U_i} = (1 - \rho /{\rho _0}){v_i} + \rho {u_i}/{\rho _0}$， ${u_i}$为流体在 $i$方向的速度分量， ${v_i}$为控制表面在 $i$方向的速度分量，n表示平面法向向量；n_i为n在i方向的分量； $p_{\rm{L}}^\prime$为偶极子项，也称为荷载噪声，

(4) $p_{\rm{L}}^\prime (x,t) = \frac{1}{{4{\text{π}} }}\left( {\left( { - \frac{\partial }{{\partial {x_i}}}} \right)\int_S {{{\left[ {\frac{{{L_i}}}{{r(1 - {M_r})}}} \right]}_{{\rm{ret}}}}} {\rm{d}}S} \right),$

${L_i} = {P_{ij}}{n_i} + \rho {u_i}({u_n} - {v_n})$， ${P_{ij}} = (p - {p_0}){\delta _{ij}} - {\sigma _{ij}}$； $p_{\rm{Q}}^\prime$表示四极子项，在船舶螺旋桨噪声问题中一般忽略不作讨论.

对4~5个旋转周期的声压时域结果作FFT，可得声压级频谱图. 对声压时域结果进行均方根统计，采用下式可以求得总声压级（sound pressure level, L_p）：

(5) $L_p = 10\lg\; {\left( {{{p_{{\rm{rms}}}}}} \left/{{{{p_{{\rm{ref}}}}}}}\right. \right)^2}.$

式中： ${p_{{\rm{rms}}}}$为声压的均方根； ${p_{{\rm{ref}}}}$为参考声压，水中取 $1\;{\text{µ}} {\rm{Pa}}$.

图4给出不同进速系数下桨叶的推力系数、扭矩系数和效率模拟结果与前人实验结果的对比. 图中，k_q,exp、 ${k_{{\rm{q,s}}{{\rm{i}}\rm{m}}}}$分别为扭矩系数的实验值和模拟值，k_t,exp、 ${k_{{\rm{t,s}}{{\rm{i}}\rm{m}}}}$分别为推力系数的实验值和模拟值，η为效率，η_exp、 ${\eta _{{\rm{s}}{{\rm{i}}\rm{m}}}}$分别为效率的实验值和模拟值。从图4可以看出，数值模拟结果与实验值^[20]吻合较好.

按照0.7倍桨叶半径的圆周截取螺旋桨表面截线，与实验测得的压力系数C_p对比，如图5所示. 图中，C为弦长。从图5可知，桨叶表面压力系数分布与实验结果吻合较好.

综上所述，桨叶的流场数值模拟不管是宏观结果还是局部压力结果，都能够较好地和实验结果吻合，验证了流场计算的可靠性.

按照图1（a）的计算域和图2（a）的网格划分方法，对桨叶模型在J=0.833, n=10 r/s的工况进行计算. 桨叶在转动方向旋转1°所需的时间为 $2.7 \times $ $ {10^{ - 4}}$ s，因此声学计算的时间步长取 $2 \times {10^{ - 4}}$ s. 湍流模型选择DES k-omega，采用0.001 s的粗时间步长计算1000步后，再用 $2 \times {10^{ - 4}}$ s的细时间步长计算0.4 s，流场基本达到稳定. 调用FW-H on the fly模型，监测4个旋转周期声压结果，开展频谱分析.

在该工况下，垂直于桨叶水平轴1倍半径处的监测点噪声声压级频谱如图6所示. 图中，f为频率。通过对比发现，现有结果与龚京风等^[7]的模拟结果较接近.

因测点噪声频谱具有一定的对称性，为了便于比较，图7给出轴向平面上4个典型测点的偶极子噪声声压级频谱. 可以看出，除了水平方向的测点，叶片通过频率（BPF）及谐波上均出现了显著高于宽带噪声声压级的线谱噪声. BPF及其谐波f_n的定义如下：

(6) ${f_{\rm{n}}} = n{\rm{BPF}},\;{\rm{BPF}} = {B_0}{\varOmega _0}.$

式中：n为谐波次数， ${B_0}$为叶片数量， ${\varOmega _0}$为叶片角速度（以r/s为单位）.

偶极子线谱噪声是由叶片周期性非定常荷载所产生的. 叶片与湍流相互作用导致的表面力脉动构成低频段的宽带噪声，高频段的宽带噪声和边界层分离密切相关，但边界层流动复杂、难以预测的特点导致高频段宽带噪声难以预测准确^[21].

图8给出轴向平面内4个测点的单极子噪声声压级频谱. 轴向方向上的A0测点在低频段的线谱特征不明显，其他方向上的频谱在BPF及谐波上会出现显著的线谱噪声. 这是由于桨叶旋转引起的周期性排水体积主要在径向方向，在轴向方向上排水体积周期性振荡特征不显著，因此除轴向上的测点都能够捕捉到显著的线谱噪声.

利用式（5）对各个方向测点4个旋转周期的声压结果进行统计，可得总频段声压级，分析单极子和偶极子噪声声压级的指向性分布. 图9、10分别给出单极子、偶极子在轴向平面和径向平面内的分布特征. 图中，L_pM为单极子噪声声压级，L_pD为偶极子噪声声压级. 在轴向平面内，单极子噪声声压级呈现出8字形的指向性，偶极子噪声声压级呈现出 $\infty $型的指向性，这点特征在文献[22]中有描述. 其中，偶极子噪声声压级在水平方向最大可以达到约130 dB，单极子噪声声压级的最大值仅为100 dB，因此偶极子噪声为主要的噪声来源. 只有在特定的方向上，单极子噪声声压级对总噪声声压级有较大贡献. 在径向平面内，偶极子噪声声压级和单极子噪声声压级在各个方向上均呈现出均匀的分布特性.

表2给出桨叶在部分测点上噪声声压级的组成成分. 单极子噪声声压级在垂直方向上更大，偶极子噪声声压级在水平方向上更大. 总体来看，轴向方向噪声主要来源于偶极子噪声，在垂直方向上，偶极子声压级和单极子声压级较接近，两者对总噪声声压级的贡献都应考虑.

suboff单艇模型的计算域如图1（b）所示，图11给出suboff上表面截线的压力系数与实验值^[12]的对比. 可知，现有结果可以和实验值吻合.

表3给出不同的网格划分方案和湍流模型的模拟结果比较. 表中，F_d为总阻力，E为误差。从表3可以看出，采用DES模型和2种不同的RANS模型所得到的suboff表面阻力均能够较好地和实验结果吻合.

图12给出在x/L=0.978的桨盘面处，半径取0.0635 m为圆周的截线上x方向的无量纲速度分布. 可以看出，DES模型较RANS模型能够和实验结果^[11]吻合，RANS模型显现出空间脉动量估计过小的缺陷.

按照图1（b）、2（c）的前处理方法，桨放置在艇模型尾部，研究艇尾流场对桨叶噪声声压级的影响. 入口来流速度U设定为3.051 m/s，转子转速为10 r/s. 湍流模型选择DES SST k-omega，在粗步长计算稳定后调用FW-H模型，频谱分析选用4个周期的结果.

在艇尾部伴流场作用下，桨叶辐射的噪声声压级频谱发生了明显变化，如图13所示，与均匀流场的偶极子噪声声压级频谱图相比，主要有以下2点差异. 1） 在伴流场作用下，各个方向测点的低频线谱噪声声压级趋于一致，因此特定频率下的指向性特征不再明显，在均匀流场下各方向监测点的低频线谱噪声声压级差异更显著. 2） 从频段总声压级上看，伴流场使得低频段宽带噪声声压级有显著的提升，对总偶极子噪声声压级有一定的增加.

线谱噪声声压级的指向性在艇尾部伴流场作用下发生显著的变化，可以从桨叶表面的激振力得到解释，图14分别给出均匀流和伴流场下桨叶表面在轴向（x方向）和径向（y方向）的激振力的变化. 图中，F_x、F_y分别为轴向和径向的激振力。可以看出，在均匀流下，桨叶激振力在径向波动较小，在轴向有较大的波动，2个方向上的规则周期性都不显著；在伴流场下，轴向激振力的均值有一定增加，因此对偶极子噪声的低频宽带噪声有一定提升，同时桨叶在径向和轴向上的激振力都有较明显的周期性，轴向和径向力振幅趋于一致，使得在各个方向测点上的线谱声压级趋于均匀.

对比图15、8的单极子噪声声压级频谱结果可知，在伴流场和均匀流下，单极子噪声声压级频谱呈现出相似的分布特性：轴向测点的单极子声压级较小，垂直方向的单极子噪声的线谱噪声更大. 除了水平方向测点，其他方向上测点的单极子噪声声压级较大，且频谱中出现了对应于叶频及谐波的线谱噪声，反映了除水平方向上测点由于桨叶在径向方向上的排水体积具有显著的周期性，尾流方向测点因为不在桨叶周期性排水方向，单极子噪声的声压级较小.

为了对比尾部伴流场的影响，声源辐射面只选取桨叶面. 在尾部伴流场作用下，单极子噪声声压级和偶极子噪声声压级的指向性特征没有本质改变，分别呈现出8字形和∞形分布，如图16所示.

表4给出伴流场下测点噪声声压级的定量结果。可以看出，伴流场下的单极子噪声声压级基本保持不变，但是偶极子噪声声压级的增加较明显，在垂直方向的测点增加尤为显著，如A9测点偶极子声压级增加5 dB.

从式（3）可知，使用不可穿透面方法的单极子噪声声压级只取决于桨叶的几何特征和运动特性，因此在桨叶转速不变的情况下不会发生变化，该特征导致不可穿透面无法用于计算桨叶空化引起的单极子噪声声压级. 偶极子噪声主要来源于桨叶表面的脉动力，由于艇尾部脱落的涡结构与桨叶的相互作用，使得桨叶轴向压力表面脉动得到增强，因此偶极子噪声声压级有了一定的增加，伴流场下桨叶轴向和径向的激振力的周期性和振幅趋于一致，偶极子线谱噪声声压级的指向性趋于均匀分布.

（1）采用DES和RANS方法，均能对艇和桨叶的宏观流动特征进行准确刻画，但是DES模型在空间脉动量的计算上相较于RANS模型更准确.

（2）在可靠的网格划分和湍流模型方案下，suboff和DTMB4119桨叶模型的水动力学数值模拟结果和实验结果吻合较好.

（3）桨叶在轴向平面内的单极子噪声声压级呈现出8字形指向性分布特征，偶极子噪声声压级呈现出∞型分布；在径向平面内桨叶噪声各个方向上趋于一致. 不管是在均匀来流还是在非均匀流场中该特征都有体现.

（4）在测点噪声频谱中，低频段可以观察到显著的线谱噪声，这些线谱噪声的频率恰好对应于桨叶BPF及谐波频率.

（5）在伴流场作用下，各个方向测点的线谱噪声声压级趋于一致，线谱指向性相较于均匀流场的桨叶不再明显.

（6）艇尾部伴流场增加了桨叶表面的激振力，提升了偶极子声压级，由于不可穿透面方法的简化，单极子声压级不会随着伴流场发生改变.