0 引言

对于超声速飞行器的设计，层流-湍流转捩是必须要考虑的问题，也是流体力学研究的热点。在超声速飞行条件下，飞行器钝头体附近产生弓形激波。当气流穿过前缘附近激波角变化较大的弓形激波时，激波后的熵值变化较大，形成一层横向熵梯度较大的区域，称之为熵层。气体向下游流动，将包围着飞行器表面，从而增加表面气动热流^[1]。因此研究熵层的流动稳定性，对于飞行器气动热防护的设计具有实际意义。

早期国外对钝体模型的实验研究发现，熵层对边界层稳定性的影响具有重要作用。Stetson^[2]通过实验研究发现，在边界层外存在增长着的扰动，这表明熵层存在不稳定性。Lysenko^[3]通过钝板实验也测量出熵层不稳定模态，同时指出这个不稳定模态影响着边界层内的第一模态。为了证实Stetson实验等关于熵层不稳定性的存在，Dietz & Hein^[4]利用数值方法计算钝板模型，指出熵层有较小的增长率，其扰动演化与文中^[4]Aerodynamisches Institute的实验结果吻合。Fedorov^[5]根据Yakura^[6]渐近法分析得到了熵层基本流方程表达式，在此基础上进行线性稳定性分析，得出同样的结果。

早在实验^{[3, 7-8]}研究钝度对边界层转捩的影响时发现，在小钝度范围内增加钝度，转捩位置往下游移动；当钝度增加到一定大小时，转捩位置反而提前，这种现象称之为“transition reversal”。然而，Rosenboom & Hein^[9]和Lei & Zhong^[10]在计算N值时并没有发现这种现象，且钝度不断增加后，N值越小。Balacumar^[11]和Kara et al^[12]研究不同钝度的钝板、钝楔和钝锥的慢声波感受性时，得到的感受性系数亦是随着钝度的增加而始终减小。Fedorov^[13]认为这种“transition reversal”的现象可能与熵层不稳定性有关。究竟是否存在“transition reversal”的现象是值得研究的。

本文将利用直接数值模拟计算马赫数为4.5、前缘为圆头的0°迎角钝板，并在此基础上进行线性稳定性分析，研究超声速绕钝板的熵层不稳定性。

1 控制方程和计算方法 1.1 控制方程

0°迎角超声速绕钝板流动是二维流动问题，直角坐标系(x, y)下二维可压缩流动N-S的无量纲守恒型方程为：

$ \frac{{\partial U}}{{\partial t}} + \frac{{\partial E}}{{\partial x}} + \frac{{\partial F}}{{\partial y}} + \frac{{\partial {E_v}}}{{\partial x}} + \frac{{\partial {F_v}}}{{\partial y}} = 0 $

(1)

其中，$\;\;\;\;\;\;\;U = {[\rho, \rho u, \rho v, \rho {e_s}]^{\rm{T}}}, $

总内能$\;\;\;\;\;\;\;{e_s} = \frac{p}{{\rho \left( {\gamma - 1} \right)}} + \frac{1}{2}({u^2} + {v^2}), $

状态方程$\;\;\;p = \frac{{\rho T}}{{\gamma M_\infty ^2}}{\rm{。}}$

ρ、u、v、T、p分别为无量纲的密度、流向速度、法向速度、温度和压强，无量纲特征尺度分别为ρ_∞、u_∞、u_∞、T_∞和ρ_∞u_∞²，下标“∞”表示来流参数。长度的无量纲尺度为圆头半径r。M_∞表示来流马赫数，γ表示比热比。E和F表示对流项，E_v和F_v表示粘性项。图 1为坐标系和计算模型示意图，计算时采用贴体坐标系(ξ, η)，可对式(1)进行坐标变换得到(ξ, η)下的N-S方程，形式如下：

$ \frac{{\partial \bar U}}{{\partial t}} + \frac{{\partial \bar E}}{{\partial \xi }} + \frac{{\partial \bar F}}{{\partial \eta }} = \frac{{\partial {{\bar E}_v}}}{{\partial \xi }} + \frac{{\partial {{\bar F}_v}}}{{\partial \eta }} + \bar M $

(2)

其中，

$ \begin{array}{l} \bar U = JU, \bar M = JM, \bar E = J({\xi _x}E + {\xi _y}F), \\ {{\bar E}_v} = J({\xi _x}{E_v} + {\xi _y}{F_v}), \bar F = J({\eta _x}E + {\eta _y}F), \\ {{\bar F}_v} = J({\eta _x}{E_v} + {\eta _y}{F_v}), J = \partial \left( {\xi, \eta } \right)/\partial \left( {x, y} \right){\rm{。}} \end{array} $

1.2 计算参数和边界条件

本文选择马赫数为4.5、前缘为圆头的0°迎角钝板，可参考图 1的示意图。气体为理想气体，黏度和热传导系数均采用Sutherland公式，流动条件见表 1。

对于直接数值模拟，计算模型可分为基本流计算和扰动演化。不同的计算模型边界条件也将不同。在对称面处，两个模型均采用对称边界条件；在出口处，各变量由上游三个点进行二次外推；在壁面处，基本流计算时采用无滑移绝热壁，扰动演化时速度脉动为0，温度脉动的法向导数为0；在远场，基本流计算时上边界等于来流值，扰动演化时上边界各变量的脉动值均为0。

1.3 数值方法

对于直接数值模拟，基本流计算模型切换到扰动演化计算模型时，为了避免切换过程中引入数值扰动，需要保证通量分裂方法和时间离散格式相同，以及边界层和熵层区域的空间离散格式相同。在这两个计算模型下，通量分裂方法均选择AUSM plus分裂方法，时间离散格式均选择三步Runge-Kutta格式。但对于空间离散，控制方程的对流项所选择的计算格式将不同。计算基本流时，对流项在激波附近采用五阶WENO格式(Jiang & Shu^[14])，其他区域包括边界层和熵层区域采用五阶迎风格式；而进行扰动演化时，对流项统一采用五阶迎风格式。

本文的计算无量纲长度为1500，即7.5 m。在进行基本流计算时，计算域的流向网格和法向网格均为不均匀网格。对于法向网格分布，壁面附近网格较密，沿法向向上逐渐变稀，并且保证边界层区域内约包含101个法向点，熵层中广义拐点(广义拐点下文将介绍)附近的区域内约包含31个法向点。对于流向网格分布，头部区域较密，沿下游逐渐变稀。这样计算域的网格点数为1001×301。在进行扰动演化时，法向网格分布与上述一致，但流向网格为均匀分布，并保证一个扰动波波长大小的区域内至少包含21个流向点。这样计算域的网格点数为3001×301。数值计算程序由本课题组提供，其正确性可参见文献[15]。

2 计算结果 2.1 熵层和边界层

图 2给出直接数值模拟得到的温度云图，同时给出激波、熵层外缘和边界层外缘的位置曲线，图 2(a)给出的是流向区域为x=0~100的云图，图 2(b)给出的是流向区域为x=100~1000的云图，其中区域为x=400~430的云图是红色矩形的放大图，粉色曲线表示激波位置，黄色曲线表示熵层外缘位置，黑色实线表示边界层外缘位置。边界层外缘与壁面之间为边界层区域，边界层外缘与熵层外缘之间为熵层区域。

边界层外缘位置可通过总焓h₀参数确定。根据Kimmel & Klein^[16]的经验，总焓由壁面开始过峰值后等于来流总焓的1.005倍，则为该流向位置的边界层外缘位置。例如，在x=200处，如图 3(a)所示，给出总焓h₀沿垂直壁面法向y的变化，黑点即为边界层外缘位置。对于熵层外缘的位置，可根据熵变的大小来确定，Laible^[17]选择的经验方法是基于当地激波后熵变Δs的大小来判断，这样在确定每个位置的熵层外缘时需计算该位置激波后的熵变，增加了工作量。本文采用新的简单的经验方法，是基于某一固定值来确定。即，熵层外缘满足沿法向的熵变Δs等于驻点处熵变的0.01倍。某一位置的无量纲熵变Δs为：

$ \Delta s = \frac{\gamma }{{\gamma - 1}}{\rm{ln}}\left( {\frac{T}{{{T_\infty }}}} \right) - {\rm{ln}}\left( {\frac{p}{{{p_\infty }}}} \right) $

图 3(b)中给出熵变Δs沿垂直壁面法向y的变化，并根据上述方法确定出熵层外缘的位置。如此钝板流场中的熵层区域和边界层区域也就可以确定了，如图 2所示。需要说明的是，在图 2的x=0~1000流向范围内，没有明显看出熵层外缘与边界层外缘相交的趋势，换句话说，熵层外缘与边界层外缘夹角很小，熵层区域的特征仍然显著，熵层还没被“熵吞”。在这段区域内，熵层外缘与边界层外缘近乎平行，熵层流场在下游满足平行流假设。

2.2 熵层不稳定性

图 4给出x=200的流向速度和ρ∂u/∂y曲线，可以看出ρ∂u/∂y曲线出现两个极大值，故存在两个广义拐点，分别位于边界层和熵层，其中熵层区域的广义拐点记为GIP-Entropy。根据无粘不稳定理论^[18]，除了边界层中的不稳定模态，在基本流场中的熵层区域也有无粘不稳定模态。

为验证理论，针对熵层的广义拐点，进行线性稳定性理论(LST)分析。线性稳定性理论能够描述扰动的参数特征。它基于小扰动、平行流假设，忽略扰动方程的非线性项和基本量沿流向的变化，小扰动q′形式如下：

$ \begin{array}{l} q\prime \left( {x, y, z, t} \right) = {A_0}\hat q\left( y \right){{\rm{e}}^{{\rm{i}}\left( {ax + \beta z - \omega t} \right)}} + {\rm{c}}.{\rm{c}}.\\ q = \left\{ {\rho, u, v, T, p} \right\} \end{array} $

(3)

其中A₀表示扰动波的初始幅值，ω表示圆频率，α和β表示流向波数和展向波数，c.c.表示共轭项。将式(3)代入简化的扰动方程，结合齐次边界条件，得到扰动波的色散关系式，即f(α, β, ω, Re)=0，以及${\hat q}$的解。对于空间模式的扰动，ω为实数，α和β为复数，Re表示雷诺数。${\hat q}$也是复数，形式为$\hat q = {{\hat q}_r} + {{\hat q}_i}$，其模为$\left| {\hat q} \right| = \sqrt {\hat q_r^2 + \hat q_i^2} $。LST分析程序由本课题组提供，其正确性可参见文献[15]。

图 5给出二维扰动的空间增长率-α_i及相速度随圆频率的变化，可以看出中性解的上支界点的无量纲圆频率约为0.13，对应无量纲相速度约为0.92，与图 4中广义拐点GIP-Entropy对应位置的流向速度相等，这表明该模态为熵层引起的不稳定模态。图 6给出该不稳定模态的特征函数幅值分布，可以发现特征函数幅值的峰值对应的无量纲法向位置约为5.16，与图 4中熵层中广义拐点的位置相等，进一步证明了熵层中存在不稳定模态。

图 7给出x=200处不同圆频率下无量纲流向速度和温度的特征函数幅值分布，七个圆频率表示在图 5中。图 7(a)可以看出，较小圆频率时流向速度特征函数峰值在熵层广义拐点上，随着圆频率的增加，该峰值幅值逐渐减小，同时更靠近壁面的位置的幅值开始增加，并超过原峰值。可以发现，此时的峰值位置是图 4中Boundary-layer edge和GIP-Entropy之间的极小值对应的位置。而图 7(b)可以看出，不同频率下温度特征函数幅值的峰值位置均在熵层广义拐点GIP-Entropy上。

图 8和图 9给出x=200的中性曲线及β_r=0、0.04、0.08、0.12时空间增长率随圆频率的变化，可以看出熵层存在着三维不稳定波，但二维不稳定波的空间增长率最大，因此在研究熵层扰动对边界层的影响时，选择熵层二维不稳定波。

2.3 熵层扰动的演化

图 10给出熵层二维中性曲线沿流向的变化，可以看出，中性曲线是半开放半封闭的。相比于边界层内中性曲线的左封闭右开放的形状，熵层的形状则相反，即左开放右封闭。换句话说，熵层扰动在上游是不稳定的，而到下游则是稳定的。图 11给出无量纲圆频率为0.07的空间增长率沿流向的变化，可以看出空间增长率沿流向单调递减，且约在x=1000之后空间增长率为负，扰动开始衰减。另外注意到增长率为万分之一的数量级，扰动幅值变化很小。

为研究熵层扰动沿流向的演化过程，在x=400处加入通过稳定性分析得到的熵层扰动，进行扰动的数值模拟。图 12给出密度脉动的云图，可以看出扰动主要分布于熵层区域，边界层内扰动未见到明显变化。

图 13给出数值模拟的速度脉动的幅值曲线沿流向的变化，同时还给出线性PSE和LST的结果。可以发现约在x=1000之前，扰动幅值缓慢上升，在临界位置之后缓慢下降。而且线性PSE(线性PSE的程序验证参见文献[19])和LST的结果也与数值结果吻合得很好，由于LST对方程组作了平行流假设，说明了熵层内流动在下游具有很好的平行性。

3 结论

本文以前缘为圆头的0°迎角钝板为模型，通过直接数值模拟(DNS)计算基本流场，在此基础上利用线性稳定性(LST)分析研究了熵层中的不稳定性，同时根据经验方法确定了熵层区域和边界层区域，得到结论如下：

1) 熵层中存在广义拐点，通过稳定性分析表明熵层中存在不稳定扰动。另外，熵层存在三维不稳定波，但二维不稳定波增长率最大。

2) 扰动圆频率较小时流向速度特征函数峰值在熵层广义拐点上，随着圆频率的增加，该位置的峰值大小逐渐减小，同时边界层外缘与熵层广义拐点之间的某一位置的幅值开始增加，并超过原峰值大小。

3) 熵层不稳定模态集中在前缘附近，且空间增长率很小。利用直接数值方法模拟熵层扰动的演化过程，扰动在前缘附近幅值缓慢上升，之后缓慢下降。在距离前缘较远的范围内DNS结果与LST结果吻合得很好，表明熵层流场在下游具有很好的平行性。

需要指出的是，究竟是否存在“transition reversal”的现象还需要进一步研究。

致谢:

本文所使用的DNS程序由赵磊博士提供，PSE程序由张存波博士提供，在这里特别感谢。