气动噪声预测中非紧致格林函数的数值计算方法

宋亚辉，刘秋洪，蔡晋生

（西北工业大学 翼型叶栅空气动力学国防科技重点实验室，西安 710072）

对于非紧致结构(即边界特征尺寸大于声波波长或与之相当)，非定常流动激发的气动噪声在向空间辐射的同时，会在结构壁面上发生反射和散射，导致气动噪声的实际声场与自由空间下的辐射声场存在明显不同[1,2]，因而需要考虑非紧致边界对气动噪声的影响。当直接采用声模拟理论对非紧致结构气动噪声进行计算时，就需要寻求满足边界散射条件的非紧致格林函数。

非紧致格林函数的求解有理论解析法和数值计算法两种。对于形状较规则的结构，非紧致格林函数可通过理论解析法求得，如Gloerfelt[3]通过理论推导得到了圆柱空间非紧致格林函数的解析表达式。遗憾的是，Gloerfelt推导的非紧致格林函数的二阶偏导数存在错误，文献[4]虽然指出了这一错误的存在，但也没有给出正确的结果。对于外形复杂的结构，难以得到非紧致格林函数的解析表达式，需要采用数值方法求解。Hu和Jones等[5―7]发展了一种基于基函数的数值计算非紧致格林函数的方法，但是该方法计算过程繁琐，不利于程序开发。Takaishi等[8]基于边界元思想提出了一种求解非紧致格林函数的边界积分方法，该方法简明且易于实现数值计算，但文献[8]只计算了非紧致格林函数的1阶偏导数，当应用FW-H方程[9]计算噪声时，还需要进一步计算非紧致格林函数的2阶偏导数。

本文借鉴文献[8]的思想，探讨了非紧致格林函数的数值计算方法，并针对圆柱这一简单外形，采用数值方法计算了非紧致格林函数及其1阶、2阶偏导数，并与理论解析结果进行了对比，检验了数值方法的可靠性。

1 非紧致格林函数

为了简化问题，本文考虑高雷诺数、低马赫数、边界静止时频域的FW-H方程[3]

式中x和y分别表示观察点和源点位置，Ω是声源区域，S是边界，p表示单位面积固体表面作用在流体上的作用力，ρa表示声学密度脉动，ω表示圆频率，Tij为Lighthill应力张量，G(x,y,ω)是频域的格林函数，满足

对边界为声学硬边界的非紧致结构，格林函数G(x,y,ω)在边界上满足

式中n是边界S的单位法向量，正方向由边界指向声源区域。故式(1)化为

式(4)即是采用非紧致格林函数所得到的考虑边界散射的FW-H方程。

对于边界S上的点z，根据格林定理，有

式中，G0(z,y,ω)是自由空间格林函数，满足

将式(2)、(3)、(6)代入式(5)，并利用δ函数的性质，可得G(x,y,ω)的最终表达式为

2 计算方法及结果

为简化计算，选择二维圆柱算例，采用数值计算方法和理论解析方法计算非紧致格林函数及其偏导数，检验计算方法的可靠性。

2.1 非紧致格林函数的数值计算方法

为了便于区分符号，用GN(x,y,ω)表示式(7)中的非紧致格林函数。将流体中的源点y趋近于边界上的点zp，对式(7)有

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

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

通过数值方法求解线性方程组(11)，将结果代入式(7)、(8)和(9)等的等号右端第二项，即可实现对非紧致格林函数及其偏导数的数值计算。

2.2 非紧致格林函数的理论解析方法

二维圆柱，非紧致格林函数的解析表达式易于求得，可用于对上述数值计算方法的验证。文献[3]和文献[4]给出了圆柱空间非紧致格林函数及其1阶、2阶偏导数的解析表达式，但2阶偏导数均存在错误。本文略去推导过程，直接给出非紧致格林函数2阶偏导数的正确结果。

图1所示坐标系，当rx≥ry时，非紧致格林函数2阶偏导数的正确解析表达式为

图1 坐标系及各变量的意义Fig.1 Notations and coordinate systems

2.3 计算结果

选择圆柱直径a=0.05m，点源位于(2a,0°)，取波数k=40(这时声波的波长λ=0.157m≈3a，声波波长与边界特征尺寸相当，满足非紧致声学边界条件)。取观察点位于直径120a的圆上和位于距圆柱圆心10a至200a的圆环区域，采用数值方法求解式(7)至(9)得到数值结果(Numerical result)，并与解析结果(Analytic result)进行对比。

图2是观察点位于直径为120a圆上时非紧致格林函数的数值结果和解析结果的幅值指向性曲线对比，图中角度坐标的取值原则和图1所示相同。从图中可看出数值结果与解析结果吻合良好，由于非紧致边界的存在，导致格林函数及其偏导数的空间分布出现了“瓣状”。

图2 非紧致格林函数的数值结果和解析结果指向性曲线Fig.2 Directivity plots of numerical result and analytic result of non-compact Green’s function

图3是观察点位于距圆柱圆心10a至200a圆环区域时非紧致格林函数二阶偏导数的数值结果和解析结果的实部、虚部云图对比。从图中可以看出，数值结果和解析结果基本一致，由于非紧致边界的存在，使得格林函数的值在观察区域随着距离的增加在某些方向上出现一定规律的波动。

图4是观察点位于直径为120a圆上、不同波数时，非紧致格林函数二阶偏导数的自由入射项(Incident value)、散射项(Scattered value)及其两项之和(Total value)的对比。从图中可以看出，非紧致边界的存在使格林函数的值在不同方向上出现了波动，且随着波数的增大，非紧致边界的散射的影响越来越明显。

图3 非紧致格林函数的数值结果和解析结果的云图Fig.3 Contour plots of numerical result and analytic result of non-compact Green’s function

3 结语

基于边界元思想，对气动噪声预测中非紧致格林函数的数值计算方法进行了初步研究。该方法适用于任意外形结构，易于程序开发，可有效节约计算资源。

为验证该数值计算方法的可靠性，应用该方法计算了二维圆柱的非紧致格林函数及其1阶、2阶偏导数，数值计算结果与理论解析结果一致，表明该数值计算方法是正确的。另一方面，此数值计算方法将非紧致格林函数分为散射项和入射项两部分，可定量分析壁面散射对声场的影响。算例计算结果表明，随着波数的增加，非紧致声学边界的散射作用越来越明显。

笔者今后将采用本文非紧致格林函数的数值计算方法结合式(4)计算非紧致结构的气动噪声。

图4 不同波数时 |∂2GN(x,y,ω)∂y1∂y2 |的自由入射项、散射项及其两项之和的对比Fig.4 Comparison of the incident value,scattered value and total value of |∂2GN(x,y,ω)∂y1∂y2|

[1]Powell A.Aerodynamic noise and the plane boundary[J].The Journal of the Acoustical Society of America,1960,32(8):982-990.

[2]Ffowcs 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.

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

[4]毛义军.低马赫数流动中声比拟理论的拓展及其应用研究[D].西安：西安交通大学能源与动力学院，2010，15-17.

[5]Hu F Q,Guo Y P,Jones A D.On the computation and application of exact Green’s function in acoustic analogy[C].11th AIAA/CEASAeroacousticsConference,Monterey,California,May 2005,on CD-ROM.

[6]Jones A D,Hu F Q.An investigation of spectral collocation boundary element method for the computation of exact Green’s functions in acoustic analogy[C].12 th AIAA/CEASAeroacousticsConference,Cambridge,Massachusetts,May 2006,on CD-ROM.

[7]Jones A D,Hu F Q.A three-dimensional time-domain boundary element method for the computation of exact Green’s functions in acoustic analogy[C].13 th AIAA/CEASAeroacousticsConference,Cambridge,Massachusetts,May 2007,on CD-ROM.

[8]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.

[9]Ffowcs Williams,J.E.and Hawkings,D.L.Sound generation by turbulence and surface in arbritary motion[C].Proceedings of the royal society of London,series A,mathematics and physical sciences,1969,264:321-342.