朱志斌, 冯 峰, 沈 清

(中国航天空气动力技术研究院, 北京 100074)

引 言

边界层转捩会对飞行器摩擦阻力、 热流密度和流动分离再附特征等产生显著影响, 是高超声速飞行器研制中必须关注的气动现象. 但是由于影响流动转捩的因素众多, 且转捩过程复杂, 对边界层转捩特性和机制认识还十分有限, 严重制约了高超声速气动技术的进一步发展.

巡航、 滑翔类高超声速飞行器普遍采用扁平、 后掠的高升阻比几何构型, 受压力梯度和后掠角的共同作用, 飞行器表面边界层呈现显著的三维效应, 且边界层转捩通常受横流速度拐点导致的无黏失稳机制主导[1]. 这显著区别于传统平板、 圆锥等简单几何外形诱导的高超声速二维边界层流动, 二维边界层转捩一般由速度剖面上广义速度拐点引起的第1模态、 第2模态等高阶模态流向行波不稳定性引起[2]. 目前, 对高超声速边界层横流转捩过程和机理的认识不足, 导致对横流转捩的精确预测难度较大, 也难以准确解释横流失稳影响因素及作用规律[3].

美国和澳大利亚联合开展的高超声速国际飞行研究实验(Hypersonic International Flight Research Experimentation, HIFiRE)项目中设计的第5次飞行试验(HIFiRE5)专门针对高超声速边界层横流转捩问题进行研究. HIFiRE5项目以横纵比2∶1的椭圆锥构型为飞行载荷, 开展了大量的风洞实验、 数值分析工作,并于2016年5月成功完成飞行试验. 飞行试验前, Holden等[4]在LENS1风洞中, 采用全尺寸模型, 通过温敏漆测热技术, 观测到了起始于椭圆锥模型中心线及侧方下游的转捩形态. Berger等[5]采用全局磷光热成像技术对38.1%缩比模型表面热流进行了测试, 对比分析了迎角、 侧滑角及Reynolds数对热流分布的影响. Juliano等[6-7]以相同缩比HIFiRE5椭圆锥模型为对象, 在普渡大学Mach数为6的静音风洞中(Boeing/AFOSR Mach-6 Quiet Tunnel, BAM6QT), 采用温敏漆技术对模型表面温度进行整体测量, 研究了来流噪声、 壁面粗糙度、 迎角和Reynolds数对边界层转捩的影响. Borg等[8-10]通过压力传感器、 温敏漆和油流显示技术, 观测到HIFiRE5缩比模型在静音条件下的定常横流涡和行波扰动; 而在噪声条件下, 只通过油流图像观测到定常横流涡. 飞行试验由热流判定的转捩位置与前期地面试验结果整体定性一致, 如HIFiRE5b飞行器表面边界层出现观察到3个叶片状的转捩线形态, 转捩起始位置分别位于飞行器中心线、 侧前缘以及二者之间的中部区域[11].

HIFiRE5转捩预测的数值方法主要基于流动稳定性分析. 线性稳定性理论应用于三维外形的固有缺陷, Choudhari等[12]采用基于准平行性假设和曲率效应结合的稳定性分析方法, 对典型飞行和试验工况下HIRiRE5边界层转捩特征进行了预测. Gosse等[13]采用抛物化稳定性分析方法, 对HIFiRE5飞行弹道点下对称面的边界层流动进行了转捩预测分析. 研究发现, 侧缘平面边界层可能会出现第2模态的增长, 而中心线处流动受两侧旋涡影响, 呈现出更加复杂的不稳定性增长模式. Li等[14]采用经典的稳定性分析方法, 结合风洞试验数据, 对HIFiRE5缩比模型边界层转捩现象进行了分析, 并指出需要依靠对三维边界层转捩现象的直接数值模拟来改进经典稳定性理论在实际应用中的基础假设.

近期, 国外已有部分研究者采用非定常精细数值模拟方法对HIFiRE5风洞试验模型边界层转捩现象开展研究. Dinzl等[15]对风洞试验模型的定常横流涡进行了直接数值模拟, 采用分布式粗糙度来激发表面最不稳定扰动波, 并发现高热流条带是由流向速度扰动冲击壁面引起. Tufts等[16]采用高保真非定常精细模拟方法研究了HIFiRE5风洞试验模型边界层转捩现象, 利用平面声波线性叠加方法引入扰动, 通过调整扰动幅值复现了静音及噪声来流条件下的模型热流分布特征, 并获得了压力脉动能谱. 上述研究展现了精细数值模拟方法在高超声速边界层横流转捩问题研究中的良好应用潜力.

本文采用大涡模拟方法对椭圆锥边界层横流转捩现象进行了细致模拟分析, 深入认识高超声速边界层横流失稳转捩特征及机制, 为高超声速飞行器边界层转捩准确预测及气动设计奠定了基础.

1 计算模型及工况

以普渡大学静音风洞开展的椭圆锥风洞试验[6-7]为数值研究的参考对象, 即模型为38.1%缩比的HIFiRE5外形, 模型长度328 mm, 见图1. 风洞试验工况参数如表1所示, 本文只考虑0°迎角. 试验通过调节泄流槽的打开与闭合实现不同来流噪声条件, 在静来流条件湍流度低于0.05%, 而在噪声条件下达到约3%.

图1 椭圆锥试验模型[6-7]Fig. 1 Experimental model of the elliptic cone[6-7]

表1 风洞试验工况

2 流场模拟方法

2.1 控制方程及数值格式

数值计算采用了隐式大涡模拟方法[17-18], 即针对Favre滤波三维可压缩Navier-Stokes方程[19], 基于亚格子模型对不可解尺度流动建模求解假设, 依靠数值格式的耗散特征来抑制亚格子湍流动能积聚.

为精细刻画流场小尺度湍流结构, 同时无振荡捕捉激波间断, 流通量离散采用通量限制型偏心紧致格式[20], 并通过与TVD格式结合, 保证算法对激波间断的无振荡捕捉能力. 黏性项通过4阶中心差分计算. 时间推进采用2阶显式Runge-Kutta方法. 该大涡模拟方法已多次在高超声速复杂非定常流动中得到验证、 应用[21-22].

2.2 计算网格及边界条件

基于椭圆锥模型横、 纵向对称性, 采用1/4椭圆锥模型为对象进行计算建模, 以降低计算花费. 使用多块结构化网格划分计算域, 见图2, 流向计算域为82 mm≤x≤328 mm, 以便于在模型适当远的位置x=82 mm施加人工入口扰动, 相同的方法已在文献[15-16]中采用. 空间网格分布与脱体激波位置相匹配呈扁平型. 法向网格在近壁区进行加密, 壁面法向第1层网格尺度为5×10-4mm. 展向网格均匀分布, 流向网格尺度则逐渐增大. 流向、 法向和展向网格数分别为906×161×301, 对应的网格单元总量为4.344×107.

图2 计算网格示意图Fig. 2 Schematic of computational grids

大涡模拟以常规2阶Roe格式求解得到的全模型定常层流解为初始流场进行非定常推进计算. 计算域流向入口从全模型层流流场中截取得到. 为模拟真实情况下来流扰动和壁面粗糙度等对流场的干扰作用, 在入口处引入三维速度小扰动, 其形式与三维最不稳定T-S波相似, 曾在超声速喷流混合流场模拟中得到应用验证[23], 具体为

u′,v′,w′=AG(η)[cos(βθ±ωt)+ cos(0.5βθ±0.5ωt)]

其中,A为扰动幅值,G(η)为壁面距离η的函数, 形式为

扰动参数ω=2πf, 扰动频率f由当地边界层外缘速度Ue和边界层厚度δ确定, 即

f=Ue/(10δ)

θ为无量纲展向位置参数, 范围为0～1, 展向扰动取20个周期, 因此频率参数β为20×2π.

法向远场边界设置为来流条件. 物面采用等温无滑移边界条件, 壁温设为300 K. 展向对称面采用对称边界条件. 出口边界采用外插方法设置, 并对网格作缓冲层处理. 大涡模拟非定常计算无量纲时间步长为5×10-7, 特征长度L=1 m, 时间步长特征参考量为L/U∞, 待入口扰动波流出出口边界2倍周期的时间即流场充分发展后作流场统计, 统计步数为1×105.

2.3 数值结果验证

基于试验来流噪声条件, 入口扰动幅值分别设为0.05%和3%. 在核心计算区域内, 无量纲网格尺度分别为y+≈0.2,Δx+≈10,Δz+≈5, 已接近满足DNS网格分辨率要求[24], 即y+<1, 10≤Δx+≤20, 5≤Δz+≤10, 表明网格尺度达到大涡模拟要求.

图3给出了不同扰动幅值下数值结果与文献[6-7]中风洞试验结果的椭圆锥时均热流云图对比. 当小扰动幅值A=0.05%时, 见图3(a), 侧前缘和中心线间的中部区域出现两组细长条高热流条带结构. 在湍流度0.05%的试验静来流条件下, 见图3(c), 模型表面出现热流条带结构, 且中部区域高热流条带形态及位置均与数值模拟小扰动幅值A=0.05%的结果相符. 当扰动幅值为A=3%时, 见图3(b), 时均热流分布反映了边界层流动的横流转捩形态特征, 整体模型表面出现多个峰状转捩阵面, 中心线和侧前缘间的中部区域出现两个转捩阵面. 对应试验噪声来流条件下, 见图3(d), 椭圆锥模型边界层表面出现双肺叶状的转捩形态, 与计算结果对比可发现, 大扰动幅值计算得到热流分布反映了与试验观测数据相似的转捩形态特征. 此外, 风洞试验和时均计算结果在模型侧前缘区域均没有出现转捩形态特征. 整体而言, 本文大涡模拟能够有效预测椭圆锥边界层失稳转捩现象, 且获得了丰富的流场信息, 有助于获得高超声速边界层流动横流失稳特征的机理性认识.

(a) Numerical simulation: A=0.05%

(b) Numerical simulation: A=3%

(c) Experiment: quiet inflow

(d) Experiment: noisy inflow图3 数值结果与风洞试验热流数据对比Fig. 3 Heat-flux comparison between numerical results and wind tunnel experimental data

3 计算结果分析

3.1 流场结构特征

不加入口扰动计算得到的空间流场和壁面热流分布如图4所示. 可以看到在高超声速来流条件下, 椭圆锥模型不同周向角位置激波强度不同, 侧缘区域压强显著大于中心线区域, 形成的压力梯度使得流动由再附线向中心线流动. 在模型侧缘和中心线间的中部区域, 出现多条平行于来流方向的高热流条带.

图4 无扰动流场Fig. 4 Flowfield without inlet disturbance

图5展示了x=200 mm处流向截面流场. 可以看到流动从侧缘向中心线汇聚, 使得边界层从再附线到中心线逐渐增厚, 并在中心线附近形成蘑菇状形态的反向旋转旋涡结构. 在中心线两侧区域, 边界层内、 外部展向速度和流向涡量方向相反, 表明边界层外缘存在横向流动剪切作用.

(a) Density

(b) X component velocity

(c) Z component velocity

(d) X component vorticity

(e) Temperature图5 x=200 mm处流向截面流场Fig. 5 Flowfield of streamwise section at x=200 mm

图6给出了中心线附近不同展向位置处的速度分量对比. 从中可发现, 在0≤z≤8 mm范围内, 流向速度分量呈S分布, 展向速度分量呈反S分布. 在z=4 mm处速度剖面曲线扭曲最为严重. 而当z>8 mm时, 流向速度剖面恢复至常规边界层速度型, 展向速度剖面呈钩状分布, 出现速度拐点. 如图7所示, 中心线区域外侧的速度剖面显示, 模型中部区域展向速度剖面仍存在速度拐点, 而模型侧前缘处当地展向速度则接近为0, 流向速度均保持类似二维边界层分布形态. 根据速度剖面呈现的流动及潜在不稳定性特征, 可将椭圆锥边界层流场大致分为3个不同的区域, 即中心线区域、 侧缘区域以及中心线和侧缘间的中部区域. 在中心线区域, 流动汇聚产生流向涡, 流向、 展向速度剖面均呈明显扭曲形态, 流动非常容易失稳. 在中心线与侧缘之间的中部区域, 流向速度体现为黏性不稳定性, 但展向速度剖面仍存在拐点, 潜在表现为横流不稳定性. 在侧缘区域, 流动只表现为流向速度黏性不稳定性.

(a) X component velocity

(b) Z component velocity图6 中心线附近速度剖面Fig. 6 Velocity profile near the centerline

(a) z=25 mm

(b) y=0 mm图7 中心线外侧速度剖面Fig. 7 Velocity profile outside the centerline

3.2 扰动幅值影响

为研究横流转捩对扰动幅值的敏感性, 开展了不同扰动幅值下椭圆锥边界层转捩现象计算分析.图8展示了大涡模拟得到的瞬态涡系结构, 以速度梯度第2不变量等值面展示(Q2=0.01).由图可见, 不同扰动幅值下流场涡系结构具有明显差异. 当扰动幅值非常小(A≤0.05%)时, 中心线区域出现大尺度的流向涡结构, 在中心线和侧前缘间的中部区域出现若干条平行于来流方向的流向涡. 当扰动幅值大于0.3%后, 流场中出现显著的横流失稳转捩现象, 且扰动幅值越大, 流场涡系结构愈加丰富. 从瞬时涡系结构的分布形态看, 中心线区域涡结构尺度最大, 呈现出明显的失稳、 转捩及湍流发展演化特征. 侧前缘和中心线间的中部区域按转捩形态, 可细分为靠近中心线侧和靠近外缘侧两个子区域. 在较小的扰动幅值下(A=0.3%), 两个子区域易于区分, 靠近中心线侧的涡系结构尺度相对较大、 影响区域较小, 而靠近外缘侧区域的涡结构尺度小、 影响区域较大. 较大的扰动幅值下(A=3%), 中心线和中部区域间的涡系结构在展向逐步融合.

(a) A=0%

(b) A=0.05%

(c) A=0.3%

(d) A=3%图8 不同扰动幅值下瞬态涡系结构(Q2=0.01)Fig. 8 Instantaneous vortex structures under different disturbance amplitudes (Q2=0.01)

通过不同扰动幅值下瞬时和时均热流云图, 可进一步分析边界层流动失稳转捩特征. 见图9, 当扰动幅值非常小时(A≤0.05%), 微小的扰动也会迅速引起中心线区域旋涡结构的发展演化, 使得壁面产生对应的局部高热流分布形态. 而侧前缘和中心线间的中部区域出现两组细长条高热流条带结构, 对应展向坐标范围分别为20～25 mm和40～60 mm. 这两类热流条带瞬态值和时均值没有明显差异, 表明对应的流向涡结构是定常的. 当扰动幅值大于0.3%时, 热流结果反映出流场中存在显著的非定常流动结构, 在中心线区域流场中存在大尺度涡结构, 中心线和侧缘间的中部区域出现流向条带分布. 时均热流分布与小扰动幅值下定常流向涡结构位置相对应. 模型表面中心线处转捩起始位置最为靠前. 中心线和侧前缘间的中部区域出现两个转捩阵面, 其中靠近侧前缘一侧转捩形态随扰动幅值变化相对较小, 而靠近中心线一侧的转捩起始位置随扰动幅值增大迅速前移.

瞬态热流分布清晰地显示了流动失稳特征, 从中可发现, 在时均转捩阵面之前, 流场中存在斜向的条纹结构, 条纹结构方向与壁面摩擦力线一致, 表明对应的流动结构存在于边界层底层, 且尺度较小. 转捩后的热流条带整体形态平行于来流方向, 表明边界层流动由外缘处的大尺度流动结构主导. 中心线和侧缘之间的中部区域靠近中心线侧的流动受到中心线大尺度涡和横流扰动的共同影响, 发生失稳转捩位置相对更靠前, 产生的流动结构尺度也更大. 在侧前缘区域, 不同扰动幅值下的热流分布均没有出现转捩形态特征.

(a) Instantaneous: A=0%

(b) Time-averaged: A=0%

(c) Instantaneous: A=0.05%

(d) Time-averaged: A=0.05%

(e) Instantaneous: A=0.3%

(f) Time-averaged: A=0.3%

(g) Instantaneous: A=3%

(h) Time-averaged: A=3%

3.3 横流转捩特征

为认识静来流及噪声来流条件下高超声速边界层横流转捩特征,图10显示了0.05%和3%扰动幅值下流向截面的瞬时密度分布. 可以发现小幅值扰动下, 随流动向下游发展, 中心线处汇聚卷起的旋涡逐渐发生变形演化, 侧缘和中心线间的中部区域在靠近壁面处出现两组独立的涡结构. 靠近中心线一侧的涡数量少、 尺度相对较大, 并且其形成与主涡存在关联, 靠近侧缘的涡数量相对较多, 并且在靠近侧缘肩部位置逐渐减弱. 大扰动幅值下, 中心线处主涡结构失稳产生大尺度的流动结构, 中心线和侧缘间的中部区域非定常流动结构尺度相对较小, 且更贴近壁面.

(a) A=0.05%

(b) A=3%图10 流向截面瞬态密度Fig. 10 Instantaneous density at streamwise sections

图11给出了0.05%和3%来流扰动幅值下不同展向位置截面瞬态密度分布.A=0.05%小幅值扰动下, 模型后部中心线区域(z=0 mm)涡外缘位置出现了微弱的流场扰动, 而在展向中部(z=25 mm)及靠近侧缘处(z=50 mm), 没有出现扰动形态, 表明边界层流动处于稳定状态. 而A=3%扰动幅值下, 瞬态流场密度分布形态, 椭圆锥边界层流动向下游发展过程中出现了明显的失稳转捩现象.

(a) A=0.05%: z=0 mm

(b) A=3%: z=0 mm

(c) A=0.05%: z=20 mm

(d) A=3%: z=20 mm

(e) A=0.05%: z=50 mm

(f) A=3%: z=50 mm图11 展向截面瞬态密度Fig. 11 Instantaneous density at spanwise sections

图12对比了不同扰动幅值下的密度脉动均方根和湍动能分布, 以展现二者湍流发展差异. 在较小的入口扰动幅值下, 只有在中心线区域流场出现较为微弱的流场脉动, 中心线以外流场可认为是定常的. 而在入口扰动幅值较大时, 流场非定常脉动特征显著, 其中边界层外缘处脉动幅值较大. 结合瞬态流场涡系结构特征, 可确定较小扰动幅值下, 中心线和侧缘间的横流涡结构是定常的, 来流扰动提高会促发定常横流涡的二次失稳, 进而产生横流失稳转捩现象.

(a) Root mean square of density fluctuation

(b) Turbulent kinetic energy图12 不同扰动幅值密度脉动均方根、 湍动能对比Fig. 12 Comparison of root mean square of density fluctuation and turbulent kinetic energy between different disturbance amplitudes

3.4 非线性动力学分析

采用压强频谱分析、 相空间以及Lyapunov指数等方法对三维边界层的失稳转捩非线性动力学特性进行分析.

非定常统计过程中在椭圆锥表面布置了12个探测点, 其位置如图13所示. 探测点1～4位于中心对称线, 探测点5～8和9～12分别位于不同展向位置处.

图13 探测点位置示意图Fig. 13 Schematic of probe positions

图14为通过快速Fourier变换得到的压力脉动频谱曲线.图14(a)显示探测点1处脉动幅值较小, 对应于该点处于层流流动状态. 探测点2在低频200～400 kHz范围内存在较大脉动幅值. 而探测点3和4呈现出显著的宽频脉动特征, 频率范围约为200～1 000 kHz, 可推断当地流动已达到湍流状态.图14(b)和(c)脉动频谱对比发现, 在侧缘处(探测点12), 压强脉动幅值微小, 流场脉动微弱, 而在接近侧缘位置(探测点8和11), 压力开始出现较强震荡, 表明扰动进入发展放大阶段. 此外, 在转捩发展的起始位置处(探测点6和7), 压力脉动具有较低的特征频率频谱, 而在湍流区域(探测点5, 9, 10), 压力在宽频范围内较连续光滑地过渡.

(a) Probe 1～4

(b) Probe 5～8

(c) Probe 9～12

图15为各探测点的压强脉动相轨迹图, 其中横轴坐标为压强, 纵轴为其对时间的导数. 相轨迹曲线的拓扑结构可直观展示边界层转捩的非线性演化特征. 在探测点1, 压强脉动相图的范围较小, 可近似认为是一个不动点, 即该点流动接近于稳定状态. 在探测点2处, 相图轨迹在多个环形线圈间变换, 表明流动出现非线性演化及分叉, 呈现出向混沌状态转变趋势. 探测点3和4形成内外侧线圈缠绕结构的相图轨迹, 显示了流场中存在不同尺度、 多周期的非定常运动. 探测点6的相图轨迹形态与探测点2相似, 但压强导数幅值相对较小, 结合转捩演化过程结果可发现, 二者均处于转捩的起始区域. 探测点5和7分别处于模型中部内外侧横流转捩发展区域, 形成围绕多个中心的线圈缠绕相图轨迹形态. 探测点8的相图曲线围绕单个空心的多重线圈形态, 且相图范围较小, 可推断流场中存在多周期的小尺度流动结构. 探测点9和10处于全湍流区域, 其相图轨迹也与探测点3和4相似. 探测点11靠近模型侧缘, 压强导数相对较小, 说明对应的流动结构的尺度较小. 探测点12的相图范围较小, 且呈现规则的多重线圈, 表明此侧缘处流动结构还未开始非线性演化.

(a) Probe 1

(b) Probe 2

(c) Probe 3

(d) Probe 4

(e) Probe 5

(f) Probe 6

(g) Probe 7

(h) Probe 8

(i) Probe 9

(j) Probe 10

(k) Probe 11

(l) Probe 12

Lyapunov指数可定量评价动力系统非线性特征, 表征了系统在相空间中相邻轨道间收敛或发散的平均指数律. 最大Lyapunov指数大于零时, 意味着在系统相空间中, 初始两条轨线间的差别会随着时间的演化而呈指数律的增加以致达到动力学混沌状态[25]. 采用Wolf等[26]提出的基于轨道跟踪的Lyapunov指数计算方法, 从各探测点的压强时间序列提取了最大Lyapunov指数. 如图16所示, 在中心线处(探测点1～4, 实线所示), 最大Lyapunov指数(Lyapunav exponent, LE1)由接近零急剧增大, 表明中心线处流动扰动经历强烈的非线性增长, 流动出现显著的混沌特征. 最大Lyapunov指数沿展向的分布(虚线及点划线)显示, 模型中部区域压强脉动的LE1值明显低于中心线处. 靠近中心线的内侧(探测点6和9)相对较大的Lyapunov指数值表明当地流动存在较强非线性作用, 而在靠近侧缘的外侧区域, Lyapunov指数值较小, 表明当地流动非线性特征较弱.

图16 压强时间序列最大Lyapunov指数Fig. 16 Maximal Lyapunov exponent of pressure time series

4 研究结论

本文采用隐式大涡模拟方法, 通过在计算域入口引入三维速度扰动, 对不同扰动幅值下的HIFiRE5风洞试验椭圆锥模型的横流失稳转捩现象开展了细致数值研究, 得到如下结论:

(1)大涡模拟精细模拟了高超声速边界层横流失稳转捩的非定常演化过程, 预测了静音及噪声来流条件下边界层转捩的形态特征, 获得了与试验数据接近的热流分布形态.

(2)椭圆锥流场中心线流动汇聚形成的流向涡结构非常容易失稳, 在中心线及侧缘间的中部区域存在显著的横流不稳定性, 两种转捩机制共同影响三维边界层转捩过程.

(3)扰动幅值对椭圆锥三维边界层转捩影响显著. 在静来流条件下, 横流区域出现两组独立的定常横流涡结构; 噪声来流条件下, 中心线主涡和中部横流涡发生失稳转捩, 在模型表面形成多峰状的转捩阵面.

(4)探测点压强脉动的非线性动力学分析, 展示出三维边界层发生失稳转捩的非线性演化特征, 深化解释了模型表面边界层内周期性小扰动发展至混沌的非线性动力学演化过程.