0 引言

高超声速三维边界层的转捩预测是飞行器设计中必须考虑的重要问题，特别是对长时间在大气中飞行的高速飞行器。由于湍流状态下壁面摩阻和热流比层流下大得多，只有准确地预测出转捩位置，才能准确计算阻力和表面热流，适当设计热防护措施。实际中飞行器飞行时常常有后掠角，这时在边界层内与外流流线垂直的方向上由于压力梯度会产生二次流动，称为横流。横流失稳是触发三维边界层转捩至湍流的主要原因。然而，目前对横流失稳触发转捩发生的机理，特别是对高超声速三维边界层^[1-4]，还并不清楚。

关于横流不稳定性的研究，早期主要针对低速(不可压或亚声速、个别低超声速)后掠翼和后掠平板的边界层。国际上最具代表性的是德宇航Bippes及其合作者、以及美国Saric及其合作者的工作，比较全面的综述可见文献[3-4]。国内也开展了不少关于横流感受性^[5]、横流稳定性^[6-8]以及横流转捩预测模型的研究^[9-10]。目前，已经有了比较一致的认识，即有横流的边界层中可以存在两种由横流不稳定性导致的流动——定常横流涡和低频的横流行进波。定常涡对壁面粗糙度非常敏感，很容易被激发。激发出的定常涡经历一定程度的线性增长后，在下游会由于非线性作用幅值出现饱和。饱和的定常涡导致平均流发生修正，使得在修正的平均流上高频扰动再次失稳，即出现二次失稳现象，促使转捩发生^[11]。横流转捩中定常涡和行进波的非线性饱和过程很长，在此过程中，扰动幅值变化非常缓慢，转捩发生的位置与开始发生非线性饱和的位置有相当长的距离。因此，若仍采用目前普遍用于Mack模态转捩的预测方法，即以首次不稳定波(横流定常涡或行进波)的幅值为基础作为转捩判据是不可行的^[11-14]。若以此为判据，任何在输入扰动和计算扰动演化过程中引入的很小的误差，也会给转捩位置的预测带来很大的影响。因此，横流失稳导致的转捩预测相比Mack模态失稳而言，更依赖于对转捩机理的认识和理解。

由于横流转捩是通过定常涡的二次失稳导致的，因此横流的二次不稳定性受到了广泛的关注^{[10, 13-15]}，研究大多集中在后掠翼上。尽管在实验方面已经有了不少有价值的结果^{[3-4, 15-16]}，而在理论和数值研究方面，大量的工作集中在横流首次失稳模态的非线性作用、二次失稳研究上，真正将其与转捩和湍流联系起来的却不多，且主要针对的是亚声速流动^[17-19]。

本文针对高超声速后掠钝板的三维边界层，采用直接数值模拟方法(DNS)和全局不稳定性(Bi-global)的分析方法，研究了横流定常涡模态的首次失稳和二次失稳；采用DNS计算了高超声速三维边界层横流失稳转捩至湍流的过程，分析了边界层转捩过程的机理。

1 基本流场 1.1 计算工况

本文选取半无限长后掠钝板为研究对象，其前缘为一半径为3.5 cm的圆柱形，如图 1所示。x、y、z分别为流向、法向、展向。飞行马赫数为6，来流迎角为0°，后掠角为45°。飞行高度为30 km，来流温度为226.5 K，以头部半径为特征长度的雷诺数为7.91×10⁴。

1.2 计算方法和边界条件

直角坐标系下三维N-S方程的守恒型控制方程为：

$ \frac{\partial \boldsymbol{U}}{\partial t}+\frac{\partial \boldsymbol{E}}{\partial x}+\frac{\partial \boldsymbol{F}}{\partial y}+\frac{\partial \boldsymbol{G}}{\partial z}=\frac{\partial \boldsymbol{E}_{v}}{\partial x}+\frac{\partial \boldsymbol{F}_{v}}{\partial y}+\frac{\partial \boldsymbol{G}_{v}}{\partial z} $

(1)

其中：U为守恒型变量，U=[ρ, ρu, ρv, ρw, ρ(c_vT+$\frac{1}{2} $(u²+v²+w²))]^T，ρ、u、v、w、T分别表示密度、流向速度、法向速度、展向速度和温度，c_v为定容比热；E、F、G分别为流向、法向和展向的无黏通量；E_v、F_v、G_v为相应的黏性通量。N-S方程中具体系数矩阵见文献[20]。

黏性系数采用Sutherland公式计算，具体为$ \mu=\mu_{\infty}\left(\frac{T}{T_{\infty}}\right)^{3 / 2}\left(\frac{T_{\infty}+110.4}{T+110.4}\right)$，下标∞表示来流量。热传导系数κ由$ \kappa=\frac{\mu c_{p \infty}}{P r}$给出，Pr=0.72，c_p∞为来流定压比热。

由于模型是二维的，展向无限长，因此∂/∂z=0。但由于流动存在后掠角，因此是三维的，即速度分量包含u、v和w。

采用基于有限差分的直接数值模拟程序计算基本流场。对流项离散采用五阶WENO格式，在边界处需进行降阶处理^[22]。对于边界点，采用三点二阶的偏心差分格式；对于次边界点，采用四点三阶的偏心差分格式。黏性项离散采用六阶中心差分格式，时间推进采用四步四阶Runge-Kutta法。壁面采用无滑移和绝热边界条件，远场条件给定自由来流。头部对称面处采用对称边界条件，出口采用线性外推边界条件。本文采用的基本流程序代码已成功地用于后掠钝板的基本流计算中，可见文献[21]。

由于横流涡扰动和二次失稳扰动是三维的，因此，对于定常横流涡的演化，以及非定常二次失稳模态的演化，均采用三维N-S方程。

扰动演化的计算中，入口给定具体的扰动波形式，出口采用嵌边区^[23]，壁面采用无滑移边界条件，外边界采用远场边界条件，展向采用周期边界条件。数值离散格式采用与基本流计算中相同的格式。

1.3 基本流场

计算域如图 2所示。流向计算域为100个头部半径(3.5m)，法向计算域包含激波。流向和法向网格数分别为1001和301，CFL数取为0.4。

求解二维的N-S方程直至所有网格点上的变量达到定常状态。将速度在边界层外缘的势流方向和垂直于势流的方向进行分解，得到势流方向速度分量u_s和横流速度分量u_c的分布，如图 3所示。可以看到，流动中横流的量级为10^-2。从头部到下游，横流强度逐渐减弱。

2 首次失稳——横流定常涡的演化

首先对基本流进行线性稳定性(LST)分析。在平行流假设下，小扰动可以写成如下形式：

$ \varphi(x, y, z, t)=\hat{\varphi}(y) \mathrm{e}^{\mathrm{i}(a x+\beta z-\omega t)} $

(2)

其中$ \hat{\varphi}(y)$为扰动的形函数，扰动在流向和展向和时间方向具有波的形式。考虑空间模式，即扰动沿空间增长，那么频率ω为实数，α和β为复数，其虚部的相反数分别表示扰动波在流向和展向的增长率。将上式带入到线性化的N-S方程的扰动形式中，并忽略基本流量的流向导数及法向速度，得到准平行流假设下的线性稳定性方程:

$ L(\hat{\varphi})=0 $

(3)

在壁面和远场扰动为零。这样方程结合齐次边界条件，构成特征值问题。通过求解可以得到边界层扰动波的色散关系和扰动的特征函数。

横流有两种失稳模态，一种是定常的横流涡，另一种是低频的横流行进波。在高空低背景扰动情况下，定常横流涡被认为是三维边界层中占据主导作用的失稳模态^[3]。因此，在首次失稳中，我们暂时只关注定常横流涡。本文选取β₀=3这一定常涡扰动，来计算其在流场中的演化情况。

由于基本流中引入的是三维扰动，因此首先将基本流在展向延拓成一个波长，即L_z=2π/β₀，然后在计算域的入口引入扰动来研究首次失稳扰动演化，扰动形式为:

$ \varphi(y, z)=A_{0} \hat{\varphi}(y) \mathrm{e}^{\mathrm{i} \beta_{0} z}+\mathrm{c.c} $

(4)

其中，A₀代表入口扰动幅值，取为10^-3，$\hat{\varphi}(y) $代表线性稳定性分析得到的β₀扰动的特征函数。

用DNS计算这一扰动在流场中的演化，直至流场定常。图 4给出了定常横流涡的流向速度云图。可以看出在上游，流向速度在展向分布比较均匀，随着扰动向下游演化，壁面附近的流体逐渐被卷起，形成涡结构。

将DNS得到的扰动场沿展向做傅里叶变换，可以得到不同展向波数的各阶谱的速度扰动幅值沿流向的变化，如图 5所示。图中(0, 1)表示加入的定常基本波(β₀=3)，(0, 2)和(0, 3)分别表示2倍展向波数和3倍展向波数的波。可以看出，在上游基本波的幅值快速增长，在x=20处由于幅值较大开始非线性作用产生基本流修正(0, 0)和二次谐波(0, 2)。当基本流修正增长到幅值约为0.1左右时，(0, 1)波幅值达到饱和，随后基本流修正也开始出现饱和。再向下游演化，扰动幅值不再增大，扰动维持在接近饱和的状态。

由于定常涡的非线性作用会产生基本流修正，修正后基本流的失稳特征将会出现变化。因此我们将图 5中得到的基本流修正项(0, 0)加入到基本流中，再次分析稳定性特征。图 6给出了新的基本流在不同流向位置的失稳区内增长率等值线图。同时，作为对比，也给出了原始基本流的稳定性分析结果。可以看出，在x=20的地方，新的基本流增长率分布与原始基本流一致，这是由于基本流修正的幅值还比较小，对稳定性影响不大。在下游x=30的位置，相比原始基本流，增长率以及失稳区域明显减小；继续向下游，在x=40处新的基本流几乎不再失稳。说明由于基本流修正的作用，流场不稳定区域减小，这意味着忽略扰动展向的变化，只考虑(0, 0)基本流修正，只会得到更稳定的基本流。但是转捩的发生依赖于流场在更大范围内的失稳，说明具有展向尺度的二次失稳扰动可能是促进转捩的关键因素。

3 二次失稳及转捩至湍流的过程

首先对首次失稳的流场进行二维全局稳定性分析。全局稳定性分析方法可以研究基本流在两个方向是快变而在另一个方向是慢变的三维流场的稳定性问题。扰动可以写成如下形式：

$ \varphi(x, y, z, t)=\hat{\varphi}(y, z) \mathrm{e}^{\mathrm{i}(a x-w t)} $

(5)

其中$ \hat{\varphi}(y, z)$为扰动的形函数，扰动在流向和时间方向具有波的形式。与LST类似，经过推导可以得到关于$ \hat{\varphi}(y, z)$的二维特征值问题。具体公式见文献[21]。

对首次失稳后流场的某一个流向站位的剖面作全局稳定性分析，可以找到不稳定的模态，称为二次失稳模态。

在x=36处，即首次失稳波接近饱和(幅值约为0.2)位置处进行全局不稳定分析，可以找到两个二次模态，分别对应于z模态和y模态^[13]。图 7(a)(b)分别给出了z模态和y模态的特征函数$|\hat{\rho}| $。对于z模态，其最大值的位置在法向y≈1.5，即二次失稳扰动发生在饱和涡结构的侧部，这一类型的模态是由流场中展向速度剪切导致的。而对于y模态，其最大值的位置在法向y≈1.6，对应饱和涡结构的顶部位置，这一类型的模态是由流场中法向速度剪切导致的。图 8给出了x=36处，时间模式的全局稳定性分析得到的z模态和y模态的增长率随频率的变化。显然，z模态具有更大的增长率，应是更容易触发转捩的二次扰动。因此，在本文后面的计算中，只考虑更不稳定的z模态。

以x=36为计算域入口，引入全局稳定性分析得到的两个不同频率(ω=5和ω=4)的二次失稳模态。为了使转捩尽快发生，它们的初始幅值都选为A₀=0.01。使用DNS计算二次失稳扰动演化，其中流向计算域为x=36~96，展向计算域仍为一个基本波的波长。流向、法向和展向计算域内的网格点数分别为2401、151和61。

图 9给出了计算得到的流场瞬时速度等值面。可以看到在上游靠近入口的位置，二次失稳扰动幅值很小，流场的涡结构主要呈现首次失稳后接近饱和的定常涡结构。随着扰动向下游演化，高频的二次失稳扰动增长起来，饱和涡结构开始扭曲，继续向下游演化，规则的涡结构破碎，转捩完成，最终形成湍流。

图 10给出了不同流向位置的二次失稳扰动速度云图以及首次失稳横流涡速度等值线图。从图中可以看出在x < 51的范围内，二次失稳扰动的形状与首次失稳横流涡的速度等值线形状相近，说明首次失稳流场的剪切导致流动出现二次失稳。而在x>51的下游，明显观察到二次失稳扰动形状逐渐变得不规则，说明流场中不同频率(展向波数)的高次谐波增长起来，流场开始出现转捩。

为了进一步分析流场的转捩过程，首先分析二次失稳扰动各阶谱的演化特征。图 11给出了傅里叶分析得到的不同频率的扰动幅值沿流向的变化，图中ω=4和ω=5的波表示入口加入的两个基本波，无量纲频率ω=0、1、2、3的波表示由基本波相互非线性作用激发的扰动，更高频率的扰动在图中用灰线表示。可以看到，在x < 44的范围内，入口加入的基本波快速增长，同时由于它们的幅值较大，非线性作用得到的ω=0和ω=1的扰动被激发并且快速增长。继续向下游演化，其它频率的扰动也被激发，并且它们都会经历一个快速增长阶段，最终在x=55附近各阶谱的幅值趋于稳定，其中频率为0的谱幅值最大，并且频率越大，幅值越小。

图 12给出了壁面摩擦系数(C_f)沿流向的变化，其中紫线是层流解的结果，红线是采用SST模型得到的湍流剖面的结果，而黄线是首次失稳后横流涡的壁面摩擦系数，黑色实线是二次失稳后的壁面摩擦系数。可以看到首次失稳后壁面摩擦系数有一定的抬升，但是抬升量不大，并没有达到湍流的壁面摩擦系数值。说明首次失稳幅值增大到一定程度对平均流产生修正，导致壁面摩阻上升，但是由于首次失稳波已经饱和，其本身无法再使得壁面剪切增大，使得C_f进一步抬升，因此转捩并未发生。而加入二次失稳扰动之后，壁面摩擦系数在40 < x < 60的范围内快速抬升，并达到湍流剖面的摩擦系数，说明在x=40下游，流场开始转捩，并逐步演化为湍流。这是由于引入高频的二次失稳波后，在非线性作用下，各阶扰动经历快速的增长(见图 11)，不断从平均流中吸取能量，从而进一步修正平均流使得C_f更加快速抬升。这样的结果，导致壁面摩擦系数超过了相应的湍流的估计值，因而产生了“过冲”^[24]的现象。当各阶扰动的幅值达到一定程度，不再继续增长，“过冲”现象逐渐减弱或消失，C_f曲线会逐渐趋向于充分发展湍流的估计值。同时在C_f曲线快速抬升的阶段，低频扰动快速增长并且幅值最大，说明低频扰动的增长主导平均流修正及C_f曲线的抬升。

将展向计算域拓展为两个基本波波长(2L_z)，其他条件保持不变，计算了二次失稳扰动的演化，图 12中黑色虚线给出了壁面摩擦系数(C_f)沿流向的变化，与展向计算域为一个基本波波长(黑色实线)结果相比，二者几乎一致，特别是在C_f曲线抬升阶段，二者完全重合，说明展向计算域对横流涡转捩过程影响不大。

图 13给出了经过Van Driest变换后的平均速度剖面的分布，图中虚线分别表示壁面律(u⁺=y⁺)和对数律(u⁺=2.5lny⁺+4)。可以看出在61 < x < 91的范围内，平均速度剖面与对数率和壁面率符合得很好，而且随着向下游演化，剖面越来越接近于对数律，说明剖面逐渐向完全湍流发展。

4 结论

本文以马赫数为6的后掠钝板边界层为研究对象，采用直接数值模拟方法研究了横流定常涡从首次失稳、二次失稳到转捩发生的过程。得到以下结论：

1) 横流定常涡的非线性作用引起平均流修正，可使壁面摩擦系数曲线有一定程度的抬升，但是横流涡饱和后壁面摩擦系数不再增大，转捩不会发生。

2) 在横流涡非线性饱和(幅值约为0.2)的基础上，发生二次失稳，产生高频的不稳定波。二次失稳波由于非线性作用产生的低频谐波及定常涡迅速增长，促使壁面摩擦系数急剧抬升，同时首次失稳的饱和横流涡结构破碎，最终促使转捩发生。

致谢: 此项工作得到天津大学周恒院士的关心和指导，天津大学赵磊博士提供了DNS程序及全局稳定性分析程序，并在国家超算天津中心“天河一号”完成计算，特此致谢。