V 字形钝前缘激波反射迟滞现象

2022-08-30张志雨钱战森李祝飞

南京航空航天大学学报 2022年4期

王军，张志雨，刘愿，钱战森，李祝飞

（1.中国科学技术大学近代力学系，合肥 230027；2.中国航空工业空气动力研究院，沈阳 110034）

吸气式宽速域飞行器从地面起飞逐渐加速至高超声速以及返航过程中，经历了多个气动环境迥异的飞行阶段，给气动外形设计提出了巨大挑战［1］。特别是，宽速域飞行条件下的激波干扰更加复杂多变，甚至与飞行阶段的历史效应紧密相关，认识其变化规律对工程设计尤为重要。

三维内转式进气道［2-3］具有较好的宽速域适应能力，受到广泛关注，这类进气道在V 字形唇口根部位置产生的复杂激波干扰，容易导致严酷的气动力/热问题［4-6］。为了深入认识V 字形唇口部位的复杂流动，肖丰收等［7-8］最早提炼出由半径比R/r（根部倒圆半径R和前缘钝化半径r）和半扩张角β表征的V 字形钝前缘模型，在来流马赫数Ma∞=6条件下，发现不同几何参数下V 字形根部主要产生异侧激波规则反射（Regular reflection，RR）、马赫反射（Mach reflection，MR）以及同侧激波规则反射（Regular reflection from the same family，sRR）3 种类型。进一步的研究表明，V 字形钝前缘的激波反射类型对流场非定常振荡特性［9-13］、壁面气动力/热载荷［13-17］、下游流场演化［18-21］以及热防护优化设计［22-23］等都起着关键作用。肖丰收等［7］和蒙泽威等［14］对V 字形钝前缘的研究表明，反射类型在不同来流马赫数下的变化明显，激波干扰结构在宽速域条件下将变得更加复杂。值得注意的是，当飞行器经历不同速域时，可能产生激波反射类型双解以及转变迟滞现象，一旦出现双解，极容易引发流场结构突变，并给壁面气动载荷预测带来不确定性。张志雨［13］通过数值模拟初步证实了存在V 字形钝前缘激波反射迟滞现象，但对迟滞过程中反射类型转变机制的认知，仍不明晰。

本文采用数值模拟方法，并辅以风洞实验，通过改变Ma∞，研究V 字形钝前缘激波反射迟滞现象，以期阐明迟滞转变机制以及迟滞现象对壁面气动载荷影响，为内转式进气道V 字形唇口设计提供有价值的参考。

1 模型与方法

1.1 V 字形钝前缘模型

V 字形钝前缘模型如图1 所示，由两侧斜掠直前缘和根部倒圆区域组成，模型半扩张角β=18°，倒圆半径R=15 mm，前缘钝化半径r=15 mm，以及斜掠直前缘长L=21r。其中，x方向为流向，y方向为横向，z方向为展向，φ为V 字形根部倒圆区域的周向角。

图1 V 字形钝前缘模型示意图Fig.1 Schematic of the V-shaped blunt leading edge model

1.2 风洞实验方法

实验在中国空气动力研究与发展中心的FL-31 下吹式常规高超声速风洞［13］中进行，来流马赫数Ma∞=6，静温T∞=55.5 K，静压p∞=633.4 Pa。此外，风洞配备了高速纹影系统来捕捉激波干扰结构，纹影的拍摄速率为16 kHz，曝光时间为62.5 μs。

1.3 数值模拟方法

数值求解基于雷诺平均的三维Navier-Stokes控制方程，无黏通量采用Roe 格式进行差分分裂，对流项采用二阶迎风格式离散，黏性项采用二阶中心差分格式离散，使用k-ωSST 湍流模型，空气采用量热完全气体假设，分子黏性系数由Sutherland公式计算。计算域、边界条件以及驻点附近网格分布如图2 所示，鉴于来流和几何的对称性，为提升计算效率，采用1/4 计算域进行模拟，来流采用压力远场条件、压力出口、对称边界（y=0 和z=0 平面）和无滑移绝热固壁边界条件。为了便于同本文风洞实验结果对照，来流条件参照FL-31 风洞参数进行设置。

图2 计算域、边界条件和驻点附近网格Fig.2 Computational domain, boundary conditions and surface mesh near the stagnation point

计算域采用结构化六面体网格离散，并在近壁面区域进行了网格加密。以残差下降4 个数量级，或者残差不再变化同时驻点位置流场参数保持稳定作为收敛判据。为了进行网格无关性验证，如图3 所示，采用了4 套不同疏密程度的网格，提取壁面中心线上的压力进行对比，其中压力p采用相同来流条件下正激波波后的总压p0进行无量纲化。可以看出，采用第3 套和第4 套网格获得的压力结果几乎重合，最大的p/p0相差不超过0.4%，表明本文采用的数值模拟方法能够满足需求。因此，采用第3 套网格（ζ × ξ × η=450×300×90）进行后续研究，总网格量约为1 200 万个，壁面网格第一层高度为1 μm，保持壁面y+＜1。上述模拟方法，在笔者前期的研究中已得到广泛地验证［8，16，18-21］，在本文2.1 节将进一步结合风洞实验考核该方法的可靠性。

图3 壁面中心线压力对比Fig.3 Comparison of surface pressure along the centerline

在计算过程中，首先，分别在Ma∞=5.7 和6.7条件下，对应求解得到RR 和MR 的稳定流场。接着，在收敛解基础上，改变0.1 个马赫数继续计算，直至激波干扰类型发生转变，进而得到反射类型转变临界Ma∞。本文中，p∞和T∞始终保持和FL-31风洞来流参数一致，Ma∞逐步增大过程为5.7→6.5，Ma∞逐步减小过程为6.7→5.9。研究发现，迟滞转变区间约为Ma∞=5.9～6.5，第2 节将重点对此区间中V 字形钝前缘激波反射双解现象以及转变迟滞过程进行介绍。

2 结果与讨论

针对R/r=1、β=18°的V 字形钝前缘构型，首先介绍增大和减小Ma∞过程中激波反射类型演变过程；然后，分析导致激波反射迟滞现象的根源，并基于流动结构特征建立转变边界；最后，评估迟滞现象对壁面压力特性影响。

2.1 改变Ma∞引起的激波反射迟滞现象

图4 和图5 分别给出了迟滞区间中Ma∞=5.9和6.5 时，V 字形钝前缘根部反射类型对应为RR和MR 的三维流场。图6 给出了Ma∞在增大和减小过程中x-z对称面流场的马赫数云图，并叠加了压力等值线，为了揭示迟滞现象对壁面压力分布影响，对应补充了壁面无量纲p/p0云图。

初始Ma∞=5.9 时，如图4（a）和图4（b）所示，直前缘脱体激波（Detached shock，DS）在V 字形根部倒圆区域相交，发生了同侧激波规则反射RR，在DS 相交点（Intersection point，IP）产生两道透射激波（Transmitted shock，TS）。如图4（b）和图4（c）中流线所示，由于V 字形根部倒圆区对气流的几何约束，来流通过激波偏折并汇聚到倒圆区根部中心位置，气流滞止导致压力显著提升，沿直前缘方向产生较大的逆压梯度，诱导驻点附近气流向两侧偏折，形成大范围的流动分离，分离区内产生大尺度涡和分离激波（Separation shock，SS）。两侧的SS 在驻点上游与TS 发生规则反射后继续相交于点IP′。以Ma∞=5.9 的流场为初场并逐渐增大Ma∞，如图6（a～d）所示，当Ma∞=5.9→6.4 时，激波反射类型一直为RR，DS 的脱体距离随Ma∞增大而减小，TS 长度逐渐缩短，交点IP 向下游流场移动并不断靠近IP′。

图4 Ma∞= 5.9 时RR 反射类型三维流场Fig.4 Three-dimensional flow field of RR at Ma∞= 5.9

当Ma∞=6.5 时，原有的RR 平衡状态遭到破坏，激波干扰结构通过调整后重新达到稳定格局。如图5（a）和图5（b）所示，DS 与马赫杆（Mach stem，MS）相交，发生了马赫反射MR，由三波点（Triple point，TP）发出的TS 以及剪切层Σ，入射到直前缘与倒圆区域相连接的位置。由于TS 前后存在较大的逆压梯度，导致出现小范围的流动分离，并产生SS。SS 与TS 直接相交，形成规则反射。从图5（b）和图5（c）中流线可以看出，包裹在剪切层中的气流由两侧向中心汇聚，并在驻点附近对撞后向上游偏折，使得MS 下游形成大尺度的反转涡对（Counter-rotating vortex pair，CVP），MS 后的流动分为两部分，一部分为对撞导致的逆流，另一部分是直接穿过MS 的气流，两部分气流最后都沿横向溢流。以Ma∞=6.5 的流场为初场并逐渐减小Ma∞，如图6（e～h）所示，当Ma∞=6.5→6.0 时，激波反射结构一直为MR，Ma∞减小直接导致DS脱体距离增大，MS 长度逐渐减小，TP 向下游移动，两侧的分离区范围逐渐增大。当Ma∞=5.9 时，如图6（a）所示，MS 突然消失，激波反射类型由MR转变为RR。

图5 Ma∞= 6.5 时MR 反射类型三维流场Fig.5 Three-dimensional flow field of MR at Ma∞= 6.5

图6 改变Ma∞引起的激波反射迟滞现象Fig.6 Shock reflection hysteresis phenomena by variation in Ma∞

图7 给出了来流Ma∞=6 时FL-31 风洞的实验纹影和数值纹影。由于风洞实验中存在扰动等因素，在实验中也观察到V 字形钝前缘激波反射存在双解现象，相应的两种波系干扰结构均与数值模拟吻合良好。

图7 风洞实验和数值模拟中的激波反射双解现象Fig.7 Dual-solution phenomena in the wind tunnel experiments and numerical simulations

前文的流场分析已经表明，当来流Ma∞=5.9～6.5 时，Ma∞增大或减小过程中，来流条件和几何构型相同，激波反射类型却分别对应为RR 和MR，其双解特性与经典二维理论认识的差异，以及对应的转变机制如何，将在2.2 节进行分析。

2.2 RR-MR 迟滞转变机制

由于x-z对称面流场具有准二维特征，并且能够体现激波结构的主要特征，利用二维理论可以近似分析该平面内的流场参数变化。图8 给出了来流Ma∞=6 条件V 字形钝前缘x-z对称面主激波结构对应的激波极曲线，考虑到流动对称特性，图中只展示了其中一侧分支。从图8 可以看出，来流首先经过激波DS 压缩之后，对应由状态（1）点转变为状态（2）点，当RR 时，气流通过TS 后对应于激波极曲线上的状态（3，4）点，当MR 时，气流通过TS 后对应于激波极曲线上的状态（3′，4′）点。

结合图8 中的反射类型示意图，理论上可以将MR 分为3 类［33］：当状态（3′，4′）的净偏折角为正时，对应为直接马赫反射（Direct-Mach reflection，DiMR）；当状态（3′，4′）的净偏折角为零时，对应为固定马赫反射（Stationary-Mach reflection，StMR）；当状态（3′，4′）的净偏折角为负时，对应为逆马赫反射（Inverse-Mach reflection，InMR），在前人的研究中InMR 被认为是定常激波反射中的反常结构［34］。经典二维激波反射类型转变理论［30］认为双解存在于von Neumann 准则θvN和脱体准则θD之间，双解中的MR 通常为DiMR 或StMR。然而，V 字形钝前缘的双解位于von Neumann 转变边界之下，特别地，如图6（e～h）和图7 所示，V 字形钝前缘的MR 结构中，由三波点TP 发出的滑移线Σ 呈扩张状，状态（3′，4′）点的气流的净偏折角为负，对应出现了InMR 结构。实际上，本文构型中InMR 中的MS 之所以能稳定存在，主要是依靠其下游的反转涡对CVP 维持。

图8 Ma∞=6 时双解区激波极曲线Fig.8 Shock polar diagram for the double solution at Ma∞=6

V 字形钝前缘双解区内的InMR 反射类型以及迟滞转变特性，都与经典二维理论认识存在显著差异，有必要探究其激波反射迟滞现象的转变机制。图9 给出了激波反射类型分别为RR 和MR 时的x-z对称面的流场结构示意图。2.1 节研究表明，流场中主激波结构上的特征位置（如IP、IP′、TP 和CVP）与激波干扰类型的转变密切相关，对这些特征位置的定量分析，有助于理解迟滞转变机制。为了便于下文对这些特征位置的分析，将特征点到驻点S 的水平距离定义为参数d。图10 给出了特征点位置参数d随Ma∞变化，其中，散点符号表示数值结果数据，曲线为理论分析或由散点通过二次多项式拟合得到的结果。

图9 不同激波反射类型示意图Fig.9 Sketches of different shock reflection structures

图10 迟滞过程中不同反射类型特征点位置变化Fig.10 Positions of characteristic points of different shock reflections during the hysteresis process

当激波反射类型为RR 时，如图9（a）所示，两侧DS 的交于点IP，两侧SS 与TS 分别干扰后进一步相交于点IP′。如图10 所示，当增大Ma∞时，DS脱体距离δ逐渐减小，导致点IP 向下游移动，即点IP 到驻点S 的距离dIP，S逐渐减小；点IP′的位置主要受大分离区范围以及SS 影响，数值模拟结果显示大分离区范围以及SS 随Ma∞变化并不明显，如式（1）和图10 所示，可以将点IP′到驻点S 的距离dIP′，S视为定值。当点IP 和点IP′重合时，即满足dIP，S=dIP′，S，波系结构将不再稳定，将会引起激波干扰类型从RR 向MR 转变。因此，确定点IP 的位置，是准确预测RR 向MR 转变的关键。

由几何关系可知，点IP 到驻点S 的距离dIP，S主要与DS 的脱体距离δ有关，dIP，S和δ满足式（2）关系。为了求解δ，可以将斜掠直前缘看作一段后掠圆柱，DS 将沿其流向逐步发展，当直前缘长度L足够长时，DS 最终达到充分发展状态，即激波的脱体高度δ到达定值δf。Zhang 等［32］发现δ主要与垂直于后掠圆柱方向的气流分量Ma∞sinβ有关，如式（3）所示，根据无黏圆柱脱体激波高度近似理论［35］，可以求解得到充分发展状态的脱体距离δf。由于本文采用的构型，直前缘长度限制，δ未达到充分发展状态，但基于δf变化规律，根据数值模拟结果进行修正，可以近似获得δ，如式（4）和图10 所示。联立式（1～4），进而可获得式（5）表示的转变边界fRR→MR（Ma∞）。在图10 中，dIP，S和dIP′，S变化曲线在Ma∞= 6.48 时相交，从理论上确定了临界转变边界fRR→MR位置，而数值计算得到的fRR→MR位置约为Ma∞= 6.5，两者偏差仅约0.3%。

当激波反射类型为MR 时，如图9（b）所示，DS、TS 和MS 交于TP，此时，汇聚到V 字形根部气流对撞后向上游偏折，进而产生CVP。当减小Ma∞时，δ逐渐增大，TP 逐渐向下游移动，同时MS长度减小。由图6 可知，MS 长度减小至一定值后突然消失，反射类型由InMR 转变至RR。实际上，InMR 中MS 产生或消失都与CVP 存在很强的关联。本文采用对称面上逆流最高位置（鞍点，见图5（c）至驻点S 的距离dCVP，S表征CVP 尺度变化，采用TP 至驻点S 的距离dTP，S表征MS 变化，如图10所示，dTP，S和dCVP，S变化趋势基本一致，这进一步印证了MS 主要受CVP 尺度变化影响。

结合图10 中dCVP，S、dTP，S和dIP，S变化过程，不难推断出，InMR 转变至RR 的充分条件为，CVP 尺度减小至不再影响DS 直接相交，即dCVP，S=dIP，S，此时MS 将消失，反射类型转变。InMR 中dCVP，S和dTP，S均难以理论获得，从数值模拟结果中，可以拟合获得dCVP，S和dTP，S随Ma∞的变化规律，如式（6）和图10 所示，其中，dTP，S可以近似看作dCVP，S与CVP 导致的MS 脱体距离之和。联立式（2～4）和式（6），进而可获得式（7）表示的转变边界fMR→RR（Ma∞）。图10 中，dCVP，S和dIP，S变化曲线在Ma∞=5.85 时相交，从理论上确定了转变边界fMR→RR位置；而数值计算得到的fMR→RR位置约为Ma∞=5.9，两者偏差仅约0.8%。

以上分析表明，迟滞现象将V 字形钝前缘激波干扰问题变得更加复杂，反射类型双解的出现不仅使流场结构、流场参数明显改变，而且不同反射类型对应产生的透射激波、剪切层和超声速射流等与壁面相互作用，将会导致其气动载荷特性出现极大的差异。

2.3 迟滞现象对壁面压力特性影响

图11 给出了不同Ma∞条件下RR 时V 字形根部倒圆段壁面中心线上p/p0分布和局部压力峰值附近流场，结合图6（a～d）可以看出，复杂激波干扰导致其壁面压力分布极不均匀。当φ=0°～12°时，流场上、下两侧的脱体激波DS 之间，以及分离激波SS 之间的相互干扰，使得经过多道激波压缩后的气流从两侧向流场中心附近汇集，超声速气流直接冲击驻点，形成滞止激波BS，导致驻点位置处出现极其严酷的中心压力峰值（Central peak）。当Ma∞= 6.4 时，BS 形态受到上游波系影响发生变化，导致中心压力峰值分布一定程度改变。当φ=12°～72°时，由于处于大分离区位置，壁面p/p0相对较低。

图11 RR 反射类型壁面压力分布Fig.11 Surface pressure along the centerline of RR

图12 给出了不同Ma∞条件下MR 时壁面中心线上p/p0分布和局部压力峰值附近流场。结合图6（e～h）可以看出，由TP 发出的TS 和Σ 入射倒圆段的壁面，TS 与近壁面附近的SS 发生规则反射，产生透射激波，在下游气流再附，产生再附激波。因此，从倒圆区域的起始位置处（φ=72°）开始，p/p0开始迅速升高，在φ=36°～42°位置的p/p0出现最外侧峰值（Outermost peak）。进一步地，包裹在剪切层中的气流由两侧向中心汇聚，在驻点附近发生对撞，在φ=0°位置的p/p0出现中心峰值（Central peak）。

图12 MR 反射类型壁面压力分布Fig.12 Surface pressure along the centerline of MR

图13 给出了不同Ma∞条件下RR 和MR 时的壁面不同位置压力峰值变化。当Ma∞增大或减小时，中线峰值在Ma∞=5.9～6.5 范围出现迟滞环。当RR 时，中心峰值随着Ma∞升高而逐渐增大，在Ma∞=6.4 时有所减小；当MR 时，中心峰值和最外侧峰值都随着Ma∞降低而逐渐减小。由于不同激波反射类型对应的压力峰值产生机制存在差异，RR 工况对应中心压力峰值约为MR 工况对应中心压力峰值的2～3 倍。此外，根据前期将壁面条件设置为等温壁的研究结果［13，16-17］，RR 工况对应的最大热流值也远高于MR 工况对应的最大热流值，更容易造成烧蚀破坏。