0 引　言

噪声级是衡量水下航行器（潜艇）的一项重要性能参数。潜艇的噪声分为3类：机械振动噪声、螺旋桨噪声和流噪声。水下航行器噪声产生机理十分复杂，不同部位受到三大噪声源的影响也不同，如流噪声对首部声呐处自噪声的影响较大。对于潜艇装置的声呐，控制流噪声造成的自噪声意义重大。据调查，声呐平台自噪声降低5 dB，本艇声呐探测距离增加60%，探测目标的海区面积为原来的3倍左右^[1]。虽然流噪声的数值比机械振动噪声要低很多，但是随着潜艇速度的不断提高，其在潜艇总噪声中的比重也迅速增加，从潜艇航行隐蔽性而言，降低流噪声也十分重要。

目前，工程上对于流噪声的研究主要分为试验测量和数值仿真。试验研究具有花费大、对试验仪器测量精度要求高、等比例模型试验难以完成等缺点，极大限制了工程上广泛借助试验方法来对水下航行器的流噪声进行研究。随着CFD技术以及计算机资源的极大提升，数值模拟成为可能。数值模拟虽不能完全代替试验测量，但是作为试验测量的补充，在对比分析、了解噪声产生机理、优化选型、为试验减少工况等方面起到了重要作用。数值模拟在一些方面也具有试验不可比拟的优点，例如在数值仿真过程中可以在任意场点上布置虚拟接收器，有效避免了试验中水听器安装困难的问题，也避免了试验中声接收器与流场之间的相互干扰问题。此外，试验中环境噪声或水洞洞壁的声透射等问题在数值模拟中也不存在。

本文首先通过与已有文献计算结果进行对比，说明LES和ACTRAN声学软件联合计算流噪声这种方法满足一般工程的精度要求。随后对2种不同头部线型的水下航行器在不同航速、不同攻角、不同测试点的流噪声进行研究，旨在分析头部线型、来流速度和测试点位置对流噪声的影响，然后与黄桥高^[2]水下航行器流噪声特性水洞试验结果进行对比，本文数值模拟得到的结论与水洞试验相同，进一步说明本文所用的流噪声计算方法可以用于流噪声的对比分析，最后根据数分析给出关于潜艇首部的线型设计和声呐安装位置的一些建议。

1 计算模型 1.1 头部线型说明

本文选择双参数椭圆曲线作为头部线型，双参数椭圆线型的表达式为：

${y^m} = 1 - {\left( {1 - x} \right)^n}\text{，}0 \leqslant x \leqslant 1\text{，}$

(1)

根据式（1），水下航行器头部线型方程为：

$y = \frac{{{D_f}}}{2} + \frac{{D - {D_f}}}{2}{\left[ {1 - {{\left( {1 - \frac{x}{{{L_h}}}} \right)}^n}} \right]^m},0 \leqslant x \leqslant {L_h}\text{，}$

(2)

式中：x为航行器物理横坐标；y为航行器物理纵坐标；D为航行器的最大直径；D_f为航行器前端面直径；L_h为航行器头部长度；m，n为头部参数^{[2, 3]}。

本文选择如图1所示，2种头部线型。

1）D=70 mm，D_f=42 mm，L_h=39 mm，m=2，n=1.85；

2）D=70 mm，D_f=44 mm，L_h=36 mm，m=2.3，n=1.98。

1.2 计算域与计算网格

本文虽然只研究头部线型的噪声性能，但是建模时不能做截断处理只建立头部的线型，为了避免截断面后会出现回流的影响，这里将模型平滑延长，模型总长为700 mm，如图2所示。根据水下航行器的实际尺寸，综合考虑计算量、以及所需的流场详细信息的准确性。流场域的大小满足：入口在航器首部1倍长度处，出口在航器尾部2倍长度处^[4]，流场域选择圆柱域，半径为航行器半径的10倍，见图3。为了节省计算资源，网格划分在航器周围采用外O-BLOCK的技术形成加密区，以达到尽可能获得多的航行器头部流场流动信息的目的。整个模型采用结构化的六面体网格，在网格生成过程中，分布应该稀疏合理。靠近边界层，网格太稀不利于捕捉到流动细节；物面附近参数变化梯度大，网格太稀则与艇体实际形状误差大；离艇体较远处，参数变化小，网格太密则浪费计算资源。图4为头部的局部网格示意图，可见航行器壁面以及头部的网格尺寸很小，以便捕捉详细的边界层流场信息。计算域网格总数为453 040，节点总数为443 526，为了准确捕捉壁面信息y⁺取50。

2 计算原理 2.1 湍流模型

对于流噪声的预报，其不可压缩粘性湍流流动优先采用大涡模拟（LES）。LES的基本思想可概括为：用瞬时的N-S方程直接模拟湍流中的大尺度涡，对于小尺度涡的处理采用近似的模型来考虑。相对于RANS方法对控制方程中的时间量做平均得到脉动量的时均值的方法，LES可以为噪声计算提供更详细的更全面的流场信息，从而使得LES更适合于声学计算。要实现LES必须完成以下2个重要环节的工作：一是建立一种数学滤波函数；二是加入被滤掉小涡对大涡运动的影响，LES中通过加入附加应力项实现，该应力项称为亚格子尺度应力。

LES方法中，将任意流体变量U分解为大尺度的平均分量 $\overline U $ 和小尺度分量 $U'$ ，即

$U = \overline U + U'\text{。}$

(3)

空间中位于x₀处的滤波后的变量定义为

$\overline U \left( {{x_0},t} \right) = \int {U\left( {x,t} \right)G\left( {{x_0},x,\Delta } \right)} \text{。}$

(4)

式中：G为滤波函数。滤波后参量的结构和大小由滤波函数决定，而滤波函数依赖于 ${x_0} - x$ 以及滤波宽度 $\Delta = {\left( {{\Delta _1}{\Delta _2}{\Delta _3}} \right)^{1/3}}$ ， ${\Delta _i}$ 其中为第i个空间坐标中的滤波宽度。

亚格子尺度模型采用Smagorinsky模型，该模型是基于均衡假设，即小尺度参量将从大尺度参量吸收的全部能量完全地、瞬时地耗散干净。根据Smagorinsky的基本SGS模型，SGS应力具有下面的形式：

${\tau _{ij}} - \frac{{{\tau _{kk}}{\delta _{ij}}}}{3} = - 2{\mu _t}{\overline S _{ij}}\text{。}$

(5)

式中：μ_t是亚格子尺度的湍动粘度；亚格子尺度应力τ_kk的各向同性部分既可以通过建模添加到滤波后的压力函数中，也可以忽略不计^[5]。

亚格子尺度的湍动粘度为

${\mu _t} = {\left( {{C_s}\Delta } \right)^2}\left| {\overline S } \right|\text{。}$

(6)

式中： $\left| {\overline S } \right| = \sqrt {2{{\overline S }_{ij}}{{\overline S }_{ij}}} $ ；C_s为Smagorinsky常数，通常值为0.1左右时可以得到最佳的结果。

2.2 流噪声方程

从声学基本概念来讲，声学其实是流场的不平衡扰动。通过紧致声源的假设，1952年Lighthill^{[6, 7]}提出声学类比理论，使噪声的求解成为可能。该理论将声学计算和流场计算结合在一起，通过求解流场信息进而求解流噪声。Lighthill方程是由N-S方程导出，得到密度波动方程：

$\frac{{{\partial ^2}\rho '}}{{\partial {t^2}}} - c_0^2{\nabla ^2}\rho ' = \frac{{{\partial ^2}{T_{ij}}}}{{\partial {x_i}\partial {x_j}}}\text{。}$

(7)

式中：c₀为流体介质中等熵条件下的声速； $\rho ' = \rho - {\rho _0}$ ，其中ρ和ρ₀分别是扰动与未扰动时流体的密度；T_ij为Lighthill应力张量，可表示为

${T_{ij}} = \rho {u_i}{u_j} + {P_{ij}} - c_0^2\left( {\rho - {\rho _0}} \right){\delta _{ij}}\text{。}$

(8)

本文所用的ACTRAN对流噪声的预报是基于Ligthhill声类比的方法，其对波动方程进行了改进如下：

$\frac{{{\omega ^2}}}{{{c^2}}}\psi + \frac{{{\partial ^2}\psi }}{{\partial {x_i}\partial {x_i}}} = \frac{1}{{i\omega }}\frac{{{\partial ^2}{T_{ij}}}}{{\partial {x_i}\partial {x_i}}}\text{，}$

(9)

Ligthhill应力张量则改写为式：

${T_{ij}} = \overline \rho {\overline u _i}{\overline u _j} + {\overline \tau _{ij}} + \left( {\overline p - {c^2}\overline \rho } \right){\delta _{ij}}\text{，}$

(10)

Ligthhill声类比方程的弱变分方程为：

$\begin{split}& \int_\varOmega {\displaystyle\frac{{{\Omega ^2}}}{{{c^2}}}} \psi \delta \psi {\rm d}\varOmega - \int_\varOmega {\displaystyle\frac{{\partial \psi }}{{\partial {x_i}}}\displaystyle\frac{{\partial \delta \psi }}{{\partial {x_i}}}{\rm d}\varOmega } = \\& \int_\varOmega {\displaystyle\frac{i}{\varOmega }\displaystyle\frac{{\partial \delta \psi }}{{\partial {x_i}}}\displaystyle\frac{{\partial {T_{ij}}}}{{\partial {x_i}}}{\rm d}\varOmega }- \\& \int_\Gamma {\displaystyle\frac{{\delta \psi }}{{i\varOmega }}\displaystyle\frac{\partial }{{\partial {x_i}}}\left( {{c^2}\rho {\delta _{ij}} + {T_{ij}}} \right){n_i}{\rm d}\Gamma } \text{。}\end{split}$

(11)

其中：方程右端第1项表征体声源对声场的贡献，第2项表示面声源对声场的贡献。体声源包括3个部分：对流输运的散度、熵源项、粘性应力的散度。ACTRAN在计算时，忽略了熵源项和粘性应力。在ACTRAN计算中，壁面边界默认为刚性壁面，此时面源项等于0，对应与式（11）中右端第2项为0。

3 流噪声计算方法验证

本文流噪声的计算首先采用LES方法获得三维流场的瞬态信息，然后基于Lighthill的声类比理论对流噪声进行模拟。本文对研究数据较丰富的“大青花鱼”简化潜艇模型的三维流场及流噪声进行了研究，得到了x方向声压总级变化特性对比。“大青花鱼”简化模型尺度：艇长L=3.2 m，型宽B=0.4 m。流噪声对比计算结果见表1，3种结果流噪声的辐射特性基本一致，在沿长度方向上呈逐渐递减的规律。对于离艇体较近的区域声压2种不同方法（BEM和Actran）计算结果存在一定差别，而当艇体距离不断增大时，3个结果具有很高的吻合度。本文计算结果与孟堃宇^[8] 、刘继明^[9]研究数据进行对比可知，除了在艇体表面较近处，本文采用的计算结果满足工程要求的精度（误差少于3 dB）；而对于流噪声的计算，ACTRAN采用的是有限元和无限元来模拟声源，而Sysnoise是采用边界元，从所用方法来说，ACTRAN得到的近场结果更加准确。综上讨论，本文选用对近场和远场流噪声模拟具有更高精度的ACTRAN。

4 计算流程及数值结果 4.1 流场计算

本文划分好网格后，将其导入Fluent中进行流场计算。流场数值模拟的计算步骤为：首先给定来流速度，应用RNG k-ε进行定常计算。为提高计算精度，计算时对动量、湍动能的处理采用了二阶迎风差分格式，压力速度耦合采用了SIMPLE算法。待定常计算充分收敛之后，通过阻力性能来验证流场模拟的精度；在比较满意的条件下以此定常解作为初值进行LES非定常计算。LES采用Smagorinsky-lilly粘性模型，压力速度耦合选择SIMPLEC算法，压力空间离散采用PRESTO!方法，动量的处理采用有界中心差分（Bounded Central Differencing）。下面仅对第1种头部线型在8 m/s时的流场计算进行分析。

根据桑海公式，对于Re在1.0E+06到1.0E+09之间的估算公式：

${C_f} = \frac{{0.4631}}{{{{\left( {\lg {\mathop{ Re}\nolimits} } \right)}^{2.6}}}}\text{。}$

(12)

将 ${\mathop{Re}\nolimits} = \frac{{\rho \nu L}}{\gamma } \approx \frac{{1000 \times 8 \times 0.7}}{{0.001003}} = 5.60{\rm E} + 06$ 代入上式，计算得： ${C_f} = 3.20{\rm E} - 03$ 。从粘阻系数可以看出，流场计算具有较高的精度。

4.2 流噪声计算

在非定常稳定后，保存5 000步流场结果，时间步长为1.0E–04 s，取该值的原因为：LES需要满足CFL在1附近的条件，其值的大小由最小求解涡的时间尺度来决定，若时间步长取太短固然能使CFL数小于1但会造成迭代步数增多，消耗大量时间的同时也会积累大量的舍入误差，综合考虑计算资源，所以时间步长取为1.0E–04 s。根据采样频率及采样定理，转换得到的频域流场信息分辨率为2 Hz，带宽是5 kHz。由于声学计算软件ACTRAN基于有限元、无限元的方法，则离散分为2个部分，一个是远场的无限元区域，一个是近场的有限元区域；无限元不需要划分网格，在ACTRAN里直接设置其阶数即可（一般对于近场网格比较密的情况下，选择默认设置为5）。这里ICEM的网格划分只针对有限元区域，声场中网格大小的划分和分析频率密切相关。分析频率越高，波长越短，按最高频率计算波长，要求每个波长范围内至少要有6～8个线型单元。其计算模型设置如图6所示：在声源域，将有限元网格灰色区域的外表面设置为无限元基面，潜艇表面采用默认设置（隔声面）。布置2个接收点，A点位于头部驻点处，B点位于头部端面15 mm半径处。

利用Fluent进行三维流场分析，通过ACTRAN中的ICFD模块，将流体计算结果直接导入到ACTRAN中求解流噪声辐射问题。ICFD为ACTRAN的中将CFD结果转化为声源的模块，首先将流体计算中速度、压强等信息通过积分算法插值到声学网格上，得到声学模型的时间历程载荷，再利用其内部的傅里叶变换算法将时间历程载荷转化为频谱载荷，最后利用ACTRAN进行求解。

图7给出了模型1和模型2当来流无攻角时，在2个不同速度下，2个不同接收点处的流噪声声压级特性曲线。

表3是根据声压级特征曲线计算出的声压总级值（单位为dB）。计算公式如下：

$OSPL = 10\lg \left( {\Delta f\sum\limits_{f = f0}^{{f_n}} {{{10}^{0.1{L_p}\left( f \right)}}} } \right)\text{。}$

(13)

式中：f对应声压级特性曲线x轴的各个频率值， ${L_p}\left( f \right)$ 为频率值对应的声压级。

1）速度对流噪声的影响分析

图7给出了模型1和模型2在来流无攻角时不同速度下2个测试点的流噪声声压级频谱特性曲线。从图中可以看出，流噪声声压级随着速度的增大而增大。声学能量主要集中低频下，随着频率的增大，流噪声的声压级具有下降的趋势。分析原因普遍认为，转捩区的声辐射是构成水下航行器流噪声的主要组成部分，随着来流速度的增加，转捩点位置前移是导致辐射到相应位置在噪声增大的原因。根据表3，模型1和模型2不同攻角情况下，在10 m/s时的声压总级比8 m/s时的声压总级有不同程度的升高，最高不到8 dB。

2）头部线型对流噪声的影响分析

根据表3，对比模型1和模型2的声压总级可以发现，相同速度下模型2的噪声总声压级要比模型1大，分析原因认为，首部流噪声主要源于流场转捩区的声辐射，模型1的头部端面直径小、且头部曲线段长度长，这种线型使得最小压力点和边界层转捩点（湍流起始点）的位置后移，这样辐射到模型头部相应监测点处的噪声要小。

3）测试点位置对流噪声的影响分析

表3给出了模型1在0°，2°，4°攻角和模型2在0°攻角时不同测试点在不同速度下的流噪声总声压级。头部端面驻点处的声压总级比15 mm半径处的声压总级大约小11～17 dB，这里与黄桥高试验在差值大小上有一定的差别，这里分析原因为：黄桥高试验统计频率范围为15～35 kHz，转捩区的辐射声波入射到模型头部端面时，声波频率高的时候，声波的波长比较短，则需要考虑声波传播的绕射损失。而本文由于考虑到数值模拟分析最高频率与声波波长有关，在生成声学计算网格需要划分很密，超出计算机的计算能力，本文分析的最高频率选为5 kHz，这里转捩区的声传播过程中绕射损失要小很多。相比头部端面其他位置，驻点处的流噪声声压总级最小。

5 结　语

本文在不同速度、不同攻角下，分别对2种不同头部线型的水下航行器流噪声进行了数值仿真，对计算结果分析可得出以下结论：

1）在不同速度下，声学能量主要集中低频下。随着频率的增大，流噪声的声压级具有下降的趋势。对于同一线型的头部，流噪声随着速度的增大，呈现出升高的趋势。不同线型的头部在同一速度下的流噪声大小不同，就本文分析的2种线型而言，线型1具有较好的噪声性能，同时结合流噪声声源主要来自于转捩区的声辐射这一共识，选择头部端面面积小而头部较长的线型具有更好的噪声性能。

2）对比各个速度各个线型A接收点（驻点）的流噪声明显小于B接收点的流噪声，那么在潜艇首部声呐的布置上，从降低自噪声的出发点来看，尽量把声呐布置在驻点附近，有利于提升声呐的探测效能。

3）将数值仿真结果与试验结果进行对比，说明本文计算流噪声所用的方法在精度上具有一定的说服力。将这种数值仿真方法用于试验前的优化设计可以减少试验方案，从而使设计者跳出考虑试验高花费的限制，为大胆尝试各种不同的线型提供了便利。