基于K-FWH声比拟方法的串列双圆柱气动噪声研究

陈武，周毅

(南京理工大学 能源与动力工程学院，南京 210094）

飞机的气动噪声问题是航空航天工程领域中不可忽视的研究内容之一。飞机的气动噪声主要分为发动机噪声和机体噪声2类。随着高新技术的发展与应用，发动机噪声被大大降低，因此机体噪声成为飞机噪声的主要噪声［1］。在机体噪声中，飞机起落架和增升装置是主要的气动噪声声源。Chow等［2］对空客A340进行了实验研究，结果发现起落架噪声比襟翼噪声高6 dB。因此，研究飞机起落架的气动噪声对于飞机降噪问题具有重要意义。由于飞机起落架的结构十分复杂，在研究中通常将其简化为串列双圆柱模型［3］。相比较于单圆柱绕流，串列双圆柱的流场更加复杂，流场中不仅存在流体分离、涡脱落过程，还存在上下游圆柱间的流动干扰等非定常流动现象。上下游圆柱产生涡脱落导致其受力发生剧烈变化，从而引起气动噪声。因此，对串列双圆柱绕流气动噪声特性的研究具有重要的学术价值和实际的工程意义。

国内外学者对串列双圆柱绕流问题做了大量的实验和数值模拟研究。美国国家航空航天局的Langley研究中心［4-6］在QFF和BART，风洞中进行了大量的空气动力学实验和声学实验，这些实验研究成为后来研究双圆柱绕流噪声的空气动力学特性和气动噪声特性的基石，大量的数值模拟研究都以此为标准进行后续研究。随着计算机性能的提升和并行数值方法的发展，计算流体力学(CFD）结合FWH(Ffowcs Williams Hawkings）方程的声比拟方法渐渐成熟，成为当前计算流体气动噪声最有应用前景的方法之一［7］。

国外研究人员对串列双圆柱噪声问题进行了一系列数值仿真研究。Lockard等［3］应用基于有限体积法的三维求解器CFL3D，对串列双圆柱绕流进行了数值模拟计算，结合FWH方程计算了远场噪声并与实验结果进行了对比。Brès等［8］使用离散玻尔兹曼方法与FWH方程相结合的方法，对串联圆柱体的远场气动噪声进行了预测。Papaioannou等［9］比较了串列双圆柱的二维模拟与三维模拟的区别，结果表明，三维模拟能更好地模拟尾迹的涡量场、圆柱表面受力及涡脱落频率，得到了与实验更接近的临界间距比。

与此同时，国内学者对串列双圆柱噪声问题也开展了相关研究。刘敏等［10］对不同间距比的串列圆柱扰流问题进行了大涡模拟，并采用Farassat-1A方程计算了相应的载荷噪声。赵良举等［11］对二维横掠串列双圆柱绕流噪声问题进行大涡模拟，研究了不同流速、直径、间距比对绕流流场及其产生的气动噪声的影响。龙双丽等［12］采用大涡模拟与声比拟相结合的混合方法，预测了不同雷诺数下二维圆柱绕流的远场气动噪声。余雷等［13］采用基于非线性k-ε湍流模型的限制数值尺度(Lim ited Numerical Scales，LNS）方法对双圆柱绕流进行了数值模拟，将计算结果与相应的实验结果进行了对比分析，结果两者吻合良好。宁方立等［14］采用大涡模拟与FWH方程相结合的方法对串列双圆柱绕流进行三维数值仿真，并对下游不同直径的情况进行对比研究。高威和陈国勇［15］采用改进型延迟分离涡模拟方法对串列双圆柱的绕流流场进行数值模拟，比较了不同来流速度、圆柱间距比和圆柱直径的双圆柱绕流的气动噪声特性。有些学者采用其他方法对串列双圆柱进行研究。周凯等［16］基于格子玻尔兹曼方法，对二维静止串列双圆柱绕流进行了数值模拟，并将模拟结果与已有研究结果进行了对比分析。Du等［17］采用流声分解方法对低雷诺数(Re＝200）下的串列和并列双圆柱绕流噪声问题进行了研究，着重分析了不同圆心距和不同流动形式下的串列、并列双圆柱绕流所引起的声场分布。此外，还有学者对飞机中的高升力体和工业生产中的水泵进行了研究［18-21］。葛明明等［18］发展了高精度隐式大涡模拟模型，并结合FWH声比拟方法对30P30N高升力体气动噪声问题进行了数值模拟，得到了更接近实验值的模拟结果，准确预测了缝翼低频窄带噪声，并得到了合理的噪声辐射指向性分布。蔡晓彤等［20］基于CFD理论与Lighthill声比拟理论，对潜水排污泵的流场和声场进行数值模拟，并对不同工况下潜水排污泵的内部压力分布特性进行分析，同时探讨了内声场和外声场的噪声产生原因及分布传播特性，结果表明，采用SST模型得到的性能曲线最接近试验结果，隔舌处是主要噪声源。

串列圆柱气动噪声问题属于典型的气动噪声测试验证问题，前人的研究大多集中在阐明不同间距比、不同来流速度及上下游不同圆柱直径的影响，而对上下游柱体涡脱落特性(阻力系数、升力系数、相位变化）及其对远场噪声的影响这一重要问题研究较少。本文对某一标准算例开展深入研究，利用相关性分析和时频分析对升阻力和噪声特性物理关联进行进一步解析。本文采用大涡模拟和基于Kirchhoff积分的K-FWH 方程［22］相结合的方法，对串列双圆柱绕流流场及其声场进行三维数值模拟研究。首先，通过大涡模拟求解不可压Navier-Stokes方程得到非定常流场，即近场声源。然后，以圆柱表面为积分面，结合K-FWH方程时域解进行积分求解，得到观测点的声压时域值，即远场辐射噪声。在此基础上，验证了标准算例的气动特性和噪声特性，证实了流动数值模拟方法及噪声简化算法的精确性。其次，在时域和频域对上下游圆柱的升阻力系数进行分析，并对上下游圆柱升阻力系数进行相关性分析。最后，对上下游圆柱产生的声压及其分解项进行相应的分析。

1 数值方法

1.1 流场计算

本文采用大涡模拟方法，利用壁面自适应局部涡黏(Wall-Adapting Local Eddy-Viscosity）亚格子模型，模拟其对大尺度流动的影响。WALE湍流模型是在湍流结构的运动和动力学性质的基础上，将转动张量包含在模型中构造而成的，优点是：该模型能正确反映出近壁区域涡黏系数与垂直壁面距离的三次方成正比的性质，且不含关于边界几何尺寸的参数，易于应用于复杂的湍流模拟。该模型表达式如下［23］：

1.2 声场计算

通过引入广义函数，按照波动方程的形式，将Navier-Stokes方程重新整理成FWH方程［24］，FWH方程是声比拟方法中最常用的形式，表达式如下：

式中：F为积分表面；ui为xi方向的速度分量；un为F＝0面上的法向速度；vn为物面速度的法向分量；δ(F）为狄拉克函数；p′为远场瞬时声压；a0为远场声速；H(F）为赫维赛德阶跃函数，当F＝0，H(F）为物体运动边界控制面方程，F＜0表示无扰动流体的空间区域，H(F）＝0，F＞0表示物面外受到流体扰动的空间区域，H(F）＝1；p为当地壁面压力。

Tij为Lighthill应力张量

Pij为表面载荷

式(2）右边第一项为由流体本身内应力造成的四极子源项Q，第二项为由于壁面压力的存在而引起的偶极子源项D，第三项为由当地的流体密度发生变化而引起的单极子源项M。本文为不可压缩流体，故不考虑单极子源项M。本文流动为低亚声速流动，四极子源项Q 相对于偶极子源项D 非常小，可以忽略不计，故只考虑偶极子源项D。本文采用的K-FWH声压时域解公式是在FWH 方程时域解的基础上结合Kirchhoff方法，该公式的主要优点是不需要知道表面压力的法向导数。此外，它也适用于具有可渗透表面的物体，而经典的FWH方程在这种情况下是无效的。这种新的边界积分公式同时也简化了声压的计算，具体参考文献［22］。可以将远场瞬时声压p′的时域解公式写成如下形式：

根据实际情况只考虑偶极子源项D，简化公式可以得到如下表达式：

式(6）右边第一项为压力项Dp，与壁面压力直接相关，第二项为压力的时间导数项Dt，与壁面压力的时间导数相关。下文将对式(6）中的简化假设进行验证。

2 计算模型与方法验证

2.1 计算模型

参照文献［25］的标准算例验证过程，串列双圆柱模型由2个直径相同的圆柱体沿着来流方向串列排布。圆柱的直径D＝0.057 15 m，2个圆柱中心点之间的间距比L／D＝3.7，来流速度Uin＝43.4 m／s，基于圆柱直径和来流速度的雷诺数为Rein＝UinD／ν＝1.66×105，ν为运动黏度。

为了捕捉到重要的流动特征，根据文献［25］的建议，在建模过程中取圆柱的展向长度为3D。选取上游圆柱体的中心作为坐标系原点，计算区域为：－10D≤X≤25D，－10D≤Y≤10D，0≤Z≤3D。对圆柱周围进行加密处理，壁面第1层网格Y＋≈2，X-Y平面上的计算域网格如图1所示。展向设置为周期性边界条件，圆柱表面为无滑移壁面边界条件，其余边界为远场边界条件。空间离散为二阶迎风格式，时间离散为二阶隐式格式，时间步长为Δt＝1×10－5s，下文展示的时间平均结果和噪声结果均在非定常流动呈现相对稳定的周期性特征之后开始统计，采样统计时间为70 000个时间步长(对应总的时间为0.7 s），对应约126个涡脱落周期。

图1 X-Y平面上的计算网格Fig.1 Computational grids of X-Y plane

2.2 数值模拟及噪声计算验证

为了验证网格收敛性，本文采用3种精细程度的网格(算例1、算例2、算例3）进行计算，对应的网格数分别为109万、238万、665万，将计算得到的上下游圆柱的平均阻力系数和涡脱落频率与实验结果进行对比(上标1表示上游圆柱，2表示下游圆柱）。上下游圆柱的气动特性计算结果如表1所示。从阻力计算结果来看，3个算例的上下游圆柱平均阻力系数均在前人模拟［25］的最大值与最小值区间内，涡脱落频率也均在最大最小值范围内。理论上网格的加密会使结果更接近真实情况，验证算例中随着网格密度的增加，结果变化并不明显，且均在合理范围内。表1说明在误差允许范围内，精细程度为算例3的网格已满足计算要求。因此，除特殊标明外，下文展示均为采用算例3(网格数为665万）的计算结果。

表1 上下游圆柱气动特性统计结果Table 1 Statistic results of aerodynamic properties for upstream and downstream cylinders

为了进一步验证流场数值模拟方法的精确性，本文计算得到上下游圆柱表面时均压力系数Cp及时均脉动压力系数Cp′和观测点处(9.11D，32.49D，1.50D）的噪声功率谱密度(PSD），并与实验数据［25］和前人数值模拟数据［13］进行对比。如图2所示，上下游圆柱的表面时均压力系数均与实验结果十分接近，其中，A为沿着圆柱一圈的度数，一圈为360°。

图2 上下游圆柱表面时均压力系数Cp分布Fig.2 Distribution of time-averaged pressure coefficient Cp on the surface of upstream and downstream cylinders

如图3所示，计算得到的上游圆柱时均脉动压力系数Cp′与实验数据趋势相同但高于实验值，但与前人模拟数据［13］接近，下游圆柱Cp′与实验数据［25］吻合较好，上游圆柱Cp′取值较小，因此前人模拟结果与本文模拟结果均与实验结果存在误差，但误差在可接受范围内。

图3 上下游圆柱表面时均脉动压力系数Cp′分布Fig.3 Distribution of rootmean square Cp′of pressure coefficient disturbance on the surface of upstream and downstream cylinders

图4所示为观测点(9.11D，32.49D，1.50D）处的噪声功率谱密度，表明噪声功率谱密度的主峰频率和峰值、次级峰频率和峰值及低频分布均与实验相吻合。由于采用的数值模拟方法为大涡模拟，很难准确捕捉流场中小尺度湍流运动。而噪声功率谱密度的高频分布与流场的小尺度涡运动密切相关，因此图4中噪声功率谱密度的高频分量略低于实验数值。图4为算例2和算例3(网格数分别为238万和665万）所预测的声压级功率谱密度，可以看出两者存在细微差别，在主峰频率上相差约3%，不同频率的功率谱密度分布一致。图4进一步说明了利用本文方法及网格布置能够准确描绘本文重点关注对象——远场噪声特性。

图4 观测点处的声压级噪声功率谱密度Fig.4 Power spectral density of sound pressure level at observation point

图5给出了位于r／D＝33.74处的声压级指向分布图及观测点(9.11D，32.49D，1.50D）的主频峰值大小，其中r是与观测点相同X-Y平面上离坐标原点的距离。声压级由噪声功率谱密度主频峰值计算得到。图中显示声压级成8字形分布，说明了计算结果具有显著的偶极子声源指向特性，也验证了本文噪声辐射指向性分布的合理性。图5表明，算例3(网格数为665万）与算例2(网格数为238万）基本指向性相同，但在0°和180°这2个方向存在着约6%的误差，算例3的偶极子特性相比较于算例2更加明显，但对其他方向的指向性分布影响较小。

图5 远场声压级指向分布图(r／D＝33.74）Fig.5 Directivity of Sound Pressure Level at r／D＝33.74

3 计算结果与分析

在上述验证工作的基础上，对Uin＝43.4 m／s，圆柱直径为D，间距比L／D＝3.7的上下圆柱涡脱落流场特性及声场特性进一步研究。

利用速度梯度张量第二不变量Q对流场进行后处理计算，如图6所示，其为Q值为1 000的瞬时等值面图。从图中可以看出，在上下游圆柱后方流体均出现回流(图6中蓝色部分）。流体流经上游圆柱，形成了涡脱落，从上游脱落的涡再作用在下游圆柱上。因而下游圆柱后部有着更为复杂的涡产生机制。另外，可以看见有流向肋条状涡结构附着在下游圆柱外表面上。

图6 Q值为1 000的瞬时等值面图Fig.6 Instantaneous isosurfaces with Q＝1 000 and corresponding stream wise velocity

图6表明，上游圆柱产生的涡结构脱落直接作用在下游圆柱表面上，并在下游形成新的涡脱落。为了进一步研究上下柱体间涡脱落关联，本文基于皮尔逊相关性系数R对上下游圆柱的升阻力系数进行相关性统计。相关性系数R的取值范围为－1～1，绝对值越大代表相关性越大。研究结果表明，上下游圆柱阻力系数间的相关系数Rd为－0.343 7，升力系数间的相关系数Rl为－0.580 5，说明上下游圆柱的升阻力在时域上呈现一定的负相关性，即上下游圆柱的大尺度涡呈现一定程度上的反相位脱落。

如图7所示，图7(a）中水平虚线为沿着流向的中心线，图7(b）中垂直虚线与下游圆柱前后壁面相切。由图7可知，上游圆柱前方的高压区明显高于下游圆柱。这是因为2个圆柱之间存在一个显著的负压回流区(速度回流区见图6(a）、(b））。此外，本文中2个圆柱的间距比仅为L／D＝3.7，两柱体间流场并未完全恢复，从而上游圆柱的前后压力差大于下游圆柱，即上游的平均阻力系数大于下游的平均阻力系数。此结果与表1中的结果相同。故可以推测随着间距比的增大，上下游圆柱阻力系数趋于相等。

图7 流场平均压力分布Fig.7 Pressure distribution of flow field

图8 上下游圆柱阻力系数及升力系数随时间变化Fig.8 Time history of drag coefficients and lift coefficients on the surface of upstream and downstream cylinders

通过快速傅里叶变换，将时域数据转换为频域，图9为上下游圆柱阻力系数(和）与升力系数(和）频谱。从图中可以看出，虽然上游圆柱涡脱落对下游圆柱涡脱落过程有着显著影响(上下游圆柱升阻力系数成负相关性），但是上下圆柱的涡脱落频率相同，均为187.6 Hz。升力系数的主频与之前验证的噪声功率谱密度中的主频率相同，而阻力系数的主频为2倍的升力系数主频。这是因为一次完整的涡脱落过程包含纵向运动(Y方向）及由此导致的流向运动(X方向）。其中，一个周期的脱落包含2次纵向运动和1次流向运动。因此，阻力系数(流向运动）的主频为升力系数(纵向方向）主频的2倍关系。由图9还可以看到，下游圆柱的阻力系数频谱和升力系数频谱均在上游圆柱的上方，说明下游圆柱导致的扰动在不同频率上均比上游大。可以进一步推测，与圆柱表面压力相关的观测点声压主要是由下游圆柱产生的声压主导。另外，在下游升力系数频谱上，可以观察到次级峰的存在，但阻力系数的次级峰不明显。由此可以推测上游脱落的涡旋主要作用在下游圆柱上下表面(见图6），流向运动影响较小，进而只在下游升力系数产生了次级峰的现象。

图9 上下游圆柱阻力系数及升力系数频谱Fig.9 Frequency spectra of drag coefficients and lift coefficients on the surface of upstream and downstream cylinders

图8(b）表明升力系数时程曲线呈现规则变化，且上下游圆柱的脱落频率相同(f＝187.6 Hz），故可采取希尔伯特变换进一步评估上下游升力系数(和）之间的相位差演化。图10为上下游圆柱升力系数相位差的时间演化图，对应图8(b）中的时间范围，其中的相位差，°，两者相位完全相同，当°，两者相位完全反相。从图10中可以看到，上下游圆柱升力系数的相位差在选定时间范围内由0°～180°变化，且相位差在90°～180°的时间更长。

图10 上下游圆柱升力系数相位差的时间演化Fig.10 Time history of lift coefficient phase difference between upstream and downstream cylinders

如图11概率密度函数(PDF）分布所示，上下游圆柱升力系数的相位差概率密度主要分布在90°～180°，小于60°的概率相对较小，这与图8中初步观察得到的结论相符合。升力系数与流体的涡脱落密切相关，图11得到的结果表明在串列双圆柱绕流中，上下游涡脱落同相位(即涡同时向上或向下脱落）的概率较小，而呈现反向脱落(相位差为90°～180°）的概率相对较大。本节前文中对上下游圆柱升力系数的相关性进行计算，升力系数间的相关系数Rl为－0.580 5(即上下游圆柱的大尺度涡呈现一定程度上的反相位脱落），图11对此结果进一步作出了解释，说明上下游圆柱的涡脱落呈现一定相位差，相位差在90°～180°，且相位差接近180°的概率最大。

图11 上下游圆柱升力系数相位差的概率密度分布Fig.11 Probability distribution of lift coefficient phase difference between upstream and downstream cylinders

根据式(6）可将声压分解为两部分：一部分是基于圆柱壁面瞬时压力的压力项Dp，另一部分是基于其瞬时压力的时间导数项Dt。图12为观测点(9.11D，32.49D，1.50D）处上下游圆柱的声压P及其分解项时程曲线。从图12中可以看出，瞬时压力的时间导数项Dt是主要的声压组成部分，与声压Dp＋Dt基本一致。下游圆柱产生的声压值明显大于上游圆柱，可以从时域上判断下游圆柱为总声压的主要产生者，从而验证了之前的推测。

图12 上下游圆柱声压及分解项随时间变化Fig.12 Time history of acoustic pressure and its decompositions on the surface of upstream and downstream cylinders

为了进一步揭示串列双圆柱绕流声场特性，通过快速傅里叶变换，将声压时域信号转换为频域信号。图13为上下游圆柱声压分解项的噪声功率谱密度。如图13所示，不论是Dp还是Dt，下游圆柱的功率谱密度均大于上游圆柱。由于总声压是上下游圆柱产生的声压的叠加，可以证明，下游圆柱产生的能量在串列双圆柱绕流噪声中占主导地位。图13中时间导数项Dt在整个频段上占主要地位，而压力项Dp则主要与低频分布有关。由此可见，时间导数项Dt有更加剧烈的波动，产生更高频率的能量。从图13中可以观察到下游的压力项和时间导数项均出现次级峰的现象。进一步说明了上游圆柱脱落的涡直接作用在下游圆柱上，进而导致了下游圆柱频谱产生次级峰的现象。

图13 上下游圆柱声压分解项的声压级噪声功率谱密度Fig.13 Sound pressure level noise power spectrum of sound pressure decomposition on the surface of upstream and downstream cylinders

基于皮尔逊相关性系数R对偶极子源项D及其分解项进行了相关性分析。表2为上下游圆柱声压(D1，D2）及其分解项之间的相关性系数统计。从表中可以看到，上下游各自的声压D1、D2与其时间导数项、的相关性系数R均在0.98以上，进一步说明Dt为D 的主要组成部分。这与前文的结论相符。压力项Dp与声压D 的相关性系数R为0.15左右，说明Dp与D 有一定的相关性。另外，本文计算了上下游圆柱对于选定观测点的总声压Da与上下游声压的相关性系数R。总声压与下游声压的相关性系数为0.936 5，与上游声压仅为－0.149 4，验证了在此观测点处下游声压在总声压中占主导地位。

表2 相关项之间的相关系数Table 2 Correlation coefficients between related term s

从图12中可以看出，上下游圆柱产生的声压在时域上呈现一定的规则变化，且从图13中得到上下游声压主频相同(f＝187.6 Hz），故也采用希尔伯特变换进一步评估上下游声压(D1和D2）之间的相位差演化。图14为上下游圆柱声压相位差的时间演化，对应图12(b）中的时间范围，其 中为D1和D2相位差，≤180°。当°，两者相位完全相同，当＝180°，两者完全反相。

图14 上下游圆柱声压相位差的时间演化Fig.14 Time history of acoustic pressure phase difference between upstream and downstream cylinders

从图14中可以看出，在相同演化时间段中，图14的相位差随时间演化规律与图10十分类似，上下游圆柱声压相位差随着升力系数相位差同步变化，两者变化曲线大致相同，同样可以看出上下游圆柱声压相位差在90°～180°的时间分布更长。本文进一步统计了上下游圆柱声压相位差的概率分布情况，如图15所示。由图15可以得到，声压相位差概率密度分布基本与升力系数相位差概率密度分布一致，上下游声压的相位差分布在90°～180°的概率较大。此结果说明了上下游圆柱产生的声压与各自圆柱的升力系数(即大尺度涡脱落特性）密切相关，相位差的时间演化规律和概率密度分布大致相同。

图15 上下游圆柱声压相位差的概率密度分布Fig.15 Probability distribution of acoustic pressure phase difference between upstream and downstream cylinders

图16 基于上下游圆柱声压相位差的选定观测点声压脉动值条件平均分布Fig.16 Conditional average of acoustic pressure pulsation at selected observation points based on acoustic pressure phase difference between upstream and downstream cylinders

4 结论

1）证实了WALE大涡模拟模型结合基于KFWH方程的声比拟方法能够较好地捕捉到串列双圆柱的流动过程及其噪声特性。上下游圆柱的涡脱落频率相同，但上游圆柱的阻力系数大于下游圆柱。上下游圆柱的升阻力系数为负相关，大尺度涡呈现一定程度的反相位脱落，且下游的升阻力系数比上游波动更加剧烈。

2）串列双圆柱的气动噪声声压主要由两部分组成：一部分为圆柱表面的瞬时压力项Dp，另一部分为瞬时压力的时间导数项Dt。其中瞬时压力的时间导数项Dt是主要的声压组成部分，与总声压Da基本一致。

3）远场观测点的瞬时声压由上下游圆柱各自产生的声压叠加产生，下游圆柱产生的噪声在串列双圆柱绕流噪声中占主导地位。时间导数项Dt与噪声功率谱密度的低频和高频分布均有关，而瞬时压力项Dp主要与低频分布有关。声压的高频能量主要由时间导数项Dt决定。由于上游涡脱落对下游圆柱的涡脱落的影响，导致下游升力系数频谱及观测点总噪声功率谱密度呈现次级峰的现象。

4）上下游圆柱产生的声压与上下游的升力系数(即大尺度涡脱落特性）密切相关，相位差的时间演化规律和概率分布大致相同。观测点的声压脉动值不受上下游声压的相位差影响，即条件平均＞约为一个固定值。

笔者团队前期对规则及分形布置柱体绕流场中涡脱落特性开展了直接数值模拟研究，研究表明，在相同阻塞率下，规则布置中相邻柱体后涡脱落呈现“反向相位锁定”(Antiphase Locking）特性，而分形布置绕流场中存在多尺度的拟序结构，相邻柱体后涡脱落不存在显著相位锁定。本文研究表明绕流场中噪声强度与柱体表面受力情况(柱体涡脱落特性）直接相关，能否利用分形布置柱体实现绕流场中噪声的被动控制是未来值得研究的方向。