多孔壁面对高速边界层最优增长条带二次失稳的影响规律

王宇天，刘建新，王晓坤，李晓明，

1.空军航空大学，长春 130022

2.天津大学 机械工程学院，天津 300072

高超声速飞行器具有重要的军事和经济价值，是当今各航空航天大国的研究热点之一。在临近空间飞行雷诺数范围内，边界层通常要经历由层流到湍流的转捩过程。由于飞行器摩擦阻力、热防护、发动机性能等均与边界层流动状态密切相关，因此边界层转捩是飞行器设计中必须要解决的基础理论问题，具有重要的工程应用价值与学术意义［1-4］。

边界层转捩描述的是扰动产生与失稳演化直至突变为湍流（Breakdown）的全过程，其中扰动的产生由感受性机制决定，失稳演化则是流动稳定性理论的范畴。线性稳定性理论（Linear Stability Theory，LST）、抛物化稳定性方程（Linear Parabolized Stability Equations，LPSE）等可以预测模态扰动的线性发展过程，但人们在Poiseuille、Couette、Blasius 边界层等流动实验中发现，扰动在中性曲线以外的稳定区域依然出现了较大增长，甚至转捩。由于扰动间的非线性作用只是扰动能量的重新分配，并不能产生扰动能量的净增加，所以必然存在一种不同于模态增长的线性增长机制，通常称为非模态或瞬态增长。通过数学分析可知，瞬态增长由线性扰动方程算子矩阵的非正规性所致［5］。由于非正规矩阵的特征向量并不彼此正交，因此即使所有特征向量均为衰减模态（主要为连续模态），也可通过一定的线性组合产生增长扰动。而在物理上，瞬态增长通常由Lift-up 机制［6］进行解释：在流向涡的作用下，壁面附近的低速流体被抛向边界层外，并将边界层外的高速流体带向壁面，进而形成了流向速度在展向呈周期性分布的条带结构。

瞬态增长中各衰减模态的线性组合系数，即幅值与相位信息，同样由感受性过程决定，所以瞬态增长率受外部扰动情况（如高幅值自由流扰动、粗糙元等）的影响。但在特定的流场参数条件下，可以通过最优扰动理论进行求解，得到在有限时间或空间范围内实现最大瞬态增长的最优扰动。已有的研究表明［7-8］，从不可压到高超声速边界层，最优扰动通常为定常流向涡的形式，其通过非线性发展增长至较大幅值的条带结构时，新的平均流剖面在展向上将出现拐点，从而发生二次失稳触发转捩。因此从转捩预测的角度上看，这是一种更为保守的亚临界转捩触发阈值。

瞬态增长理论完善了传统的线性稳定性理论框架，可以对超越了模态增长理论的旁路转捩中一些亚临界转捩现象进行定性解释，如经典的Poiseuille 管道流动［9］。除此之外，瞬态增长理论的一个重要应用是再入返回舱等大钝度飞行器的“钝头体悖论”问题，即发生在弓形激波后亚声速区的转捩现象。根据线性稳定理论，此处的顺压梯度作用将抑制Tollmien-Schlichting（T-S）波的增长，而凸面曲率也将抑制Görtler 涡，同时横流速度较小使得横流模态可以忽略，因此人们在很长一段时间内无法找到实验中出现转捩的原因，并且近期在李强等［10］的顿头平板转捩实验中也出现了这一现象。近几年来，NASA 课题组Paredes 等［11］对最优扰动进行了细致的研究，发现在驻点附近存在较强的瞬态增长，并认为壁面微小粗糙度便可能成为诱导转捩的主要原因。同时瞬态增长理论也同样能够对钝头体“转捩反转”现象［12］进行定性解释［13-14］，计算发现在钝头与锥身连接处附近存在较强的瞬态增长，并且增长率随着钝度增大而增大，因此能够在Mack 第一、第二和熵层模态增长率较小的情况下触发转捩。尽管连接处的微小粗糙元不一定能够恰好激发最优扰动，但也为这一转捩过程预测提供了一种科学的阈值指导，并且还需要进一步的实验验证。

随着对转捩机理的深入认识，各种延迟转捩技术应运而生［15］，其中多孔壁面能够有效抑制Mack 第二模态，可作为一种推迟转捩的被动控制方式［16］。目前，微孔结构可在高超声速飞行器表面涂层材料的制备工艺中自然形成，因此在工程领域受到了广泛关注。但已有的研究均针对二维边界层的局部稳定性问题，对于三维边界层二次失稳模态的控制效果还不明确。本文通过建立最优扰动的求解方法，并以最优增长条带形成的三维边界层为研究对象，研究多孔壁面对全局稳定性（Bi-Global）的影响规律，评估其在亚临界转捩中的控制效果，以期为多孔涂层的布置方案提供一定的科学指导。

1 数值方法及验证

本文以超声速绝热平板边界层流动为研究对象，基本控制方程为直角坐标系下的可压缩Navier-Stokes（N-S）方程，选取以下特征量与特征尺度对物理与时空变量进行无量纲化：

式中：x、y、z分别为流向、法向与展向坐标；u、v、w分别为流向、法向与展向速度；t为时间；ρ为密度；T为温度；p为压力；上标“*”代表有量纲物理量；下标“∞”代表自由来流值。特征长度取0.001 m，故参考雷诺数定义为。黏性系数μ由Sutherland 定律求得，传热系数κ则可通过普朗特数Pr得到：

流动稳定性分析采用原始变量，在非守恒型N-S 方程的基础上，通过理想气体状态方程p=消去压力项p。将瞬时物理量q=[ρ，u，v，w，T]T表示为基本流与扰动量之和：q=+q′，将其代入N-S 方程后减去基本流满足的方程部分，并忽略非线性作用项，便可得到如下紧致形式的线性扰动方程：

式中：系数矩阵Γ、A、B、C、D、Vxx、Vyy、Vzz、Vxy、Vxz、Vyz的具体表达式可参见文献［17］，并以可压缩边界层相似解作为基本流动［18］。为了理论方法介绍的完整性，LPSE 及其伴随方程的表达式详见附录A［19］。

1.1 最优扰动方法

传统的LST 通过求解常微分方程特征值问题，得到的是扰动在时间与空间趋近于无穷时的渐进演化行为。而瞬态增长则通过求解常微分方程的初值问题，得到扰动在有限时间或空间内的代数增长。目前，最优扰动有2 种求解方法：一种是在平行流假设条件下，基于LST 方程的奇异值分析［20-21］；另一种是考虑边界层弱非平行性的，基于伴随方程的优化方法。基于LST 的平行流最优扰动分析只能定性地给出最优扰动的展向波数范围，而对于扰动放大倍数G的计算可能存在较大误差，因此本文将在考虑边界层弱非平行性的LPSE 方程的基础上，利用Lagrange 乘数法求解扰动增长的数学最优化问题［22-23］。

最优增长条带在流向上的尺度较大，因此扰动幅值与相位信息的变化可以全部反映在形状函数中，即可使式（A1）中α=0。采用Hanifi等［20］建议的扰动能量范数E的计算准则，以内积的形式给出：

设入口x0处最优初始扰动为，以某一指定位置x1处的扰动的能量范数的增长倍数作为目标函数：

其中：∇代表取梯度。为使式（8）恒满足δF=0，式（8）中两项内积需同时为0。第1 项自动满足，即求解齐次LPSE 方程（式（A2））：

式（8）中第2 项内积可利用内积定义（式（A5））与分步积分将其变换为

式中：L†=0即为齐次伴随LPSE方程，此时分步积分产生的边界项仅保留式（10）中的第二、第三项。式（10）进一步变换为

要使式（11）恒等式成立，等号左边两项需同时为0，这样便可以得到在x0与x1处初始扰动与伴随向量的关系式为

对于线性问题，c0、c1为归一化系数使得。在壁面处为奇异矩阵，因此可由法向动量方程求得。在入口位置x0处，式（12）中的流向扰动速度优化条件为

由于伴随密度扰动在壁面处无预设齐次条件，将导致扰动速度在壁面处非零，本文采用Tempelmann 等［23］的建议方法，直接假设，下标“w”代表在壁面处的值，该方法与Zuccher 等［24］的插值方法具有相同的效果。基于以上优化条件，便可通过迭代求解得到最优扰动，具体步骤为

步骤1选取任意满足齐次边界条件的扰动作为初始扰动，可利用LST 算子进行当地求解。

步骤2利用LPSE 方程（式（A2））从x0→x1推进求解得到。

步骤3由式（13）得到伴随向量，利用伴随LPSE 方程从x1→x0推进求解得到。

步骤4由式（12）与式（14）得到新的初始扰动0。

重复步骤2～步骤4，直至目标函数J()=G收敛。计算结果表明该优化系统具有吸引子，通常迭代3～4 步便可以收敛至10-4。

以x1=1 m 为优化目标位置，在Ma=3.0 工况下，当展向波数β=0.25 时取不同的初始位置x0，或当x0=0.36 m 时取不同展向波数β，最优扰动增长倍数的计算结果如图1 所示，其中纵坐标以来流单位雷诺数Re∞做归一化处理。通过与Paredes 等［8］的结果对比验证了程序的正确性，并且可以看出最优扰动增长倍数与优化初始位置的选取密切相关。

图1 最优扰动增长倍数的验证Fig.1 Validation of optimal disturbance amplification

1.2 全局稳定性分析

不同于局部LST 分析，对于基本流在展向z同样是快变，而流向x仍是慢变的三维边界层，需采用全局稳定性分析（Bi-Global）。此时扰动可以写为如下行进波的形式：

方程在法向上同样采用10 阶中心差分格式，展向上采用适用于周期性边界条件的Fourier 谱方法离散，结合齐次边界条件（式（A4）），便可通过求解广义特征值问题得到二次失稳扰动的空间增长率。对于大型稀疏矩阵，本文采用隐式自启动Arnoldi 方法。

对于多孔壁面，由于孔隙直径一般为几十微米量级，孔隙间的相互作用较小，因此其诱导的流向与横向扰动速度可忽略不计，但对于法向速度与温度扰动具有半透射作用。Fedorov 等［25］通过引入声导纳系数将其与压力扰动相关联，扰动在壁面处的边界条件变为

式中：Ac、Bc为复导纳系数，其与涂层材料、平均流以及扰动特征相关。由于涂层厚度（孔隙深度）h通常远大于孔隙直径2r，可将其看作细长管。基于声波在细长管内的传播理论，推导出复导纳系数的计算公式为

式中：ϕ为开孔率，当孔隙截面为圆形时，易知最大开孔率ϕ→π/4；下标“w”代表基本流动在壁面处的值；Jn则为n阶第一类Bessel 函数。

已有的直接数值模拟（Direct Numerical Simulation，DNS）与实验结果均证明该理论边界条件的合理性［26］。本文选取Fedorov 等［25］针对二维边界层的计算结果进行验证，流场条件为：Ma∞=6.0，=243.9 K，=1.4，计算结果如图2 所示，验证了程序的正确性。

图2 多孔壁面条件下第二模态增长率的验证Fig.2 Validation of growth rate of the second modes over porous wall

1.3 直接数值模拟

对于最优增长扰动的非线性演化过程，本文采用高精度有限差分方法求解N-S 方程得到。计算程序在中国科学院李新亮等的开源程序OpenCFD-SC［27］的基础上发展而来。当采用层流边界层相似解作为基本流时，由于相似解并不是N-S 方程的稳态解，因此需要在N-S 方程中引入人工源项以维持基本流保持不变［28］。在对无黏通量进行空间离散时，采用截断误差最小的Steger-Warming 通量分裂格式得到正负通量，再采用7 阶迎风格式进行离散，黏性通量采用8 阶中心差分格式进行离散。时间推进采用具有TVD（Total Variation Diminishing）性质的3 步3阶Runge-Kutta 方法。在壁面处，速度与温度采用无滑移等温条件，并假设∂p/∂y=0 作为壁面压力条件。

2 计算结果与分析

已有的研究表明，多孔壁面通常只对无黏的Mack 第二模态起稳定作用。但根据3 层结构理论［29］，满足的Mack 第一模态同样为无黏模态，但多孔壁面却起促进作用。为了研究多孔壁面在不同中高马赫数范围内对二次失稳模态的影响规律，本文将研究来流马赫数为3.0 与6.0 这2 种情况，来流总温=333 K，Pr=0.7，单位雷诺数Re∞=106m-1。

下面将首先利用伴随方法得到线性最优扰动，然后通过DNS 计算不同初始幅值最优扰动的非线增长过程得到最优增长条带，再以某一流向位置的y-z截面作为新的三维边界层剖面进行二次失稳分析。

由于条带基本流在展向上具有周期性，对于基频（Fundamental）或亚谐型（Subharmonic）二次失稳模态，可以根据其对称或反对称性质分为奇（Sinuous 型）、偶（Varicose 型）模态，二者具有不同的增长率。偶模态的扰动分量在一个周期内呈对称分布，为反对称分布，而奇模态的对称性与之相反。本文在实际计算中，为更好地区分奇/偶模态，展向上采用对称离散方式，即只需保留一半的展向网格即可。通过网格分辨率验证，Ny×Nz=201×64 已满足计算需求。

2.1 最优扰动

图3 分别为Ma∞=3.0 与Ma∞=6.0 流场条件下，不同起始位置的最优扰动增长倍数随展向波数β的变化规律。可以看出在高超声速条件下，最大瞬态增长倍数相对降低，并且产生最大瞬态增长的展向波数β减小。

图3 初始位置x0 与展向波数β 对最优扰动增长倍数的影响（x1=1 m）Fig.3 Effect of initial position x0 and spanwise wavenumber β on optimal disturbance amplification（x1=1 m）

由图3 可以看出，选取不同的初始位置将改变扰动的最优增长倍数。从感受性的角度看，由于前缘区域的边界层较薄且具有较强的非平行性，因此自由来流涡扰动或壁面粗糙度具有更大的感受性系数，因此由较高幅值环境扰动所引起的瞬态增长通常发生在中性曲线下支前。为不失研究对象的一般性，本文选取靠近前缘的初始位置，即x0=0.09 m，同时在Ma∞=3.0 与Ma∞=6.0 流场条件下，分别选取展向波数β=0.25 与β=0.10 的最优扰动作为初始扰动，初始与终点位置处扰动的形状函数如图4 所示。可以看出初始最优扰动具有流向涡的形式，也就是展向与法向扰动速度较大，其增长由Lift-up 机制主导。并且，随着马赫数增大至高超声速，最优初始扰动的温度和密度扰动分量与横向和法向扰动速度具有相同量级。当流向涡发展为条带时，温度与密度扰动将超过流向速度扰动成为扰动的主要分量，这与定常Görtler 模态相似［30］。

图4 初始与终点位置处的最优扰动形状函数（x0=0.09 m，x1=1 m）Fig.4 Shape function of optimal disturbances at initial and final locations（x0=0.09 m，x1=1 m）

2.2 马赫数3.0 工况

当最优扰动初始能量（式（4））E0=0.000 4时，x=0.2 m 处的截面条带如图5 所示。图中深灰色线为流向速度等值线，其大小从0.1 依次递增0.1 至1.0，可以看出此时展向流场已具有一定的非平行。

图5 x=0.2 m 处的流向速度等值线与梯度云图（E0=0.000 4）Fig.5 Contours of isolines and gradient of streamwise velocity at x=0.2 m（E0=0.000 4）

多孔壁面条件设为：开孔率ϕ=0.5，孔隙半径r=0.1 mm，约为当地边界层厚度的2%。研究表明［26］：当孔隙深度大于0.8 倍的当地边界层厚度时，控制效果基本不再随孔隙深度而改变。设置孔隙深度h=2 mm，在本文中的超声速条件下|Λh|＞3，根据式（19）中双曲正切函数tanh 的性质可知tanh(Λh)→1，因此其与Λh→∞时等效。图6 为基频与亚谐频失稳模态的增长率-αi随角频率ω的变化规律，以及一维基本流的LST增长率。当最优扰动条带的幅值较小时，二次失稳模态的增长率未超过第一模态，并且基频偶模态（FV）增长率高于基频奇模态（FS），而对于亚谐型失稳则是奇模态（SS）高于偶模态（SV）。与第一模态类似，多孔壁面促使二次失稳模态增长率变大，并且对亚谐失稳的作用更强。

图6 有无多孔壁面作用下的二次失稳模态增长率（E0=0.000 4）Fig.6 Growth rates of secondary instability modes with and without porous wall（E0=0.000 4）

图7 ω=0.06 时二次失稳模态流向扰动速度模值Fig.7 Modulus of streamwise disturbance velocity of secondary instability at ω=0.06

当提高最优扰动的初始能量至E0=0.01时，x=0.2 m 处截面条带如图8 所示。此时条带基本流进一步扭曲，展向剪切增大了10 倍以上。并且可以注意到等值线在中心处出现下凹区，这在Görtler 条带［30］中并未出现，将产生新的二次失稳模态。

图8 x=0.2 m 处的流向速度等值线与梯度云图（E0=0.01）Fig.8 Contours of isolines and gradient of streamwise velocity at x=0.2 m（E0=0.01）

图9 为二次失稳模态的增长率-αi随角频率ω的变化规律，此时二次失稳模态的增长率远高于第一模态，将成为促发转捩的主要因素。与不可压Blasius 边界层［31］最优扰动条带及Görtler涡［30］的二次失稳相同，亚谐模态与基频模态具有相近的增长率，且低于基频模态，因此在转捩预测时通常更关心基频失稳。与图6 相似，基频奇模态（FS）增长率高于基频偶模态（FV），而亚谐型失稳模态则与之相反。

图9 有无多孔壁面作用下的二次失稳模态增长率（E0=0.01）Fig.9 Growth rates of secondary instability modes with and without porous wall（E0=0.01）

图10 ω=0.085 时二次失稳模态流向速度扰动实部Fig.10 Real parts of streamwise disturbance velocity of secondary instability modes at ω=0.085

2.3 马赫数6.0 工况

在进行Bi-Global 分析之前，首先对高超声速二维基本流动进行稳定性分析。图11 为x=0.2 m 处LST 分析得到的Mack 第二模态增长率-αi随角频率ω的变化规律。多孔壁面条件设为：开孔率ϕ=0.5，孔隙半径r=0.15 mm，约为当地边界层厚度的1.5%，孔隙深度h=2 mm。从图中可以看出：多孔壁面对低频模态（即第一模态）有促进作用，而对高频模态（即第二模态）有明显的抑制作用，并且对于三维模态具有相同的控制效果。在该工况下，第二模态由慢模态［32］（即第一模态）发展而来，图12 为二维快/慢模态相速度c=ω/αr随角频率ω的变化规律，多孔壁面对相速度的影响整体较小，并且对慢模态的影响主要体现在第一模态频率区，对第二模态无明显影响。快/慢模态的同步点频率［32］约为0.237，而图11 中控制效果的转折点频率约为0.22，二者较为接近，并且对于三维模态具有相同的规律。

图11 x=0.2 m 处LST 分析得到的Mack 模态增长率Fig.11 Growth rates of Mack mode analyzed by LST at x=0.2 m

图12 x=0.2 m 处LST 分析得到的快/慢模态相速度Fig.12 Phase velocity of fast/slow modes analyzed by LST at x=0.2 m

当扰动初始能量E0=0.004 时，x=0.2 m处截面条带如图13 所示，图中深灰色线为流向速度等值线，其大小从0.1 依次递增0.1 至1.0，流向速度法向梯度与展向剪切与Ma∞=3.0 工况具有相似的分布（图5）。

图13 x=0.2 m 处的流向速度等值线与梯度云图(E0=0.004)Fig.13 Contours of isolines and gradient of streamwise velocity at x=0.2 m (E0=0.004)

高超声速条带的二次失稳情况更加复杂，由于基频失稳通常在转捩中占主导地位，下面将重点研究基频失稳。该流场条件下，共出现2 个奇模态（S1，S2），3 个偶模态（V1，V2，V3）。光滑壁面条件下，以流向扰动速度模值||的最大值做归一化处理，图14 为二次失稳模态的流向扰动速度模值||及温度扰动模值||的分布云图。对于高超声速流动，温度扰动幅值远高于速度扰动，并且速度扰动峰值靠近壁面区，而温度扰动峰值位于靠近边界层外缘的临界层。

图14 ω=0.28 时二次失稳模态流向速度与温度扰动模值Fig.14 Modulus of streamwise velocity and temperature disturbances of secondary instability modes at ω=0.28

图15 为二次失稳基频奇/偶模态的增长率-αi随角频率ω的变化规律。V1 模态增长率高于S1 模态，是二次失稳的主导模态，而S2 与V3增长率较小可以忽略。V2 模态是一个较为特殊的模态，当ω＜0.27 时，其增长率明显小于V1 模态与S1 模态，但其高频区的增长率较高，可能成为转捩中的最危险的模态。当采用多孔壁面边界层条件时，其影响规律与Mack 模态相似，即只对高频二次失稳模态起抑制作用，而对低频二次失稳模态起促进作用。S1 模态与V1 模态的转折频率分别为0.215 与0.19，这与在图11 中Mack模态的转折频率相近。V2 模态的转折频率为0.28，远高于其他模态，但V2 的快速增长只能出现在高频区域，相对的低频区域即使有促进作用，由于其增长率远小于S1 模态与V1 模态，因此不会对转捩产生影响。

图15 有无多孔壁面作用下的二次失稳模态增长率（Ma∞=6.0）Fig.15 Growth rates of secondary instability modes with and without porous wall（Ma∞=6.0）

综合以上研究结果可以看出，虽然二次失稳模态在本质上都属于无黏失稳，特别是对于奇模态，其应是由展向剪切主导的失稳模态，但在不同的流场条件下，多孔壁面对其影响规律却不尽相同，并与无条带时的局部稳定性特征密切相关。Song 等［33］针对Görtler 涡的二次失稳问题，证明了二次失稳模态可以由上游的Mack 第一、第二模态演化而来，由此可以推断二次失稳模态也将保留二维边界层中Mack 模态的稳定性特征。但在物理空间中，上游某一特定频率的模态扰动演化到下游也应为一特定扰动模态，但在Bi-Global 分析中却会得到多个二次失稳模态，这其中的关联机制还尚不清楚，需要借助模态分解技术［28］做进一步详细研究。

3 结论

本文针对瞬态增长理论，利用Lagrange 乘数法与变分法，推导出以伴随抛物化稳定性方程为约束条件的优化系统，建立并实现了最优扰动的数值求解方法。以可压缩平板边界层流动为研究对象，针对Mack 第一、第二模态分别主导转捩的超声速（Ma∞=3.0）与高超声速（Ma∞=6.0）2 种工况开展研究。计算表明：最优扰动具有流向涡的形式，并通过Lift-up［6］增长机制发展为条带结构，并且在高超声速条件下，其与定常Görtler 模态相似［30］，温度与密度扰动将超过速度扰动成为扰动的主要分量。

以有限幅值最优扰动非线性发展形成的三维边界层为研究对象，并利用Fedorov 等［25］建立的等效边界条件，研究了多孔壁面对最优增长条带二次失稳模态的影响规律。结果表明：与第一模态相同，多孔壁面对超声速条带的二次失稳扰动只起促进作用。对于高超声速条带，多孔壁面对第一模态频率范围内的二次失稳扰动为促进作用，对第二模态频率范围内的二次失稳扰动起抑制作用，并且转折频率接近局部快/慢模态的同步频率。虽然目前还无法对其中的物理机制做出解释，但对工程应用中多孔涂层的布置方案具有一定的指导意义。

附录A：

基于边界层慢变假设，将扰动分解为快变的波数函数和慢变的形状函数两部分：

扰动在边界处满足Dirichlet 边界条件：

在法向上采用10 阶中心差分离散后，结合齐次边界层条件（式（A4）），此时便将扰动方程求解由LST 的特征值问题转变为偏微分方程的初边值问题，在流向上采用隐式欧拉格式推进求解即可得到扰动的空间发展过程，详细计算方法与残余椭圆性问题可参见文献［19］。

伴随LPSE 方程可以通过Green-Lagrange恒等式推导而来，首先定义向量内积为

式中：上标“H”表示共轭转置。本文以上标“†”表示相应的伴随变量，将伴随向量与齐次LPSE 方程（式（A2））作内积，通过分步积分可以得到如下Green-Lagrange 恒等式：

式中：L†=0 即为伴随LPSE 方程；为分步积分产生的边界项，伴随向量同样满足齐次边界条件（式（A4））。伴随算子L†为

伴随LPSE 方程求解采用与LPSE 相同的数值方法，只是改为由下游向上游推进求解。