1 引言
火箭发动机通过稳定燃烧将自携的推进剂化学能转化为工质内能，高速排出产生反向推力，并产生巨大的射流噪声。大型火箭发动机的燃气射流噪声级别可达170 dB，将对火箭发射结构产生冲击破坏并引发振动破坏，同时引起发射人员听觉功能下降[1]。为应对严峻的噪声威胁，各国对射流噪声进行了长期且广泛的研究。19世纪60年代，Lighthill由流体力学基本方程出发，提出了Lighthill应力张量，建立了声类比理论，标志着气动声学领域的开篇，为气动声学研究提供了理论依据[2]；而后Curler 率先考虑固体表面对流动噪声的影响，成功解释了圆柱旋涡脱落引发的噪声问题，但Curler理论并未涉及运动固体于流体中的边界发声问题[3]；而Ffowcs Williams 和Hawking 应用广义Green 函数方法，提出了FW-H模型，解决了上述Curler理论缺陷，推动力声类比理论的进展及应用[4]。理论建立后，国外许多研究人员进行了应用研究。Andersson[5]首次利用有限体积法对马赫数为0.75的射流外流场及声场进行了LES数值模拟，并被试验结果验证，证明了有限体积法在噪声计算中的可行性。Evren等[6]利用大涡模拟研究了火箭发射过程中，焰流垂直冲击平板的声学特性，并验证了一种空隙隔板的降噪方法。Wang等[7]利用大涡模拟和FW-H方法预测了距机翼的不同距离航空发动机的噪声，与试验结果几乎完全一致。而国内自20世纪70年代也针对火箭发动机的射流噪声问题进行了各种研究。马大猷等[8-9]提出了在声速和亚声速条件下，声功率与气室超压的2.3次方成比例经验公式。张磊等[10]通过数值模拟与试验，研究了喷管尺寸对速燃固体火箭射流噪声的影响。结果表明，喷管尺寸增大声压级增大，出口马赫数相同的喷管声压级比较接近。何泽鹏等[11]对亚音速自由射流噪声进行了模拟研究。结果表明，噪声受射流出口流速、理想膨胀状态、出口马赫数及雷诺数的影响，其主要沿射流轴线和射流剪切边界层2个方向传播。以上研究提出了CFD(computational fluid dynamics)和CAA(computational aero acoustics)方法的理论，并进行了应用验证。
本文针对小型加长喷管模拟发动机，采用三维瞬态CFD/CAA的方法，应用大涡模拟和FW-H表面积分理论，分析了在总温3 430 K，总压1 MPa的工况下，射流噪声场总声压级和频域特性的空间分布。
2 物理模型与计算方法
2.1 物理模型
理想喷管射流场结构如图1所示[12]。火箭发动机产生的高温高压高速燃气由喷管喷出后以湍流状态向下游传播，在接触大气后向四周进行扩散，与周围介质进行产生激烈的热能和动能交换，随着周围介质的不断卷入，形成射流速度核心区，核心区侧面是燃气射流与周围静止大气之间湍流剪切运动形成的含涡结构的射流剪切边界层，剪切流诱导产生气动噪声[10]。
图1 典型的喷管射流场结构示意图
Fig.1 Schematic of the structure of fluidic injection
图2给出了本文研究的圆管加长的收缩-扩张喷管(拉瓦尔喷管)流场噪声计算域，共被划分为3个区域：火箭发动机喷管区域、近场声源区域、远场流动区域。将坐标原点置于燃烧室入口圆心处，设定喷管出口直径为基准长度D，喷管外流域长度和宽度分别取80D和30D，声源面宽度为8D，长度为50D。喷管具体尺寸如下：燃烧室入口半径为23 mm，收缩段长度为44 mm，喉部半径为2.5 mm，扩张段长度为13 mm，喷管出口直径为12 mm，加长喷管长度为150 mm。
图2 物理模型示意图
Fig.2 Geometry model
2.2 计算方法
CFD计算通常采用RANS(reynolds average navier-stokes)和LES湍流模型，但由于噪声模拟需要计算方法能够捕捉到微小的压力脉动，而RANS方法对湍流脉动会进行时间平均处理，直接导致气动噪声计算结果不准确[15-16]。因此，本文选用LES模型对计算域开展三维瞬态模拟，其核心思想是对于相对尺度较大的涡使用Navier-Stokes方程进行求解，而对于相对尺度较小的涡则通过建立亚格子模型(subgrid scale model,SGS)进行求解。SGS选用Smagorinshy涡粘模型。
为捕获声场数据，基于计算流场域中设定的声源面，选用FW-H表面积分法进行声学模拟求解，将流场的压力脉动数据提取传递为远场噪声的声源。FW-H方法采用无限制、无粘波动传播，在波传递到表面的过程中减少了数值上的色散和耗散，其声源面是一个封闭的表面[14]。FW-H方程表达式如下：
其中Ui=[1-(ρ/ρ0)]vi+(ρui/ρ0)；Li=Pijnj+ρui(un-vn)；p′为远场声压；u、ρ、Pij分别为流体的速度、密度、应力张量；v为物面速度；c为声速；Tij为Lighthill应力张量，满足
为Heaviside函数；δ(f)为狄拉克(Dirac Delta)函数。
FW-H方程右边分别代表四级子、偶极子和单极子声辐射。面积分为单极子、偶极子和积分面内的四极子声源相对应。四极子声源是由Heaviside函数决定的，为体声源。当高速燃气喷射到静止空气中时，湍流产生大量的旋涡，方程右边第一项中的应力张量发生变化，会产生四级子辐射源。
2.3 网格划分及初始边界条件
本文对流场计算域采用分块结构化网格划分。为保证计算收敛性以及有效捕捉射流核心区及剪切层压力脉动，对近壁面处和射流核心区进行网格加密处理。壁面网格的厚度为0.01 mm，无量纲数y+≈1。计算域的总网格数量为1 400万，计算域网格如图3所示。
图3 计算区域网格划分和边界设置示意图
Fig.3 The grid division and boundary setting
由于喷管出口做复杂的湍流流动，为了简化求解，进行以下合理的假设：燃气射流不进行化学反应；将燃气射流设定为可压缩的理想气体；周围固壁对声波无反射作用。结合实际工况设置入口边界条件为：燃烧室左部为压力入口Inlet 1，总压p0设置为1 MPa，总温T0设置为3 430 K；左侧柱面入口边界Inlet 2设置为压力入口p1设置为111 325 Pa；设置柱面边界Inlet 3为压力远场边界，p2为101 325 Pa；设置远场边界为压力出口Outlet 2，压力p3设置为101 325 Pa，温度设置为300 K；壁面Wall采用无滑移、绝热壁面条件，壁面函数设置为标准壁面函数；内部声源面Ac-Source设置为Interface。
图4为噪声采集点分布图，图中各点分布在不同的射流轴线夹角θ以及距喷管出口距离r处。
图4 噪声采集点位置分布示意图
Fig.4 Sketch map of noise receivers
3 流场结果及分析
3.1 稳态流场特性
为验证网格加密的有效性和减少数值求解的收敛时耗，在进行瞬态计算前，利用Standard k-ε模型进行了稳态计算，将稳态计算结果作为瞬态流场计算的初始条件。
图5为稳态计算流场结果。在图5(a)中，燃气通过喷管后的温度在2 000 K以下，射流核心区存在交替的低温区，根据流体动力学知识可知该位置处于膨胀波区，核心区的温度梯度较大，而射流外部区域相对较小。图5(b)中，射流在加长喷管内保持1 800 m/s以上的流速，与图1射流结构布局基本一致。速度和温度核心区均以小扩张角的喇叭形区域内分布，未溢出加密的柱形区域，验证了网格加密区域的空间有效性。
图5 稳态计算流场分布云图
Fig.5 Steady jet flow field result
3.2 瞬态流场特性
基于前述稳态结果，应用LES Smagorinshy模型进行了瞬态流场计算，时间步长设置为10-5 s。为弱化稳态结果的影响，设置时间步n=4 000，使得nΔtUj/xf>7.5，Uj为喷管出口速度，xf为整体流场长度[15]。经过6 000核时计算得到如下流场分布结果。
由图6(a)可以看出：燃气射流的静压在流经扩张段后基本稳定在小范围内波动，说明压力的传播容易耗散。而瞬态温度场如图6(b)和速度场如图6(c)存在明显的射流核心区，但瞬态流场并集中在射流中心轴线区域，而在射流核心区的外缘明显可见涡的形成、发展和耗散过程。
4 射流声学结果与分析
脉动压力p*是流场压力的时间导数，也是评判噪声传播的重要依据，根据图7的p*分布情况，可以观察到波向远场传播，而振幅在径向方向衰减。这是球面扩散和数值色散的结果，这表明需要一个声学求解器来求解远场噪声水平。
图6 瞬态流场云图
Fig.6 Transient jet flow field result
图7 脉动压力云图
Fig.7 Fluctuant pressure
4.1 声压级分布分析
声压则定义为p′=pr-pm，式中pr为瞬时压强；pm为平均压强。声压级(sound pressure level,SPL)定义为
SPL=20lg(p′/pref)
式中: p′为当地的压力扰动；pref为参考压强，对于空气，pref=2.0×10-5 Pa。
而总声压级(overall sound pressure level,OASPL)则定义为如下公式。
表1为各采集点的噪声总声压级，由表中数据可知，各接收点的OASPL大为123.6 dB，在θ为30°，r为200 mm处，最小为102.7 dB位于在θ为90°，r为600 mm处，总声压级均在100 dB以上。
由图8总声压级-角度曲线可知，随着射流轴线夹角θ的增大，OASPL总体呈现下降趋势，且折线斜率越来越大，并在50°～60°达到最大，而在靠近90°折线变缓，斜率变小，总声压级下降变缓。在r<400 mm的近喷口区域，噪声在射流轴线方向上的传播占据主导了地位，因此随θ增大，OASPL随之变小，近轴线方向OASPL越大。但在r>400 mm的接收点簇在θ<40°方向OASPL逐渐增大，而后又变小，最终导致在40°方向时OASPL最大。这是因为噪声沿剪切层法向和射流轴线方向的传播耦合作用。综上OASPL在空间上均表现出指向性。
表1 总声压级分布表(dB)
Table 1 The OASPL of each receiver
角度θ/(°)距离r/mm20030040050060030123.6121.1117.5114.4112.440120.4119.6118.1116.1114.150117.3116.0115.4114.5113.460114.4111.8110.7110.2109.870112.6109.3107.4106.3105.480110.4108.1106.2104.5103.390123.3106.4105.1103.6102.7
图8 总声压级-角度曲线
Fig.8 OASPL for different receiver angles
图9为总声压级-距离曲线，由图9可知，随着r的增大各个方向的总声压级曲线都呈下降趋势，越靠近喷管出口总声压级越大，这是因为球面波传播到更大的空间，能量被分散给更多的介质，使得接收点的总声压级减小；而图9中各段折线的斜率逐渐减小，对于越大的距离，对于增加相同的半径，包络面周长的相对变化变小。因此声压级减小变缓。
4.2 频域分布分析
在射流噪声的研究中，除总声压级，噪声的频率特性也是重点分析的内容。在得到较为稳定的瞬态流场解后，使用FW-H声学方程同步对流场数据进行转换，得到喷管的声学特征分布。通过FFT(fast fourier transform)对采集到的数据进行频域特性分析。图10所示为振幅-频率偏移曲线，在等距离变角度采集点的噪声频谱分布呈现出一致性，振幅分布主要在3 000～7 000 Hz范围内，峰值频率分布介于4 000～5 000 Hz之间，而振幅也随角度的增大而减小。
图9 总声压级-距离曲线
Fig.9 OASPL for different receiver distances
图10 振幅-频率偏移曲线
Fig.10 Amplitude-frequency migration diagram
图11为80°方向上距喷管出口不同距离的得到的频谱数据偏移曲线，由图中可以看出频率分布具有宽频性，而在间隔100 mm的接收点的声压频谱分布是：随着频率的增加声压逐渐缓慢增大，而后缓慢减小，最大声压的频率分布在4 500 Hz附近。
图11 声压-频率偏移曲线
Fig.11 SPL-frequency offset diagram
表2给出各个采集点的声压峰值频率。由表2可知：空间各采集点声压峰值频率最小为4 085 Hz和最大为4 875 Hz，极差为790 Hz，峰值频率的平均值为4 401.4 Hz，相对偏差为18.0%。对于整体的求解频率范围相对偏差为7.9%。对于整个宽频域而言，声压峰值频率表现出了明显空间一致性，并没有出现明显的指向性。
表2 各采集点声压峰值频率
Table 2 The peak frequency for receiver SPL
角度θ/(°)距离r/mm200300400500600304 6874 1854 1764 8754 86940 4 6644 1804 1664 0854 180504 1734 1674 1714 6794 166604 2464 2584 2364 2444 175704 6624 4944 8244 8304 516804 2514 4144 3314 5974 497904 1744 3884 4144 6594 416
5 结论
本文以额定工况下的小型模拟火箭发动机为研究对象，采用三维瞬态CFD/CAA的方法，基于LES模型和FW-H表面积分理论，对喷管射流噪声场总声压级和频域特性的空间分布开展了全面分析，得到如下结论：
1) 压力的时间导数云图可以表征噪声的传播方向，但由于球面扩散和数值色散，不能清晰表征远场噪声；
2) OASPL随着距离r的增大而减小，且距离越大减幅越小；
3) OASPL在空间上表现出指向性，在r<400 mm时，OASPL越靠近射流轴线越大，在r>400 mm时，OASPL在40°方向上最大；
4) 噪声频率分布具有宽频性，而声压峰值频率分布具有空间一致性，在各个采集点均在4 000～5 000 Hz。
[1] 陈劲松，何冠杰，吴新跃，等.单喷管火箭自由喷流噪声数值模拟[J].宇航学报，2020,41(03):371-378.
Chen J S，He G J，Wu X Y，et al.Numerical simulation on free jet noise of single-nozzle rocket[J].Journal of Astronautics，2020,41(03):371-378.
[2] Lighthill M J.On sound generated aerodynamically Ⅰ.general theory[J].Proceedings of Royal Society of London Ser A，1951,211(1107):564-587.
[3] Curle，N.The generation of sound by aerodynamic means[J].Aeronautical Journal，1961,65(611):724-728.
[4] Williams J E F，Hawkings D L.Sound generation by turbulence and surfaces in arbitrary motion[J].Philosophical Transactions of the Royal Society of London.Series A，Mathematical and Physical Sciences(1934—1990)，1969,264(1151):539-546.
[5] Andersson N，Eriksson L E，Davidson L.Investigation of an isothermal Mach 0.75 jet and its radiated sound using large-eddy simulation and Kirchhoff surface integration[J].International Journal of Heat & Fluid Flow，2005,26(03):393-410.
[6] Yenigelen E，Morris P J.Numerical investigation of a noise reduction strategy for rocket launch vehicles[C].AIAA Aviation 2020 Forum，2020.
[7] Wang Z N，Proenca A，Lawrence J，et al.Large-eddy-simulation prediction of an installed jet flow and noise with experimental validation[J].AIAA Journal，2020,58(02):1-10.
[8] 马大猷，李沛滋，戴根华，等.小孔喷注噪声和小孔消声器[J].中国科学，1977(05):445-455.
Ma D Y，Li P C，Dai G H，et al.Hole jet noise and hole muffler[J].Science China,1977(05):445-455.
[9] 马大猷，李沛滋，戴根华，等.湍流喷注噪声的压力关系[J].物理学报，1978(02):122-125.
Ma D Y，Li P C，Dai G H，et al.Pressure relations of turbulent injection noise[J].Acta Physica Sinica，1978(02):122-125.
[10] 张磊，阮文俊，王浩等.固体火箭发动机燃气射流流场和声场数值计算[J].固体火箭技术，2015(02):198-202.
Zhang L，Ruan W J，Wang H，et al.Numerical calculation of gas jet flow and acoustics fields for solid rocket motor[J].Journal of Solid Rocket Technology，2015(02):198-202.
[11] 何泽鹏，傅德彬，杜文强.亚音速自由射流噪声的模拟与研究[C]//第十届全国流体力学学术会议论文集.杭州，2018.
He Z P，Fu D B，Du W Q.Simulation and study of subsonic free jet noise[C]//The 10th National Conference on Fluid Mechanics.Hangzhou,2018.
[12] Cuppoletti D R，Malla B，Gutmark E J，et al.The response of supersonic jet noise components to fluidic injection parameters[C].AIAA/CEAS Aeroacoustics Conference，2013.
[13] 余雷，宋文萍，韩忠华，等.基于混合RANS/LES方法与FW-H方程的气动声学计算研究[J].航空学报，2013,34(08):1795-1805.
Yu L，Song W P，Han Z Z，et al.Aeroacoustics noise prediction using hybrid rans/les method and FW-H equation[J].Acta Aeronautica et Astronautica Sinica，2013,34(08):1795-1805.
[14] 李爱琴，李军伟，王宁飞.火箭发动机自由射流噪声模拟过程[J].计算机仿真，2014,31(12):64-68.
Li A Q，Li J W，Wang N F.Rocket engine free jet noise simulation process[J].Computer Simulation，2014,31(12):64-68.
[15] 朱余佳.裂隙性质对含裂隙基坑土坡渗流特性影响数值模拟研究[J].重庆工商大学学报(自然科学版),2021,38(01):105-112.
Zhu Y J.Numerical simulation study on the influence of fracture properties on seepage characteristics of soil slope with fissures[J].Journal of Chongqing Technology and Business University(Natural Science Edition),2021,38(01):105-112.
[16] 张毅，屠珊，洪振瀚，等.喷管入口气温影响下的射流流场及声场特性[J].西安交通大学学报，2020,54(02):175-181.
Zhang Y，Tu S，Hong Z H，et al.Subsonic jet flow field and sound field characteristics influenced by nozzle inlet air temperature[J].Journal of Xi’an Jiaotong University，2020,54(02):175-181.