横向矩形微槽抑制高超声速第二模态扰动波的参数化研究*

2022-10-16刘勇涂国华向星皓李晓虎郭启龙万兵兵

物理学报 2022年19期

刘勇涂国华向星皓李晓虎郭启龙万兵兵

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

针对高超声速飞行器边界层转捩控制问题,以马赫数6 平板边界层的第二模态扰动波为研究对象,采用线性稳定性理论(LST)和直接数值模拟(DNS)分别开展了离散模态的同步模式研究和大尺寸(0.4 mm 宽)横向矩形微槽开槽位置对第二模态扰动波的控制作用研究.LST 分析表明: 涡/熵波会导致Mack 第二模态和“mode I”模态(通常来源于快声波)的分支类型发生改变.通过DNS 发现,开槽表面对基本流的影响程度与边界层流向位置(或厚度)相关,随着开槽位置后移(边界层厚度增加),开槽表面对基本流动的影响程度减弱,摩擦阻力系数和压差阻力系数也逐渐减小.DNS 结果还表明,位于快/慢模态同步区间之前的开槽工况对第二模态扰动波依然有抑制效果,这与文献中关于小尺寸(微米量级)微孔隙位置对第二模态控制作用的结论不同,同时发现,当矩形微槽布置在最大增长率区间范围内或快/慢模态同步区间位置时,对第二模态扰动的抑制效果最佳.

1 引言

通用航空飞行器(CAV)是满足许多国家愿景和军事计划的新型运载工具,其面临的一个主要困难是: 在承载有效载荷下,如何从太空中以高超声速再入大气层而不被烧毁[1].之前,大量研究表明飞行器高速飞行时边界层由层流到湍流转捩会导致壁面摩阻和热流成倍增加[2],有效载荷相差可高达一倍[3],严重影响飞行安全和飞行效率.因此,高超声速飞行器迫切需要延迟边界层转捩,尽可能保持飞行器的绕流为层流状态.

但是,转捩是一个复杂的物理演变过程,它强力地依赖基本流和额外的扰动.转捩的具体过程可以根据环境扰动的量级分为五个路径[4,5],这5 条路径中有多种失稳形式,比如模态增长(第一、二模态、横流失稳等)、瞬态增长、多参数不稳定、旁路转捩等.当飞行器在大气中飞行时,转捩物理过程主要包括感受性、边界层内扰动波线性失稳和扰动波非线性作用三个阶段.在线性失稳区(即第二个阶段),扰动波的幅值通常增加数百倍(风洞工况)乃至数百万倍(飞行工况),空间区域可占层流区的50%以上,因而,抑制扰动波的线性增长可有效地延迟边界层转捩,从而减小飞行器表面摩擦阻力和热流.高超声速飞行器与低速飞行器相比,转捩机制更为复杂,最为典型的差别就是高超声速边界层转捩过程中除了存在低速边界层常见的第一模态、横流失稳等转捩机制,还存在第二模态和其它更高阶模态,它们的扰动频率较高,表现出声学性质的无黏不稳定性.对于大气中飞行的高超声速飞行器,由于壁温比(壁温与总温之比)低,第一模态扰动会被抑制,但是会促进第二模态的发展.因此,在高超声速飞行条件下,第二模态是导致转捩的主要原因之一,它普遍存在于轴对称体(比如导弹)、平面体(比如超燃冲压发动机入口前的压缩面)和通用飞行器的附着区(比如机翼前缘).因此,抑制第二模态能够有效延迟边界层转捩,维持飞行器表面的层流状态,使摩擦阻力和热流降为湍流的1/5 至1/3,实现飞行器降热减阻和增程增载的目的.

微孔隙表面在对基本流场没有显著影响的情况下,能通过吸声机制有效地抑制第二模态扰动的发展,是最有可能投入工程实践的技术之一.Malmuth 等[6]最早发现微孔隙表面能够抑制第二模态扰动的增长,并从无黏扰动的特征值问题出发给出了一个简单的模型来预测微孔隙表面对第二模态扰动的抑制作用.随后,Fedorov 等[7]和Rasheed等[8]对具有等间距圆柱形盲孔的多孔涂层覆盖的高超声速边界层进行了第二模态扰动波的稳定性分析,发现多孔涂层对高频扰动声波能量的吸收使第二模态扰动的幅值增长率大幅度减小,并通过尖锥的高超声速风洞实验证实了微孔隙能显著延迟转捩.Chokani 等[9]采用双谱分析的方法发现多孔涂层能显著减弱第二模态的次谐波和谐波共振现象,非线性相互作用被大大降低.Egorov 等[10,11]采用直接数值模拟(DNS)对平板、尖锥和压缩拐角表面的多孔涂层感受性问题进行了相关研究,研究发现多孔涂层对与感受性相关的边界层模态初始振幅和声波扰动影响较小,但能强烈抑制第二模态扰动幅值.Dong 和Li[12]采用DNS 研究了第二模态和自由流声波在凹腔表面的局部感受性,发现当凹腔较浅时,局部散射和感受性效应随凹腔深度的增加而增加,对于深凹腔则相反.Long 等[13]采用Doak 动量势理论将动量密度分解为涡、声和热分量来研究多孔表面的能量输运,他们发现多孔涂层表面的正涡源被显著抑制,较少的旋涡能量被传输到临界层,由于波动总焓向外传输中的能量损失,声能最终耗尽,声辐射从多孔涂层表面消失.

文献研究还发现,影响多孔涂层对第二模态扰动控制效果的主要因素有: 孔隙形状参数[7,14-19](孔隙率和宽深比)、涂层分布位置[20-22]和壁面温度[23]等.Brès 等[24]发现当多孔涂层表面出现较大的宽深比和高的孔隙率时,会在涂层表面出现新的不稳定模态,这种不稳定性与腔体共振相关,可能比第二模态更不稳定.基于第二模态的声学特性,Zhao等[25,26]设计了近零阻抗孔隙结构和反射方向可控的超表面,其通过广义斯涅尔定律对扰动波的反射方向进行调控,使入射声波和反射声波在表面处相位相反,强度相互抵消,实现抑制第二模态的目的.在超声速边界层内,慢模态来源于慢声波,快模态来源于快声波或涡/熵波,当快模态和慢模态相速度同步时,会引发相应的模态转换机制[5,27],同时,利用基于空间双正交分解系统的多模态分解技术[28],可以将扰动场分解为简正模态,从而对不同模态扰动演化分别开展物理机制研究.快模态和慢模态相速度同步点与多孔涂层分布位置间的相对关系,对第二模态的控制作用有紧密联系,当多孔涂层布置在同步点下游时,第二模态扰动被抑制,反之,则可能被促进[20-22].多孔涂层表面结构尺寸通常为微米量级,不易于加工,于是郭启龙等[17]开展了较大尺寸(0.1 mm 量级)的矩形微槽数值模拟,发现大尺寸的横向微槽同样对第二模态抑制效果,随后,Liu 等[29]通过数值模拟发现横向矩形微槽的微小局部加工形变有可能会增强对第二模态扰动的抑制作用.郭启龙等[17]和Liu 等[29]给出的微槽道附近的流动图像显示槽道内出现了回流区,槽道外出现了声波散射,因此较大孔隙的多孔涂层与微米孔隙的多孔涂层相比,除了具备吸声机制外,可能还存在剪切耗散机制和声波散射机制,而文献中常用的声阻抗边界模型尚不能充分体现这两种机制,需要采用高精度DNS 等更精确的方法来研究两种机制对第二模态扰动的影响.

目前,在高超声速条件下多孔涂层分布位置对第二模态扰动的影响研究大都基于微米量级的孔或槽,对于大尺寸(毫米量级)矩形微槽开槽位置对第二模态扰动的影响没有相关结论.本文在郭启龙等[17]和Liu 等[29]的工作基础上首次采用线性稳定性理论(LST)和DNS 分别开展了第二模态扰动的流向演化与横向矩形微槽的开槽位置研究.针对马赫6的空间发展高超声速平板边界层,对比了微槽道平板与光滑平板的表面边界层和第二模态扰动波的演化情况,获得了快慢模态同步区间与开槽位置相对关系对第二模态的影响规律.

2 数值方法

2.1 空间线性稳定性分析

在进行矩形微槽转捩控制研究之前,首先需要对边界层稳定性特征和转捩特性有足够的认识.线性稳定性理论是最常用的研究高超声速边界层失稳特性的理论之一,能很好的分析外来扰动是如何在边界层中演化发展.线性稳定性理论在线化小扰动假设条件下,将瞬态变量q分解为基本量q0和扰动量q′.对于二维平板边界层,由于流动只在法向变化剧烈,在流向是慢变的,基于局部平行假设,扰动形式可以写成:

其中(y) 表征的是扰动的形状函数; c.c.表示共轭复数;对于空间模式,α是复数,用α=αr+iαi来表示,αr是流向波数,当地相速度被表示为c=ωr/αr,αi指扰动的增长率,如果αi＜0,则流动是不稳定的.值得注意的是,LST 由于引入平行流假设,其所预测的扰动增长率与DNS 计算结果有轻微出入[30-33].

本文选择Wartemann 等[34]的平板稳定性分析工况对现有的LST 程序进行验证,在空间模式下,通过与文献对比增长率 −αi随频率ωr的变化曲线,验证程序的准确性已满足计算要求(如图1 所示).

图1 增长率 −αi 随频率 ωr的变化Fig.1.Growth rate −αi versus frequency ωr .

2.2 直接数值模拟

文献通过声阻抗边界与LST 相结合的方法来研究微槽道对边界层稳定性的影响,如Fedorov 等[7]和Zhao 等[22]的做法,但是这种方法忽略微槽道的内部流动,存在较大误差.采用高精度算法直接求解二维Navier-Stokes 方程(即DNS)是一种更精确的方法.DNS的控制方程是曲线坐标系下的Navier-Stokes 方程,它是通过对笛卡尔直角坐标系(x,y,t)下的Navier-Stokes 方程进行坐标变换和无量纲化处理后获得.二维任意曲线坐标系(ξ,η,τ)下的控制方程可表示为

这里J表示雅可比矩阵,其可被表示为

其他变量分别被表示如下:

本文采用高精度有限差分格式对控制方程进行离散,无黏项采用4 阶WENN 格式[30],黏性项采用4 阶中心差分格式.对于带有矩形微槽的平板边界层基本流计算时,可看成定常问题来求解[7,17];当平板前缘施加第二模态扰动时,采用非定常时间离散格式来求解控制方程.定常问题计算的时间离散采用LUSGS 格式,非定常问题计算的时间离散采用3 阶Runge-Kutta 格式.

3 高超声速边界层离散模态的同步模式

本文选取的来流参数与文献[16-18,29]相同,如表1 所列.基于表1的来流参数,使用Blausis解程序生成平板边界层基本流场,采用线性稳定性理论(LST)来研究高超声速边界层内快慢模态和涡/熵离散模态之间的相互作用.高超声速平板边界层感受性问题在之前已经得到广泛的关注,先前的研究已经发现不同的分支模式和模态转换有着紧密的联系,快模态和慢模态之间的同步会导致不同的模态分支拓扑结构[35-37].首先考察第二模态在不同流向位置上的稳定性特征,分别选取流向50,200 和1000 mm 处来研究模态转换机制随流向位置的变化.图2 给出三个不同流向站位的离散谱信息,其中图的左半部分为扰动波相速度c随频率ωr的变化,右半部分为扰动增长率 −αi随频率ωr的变化,图中相速度1±1/M∞分别表示慢声波和快声波,相速度1 表示熵、涡波.随着流向位置的后移,快声波激发出一系列的快模态,记为mode I,mode II···;相应的慢声波激发出一个慢模态.从图中可发现,慢模态和快模态的模态分支类型随流向位置的不同发生了变化.根据Mack 模态(也称第二模态)和mode I 相速度的分支类型(“P”表示两者平行,“X”表示两者交叉,“S”表示Mack 模态与慢模态属于同一支,“F”表示Mack 模态与快模态属于同一支),发现图2(a)、图2(b)、图2(c)分别属于“X-S”型、“P-S”型,和“X-F”型.即,在上游Mack 模态属于前缘慢模态在高频区域的分支,并与一系列快模态发生同步共振(“X-S”型);随着向下游移动,慢模态分支在同步点之后产生跳跃,模态分支类型由“X-S”型转变为“P-S”型,同时,同步点也拓展为一个小的同步区间(棕色方框内);继续下游,Mack 模态从慢模态分支跳转到快模态分支,并且其总能与mode I产生一个小的同步区间,此时,模态分支类型从“P-S”型跳转为“X-F”型分支.

图2 马赫数6 平板边界层离散模态的流向演化;(a) x=50 mm,(b) x=200 mm,(c) x=1000 mmFig.2.Streamwise evolution of discrete modes of the Ma6 flat plate boundary layer: (a) x=50 mm,(b) x=200 mm,(c)x=1000mm.

表1 来流参数设置Table 1.Free stream parameter setting.

为了进一步解释Mack 模态和mode I 模态转换机制,图3 绘制了高超声速边界层内三个流向站位(50,200 和1000 mm)的离散谱和连续谱分布.在高超声速平板前缘,涡波和熵波对Mack 模态和mode I 影响较大,Mack 模态被涡/熵波“阻隔”,相速度小于1,属于mode S 分支.mode I 被涡/熵波“吸引”,相速度从快声波穿过涡/熵波与Mack模态产生同步.随着流向位置的后移,涡波和熵波逐渐消失,Mack 模态迅速跳转为快模态,而mode I变为慢模态,两者又继续产生同步现象.

图3 不同流向位置的离散谱和连续谱分 (a) x=50 mm;(b) x=200 mm;(c) x=1000 mmFig.3.Discrete and continuum spectrum at different streamwise locations: (a) x=50 mm;(b) x=200 mm;(c) x=1000 mm.

4 数值实验

4.1 线性稳定性分析

基于表1 来流条件得到光滑平板的层流场,采用LST 对矩形微槽进行辅助设计.图4 给出了Mack 模态中性曲线和不同频率的N值曲线,本文选取400 kHz为第二模态特征频率,其对应无量纲频率ωr约为1.42.图4(a)给出了三套网格下的计算得到的中性曲线,可见它们基本完全重合,验证了计算结果的网格无关性,同时也得到400 kHz 扰动波的失稳区间为[158 mm,260 mm].由图4(b)可见,取400 kHz 扰动波的最大N值对应的流向位置为260 mm.图5 给出该扰动波的流向速度和温度扰动型函数.图6(a)和图6(b)分别给出400 kHz 快/慢模态增长率和相速度沿流向的变化.由图6(a)发现慢模态对应的增长率(−αi)大于0,所以400 kHz的Mack 模态扰动波来源于慢模态,由图6(b)发现400 kHz 快慢模态对应的同步区间范围大概在流向x/L0∈[154,261] 处.

图4 Mack 模态的线性稳定性分析 (a) 中性曲线;(b) N 值曲线Fig.4.Linear stability analysis of Mack modes: (a) Neutral curves;(b) N-value curves.

图5 扰动型函数 (a) 流向速度实部和虚部;(b) 温度实部和虚部Fig.5.Perturbation shape function: (a) Streamwise velocity real and imaginary parts;(b) temperature real and imaginary parts.

图6 400 kHz的Mack 模态沿流向发展情况 (a) 增长率 −αi ;(b) 相速度cFig.6.Development of the 400 kHz Mack mode perturbation along the flow direction: (a) Growth rate −αi ;(b) phase velocity c .

4.2 DNS 方案

为了研究开槽位置xl对矩形微槽抑制Mack模态的影响,设置如图7 所示的物理模型,其中A1,A2,A3,A4,A5 分别表示五个不同开槽工况,其中A1 布置在同步区间之前,A2,A3,A4 布置在同步区间之内,A5 布置在同步区间之后,此外,还有未开槽的平板工况A0 作为对照组,具体的开槽位置布置见表2.DNS的流向计算区域为100 至300 mm,法向计算区域为0 至54 mm(约为边界层厚度的30 倍),由于A4,A5 工况开槽区靠近超声速出口,为了避免扰动反射对流场的影响和透射系数能恢复到定值,选择将A4,A5 工况流向计算域加长至400 mm.A1—A3 工况流向布置2600 个网格点,对于加长的A4 和A5 工况流向布置3100个网格点,法向布置301 个网格点,法向网格以“tanh”函数沿着边界层向外增长,其中第一层网格间距为0.0008.流向布置20 个微槽,微槽宽0.4 mm,槽深1.392 mm,开槽率0.5,开槽区总长16 mm保持不变,单个微槽采用20×30 个网格单元来识别槽内流动,微槽内流向采用均匀分布,法向以“tanh”函数向下增长,第一层网格间距为0.0008.边界设置为: 左边界采用高超声速入口条件,右边界超声速无反射出口条件,上边界采用远场条件,下边界采用无滑移等温壁面条件.本文数值模拟采取的方案是: 先利用高精度数值模拟软件计算当前来流条件下平板边界层的基本流场,通过截取流向100 mm 站位的剖面作为入口条件,然后计算开槽工况的基本流场,待基本流场计算收敛后,在入口施加400 kHz的第二模态扰动波来研究其在微槽表面的演化.

图7 开槽位置Fig.7.Grooving locations.

表2 开槽位置参数Table 2.Grooving location parameters.

4.3 网格收敛性验证

前文已经开展了线性稳定性分析的网格收敛性验证,此处补充DNS的网格收敛性测试.选择A1 工况,对流向和法向网格加密一倍并保证网格分布同比例变化,时间步长减小一半,因此网格加密后相同区域的计算量将增加8 倍,为了减少计算量,网格加密后的流向计算域从[100 mm,300 mm]缩减为[100 mm,150 mm].图8 对比了标准网格和加密网格计算的基本流在流向140 mm 处的无量纲流向速度和温度剖面,从图中发现标准网格和加密网格基本流剖面符合得很好.图9 和图10 分别对比了A1 工况两套网格在施加400 kHz 第二模态扰动的瞬态流场和扰动幅值,可发现两套网格下的瞬态流场和扰动幅值基本一致,表明计算结果是网格无关的,具有相当高的可信度.

图8 流向140 mm 处基本流剖面 (a) 无量纲流向速度;(b) 无量纲温度Fig.8.Basic flow profile at the 140 mm streamwise location: (a) Dimensionless streamwise velocity;(b) dimensionless temperature.

图9 施加第二模态扰动后的瞬态压力云图 (a) 标准网格;(b) 加密网格Fig.9.Transient pressure contours after imposing the second-mode perturbation: (a) Standard grid;(b) fine grids.

图10 标准网格和加密网格上400 kHz 第二模态扰动幅值沿流向演化.Fig.10.Evolution of the 400 kHz second-mode amplitude in the standard grid and the fine grid.

5 结果与讨论

5.1 开槽位置对基本流的影响

图11 给出了基本流的压力云图,从图中发现高超声速平板前缘的压力明显大于平板后缘,矩形微槽对基本流场的影响取决于背景流场的强弱,其在边界层内产生压缩/膨胀波系结构并以tan−1(1/Ma)倾角向边界层外势流区辐射.在壁面以上(y≥0),平均流修正参量可以定义为Δφ=φ −φref,其中φ为施加微槽的基本流,φref表示光滑平板的基本流.图12(a)和图12(b)分别给出了A1—A5 工况首槽和尾槽中心线处基本流压力修正量 Δp的分布.通过对比发现随着开槽位置的后移首槽和尾槽中心线处压力修正量 Δp有减弱趋势,主要原因是边界层沿流向逐渐增厚.为了进一步研究开槽位置对基本流场的影响,图13 给出了5 个工况流场摩擦阻力系数Cdf、压差阻力系数Cdp和总的阻力系数Cdx(包括摩擦阻力系数Cdf和压差阻力系数Cdp)随开槽起始位置xl的变化.从图13 中发现,随着开槽位置的后移Cdx,Cdf和Cdp逐渐减小,主要原因仍然与边界层厚度有关,边界层沿流增厚后,法向梯度减小,从而摩阻减小.由于微槽道对压力的修正沿流向逐渐减弱,因而微槽后移时,压差阻力也逐渐减弱.光滑平板的压差阻力等于零,但是,由于每个微槽的左侧壁面和右侧壁面的压力出现微弱差异(微弱的逆压梯度),导致微槽道表面的压差阻力不为零,其大小比摩擦阻力低三个量级.

图11 基本流压力云图Fig.11.Pressure contours of the base flow.

图12 A1—A5 工况基本流修正压力剖面分布 (a) 首槽中心线;(b) 尾槽中心线Fig.12.Pressure correction of the A1—A5 basic flow: (a) First groove centerline;(b) tail groove centerline.

图13 阻力系数(Cdx,Cdf,Cdp)随开槽位置 xl的变化Fig.13.Drag coefficients (Cdx,Cdf,Cdp) versus grooving locations xl .

5.2 开槽位置对第二模态扰动的影响

在平板入口处添加第二模态扰动波,图14 给出了五个开槽工况瞬态脉动压p′云图,可发现第二模态扰动具有流向行波特征,扰动波在向下游的传播过程中增长,在开槽区域,压力脉动被抑制.图15 给出A0—A5 在y=0(壁面)的脉动压力沿流向的分布,与光滑平板相比,开槽能够有效降低下游压力脉动的幅值,其中A3 工况的幅值抑制率最高,达到59%,A1 工况的幅值抑制率最低,接近35.8%.

图14 脉动压力云图 (a)—(f) A0—A5 工况Fig.14.Pulsating pressure: (a)—(f) A0—A5 cases in turn.

图15 壁面脉动压力沿流向分布 (a)—(e) A1—A5 工况壁面脉动压力Fig.15.Distributions of the wall pressure fluctuation: (a)—(e) The A1—A5 cases.

为了比较开槽位置对第二模态扰动波控制效果的影响,采用(6)式来提取扰动场的第二模态扰动幅值.式中u′(x,y,t) 表示通过数值模拟得到的扰动速度,F[·] 表示Fourier 变换,(x,y,f)为模态的本征函数.ε为扰动的初始幅值,本文取0.0001,max[·]表示提取扰动法向最大值,采用(7)式来对扰动幅值进行归一化.

图16 给出400 kHz 第二模态扰动波的归一化幅值A沿流向的变化,其中图16(a)和图16(b)分别表示标准计算域的归一化扰动幅值和在加长区的归一化扰动幅值.从图16(a)可以发现在x/L0∈[100,115]处扰动波有轻微的增长现象,这主要是因为入口添加的第二模态扰动波是在基本流的法向速度等于零(平行流)的假设下获得的,而真实流动具有微小的法向速度,因此入口存在一个调制过程,同样的现象在类似的文献中也有出现.图16(a)还显示,在x/L0∈[170,272]的区域内,扰动波是增长的,增长区域与LST 得到的中性曲线包含的区域(图4(a))相比延后了12 mm,导致该差异的主要原因可能是LST 采用了平行流假设,而真实流动是非平行的,但是LST 与DNS的结果偏差比较小,比如失稳位置偏差小于7%、失稳区域长度偏差约为0,因而相互验证了对方的正确性.

图16 400 kHz 第二模态沿流向的发展 (a) A0—A5 工况流向[100,300] mm 扰动幅值;(b) A0,A4,A5 工况在加长区[300,400] mm 扰动幅值Fig.16.Development of the 400 kHz second mode along the streamwise direction: (a) A0—A5 cases flow direction[100,300] mm disturbance amplitude;(b) A0,A4,and A5 cases in the extended area of [300,400] mm.

开槽平板与光滑平板相比,在开槽区域,扰动幅值出现一定波动,但是总体趋势是减小,在槽区结束处,压力脉动比光滑平板小.在槽的下游,即使表面没有微槽道,压力脉动的幅值仍然比光滑平板小,但增长趋势与光滑平板一致,表明在微槽道对脉动的影响可以持续至下游区域,对基本流的影响逐渐减弱.

为了定量地描述矩形微槽对边界层扰动的影响,本文采用开槽结束后60 mm 处的扰动幅值与对应站位平板幅值之比作为矩形微槽的透射系数T,如(8)式所示,其中xe表示开槽结束位置.图17给出了透射系数T 随开槽起始位置xl分布,从图中发现A3 工况的透射系数相较于其他工况是最小的.从图15—图17 都可以看出,A3 工况对第二模态扰动波的抑制效果最强,其背后的物理本质是A3 工况的开槽位置正好是第二模态扰动波增长率最大的位置(见图6(a)),也是快/慢模态同步区域的中心位置.矩形槽道一方面通过吸波机制抑制第二模态扰动波增长,另一方面也减弱了快/慢模态的同步强度,根据Fedorov[5]的结论,快/慢模态同步是产生第二模态扰动波的主要机制,因而,减弱同步强度有利于抑制第二模态扰动波.