含激波、旋涡和声波的复杂多尺度流动数值模拟研究

2018-06-29 11:06张树海王益民
空气动力学学报 2018年3期
关键词:旋涡激波插值

张树海, 李 虎, 王益民

(中国空气动力研究与发展中心 空气动力学国家重点实验室, 四川 绵阳 621000)

0 引 言

激波和旋涡是可压缩流动的基本结构。在许多流动中,比如超声速混合层、超声速射流和燃烧等,激波与旋涡共存,它们之间会发生相互作用,如激波与激波之间、激波与旋涡之间以及旋涡与旋涡之间的相互作用,这种作用,不仅会引起激波和旋涡的变形,还会产生噪声。因此,吸引了很多学者开展理论、实验和数值模拟研究。

1955年,Hollingsworth 和Richards[1]首先研究了激波与旋涡的相互作用,他们利用激波管产生平面激波,该激波扫过一个具有一定攻角的机翼后,产生一个旋涡,激波碰到激波管端部的固壁后反射回来,与旋涡发生干扰。该实验发现,激波与旋涡干扰后,伴随着激波与旋涡变形会产生声波。Dosanjh 和Week[2]对此问题进行了较细致的实验,测量了第一道声波的强度,成为后来实验研究和数值模拟研究的一个标准。Barbosa 和Skews[3]设计了一个分叉型激波管,激波管所产生的激波被分成两支,一支通过激波管内的拐角产生旋涡,另一支直接与旋涡发生干扰,研究发现,激波与旋涡干扰后,在旋涡中心附近一个很小的区域,出现压力脉冲区,其脉冲压力是周围压力的两倍以上。

为了解释声波的产生机理,早在20世纪60年代,Ribner[4]就发展了线性分析理论。他利用傅立叶变换将旋涡分解成平面正弦剪切波,这种单正弦波与激波相互作用后产生声波,经过合成后得到的激波与旋涡干扰产生的声波在垂直于激波的方向反对称分

布,其理论结果与Dosanjh和Week[2]的测试结果的反对称部分一致。后来,Week和Dosanjh[5]发展了Lightill[6]的声比拟理论,将周向压力分布展开成四极、两极和单极声源之和。单极的作用使得声波偏离了单纯反对称态,他们所得到的声波与其实验一致。

到20世纪90年代,伴随着高阶精度捕捉激波格式的成熟,采用数值模拟研究激波和旋涡干扰并捕捉声波已成为研究声波产生机理的重要手段。Ellzey等[7-9]采用四阶差分格式,通过求解非定常可压缩Euler方程,系统地模拟了激波与旋涡的干扰,得到的第一道声波强度与实验测得的声波强度基本一致。Inoue等[10-12]利用六阶紧致格式直接模拟了不同强度激波与旋涡相互作用及声波的产生过程,给出了弱激波与旋涡相互作用声波的产生及传播机理。他们通过精细的模拟发现,当激波碰到旋涡的瞬间,旋涡边缘发生局部压缩和膨胀,形成双极的第一道声波,随着激波与旋涡干扰的继续,出现第二组压缩和膨胀区,此压缩和膨胀区与前面出现的压缩和膨胀区域交替出现,形成第一道完整的声波,第一道声波演变成四极波。激波通过旋涡后,在旋涡两侧,形成两道反射激波,反射激波后伴随第二道和第三道声波。Grasso和Pirrozli[13-16]也系统地模拟了较大范围激波与旋涡的干扰问题,并根据激波变形的形状,将激波与旋涡的干扰分成三类:弱干扰,具有规则反射的强干扰,具有马赫反射的强干扰。

除了激波与旋涡之间的相互作用外,在包含激波和旋涡的复杂流场中,还存在旋涡与旋涡之间的相互作用,两个或多个旋涡之间的相互作用也是噪声产生的重要机制。比如,两个同向旋转的旋涡的合并过程,是剪切层噪声的重要部分。

这种包含激波、旋涡的复杂流场的非定常运动产生的噪声是激波噪声。目前,虽然气动声学理论已有相当的发展,但是关于激波噪声的产生机制认识还很不够。伴随高阶精度数值格式和大型并行计算机的快速发展,数值模拟已成为研究激波噪声产生机制的重要手段,但是激波噪声的计算仍然是计算气动声学所面临的相当严峻的挑战,具体表现在:

(1) 流场的强非定常特性。激波噪声是激波与其它流场结构相互作用产生的,其中存在激波变形和产生过程,与激波相互作用的流场结构也随之发生变化。比如,在激波与旋涡的相互作用过程中[17-20],入射激波会变形并产生反射激波,在旋涡中心附近还产生小激波串,同时旋涡被压扁,而且出现扭转运动,当激波强度达到一定程度之后,会引起旋涡的破碎。在激波与剪切层相互作用过程中,存在有激波的局部震荡[21]。因此,传统的定常流体动力学问题计算方法或是大尺度问题的非定常方法已不再适用。尽管已发展了多种高阶精度Runge-Kutta方法[22],但大多Runge-Kutta方法是显式格式,由于稳定性条件的限制,允许的时间步长非常小,对包含边界层的流动计算效率非常低,因此仍需要发展时间离散精度和效率更高的计算方法。

(2) 流场的多尺度特性。激波噪声是流场的一部分,但是在很多情况下,它的扰动尺度比流体动力学尺度小很多[23]。比如在马赫数为1.5的喷流中,距喷流直径40倍处测得的噪声强度是124 dB,声场的脉动速度与喷流的速度之比是1.5×10-4。特别地,当流场中存在强激波时,声波引起的物理量脉动比激波引起的变化尺度差距更大。这就要求不仅能精确计算宏观尺度的流体动力学量,能够捕捉激波,还能精确计算微弱尺度的声波,对计算方法和计算网格的要求远比动力学计算高。

(3) 模拟激波噪声所需的计算区域远大于模拟气动力所需的计算区域。研究激波噪声既关心近场流体动力学的变化,挖掘噪声产生的机制,又关心远场激波噪声的传播特性,计算区域边界与声源的距离非常远,这是与传统计算流体动力学的鲜明差异。

(4) 激波与声波的计算对格式需求存在矛盾的因素。捕捉激波需要格式具有一定的耗散,而声波是近乎无色散、无耗散的等熵运动,声波的模拟需要高分辨率和无耗散格式。因此,现在还不存在一种较为理想的高阶精度、高分辨率、低耗散的激波捕捉格式。目前,声波的计算普遍采用线性紧致格式,精度和分辨率非常高,耗散性小,适合于声波的计算。但是线性格式不能捕捉流场中的强激波,而捕捉激波普遍采用WENO格式[24]。尽管WENO格式的鲁棒性非常好,可以模拟很强的激波,但同所有的捕捉激波格式一样,在激波下游会产生微弱的非物理波动。虽然这种波动不影响流场的整体特性,但其强度常常会比流场中的噪声高。另外,虽然可以设计很高精度的WENO格式[25],但格式的耗散性比较大,对声波的分辨率不满意。尽管已有很多种简单的混合格式[26]、优化WENO格式[27]、非线性加权紧致格式[28-29]以及DRP格式[30],试图提高捕捉激波类格式的分辨率,但实际应用比WENO格式并没有明显的优势,因此目前还不存在一种既能捕捉激波又适用于声波计算的较理想数值方法。

针对这些问题,我们近年开展了较为系统性的研究,解决了计算激波噪声的一些方法上存在的问题,对包含有激波、旋涡和声波的一系列问题进行了直接数值模拟,揭示了激波噪声产生机制。本文对这些研究做简要回顾。

1 针对激波噪声计算的数值方法的研究

1.1 WENO格式的不收敛问题[31-32]

WENO格式是一种高阶精度捕捉激波的格式[24],被广泛用于工程实际和科学计算。但是像所有其它的激波捕捉格式一样,用WENO格式计算定常激波时,在激波下游会出现微弱的非物理波动,导致残差不能收敛到满意程度。虽然这种波动不会影响整体的气动特性,但是至少会造成两方面的影响,第一,残差是结束计算的重要判据,残差无法收敛到很小的值会使计算无法判断何时停止,只能靠经验结束进程;第二,虽然激波后的非物理波动是微弱的,不影响气动力等宏观流动特性,但是对于激波噪声,这种微弱的非物理波动则非常大,有时比激波噪声的幅度高几个数量级,因此计算气动噪声问题,应该消除这种微弱的非物理波动。经过系统的数学物理分析和大量的数值实验,我们提出了两种方法消除这种微弱的非物理波动[31-32]。

1) 针对五阶WENO格式,我们提出了一种新的光滑测试因子如下:

(1)

具体形式为:

(2)

采用这种光滑测试因子后,消除了激波下游的微弱非物理波动,对典型问题计算可以收敛到机器零。

2) 前面的方法只适合于五阶精度的格式,对更高阶精度格式不适用,而且会引起格式在高阶极值点附近降阶。经过大量的数值实验,我们发现WENO格式在特征投影过程中,采用Roe平均方法计算半点上的物理量确定特征值和特征矩阵,是影响WENO格式不收敛的重要原因。为此,我们提出采用迎风插值方法代替原来的Roe平均方法。即:

具体的迎风插值公式如下:

一阶插值:U(1)=UiU(2)=Ui+1

(3)

(4)

最优加权插值:

(5)

(6)

最优插值精度与WENO重构过程精度一致,如五阶精度的WENO格式采用五阶插值,七阶精度的WENO格式采用七阶插值,九阶精度的WENO格式采用九阶插值。图1和图2是采用两种方法对马赫数2的一维定常激波计算得到的密度分布和收敛历程,图1中ResA为残差的平均值。可以看到,改进的WENO格式可以收敛到机器零。值得指出的是,采用最优加权插值方法求解半点上的值,使WENO格式具有自相似特性,并一致高阶精度,详见文献[31-32]。

1.2 一种类谱分辨率的高阶精度紧致格式[34-35]

WENO格式是一种高阶精度捕捉激波格式,虽然格式的精度可以很高,但是对声波的分辨率并不理想,而且耗散性也比较大,采用WENO格式模拟湍流结构、激波噪声等复杂多尺度流动并不理想。为了提高格式的分辨率,降低耗散,并保持捕捉激波的鲁棒性,我们将WENO格式的设计思想应用到紧致格式中,设计了一种类谱分辨率的高阶精度紧致格式。这一格式是基于Lele的半点型线性紧致格式[33]设计的,Lele的半点型线性紧致格式形式是:

(7)

该格式左端是网格点上的函数导数值,右端是半点的函数值,它最大的优点是分辨率非常高,接近谱方法分辨率。但是右端半点的函数值是未知的,Lele[33]设计了一种线性紧致插值的方法求得半点的值。图3是Lele的半点型线性紧致格式的设计模板,可以看出,该格式有两个缺点:(1)设计模板包括了整点和半点的函数值,但Lele仅用了部分信息,格式没有达到该模板上的最优精度;(2)半点上的值需要插值方法求得,任何插值方法都存在传递误差,导致格式对高波数的分辨率大幅度下降。

为克服Lele格式的缺点,我们设计了一种类谱分辨率的高阶精度紧致格式。在网格节点上,我们采用如下格式计算其一阶导数:

(8)

采用同样的格式计算半点上的一阶导数值:

(9)

求得网格点和半点处的一阶导数之后,通过求解相应的控制方程,即可求得下一时刻的网格点和半点处函数值。

与Lele的格式相比,式(8)、式(9)给出的数值格式有两个优点:1) 格式的精度有一定的提高,Lele格式最高达到十阶精度,我们设计的格式可以达到十四阶;2) 网格点和半点上的函数值都是通过同一格式计算得到,消除了原格式紧致插值引起的传递误差,保证了计算过程的精度和分辨率。图4给出了六阶和八阶格式修正波数及其与 Lele 格式的比较,相比之下我们改进的格式分辨率与精确值相差无几,详见文献[34]。

线性格式不能捕捉流场中的强激波,为了使格式具有捕捉强激波的功能,我们将WENO重构思想应用到插值技术,发展了高阶WENO插值方法,并将该方法应用于式(8)、式(9)右端,初步发展了一种非线性类谱分辨率的紧致格式。具体的WENO插值方法是:在格式的设计模板S2r-1={xi-r+1,...,xi+r-1}上,采用插值函数:

(10)

可以得到半点上的通量值:

(11)

与WENO格式重构类似,将模板S2r-1分成三个子模板,在每个子模板上,存在r阶插值:

(12)

采用加权思路,可以得到半点上的加权插值:

(13)

在式(8)的右端,采用如下混合插值方法替代原来的精确值fi+1/2,即:

(14)

而在式(9)中,采用类似的混合插值方法。得到的非线性格式,能够保证格式的精度,具有较好的捕捉激波能力,同时格式的分辨率比现有的WENO明显提高,耗散性也大幅度降低,详见文献[35]。

2 包含激波、旋涡和声波的复杂多尺度流动数值模拟

2.1 两个同向旋转的Gaussian涡合并过程[20]

Gaussian涡是一种典型的旋涡,在一定的条件下,两个同向旋转的Gaussian涡会合并[36],对于合并条件,已有很多研究[37]。1995年,Mitchell等[36]采用直接数值模拟方法,研究了两个同向旋转的

(a)r=λ/2

(b)r=2λ

图7两个同向旋转的Gaussian涡合并过程中
产生的声波及其比较
Fig.7Far-fieldpressuretracesatr=λ/2andr=2λandthecomparisonbetweenourdirectnumericalsimulation(DNS)andthatbytheMöhring’sequation

值得指出的是,我们的结果与Mitchell的结果并不一致,其原因可能是由于初始参数的微弱不一致性,导致两个旋涡合并时间的差异。为了验证我们结果的正确性,我们进行了网格收敛性研究,初始条件的影响以及采用不同方法对此问题进行了模拟,结果一致,详见文献[20]。

2.2 两个Taylor旋涡的相互作用[20]

Taylor涡也是一种典型的旋涡,与Gaussian涡结构不同,Taylor涡具有双层结构。为了深入了解涡声的产生机理,我们系统研究了两个Taylor旋涡的相互作用,并根据旋涡特性,将两个Taylor旋涡的相互作用分成四类:(1) 两个反向旋转的Taylor涡的相互作用;(2) 两个同向旋转Taylor涡的相互作用;(3) 两个同向旋转Taylor涡的合并;(4) 两个尺度或强度相差较大的旋涡的相互作用。

2.2.1 两个反向旋转的旋涡的相互作用

Taylor涡具有双层结构,内层外层涡量符号相反。如果两个反向旋转的旋涡距离足够近,它们之间会发生耦合,结果导致内层外层分离,两个涡核逐渐靠拢,形成一个涡偶极子,两个Taylor涡的外层逐渐靠拢,形成一个反向运动的涡偶极子,如图8所示,是强度为Mv=0.5、初始距离为4倍旋涡半径的两个反向涡相互作用的涡量演化历程。详见文献[20]。

2.2.2 两个同向旋转Taylor旋涡的相互作用

图9是两个强度为Mv=0.5、初始旋涡间距是4倍旋涡半径的两个同向旋转的Taylor旋涡相互作用过程的涡量演化历程。这种相互作用存在两种特性,1) 两个同向旋转的旋涡将演化成两个不对称的旋涡区,每一个涡区存在三个涡核结构;2) 与两个反向旋转的旋涡相互作用具有本质区别的是,每一个涡区的两个涡核来自同一个初始旋涡,强的涡核来自初始旋涡的内层,而初始旋涡的外层逐渐从初始旋涡中分离形成围绕涡核的两个独立的涡核。

2.2.3 两个同向旋转Taylor旋涡的合并过程

与两个同向旋转的Gaussian旋涡相似,在一定条件下两个同向旋转的Taylor旋涡也会合并。如图10所示,是两个强度为Mv=0.5的同向旋转的Taylor旋涡的合并过程,两个旋涡的初始距离是旋涡半径的2倍。两个涡核合并形成主涡核,而在涡核的外部,形成两个旋转臂结构,两个涡的外层逐渐集中于旋转臂附近,形成一个三核涡结构。

2.2.4 两个强度或尺度相差较大的旋涡的相互作用

两个强度相差很大的旋涡的耦合过程是最令人感兴趣的。1) 两个强度或尺度相差较大的旋涡的相互作用会产生多核涡结构;2) 强度或尺度的较大差异,导致它们之间相互作用所产生的多核涡结构,具有非对称特性。这种相互作用包括有两种典型状态,一是强度相差较大的两个Taylor涡的相互作用,另一是尺度相差较大的Taylor涡的相互作用。

图11是两个强度相差较大的Taylor涡的相互作用过程,其中强旋涡的强度为Mvu=-0.8,顺时针旋转,弱旋涡的强度为Mvd=0.25,逆时针旋转,它们之间的初始距离为半径的4倍。可以看到,与两个等强度的反向旋涡相互作用类似,两个强度相差较大的Taylor涡相互作用会产生两个涡偶极子。两个旋涡的涡核逐渐靠拢,形成一个强的涡偶极子,两个旋涡的外层逐渐从原旋涡中分离,形成一个弱的涡偶极子。由于其强度相差较大,两个涡偶极子都呈现非对称特性,导致弱的涡核围绕强涡核旋转,较弱的涡偶极子围绕强涡偶极子旋转。

图12是两个尺度相差较大的反向旋涡相互作用的过程,其中两个旋涡的强度都是0.5,一个逆时针旋转,一个顺时针旋转。两个旋涡的尺寸相差较大,上面的旋涡半径为1, 而下面的旋涡半径为0.2, 两个旋涡的初始距离为2.2。可以看到,小旋涡的涡核穿过大旋涡的外层,进入到大旋涡的涡核外围,并围绕大旋涡的涡核不断旋转,最终形成一个椭圆形的涡核,而旋涡的外层形成两个涡核结构,围绕大的涡核旋转,形成一个三核涡结构,并持续做蛙跳运动。

2.2.5 两个Taylor涡相互作用过程中声波产生机理分析

两个Taylor涡的相互作用包括丰富的动力学特性,这种过程会产生噪声。图13和图14 是两个典型状态两个Taylor涡相互作用的声压Δp=(p-p0)/p0等值线和其径向分布。旋涡的强度均是0.5,两个旋涡初始距离都是半径的4倍。其中图13是两个反向Taylor涡相互作用所产生的噪声,图14是两个同向Taylor涡相互作用所产生的噪声。可以看到,两种旋涡相互作用产生的噪声具有很大区别,两个反向旋涡相互作用仅产生几个声脉冲,而两个同向旋转的Taylor涡相互作用,与两个同向旋转的Gaussian涡相互作用很类似,会产生持续的强噪声。

以两个尺度相差较大的旋涡相互作用为例,分析噪声产生机理。我们计算了Lamb量·(ρω×u),它代表着噪声产生的源。如图15所示,是Lamb量的等值线演化过程。可以看到,旋涡相互作用具有三个明显的运动:1) 初始圆形旋涡被扭曲,弱旋涡被撕扯成两个涡带,他们与强旋涡继续相互作用;2) 尽管强旋涡的涡心位置几乎不变,但整个涡核被压成椭圆形;3) 两个涡带与强旋涡的涡核之间相互作用产生一种蛙跳式运动。图16是监测点(x,y)=(100,0)和(0,100)的声压随时间变化曲线,压力峰值P1、P2、P3、P4、P5、P6和P7位于时间轴103、121、139、163、190、223和262。考虑到声从涡心附近产生声波到达监测点所需要的时间,它们对应于长轴在y轴的方向,而声波谷M1、M2、M3、M4、M5、M6和M7位于时间轴114、130、151、176、206、241和285,它们对应于长轴位于x的方向。这意味着,噪声主要由旋涡的蛙跳运动产生。

2.3 激波与旋涡的相互作用[17-19]

激波与旋涡相互作用是激波噪声的一个简化模型,通过研究激波与旋涡相互作用,可以深入了解激波噪声产生机制。通过对激波与旋涡相互作用系统的模拟,我们发现激波与强旋涡相互作用存在多级干扰。第一级干扰是入射激波和初始旋涡的相互作用,第二级干扰是反射激波和变形旋涡的相互作用,第三级干扰是涡心附近产生的小激波串和变形旋涡的干扰,如图17所示。与此同时,旋涡的变形也存在多级特征。在第一级干扰中,初始圆形旋涡被压成椭圆形;在第二级干扰中,椭圆形旋涡被压回成圆形旋涡;在第三级干扰中,圆形旋涡又被压缩成椭圆形。此外,这种多级干扰会引起旋涡在水平线附近做振荡运动。根据激波与强旋涡相互作用过程中的多级干扰过程,当时我们预测,这种多级干扰过程会产生更为复杂的声波。这一预测被Chatterjee和Vijayaraj[39]的研究所证实,如图18所示。此外,我们还系统研究了激波与涡列相互作用过程流场结构和声波产生机制[17-19]。

2.4 激波与剪切层相互作用[40]

激波与剪切层相互作用是一个典型的问题,也可以看做是研究超声速喷流激波噪声的一个简化模型。为了了解超声速喷流激波噪声的产生机制,我们设计了一个激波与剪切层相互作用的模型。如图19所示,一道斜激波打到一个剪切层上,下边界是一个反射固壁,透射激波打到固壁上后,形成反射激波,反射激波再次与剪切层相互作用发生反射,形成类似于超声速喷流的激波串结构。通过较为系统的数值模拟,我们揭示了两种噪声产生机制,一是激波与旋涡相互作用,发生在图20所示的区域1, 另一是激波泄露机制,发生在图20所示的区域2。

图21是第一个区域入射激波与剪切层中旋涡相互作用过程的一个演化周期中几个典型时刻的胀量图,其中红色虚线标记的是剪切层下层流体中被追踪的两道声波。我们看到由于流体以超声速向下游传播,声波也随着流体向下游方向对流,并且声波的辐射半径在逐渐增大。当t=91.69时,产生了第二道声波,可以明确地看到声波的产生源在涡与小激波相互作用的位置,这和张树海等[17-19]的激波与旋涡相互作用的研究结果相同。剪切层的涡列穿过激波,每个涡与激波发生干扰时,根据这种激波噪声产生和传播的机制,在噪声源处向外辐射出一道道的声波,同时我们可以看到当声波传播到下壁面时声波反射回内场,从而声波再次与激波和剪切层相互干扰。

图22是第二个区域反射激波与剪切层相互作用过程中一个演化周期中几个典型时刻的胀量图,其中红色虚线标记的是剪切层下层流体中的一道声波。我们看到由于流体以超声速向下游传播,声波也随着流体向下游方向对流,并且声波的辐射半径在逐渐增大。当t=91.69时,噪声源处产生了第二道声波,可以明确地看到声波的产生源在旋涡之间的鞍点处,即辫子区的位置,这和Manning等[41-43]的激波泄漏机制完全相同。在剪切层的涡列穿过激波时,涡对与激波发生干扰,在涡对之间的鞍点处激波透过剪切层,同时在鞍点处激波被辫子区的旋臂所阻碍,激波与辫子区的相互干扰使得激波强度变弱,部分激波穿过剪切层,还有一部分激波没有穿过剪切层,而是在鞍点处以声波的形式发生泄漏。这说明了强激波与剪切层相互作用中也存在激波泄漏机制,在鞍点处激波泄漏并向外辐射声波。根据这种激波噪声产生和传播的机制,当剪切层穿过激波时,在噪声源处向外辐射出一道道的声波,同时我们可以看到当声波传播到下壁面时声波反射回内场,从而声波再次与激波和剪切层相互干扰。

2.5 轴对称超声速喷流[44]

超声速喷流噪声是一种典型的激波噪声,与很多工程密切相关。采用五阶WENO格式直接求解轴对称Navier-Stokes方程,数值模拟了射流马赫数(完全膨胀)为Mj=1.19(啸声呈现轴对称模态,三维效应不明显)的轴对称欠膨胀射流,其对应的射流声速马赫数为Ma=1.05,雷诺数为Re=6.216×105,喷管直径为D=25.4 mm,喷管厚度为0.4D,环境温度为288.15 K。

图23给出了喷管唇口壁面[0.0,0.642D]和喷管出口平面[0.0,2.0D]上压力监测点处噪声信号的频谱信息,包括频率和声压级。结果显示,在大于5000 Hz的频率区间内,有4个尖锐的声压级突起:前两个声压级突起的频率为6567 Hz(128 dB)和8637 Hz(123 dB),分别是啸声A1模态和啸声A2模态;另外两个声压级突起的频率为12087 Hz(125 dB)和14310 Hz(123 dB),分别是啸声B模态的谐频和啸声A0模态。

(a) 喷管唇口壁面

(b) 喷管出口平面

图23监测点处压力信号的谱分析
Fig.23Soundpressurelevel

综合流场结构和声场信息可以得到的啸声模态的图像及其产生位置,其中流场结构用数值纹影(密度梯度的模)表示,声场由胀量场表示。根据啸声频率可以得到啸声的波长,利用波长可以在胀量场中清晰地分辨出啸声A1模态和A2模态,相应的波长分别为λ=2.04D和λ=1.51D。从图24中可以清晰地分辨出向上游传播的啸声,并且啸声产生位置都是在第三个激波栅格和第五个激波栅格之间的区域。

2.6 声波穿过激波的畸变过程

声波产生在近流场区域,流体介质的非定常运动产生的噪声经过复杂流场结构传播达到远场。当噪声穿过复杂流场结构过程中,会与流场相互作用,发生畸变。其中,声波穿过激波是一个典型问题,我们采用直接数值模拟方法,研究了声波穿过激波的畸变特性。

初始条件为:

(15)

(16)

(17)

计算区域为x∈[-10,10],马赫数M=2.0,声波幅值ε=1.0d-5,频率ω=0.6π,比热比取γ=1.4,计算结果如图25所示,由图可以看出,声波穿过马赫数M=2.0的激波,振幅放大了约3.64倍多。

2.7 声波与旋涡相互作用

当声波穿过旋涡时,声波的幅值和频率都会发生变化。采用我们发展的线性紧致格式,通过求解Navier-Stokes方程,我们研究了平面声波穿过等熵涡的散射特性。

图26是计算模型,平面声波与一个二维Taylor涡相互作用。计算区域的左端入射声波为:

(18)

二维Taylor涡置于计算区域的中心。图27是入射声波频率ω=0.6π与强度为Mv=0.25的Taylor涡相互作用的瞬态脉动压力场云图。可以看到,旋涡对声波的影响非常明显。在旋涡下游,形成两条强噪声干扰条带,在强条带外侧,有二次条带,由于旋涡以逆时针方向旋转,干扰条带关于y=0不对称。图28是散射声压的均方根(RMS)等值线图,图29是离涡心r=10的圆周上散射声压的指向性分布。

3 结束语

近年来,针对激波噪声计算的需求,我们开展了数值方法研究,通过对典型问题的直接数值模拟,揭示了激波噪声产生机理。

(1) 在数值方法方面,针对高阶精度捕捉激波的WENO格式不收敛问题,提出了一种新的光滑测试因子,发展了一种迎风插值技术计算半点上物理量,用以确定WENO重构过程中特征投影所需物理量。所改进的WENO格式,消除了激波下游的微弱非物理波动,对典型问题计算可以收敛到机器零。发展了一种类谱方法的紧致格式,这一格式包括线性和非线性两部分,其中线性格式具有高阶精度、高分辨率和低耗散特性,是模拟低速气动声学问题的理想算法。非线性格式捕捉激波能力与WENO格式相当,分辨率明显提高,耗散性大幅度降低。

(2) 采用直接数值模拟方法,系统研究了两个旋涡相互作用、激波与旋涡/涡列相互作用、激波与剪切层相互作用以及超声速喷流激波噪声。发现激波与强旋涡相互作用具有多级干扰特性,提出每级干扰都有自己的声波产生。

当然,对激波噪声的认识还相当有限,有许多问题有待进一步解决。

参 考 文 献:

[1]Hollingsworth M A, Richards E J. On the sound generated by the interaction of a vortex and a shock wave[R]. British Aeronautical Research Council, Fluid Motion Sub-Committee, 18257, FM2371, 1956.

[2]Dosanjh D S, Weeks T M. Interaction of a starting vortex as well as a vortex street with a travelling shock wave[J]. AIAA J, 1965, (3): 216-223.

[3]Barbosa F J, Skews B W. Shock wave interaction with a spiral vortex[J]. Phys Fluids, 2001, 13: 3049-3060.

[4]Ribner H S. The sound generated by interaction of a single vortex with a shock wave[R]. UTIA Report No. 61, 1959.

[5]Weeks T M, Dosanjh D S. Sound generation by shock vortex interaction[J]. AIAA J, 1967, (5): 660-669.

[6]Lighthill M J. On sound generated aerodynamically I: general theory[J]. Proceedings of the Royal Society of London, 1952, 211: 564-587.

[7]Ellzey J L, Henneke M R, Picone J M, et al. The interaction of a shock with a vortex: shock distortion and the production of acoustic waves[J]. Phys Fluids, 1995, 7: 172-184.

[8]Ellzey J L, Henneke M R. The shock-vortex interaction: the origins of the acoustic wave[J]. Fluid Dyn Res, 1997, 21: 171-184.

[9]Ellzey J L, Henneke M R. The acoustic wave from a shock-vortex interaction: comparison between theory and computation[J]. Fluid Dyn Res, 2000, 27: 53-64.

[10]Inoue O, Hattori Y. Sound generation by shock-vortex interaction[J]. J Fluid Mech, 1999, 380: 81-116.

[11]Inoue O. Propagation of sound generated by weak-shock interaction[J]. Phys Fluids, 2000, 12: 1258-1261.

[12]Inoue O, Takahashi T, Hatakeyama N. Separation of reflected shock waves due to secondary interaction with vortices: Another mechanism of sound generation[J]. Phys Fluids, 2002, 14: 3733-3744.

[13]Grasso F, Pirozzoli S. Shock-wave thermal inhomogeneity week shock-vortex interaction in a mixing zone[J]. AIAA J, 1999, 33: 1797-1802.

[14]Grasso F, Pirozzoli S. Shock-wave-vortex interactions: shock and vortex deformations, and sound production[J]. Theroet Comput Fluid Dyn, 2000, 13: 421-456.

[15]Grasso F, Pirozzoli S. Simulations and analysis of the coupling process of compressible vortex pairs: Free evolution and shock induced coupling[J]. Phys Fluids, 2001, 13: 1343-1366.

[16]Pirozzoli S, Grasso F, D’Andrea A. Interaction of a shock wave with two counter-rotating vortices: Shock dynamics and sound production[J]. Phys Fluids, 2001, 13: 3460-3481.

[17]Zhang S, Zhang Y T, Shu C W. Multistage interaction of a shock wave and a strong vortex[J]. Phys Fluids, 2005, 17: 116101.

[18]Zhang S, Zhang Y T, Shu C W. Interaction of an oblique shock wave with a pair of parallelvortices: Shock dynamics and mechanism of sound generation[J]. Phys Fluids, 2006, 18(12): 126101.

[19]Zhang S, Jiang S, Yong-Tao Zhang, et al. The mechanism of sound generation in the interaction between a shock wave and two counter-rotating vortices[J]. Phys Fluids, 2009, 21: 076101.

[20]Zhang S, Li H, Liu X, et al. Classification and sound generation of two-dimensional interaction of two Taylor vortices[J]. Physics of Fluids, 2013, 25: 056103.

[21]Lui C V. C. M. A numerical investigation of shock-associated noise[D]. Stanford University, 2003.

[22]Gottlib S, Shu C W. Total variation diminishing Runge-Kuttaschemes[J]. Mathematics of Computation, 1998, 67: 73-85.

[23]Tim Colonius, SanjivaK Lele. Computational aeroacoustics: progress on nonlinear problems of sound generation[J]. Progress in Aerospace Sciences, 2004, 40: 345-416.

[24]Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126: 202-228.

[25]Balsara D S, Shu C W. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order accuracy[J]. Journal of Computational Physics, 2000, 160: 405-452.

[26]Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction[J]. J Comput Phys, 2002, 178: 81-117.

[27]Wang Z J, Chen R F. Optimized weighted essentially non-oscillatory schemes for linear waves with discontinuity[J]. J Comput Phys, 2001, 174: 381-404.

[28]Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes[J]. J Comput Phys, 2000, 165: 22-44.

[29]Zhang S, Jiang S, Shu C W. Development of nonlinear weighted compact schemes with increasingly high order accuracy[J]. J Comput Phys, 2008, 227: 7294-7321.

[30]Tam C K W, Webb J C. Dispersion-Relation-Preserving finite difference schemes for computational acoustics[J]. J Comput Phys, 1993, 107: 262-281.

[31]Zhang S, Shu C W. A new smoothness indicator for the WENO schemes and its effect on the convergence to steady state solutions[J]. Journal of Scientific Computing, 2007, 31: 273-305.

[32]Zhang S, Jiang S, Shu C W. Improvement of convergence to steady state solutions of Euler equations with the WENO schemes[J]. Journal of Scientific Computing, 2011, 47: 216-238.

[33]Lele S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992, 103: 16-42.

[34]Liu X, Zhang S, Zhang H, et al. A new class of central compact schemes with spectral-like resolution I: Linear schemes[J]. Journal of Computational Phys, 2013, 248: 235-256.

[35]Liu X, Zhang S, Zhang H, et al. A new class of central compact schemes with spectral-like resolution II: Hybrid weighted nonlinear schemes[J]. Journal of Computational Phys, 2015, 284: 133-154.

[36]Mitchell B E, Lele S K, Moin P. Direct computation of the sound from a compressible co-rotating vortex pair[J]. Journal of Fluid Mechanics, 1995, 285: 181-202.

[37]Cerretelli and Williamson C H K. The physical mechanism for vortex merging[J]. J Fluid Mech, 2003, 475: 41-77.

[38]Mohring W. On vortex sound at low Mach number[J]. J Fluid Mech, 1978, 85: 685-691.

[39]Chatterjee A, Vijayaraj S. Multiple sound generation in interaction of shock wave with strong vortex[J]. AIAA Journal, 2008, 46: 2558-2567.

[40]Liu X, Zhang S. Direct numerical simulation of the interaction of 2D shock wave and shear layer[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45: 61-75. (in Chinese)刘旭亮, 张树海. 二维激波与剪切层相互作用的直接数值模拟研究[J]. 力学学报, 2013, 45: 61-75.

[41]Manning T A. A numerical investigation of sound generation in supersonic jet screech[D]. Stanford University, 1999.

[42]Suzuki T, Lele S K. Shock leakage through an unsteady vortex-laden mixing layer: application to jet screech[J]. Journal of Fluid Mechanics, 2003, 490: 139-167[43]Lui C C M. A numerical investigation of shock-associated noise[D]. Stanford University, 2003.

[44]Li H. Numerical study of shock-associated noise in axisymmetric supersonic jet[D]. China Aerodynamics Research and Development Center, 2016. (in Chinese)李虎. 轴对称超声速射流激波噪声数值模拟研究[D]. 中国空气动力研究与发展中心, 2016.

猜你喜欢
旋涡激波插值
滑动式Lagrange与Chebyshev插值方法对BDS精密星历内插及其精度分析
面向三维激波问题的装配方法
二元Barycentric-Newton混合有理插值
一种基于聚类分析的二维激波模式识别算法
基于HIFiRE-2超燃发动机内流道的激波边界层干扰分析
大班科学活动:神秘的旋涡
山间湖
基于pade逼近的重心有理混合插值新方法
斜激波入射V形钝前缘溢流口激波干扰研究
为领导干部荐书