激波/湍流边界层干扰压力脉动特性数值研究1)

2021-11-09童福林段俊亦周桂宇李新亮

力学学报 2021年7期

关键词：边界层外层激波

童福林段俊亦周桂宇李新亮

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

† (中国科学院力学研究所高温气体动力学国家重点实验室,北京 100190)

** (中国空气动力研究与发展中心计算空气动力研究所,四川绵阳 621000)

†† (中国科学院大学工程科学学院,北京 100049)

引言

激波与湍流边界层的相互作用问题广泛存在于各类高速飞行器中.已有研究表明[1-2],在强激波作用下,干扰区内压力脉动峰值急剧升高,分离激波的非定常运动导致物面脉动压力低频能量的急剧增强,这将使得飞行器结构出现振动疲劳问题,进而严重影响飞行安全.因此,深入研究激波/湍流边界层干扰区内压力脉动特性对认识和理解激波与湍流边界层的相互作用机制非常重要.

自20 世纪70 年代以来,国内外学者一直致力于对该问题的研究.早期研究主要以风洞试验为主,受限于实验测量技术,风洞试验研究多局限于物面压力脉动信号,并以此获得干扰区内典型流动结构的特征信息,如分离激波非定常运动特性、湍流结构、速度场等.Settle 等[3]研究了激波强度对压缩拐角物面压力脉动强度的影响规律.Ringuette 等[4]研究表明,低雷诺数下,干扰区内物面压力脉动强度降低而间歇性增强.在分离激波低频振荡方面,Erengil等[5]采用条件取样和变窗口系综平均技术分析了物面压力分布与分离激波非定常运动的相关性.结果表明,分离泡内的低频压力脉动与激波位置密切相关:随着激波往下游移动,分离泡内压力升高,反之则降低.Dolling 等[6]比较了不同构型诱导的激波/湍流边界层干扰区物面压力脉动特性,研究发现,分离激波的低频振荡机制主要来源于下游流场的局部或全局脉动,其过零频率与上游湍流边界层的特征频率无关.此外,Selig 等[7]基于物面压力与边界层质量通量的时空关联分析,给出了干扰区内湍流结构的演化规律.Piponniau 等[8]采用本征正交分解和线性随机估计方法,通过实验测得的物面压力脉动信息对速度场进行了重构.

由于风洞试验能够得到的压力信息极为有限,为此国内外学者针对干扰区内压力脉动特性也开展了大量的高精度数值模拟研究.特别是近些年,直接数值模拟方法(direct numerical simulation,DNS)已经在激波/湍流边界层干扰问题方面取得了长足进展,其中Pirozzoli 等[9-11]的工作最具代表性.Pirozzoli等[9]分析比较了干扰区内物面压力脉动与剪切层压力脉动的差异,并基于声共振机制成功解释了分离激波的非定常运动现象.Bernardini 等[10]着重探讨了干扰区内物面压力脉动特性的演化规律,如脉动强度、时空关联、对流速度和频谱等,发现物面压力脉动对流速度明显降低而其展向尺度则显著增强.随后,Volpiani 等[11]进一步考察了壁面温度对物面压力脉动特性的影响规律.与此同时,国内童福林等针对物面压力脉动特性也开展了较为详细的直接数值模拟研究,归纳了激波强度[12]、壁面温度[13]、膨胀效应[14-15]等因素的影响规律.

尽管在激波/湍流边界层干扰区压力脉动特性研究方面,国内外已开展了大量的风洞试验和数值研究工作,但是这些研究大多关注物面压力脉动特性及其与分离激波非定常运动的相互关系方面,而关于边界层内压力脉动在干扰区内演化机制的认识仍然是很有限的.实际上,Duan 等[16]的研究结果表明,边界层内压力脉动结构特征与物面压力脉动完全不同,其特征频率沿法向将降低而特征尺度将急剧增大.然而,目前学术界对边界层内外层大尺度压力脉动结构与分离激波的相互作用机制仍然缺乏充分的认识,亟待开展相关研究工作.

本文采用直接数值模拟方法研究M∞=2.25,Reθ=3567 来流下激波角β=33.2°的入射激波/平板湍流边界层干扰区内压力脉动特性.为了便于比较和验证结果,计算参数的选取与Dupont 等[17]的实验和Fang 等[18]的DNS 相近.计算结果与前人的DNS 结果及风洞实验数据进行了仔细验证,着重探讨干扰区外层压力脉动与物面压力脉动演化机制的差异,特别是激波干扰对脉动强度、功率谱、特征尺度及时空关联特性等的影响规律.

1 DNS 计算

控制方程为曲线坐标系(τ,ξ,η,ζ)下的三维可压缩无量纲Navier-Stokes 方程组

其中J−1为直角坐标系(x,y,z)变换为曲线坐标系(ξ,η,ζ) 的Jacobian 矩阵,U为守恒变量,F,G和H分别为(x,y,z)方向上对应的无黏通量,Fv,Gv和Hv分别为(x,y,z)方向上对应的黏性通量,具体表达式参见文献[19].本文采用的解算器是高精度有限差分软件OPENCFD-SC,该求解器目前已在压缩拐角[19-20]、膨胀角入射激波/湍流边界层干扰[14-15,21]等复杂流动问题中得到广泛应用和验证确认,可以保证计算结果的可靠性.

对流项计算时,采用Martin 等[22]优化构造的WENO_SYMBO_LMT 格式以及Steger-Warming流通量分裂方法,这样能够在极大抑制激波间断区的强数值振荡同时也保证了对边界层内不同尺度湍流脉动结构的高分辨率捕捉.黏性项采用八阶中心差分格式进行离散,时间方向上推进采用三阶Runge-Kutta 方法.

图1 为计算模型示意图,来流方向为图中从左往右,来流参数如下:马赫数M∞=2.25,基于单位长度来流雷诺数Re∞=2.5×104mm−1,静温T∞=169.44 K.模型长度Lx×Ly×Lz=137.6 mm×12.7 mm×4.4 mm,这里Lx,Ly和Lz分别为流向、法向和展向长度,坐标系原点取为计算域入口处.采用与文献[21]类似的方法生成入射激波,在计算域上边界xin=82.2 mm左右两侧分别设置为波前和波后参数,即:在xin左侧流动参数取为自由来流参数,而在xin右侧则按照激波关系式Rankine-Hugoniot 给出.这里激波角β取为33.2°,图中入射激波在壁面上的名义入射点为xs=101.6 mm.

图1 计算模型Fig.1 Computational model

计算网格点数为Nx×Ny×Nz=3700×300×250 (流向×法向×展向).图2 为计算网格示意图,为了便于展示,流向间隔了10 个点,法向间隔了5 个点.网格生成采用代数解析方法,这里流向网格点采用三段分布:上游转捩区0 mm

图2 计算网格示意图Fig.2 Sketch of the computational grid

边界条件设置如下:计算域入口处取为相同来流下的层流解(laminar inflow),具体分布如图3 所示.壁面取为无滑移和等温壁,壁温Tw=321.9 K,上边界采用简单无反射边界条件.出口处采用超声速出口无反射条件及缓冲区,以便消除扰动波反射带来的影响.计算域展向两侧采用周期性边界条件.

图3 入口处层流剖面Fig.3 Laminar profile at inlet

计算时,通过在上游壁面添加多频吹吸扰动从而转捩生成充分发展湍流边界层,进而与下游入射激波产生相互作用.这里,法向扰动速度Vbs的具体表达式如下[18]

DNS 计算时取lmax=10,mmax=5,扰动幅值A=0.2,扰动频率ϖ=0.628U∞/δ,这里δ表示参考点处边界层厚度,下文类似.扰动带起始点分别为xa=7.62 mm 和xb=20.32 mm.式中相位差φl和φm取为0 到1 的随机数.

表1 给出了参考点xref处湍流边界层参数与Fang 等[18]DNS 结果比较情况,表中Reδ,Reδ*,Reθ分别为基于边界层厚度、位移厚度和动量厚度的雷诺数,Cf为摩阻系数.两者较为接近,这也表明本文计算结果的准确性.为了进一步验证本文DNS 结果的可靠性,还与以往数值模拟结果和风洞实验数据进行了比较,包括湍流边界层的速度剖面、脉动强度,物面压力频谱以及平均压力分布.

表1 参考点xref 处湍流边界层参数Table 1 Turbulent boundary layer parameters at xref

首先,图4～图6 分别给出了上游湍流边界层在参考点xref处的统计特性.图4 为平均流向速度剖面,定义为

图4 平均流向速度剖面Fig.4 Mean streamwise velocity profile at xref

式中ρ和u分别为密度和流向速度,上标 − 表示平均量.需要特别说明的是,本文中平均指的是时间和展向平均.速度分布在黏性底层y+<10 符合线性律,在30

图5 湍流脉动强度Fig.5 Turbulence intensities at xref

图6 物面压力功率谱Fig.6 Power spectral density for wall pressure at xref

此外,图7 还给出了干扰区内物面压力分布情况,其中物面压力P*和流向坐标x*分别定义为

图7 物面压力分布Fig.7 Distribution of wall pressure

式中,P1是反射激波的波后压力值,xr*为转换后的反射激波脚对应的平均流向位置,这里采用与Fang等[18]类似的方法,即(P1+P∞)/2 对应的流向位置.图中红色虚线代表压力分布的无黏解.可见,计算得到的压力分布与Fang 等[18]的DNS 结果及Dupont 等[17]的实验结果相符较好.在黏性干扰作用下,干扰区内压力梯度变化较为缓和且在下游逐步趋近于无黏解.

2 脉动强度

图8 给出了上游湍流边界层(turbulent boundary layer,TBL)压力脉动强度沿法向分布情况.可见,峰值位置出现在近壁区,约为y/δ=0.04,其量值沿物面法向急剧降低,在边界层外缘,约为近壁峰值的28%.采用当地物面剪切应力τw进行无量纲化后,计算得到的脉动分布与Bernardini 等[28]和Duan 等[16]的DNS 结果符合较好.

图8 上游湍流边界层压力脉动强度分布Fig.8 Distribution of pressure fluctuation intensity for the incoming TBL

从图9 可以看出,与上游湍流边界层相比,干扰区内压力脉动强度呈现截然不同的分布规律.从量值上来看,脉动峰值呈现先急剧升高,这主要是激波的增强作用,随后在下游再附区伴随着扰动边界层的恢复则呈现缓慢衰减的流向演化趋势.

图9 干扰区压力脉动强度云图Fig.9 Contour of pressure fluctuation intensity in the interaction region

为了定量描述这一现象,图10 给出了物面压力脉动强度沿流向的分布情况.计算得到脉动强度与Fang 等[20]的DNS 结果符合良好,两者在整体往下偏移0.02 的情况下与Dupont 等[17]的实验数据也较为吻合.这主要是由于压力传感器截止频率太小导致实验数据明显偏低[17].图中符号S和R分别表示分离点和再附点位置,通过平均摩阻的过零点确定.可见,分离区内物面压力脉动急剧增强,其峰值约为上游的4.3 倍,随后在下游再附区脉动峰值缓慢降低,在x*−xr*=12 时,仍约为上游的2.2 倍.另外,从峰值位置来看,其在分离区呈现逐渐远离壁面的态势,这与激波诱导生成分离区剪切层相关,如图9中平均音速线所示(粉色点划线).强脉动主要出现在干扰区内压缩波系、剪切层附近等强剪切区域.此后,在下游再附边界层,其峰值位置仍集中在外层区域.

图10 物面压力脉动强度分布Fig.10 Distribution of wall pressure fluctuation intensity

3 功率谱

沿外层y/δ=0.8 和物面分别选取5 个流向站位进行对比分析,其中站位1 位于上游湍流边界层内,站位2 位于分离区内,站位3～ 5 位于下游再附边界层内,各站位具体位置见图9.如无特别说明,下文中外层均为y/δ=0.8 处.

图11 比较了上游湍流边界层站位1 压力脉动的预乘功率谱ωφ(ω),其中ω和φ(ω)分别为角频率和功率谱密度.功率谱密度的计算采用Pwelch 方法和Hamming 窗函数,同时对分段数据按50%重叠处理.为了便于定量比较,图中还采用当地脉动强度的平方进行了归一化处理.显而易见,外层压力脉动的时间尺度明显要大于物面处,前者峰值频率约为ωδ/U∞=3.0,后者约为ωδ/U∞=12.3,这与Duan 等[13]的研究结论是一致的.这是由于外层压力脉动主要表征为大尺度脉动结构.

图11 站位1 压力脉动预乘功率谱Fig.11 Pre-multiplied power spectral density for pressure fluctuations at station 1

图12 分别给出了干扰区外层和物面压力脉动预乘功率谱的变化情况.从两者的比较来看,以下几个方面值得特别关注.首先,相较于上游站位1,分离区站位2 处两者都出现了峰值频率向更低频区偏移的现象,这与激波的低频振荡运动密切相关.

图12 干扰区内不同站位压力脉动预乘功率谱Fig.12 Pre-multiplied power spectral density for pressure fluctuations at various streamwise locations

如图12(a)和图12(b)所示,外层和物面压力峰值频率分别约为ωδ/U∞=0.7 和ωδ/U∞=4.4.以往研究表明[29],在分离激波低频振荡的作用下,分离区内物面压力脉动的时间尺度通常比上游湍流边界层低1～ 2 个数量级.本文计算结果进一步验证了该结论,同时结果还表明,外层压力脉动也同样符合这一规律.在功率谱的定性分布方面,显然,外层压力脉动的变化较物面压力更为剧烈,这很可能是由于外层站位非常靠近分离激波,受其非定常运动的影响更为直接,而物面压力处于分离泡内(见图9),影响则要弱得多.

从再附区站位3～ 5 的结果来看,外层压力脉动呈现快速恢复的态势,且峰值频率向高频区偏移,逐步趋近于上游站位1,这与物面压力脉动的演化过程则完全不同.从图12(a)中可以清楚看到,再附边界层内物面压力脉动的低频能量占比仍相对较高,峰值频率依然维持在较低频区,并没有出现往高频区恢复的现象.造成两者该差异的原因可以认为主要有两方面:一方面,上游激波诱导分离形成了剪切层,分离泡上方剪切层内不稳定涡结构沿流向急剧发展,并进入到下游再附边界层内,这些不稳定结构极大地抑制了物面压力脉动结构的恢复;另一方面,外层压力脉动法向位置相对较高,尽管在站位2 处激波干扰对其影响显著,但是过激波后外层压力脉动受近壁区流动结构的影响则要小得多,主要表征为外层大尺度脉动结构的破碎过程,更为详细的定量分析可见下一节的两点相关统计结果.

4 两点相关

通过两点相关性考察干扰区外层压力脉动拟序结构长度尺度的变化规律,且定量分析比较其与物面压力脉动结构尺度的差异.两点相关系数Cpp(Δx,Δz)定义为

式中p′为雷诺平均后得到的压力脉动,Δx和Δz分别代表流向间距和展向间距,x0为特征点流向位置,yref分别取物面和外层,上横线代表时间和展向平均.

图13 分别给出了上游湍流边界层外层和物面脉动压力两点相关系数的分布情况.以往的零压力梯度平板结果表明[28],小间距时脉动压力相关性曲线可近似为圆形,这是小尺度压力脉动结构的各向同性所决定的;大间距时,大尺度脉动结构的各向异性导致相关性曲线在展向出现了急剧拉伸,近似为椭圆形分布,表征了展向拉长的大尺度结构.显然,图13(a)中的结果也完全符合该规律.然而,对于外层压力脉动,如图13(b)所示,整体来看,外层脉动结构的流向尺度和展向尺度都明显大于物面压力脉动.这里以相关性系数0.3 为参考确定其特征长度,外层结构的流向和展向尺度分别约为0.7δ和1.0δ,而物面压力脉动则分别约为0.14δ和0.2δ.此外,尽管外层相关性等值线小间距时仍趋近于圆形,但随着间距增大,等值线在流向和展向都被急剧拉伸,展向明显强于流向,这表明外层压力脉动拟序结构仍以展向大尺度结构为主,但其特征尺度要大得多.

图13 站位1 压力脉动两点相关系数Fig.13 Two-point correlation coefficient of pressure fluctuations at station 1

图14(a)和图14(b)分别给出了站位2 物面和外层压力脉动的两点相关等值线,其中等值线从外到内依次为0.3,0.5,0.7 和0.9.为了便于比较说明,图中也给出了站位1 的结果.从分布规律上来看,等值线形态基本类似于上游湍流边界层,仍然以展向拉伸的椭圆形分布为主,但整体呈现向外扩张的趋势,这表明激波干扰作用增大了外层和物面压力脉动结构尺度,特别是在展向.需要特别关注的是,从图14(b)中可以清楚看到,同样以等值线0.3 为参考,站位1 处的展向尺度约为1.0δ,而站位2 处增大到约2.0δ,而流向尺度两者基本维持在0.7δ.这很可能是由于站位2 外层位于入射激波附近(见图9),激波的展向结构增强了压力脉动的展向相关性.更为细致的定量分析还有待深入研究.

图14 站位2 两点相关系数(虚线:站位1)Fig.14 Two-point correlation coefficient of pressure fluctuations at station 2(dashed lines:station 1)

如图15 所示,再附区站位3～ 5 外层和物面压力脉动的两点相关等值线与上游湍流边界层也基本类似.小间距时,曲线接近于圆形,这表明小尺度结构并没有发生实质改变,仍以各向同性为主.从大尺度结构的演化历程来看,两者差异显著.研究发现,大间距时,图15(a)～图15(c)中等值线以向外扩张为主,这表明再附区内物面压力大尺度脉动结构的特征尺度仍在增大,而图15(d)～图15(f)等值线则以向内收缩为主,对应的是外层压力大尺度脉动结构特征尺度的减小.

图15 站位3～ 5 两点相关系数(上:物面;下:y/δ=0.8;虚线:站位1)Fig.15 Two-point correlation coefficient at station 3～ 5 (top:wall;down:y/δ=0.8;dashed lines:station 1)

为了定量描述该现象,图16 进一步给出了干扰区内物面及外层压力脉动积分尺度的分布.积分尺度可以用来表征脉动场最大尺度结构的特征尺度,这里流向和展向积分尺度分别定义为[10]

可见,外层压力脉动流向和展向积分尺度远大于物面压力脉动,这与Duan 等[16]的研究结果是一致的.从整体变化规律上来看,外层压力Λx和Λz在站位2 处急剧增大,随后在站位3～ 5 则又缓慢减小,而物面压力的Λx和Λz整体上呈现逐步增大的趋势,这与图15 中相关性分析结果是相符的.从量值变化上来看,以图16(b)中的展向积分尺度Λz为例,对于外层压力,在激波的增强作用下,站位1 处Λz=0.73δ增大到站位2 处Λz=1.17δ,随后在再附区又缓慢衰减到站位5 的Λz=1.02δ;而物面压力,从站位1 处Λz=0.26δ逐步增大到站位5 处的Λz=0.82δ,后者约为前者的3.15 倍.可以看到,在站位5 处,此时再附区内物面压力脉动结构的展向和流向积分尺度已非常接近于外层的积分尺度.综上所述,可以认为这是由于外层压力脉动大尺度结构在过激波后呈现快速破碎的过程,而物面上大尺度结构仍在缓慢增长,这也解释了图12 中物面压力脉动低频能量的相对升高而外层压力脉动峰值频率往高频区的偏移现象.

图16 压力脉动积分尺度分布Fig.16 Distribution of integral scale for pressure fluctuations

5时空关联

进一步分析干扰区外层和物面压力脉动的时空关联特性,从而得到向下游传播的压力波的对流速度.

图17 分别给出了站位1 处物面和外层压力的时空关联云图,图中时空关联系数Rpp(Δx,Δt)定义为

式中Δt为延迟时间.以往研究结果表明[28,30],由于脉动结构往下游传播,物面压力时空关联云图呈现倾斜的类椭圆形分布,如图17(a)所示,本文计算结果也完全符合该结论.相较于物面压力,从图17(b)中可以看到,外层压力时空关联云图的椭圆率和倾斜角都明显增大.

图17 站位1 压力脉动时空关联系数Fig.17 Space-time correlation coefficient for pressure fluctuations at station 1

图18 和图19 分别给出了分离区站位2 和再附区站位3～ 5 压力脉动时空关联系数与站位1 (黑色虚线)的比较情况.如图18(a)所示,可以发现,相较于上游湍流边界层,分离区物面脉动压力时空关联曲线更为聚集,同时倾斜角也急剧减小,但整体形态仍满足椭圆形(只是等值线0.1 较为杂乱).值得注意的是,对于站位2 外层脉动压力,显然,图18(b)中等值线形态变化非常剧烈,已完全不符合椭圆形分布,特别是等值线0.1 和0.3,在时间方向上可近似为直线分布.造成这个现象的原因很可能是由于站位2 外层处于强激波区(见图9),在分离激波的强干扰下,此时泰勒冻结假设完全失效,脉动压力在时空上表征为强非线性.此外,从图19 中还可以看到,无论是再附区外层还是物面,其脉动压力时空关联等值线整体形态与上游站位1 都较为类似.从站位3～站位5,强相关等值线被急剧拉长压缩,且倾斜角增大,呈现逐步逼近上游站位1 的演化过程.从定性比较的角度来看,外层脉动压力时空关联曲线的恢复速度明显要快于物面.

图18 站位2 压力脉动时空关联系数(虚线:站位1)Fig.18 Space-time correlation coefficient for pressure fluctuations at station 2 (dashed lines:station 1)

图19 站位3～ 5 压力脉动时空关联系数(上:物面;下:y/δ=0.8;虚线:站位1)Fig.19 Space-time correlation coefficient at station 3～ 5 (top:wall;down:y/δ=0.8;dashed lines:station 1)

压力脉动对流速度的计算方法如下[16,31]:首先给定延迟时间Δt,随后依据时空关联曲线通过

获得流向间距rx,进而得到对流速度Uc=rx/Δt.采用上述方法计算得到了上游湍流边界层物面和外层脉动压力对流速度随延迟时间的分布情况.如图20 和图21 所示.总体来看,外层脉动压力对流速度要远大于物面脉动压力,这与Duan 等[16]的研究结果也较为符合.另外,对于外层脉动压力,可见,其对流速度对延迟时间的依赖性较小,基本维持在Uc=0.94U∞,与当地平均流向速度接近.然而,物面脉动压力对流速度随延迟时间则变化剧烈,在小延迟时间(小尺度脉动)时,对流速度约为Uc=0.56U∞,而在大延迟时间(大尺度脉动) 时,其增大到Uc=0.73U∞.已有研究表明[13],这是由于物面大尺度脉动压力与外层大尺度结构相关,而小尺度脉动压力则主要受近壁区结构影响,故前者对流速度明显大于后者.

从图20 和图21 中还可以看到,分离区物面脉动压力对流速度降低,约为Uc=0.2U∞,随后下游再附区则又逐渐增大,但仍小于上游湍流边界层.Na等[32]在强逆压力梯度平板的DNS 结果中也发现了类似规律,本文的计算结果进一步证实了该现象在入射激波/湍流边界层干扰区内的适用性.从图20(b)中可以看到,采用当地外层速度Ue对Uc进行归一化后,分离区和再附区各站位处物面脉动压力的对流速度分布均吻合较好,约为Uc=0.25Ue,与Bernardini 等[10]在跨声速激波/湍流边界层干扰下的研究结论也较为接近,但后者为初始分离(incipient separation)工况.如图21 所示,干扰区外层不同站位处脉动压力对流速度分布与上游结果较为类似,其对延迟时间的依赖性依然较弱,但其量值减小到约为Uc=0.4Ue.结合之前的功率谱和两点相关分析,可以推测外层大尺度脉动结构特征尺度的减小可能是造成该现象的重要因素.