万振华，周 林，孙德军

（中国科学技术大学近代力学系，安徽 合肥230027）

0 引 言

一般认为长高比L／D＜6的方腔流动为开式流动，其在航空、航天以及其它工业领域应用广泛，如起落架舱、武器舱、敞篷汽车等。开式方腔流动会产生强烈的振荡，有可能引发结构疲劳损伤和噪声污染。振荡的诱因、幅值、抑制方法等都很值得研究，涉及流动失稳、波涡非线性作用、自持振荡等许多理论问题。故20世纪50年代以来，方腔流致振荡和发声问题就已经引起了众多研究者的兴趣［1］。Krishnamurty［2］最早从实验得到方腔声场结构，Rossiter［3］总结了预测振荡频率的半经验公式，Tam和Block［4］引入线性稳定性分析提出了新的频率预测模型，Rowley等［1，5－6］对方腔的空间稳定性、降维模型、理论模型、控制等问题进行了一系列研究，Shieh［7］等对高Reynolds数下不同构型方腔进行了湍流数值模拟研究。衣云峰［8］研究了低亚声速突出和陷落式圆柱空腔流动，并归纳了频率公式。罗柏华等［9］实验研究了纯音激励方式的噪声抑制。此外，罗柏华［10］、侯中喜等［11］、李晓东等［12］等使用不同数值方法都成功模拟出了方腔流场的振荡、声波传播性质。

虽然国内外研究内容已经很多，但仍然有很多问题有待深入研究，如目前对方腔的振荡存在两种认识［1，6］：一种认为它是自持的，另外一种认为它是轻微衰减的系统，所以具体振荡性质需由实际情况具体分析。关于方腔振荡的低频成分也有几种不同产生方式，如涡－后拐角撞击方式变化导致低频产生［13］，湍流间歇结构导致模态切换［14］引入低频等。此外，目前对回流区与剪切层相互作用研究较少。

本文采用直接数值模拟方法，研究了低Reynolds数、亚声速开式方腔流动的振荡及噪声问题，特别是回流区与剪切层相互作用及其与低频的关系。首先，研究了来流边界层厚度的影响，探讨了所涉及参数范围内的方腔流动振荡机制；针对来流边界层变薄后出现的低频成分，从回流区与剪切层相互作用出发分析了其产生的机制。其次，采用本征正交分解（POD）方法，分析了回流区与剪切层相互作用对剪切层涡结构的影响，以及涡结构改变与涡－后拐角撞击方式的关联。最后，计算得到了方腔的声场结构，通过胀量场的分析得到了声波产生和传播的过程。

1 数值方法

1.1 计算方法

采用直角坐标系下二维可压缩Navier－Stokes方程：

考虑到捕捉噪声，所以空间使用了的6阶紧致格式离散［15］，时间方向使用低耗散和低色散的Runge－Kutta格式［16］，粘性项使用的6阶中心差分格式。方程以远场声速a∞进行无量纲化。Reynolds数定义为Re＝ρ∞a∞D／μ∞，取Prandtl数Pr＝0.72。粘性系数以Sutherland公式给出，其中远场温度T∞＝300K。

计算模型如图1所示，方腔的高度定义为D，长度为L，方腔前端延伸到－6D，尾端延伸到13D，上端延伸到15D。边界条件设定如下：方腔上部左、右端分别为无反射亚声速出、入口，上方设定为无反射出口，固壁边界设为等温壁［17］。在方腔上部三个方向施加海绵层边界条件［18］以抑制边界反射。另外，采用了10阶中心滤波抑制数值伪振荡。初始条件给定如下：腔体上部以基本流速度型进行初始化；腔体内部速度设为0。密度和压力设为远场条件。

图1 方腔设定与计算域的示意图Fig.1 Schematic diagram of cavity configuration and computational domain

1.2 计算参数选择

本文主要研究高亚声速低Reynolds数（相对于湍流发生的Reynolds数）方腔流动振荡和发声问题。选择典型的方腔长高比为L／D＝2，Re＝2500，同时也便于和已有实验结果进行比对。来流为层流状态，研究了不同来流马赫数M、来流边界层动量厚度θ的影响，具体计算参数如表1所示。

表1 算例列表Table 1 Computational cases

1.3 网格收敛性验证

本文使用的程序由文献［19－20］中程序发展而来，涉及的空间格式、边界条件等均已很好验证，如采用本文中的数值算法计算文中［21］的Taylor－Green涡算例，条件设定与文献［21］相同，在13×13的网格上求解得到了Taylor－Green涡速度随时间的衰减（如图2），无量纲时间t＝1／（8π2）时的计算解与理论解的相对误差小于0.03%。考虑到声学问题对网格分辨率要求很高，所以为了确保文中使用的网格满足要求，本文做了网格收敛性分析。计算采用单向拉伸网格，针对L／D＝2的情形，文中使用的计算网格为：腔外网格为918×416，腔内网格为262×118，总网格点数约为40万，与文献［5］网格数相当，网格在固壁附近和唇口处最密，唇口两端靠近壁面（x方向）的网格间距Δx／D≈0.005，在－0.1＜y／D＜0.1区间内集中了32层网格，最小的网格间距Δy／D≈0.005。再考虑一套加密的网格：腔外网格为1118×510，腔内网格为308×152，总网格点数约为60万。比较算例1在两套网格下后拐角位置的密度随时间的变化。如图3（a）所示，二套网格下密度的变化基本相同，相应峰值差别很小，且相位基本保持同步。为了更精确地比较频谱误差，图3（b）给出了前200个无量纲时间密度的FFT变换后在谱空间上幅值。从图中可以看出在各种频率下，二者幅值相差很小，如在主频f＝0.216上幅值的误差小于1×10－3，故说明当前选择的计算网格密度已经足够。

图2 Taylor－Green涡的速度随时间衰减Fig.2 The decaying of velocity for a Taylor－Green vortex

图3 网格收敛性验证Fig.3 Validation of mesh convergence

2 计算结果及讨论

2.1 流场的振荡现象

目前，对方腔流场振荡机制的认识存在两种观点。一种观点认为，振荡是自持的：腔体内外流体的剪切诱导Kelvin－Helmholtz不稳定性，产生涡卷起；涡卷起后向下游运动，进而与后拐角撞击压缩产生往前传的声波；声波到达方腔前端会再次激发剪切层失稳，使新涡卷起，完成自持振荡过程。所以，在无任何外部持续激励时，方腔流动仍会持续地振荡。此时方腔流动的动力学系统可认为是一个关于不稳定平衡点（Navier－Stokes方程定常解）的稳定的极限环［5－6］，振荡的幅值取决于非线性效应，如剪切层失稳中扰动增长达到饱和。还有一种观点认为，方腔流动是稳定的、轻微衰减的系统，只有在外界持续激励下（如边界层内湍流脉动，壁面粗糙引起扰动等）才能维持振荡［1，6］；撤除激励，振荡现象也将消失。

本文的计算中没有施加任何外部激励，在所计算的参数范围内得到了流场的振荡现象，计算结果支持第一种观点，即属于自持振荡。振荡的强弱可以采用声压级来表征（图4），结果表明，边界层动量厚度θ越小，则方腔底部声压级越高，即振荡幅值越大。这与剪切层的不稳定性特性相吻合，即剪切层θ减小导致扰动的空间增长率增加。所以，振荡幅值与θ的变化关系主要与剪切层的空间不稳定性增长率有关，这种剪切层的对流不稳定性与下游方腔中的反馈机制相结合，仍然可以形成自持的振荡，而无需来自外界持续扰动的激励。

Rossiter［3］于1964年提出了预测振荡频率的半经验公式，其形式如下：

图4 方腔底部的声压级（M＝0.6）Fig.4 Sound pressure levels of cavity floor at M＝0.6

其中，κ＝0.57是一个经验常数，α＝0.25为涡撞击后拐角到产生声波的相位滞后，n为剪切层中涡的数目，对应模态的阶数。n＝1，2的模态分别称为Rossiter I、II模态。图5给出了不同马赫数下得到的振荡频率，在本文的计算参数范围内只得到了Rossiter II模态，这与Krishnamurty［2］的实验结果是基本符合的。本文的频率与文献［5］中Rossiter II模态的频率接近，但都略低于实验［2］得到的频率，原因可能是实验中的Reynolds数更高和来流边界层厚度不同造成的。此外，实验中同一马赫数下出现多个频率，这可能是由于来流边界层厚度不同造成的频率偏移。在本文中这种偏移很小，而在湍流模拟中这种频率偏移现象明显［12］，原因可能是湍流情况下频率不再是明确的大尺度涡结构对流主导，导致对边界层厚度改变更加敏感。

图5 不同马赫数下的振荡频率Fig.5 The oscillation frequencies at different Mach numbers

当M＝0.6，由图6（a）和图6（c）中的相图可知，算例1的振荡是周期的；随着来流边界层厚度的减小，算例2则为准周期的。相应地，从图6（b）和图6（d）可以看出，算例1中只存在对应Rossiter II模态的单一频率St＝0.72；而算例2中存在两个频率，St1＝0.72的基频能谱占主导地位，仍和Rossiter II模态一致，但同时出现了低频成分St2＝0.36，视为亚谐频。虽然亚谐频非常靠近Rossiter I模态的频率，但其不符合Rossiter关于I模态中只有单涡的定义，故我们认为这个低频不是I模态。它的产生与Rockwell等的实验结果［13］相似，Rockwell等认为这是因剪切层中涡－后拐角撞击方式变化造成的。Rockwell的实验［13］是在水洞中进行的，本文计算的低Reynolds数、高亚音速情况的气体流动，结果表明，涡－后拐角撞击不稳定仍然存在，它们是回流区与剪切层相互作用形式不稳定的结果。

图7（a）给出了算例1的时均流场结构，从图中可以清晰地区分几个主要区域：主回流区A，A区诱导的主反向流动区域B，次回流区C，A区诱导的次

反向流动区D。回流区流动可以改变剪切层中涡心轨迹和横向运动速度。图7（b）为涡心运动轨迹，可以将其大致分为4个区域：I区为自由发展区，涡心横向运动速度以及纵向位置没有改变。II区涡受次回流区C后端牵引运动方向会向方腔内偏斜，同时由于越靠近方腔横向对流速度越小，涡横向运动速度会减小。III区涡进入主回流区，涡会被主回流区A前端牵引上升，上升越高横向对流速度越大，涡的横向运动速度越快。IV区涡受A后端牵引位置逐渐下降，涡横向运动速度减小。算例1中的剪切层涡都经过上述4个区域，最终与后拐角发生稳定的撞击，即每次涡－后拐角的撞击方式都相同，如图7（c）所示，整个流场也呈Rossiter II模态形式的周期性振荡，振荡周期由相邻两次涡－后拐角撞击的时间间隔决定。

图6 （1.9D，0）位置的相图与功率谱Fig.6 The phase diagram and power spectrum of perturbed pressure at the location（1.9D，0）

当来流边界层厚度减小（算例2），由于剪切层的不稳定性增强，剪切层与腔内回流区的相互作用也明显增强，导致剪切层涡呈高低交替出现，它们分别走过高、低两条轨迹，每条轨迹也可大致分成上述的4个区域，如图7（d）。由于此时的流动是准周期的，故各个周期内涡的运动轨迹还有细微的差异，但都接近于图7（d）所示的高、低涡运动轨迹。

剪切层也可以通过涡量输运的方式改变回流区的运动。现仍结合算例2进行分析。如果剪切层中前一个涡走低轨道，那么它将与后拐角发生图7（e）形式的撞击，其特点是撞击后大部分涡量进入回流区，使得主回流相对变强，此即Rockwell等［13］实验中的PC（Partial Clipping）形式的撞击；主回流的加剧会使得剪切层中下一个涡走高轨道，它与后拐角将发生图7（f）形式的撞击，其特点是撞击后只有少量涡量进入回流区，主回流相对变弱，此即Rockwell等［13］实验中的PE（Partial Escaping）形式的撞击，主回流减弱导致下一个涡将走低轨道。如此周而复始，形成准周期流动。由于此时两次相邻涡的撞击不再是周期的，只有完成了一个PC到PE形式涡－后拐角撞击过程才构成一个主周期，所以，这种撞击方式的切换导致了新的低频（亚谐频）成分的产生。相对应的是，算例1中的涡撞击形式介于PC与PE之间，且每次都相同，形成周期性流动，因此没有低频成分出现。

值得指出的是，前述Rossiter频率预测公式（式2）是假设恒定的剪切层涡对流速度得到的，而回流区与剪切层的相互作用会改变涡横向运动速度，会影响Rossiter公式的准确性。

算例3和算例4的结果类似，当来流边界层较厚时，是没有低频的Rossiter II模态周期性流动，当来流边界层变薄，会有低频成分出现，不再赘述。

2.2 本征正交分解（POD）分析

为了进一步明晰回流区与剪切层相互作用对涡模态的影响，采用了等熵可压缩的POD方法［22］进行了分析。由表2可知，算例2的能谱比算例1更宽，表明薄来流边界层的相互作用更强。算例1中前两个模态就捕捉了超过90%的能量，算例2中大部分能量需要前4个模态才能捕捉。算例1中的模态1和2的动力学结构基本类似，故图8只给出了模态1和3。算例2中模态1和2，模态3和4的结构类似，因此图9中只给出了模态1、3和5。可以看出，低阶模态中，剪切层中均存在2个大涡结构，对应于功率谱中Rossiter II模态，故低频成分不对应于Rossiter I模态。图8（b）和图9（c）中涡数目增多，涡结构变小，对应于功率谱中高频成分，但能量很小。图8（a）和图9（a）的涡结构类似，涡纵向结构比较宽，这种涡结构倾向于PC形式的撞击，图9（b）中的涡结构则不同，比较扁平，更容易形成PE形式的撞击。比较图8和图9还可以看出，算例2主回流区中涡结构要比算例1的更加复杂和紊乱。总体而言，POD模态验证了低频成分不是Rossiter I模态，回流区与剪切层的相互作用会影响涡结构的形态，从而改变涡－后拐角的撞击方式。

图8 算例1的POD模态的涡量场ω∈（－5，5）Fig.8 Vorticity filed of different POD modes for Case 1

图9 算例2的POD模态的涡量场ω∈（－5，5）Fig.9 Vorticity field of different POD modes for case 2

表2 能量分数Table 2 Energy fractions

2.3 声波结构对比及分析

本文中的声波结构通过胀量场和密度梯度两种方式分别进行了表征。图10（a）～图10（c）为Krishnamurty［2］的实验结果，声场由密度的纹影图给出。图11（a）～图11（c）为本文DNS模拟得到的结果。实验与计算均使用层流来流边界层，且都是被Rossiter II模态主导，所不同的是实验的Reynolds数高出一个量级。由二者结果比较可知，虽然计算中Reynolds数较低，但二者声场的定性结构，相位关系都基本相符。此外，从实验和计算的结果均可看出，随着马赫数的增加，方腔外声波的传播方向会发生明显变化，且两道声波之间相位差会变大，这些差别主要来自于声波传播过程中的对流效应，波阵面的形态取决于波速和当地流动速度的叠加（见图12）。故在M＝0.6时，声波传播方向与流动夹角约为135°角，随着马赫数增加这个夹角逐渐减小。

图12描述了声波产生和传播的具体过程：首先，涡－后拐角的作用产生声波，声波分别在腔外和腔内传播（表示为波阵面I和II）。其次，因腔内对流速度较小，波阵面II传播快于I，二者产生相位差。再次，腔内的波阵面II到达前壁面，部分透过剪切层辐射出方腔，剩下部分发生反射。最后，可观察到腔外有两道声波I、II，腔内有II的反射声波。声波I、II在腔外相位差大小被腔体内外的对流速度的差异决定，故在算例1中声波I、II相位差较小，而在算例4中声波I、II相位差明显。图13以胀量场的形式给出了算例4中声波传播的具体过程，进一步验证了图12的描述的过程。值得注意的是，本文的计算很好地验证了Tam等人［4］在理论建模时关于声波传播的描述，不同是计算中未能看到明显的声波II在腔体下壁面的反射声波，可能原因是反射声波相对其它声波强度较弱，在旋涡中又易耗散，造成计算中不容易捕获。

图12 波阵面传播过程示意图Fig.12 Sketch of propagation of wave fronts

图13 算例4的4个不同时刻的胀量场Fig.13 4snapshots of dilatation field for case 4

3 结 论

本文采用直接数值模拟方法，研究了低Reyn－olds数下L／D＝2的方腔的亚声速流致振荡现象及其诱导的噪声。在所计算的参数范围内，方腔流场的振荡是自持的，并且由Rossiter II模态主导。振荡幅值、频率与来流边界层厚度密切相关，当来流边界层变薄，振荡幅值增大，周期振荡会转变为准周期振荡，并产生低频成分。低频成分产生的根本原因是回流区与剪切层的相互作用。这种相互作用会改变剪切层涡的结构、运动轨迹及速度，进而导致了不同的涡与后拐角撞击方式；而撞击方式的切换构成了低频成分。采用本征正交分解（POD）方法，分析了不同振荡形式所对应的本征模态涡结构。对于没有低频成分的周期振荡情况，只发现一种能量集中的本征模态；而对于出现低频成分的准周期振荡情况，发现了两种能量较为集中的本征模态，它们反映了回流区与剪切层相互作用对涡结构的影响，并与后拐角的不同撞击方式有关。声场的计算表明，涡－后拐角的撞击作用会产生声波，并分别在腔外和腔内传播，形成两个波阵面，它们之间的相位差与对流效应有关。

［1］ROWLEY C W，WILLIAMS D R.Dynamics and control of high－Reynolds－number flow over open cavities［J］.AnnualReviewFluidMechanics，2006，38：251－276.

［2］KRISHNAMURTY K.Sound radiation from surface cutouts in high speed flow［D］.California Institute of Technology，1956.

［3］ROSSITER J E.Wind－tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds［R］.Aeronautical Research Council Reports and Memoranda，Technical Reports 3438，1964.

［4］TAM C K W，BLOCK P J W.On the tones and pressure oscillations induced by flow over rectangular cavities［J］.JournalofFluidMechanics，1978，89（2）：373－399.

［5］ROWLEY C W，COLONIUS T，BASU A J.On self－sustained oscillations in two－dimensional compressible flow over rectangular cavities［J］.JournalofFluidMechanics，2002，455：315－346.

［6］ROWLEY C W，WILLIAMS D R，COLONIUS T，et al.Linear models for control of cavity flow oscillations［J］.JournalofFluidMechanics，2006，547：317－330.

［7］SHIEH C M.Parallel numerical simulations of subsonic，turbulent，flow－induced noise from two－and threedimensional cavities［D］.Pennsylvania State University，2000.

［8］衣云峰，何祚镛.圆柱形腔流激振荡及其耦合共振的研究［J］.声学学报，1996，21（4）：339－456.

［9］罗柏华，胡章伟.流动诱导空腔振荡及声激励抑制的实验研究［J］.南京航空航天大学学报，1996，28（3）：331－336.

［10］罗柏华.二维高亚音速空腔流激振振荡的数值模拟研究［J］.空气动力学学报，2002，20（1）：84－88.

［11］侯中喜，易仕和，王承尧.超声速开始空腔流动的数值模拟［J］.推进技术，2001，22（5）：400－403.

［12］李晓东，刘靖东，高军辉.空腔流激振发声的数值模拟研究［J］.力学学报，2006，38（5）：599－604.

［13］ROCKWELL D，KNISELY C.Vortex－edge interaction：mechanisms for generating low frequency components［J］.PhysicsofFluids，1980，23（2）：239－240.

［14］GLOERFELT C，BOGEY C，BAILLY C.Numerical evidence of mode switching in the flow－induced oscillations by a cavity ［J］.InternationalJournalof Aeroacoustics，2003，2（2）：193－218.

［15］COLONIUS T，LELE S，MOIN P.Sound generation in a mixing layer［J］.JournalofFluidMechanics，1997，330：375－409.

［16］BERLAND J，BOGEY C，BAILLY C.Low－dissipation and low－dispersion fourth－order Runge－Kutta algorithm［J］.ComputersandFluids，2006，35：1459－1463.

［17］POINSOT T J，LELE S.Boundary treatment for direct simulations of compressible viscous flows［J］.Journal ofComputationalPhysics，1992，101：104－129.

［18］FREUND J B.Proposed inflow／outflow boundary conditions for direct computation of aerodynamic sound［J］.AIAAJ.，1997，25（4）：740－742.

［19］WEI M J.Jet noise control by adjoint－based optimization［D］.University of Illinois at Urbana－Champaign，2004.

［20］周林，万振华，魏明骏，等.亚声速时间混合层的噪声研究［A］.第十三届全国分离流、旋涡和流动控制会议论文［C］，2010：220－224.

［21］REMBOLD B.Direct and large－eddy simulation of compressible rectangular jet flow［D］.Swiss Federal Institute of technology，2003.

［22］ROWELY C W，COLONIUS T，MURRAY R M.Model reduction for compressible flows using POD and Galerkin projection［J］.Physcia.D.，2004，189：115－129.