1 引言
战斗机和轰炸机内埋武器舱可压缩空腔流动噪声可超过160 dB，能够致使机载电子设备失灵、舱内结构疲劳损伤，甚至诱发武器在分离过程中非常规运动，导致投放精度降低[1]。此外，空腔流强噪声辐射还会影响商用航空运输的噪声适航性及军用战斗机、轰炸机飞行隐蔽性[2]。20世纪50年代，Krishnamurty[3]最早开始研究空腔流动现象，并指出空腔流中存在声-流耦合共振的声反馈环机制。当高速气流掠过空腔上方表面后，形成高度有序的剪切层涡列结构，该涡列与空腔后壁撞击产生反射声波扰动空腔前缘初始剪切层，当扰动波与反射压力波频率相匹配时，空腔内出现声场-流场相互耦合的声反馈环流动振荡现象。之后，Rossiter[4]发展了半经验公式对声-流共振现象提供了更完整的描述，并可预测空腔流振荡频率，在空腔流研究中得到广泛应用[5]。
针对超声速飞行条件下战斗机内埋武器分离投放可靠性问题，高保真数值技术在探索超声速空腔流机制方面展现出良好的前景。其中LES(large eddy simulation，大涡模拟)方法可以精确求解空腔内各类涡结构及其产生的噪声特性，具有很强的预测分析空腔流噪声机制的能力。Sinha等[6]和Smith等[7]较早采用三维LES方法分别计算了Mach 2.0和Mach 1.5超声速开式空腔流动噪声。其中，Smith等[7]的LES结果显示压力频谱峰值与实验结果一致。Rizzetta和Visbal[8]采用了基于动态Smagorinsky亚格子模型的LES技术求解了长高比L/D=5，Mach数1.19，Reynolds数 200 000的超声速方腔流动，并应用了高精度紧致格式和大量计算网格，获得的方腔底部多个位置上压力谱与实验结果在幅值和特征频率上趋势一致。但上述研究缺乏细致的流场动力学分析，近期在深化空腔流机理认知以改进噪声控制的需求牵引下，人们再次开始关注展向受限条件下空腔流动基本特性。Brès和Colonius[9]使用DNS(direct numerical simulation，直接数值模拟)方法和整体稳定性理论分析了低Reynolds数1500、亚声速0.325空腔流动展向稳定性，并发现三维空腔模态具有一个与方腔宽度相当的展向波长，其频率比二维Rossiter模态低一个量级。Sun等[10]完成了一个有限宽度超声速开式方腔流动LES计算分析，研究了侧壁效应、自由来流Mach数、Reynolds数等对空腔流动的影响。整体而言，目前对可压缩空腔流动及噪声特性及机理细节认识仍不足。本文采用基于计算气动声学技术的LES方法计算分析三维Mach 1.4和Mach 0.6超/亚声速来流条件下可压缩开式方腔流动。针对空腔流动中剪切层演化、多尺度涡结构、声反馈、可压缩效应等基本动力学过程开展分析，以深化可压缩空腔流动及噪声现象认识和机制理解。
2 数值方法
2.1 控制方程及SGS模型
为实现LES方法，采用无量纲化三维可压缩Favre滤波Navier-Stokes方程
(1)
(2)
(3)
其中
分别为Favre滤波密度、速度、压力、总能、热通量和黏性应力张量。以方腔深度D*，来流密度ρ∞*，来流声速a∞*，来流温度T∞*，来流分子黏性系数μ∞*等作为特征参考量对方程(1)-(3)进行无量纲化处理。对于理想气体，状态方程的无量纲形式可写为
(4)
其中比热比
为Favre滤波温度。方程(1)-(3)中的亚网格尺度(subgrid-scale，SGS)项通过经典Smagorinsky模型建模，SGS黏性应力张量为
(5)
其中
模型系数CS=0.012，CI=0.006 6，湍流Prandtl数也取为常值Prt=0.9。Δ是LES的滤波尺度，与当地网格间距相关。δij为Kronecker符号。SGS热通量建模
(6)
其中cp是等压比热，湍流Prandtl数Prt=0.9。此外，对于SGS黏性扩散项
以及由SGS黏性应力和SGS热通量引起的
和
由于它们在方腔流LES计算中影响远低于SGS热通量和黏性应力张量影响而被忽略。
2.2 计算格式
控制方程采用基于计算气动声学的高精度数值格式进行离散，空间格式使用Tam等[11-12]提出的低色散低耗散误差7点4阶DRP格式，时间推进采用Berland等[13]提出的优化低存储6步4阶Runge-Kutta格式。由于DRP数值格式不能精确离散高波数短波，为避免数值高频波导致Gibbs振荡，引入Tam[11-12]构造的人工选择性阻尼项对非物理数值短波进行滤波。选择性阻尼项、黏性项及SGS项在计算中均在Runge-Kutta最后一个子步引入。对于含强湍流脉动和激波间断的超声速空腔流场，采用Bogey等[14]提出的自适应空间滤波法对守恒变量滤波以捕捉激波并维持计算稳健性。
为了尽量避免边界处出现非物理声波反射，对于亚声速空腔流动，远场辐射边界条件通过线化Euler方程在远场的渐近解构造[12]。对于超声速空腔流动，来流和远场给定超声速自由来流条件，出口边界使用了特征无反射边界条件。空腔固壁采用无滑移
零压力梯度
及绝热壁
边界设置，其中n表示固壁法向。为模拟风洞实验中的侧壁效应，计算模型侧壁远场应用了无穿透滑移壁边界
在超声速和亚声速空腔流计算域的下游出口处均设置海绵吸声区，以避免大尺度涡穿过下游出口边界产生非物理声反射。
2.3 模型设置
计算模型为一个有限宽度的开式方腔，如图1所示，方腔长深比L/D=6，宽深比W/D=2，D、L、W分别为方腔深度、长度和宽度。以方腔深度为特征长度的Reynolds数Re=10 000，来流边界层厚度为δ=0.167。
图1 方腔模型及计算网格示意图
Fig.1 Sketch of cavity model and grid
亚/超声速方腔流动使用相同Cartesian网格，900万个网格点求解空腔流动及噪声。空腔内设置199×103×91个网格点，计算域范围为x1=0～6，x2=-1～0，x3=-1～1。空腔外计算域339×141×183个网格点，计算域范围为x1=-3～13，x2=-1～9，x3=-2.5～2.5。壁面附近对网格进行加密，各方向最小网格尺度为Δx1=0.002，Δx2=0.001，Δx3=0.002。超声速方腔算例进行500 000次迭代，时间步长Δt=0.000 5，流场统计在第100 000步后启动。亚声速方腔算例进行250 000次迭代，时间步长Δt=0.000 5。为加速计算，LES程序采用MPI方法，使用84个进程，每网格点每次迭代约花费0.5 μs CPU计算时间，一次计算花费约480 h。
2.4 Rossiter空腔频率预测模型
Rossiter[4]于1964年提出了预测空腔流噪声特征频率的半经验公式，或称Rossiter模态预测，其形式如下
(7)
式中: fm为第m个模态的频率，u、Ma为自由来流速度和Mach数， γ为比热比，α为相位移。κν为大尺度涡结构对流速度与自由来流速度的比值，是依赖于空腔几何形状和测试条件常数。本文LES计算中在方腔底部(3，-1，0)和中心(3，-0.5，0)设置两个测压点，并将计算压力脉动值与Rossiter预测结果对比，相互验证并确认空腔噪声特征频率。
3 验证与校核
将计算统计得到的空腔流时均流向速度和湍流脉动速度与Dudley和Ukeiley[15]的实验结果和DES数值结果进行对比，以验证LES结果的准确性。图2展示了相同来流条件下本文LES计算和Dudley和Ukeiley的实验测试的方腔 x3=0截面上平均流向速度，两图等值云图范围均在u/U∞=-0.25～1之间，增量均为0.25。图中LES和实验获得的超声速空腔流流向速度沿空间分布完全一致，均展现了明显的剪切流动特征，即方腔上方剪切层沿下游逐渐增厚，受方腔后壁遮挡而浸入腔内的动力学现象，空腔前部及底部平均流向速度相对较低。
图2 时均流向速度场云图
Fig.2 Flow field of time-averaged streamwise velocity
图3展示了方腔x3=0截面上流向脉动速度均方根，LES结果与文献中DES结果[15]采用等值云图范围一致，均为〈u1rms〉/U∞=0～0.25，增量0.05。计算与实验结果均显示了湍流脉动速度沿剪切层向空腔后部逐渐增强，在方腔后部剪切层区域湍流脉动达到最大值的。受方腔后缘诱导，后缘表面也形成一个小范围的强湍流区。值得注意的是，在靠近空腔后底部小的湍流区域，文献DES获得的湍流脉动速度强度高于本文LES计算统计结果。这是由于DES计算的底壁面边界层较薄，易于在腔内剪切层于后壁回流诱导的逆压梯度下，形成较大范围湍流分流区，致使后底壁区呈现较LES结果高的湍流脉动强度。
图4展示了方腔内7个流向站位上LES与实验测量获得的时均流向速度型和流向脉动速度均方根。如图所示，二者很好地吻合，表明LES计算准确地捕捉到了剪切层和方腔流动对流场平均速度及湍流脉动的影响。在空腔下游第六和第七个站位处，数值和实验的时均流向速度和流向脉动速度均方根在靠近空腔底壁面时出现最大约10%的相对误差。对比本文LES及文献DES结果[15]，发现在空腔后底壁面区域，LES计算湍流脉动速度强度略低于实验测试结果，而DES结果则略高于实验结果。可见未来仍需不断深入探索以完善更精确的空腔流动数值分析方法。
图3 时均流向脉动速度均方根云图
Fig.3 Flow field of time-average root mean square of streamwise velocity fluctuation
图4 空腔内不同流向站位时均流向速度和 流向脉动速度均方根示意图
Fig.4 Time-average streamwise velocity and root mean square of streamwise velocity fluctuation in cavity
4 LES结果与分析
4.1 瞬态流场
LES方法可精确计算空腔流动涡、声、激波等全部流场细节，便于对复杂振荡空腔流动的演化特征和机制开展分析探索。如图5，使用Q准则等值面(Q=3)和压力云图展示了Mach 1.4超声速开式空腔瞬态流场。来流进入计算域后，在壁面诱导的气流压缩效应下计算域前缘形成一道激波。气流经过空腔前缘后，与腔内气体混合使得边界层向剪切层转化，且由于流动不稳定性剪切层卷起形成展向涡结构。涡结构在向下游运动过程中沿展向失稳，产生大量Λ涡及其他更小尺度的结构，在高速来流裹挟下继续向下游对流。在空腔后缘，剪切层被撞击产生大量小尺度湍流涡结构，一部分侵入空腔后部与腔内流体急剧混合，另一部分溢出空腔沿壁面继续向下游流出计算域。剪切层撞击空腔后缘还形成向上游传播的声波，声波辐射出空腔后被超声速来流的对流作用下向下游远场传播，在腔内声波穿过亚声速流动辐射至空腔前壁面。
图5 瞬态Q准则等值面(Q=3)及压力云图
Fig.5 Instantaneous Q-criterion isosurface(Q=3) and pressure field
图6展示了空腔中心x1x2截面(x3=0)上4个连续时刻密度梯度幅值云图，由图可见复杂超声速空腔流场动力学演化。首先，图中清晰地展示了剪切层大尺度结构的运动发展过程。在t=139时刻大尺度结构形成并向下游对流，t=139.5时刻该大尺度结构向方腔后缘靠近，同时新的涡结构卷起，t=140时刻原大尺度结构继续逼近方腔后缘且头部结构开始变形，新的涡结构显著增强，t=140.5原大尺度结构与方腔后缘撞击，新大尺度涡结构形成并向下游对流，进而完成了一次大尺度涡结构演化周期。同时，如图6(c)所示，可观察到剪切层附近的亚声速流动区域内形成向上游辐射的声波，其波长与大尺度涡结构相当。这证实了空腔振荡流动中剪切层与空腔后壁面撞击产生向空腔前方传播的声波，并扰动空腔前缘剪切层，形成自持振荡的声反馈环现象。此外，超声速来流在开式空腔诱导下在空腔前缘也形成一道激波，如图6(a)所示，由于与剪切层相互作用，其始终处于振荡状态，并反向作用影响剪切层的发展演化。整体而言，三维超声速开式空腔流动的声反馈环特性较为显著。
图6 密度梯度幅值(纹影图)演化云图
Fig.6 Evolution of density gradient magnitude(schlieren image)
4.2 空腔压力脉动
在空腔底部壁面中心P1点(3，-1，0)和空腔空间中心P2点(3，-0.5，0)两个典型位置设置测压点以观测腔内压力变化。如图7(a)所示，P2点压力脉动幅值略高于P1点，这是由于湍流在空腔空间中的脉动强度强于空腔底壁所致。使用了Fourier变换计算两点上压力频谱，由于压力信号计算及存储时间有限，对原始压力脉动信号做Hamming窗滤波处理以方便识别主频信息。
图7 空腔内P1、P2两测点的压力脉动和频谱曲线
Fig.7 Pressure fluctuation and frequency spectrum of P1 and P2 measuring points in the cavity
从图7(b)可见，压力脉动频谱具有多个峰值，其中前3个最优峰值频率被标记为f1=0.88、f2=1.88和f3=2.84(圆频率)。两观测点上前2个主导频率一致，但第3个主导频率在空腔中心位置较为显著，展示了腔内复杂的非线性涡系运动特性。结合第4.1节流场演化特征，分析可知主频f2与声反馈导致的剪切层涡运动相关，而频率f1与空腔后部运动较缓慢的大尺度涡运动相关，据信f3应是2种运动交互作用衍生频率。值得说明的是，相似的脉动主频在Dudley和Ukeiley[15]和Sheta等[16]的工作中也被发现过，表明超声速空腔流存在固有动力学特性。
此外，如图7(b)所示，虚线为Rossiter公式[4]预测的空腔Rossiter模态，其中设定常参数κν=0.66，相位移由经验公式α=0.062×L/D计算得到，本文方腔模型α=0.397。LES计算得到的3个压力脉动主导频率与Rossiter预测频率非常吻合，证明Rossiter公式对超声速空腔流动共振模态也具有预测能力。
4.3 压缩性效应
为明确压缩效应对空腔流动的影响，对Mach 0.6亚声速空腔流进行LES计算及对比分析。值得说明的是，除来流速度外，该算例其余来流参数以及计算网格均与Mach 1.4超声速空腔流一致。图8展示了亚声速空腔流动密度梯度幅值云图随时间的演化。由图8(a)可见，该流动仍然受剪切层涡结构运动主导，在空腔前缘位置卷起一列规则小涡。约在x1=2～5站位，细小涡积聚形成大尺度涡团结构。如图8(b)，包含大量小涡的大尺度涡团开始向空腔后缘撞击，破碎成大量的小湍流涡，同时形成向上游传播的声波。声波在空腔中沿腔内向上游传播，则与初始剪切层耦合形成声-流共振状态。声波向腔外辐射，形成上游远场占优的强声场，这与超声速空腔流中声波无法向上游传播不同。另外，由于没有激波出现，缺少了激波与剪切层相互干扰现象。
图8 亚声速(Ma=0.6)空腔流动密度梯度幅值云图
Fig.8 Density gradient magnitude of subsonic(Ma=0.6) cavity flow
图9为亚、超声速空腔流在空腔中心P2点(3，-0.5，0)的压力脉动及频谱。由图9(a)可见亚声速空腔流内压力脉动远小于超声速情形。若基于标准大气压条件将LES计算腔内无量纲压力转换为有量纲压力，即大气密度
环境声速a∞*=340.3 m/s，根据p*=p×ρ∞*(a∞*)2，并使用声压级公式Lp=20log(prms/pref)，其中prms为压力脉动均方根，参考声压pref=2×10-5 Pa。可计算亚声速空腔中心整体声压级为165.58 dB，而超声速空腔对应为179.11 dB，即低出约13.53 dB，可见压缩效应对空腔流动的噪声环境影响显著。此外，图9(b)频谱显示，亚声速空腔流具有一个明显占优的主频f=0.59，小于超声速空腔流的三个主导频率f1=0.88、f2=1.88和f3=2.84，意味着亚声速空腔流主要涡结构运动相对缓慢。与图8流场结合分析，发现亚声速空腔主频f=0.59由空腔声-流耦合现象导致的剪切层大尺度涡团结构运动引起，该频率也是声辐射的主频。整体而言，压缩性效应对空腔流动振荡及腔内噪声环境影响显著。
图9 亚/超声速空腔流压力脉动及其频谱曲线
Fig.9 Comparison of subsonic/supersonic cavity flow for pressure fluctuation and its frequency spectrum
5 结论
1) 超声速(Mach 1.4)空腔流动LES数值结果与文献中实验及数值结果对比，表明基于低色散低耗散计算气动声学技术的LES方法对空腔流计算的精确性、鲁棒性，解决了常规CFD方法或商用软件计算精度低不易捕捉噪声及边界非物理声反射等问题。
2) 超声速来流条件下，剪切层大尺度涡结构与空腔后壁发生撞击并产生经腔内向前传播的反射声波，该声波与空腔前缘初始剪切层相互作用，导致空腔内形成声-流耦合自持振荡的声反馈环现象。基于该认识，可设置专门装置在超声速战斗机武器舱等工程应用中破坏来流剪切层内大尺度涡结构的一致性，以削弱声-流耦合现象，降低腔内噪声水平。
3) 超声速空腔流展现出复杂非线性涡系运动导致多个特征频率。观测发现前3个频率较为显著，与流场演化历程结合分析可知第2个主频f2最优，其与声反馈导致的剪切层涡运动相关，而第1个主频f1与空腔后部较缓慢的大尺度涡运动相关，第3个主频f3是前2种运动相互作用衍生。
4) 当来流压缩效应减弱时，亚声速(Mach 0.6)空腔流动剪切层涡团与后缘撞击形成的声-流共振现象与超声速情形基本一致，但亚声速空腔流动的压力脉动强度明显减弱，且流动振荡主导频率减小，表明腔内涡运动缓慢。
[1] 李菁，陈帮，胡俊香.三维内埋式航弹与载机分离非定常流场数值模拟[J].兵器装备工程学报，2019，40(4):119-123.
LI Jing，CHEN Bang，HU Junxiang.CFD numerical simulation of the interference flow fluid of three-dimensional airplane and aviation bomb[J].Journal of Ordnance Equipment Engineering，2019，40(4):119-123.
[2] 王文益， 陈晨.配备GPS/INS 组合导航系统的无人机诱捕方法[J] .兵器装备工程学报， 2020， 41(11):212-217．
WANG W Y， CHEN C.Method of Trapping UAV Equipped with GPS /INS Integrated Navigation System[J] .Journal of Ordnance Equipment Engineering， 2020， 41(11):212-217．
[3] KRISHNAMURTY K.Acoustic radiation from two-dimensional rectangular cutouts in aerodynamic surfaces[R].Tech.Rep.CR-3487，NACA，1955.
[4] ROSSITER J E.Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[R].Tech.Rep.3438，Aeronautical Research Council Reports and Memoranda，1964.
[5] 叶坤，叶正寅，武洁，等.DMD和POD对超燃冲压发动机凹腔流动的稳定性分析[J].气体物理，2016，1(5):39-51.
YE Kun，YE Zhengyin，WU Jie，et al.Stability analysis of scramjet open cavity flow base on POD and DMD method[J].Physics of Gases，2016，1(5):39-51.
[6] SINHA N，ARUNAJATESAN S，UKEILEY L S.High fidelity simulation of weapons bay aeroacoustics and active flow control[C].AIAA Paper 2000-1968，2000.
[7] SMITH B R，JORDAN J R，BENDER E E，et al.Computational simulation of active control of cavity acoustics[C].AIAA Paper 200-1927，2000.
[8] RIZZETTA D P，VISBAL M R.Large-eddy simulation of supersonic cavity flowfields including flow control[C].AIAA Paper 2002-2853，2002.
[9] BR
S G A，COLONIUS T.Three-dimensional instabilities in compressible flow over open cavities[J].Journal of Fluid Mechanics，2008，599:309-339.
[10] SUN Y Y，ZHANG Y，TAIRA K，et al.Width and sidewall effects on high speed cavity flows[C].AIAA Paper 2016-1343，2016.
[11] TAM C K W，WEBB J C.Dispersion-relation-preserving finite difference schemes for computational acoustics[J].Journal of Computational Physics，1993，107(2):262-281.
[12] Tam CKW，Dong Z.Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[J].Journal of Computational Acoustics，1996，4(2):175-201.
[13] BERLAND J，BOGEY C，BAILLY C.Optimized explicit schemes:matching and boundary schemes and 4th-order runge-kutta algorithm[C].AIAA Paper 2004-2814，2004.
[14] BOGEY C，DE Cacqueray N，BAILLY C.A shock-capturing methodology based on adaptive spatial filtering for high-order non-linear computations[J].Journal of Computational Physics，2009，228:1447-1465.
[15] DUDLEY J G，UKEILEY L.Detached eddy simulation of a supersonic cavity flow with and without passive flow control[C].AIAA Paper 2011-3844，2011.
[16] SHETA E F，HARRIS R E，LUKE E A，et al.Hybrid RANS/LES acoustics prediction in supersonic weapons cavity[C].AIAA Paper 2015-0009，2015.