亚音速圆口喷流的大涡模拟

林 健, 赖焕新

(华东理工大学承压系统与安全教育部重点实验室,上海 200237)



亚音速圆口喷流的大涡模拟

林 健, 赖焕新

(华东理工大学承压系统与安全教育部重点实验室,上海 200237)

采用大涡模拟方法研究了马赫数M=0.75,雷诺数Re=50 000的圆口可压缩喷流问题。在与实验数据比对验证的基础上,分析了喷流剪切层不稳定性发展的两个主要阶段,着重考察了剪切层中轴向速度脉动u′、径向速度脉动v′、压力脉动p′在空间任意两点上的时空相关系数随时间和空间发展分布的特点,并由此分析剪切层中旋涡脉动的迁移特性,为进一步研究与控制喷流噪声提供基础。

亚音速; 喷管喷流; 大涡模拟; 剪切层; 相关性分析

亚音速喷流噪声是现代民航客机中主要的噪声来源之一。随着世界各国对民航客机噪声适航标准与要求的日益苛刻,进一步深入研究喷流流场中的涡运动与发声过程,从而为研究和改善大型民航客机噪声水平提供可能,依然是当今民用航空噪声控制的必由之路。

自从Lighthill在1952年创立气动声学[1-2]这一学科以来,喷流噪声一直是气动声学的基本问题,声源机理的解析是问题的关键。Lilley[3]总结了自Lighthill至上世纪末共50多年中喷流噪声预测的Lighthill方程及各种版本的声比拟理论方程,指出这些方程中的源项本身都是未知量,必须结合声源流场的准确预测。Tam[4]指出具有物理意义的声源是流场中的旋涡,而声比拟中的所谓四极子、偶极子等仅是用于数学上重构某一给定声场的等价声源,它们并不具有物理意义。Tam[4-5]在总结分析了大量实验数据之后,进一步提出喷流中存在两类噪声源,前者是大尺度涡结构,后者是湍流的小尺度脉动,它们各自可用一种相似频谱来描述。但由于实验测量的困难,喷流流场中的涡运动迄今并未被完全清楚地认识。例如对于喷管出口下游处极薄的剪切层尚无有效的测量方法,而该区域的流动状态对于喷流噪声又非常关键。

本文选取Jordan[6]和Andersson[7]研究的亚音速喷管喷流为对象,采用大涡模拟(LES)方法进行计算分析,在方法验证的基础上,着重考察分析剪切层中不稳定性发展过程。在此基础上,分析了两点两时刻流场中脉动速度和脉动压力的相关系数,讨论它们随时间和空间发展分布的特点,为下一步分析亚音速喷流的噪声场提供基础。

1 亚音速喷管喷流计算

1.1 物理模型及网格划分

本文计算的喷管几何尺寸与流动参数来源于文献[6-7]。其中喷管出口直径Dj为50 mm,出口马赫数为0.75。计算域如图1所示,喷管外长度取为90Dj,上游宽度为20Dj,下游宽度为40Dj。整个计算区域均采用结构化网格,节点总数为1 640 854。对剪切层区域进行加密处理,喷管内部采用O型网格划分,如图1所示。

1.2 控制方程与边界条件

为对喷管中的可压缩流动进行大涡模拟,采用Favre[8]过滤得到的LES控制方程组:

连续性方程:

(1)

(2)

(3)

图1 计算域与网格划分Fig.1 Computational domain and mesh

(4)

(5)

(6)

为使方程组封闭补充状态方程:

(7)

式中R=287.1 J/(kg·K)。

本文选用动态Smagorinsky-Lilly模型,计算亚格子尺度应力采用式(8):

(8)

(9)

(10)

本文先采用RNGk-ε湍流模型进行稳态计算,当整体质量流量误差小于1%时认为达到收敛,将所得结果作为LES的初始流场。计算时选用Density-Based求解器,流体介质设置为理想气体,黏度采用Sutherland模型。瞬态时间步长设置为5×10-7s,确保计算域内90%区域迭代步长CFL数小于0.5。本文从初场开始,计算时长约130(Dj/Uj),此时非稳态流动沿时间已充分发展,可排除初场的影响,在此基础上再存取数据进行统计分析。

上游喷管入口设置:给定总压147 152Pa,总温320.4K,同时给出参考静压101 325Pa; 喷管外计算域入口以及四周环境界面设置为无反射的压力边界条件,给定总压101 325Pa,总温288K和参考静压101 325Pa;下游出口设置为无发射的压力出口边界条件,给定静压101 325Pa,总温288K; 喷管固壁采用无滑移、绝热壁面边界条件。

2 计算结果及分析

2.1 数值计算与实验结果的比较

为方便表述,对本文的符号作出如下说明:〈…〉t表示时间平均; 〈…〉θ表示轴向平均; Uj为喷管出口速度; u和u′分别表示轴向速度和轴向速度脉动,v和v′分别表示径向速度和径向速度脉动; uc为中心线轴向速度; 坐标原点为喷嘴出口圆心。

图2示出了本文计算所得喷管中心线时均轴向速度,图中给出了计算结果与文献[6]实验值、文献[7]计算值的分布对比。中心线速度〈uc〉t≥0.95Uj时对应的距离为势流核长度。从中可以看出,本文结果与文献[7]计算结果基本一致,预测势流核长度均比实验值偏短,这主要与亚格子模型的选取和网格质量有关。

图3示出了在喷口下游x/Dj为1、2.5和5这3个位置处轴向速度沿径向分布对比,可以看出在x/Dj为1和2.5这两处,本文计算结果更接近实验值; 但在x/Dj=5处文献[7]计算值稍好,可能是由于本文采用的网格在下游位置比文献[7]的网格较少导致的。

图2 喷流中心线上轴向速度的分布Fig.2 Axial velocity along the centerline

图3 x/Dj为1、2.5、5处轴向速度沿径向分布Fig.3 Radial profiles of axial velocity at axial positions x/Dj=1,2.5,5

图4给出了不同轴向位置上湍流脉动的径向分布,可以看出当前结果和实验值基本吻合,整体上与文献[7]的计算结果一致,表明本文的计算方法在模拟亚音速喷流时是可行的。

2.2 断面流速分布的相似性

图5示出了在喷管自由喷流轴截面上,各个不同断面〈u〉tθ/Uj的径向分布。从图中可以看出,不同断面轴向流速的分布显示出一定程度的相似性,即轴线上流速最大,往下游发展流速逐渐减小,距离喷管出口越远流速衰减越慢。为进一步分析断面流速分布的相似特性,定义剪切层厚度δ[10]为

(11)

式(11)中r0.10、r0.95以及下文中的r0.5分别表示不同断面上轴向流速为10%〈uc〉t、95%〈uc〉t和50%〈uc〉t(即中心线轴向流速的10%、95%和50%)处对应的半径值。

图6中离散点为本文剪切层厚度的计算结果。由于数据量太大,本文仅选取了7个位置的数据结果。使用最小二乘法对这些离散点进行线性拟合,延长该拟合直线使之与x轴相交于x0,即为喷流虚源。

将流速〈u〉tθ和断面上径向坐标r分别以量纲为一坐标〈u〉tθ/〈uc〉t和(r-r0.5)/(x-x0)来表示,则所有断面上量纲为—流速均落在一定的区域之内,如图7所示。可以看出本文计算结果基本上已经体现了自由喷流断面流速分布相似性的这个重要特性。但可以看出图中量纲为—流速一致性还有待改善,这主要是受到数据样本容量的限制以及虚源定位方法准确性限制[11]所致。

图4 不同轴向位置上湍流脉动的径向分布Fig.4 Radial profiles of turbulent fluctuations at axial positions x/Dj=1,2.5 and 5

图5 不同轴向断面上〈u〉tθ/Uj的径向分布Fig.5 Radial profiles of 〈u〉tθ/Uj at different axial position

图6 剪切层厚度δ/Dj沿轴向坐标的分布Fig.6 Shear layer thickness δ/Dj with axial coordinate

图7 〈u〉tθ/〈uc〉t随(r-r0.5)/(x-x0)分布Fig.7 Variation of 〈u〉tθ/〈uc〉twith (r-r0.5)/(x-x0)

2.3 剪切层不稳定性的发展过程

剪切层不稳定性的发展过程直接影响着喷流流场的湍流特性,对于亚音速喷流噪声机理的研究具有重要的意义。许常悦等[12]通过分析涡量厚度δΩ沿剪切层的分布来对该过程进行了分析。Wan等[13]用喷管半值速度时的半径r0.5和涡量厚度δΩ对自由流喷管和旋流式喷管进行了对比研究。两者的物理意义基本相似,均表示剪切层增长率。

图8示出了喷管r0.5沿轴向坐标的分布(rj为喷管出口半径),图中曲线主要呈现出两个斜率:第1阶段(515)斜率快速增长。为进一步分析该现象,图9示出了两个时刻(t1=223Dj/Uj和t2=224Dj/Uj)涡强Ω0的瞬时分布,Ω0=‖×V‖。

图8 喷管半值速度宽r0.5沿轴向坐标的分布Fig.8 Half-width of jet r0.5 along x

从图9可以看出当515时,剪切层开始混合、扰动涡合并,导致了曲线第2个阶段斜率快速增长的变化。

2.4 时空关联分析

远场噪声的强度及能谱与声源中任意两点两个时刻的时空相关系数有密切关系[14]。如图10所示,对于剪切层中轴向速度的脉动值,点x与其下游点x+ξ处的两点时空相关系数Ruu定义如下:

(12)

式中:x为测点在流场中的轴向坐标;ξ为两个测点间的距离;τ为延迟时间。时空相关系数越大,表明两个现象相关的程度越高。当两个测点在湍流的同一个旋涡内时,湍流脉动速度的相关程度就高,对应相关系数就越接近于1。反之,两个脉动值的相关程度就低,时空相关系数会趋于零[15]。

图11示出了剪切层中r=0.5Dj,x=2.5Dj处的点与其下游x+ξ处的点两点间的时空相关系数分布。当ξ=0时,对应图中相关系数从最大值1起始的一条曲线,该曲线表示x=2.5Dj位置上t时刻与t+τ时刻的两个脉动的相关系数随τ逐步变化的关系。图中其他每一条连续曲线均代表固定点x=2.5Dj处与x+ξ(ξ以0.1Dj递增)处的时空相关系数。可以看出时空相关系数的峰值随着ξ的增加而下降,表明两点之间脉动的相关程度随着空间距离的增大而降低。当空间距离增大时,旋涡在向下游发展的过程中不断地卷入周围的流体介质,互相之间发生质量、动量和能量的交换,并伴随着能量的耗散,因此两个脉动现象之间的相关关系逐步被解除。

图9 涡强的瞬时分布Fig.9 Instantaneous vortical shown by Ω0

图10 剪切层中脉动相关系数定义的对应位置x和测点间距离ξ的示意图

Fig.10 Two-point space-time correlations in the shear layer are obtained for locationxusing a spatial separationξ

图12示出了喷流剪切层中r=0.5Dj,x=2.5Dj处的两点时空相关系数Ruu随ξ和τ的二维分布图。图中等值线表示Ruu,点划线是Ruu极大值的连线,其斜率表征了旋涡的迁移速度UA。

图11 剪切层中r=0.5Dj,x=2.5Dj处的轴向 速度脉动时空相关系数Fig.11 Two-point space-time correlation of axial velocity fluctuation at a location x in the shear layer r=0.5Dj,x=2.5Dj

图12 剪切层中r=0.5Dj,x=2.5Dj处Ruu云图Fig.12 Contours of Ruu obtained in the shear layer ofr=0.5Dj,x=2.5Dj

由于旋涡在迁移过程中发生质量、动量和能量的变化,因此理论上只有在原点处的斜率才等于当地的旋涡迁移速度UA[16]:

(13)

式中τRmax表示两点时空相关系数取极值时的τ值。本文在图11中处理时空相关系数时,ξ的增量取值为0.1Dj,因此近似采用Δξ=0.1Dj估算喷管唇线上x=2.5Dj、5Dj及10Dj3个点上的旋涡迁移速度,结果列于表1。

由表1可看出本文计算结果与文献[7]计算结果基本一致。距离喷管出口x/Dj<5内,UA趋于恒定,为0.65Uj。随着流动向下游发展和剪切层的增长变厚,使得旋涡的迁移速度降低,在x/Dj=10时UA减小为0.53Uj。计算结果略低于文献[7]的0.56Uj,可能是本文的计算网格较少造成的。

表1 剪切层喷嘴唇线(r=0.5Dj)上各点的迁移速度UATable 1 Convection velocity UA at different positions in the shear layer r=0.5Dj

在可压缩流动中,声波的传递与压力的脉动值p′密切相关,为此本文进一步对压力脉动p′、径向流速的脉动v′、以及轴向速度脉动u′之间的时空相关进行了分析。参照Ruu的定义式(12),定义以下时空相关系数Rvv、Rpp、Rup、Rvp分别如下:

图14示出了剪切层中r=0.5Dj,x/Dj为2.5、3和5这3处Rup与Rvp的对比。可以看出Rup与Rvp基本同相位,但各峰值点处Rvp数值均高于Rup,显示了p′与v′之间存在较强的相关关系。

由于Lighthill声源的主要贡献来自轴向速度,压力项对总声源贡献较小[14]。所以v′增大虽然会伴随出现p′增大,但并不会对总声源造成明显的增强,因此在考虑降噪措施时,如果某一种方法能适当地增大v′而减小u′,则可减弱主要噪声源的贡献,该方法将有可能是一种可行的降噪措施。

图13 剪切层中r=0.5Dj,x=2.5Dj处Rvv、Rpp云图Fig.13 Contours of Rvv、Rpp obtained in the shear layer (r=0.5D,x=2.5Dj)

图14 剪切层中不同位置处压强脉动与速度脉动相关系数Fig.14 Two-point space-time correlation of pressure and velocity fluctuation obtained in the shear layer

3 结 论

本文采用大涡模拟的方法对马赫数为0.75的亚音速圆口喷管喷流进行研究,得到如下结论:

(1) 亚音速喷流剪切层不稳定性的发展过程主要分为两个阶段:第1阶段扰动涡开始失稳并进行配对,剪切层厚度缓慢增加; 第2阶段剪切层开始混合、扰动涡合并,剪切层厚度快速增长。

(2) 剪切层中轴向速度脉动u′、径向速度脉动v′和压力脉动p′的时空相关分布规律类似,而p′的迁移速度比u′和v′的迁移速度要快。

(3) 剪切层中Rup与Rvp基本同相位,p′与v′之间存在较强的相关关系。如果某一种方法能适当地增大v′而减小u′,则该方法将有可能是一种可行的降噪措施。

[1] LIGHTHILL M J.On sound generated aerodyne-mically:I.General theory[J].Proceedings of the Royal Society,A:Mathematical,Physical and En-gineering Sciences,1952,211(1107):564-587.

[2] LIGHTHILL M J.On sound generated aerodyne-mically:II.Turbulence as a source of sound[J].Proceedings of the Royal Society：A.Mathematical,Physical and Engineering Sciences,1954,222(222):1-32.

[3] LILLEY G M.Jet noise and classical theory and experiments[C]//Hubbard H H.Aeroacoustics of Flight Vehicles:Theory and Practice.Hampton,Virginia:NASA,1991:211-289.

[4] TAM C K W,VISWANATHAN K,AHUJA K K,etal.The sources of jet noise:Experimental evidence[J].Journal of Fluid Mechanics,2008,615:253-292.[5] TAM C K W.Jet noise:Since 1952[J].Theoretical and Computational Fluid Dynamics,1998,10(14):393-405.

[6] JORDAN P,GERVAIS Y,VALIERE J C,etal.Final results from single point measurements[R].[s.l.]:[s.n.],2002.

[7] ANDERSSON N,ERIKSSON L E,DAVIDSON L.Largeeddy simulation of subsonic turbulent jets and their radiated sound[J].AIAA Journal,2005,43(9):1899-1912.

[8] VREMAN B.Direct and largeeddy simulation of the compressible turbulent mixing layer[D].The Netherlands:University of Twente,1995.

[9] LARCHEVEQUE L,SAGAUT P,MARY I,etal.Large-eddy simulation of a compressible flow past a deep cavity[J].Physics of Fluids,2003,15(1):193-210.

[10] ARAKERI V H,KROTHAPALLI A,SIDDAVARAM V,etal.On the use of microjets to suppress turbulence in a mach 0.9 axi-symmetric jet[J].Journal of Fluid Mechanics,2003,490:75-98.

[11] Tide P S,Babu V.Numerical predictions of noise due to subsonic jets from nozzles with and without chevrons[J].Applied Acoustics,2009,70(2):321-332.

[12] 许常悦,吴丹,孙建红.可压缩湍流圆孔射流的大涡模拟[J].南京航空航天大学学报,2010,42(5):583-587.

[13] WAN Z H,ZHOU L,YANG H H,etal.Large eddy simulation of flow development and noise generation of free and swirling jets[J].Physics of Fluids,2013,25(12):1-27.

[14] 李栋.可压缩湍流时空关联与混合层噪声的数值研究[D].北京:中国科学院大学,2013.

[15] 张强.气动声学基础[M].北京:国防工业出版社,2012.

[16] 张兆顺,崔桂香,许春晓.湍流理论与模拟[M].北京:清华大学出版社,2005.

Large Eddy Simulation of Subsonic Round Jet

LIN Jian, LAI Huan-xin

(Key Laboratory of Pressurized Systems and Safety,Ministry of Education,East China University of Science and Technology,Shanghai 200237,China)

A compressible round jet is investigated by large eddy simulation.The jet Mach number is 0.75 and the Reynolds number is 50 000.Numerical results are validated by available experimental data.Two major stages for the shear layer instability development are observed and analyzed.After this,two-point space-time correlations of axial velocity fluctuationu′,radial velocity fluctuationv′ and pressure fluctuationp′ in the shear layer are investigated.Their distribution and evolution in time and space are presented and discussed.The migration characteristics of vortices are analyzed.The study provides a base for further investigation and controlling of jet noise.

subsonic; jet flow; large eddy simulation; shear layer; correlation analysis

1006-3080(2016)05-0715-07

10.14135/j.cnki.1006-3080.2016.05.020

2016-02-23

国家自然科学基金(51576067)

林 健 (1990-),男,江苏江阴人,硕士生,主要从事流体机械气动声学的研究。E-mail:lj307585298@163.com

赖焕新,E-mail:hlai@ecust.edu.cn

V211

A