锯齿形转捩片触发高超声速进气道边界层转捩的大涡模拟

2019-10-31张红军朱志斌尚庆刘智勇沈清

航空学报 2019年10期

张红军，朱志斌，尚庆，刘智勇，沈清

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

边界层转捩对高超声速飞行器气动特性和热防护设计具有重要影响。对于吸气式发动机而言，若进气道进口处边界层为湍流，则可有效消除激波/边界层干扰导致的局部分离，增加质量捕获，减小层流边界层分离造成进气道不起动的风险。飞行试验发现，在来流低湍流度、强激波压缩、钝前缘熵层等因素作用下，高超声速进气道压缩面边界层通常保持为层流。因此，必须在前体压缩面安装人工转捩装置来获得湍流，以确保进气道的正常工作。

研究发现，钻石形或后掠斜坡形的转捩装置能够在边界层内诱导产生一系列沿流向运动的涡对，进而快速触发高超声速边界层转捩。这些转捩装置因此被称为涡流发生器，并已被广泛应用，如X-43A[1-3]、X-51A[4]和Hyfly[5]、HIFiRE[6]等高超飞行器中均采用了这种转捩装置。在涡流发生器转捩机理方面，美国兰利研究中心的Choudhari等[7-8]针对Hyper-X缩比模型前体上由斜坡型涡流发生器引起的扰动流场，通过求解二维特征值的方法对流场进行了稳定性分析，分析了强制转捩装置后3个压缩面上的扰动波模态，并采用eN方法对扰动增长率进行积分来获得转捩位置。Iyer和Mahesh[9]对平板上安装的半球形粗糙元进行了研究，将粗糙元诱发边界层转捩的机理总结为：由于粗糙元的存在，边界层内会产生三维的流动分离，进而会形成剪切层和涡结构；在粗糙元下游，不稳定的反转涡结构会不断冲击剪切层从而形成发卡涡。

国内在转捩装置的机理研究方面也开展了大量工作。中国空气动力研究与发展中心的赵慧勇等[10-11]在风洞中开展了钻石形和斜坡形涡流发生器的强制转捩试验，研究了钻石形涡流发生器触发转捩的有效高度，并应用大涡模拟方法分析了其转捩机理[12]，给出了强制转捩过程的流场结构和流向涡的失稳模式。清华大学的肖志祥等[13-16]对高超声速边界层内不同形状粗糙单元导致的强制转捩现象进行了直接数值模拟，研究发现此类强制转捩主要由粗糙元顶部的三维剪切层失稳导致，并对多种粗糙元的转捩效果进行了定量研究。朱德华等[17]通过数值模拟,从边界层稳定性和拓扑结构稳定性角度分析了钻石型粗糙元的主要转捩机理。涂国华等[18]研究发现在超声速边界层中布置悬空的细丝可促进边界层失稳,他们还针对高超声速强制湍流提出了一种湍流模型修正方法。战培国[19]归纳总结了各风洞在超燃冲压发动机前体边界层强制转捩试验中采用的主要测量和显示技术,介绍了美国开展Hyper-X前体边界层强制转捩研究风洞设备的选择依据和选用的主要风洞,分析了强制转捩装置设计过程中风洞试验研究采用的方法。

与涡流发生器不同，作者所在研究团队提出了一种呈锯齿状的转捩薄片(见图1(a)和图1(b))。由于其厚度小(约0.2 mm)，对流场的干扰很小，因此不会对进气道性能带来附加损失，具有低阻、低热流、应用方便等优点。2015年作者团队在FD-07风洞中开展了锯齿形转捩薄片在高超声速二元进气道上应用的风洞试验，试验来流马赫数Ma=5，6。由进气道波系纹影(图1(c)和图1(d))可以明显看出在不粘贴转捩片时进气道入口前存在明显的分离激波，表明进气道不起动，在粘贴转捩片后进气道入口前的分离激波消失，表明进气道顺利起动。两个马赫数的试验情况类似，文献[20]给出了试验的详细情况。由于进气道试验未能获得有效的边界层转捩信息，因此为促进锯齿形转捩片的进一步工程应用，本文针对二元进气道上的锯齿形转捩片触发高超声速边界层转捩流动开展了大涡模拟研究，以认识转捩形态特征，并揭示其触发转捩机理。

图1 锯齿形转捩薄片及进气道波系纹影[20]

1 数值方法

基于三维Favre滤波Navier-Stokes方程[21]，采用大涡模拟方法对大尺度湍流结构直接求解，通过构造亚格子模型或者依靠数值耗散作用来模化亚格子尺度流动对可解尺度的影响。本文基于亚格子模型对可解尺度流动起耗散性作用的假设，采用隐式大涡模拟方法[22](Implicit Large Eddy Simulation, ILES)，依靠数值格式的耗散特征来提供湍流动能耗散，无需添加显式亚格子模型。

对流通量离散采用特征通量限制型紧致格式[23]，该格式具有较高的精度和分辨率，并能够在精细分辨流场小尺度结构的同时光滑捕捉激波等间断。黏性通量采用二阶中心格式计算，非定常时间推进采用二阶显式Runge-Kutta方法。计算时分别以来流密度、温度和速度作为数值计算的无量纲参考量，参考长度取1 m。文献[24]对该方法的可靠性进行了验证，数值计算结果与试验结果吻合良好。

2 模型与网格

2.1 计算模型及工况

计算模型包括三楔压缩面、等熵压缩面及单楔压缩面外形，图2给出了3种外形的尺寸及转捩片构型。其中压缩面前缘半径R=0.2 mm，锯齿前缘到压缩面前缘距离L=40 mm，转捩片厚度T=0.2 mm，齿高C=4 mm，齿间夹角θ=90°，后缘长度D=2 mm。计算来流条件与开展的地面风洞试验条件一致，如表1所示。

图2 计算模型及尺寸

表1 计算条件

2.2 计算网格及边界条件

为减小计算量，取3个转捩片单元进行数值模拟分析，计算网格如图3所示。计算域为流向x∈[20,424] mm，展向z∈[0,12] mm，法向y∈[0,12] mm。网格单元总量为55 216 800，在壁面和转捩片处对网格进行加密，最小网格尺度为0.01 mm，展向网格尺度为0.13 mm，展向在第2楔面湍流区流向网格尺度为0.21 mm，对应无量纲计算网格尺度满足：

(1)

式中：ξ、η和ζ分别对应流向、法向和展向，上标“+”表示以壁面处黏性系数和壁面摩擦速度定义的归一化尺度。因此，计算网格分辨率能够满足大涡模拟的网格尺度要求[25]。为便于对比分析，无转捩片时的计算网格点密度与转捩片计算网格一致，展向范围为转捩片计算域的一半。

计算域入口剖面由相同位置处光滑壁面的二维层流流场截取得到，壁面设为无滑移等温条件，法向和流向出口外插，展向取周期条件。

3 结果分析

3.1 转捩形态特征

根据有/无转捩片情况下3级楔压缩面的大涡模拟结果分析转捩片诱导边界层转捩的形态特征。

图4给出了有/无转捩片时进气道展向对称面的瞬时流场(图中ρ、T分别为密度和温度，下标∞表示来流参数)。可以发现，无转捩片时边界层始终保持为稳定的层流；有转捩片时，流场后部区域(300 mm

图5给出了有转捩片时流场的瞬时速度梯度第二不变量等值面(u为x方向的速度，U∞为来流速度)，表示了流场瞬时涡系结构的空间发展形态。可以看到，在转捩片后方首先出现了非常微弱的展向涡结构，经过第2级压缩拐角在第3级压缩面上出现了与转捩片展向波谷位置相对应的三维流向涡，流向涡破碎后诱发大量马蹄涡结构，使流动快速进入湍流。

图4 有/无转捩片时对称面的瞬时流场

图5 瞬时涡系结构的空间发展形态

图6进一步给出了不同流向位置处的瞬时无量纲流向涡量(ωx)云图，显示了流向涡的空间发展演变过程：在第1级压缩拐角前没有涡结构出现，经过第1级压缩拐角后，在x=200 mm处出现了微弱并且规则的反向旋转涡对结构，随着流动向下游发展，在第2级压缩拐角前(x=225 mm)，涡对结构已明显增强；经过第2级压缩拐角后，涡对结构发生扭曲和变形，直至完全破碎，最终发展为混乱无规则的大尺度流动结构。

图6 不同位置的瞬时流向涡量云图

图7和图8分别为有转捩片情况下瞬时摩阻系数Cf、热流Q和壁面摩擦力线分布。可以看到，瞬时摩阻、热流分布与瞬时涡系结构相对应，在转捩前呈现明显的条带分布特征，在湍流区域分布变得无规则。此外，从摩擦力线分布可以看到，摩擦力线在分离点(x=150 mm)前保持平直，表明流动为层流状态；在分离和转捩区域内(150 mm400 mm)，摩擦力线变得无规则。

有/无转捩片时的平均摩阻系数和热流曲线对比如图10所示。摩阻曲线小于零的部分代表分离区的范围，可以看到，有转捩片时分离点位置后移，而再附点位置前移，分离区流向长度从150 mm减小为85 mm；另一方面，有转捩片时摩阻和热流在流动再附后急剧升高，经过转捩峰值后平缓变化。

图7 瞬时摩阻系数及摩擦力线

图8 瞬时热流及摩擦力线

图9 脉动动能分布

图10 平均摩阻系数和热流曲线对比

3.2 转捩机理分析

3.1节的大涡模拟结果表明，锯齿形转捩片能够触发三级楔压缩面边界层转捩，从转捩形态来看，在流向依次出现了涡对结构、条带结构和马蹄涡结构，分别对应着流动失稳的不同阶段。从流场来看，在逆压梯度作用下，两级压缩拐角处为扁平状的分离区，分离区改变了压缩拐角的流动结构，使得压缩拐角呈现出连续的凹面流动特征(如图11所示，图中ρrms为密度的均方根值)。如果在凹面边界层中施加Görtler扰动，由于离心力和壁面法向压力梯度之间的不平衡，则在流场中会出现Görtler涡， Görtler涡会进一步发展出低速高速条带，而条带的二次失稳通常被认为是导致转捩的关键因素[26-29]，这表明凹面边界层的转捩机制与平面流动存在明显差异。锯齿形转捩片产生的三维扰动虽然不是Görtler扰动，但涡对结构的发生发展过程与Görtler涡类似。

图11 分离区流动结构

与一般凹面流动不同的是，本文所研究的压缩拐角包含有流动分离，这使得对转捩机理的分析变得更加复杂。为进一步深入认识流动的曲率效应和由流动分离所产生的自由剪切层对边界层转捩的作用，参考三级楔流动中分离区的分离点位置重新设计了一等熵压缩面(图2(c))，以消除流动分离。另外，将原三级楔压缩面的后面两个压缩面去掉并将第一级延伸至相同的流向长度(图2(d))，以考察曲率效应。采用大涡模拟方法对带有相同转捩片构型的等熵压缩面、单级楔面一并进行了研究。

图12给出了带有锯齿形转捩片的单级楔面的大涡模拟结果。可以看出，对于单级楔来说，壁面边界层始终保持为层流，即在计算的流向长度范围内转捩片产生的三维扰动在没有逆压梯度的情况下未能使边界层转捩，说明逆压梯度与三维扰动同时构成了流动失稳直至转捩的必要条件，这与完全依靠涡流发生器[12]触发边界层转捩的机制明显不同。

图13、图14给出了带有锯齿形转捩片的等熵压缩面的大涡模拟结果。可以看出，对于无分离的等熵压缩面，同样出现了类似于Görtler涡的流向涡对，其发展直至转捩的过程与三楔压缩面类似，表明涡对结构是由凹面的曲率效应所产生的，流向涡对发生扭曲和变形、直至完全破碎，最终发生转捩。从图15的热流分布曲线来看，等熵压缩面在流向x=310 mm处流动完成转捩，较三级楔的情况(x=300 mm)略靠后，说明与内凹壁面相比，内凹的剪切层起到了加速流动失稳的作用，进而使转捩位置提前。

图12 单楔瞬时流场(对称面)

为验证分离区剪切层对流动失稳的放大作用，在分离区前后各选取一个位置对三楔压缩面和等熵压缩面开展线性稳定性分析。

图16给出了有/无转捩片时分离区前(x=60 mm)和分离区后(x=280 mm)两个位置处的边界层速度型。对比可发现，在x=60 mm处，转捩片的存在未对边界层速度型产生显著影响，只在靠近物面附近带来微小的变化；在x=280 mm处，有转捩片时三楔压缩面和等熵压缩面的速度剖面均存在明显的拐点，而无贴片的速度剖面不存在拐点，这直接决定了有无转捩片时流动稳定性的差别。

图13 不同位置的瞬时流向涡量云图(等熵压缩面)

图14 等熵压缩面瞬时流场

图15 等熵压缩面与三楔压缩面的热流对比

图16 边界层速度型

采用线性稳定性分析方法对以上位置的流动稳定性进行了分析，如图17所示(其中ω表示扰动波的频率，αi表示扰动波的增长率，负值表示扰动波是不稳定的)。对比发现，在x=60 mm处，有无转捩片时αi量值差别很小。而在x=280 mm处，有无转捩片时稳定性特征差异明显：无转捩片时的流动具有非常弱的第二模态不稳定性；有转捩片时，等熵压缩面的不稳定波频率范围较无转捩片时明显变宽，不稳定波最大增长率是无转捩片时的1.6倍左右；而对于三楔压缩面来说，其不稳定波频率范围则完全涵盖了前两者，不稳定波最大增长率则是等熵时的2.5倍左右，这说明分离区剪切层对改变失稳波性质、加速流动失稳起到了显著的作用。