高超声速边界层基频二次失稳条纹结构的稳定性

李玲玉，刘建新

(天津大学 机械工程学院 高速空气动力学研究室，天津 300072)

0 引言

边界层从层流向湍流的转变过程被称为边界层的转捩。转捩问题实际是湍流的起源问题，是流体力学长期关注但尚未解决的重要研究领域之一[1]。研究清楚转捩过程的机理，对于提高转捩预测的准确性和理解湍流的产生至关重要。一般认为边界层的转捩是由流动不稳定性引起的，分析流动的失稳特性和研究扰动演化的过程，是流动稳定性分析的主要工作[2]，也是理论以及工程上所关心的重要问题之一。

在高超声速边界层中，由于流动的压缩性，边界层局部会产生相对超声速区，由此产生的失稳模态被称作第二模态，它在较高来流马赫数时对流动失稳起着至关重要的作用。对于马赫数较高的情况（如Ma>4），二维的第二模态的线性增长率最大，最不稳定，可能最先增长起来[3]。但是只有二维流动是不能导致转捩的，因为湍流是三维流动的，这就说明中间过程会有三维流动的产生。Herbert的二次失稳理论是在这个过程中研究三维扰动产生原因的理论之一[3]。二次失稳的类型主要包括基频模态失稳和亚谐模态失稳[4-7]，其中前者通常被称作为K型失稳而后者则被称作为H型失稳，表现为在转捩中不同的流场分布特征。早期的不可压缩边界层的相关研究通常认为K型失稳发生在较高的来流湍流度的情况下。

直到20世纪90年代之后，人们才将二次失稳的研究对象从低速问题转移到了高超声速问题。在1992年Ng 和 Erlebacher[8]采用二次失稳理论研究了Ma=4.5时的平板边界层中第二模态的二次失稳，发现了亚谐模态的二次失稳解，但没有发现基频模态的二次失稳解。近年来，在高超声速边界层的直接数值模拟和实验研究中，人们相继发现在Mack模态主导的边界层转捩前出现了条纹结构。2010年刘建新[9]在小迎角钝锥高超声速边界层的扰动演化研究中发现：入口的等幅值的扰动波会在下游形成展向波包型分布，在更下游的流场中会出现三维扰动，扰动的定常部分为流向条纹结构，如图1所示。还发现三维扰动快速增长的机理应该属于二次失稳理论中的基频模态失稳。随后Purdue大学的Schneider小组(2011)[10]在静风洞实验中，通过热敏漆（TSP）技术，观测到了压缩锥壁面热流呈现的三维流向条纹结构，如图2所示。此后不久Yu 和 Luo（2013）[11]在平板边界层中也发现了条纹结构，研究表明这种结构是由与第二模态同频率的三维扰动快速增长引起的，类似于不可压流动中的基频模态失稳情况。Sivasubramanian 和Fasel （2016）[12]的DNS研究也表明实验中出现的条纹结构是由基频模态失稳引起的，条纹结构对应于基频共振作用过程中激发的流向涡结构。后来Christoph和Fasel （2019）[13]用直接数值模拟方法研究了马赫数为6、迎角为0°时，高超声速裙锥边界层由层流到湍流的转变过程，结果发现基频共振要强于亚谐共振。基频共振诱发了条纹结构的产生，在转捩中占主导地位。刘建新等[14]在最近研究了不同Mack模态幅值对二次失稳类型的影响。当Mack模态幅值较小时，二次失稳以亚谐模态失稳占主导；当Mack模态幅值较大时，基频模态的增长率将超过亚谐模态，成为二次失稳的主要模态。这些研究都表明，在由Mack模态主导的边界层转捩过程中，可以通过二次失稳机制激发出由基频模态主导的条纹结构。

图1 刘建新数值结果中的条纹结构[9]Fig. 1 Streaky structures in numerical results by Liu[9]

图2 Purdue大学静风洞中观测到的条纹结构[10]Fig. 2 Streaky structures observed in the quiet wind tunnel of Purdue University[10]

很多研究都发现流场在变成湍流过程中都存在比较稳定的条纹结构。从条纹结构的生成来看，很多机制都会引起或者表现为流场内条纹结构的产生。除了前面提到的边界层内由Mack模态引起的基频二次失稳可以产生条纹结构以外，最典型的，条纹结构的生成可以和边界层内失稳扰动的非线性演化至饱和过程有关，例如横流涡、哥特涡等在流场中表现为条纹结构[15-17]；三维边界层中由于边界层抬升效应也会产生典型的条纹结构，例如HiFiRE-5模型（椭圆锥）短轴处及粗糙元尾迹流动中出现的流向涡结构，可表现为流向条纹[18]；此外，当外界湍流度较大时，也可以依靠感受性机制直接在边界层内激发出条纹结构[19-20]。

针对于以上如哥特涡、横流涡等机制引发的条纹结构，人们开展了相关的二次失稳研究，并将这些二次失稳机制与转捩过程相关联。研究结果说明湍流的突变（breakdown）过程可能和条纹结构的二次失稳有关，这说明条纹结构对转捩有很大的影响，对转捩的发生起着至关重要的作用。虽然针对于Mack模态的二次失稳机制引起边界层内条纹结构的产生问题开展了一些研究工作，但对由此产生的条纹结构的二次失稳研究则关注不够。基于此，本文以Ma=4.5的可压缩平板边界层为研究对象，以Blasius相似性解为二维基本流，采用线性稳定性分析和二次稳定性分析等研究方法开展研究，旨在研究平板边界层失稳过程中基频模态产生的条纹结构的分布特征与规律，探究这种条纹结构的稳定性特征。对于此类条纹结构的稳定性研究，有助于我们更好地理解和描述条纹结构在转捩中的作用，有助于推进高超声速边界层由Mack模态引起的边界层转捩机理的研究。

1 方程及求解所用的数值方法

本节介绍本文中用到的方程及求解所用的数值方法。

1.1 线性扰动方程

设有可压缩气体流过平板，流动方向为x方向，y方向与z方向分别为法向和展向。考虑可压缩边界层，从有量纲的完全Navier-Stokes （N-S）方程出发，选取适当的特征量将其无量纲化，我们可以得到无量纲的N-S方程。将得到的无量纲N-S方程中的瞬时量设为基本流与扰动量的和：

其中，φ =(ρ,u,v,w,T)T为 瞬时量，上标 T表 示转置。φ0=(ρ0,u0,v0,w0,T0)T为 定常的基本流，φ′=(ρ′,u′,v′,w′,T′)T为扰动量,。将式（1）代入到无量纲N-S方程中，并减去定常基本流满足的N-S方程，最后略去非线性部分，可以得到线性化扰动方程：

其中，Vxx、Vyy、Vzz、Vxy、Vxz、Vyz、Γ 、A、B、C和D为系数矩阵，系数矩阵元素的具体表达式见周恒等[21]的著作。方程（2）描述了流动中小扰动的线性演化过程，对该方程求解需要给出初始条件和计算域所有边界上的条件，而且计算量很大，因此需要进行化简来降低求解的难度，下一小节介绍对该方程的简化方法。

1.2 线性稳定性方程

首先对二维基本流进行线性稳定性分析（LST）。考虑到扰动的演化相对基本流主要发生在流向、展向和法向，可以假设扰动是时间、流向和展向的周期函数，则扰动可以写成下列形式：

将式（3）代入到线性的扰动方程（2）中，并忽略基本流的流向导数和法向速度，得到准平行性假设下的边界层的线性稳定性方程：

其中，LLST是 方程的线性算子，是以α、 β 和 ω为变量的函数。可以设扰动在壁面和远场扰动为零。这样线性稳定性问题就可以被认为是一个具有齐次边界条件的特征值问题，要想有非零解，需要满足：

式（5）描述了扰动的色散信息，即波数与频率的关系。通过求解式（5）得到色散关系后，可以从式（4）中求得形函数。所以说稳定性方程的实质正是通过求解这样一个特征值问题，得到边界层扰动波的色散关系和扰动的特征函数。该方程可以使用五点四阶中心差分格式在法方向上进行离散，矩阵的特征值问题可以通过反幂法或者QZ法进行迭代求解或者直接求解。

1.3 二次失稳分析方法

经典的线性稳定性方程主要是针对二维基本流开展的稳定性研究，它所给出的失稳扰动信息又被称作首次失稳。然而，当二维基本流与首次失稳的扰动叠加时，又可以产生新的失稳现象，即二次失稳。本文采用经典的基于Floquet理论的二次失稳理论（SIT）进行二次失稳分析。根据二次失稳理论[8]可知，二次失稳分析的基本流由原基本流和有限幅值的首次失稳二维扰动组成，即：

其 中， φ0对 应于 原 基 本 流，对应于首次失稳扰动的形函数，A为首次失稳扰动的幅值。α和 ω对应于首次失稳扰动的流向波数和频率。二次失稳分析的基本流是时间t的函数，为了消除t，使得二次失稳分析的基本流定常，这里引入坐标变换：

通过这样的变换，在首次失稳扰动的行进坐标系中，新构成的基本流就可以被看作为定常的。

为了研究二次失稳，在新构成的基本流即二次失稳的基本流上叠加二次失稳扰动，则瞬时流场可表示为：

将式（10）代入到式（9）中，可以得到方程：

下面给出边界条件，在壁面y=0处：

在远场y=∞处扰动趋于零：

对于一个给定的展向波数 β2，未知的复特征值σ描述了二次失稳扰动的时间失稳。二次失稳分析既适用于时间模式，也适用于空间模式。由于本文研究的是时间模式下的情况，所以只给出了时间问题解的推导，空间问题解在这里不再赘述。类似的，该方程也可以稍加改动即可适用于求解条纹边界层的无黏失稳情况。

式（11）是一个大规模的特征值问题，我们在法向上采用四阶中心差分格式，流向采用傅里叶谱方法进行离散。特征值问题选用Arnoldi子空间迭代法进行求解，该方法是求解大规模稀疏矩阵局部特征值的迭代方法。ARPACK[22]是目前最流行的实现Arnoldi方法求解大型稀疏矩阵特征值的开源软件包。

1.4 程序验证

为了验证本文所用程序的正确性，计算了Ng 和Erlebacher 等[8]Ma=4.5平板边界层二次失稳分析的算例。首先用线性稳定性方程求解器计算了Re=10000、α=2.52、T=61.11K的工况，得到了Mack模态扰动。图3和图4分别给出了Mack模态的速度和 温度随法向坐标的变化关系。

图3 Mack模态的流向速度沿法向的分布Fig. 3 Distribution of streamwise velocity in Mack mode along the normal direction

图4 Mack模态的温度沿法向的分布Fig. 4 Distribution of Temperature in Mack mode along the normal direction

从图中不难看出，Mack模态扰动的温度最大值明显大于速度最大值，因此以Mack模态的温度最大值来定义Mack模态扰动的幅值，即这样定义可以避免负温度的出现。为了和文献[8]作比较，本文的Mack模态扰动的幅值分别取为6%、3%和2%。

在二次失稳分析中，不考虑Mack模态的非线性作用，利用上述稳定性分析得到的Mack模态扰动乘以相应的幅值加上相似性解作为二次失稳分析的基本流，采用时间模式计算了三组不同幅值的工况。来流参数为：Ma=4.5， 来流温度为T=61.11K，边界层位移厚度定义的雷诺数为Re=10000。首次失稳扰动为二维的Mack模态扰动，无量纲流向波数为α=2.52 ，无量纲频率为 ω=2.27。最后将计算结果与文献[8]结果进行对比，图5为本文计算结果和文献结果的比较,从图中可以看出本文的计算结果与文献结果吻合的很好，从而验证了本文所用的分析程序是可靠的。

图5 二次失稳分析的结果与文献[8]的结果对比图（实线为文献中结果，方块为本文结果）Fig. 5 Comparison between the secondary instability analysis results and those in Ref.[8](The solid lines are the results in literature,and the square symbols are the results in this study)

2 条纹结构的稳定性分析

2.1 首次失稳扰动幅值对二次失稳的影响

此前研究中都发现，首次失稳的扰动幅值对二次失稳的主导模态会产生较大的影响，为了研究这种影响，我们利用1.4节中通过线性稳定性方程求解器得到的Mack模态扰动作为首次失稳扰动，首次失稳扰动乘上相应的幅值A加上相似性解，作为二次失稳分析的基本流。研究了首次失稳扰动的幅值A从0.06增加到0.8时，亚谐模态和基频模态的增长率随展向波数 β2的变化关系，结果如图6和图7所示。从两张图中可以看出，随着Mack模态幅值的增加，亚谐模态和基频模态的失稳波数范围不仅扩大，其相应的增长率也在增大。亚谐模态的最大增长率从0.025增加到了0.064，基频模态的最大增长率从0.0008增加到了0.065。

图6 亚谐模态的增长率随展向波数变化的关系Fig. 6 Relationship between the growth rate of subharmonic modes and the spanwise wavenumber

图7 基频模态的增长率随展向波数变化的关系Fig. 7 Relationship between the growth rate of fundamental modes and the spanwise wavenumber

结合两张图不难看出，相对于亚谐模态来讲，基频模态的增长率随着初始Mack模态幅值的增大而增长的更快，说明首次失稳扰动的幅值对基频模态扰动的影响更显著。进而说明在首次失稳扰动小幅值的情况下，亚谐模态在流场中占主导地位；反之基频模态占主导地位，这个结论与刘建新[14]等研究的结论是相符合的。这也解释了为什么在Ng和 Erlebacher等[8]的论文中，只发现了亚谐模态的二次失稳，没有发现基频模态的二次失稳。因为文献中的首次失稳扰动幅值最大为6%，幅值等于6%属于小幅值的情况，此时流场中亚谐模态占主导地位，所以没有发现基频模态的二次失稳。

由于本文重点研究基频模态的失稳，所以围绕亚谐模态失稳的研究这里不再展开。图8给出了首次失稳扰动幅值A=0.6 ， 展向波数 β2=2.1时，基频模态的特征函数。可以看出，基频模态速度的最大值集中在y≈0.6附近。

图8 基频模态的流向速度特征函数沿法向的分布Fig. 8 Streamwise velocity eigenfunction of the fundamental mode along the normal direction

进一步的，我们可以对基频模态的特征函数进行傅里叶分析，即分析特征函数分量的值的分布特征。图9给出了最大的三个分量的分布。

图9 基频模态特征函数的傅里叶分析Fig. 9 Fourier analysis of the fundamental mode eigenfunction

2.2 二次失稳扰动幅值对条纹结构基本流的影响

通过对前人研究工作的调研，我们知道了在首次失稳扰动大幅值的情况下，基频模态在流场中占主导地位，而且随着流动向下游的发展，会产生条纹结构。为了研究条纹结构的失稳特性，我们首先要获得条纹结构。这里选择Blasius相似性解为原基本流，Mack模态作为首次失稳扰动，其幅值为0.6，基频模态扰动作为二次失稳扰动，根据一般计算的经典二次失稳结果，进行后处理，从而生成全部的流场和条纹结构的流场。这样，考虑到基频失稳扰动 σ的虚部为零，则二次失稳扰动增长起来后的流场可写为：

2.3 二次失稳扰动幅值对全部流场的影响

根据流场的性质，人们在研究中总是通过分析流动过程中涡结构的演化特点来研究扰动的发生、演化及其表现特征的。

图10给出了首次失稳扰动幅值A=0.6，展向波数 β2=2.1时，不同二次失稳扰动幅值下，全流场的流向涡量云图。可以看出流向涡量是呈对称分布的，且分布比较均匀。随着二次失稳扰动幅值B的增加，流向涡量的强度也越来越大。

图11给出了不同二次失稳扰动幅值下，全部流场的展向涡量云图，不难看出随着扰动的发展，y≈1附近的流体逐渐被卷起，形成涡结构。随着二次失稳扰动幅值B的增加，展向涡量的强度越强。

对比图10和图11两组相同二次失稳扰动幅值下的流向涡量和展向涡量云图可以发现，流向涡量的强度要显著大于展向涡量的强度，说明在流场中流向涡量占主导地位，展向涡量次之。这意味着此时流场中的动力学结构应该主要受流向涡量的制约。因此，以下研究中我们主要关心流向涡量的条纹结构，而忽略掉展向结构的影响和特征，只研究某一个流向站位的流场结构。

图10 流向涡量云图Fig. 10 Streamwise vorticity contours

图11 展向涡量云图Fig. 11 Spanwise vorticity contours

根据前面的讨论，所谓流向涡量的条纹，其主要成分为基频模态的零阶扰动，因此，此时流向涡结构即基频失稳条纹可以写成如下形式：

为了研究不同二次失稳扰动幅值对条纹结构基本流的影响，我们计算了四组不同的参数，二次失稳扰动幅值分别为0、0.2、1和2。

图12给出了两个周期内，不同二次失稳扰动幅值下的条纹结构流向速度云图。从图中不难看出，随着B的增加，条纹结构的峰值越来越大，条纹结构越来越明显。说明B对条纹结构的产生有促进作用，二次失稳扰动的幅值越大，那么对应的条纹结构的峰值也越大。

图12 不同B值下的条纹结构（流向速度）沿法向的分布Fig. 12 Streamwise velocity contours of streaky structures under different B values along the normal direction

2.4 对条纹结构的无黏稳定性分析

得到条纹结构后，我们就可以进一步对该条文结构进行稳定性分析。这里我们使用无黏稳定性方程进行分析。因为根据通常的理论，条纹结构的稳定性主要是依赖于流场的局部剪切造成的，而黏性只起到阻尼作用，因此直接使用无黏稳定性分析就可以给出稳定性分布的本质与特征。无黏的控制方程可以写为：

其中，α为流向波数，Mˆ 为相对马赫数，具体表达式可见式（17）。对应于y-z平面的拉普拉斯算子，对应压力扰动的特征函数。

为了验证我们的分析，本文计算了两个二次失稳扰动幅值为0.2和2时的线性条带流场对应的无黏稳定性的工况。图13给出了二次失稳扰动幅值为0.2时最不稳定的两支扰动的增长率随波数的变化关系，对应图中的F1和F2曲线。两支解都是随着波数的增加，增长率呈先增大后减小的趋势。从图中不难看到，两支模态的增长率分布特征与二维Blassius边界层时的分布状态十分相似。这是由于此时条纹幅值较小，流场中保留了更多的二维性的特征。图14给出了与图13相对应的两支曲线的相速度随波数的对应关系。从图中不难看到，F1模态的相速度变化不大，其值约为0.865～0.88，该值基本对应于条纹边界层广义拐点位置的速度。因此F1模态主要对应于无黏剪切模态。相对于F1，F2模态的相速度范围变化较大，其中下支点相速度接近于1，上支点相速度接近于0.88，即广义拐点的值。这种相速度的分布特征与二维边界层中第二模态扰动的分布特征十分相近。

图13 增长率随波数的变化关系Fig. 13 The relationship between the growth rate and the wavenumber

图14 相速度随波数的变化关系Fig. 14 The relationship between the phase velocity and the wavenumber

为了进一步研究两支解的性质，图15和图16分别给出了F1和F2两支曲线增长率最大时所对应的压力特征函数云图，分别为αr=0.7、ωi=0.00617668和αr=2.4、ωi=0.04775的情况。

图15 α r=0.7,ωi=0.00617668的特征函数Fig. 15 The eigenfunction for αr=0.7,ωi=0.00617668

由图15可以看到，对于F1模态，压力最大值主要集中y≈1的临界层附近，峰值点位置对应的基本流速度约为0.87，这是典型无黏剪切解特征。图16则给出了F2模态的特征函数分布，此时压力分布特征与F1大不相同，压力最大值集中在壁面附近。由图中的云图分布我们不难看到压力模态的模呈现随着法方向从大变小逐渐到零然后再变大最后逐渐在远场衰减到零的特征。这一分布特征与Mack第二模态的分布特征非常相近，是压缩模态的特征。只是此时失稳模态并非是二维模态而是三维的，这是由于流场由于存在微小的条纹，流场呈现一定的三维性，此时压缩模态受影响也使得最不稳定扰动变成了三维的。

图16 α r=2.4,ωi=0.04775的特征函数Fig. 16 The eigenfunction for αr=2.4,ωi=0.04775

图17给出了F2模态所对应的相对马赫数分布。从图中不难看到近壁区域整体上相对马赫数都小于−1，因此扰动波在此处是超声速的，此时方程（16）在局部表现为波动方程特性，从而导致压力的波动解。该特征和结论与二维基本流情况时的Mack第二模态基本相同。比较F2和F1模态的特征我们不难得出，在高超条件下，由于微小条纹的存在，流场中可以存在此时在三维边界层中类似二维情况中第一模态和第二模态的解。因此，我们将F1类的扰动称之为条纹边界层的第一模态，将F2类似的扰动称之为条纹边界层的第二模态。其中，第一模态F1失稳频率较低，增长率较小，第二模态F2增长率较大，但频率较高。

图17 Ma=4.5,cr=0.912时的相对马赫数云图Fig. 17 The relative Ma contours for Ma=4.5,cr=0.912

图18出了二次失稳扰动幅值B= 2时的增长率随波数变化的关系。可以看出，此时大幅值情况下，流场中主要存在四支解。图19为四支曲线对应的相速度分布。从图中不难看出，随着条纹幅值的增加，增长率有放大，其中第二模态的增长区域显著放大。同时，流场中失稳模态也较条纹幅值较低时更丰富。相速度分布行为随着条纹幅值的增加逐渐偏离二维边界层时的典型特征，但整体上第一模态的相速度较小，而第二模态的相速度较大。

图18 增长率随波数的变化关系Fig. 18 The relationship between the growth rate and the wavenumber

图19 相速度随波数的变化关系Fig. 19 The relationship between the phase velocity and the wavenumber

为了进一步说明问题，图20给出了F1、F3、F4最大增长率对应的压力云图，由于F1和F2曲线的压力云图类似，所以这里不再给出F2曲线对应的压力云图。图20（a）对应的是剪切解的情况，压力最大值在临界层附近。图20（b、c）对应的是压力最大值都集中在壁面附近。这与二次失稳扰动幅值为0.2时的结论类似，即第二模态类的解集中在壁面附近。其中第二模态类的解增长率较大，失稳频率范围和增长率随着条纹幅值增大而显著增大，而第一模态类的扰动相对变化范围要小。

图20 不同波数下的压力特征函数云图Fig. 20 The eigenfunction of pressure for different wavenumbers

图21和图22给出了相速度为F3和F4模态最不稳定扰动对应的相对马赫数云图。从图中可以看到，依然满足于近壁区域相对马赫数小于−1的性质，这使得扰动在该区域呈现出波动的特征，从而导致了非平凡解的产生。

图21 c r=0.8的相对马赫数云图Fig. 21 The relative Ma contours for cr=0.8

图22 c r=0.915的相对马赫数云图Fig. 22 The relative M a contours for cr=0.915

3 结论

本文采用线性稳定性分析和二次稳定性分析的方法，以高超声速平板边界层为研究对象，围绕边界层基频二次失稳条纹结构的稳定性开展了研究，得到结论如下：

1）首次失稳扰动的幅值对二次失稳类型有影响。当首次失稳扰动幅值较小时，亚谐模态占主导；反之，基频模态占主导，其特征函数在边界层内有较大的值。

2）基频模态的主要成分为流向的条纹结构，表现为流向涡。二次失稳扰动的幅值对全流场的流向涡和展向涡结构分布和强度都有影响。在相同的扰动幅值的前提下，流向涡显著强于展向涡，这表明此时流场中的动力学特征主要受流向涡结构影响。同时，首次失稳扰动幅值相同的条件下，二次失稳基频模态扰动的幅值对条纹结构（流向涡）的产生直接相关。二次失稳扰动的幅值越大，条纹结构越明显，幅值越大。

3）通过对条纹结构的无黏稳定性分析，结果表明流场中存在多失稳模态特征，且失稳特征与多种物理机制相关。其中低频失稳扰动的增长率较小，对应于拐点无黏不稳定性，而高频扰动的增长率较大，对应于由于局部相对马赫数绝对值大于1导致的压缩性的解。从解的特征解来看，条纹边界层中的失稳模态可以看作为二维边界层流动中第一模态和第二模态在三维边界层中的扩展。

综上所述，本研究通过对边界层内条纹结构的产生和无黏稳定性特征进行了初步研究，基于线性模型构造条纹基本流，研究了条纹结构中的无黏模态分布以及其物理主导机制。这推进研究高超声速边界层转捩过程中的非线性稳定性机制以及转捩中的典型流动结构的物理机理的理论认识。

这里需要额外指出的是，在本文的研究中，条纹结构的产生主要是基于线性模型构造的，并未考虑条纹产生激发的初值效应以及饱和特征，后续需要对此进行更详细的分析和研究。此外，对条纹边界层的稳定性模态在转捩过程中的意义应进一步基于直接数值模拟开展研究工作。