倪大明，张文平，明平剑

(哈尔滨工程大学动力与能源工程学院，黑龙江哈尔滨150001)

计算气动声学已经成为预测声产生和传播的普遍工具［1］，其主要分为两大类:1)是基于 Lighthill［2］声比拟分析，结合声源(单极子，双极子，四极子)求解线性波动方程;2)是求解可压缩的欧拉或N-S方程［3］.前者主要是在声源和远场建立关系，在其线性模型中反射、衍射以及一些非线性因素不被考虑，其主要用于预测声场远场，而后者考虑了非线性因素，但是由于计算资源的限制，模拟声传播很远的距离计算量是很大的，因此，其主要被用于预测声场近场.

对于高速流，Tam［4］提出了一种可压缩欧拉和可压缩N-S耦合的方法预测近场的声产生和传播.对于低速流，Hardin［5］在交错网格的基础上提出了二步法可预测声产生和传播，SHEN等在其基础上在极坐标系下改进了二步法，并将其推广到层流［6］和湍流［7］.对于均匀流场中的声传播问题，张荣欣等［8］采用切比雪夫谱元方法求解亚音速均匀流场中的辐射声场，X.Nogueira等［9］在基于非结构化网格下提出了一种高精度有限体积法可以求解亚音速和超音速均匀流场中的辐射声场.

基于非结构化网格，采用交错网格计算比较复杂，因此，本文在直角坐标系下基于非结构化同位网格应用有限体积法结合SIMPLE算法模拟低速流中的声传播问题，在界面声流速计算采用 Rhie-Chow［10］插值，在边界上采用DRP(保持色散关系)吸收边界，最后进行了算例的模拟并和声学理论分析解进行了对比分析.

1 控制方程

可压N-S方程组由连续方程、动量方程、绝热等熵的状态方程组成:

式中:ρ、u、p、c、S 分别为控制体中心的密度、速度矢量、压力、声传播的速度和动量方程源项，ρ=ρ(x，y，z，t)，u=u(x，y，z，t)，p=p(x，y，z，t)，c=c(x，y，z，t).

对于低速流，不考虑声反向散射对流场解的影响，将可压N-S方程分解为流体动力项和扰动声项，而流体动力项按照不可压流场求解，用二步法将以上基本变量进行分解，则可压变量:

式中:U是不可压速度矢量;P是不可压压力;ρ0是不可压密度;u'是声扰动速度矢量;p'是声扰动压力;ρ'是声扰动密度.

将式(4)代入方程(1)～(3)，并结合不可压N-S方程，同时令矢量f=ρu'+ρ'U，则声扰动方程的连续方程、动量方程和等熵状态方程分别为

为了消除边界反射对内部声场的影响，采用DRP 无反射边界［11］:

式中:un为可压缩法向速度，r为边界面中心到声源距离，ui'为声扰动速度分量，n对于一、二、三维问题分别取 0、1/2、1.

2 数值离散

不可压流场和声扰动场动量方程离散在时间项采用二阶前差格式，对流项采用二阶迎风格式，扩散项和压力项采用中心差分.在求解界面不可压流场流量和声扰动场流量采用Rhie-Chow插值.不可压流场求解应用明平剑［12］GTEA程序.

2.1 动量预测

根据动量方程(6)离散得到矢量f的离散方程:

式中:AP为计算单元中心矢量f的系数，Ai为计算单元邻近单元中心矢量f的系数为计算单元中心矢量f的预测值为计算单元邻近单元中心矢量f的预测值，V为计算单元体积.连续方程(5)离散得到:

式中:Σ f_f为单元界面上矢量f流量之和.

求解动量方程(9)，得到矢量f*，然后计算单元表面上的矢量f*流量:

2.2 压力修正

引入矢量f修正量f'和声压修正量pp'，则f＊＊=f*+f'，p'＊＊=p'*+pp'，同时根据动量方程(9)，有:

式(13)与(9)相减，得到矢量f修正量f'更新公式:

同时结合连续方程(5)和绝热状态方程(7)，可得

根据式(14)、(15)，可得声压力修正方程:

以上给出的SIMPLE算法以矢量形式给出，与具体坐标系及网格形式无关，适用于任意混合网格.

2.3 动量修正

声压修正后，根据式(14)更新矢量f值，其他变量按下式更新:

式中:u'i为声扰动速度分量，fi为矢量变量分量，Ui为不可压速度分量，Pref为环境参考压力，γ为比热，tref为环境参考温度，上标n+1为当前时刻计算值，n为上个时刻值，n-1为上2个时刻值.值得注意的是:若流动介质可以作为完全理想气体，则声传播速度按照式(19)更新;若流动介质为水，则声传播速度按照介质声传播速度经验公式(20)更新.

3 求解步骤

图1 系统工作流程Fig.1 Workflow of system

控制方程的求解包括不可压流场和声扰动场计算，如图1所示，其求解步骤为:

1)根据上一个迭代层得到的U*、P*，以及上一个时间层的U0、P0，求解不可压流场动量方程，得到新的U＊＊;

2)求解不可压压力修正方程，得到修正值P';

3)计算出新的P、U;

4)根据求出的P、U和上一个迭代层得到的矢量f*、声扰动压力p'*，以及上一个时间层的f0、p0'，求解声扰动动量方程，得到新的f＊＊;

5)求解声扰动压力修正方程，得到修正值p″;

6)计算出新的 p'、f;

7)根据式(17)～(20)，更新 ρ'、ui'、c;

8)判断迭代是否收敛，如果收敛，则进入下一个时间层，否则返回步骤1)继续求解;

9)判断时间是否已经达到设定值，如果达到则计算结束，反之，将目前值存入 U0、P0、p0'、f0、ρ0'、ui0'、c0，返回步骤1)，进入下一个时间层进行迭代.

4 均匀流场算例验证

为了方便研究，本文采用具有解析解的单极子点声源，其表达式如下:

式中:un'为声源法向振动速度;A为声源振幅;ω为角频率.取 A=0.006，ω =68π.

4.1 流动对平面波传播的影响

计算区域为Ω={(x，y):-25 m≤x≤25 m，-0.2 m≤y≤0.2 m}，在人工边界 x=25 m 上面采用DRP吸收边界条件，边界x=-25 m给定点声源，将计算区域划分为400个非结构化网格，时间步长为 Δt=0.000 1 s，计算2 000 步.t=0.2 s时，沿 x轴的声压型线在Ma=0和0.2的对比见图2.从图2可以看出，Ma=0时，声源波长为10 m;Ma=0.2时，声源波长为12 m，其模拟结果符合声传播的多普勒效应，与理论分析吻合.同时可以看出，边界处理较好，反射微小.

图2 t=0.2 s，Ma=0，0.2 时沿 x 轴的声压型线对比Fig.2 Comparison of sound pressure distribution of numerical results on x axis，t=0.2 s，Ma=0，0.2

4.2 流动对柱面波传播的影响

计算区域为Ω={(x，y)∶-50 m≤x≤50 m，-50 m≤y≤50 m}，在人工边界∂Ω上面采用DRP吸收边界条件，点声源位置坐标(0，0)，将计算区域划分为8 861个非结构化网格.

取流场马赫数 Ma=0，时间步长为 Δt=0.000 1 s，计算2 000步.可得到静止流场中由单极子声源产生的辐射声场云图，见图3.沿x轴的声压型线与理论解析解的对比见图4，可以看出，数值解和解析解吻合得良好.

图3 Ma=0时的声场云图Fig.3 Distributions of sound pressure at different time，Ma=0

图4 Ma=0时沿x轴的声压型线与理论解对比Fig.4 Sound pressure distributions of numerical results versus theoretics on x axis at different time，Ma=0

取流场马赫数 ，其余参数不变，可得Ma=0.2流场中单极子声源产生的声场云图见图5，从图6可以看出数值解与解析解吻合得很好，包括边界条件的处理.

图5 Ma=0.2时的声场云图Fig.5 Distributions of sound pressure at different time，Ma=0.2

图6 Ma=0.2时沿x轴的声压型线与理论解对比Fig.6 Sound pressure distribution of numerical results versus theoretics on x axis at different time，Ma=0.2

5 圆柱绕流流噪声模拟

计算区域为 Ω={(x，y)∶ -50 m≤x≤50 m，-50 m≤y≤50 m}，在人工边界∂Ω上面采用DRP吸收边界条件，将计算区域划分为25 421个非结构化网格.圆柱直径为D=1 m，流场马赫数Ma=0.2，雷诺数Re=200.时间步长为t=0.002 5 s，总计算时间为 t=3.75 s，当不可压流场计算趋于周期变化时，即当t=2.5 s时，开始计算声场，监测点设于坐标为(0 m，3 m)处.

若考虑圆柱表面有微小振动，为了方便研究，本文采用具有解析解的单极子声源，其表达式见式(21).为了将流噪声与该声源区分开，取A=0.06，ω =68π.

当t=2.75 s时，声场计算趋于稳定.取t=3 s时，半径为r=3 m圆周上的声压值，绘制成声压指向性图如图7.从图中可以看出，马赫数为0.2时，圆柱绕流流噪声的声压级处于81～123 dB之间，且沿着x轴对称分布，最大声压级处于9°和-9°位置处.由于给定声源的压力辐值5 Pa小于流噪声的压力辐值38 Pa，因此给定声源与流噪声的干涉结果对圆柱上游影响较小，主要影响区域为圆柱上游区，其使流噪声的声压级稍微有所增大.

图7 圆柱表面是否微振指向性对比(r=3 m，t=3 s)Fig.7 Directivity pattern of circular cylinder noise radiation，considering circular cylinder surface slight shock or not，r=3 m，t=0.3 s

图8 圆柱表面是否微振监测点声压随时间变化对比Fig.8 Sound pressure distribution of moniter point vs time，considering circular cylinder surface slight shock or not

图9 圆柱表面是否微振监测点声压随频率变化对比Fig.9 Sound pressure distribution of moniter point vs frequency，considering circular cylinder surface slight shock or not

观察监测点声压随时间变化如图8，从图中可以看出，单计算圆柱绕流流噪声时，声压随时间变化分布很规则，最大声压38 Pa，最小声压-26.5 Pa.通过FFT变换后，可得声压随频率变化如图9，可以看出，圆柱绕流流噪声主要集中在频率为14 Hz处，可计算出斯特劳哈尔数:

其与声学手册中风吹声斯特劳哈尔数一般为0.2相对应.

若模拟流噪声与给定声源相干涉问题，从图8中可以看出，监测点的声压随时间变化分布与单计算流噪声的声压分布有明显不同.通过FFT变换后，可得声压随频率变化如图9，可以看出，噪声主要集中在频率为14 Hz和34 Hz处，14 Hz对应圆柱绕流流噪声的主要集中频率，而34 Hz对应给定声源的频率，其计算结果与理论分析吻合，进一步验证了该方法的正确性和可靠性.

6 结束语

采用基于可压N-S方程的分解法对低速流中的声传播问题进行数值模拟，并采用基于非结构化网格的有限体积法和SIMPLE算法求解不可压流场方程的和声扰动方程，在计算界面流速和声扰动速度采用Rhie-Chow插值，在边界上采用DRP吸收边界条件消除边界反射对声场的影响.对低速流中平面波、柱面波、圆柱绕流流噪声以及流噪声与已知声波干涉传播问题进行了数值模拟，与理论分析解对比表明，这种数值模拟方法预测的声传播是稳定可靠的.由于非结构化网格的应用，该数值模拟方法可便捷求解结构复杂的低速流中的声传播问题.

［1］SHEN W Z，MICHELSEN J A，SORENSEN J N.A collocated grid finite volume method for aeroacoustic computations of low-speed flows［J］.Journal of Computational Physics，2004，196:348-366.

［2］LIGHTHILL M J.On sound generated aerodynamically.I:General theory［C］//Proceedings of the Royal Society A.London，1952.

［3］PHILIP J M，LYLE N L，ASHOK B，et al.A parallel threedimensional computational aeroacoustics method using nonlinear disturbance equations［J］.Journal of Computational Physics，1997，133:56-74.

［4］SHEN H，TAM C W.Numerical simulation of the generation of axisymmetric mode jet screech tones［J］.AIAA Journal，1998，36(10):1801-1828.

［5］HARDIN J C，POPE D S.An acoustic/viscous splitting technique for computational aeroacoustics［J］.Theoret Comput Fluid Dynamics，1994，6:323-340.

［6］SHEN W Z，SORENSEN J N.Aeroacoustic modeling of low-speed flows［J］.Theoret Comput.Fluid Dynamics，1999，13:271-289.

［7］SHEN W Z，SORENSEN J N.Aeroacoustic modeling of turbulent airfoil flows［J］.J AIAA，2001，39(6):1057-1064.

［8］张荣欣，秦国良.用切比雪夫谱元方法求解均匀流场中的声传播问题［J］.西安交通大学学报，2009，43(7):120-124.

ZHANG Rongxin，QIN Guoliang.Chebychev spectral elements method for acoustic propagation problem in a uniform mean flow［J］.Journal of Xi'an Jiaotong University，2009，43(7):120-124.

［9］NOGUEIRA X，COLOMINAS I，FELGUEROSO L C，et al.Resolution of computational aeroacoustics problems on unstructured grids with a higher-order finite volume scheme［J］.Journal of Computational and Applied Mathematics 2010，234(7):2089-2097.

［10］RHIE C M ，CHOW W L.Numerical study of the turbulent flow past an airfoil with trailing edge separation［J］.AIAA Journal，1983，21(11):1525-1532.

［11］TAM C M，WEBB J C.Dispersion-relation-preserving finite difference schemes for computational acoustics［J］.Journal of Computational Physics，1993，107:262-281.

［12］明平剑，张文平，雷国东，等.带壁面射流燃烧室燃烧数值模拟［J］.航空动力学报，2008，23(7):1205-1211.

MING Pingjian，ZHANG Wenping，LEI Guodong，et al.Numerical simulation of combustion in a wall jet combustor［J］.Journal of Aerospace Power，2008，23(7):1205-1211.