基于大涡模拟和FW-H方程的二维圆柱风噪声模拟*

沈国辉，张扬 ，姜浩，郑翀

(1.浙江大学 建筑工程学院，浙江 杭州 310058；2.国家电网温州供电公司，浙江 温州 325000)

基于大涡模拟和FW-H方程的二维圆柱风噪声模拟*

沈国辉1†，张扬1，姜浩2，郑翀2

(1.浙江大学 建筑工程学院，浙江 杭州 310058；2.国家电网温州供电公司，浙江 温州 325000)

针对不同直径二维圆柱的风噪声，采用大涡模拟和FW-H方程声类比的方法进行计算，将基准模型的计算结果与他人计算和实验结果进行对比，分析不同直径二维圆柱的远场声场特征.研究表明：基准模型的数值模拟结果与实验值非常接近，说明了本文计算方法的适用性；远场采样点的声压级频谱峰值频率随着圆柱直径的增大而减小，同根据斯特罗哈数为0.2的理论值计算的峰值涡脱频率接近；10～38 mm直径工况远场采样点的总声压级随着直径的增大基本上呈线性增加趋势；总声压级和A计权总声压级随雷诺数的增加，呈现在亚临界区内增加、而在接近临界区时减小的变化趋势.

气动噪声；二维圆柱；数值模拟；大涡模拟；声类比

Abstract：A numerical simulation using large eddy simulation and FW-H equation acoustic analogy was performed to analyze the acoustic field of 2-D circular cylinders with different diameters.The calculating results of the standard model were compared with the experimental results and the numerical results of other researchers.The far-field aeolian noise characteristics of 2-D circular cylinders with different diameters were analyzed.The Results show that the simulation results of the standard model are very close to the experiment results，indicating the rationality of the simulation.The peak frequency of sound pressure level of the far-field receiver decreases with the increase of the diameters of 2-D circular cylinders，and is close to the peak frequency calculated from the Strouhal number 0.2，which is considered to be the theoretical value.The overall sound pressure level of the far-field receiver under conditions of 2-D circular cylinders with diameters of 10～38 mm increases almost linearly with the diameter.With the increase of Reynolds number，the overall sound pressure level and a-weighted overall sound pressure level have an upward trend in the sub-critical range and a downward trend near the super-critical range.

Keywords：acoustese noise；2-D circular cylinder；numerical simulation；large eddy simulation；acoustic analogy

圆柱绕流风噪声在工程领域中大量存在，如飞机起落架、列车受电弓、高空输电线风噪声等均为圆柱或类圆柱绕流风噪声，属于气动噪声的一种.风经过圆柱时会在圆柱后方产生漩涡脱落，导致固体表面附近流体产生脉动力及湍流应力，向周围空气产生压力波动，从而形成气动噪声[1].

对于圆柱或类圆柱绕流流场特性已有相关研究[2-3]，而对于声场的研究则相对较少.对于圆柱风噪声，已有相关风洞试验研究，Fujita等[4]利用声学风洞试验研究了圆柱在雷诺数为2e5～2e6间的风噪声，King等[5]研究了圆柱在内的不同截面柱体的三维风噪声水平，Alomar等[6]研究了具有不同表面粗糙度的圆柱风噪声.鉴于声学风洞试验的高昂代价，通常采用数值模拟的手段研究圆柱风噪声问题.Brentner等[7]运用非稳态雷诺平均N-S方程进行了圆截面绕流模拟；Cox等[8]通过SST模型和声类比混合方法进行了二维圆柱绕流气动噪声模拟；龙双丽等[9]和Yang等[10]通过二维大涡模拟和声类比混合方法进行圆柱气动噪声模拟.但以上的数值计算大多采用10%的阻塞比进行计算导致结果与试验值吻合较差，同时也没有进行过不同直径二维圆柱的声场研究，亦没有进行过远处声场的反映人耳听觉特性的A计权噪声特性研究.

基于以上背景，本文针对二维圆柱的流场和声场，采用大涡模拟(LES)和FW-H方程声类比的方法进行数值模拟计算，将基准模型的计算结果与实验、其他学者结果进行对比以校验本文计算方法的合理性，分析不同直径二维圆柱远场的声压级频谱和总声压级特征，计算结果为实际工程类似截面构件风噪声水平提供参考.

1 基于LES和FW-H的计算方法

1.1 大涡模拟(LES)方法

流场的模拟采用大涡模拟进行，该方法介于直接数值模拟(DNS)和雷诺平均N-S方程(RANS)方法之间.基于空间平均通过滤波N-S方程将湍流分解为可解尺度湍流和不可解尺度湍流，其最大尺度可与平均运动的特征长度相比，最小尺度则取决于流体的粘性耗散.通过数值求解运动微分方程捕捉可解大尺度涡的运动规律，建立亚格子尺度应力模型模拟小尺度涡的湍动能及对大尺度涡的影响.该方法可以很好地模拟流场的脉动细节特性，显著降低了计算流体力学中使用直接数值模拟(DNS)时的运算量，同时获得了比非稳态雷诺平均N-S方程(URANS)如κ-ε方程等更精确的数值结果.

1.2 基于FW-H方程的声类比方法

基于FW-H方程的声类比方法主要利用CFD软件和声学软件对气动噪声进行联合仿真计算，本质是莱特希尔声类比方法.莱特希尔方程由N-S方程推导而来，经Ffowcs-Williams与Hawkings通过引入Heaviside函数H(f)及其导数Dirac函数δ(f)推导出了FW-H声学方程：

(1)

式中：c为远场声速；p为声压；t为时间；ρ∞为远场流体密度；u和v分别为流场速度和固体声源表面处速度；ni为沿声源表面外法线的单位法向量；x为选定的声场坐标轴方向；Pij为压缩应力张量；Tij为莱特希尔应力张量.下标n，i，j分别表示沿固体声源外法线方向和选定的声场坐标系i，j方向.

莱特希尔声类比法的特点是将声音的产生和传播分别计算，相比直接计算声学法计算量大幅减少.求解FW-H方程可以得到远场采样点处的声压级频谱和总声压级.基于LES和FW-H的方法流程如图1所示，由流体动力计算和声学计算两部分组成.

图1 LES及FW-H方程计算方法的流程

2 流场和声场的物理量计算

对于二维圆柱，常用于分析的参数为阻力系数CD和升力系数CL.计算公式分别为：

(2)

(3)

式中：FD和FL分别为圆柱所受到的阻力和升力；ρ为空气密度；V为风速；D为截面直径.此外，斯特罗哈数St也是一个重要的无量纲数，表征圆柱尾流涡脱落的频率特性，其定义为：

St=fsD/V

(4)

式中:fs表示涡脱的卓越频率.

(5)

式中：r为声源到测点的距离；a0为声速.根据风向可将脉动作用力分解为阻力FD和升力FL，而阻力的脉动远小于升力的脉动，且通常在垂直于风向的平面附近测量风噪声以减少风的对流影响，此时圆柱风噪声主要由升力控制.故脉动作用力Fi可表示为：

(6)

式中:L为圆柱长度.将式(6)代入式(5)，统计一定时间内的声压信息，最终可预测采样点处等效声压满足[4]：

(7)

式中:lc为声源相关长度；CLrms为升力系数有效值；λ为辐射声波长，括号中平方项反映了近场影响.

通过采样点处的声压时程可获得各频率的声压级(Sound Pressure Level，SPL)、频谱特性等.其中声压级的定义如下：

(8)

式中:p为采样点等效声压，pref为参考声压，即人的听觉下限声压，为2e-5 Pa.

通过对各个倍频段下的声压级进行累加，可得采样点处的总声压级OASPL，计算公式为：

(9)

式中：fmin为人的听觉下限频率20 Hz；fmax为人的听觉上限频率20 000 Hz.

人耳对声音强弱的主观感觉，不仅同声压级有关，还和频率和波形有关，例如总声压级为80 dB的1 000 Hz纯音在人耳听来就远比80 dB的50 Hz纯音响，但两者的声压相同.因此工程上根据等响曲线常用A，B，C和D 4种记权特性来表示[11]，即对相应频段声音信号进行折减以更好地反映人耳的听觉特性，如图2所示.其中A，B和C记权分别反映人耳对低频声压、中频声压和高频声压的响度感觉，D记权专用于飞机噪声的测量.研究表明A记权能较好地反映人耳对噪声的主观感觉.

图2 A，B，C，D计权特性曲线

3 基准模型的计算和验证

3.1 基准模型的建立

二维圆截面杆件直径D=1.9 cm在风速V=69.2 m/s(Re=89 000)时的流场和声场被广泛研究[7-10,12]，本文也针对该案例进行计算，称之为“基准模型”.“基准模型”的二维建模计算域和网格划分如图3所示，上游区域取10D(至截面圆心，下同)，下游尾流区取20D，以满足湍流充分发展的需要.上下两边界各距离截面圆心为10D，阻塞比为5%.整个流域划分为9个网格区域，上层和下层的3个区域网格对称划分.圆柱壁面附近网格划分进行加密处理，如图4所示.

图3 二维流场和网格划分

验证壁面网格划分情况的一个重要参数为壁面Y+，即经过无量纲化后表征首层壁面网格高度的参数，定义为：

(10)

式中：y为第一层网格形心高度；Ut为壁面摩擦速度；μ为流体粘度.

流场左侧为速度入流边界，右侧为压力出口边界，上下两边界设定为对称边界，圆截面为无滑移固壁边界.流场计算采用非稳态大涡模拟方法.使用Smagorinsky-Lilly亚格子模型，模型常数Cs=0.1.压力-速度耦合采用PISO算法，空间离散对单元项采用基于节点形式，对流项采用有界中心差分格式，瞬态方程取为二阶精度.时间步长取为1e-6 s.根据文献[8]的研究结果，本文计算中的声源相关长度取5D.

3.2 网格划分的参数分析

为研究近壁面网格的划分对计算结果的影响，需要进行网格划分参数分析.针对“基准模型”，采用3种第一层网格厚度(0.002 mm，0.005 mm和0.01 mm)进行计算.流场和声场的计算结果如表1所示，其中采样点位于圆心正上方128D(2.432 m)处.由表1可知，第一层网格厚度对流场和声场的结果有一定的影响，0.005 mm和0.01 mm的结果同Revell的实验值[12]相比，较0.002 mm的结果更接近.此外，0.005 mm结果同Holmes[13]根据大量实验数据总结的流场参数也接近，而采用0.01 mm时壁面网格已大多难以满足大涡模拟对Y+<1的网格划分要求，经综合考虑，本文的计算均采用0.005 mm网格划分形式.

表1 不同网格划分时标准模型的流场结果

3.3 流场的模拟结果及分析

图5为入流速度V=69.2 m/s，D=1.9 cm，壁面第一层网格厚度0.005 mm时的Y+分布情况，基本满足大涡模拟对壁面Y+<1的要求.

图5 圆柱截面壁面Y+分布

通过瞬态大涡模拟进行非稳态流场计算，图6为不同时间点的瞬态涡量图，其中深色区域显示为较大尺度涡，在尾流区可见上下交替脱落的漩涡，交替的漩涡使得圆柱截面的阻力系数和升力系数时程产生周期性波动.此外，图中还可观察到近圆柱的尾流区出现一些小尺度涡，使得尾流区并未形成低雷诺数下非常明显的卡门涡街形式，这些小尺度涡导致阻力和升力系数时程出现非卓越频率波动.

(a)t=0.094 5 s

(b)t=0.095 0 s 图6 流场涡量图

图7和图8给出了二维圆柱的阻力系数CD和升力系数CL的时程，可以发现两者的变化幅度较大，阻力系数在0.5～2.2之间波动，升力系数约在-2.2～2.0之间波动，均方根值为1.1，这种波动是尾流区上下交替脱落漩涡的直接反映.将升力系数时程经傅里叶变换得到其震荡的卓越频率为732 Hz，根据式(4)换算为斯特罗哈数St为0.201.

图7 圆柱截面的阻力系数

图8 圆柱截面的升力系数

将本文的模拟结果与其他研究结果进行比较，如表2所示，表中所有结果均采用相同的圆柱尺寸和来流速度.由表2可见：1)本文的模拟结果与实验值更加接近，CD比实验值小2.3%，St比实验值大7.4%，但与Holmes总结取值几乎相同；2)采用大涡模拟的各组计算结果比采用其他模型的计算结果更接近实验值；3)在采用大涡模拟的各组结果中，本文的计算结果较其他计算结果更接近实验值，表明本文的计算模型适用于圆柱风噪声的模拟.其原因主要为：本文的计算流域更大，本文的计算流域为上下高度20D，堵塞比为5%，而其他结果[9-10]的计算流域为上下高度10D，堵塞比为10%.

表2 标准模型的CD及St对比

3.4 声场的模拟结果及分析

采样点位于圆心正上方128D(2.432 m)处，将采样点声压时程转化为功率谱，可以发现最大声压级对应频率为752 Hz，换算为St数为0.205，同升力系数对应的St数(0.201)几乎一致，可见圆柱的气动噪声和横风向升力与漩涡脱落息息相关，其特征频率均与漩涡脱落的卓越频率非常接近.

采样点处声压级频谱如图9所示.图中还给出了Revell等[12]的实验结果，该实验结果采用与本文相同的计算工况和采样位置.由图9可知：1)本文的计算结果与实验值在频谱图上分布比较接近；2)峰值对应的St为0.205，略大于实验值0.187；3)计算的总声压级OASPL为103.2 dB，而实验结果为100 dB，非常接近；4)说明本文数值模拟的声学结果同实验值接近，计算结果可信.

图9 采样点处声压级的频谱及比较

同时从图9中还可以发现，在St<0.18时计算的声压级小于实验值，而在St>0.22时计算的声压级大于实验值.其原因可解释如下：1)实验结果包含了风洞本身的背景噪声，使得实验结果在低频段数据显著增加，从而导致实验的低频测量结果偏大,同文献[7-8]所得结果类似；2)二维声学模拟相比三维结果获得的漩涡脱落更规则，而三维结果由于存在沿展长方向的三维效应导致漩涡更混乱，因而反映到频谱则表现为实验的高频区声压级比二维计算结果小，同样结论见文献[7-8].

由式(7)可知，采样点等效声压的平方与V6成正比，代入相关参数数值可计算圆心正上方128D处的总声压级为112.7 dB，比试验值(100 dB)及模拟值(103.2 dB)偏大，相似的对比结果见文献[4].造成该结果的原因主要为：公式值只考虑了偶极子噪声的变化规律(等效声压的平方同风速6次方正相关)，实际的结果还受到其他噪声如单极子和四极子等的影响，故一般认为有效声压的平方与风速的4～6次方之间正相关[6].本例中取V5.5时式(7)的计算值为103.4 dB，同试验和模拟值接近.

4 不同直径二维圆柱的声场计算

分别进行10 mm，19 mm，31 mm，38 mm，45mm，60 mm和70mm直径二维圆柱的气动噪声数值模拟.不同直径模型的网格划分均同“基准模型”，流场大小均按照上下左右边界距离圆心为10D，10D，10D和20D进行布置.所有工况均采用大涡模拟进行流场模拟，入流速度V=30 m/s.声学计算中选取声源相关长度均为5D，采样点位于圆心正下方100 m处.

流场进入收敛阶段后，计算圆柱周围声场，获得采样点处声压时程，通过傅里叶变换进行频域分析，可获得各工况噪声的卓越频率.通常亚临界范围内圆柱绕流的St数可取0.2(本文称为理论值)[13]，可以通过式(4)反算不同圆柱直径下的涡脱卓越频率fs，将本文的模拟值和理论值绘于图10.由图10可知，采样点处的噪声卓越频率模拟值与理论值吻合较好，均随着直径增长而成反比例降低.

图10 各工况噪声卓越频率

对所得噪声频谱进行可听域(20～20 000 Hz)内1/3倍频程处理，图11为7种工况1/3倍频程未计权声压级频谱，由图可知：1)峰值声压级对应1/3倍频带随圆柱直径增加而减小，并出现明显的频带尖峰；2)在高频段(1 000 Hz以上)各工况声压级的差异小于低频段，表明不同直径圆柱表面较小尺度的涡脱落程度差异不大.

图11 1/3倍频程声压级频谱

对图11经A计权计算可获得各工况风噪声的A计权声压级，如图12所示.由图12可知：1)经过A计权后，200 Hz以下频段的声压级折减严重，30 Hz以下所有工况的A计权声压级均为0 dB；2)10 mm直径工况的低频段A声压级最小，200 Hz以内频段均为0 dB；3)38～70 mm直径工况声压级峰值被消减，导致了A记权下噪声的1/3倍频程卓越频带升高，其中60 mm直径工况由约100 Hz变为约1 600 Hz，同未计权结果差异巨大.

图12 1/3倍频程A计权声压级频谱

通过式(9)累计各1/3倍频段声压可获得总声压级OASPL，结果如图13所示，其中横坐标已转化为各直径圆柱对应的雷诺数.由图13可知：

1)从A计权和不计权总声压级结果对比看，两者的差异随着圆柱直径的增加而增大，圆柱直径10 mm时差距约为1 dB，圆柱直径70 mm时差距接近6 dB，其主要原因为随着圆柱直径增大，噪声的峰值频段向低频段迁移，由图2可知A计权的折减程度增加.

2)从不计权总声压级看，10～38 mm直径工况(对应雷诺数为2.07 e4～7.87 e4)的总声压级几乎成线性增加，38 mm直径至60 mm直径(对应雷诺数为7.87 e4～1.24 e5)增速放缓，到70 mm直径时(对应雷诺数为1.45 e5)声压级出现了下降，较前一工况降低了约4 dB，甚至小于31 mm直径工况(对应雷诺数为6.42 e4)，其原因为在较大直径时流动开始向雷诺数临界区域发展，在直径70 mm时(对应雷诺数为1.45 e5)已接近临界区域，其涡脱情况复杂，在噪声中表现为可听域频带噪声较前一工况增长不显著，甚至出现噪声减弱的情况.

3)从A计权总声压级看，其发展趋势基本和不计权声压级一致，在38 mm直径后(对应雷诺数为7.87 e4)出现A计权总声压级略微减小的情况，70 mm直径工况显著减小.

4)A计权和不计权总声压级随雷诺数的增大先增大后减小，可以预见随着圆柱直径的进一步增加，雷诺数接近临界区后两种总声压级会出现降低的趋势，文献[4]也发现在雷诺数位于(2 e5～5 e5)时未计权总声压级随着雷诺数增大而减小的趋势，与本文模拟趋势一致.

图13 各工况的总声压级和A计权总声压级

5 结 论

本文基于LES和FW-H方程方法研究不同直径二维圆柱的气动噪声特性，主要结论有：

1)本文针对“基准模型”的数值模拟获得的气动力、声学参数与实验值和经验值非常接近，说明本文计算方法适用于求解二维圆柱的气动噪声问题.

2)在采用大涡模拟计算获得二维圆柱的气动力参数中，本文的计算结果比其他结果更接近于实验值，其原因为本文采用更大的计算域(本文堵塞比为5%，其他研究为10%).

3)由7种直径二维圆柱远场采样点的声压级频谱可知：峰值频带随着圆柱直径的增大而减小，并出现明显的频带尖峰，与根据斯特罗哈数理论值0.2计算的对应卓越频率接近.

4)从总声压级角度看，10～38 mm直径二维圆柱远场采样点的总声压级随着直径的增大而基本呈线性增加，此后总声压级增速放缓，至70 mm时声压级减小，其原因为较大直径工况的雷诺数已开始接近临界区域.总声压级和A计权总声压级随雷诺数增大在亚临界区内先增大、接近临界区时出现减小趋势.

[1] LIGHTHILL J.Waves in fluids[M].Cambridge：Cambridge University Press，1978:23-40.

[2] 王汉封,邹超,王启文,等.均匀流中悬臂圆柱体气动力雷诺数效应[J].湖南大学学报：自然科学版,2015,42(5):65-71.

WANG Hanfeng，ZOU Chao，WANG Qiwen，etal.Reynolds number effects on the aerodynamic forces of a cantilevered circular cylinder in uniform flow[J].Journal of Hunan University：Natural Sciences，2015,42(5):65-71.(In Chinese)

[3] 王述良,梁枢果,邹良浩,等.基于刚性节段模型风洞试验的输电导线阻力系数研究[J].湖南大学学报：自然科学版,2016,43(3):32-40.

WANG Shuliang，LIANG Shuguo，ZOU Lianghao，etal. Study on drag coefficients of conductors based on wind tunnel tests of rigid sectional model[J].Journal of Hunan University：Natural Sciences，2016,43(3):32-40.(In Chinese)

[4] FUJITA H.The characteristics of the Aeolian tone radiated from two-dimensional cylinders[J].Fluid Dynamics Research，2010，42(1):154-168.

[5] KING W F，PFIZENMAIER E.An experimental study of sound generated by flows around cylinders of different cross-section[J].Journal of Sound and Vibration，2009，328:318-337.

[6] ALOMAR A，ANGLAND D，ZHANG X,etal. Experimental study of noise emitted by circular cylinders with large roughness[J].Journal of Sound and Vibration，2014，333:6474-6497.

[7] BRENTNER K S，COX J S，RUMSEY C L,etal.Computation of sound generated by flow over a circular cylinder：an acoustic analogy approach[R].Hampton：NASA Langley Research Center，1996:1-7.

[8] COX J S，BRENTNER K S，RUMSEY C L.Computation of vortex shedding and radiated sound for a circular cylinder：subcritical to transcritical Reynolds numbers[J].Theoretical and Computational Fluid Dynamics，1998(12)：233-253.

[9] 龙双丽,聂宏,许鑫.不同雷诺数下圆柱绕流气动噪声数值模拟[J].声学技术，2011，30(2):111-116.

LONG Shuangli，NIE Hong，XU Xin.Numerical simulation of noise induced by flow around a cylinder at different Reynolds number[J].Technical Acoustics，2011，30(2):111-116.(In Chinese)

[10] YANG D，LI J，LIU J，etal.Analysis on physical mechanism of sound generation inside cavities based on acoustic analogy method[J].Open Journal of Fluid Dynamics，2013，3：23-31.

[11] 马大猷.噪声与振动控制工程手册[M].北京：机械工业出版社，2002：718-725.

MA Dayou.Handbook of noise and vibration control engineering [M].Beijing：China Machine Press，2002：718-725.(In Chinese)

[12] REVELL J D，PRYDZ R A，HAYS A P.Experimental study of airframe noise vs drag relationship for circular cylinders[R].Burbank：Lockheed-California Corporation，1977.

[13] HOLMES J D.Windloading of structures [M].London：Spon Press，2003：83-88.

Numerical Simulation on Aeolian Noise for 2-D Circular Cylinders Based on Large-eddy Simulation and FW-H Equation

SHEN Guohui1†，ZHANG Yang1，JIANG Hao2，ZHENG Chong2

(1.College of Civil Engineering and Architecture，Zhejiang University，Hangzhou 310058，China； 2.State Grid Corporation Wenzhou Electric Power Company，Wenzhou 325000，China)

1674-2974(2017)09-0033-08

10.16339/j.cnki.hdxbzkb.2017.09.005

2016-07-21

国家自然科学基金资助项目(51178425), National Natural Science Foundation of China(51178425)；浙江省自然科学基金资助项目(Y16E080016), Natural Science Foundation of Zhejiang Province of China (Y16E080016)

沈国辉(1977—)，男，浙江台州人，浙江大学副教授，博士

†通讯联系人,E-mail: ghshen@zju.edu.cn

TU312

A