基于圆柱绕流的气动声源识别方法

王毅刚，朱朗贤，焦 燕，张 昊

（1.同济大学上海地面交通工具风洞中心，上海 201804；2.上海市地面交通工具空气动力与热环境模拟重点实验室，上海 201804）

高速运动物体的固体壁面和空气相互作用，在近壁湍流区域产生气动声源并辐射声能，从而形成噪声干扰。因此，气动声源及其声辐射和传播一直都是声学领域的研究热点。然而，目前相关研究还不充分。

对气动声源的研究，Lighthill［1-2］最早提出声比拟思想，由流体连续性方程、动量方程（N-S方程）和气体状态方程得到声比拟方程。声比拟方程左端项反映了声的传播特征，右端项为声源项，表征了声源的强度。然而，该方程难以求解复杂流场，不能直接反映声源的声辐射和传播特性。由于声比拟方程只能求解自由空间中由湍流引起的声问题，因此Curle［3］利用基尔霍夫方法将Lighthill理论进一步推广到了固体边界，并将表面偶极子声源纳入考虑范围。Ffowcs Williams等［4］进一步考虑了运动的固体边界，得到了FW-H方程。以上研究考虑了不同边界的影响，方程的声源项将声源区分为单极子、偶极子和四极子的基本声源，同时也揭示了这些声源产生的本质，但只能定性地对声源进行描述，难以对具体流动发声现象进行定量描述。虽然数值仿真技术得到了较大的发展，但是FW-H方程仍仅可用于远场噪声计算，难以求解近声场及流场内的声源及其声辐射和传播。Powell［5］提出了涡声理论，通过引入涡量概念，重新定义声源项，并揭示了涡与声的关系。涡声理论具有与FW-H方程一样的特征。随后，出现了试图描述流场内声场的研究。Bogey等［6］通过分解流场中的物理量并代入Euler方程，得到线性化的欧拉方程。Ewert等［7］在此基础上推导出了声波扰动方程（acoustic perturbation equations，APE），可以有效防止不合理涡量的产生，但由于该方程将流体的扰动作了近似分解，因此对近声场问题的求解是近似的。考虑到声学方程（FW-H方程、涡声方程）的源项反映了一定的声源强度特征，结合流体数值仿真技术，本研究团队开展了流场非定常计算，不涉及噪声的辐射和传播，直接从源项描述声源特性［8-9］。

将固体壁面与空气的相互作用所产生的气动声源看成是无数个微球形声源的组成，基于球形声源声辐射公式和非定常流场特性，建立声场、气动声源和流场之间的关系，获得声源辐射声压和流场压力梯度及脉动速度的关系，据此关系可以确定声源的大小和分布。基于圆柱绕流风洞实验和数值仿真技术，对气动声源的特性进行描述，阐述声源辐射声压和流场压力梯度及脉动速度的关系。

1 声源识别基本思想

以圆柱绕流为研究对象，其气动声源具有典型的偶极子和四极子声源特征，因此主要研究偶极子和四极子气动声源的识别方法。

1.1 单个偶极子声源的声辐射计算

假设流体流动为理想、不可压缩流体流动。理想、不可压缩流体流动的偶极子声源由流体和固体表面相互作用产生，分布于近壁面的流场区域。假设流体流动过程中产生的气动声源由无数个微小的球形声源组成，从远场看，每个球形声源为一个力点源，则该力点源的声压波动方程［10］为

式中：p、c0、t分别为声压、声速和时间；f为力点源单位体积脉动力。

对于单个力点源，方程（1）的解为

式中：x和y分别为源和观测点的位置坐标；R3为三维空间积分域；fi为i方向上的脉动力；yi为y在i方向上的值。

对于单色力源辐射的声压有

经过一系列公式推导可得

式（3）和式（4）中：ω，k分别为角频率和波数；∇x为哈密尔顿算子；f=ei fi，其中ei为i方向上的单位向量。根据式（3）将式（4）描述成非单色声源的形式，即：

式中：xi为x在i方向上的值。利用链式微分法则，对式（5）进一步求导，可以得到单个偶极子声源向远场辐射声压的计算式，如下所示：

式中：r为声源与远场测点的距离。

1.2 单个四极子声源的声辐射计算

四极子声源可以等效为2个距离很近、相位相反、强度相等的偶极子声源的叠加。采用2个大小相等、方向相反的力点源f和-f组成四极子声源，如图1所示。单个力点源所致声场如式（4）所示，则2个力点源产生的四极子单色声源声场为

图1 由一对力点源构成的四极子声源Fig.1 A quadrupole sound source composed of a pair of force point sources

令δy为2个力点源的距离，2个力点源的位置分别为y1=y-δy/2和y2=y+δy/2，2个力点源组合得到的四极子声源声场解［10］为

同样，对于非单色声源有

由式（9）得到非单色四极子声源的声压解，如下所示：

利用复合函数链式求导法则，可以将Tij(y，t- ||x-y/c0)对接收点的位置导数转换成对时间t的时间导数。因此，声源函数是一个反映湍流应力Tij(y)=ρsuiuj随时间变化的量（ui、uj分别为i和j方向上的脉动速度），具有实际物理意义。由式（10）得到单个四极子声源向远场辐射声压的计算式，如下所示：

式中：ρs为湍流区内的平均密度，近似为一个常数；ur为球形声源径向振动速度。

1.3 气动声源的确定

式（6）和式（11）分别给出单个偶极子和单个四极子气动声源的辐射声压表达式。由于流场中偶极子声源和四极子声源有无数多个，因此求2种声源的声辐射总效应十分困难。若令式（6）和式（11）中距观察点的距离r趋向于声源半径，则p就反映了该声源辐射声压的能力，但不能反映声传播特性。进一步将流场中任一微元作为声源，对声压进行求解，得到声压在流场中的大小和分布特征，即给出了声源的特征，由此建立声源识别方法。

2 圆柱绕流的远场气动声辐射

圆柱绕流具有典型的偶极子和四极子声源特征，为了呈现上述声源识别方法对声源的描述效果，开展了圆柱绕流的风洞实验和仿真计算。

2.1 风洞实验

实验在同济大学的气动-声学风洞实验平台开展。研究对象是高度为1.8 m、直径为0.1 m的光滑实心圆柱，如图2所示。该圆柱底部采用圆台结构，顶部采用半球形结构，以减少由根部和顶部流动的不稳定性所产生的附加气动噪声干扰，突显圆柱绕流产生的噪声特征。实验中来流风速为80 km·h-1。

图2 圆柱模型（单位：mm）Fig.2 Cylindrical model（unit：mm）

在喷口射流剪切层外平行于射流中心线的位置等间距地布置3个传声器测点，用于采集远场声压信息。圆柱距风洞中心线0.5 m；传声器距风洞中心线5.0 m，距地面1.2 m，间距均为1.3 m，2号传声器与圆柱的连接线垂直于风洞中线。具体布置如图3所示。采用丹麦GRAS公司的自由场传声器进行声压测量；数采及分析设备为德国HEAD Acoustics公司的SQLABШ多通道采集和分析系统。实验中采样频率为48 kHz，数据分析的频率分辨率为16 384。

图3 测点布置Fig.3 Layout of measuring points

2.2 圆柱绕流仿真计算

根据风洞实验中圆柱的几何特性建立相同的仿真计算物理模型。圆柱直径D=0.1 m，高度为18D，计算域为35D×16D×28D。圆柱中心距入口边界10D，保证来流和两侧及顶部不受边界影响，尾流流动充分发展。圆柱计算域如图4所示。

图4 圆柱计算域Fig.4 Cylindrical computational domain

选择SST（Menter）K-Omega模型进行流场定常仿真计算，采用大涡模拟（LES）方法及WALL Subgrid Scale亚格子尺度模型进行流场非定常计算。入口风速为80 km·h-1，近似认为空气不可压缩，空气密度设置为常数1.18 kg·m-3，湍流强度为1%，湍流黏性比为10。设置速度入口为入口边界条件，压力出口为出口边界条件，计算域的上端面和2个侧面为对称边界条件，其余表面为壁面边界条件。先进行大时间步长（5×10-3s）的非定常仿真计算，时长1.0 s，让流场快速稳定，最后进行小时间步长（2×10-3s）的非定常仿真计算，时长3.2 s。

2.3 仿真与实验结果对比

图5为风洞实验和仿真计算获得的3个测点声压级频谱。可以看出，峰值频率仿真值为45 Hz，比实验值高4 Hz，峰值声压级仿真值为66.0～66.8 dB（A），比实验值小1.0～1.5 dB（A），两者的变化趋势较为一致。计算结果和实验结果存在一定的误差，原因为该流动雷诺数在105量级，处于流动变化较剧烈区域，计算容易出现偏差，但该偏差没有影响流动的主要特征。因此，认为仿真计算能够满足后续研究的要求。

图5 有限长圆柱仿真与实验声压频谱Fig.5 Simulation and experimental results of sound pressure spectrogram for finite-length cylinder

3 声源特性

3.1 偶极子声源特性分析

式（5）中，f为作用于位置y处流体上的力矢量。根据数学知识可知，矢量函数f可以是一个标量的梯度∇p，也可以是一个矢量的旋度∇×u，又可以是两者之和［10］，即：

图6 时均分布Fig.6 Time-average distribution of

图7 ∇p时均分布Fig.7 Time-average distribution of∇p

图8 圆柱表面角度分布示意图Fig.8 Schematic diagram of angular distribution of cylindrical surface

图9 60°延长线上和∇p的时均分布Fig.9 Time-average distribution of and∇p on 60°extension line

图10 100°延长线上和∇p的时均分布

Fig.10 Time-average distribution ofand∇pon 100°extension line

图11 180°延长线上和∇p的时均分布

Fig.11 Time-average distribution ofand∇pon 180°extension line

图12 偶极子各源项声压级对比Fig.12 Comparison of sound pressure level for each source term of dipole

3.2 四极子声源特性

同理，根据式（11）对四极子3个源项的声辐射特性展开分析。计算单个源项在流场中的分布，如图13～15所示。从图13～15可以看出，四极子声源各项具有相似的分布特征，靠近分离点区域数值较大，与偶极子声源相比，四极子声源更远离壁面分离点，并分布于圆柱尾流区。

图13 时均分布Fig.13 Time-average distribution of

图14 时均分布Fig.14 Time-average distribution of

图15 时均分布Fig.15 Time-average distribution of

图16 60°延长线上和u2r的时均分布Fig.16 Time-average distribution of and u2r on 60°extension line

图17 100°延长线上和u2r的时均分布Fig.17 Time-average distribution of and u2r on 100°extension line

图18 180°延长线上和u2r的时均分布Fig.18 Time-average distribution of and u2r on 180°extension line

综上所述，基于式（6）和式（11）的单个声源声辐射计算方法，结合流场非定常数值计算手段，建立气动声源识别方法。进一步结合圆柱绕流声源特性分析，得到向远场声辐射的声源识别方法。

偶极子声源声压为

偶极子声源是由固体壁面和流体相互作用而在近壁面的流体中形成的，依据偶极子声源定义，其出现在壁面附近并呈现出面声源特征，如图6所示。

四极子声源声压为

由图12和图19可知，在圆柱表面附近偶极子声源强度（134 dB）和四级子声源强度（103 dB）相差31 dB，因此后者辐射的声能量可以忽略。

图19 四极子各源项声压级对比Fig.19 Comparison of sound pressure level for each source term of quadrupole

3.3 声源识别方法的有效性

式（14）和式（15）给出了偶极子声源和四级子声源的识别方法，据此给出了如图6和图13所示的圆柱绕流对应的声源分布云图。为了说明本方法的有效性，利用APE（acoustic perturbation equations）方法对圆柱绕流辐射近声场（声压级云图）进行计算。图20a为圆柱绕流的辐射声场，靠近圆柱壁面声压级最大位置表征为声源位置。因圆柱绕流偶极子声源能量远大于四级子声源，四级子声源辐射声能可以忽略，所以该位置主要对应偶极子声源。图20b为本方法确定的偶极子声源位置。2种方法声源位置相同，从一定程度上说明本方法的有效性，但前者声源区域较大，后者声源区域更小。图20a也说明了四级子声源辐射较弱，和本研究结论相同。由于目前气动声源实验（即使是Beamforming声源识别方法）较难捕捉频率低、分布复杂的声源，也较难区分偶极子和四极子声源，因此声源识别的验证工作有待进一步开展。

图20 2种方法确定的声源位置Fig.20 Sound source location determined by two methods

4 结论

（1）根据流场中气动声源辐射声压和脉动力、脉动力和压力梯度的关系，建立了偶极子气动声源辐射声压与流场压力梯度的关系式，并进一步建立了四极子气动声源辐射声压与流场脉动速度的关系式。

（3）结合圆柱绕流结果可以看出，偶极子声源和四极子声源最大值均出现在分离区及其附近，前者更靠近分离区壁面附近流场，而后者稍远离分离区壁面，位于最大偶极子声源后面。

作者贡献声明：

王毅刚：声源识别公式推导。

朱朗贤：结合声源识别方法开展数值仿真计算。

焦燕：参与公式推导和数值仿真计算。

张昊：参与数值仿真计算及实验研究。