流体激励下泵喷推进器声振耦合响应数值分析

高丹妮，余海廷，华宏星

（1.上海交通大学 振动冲击噪声研究所，上海 200240；2.高新船舶与深海装备协同创新中心，上海 200240）

泵喷推进器是一类包含前/后置定子的导管推进器，前置定子可均化和预旋来流，主要用于鱼雷中，后置定子可吸收转子产生的旋转能，多用于潜艇或水面舰船的低噪声推进[1]。然而在流体中螺旋桨的转动也会引起结构表面的脉动压力，从而产生结构振动并辐射噪声，影响舰船的作战声隐身性能。因此，研究泵喷推进器在流体激励下的声振耦合问题尤为重要。计算流体力学（CFD）是一种基于流动方程的数值计算方法，计算控制方程可以由Navier-Stokes方程（N-S方程）来描述。直接数值模拟（DNS）方法、大涡模拟（LES）方法和雷诺平均（RANS）方法是解决N-S方程的三种主要方法[2]。早在1963年，Smagorinsky J.S就提出了大涡模拟的概念和 Smagorinsky 亚格子模型[3–5]。Mauro 等人[6–8]模拟了基于耦合有限元/边界元方法的带锥形圆柱壳潜艇结构的谐响应和声学特性。Wei等人[9]利用有限元/边界元方法研究了激励力作用下潜艇的振动响应和声辐射特性。现有研究中对于泵喷推进器的流激振动和声振耦合问题的研究并不多见。本文对某泵喷推进器进行了CFD数值计算，得到了结构表面的脉动压力，并以此作为激励源施加至其水下结构声振耦合有限元/边界元模型中，得到其声辐射规律。本文也研究了不同进速系数对于水动力性能及声振耦合响应的影响。

1 基本理论

1.1 大涡模拟

N-S方程是描述湍流运动的基本控制方程[10]。DNS方法可以直接求解N-S方程，但计算量过大；LES方法能够滤掉湍流中比网格尺寸小的脉动，在一定程度上减少了计算量；RANS方法是求解时均化的N-S方程，计算量最小。

LES方法相对于DNS方法其对于时间和空间的分辨率要求较低，相对于RANS方法其能提供更多的流动信息，且能在现有计算条件下得到表征复杂湍流和高雷诺数流动信息，能精确可靠地反映湍流的运动过程。考虑到本文对于计算精度的要求及现有计算资源，对于泵喷推进器的CFD计算将采用LES方法。

LES通过空间滤波技术分离不同尺度的涡，对脉动压力影响较大的大尺度涡直接计算，对小尺度的涡进行过滤。大尺度涡可以由f(x′,t)的加权积分来描述[11]：

其中：D是整个流域的计算域，x′和x分别是滤波前后的向量，是滤波函数。

滤波之后连续性方程和不可压缩N-S方程可以表示为[11]

其中：ui是与xi有关的速度分量，是滤波之后的平均速度分量，ρ是流体密度，是亚格子雷诺应力。τij反映了大尺度涡和小尺度涡之间的相互作用，如能量、动量转化等。Smagorinsky模型，动态Smagorinsky-Lilly模型和WALE模型是较为常用的三种亚格子模型。本文使用动态Smagorinsky-Lilly模型模拟亚格子雷诺应力。

1.2 声振耦合响应

对于声振耦合问题，在低频域经常使用耦合有限元/边界元方法来求解，有限元模型一般用来描述结构动力学响应，边界元模型用来计算辐射声场。

Helmholtz积分方程可以描述为[12–13]

[A]、[B]是系数矩阵，{p}为声压向量，{vn}是结构表面法向速度。

对于结构域，将声压看作结构载荷，则声振耦合水下结构的振动方程可以描述为[12–13]

K、C及M分别为结构的刚度、阻尼及质量矩阵，{u}为位移向量，{Fs}是结构的载荷矩阵。矩阵Lc是耦合矩阵。

结合有限元/边界元方程及速度、位移关系{vn}=jω{u}等方程，可以得到耦合方程

求解上述方程，得到结构振动位移及辐射声场声压。

2 数值计算结果分析

2.1 模型描述

本泵喷推进器模型如图1所示，由导管、转子（螺旋桨）、前置定子、轴等部件组成。模型基本参数如表1所示。

表1 泵喷推进器几何参数

图1 泵喷推进器模型

建立圆柱形CFD计算域，外流场直径为3D，长为5D，如图2所示。

图2 泵喷推进器CFD计算域

将计算域分为旋转域和静止域两部分，旋转域模型如图3所示。

图3 旋转域模型

包括转子、桨毂及外流场等。其余部分为静止域，包括定子、导管、部分轴及外流道等。图4(a)、图4(b)分别显示了静止域的单流道模型（九分之一流道模型）及全流道模型。

定义进速系数[10]

图4 静止域水动力模型

其中：Va为螺旋桨进速，n为螺旋桨转速，D为螺旋桨直径。本文分别对进速系数J=0.4（工况1：Va=1.5 m/s，n=750 r/min，D=300 mm）和J=0.5（工况2：Va=1.5 m/s，n=600 r/min，D=300 mm）两个工况进行了计算。

2.2 导管脉动压力分析

流体网格的划分在ICEM中进行，并将边界层进行网格加密处理。旋转域由于螺旋桨几何的复杂性，采用非结构网格来划分，如图5所示。

静止域采用结构网格来进行划分。对单流道进行网格划分，而后对其进行周期旋转，得到全流道网格，如图6所示。旋转域和静止域通过网格交界面来进行数据传递。

图5 旋转域网格示意图

图6 静止域网格示意图

流体域参数设置如表2所示，整个计算域网格数量约为550万。设置边界条件为速度入口和压力出口，时间步长为1×10-4s，计算的频率区间为0～500 Hz。输出时间长度为0.8 s，频率分辨率为1.25 Hz。

表2 流体域参数设置

对导管内表面全部节点上的脉动压力时间序列进行傅里叶变换，得到J=0.4工况下的叶频处（50 Hz）脉动压力幅值云图如图7所示。可以看到，脉动压力在空间上呈轴对称分布，压力脉动峰值集中出现在螺旋桨所在区段即x=0.003 m附近，且其数量级为其他区域10倍以上，因此本文主要考虑该区域的流体激励力激发的结构振动及辐射声影响。

图8(a)、图8(b)分别显示了两个工况下导管内壁在x=-0.005、0.003及0.012 m处的脉动压力频谱图，此压力也将作为流体激励力施加至水下结构声振耦合系统的有限元/边界元模型中。

图7 导管内表面脉动压力分布云图

图8 导管内壁x=-0.005、0.003及0.012 m处脉动压力

从图8中可以看到，特征频率出现在叶频（Blade passing frequency，BFP）及其倍频处。工况1和工况2的平均推力系数Cd分别为0.194和0.113。对于工况1，即J=0.4，特征频率出现在50 Hz及其倍频处，对于工况2，及J=0.5，特征频率出现在40 Hz及其倍频处。频谱幅值整体趋势上由低频到高频依次减小，频率峰值主要集中在300 Hz以下的低频段，工况1的峰值大于工况2。

2.3 泵喷推进器模态分析

泵喷推进器结构的主要参数如表3所示。结构网格数量约为13.7万，如图9所示。将结构外围包裹水体，计算结构湿模态，其有限元模型如图10所示。表4给出了结构的前5阶湿模态频率，其中前3阶频率在所分析频段0～500 Hz以内。

表3 结构参数

图9 结构网格

图10 湿模态计算有限元模型

表4 结构湿模态频率

对于第1阶模态（136.55 Hz），结构的变形主要由定子的变形引起，可以认为是定子模态；第2、第3阶模态频率（305.98 Hz、306.19 Hz）极为接近，且模态振型表现出一定的轴对称性，为导管模态，如图11所示。

图11 结构前3阶湿模态

将流体的特征激励频率和结构模态频率对比如表5、表6所示。可以看到，在两个工况中，结构的第1阶湿模态频率都接近于3 BPF，第2、第3阶湿模态频率分别接近于工况1中的6 BPF和工况2中的8 BPF。

表5 特征频率与结构模态频率对照表（工况1）

表6 特征频率与结构模态频率对照表（工况2）

2.4 导管结构声振耦合响应分析

采用耦合有限元/边界元方法进行流体激励下导管振动辐射声场的计算分析。提取导管结构表面网格作为声网格，网格数量约为0.8万。图12显示了两个工况下导管结构的均方振速，可以看到其峰值出现在叶频及其倍频处以及前3阶湿模态频率处，这是因为流体激励力包含叶频及其倍频的特征频率，使得结构振动响应表现出受迫振动响应的特性；同时，当宽频流体激励力的频率接近湿模态频率时将引起结构共振，产生相应的均方振速峰值。工况1（J=0.4，750 r/min）的峰值大于工况2（J=0.5，600 r/min）的峰值。

图12 均方振速频谱

图13显示了结构振动引起的辐射声功率，可以看出辐射声功率随着频率升高有增大的趋势，且工况1的峰值大于工况2的峰值，一定程度上暗示了转子转速增大时，导管的辐射声功率也会随之增强。

图13 辐射声功率频谱

3 结语

本文针对某泵喷推进器模型，采用大涡模拟及耦合有限元/边界元数值计算方法，对不同进速系数（J=0.4和J=0.5）下导管内壁脉动压力和流体激励下的声振耦合响应进行了计算分析，得到了导管内壁脉动压力幅频特性曲线及均方振速与辐射声功率幅频特性曲线。结果表明，导管内壁脉动压力呈轴对称分布，且最大值集中在转子区域。流体激励力特征频率出现在叶频（BPF）及其倍频处，均方振速的峰值出现在叶频及其倍频处以及湿模态频率处。随着频率升高尽管脉动压力的幅值下降，声辐射效率并未下降，在较高频率处仍有较大的辐射声功率。根据以上分析结果，在进行泵喷推进器结构设计中，须着重关注结构本身的动力学特性，避免推进器湿模态与流体激励力特征频率（转子叶频及其倍频）间出现耦合，以防流激共振现象的出现。本文的研究内容可以为工程中同类的泵喷推进器的减振降噪设计提供参考。