宋亚辉，秦 浩，刘秋洪

（1.中国飞行试验研究院，西安 710089； 2.西北工业大学 航空学院，西安 710072）

非紧致边界气动噪声数值预测方法

宋亚辉1，秦 浩1，刘秋洪2

（1.中国飞行试验研究院，西安 710089； 2.西北工业大学 航空学院，西安 710072）

通过寻找满足非紧致边界的精确格林函数，发展了一种基于声模拟理论的非紧致边界气动噪声数值预测方法。该方法首先采用边界积分方法计算满足相应边界条件的精确格林函数，然后采用精确格林函数求解FW-H方程，进行远场气动噪声预测。考虑了非紧致边界对声波传播的散射作用，适用于尺寸较大且几何外形复杂的边界与非定常流动相互作用诱发的气动噪声的数值预测。以雷诺数90 000、马赫数0.2的圆柱绕流诱发的气动噪声为例，数值预测了非紧致边界条件下的气动噪声，并与采用自由空间格林函数求解FW-H的计算结果进行了对比，对该方法进行了验证。

声学；气动噪声；FW-H方程；非紧致边界；数值预测

随着计算机技术的进步和气动声学研究的深入，计算气动声学逐渐发展成为一种探索气动噪声机理、寻找噪声源位置、预测噪声大小的重要工具。计算流体力学(CFD)和声模拟理论[1](一般采用FWH方程[2]和自由空间格林函数)相结合的方法是迄今为止工程应用中最受欢迎的气动噪声预测方法[3]。其基本思想是将声场的计算分解为两步：首先，近场部分采用CFD技术进行数值计算，得到声源信息数据；然后，采用声模拟理论来求解声波从近场到远场部分的传播。声模拟理论由Lighthill[1]提出，经过研究者[2，4]不断的发展，目前被广泛用于气动噪声预测问题研究。Lighthill的声模拟方程为

式中c0为声速，p′是声压，Tij表示Lighthill湍流应力张量，定义为

式(1)解的通用表达式为[5]

式(3)即是时域的FW-H方程。其中，x和y分别表示观察点和源点位置，t和τ分别表示观察点和源点时刻，Ω是产生声源的流动区域，S是物体边界，n是物体边界上指向流动区域的法线，ph为单位面积物体边界作用在流体上的力，ρ0是密度，vn表示物体边界运动的法向速度，G(x,t|y,τ)是格林函数，一般采用自由空间格林函数是Heaviside函数。

FW-H方程实质上是Naiver-Stokes方程的恒等变形，求解气动噪声时，理论上要求声源必须采用满足低色散、低耗散的高精度方法计算流动[6]。而目前非定常流动多采用低精度的有限容积法求解，由于非定常流动激发的气动噪声在向空间传播时会在物体边界上发生反射和散射，采用自由空间格林函数直接求解FW-H方程不能完全捕捉边界对声波传播的影响，因而在流场求解精度不足够高的情况下，要求物体边界需要满足声学紧致条件(即声波波长远大于边界特征尺寸)。对于外形尺寸较大的物体或关注的声波频率较高时，边界不再满足声学紧致条件，直接应用FW-H方程会产生较大误差[7]。

针对非紧致边界问题，大量学者进行了研究，主要有以下几种方法：

(1)基于Lighthill方程变分形式的有限元法[8]，能模拟声波从近场到远场的传播过程，需要对整个求解域进行离散求解，当进行远场噪声计算时，计算量巨大，限制了在工程中的应用；

(2)基于自由空间格林函数的边界积分法[9]，该方法先求物体边界处的散射声压，然后统一计算远场的噪声，计算量小而适用于工程应用，但当前的研究一般适用于不可压流动诱发的气动噪声预测；

(3)寻找满足边界散射条件的精确格林函数(Exact/Tailored Green’s Function)法[10-12]，该方法不需要对流动进行高精度计算，简明且易于实现数值计算。规则空间的精确格林函数可以通过理论解析法[11，13]求得，而任意空间的精确格林函数则需要采用数值方法进行求解。Takaishi等[12]基于边界元思想提出了一种求解精确格林函数的边界积分方法，适用于任意物体边界形状，但Takaishi等只计算了精确格林函数的1阶偏导数，且只将该方法应用于涡声理论。本文在文献[13]中对Takaishi等的数值计算方法做了进一步研究，求解了精确格林函数的2阶偏导数，并通过理论解进行了验证。

在文献[13]的基础上，采用数值方法得到的精确格林函数来求解FW-H方程，提出了非紧致边界气动噪声数值预测方法。随后，数值预测了圆柱绕流诱发的气动噪声，对提出的方法进行了验证。

1 数值计算方法

1.1 非紧致边界气动噪声数值计算方法

实质上，求解声学波动方程的格林函数满足Helmholtz方程，而频域的自由空间格林函数是相应的Helmholtz方程的基本解。精确格林函数的表达式可以应用格林定理演绎相应的Helmholtz方程得到。在频域下，方程(1)的解也可以写成通用表达式

其中i为虚数单位，ω为圆频率表示频域下的值，C(f)是立体角函数，是频域的精确格林函数，满足

对于非紧致边界，G(x，y，ω)在边界上满足

式中n是边界S的单位法向量，正方向由边界指向声源区域。将式(6)代入式(4)，即可得到非紧致边界气动噪声数值预测方法的基本方程

假设点z为边界S上的观察点，对于精确格林函数和自由空间格林函数有

式 (9)乘 以G(x，z，ω)减 去 式 (8)乘 以G0(z，y，ω)，再将等式两边同时积分得

根据格林定理，式(10)等号左端可化为

将式(11)代入式(10)，并利用δ函数的性质，可得精确格林函数G(x,y,ω)的最终表达式为

式中

式(12)即是精确格林函数的表达式，等式右边第一项是自由入射项，意义为当边界不存在时，观察点x处的格林函数值；等式右边第二项是散射项，意义为当边界存在时，观察点x处由边界的散射产生的格林函数值。

当物体边界为静止刚性边界时，将式(12)代入式(7)得非紧致边界气动噪声数值预测方法的基本方程的最终表达式为

式(14)等式右边第一项为四极子源噪声，为体积声源直接辐射到观察点的噪声，称为直接入射噪声，等式右边第二项是由于边界存在产生的噪声，这部分噪声包含有边界与流动相互作用产生的噪声和声波传播到边界上后散射的噪声，称为边界辐射噪声。

1.2 精确格林函数数值计算方法

采用式(14)计算气动噪声关键的一步是求出式(12)所示的精确格林函数的解。将源点y无限趋近于边界上的点z¯，对式(12)有

式中Φ(z¯)是位于边界上的点z¯处的内固体角，S′是边界S除去积分奇点后的部分。将边界S离散为M个单元，式(15)离散为线性方程组如下

式中E是单位对角矩阵，H的子项满足下式

其中m=1,2,3,…M，n=1,2,3,…M。

通过数值方法求解线性方程组(16)，将结果代入式(12)的等号右端第二项，即可实现对精确格林函数的数值计算。

式(14)用到了精确格林函数的2阶偏导数，下面说明其求解方式。根据公式(12)，可得精确格林函数2阶偏导数的表达式为

对于自由空间格林函数G0(z,y,ω)，其偏导数存在理论解析解，式(18)右端自由空间格林函数各偏导数项均可由理论解析解求出，同时边界S上各点的G(x,z,ω)值也已经由式(15)求出，将这些已知结果代入式(18)的右端，即可求得精确格林函数的二阶偏导数值。

2 圆柱非紧致边界气动噪声计算算例

2.1 流动计算结果及分析

选择流动计算参数如下：圆柱直径D=0.1 m，雷诺数ReD=90 000，空气来流马赫数Ma=0.2。流场区域采用分块结构化网格进行划分，圆柱壁面第一层网格的厚度为2×10-4D，保证近壁面网格的无量纲距离y+在1的数量级，尾涡区域及圆柱附近区域做了加密处理。流动采用二维非定常不可压缩流体的雷诺时均Naiver-Stokes方程（RANS）求解，湍流模型采用k-ωSST，时间步长取Δt=5×10-5，保证一个脱落涡周期内有不少于100个时间步。

图1为圆柱壁面升力系数Cl和阻力系数Cd随时间的变化历程。可以看出，二维圆柱绕流的Cl和Cd呈典型的周期性波动。图2为对Cl和Cd进行的功率谱分析，Cl的基频为f0=157.5 Hz，相应的无量纲量斯特劳哈数St=0.231，Cl明显的谐波频率为mf0(m=3，5，7…)，Cd的基频为2f0，Cd明显的谐波频率为nf0(n=2，4，6……)。这一结果与Takaishi等[12]和Reinaldo等[14]所得结果一致。

图3为得到的圆柱壁面时均压力系数分布与Reinaldo等的计算结果及Cantwell实验结果[15]的对比。可以看出，本文结果与Reinaldo等的结果基本一致。但与实验值相比，二者在转捩点附近的压力系数都偏小，这是由于采用二维计算模型，忽略了圆柱绕流的三维展向效应所致。本文的重点在于提出和验证气动噪声计算方法，这里不再对流动做进一步的计算和验证。

图1 升力系数和阻力系数随时间变化情况

图2 升力系数和阻力系数的功率谱密度

图3 圆柱壁面平均压力系数分布

2.2 气动噪声计算结果及分析

选取26个周期的流场数据，总采样步数为3 303步，根据采样定理，声场计算频率范围选取为50 Hz～2 000 Hz，频率分辨率为6.1 Hz。选取的声压指向性观察点在以圆柱圆心为原点、半径为128D的圆上，声场分布云图观察点选取在以圆柱圆心为原点、半径为20D的圆外。

圆柱空间存在精确格林函数的理论解析解[13]，可以用来验证精确格林函数的数值计算方法。图4是采用理论解析方法和本文提出的数值方法计算精确格林函数后代入式(14)得到的声压指向性图的对比，可以看出，两种计算方法的声压一致，说明精确格林函数的数值计算方法是可靠的。

图4 数值结果和理论解析结果的声压指向性图

图5是采用本文数值计算的精确格林函数求解FW-H方程得到的声压和采用自由空间格林函数求解FW-H方程得到的声压的指向性图对比。低频情况下，当频率为f=f0和f=2f0时(对应的声波波长分别为21.6D和10.8D)，两种方法得到的结果基本一致。随着频率增高，当频率为f=3f0和f=4f0时(对应的声波波长λ分别为7.2D和5.4D)，两种方法得到的结果出现了较为明显的差异。在低频情况下，声波波长远大于圆柱边界特征尺寸，圆柱边界近似可看作满足声学紧致条件，噪声预测方法得到的结果与采用自由空间格林函数求解FW-H方程得到的结果基本一致。随着频率的增高，声波波长和圆柱边界特征尺寸越来越接近，圆柱边界已是声学非紧致边界，边界对声波传播的影响将变得明显且不能忽略，本文噪声预测方法得到的结果与采用自由空间格林函数求解FW-H方程得到的结果出现了明显的差异，表明噪声预测方法能够适用于当边界为声学非紧致边界时的气动噪声预测，考虑了边界对声波的散射效应。

图5 不同计算方法得到的声压指向性图

图6和图7分别给出了噪声预测方法得到的声压分布云图和各成分指向性图。可以看出，圆柱绕流诱发的气动噪声总声压呈明显的偶极子特征，总声压和边界辐射声压相当，直接入射声压对总声压的贡献小，说明低频时的噪声主要由边界辐射噪声组成，这与Takaishi等[12]的结果一致。

图6 总噪声声压分布云图

通过以上分析可知，边界存在时，非定常绕流诱发的气动噪声的边界辐射噪声从本质上看有两部分：一是通常认为的等效偶极子声源产生的噪声，这是边界与流动相互作用产生的；二是边界对声波传播的散射效应，这部分噪声来源于声波传播与边界的相互作用。从本质上看，偶极子源噪声和边界对声波的散射都是边界对声波传播的散射效应的外在表现，所不同的是，偶极子源噪声是流场计算所能够捕捉到的声学脉动量，噪声预测方法和采用自由空间格林函数求解FW-H方程都能预测，而非紧致边界对声波传播的散射是噪声预测方法中波动方程所能捕捉到的声学脉动量，本文噪声预测方法的意义就在于能准确预测这部分噪声，而采用自由空间格林函数求解FW-H方程则不能捕捉到这部分声学脉动量。

图7 噪声各成分指向性的对比图

3 结语

通过数值方法寻找满足非紧致边界的精确格林函数，提出了一种基于FW-H方程的非紧致边界气动噪声数值预测方法，适用于几何外形复杂的边界与非定常流动相互作用诱发的气动噪声的数值预测。通过对圆柱绕流诱发的气动噪声的数值预测，对该方法进行了验证，结果表明：

(1)通过采用圆柱边界的精确格林函数的解析解来求解提出的气动噪声数值预测基本方程，与本文数值计算方法的结果进行对比，验证了精确格林函数数值计算方法的正确性；

(2)与采用自由空间格林函数求解FW-H方程的传统方法相比，提出的气动噪声数值预测方法能够考虑非紧致边界对声波的散射效应，能更为准确的预测流动诱发的气动噪声。

[1]Lighthill M J.On sound generated aerodynamically[M].I. General Theory.Proc.Roy.London Soc.1952,211A, 1107:564-587.

[2]Ffowcs Williams J E,Hawkings D L.Sound generation by turbulence and surface in arbritary motion[J].Proceedings of the Royal Society of London,Series A,Mathematics and Physical Sciences,1969,264:321-342.

[3]Farassat F,Jay H C.Towards an airframe noise prediction methodology:Survey of current approaches[R].AIAA Paper,2006-0210.

[4]Curle N.The influence of solid boundaries upon aerodynamic sound[J].Proceedings of the Royal Society of London,Series A,Mathematics and Physical Sciences, 1955,231,1187:505-514.

[5]Goldstein M E.Aeroacoustics[M].New York:McGraw-Hill International Book Co.,1976.

[6]Hu F Q,Hussaini Y H,Manthey J L.Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics[J].Journal of Computational Physics,1996, 124:177-191.

[7]Ffowcs Williams J E,Hall L H.Aerodynamic sound generation by turbulent flow in the vicinity of a scattering half plane[J].Journal of Fluid Mechanics.1970,40(4):657-670.

[8]Ali I,Escobar M,Kaltenbacher M,et al.Time domain computation of flow induced sound[J].Computer&Fluids,2008,37:349-359.

[9]Mao Y,Qi D,Gu Y.Prediction of airfoil noise induced by low Mach number flow[J].Proc.ImechE,Part G,Journal ofAerospace Engineering,2012,226(5):561-573.

[10]Hu F Q,Guo Y P,and Jones A D.On the computation and application of exact greens function in acoustic analogy[C].In 11th CEAS/AIAA Aeroacoustics Conference. Monterey,California,AIAApaper 2005-2986,2005.

[11]Golerfelt X,Perot F,Bailly C,et al.Flow-induced cylinder noise formulated as diffraction problem for low Mach numbers[J].Journal of Sound and Vibration, 2005,287:129-151.

[12]Takaishi T,Miyazawa M,Kato C.Computational method of evaluating noncompact sound based on vortex sound theory[J].The Journal of Acoustical Society of America,2007,121:1353-1361.

[13]宋亚辉，刘秋洪，蔡晋生.气动噪声预测中精确格林函数的数值计算方法[J].噪声与振动控制，2013，33(2)：5-10.

[14]Reinaldo M O,Julio R M,Fabio S.Two and three-dimensional simulation of sound generated by flow around a circular cylinder[J].AIAAPaper,2009,3270.

[15]Cantwell B,Coles D.An experimental study of entrainment and transport in the turbulent near wake of circular cylinder[J].Fluid Mech.,1983,136:321-374.

Computational Method ofAerodynamic Noise Induced by Airflow over Non-compact Boundary

SONG Ya-hui1,QIN Hao1,LIU Qiu-hong2
(1.Chinese Flight Test Institute,Xi’an 710089,China; 2.School ofAeronautics,Northwestern Polytechnical University,Xi’an 710072,China)

A computational method was proposed to evaluate the aerodynamic noise generated by unsteady flow over non-compact boundary.Firstly,the boundary element method was applied to compute the exact Green’s Function that satisfies the specific boundary conditions.Then,the numerical results of the exact Green’s Function were used to solve the FWH equation to predict the far field noise.This numerical method can predict the radiation noise from the aerodynamic sources as well as the scattering noise from the non-compact surfaces.Especially,this approach is suitable for the boundaries with large sizes and complex geometric shapes.The present method was applied to predict aerodynamic noise generated by turbulent flows at Reynolds number of 90 000 and Mach’s number of 0.2 over a circular cylinder with non-compact boundary conditions.Results were compared with those of FW-H equations based on the free-space Green’s function to verify the effectiveness of this method.

acoustics;aerodynamic noise;FW-H equation;non-compact boundary;numerical method

O42；O355

：A

10.3969/j.issn.1006-1335.2015.02.011

1006-1355(2015)01-0052-06

2014－09－12

国家自然科学基金(11002116)；西北工业大学基础研究基金(GCKY1006)

宋亚辉（1985－），男，河南省鹿邑县人，硕士，工程师，主要研究方向：飞机结构动力学。E-mail:songyahuilym@163.com

刘秋洪，男，副教授，主要研究方向：计算流体力学和气动声学计算。E-mail:qhliu@nwpu.edu.cn