0 引言

高超声速飞行器进气道唇口入射激波与机体肩部凸拐角流动的复杂干扰包含膨胀波、激波以及边界层分离等诸多流动现象，是目前空气动力学领域的前沿课题，亟待研究^[1]。流体在经过凸拐角壁面扩张产生的顺压梯度作用之后，会经历流向速度增大、边界层变薄(拐角较大时)等膨胀松弛过程^[2]，边界层的湍流特性因此发生变化，甚至出现再层流化现象^[3-5]。与此同时，激波与边界层相互作用会在壁面产生逆压梯度促使流动分离，从而增强流场湍流强度。因此，如果激波入射到凸拐角下游，流动会先经过扩张角产生的顺压梯度作用，随后在激波入射处受到逆压梯度作用。上述因素在多数激波/边界层干扰实验研究中并未被充分考虑，导致所获得的经验公式在预测以上问题时结果出现误差^[6]。为此，Zheltovodov^[7]指出：考虑凸拐角膨胀波影响下的激波/边界层干扰对未来的工程实际问题研究具有重要意义。

早期对凸拐角壁面激波与边界层干扰问题的研究主要集中为实验研究。Chew^[8]通过实验研究了来流马赫数1.8和2.5、扩张偏转角为8°的凸拐角壁面斜激波入射流动，考察了不同激波入射角(4°、6°和8°)下的流场图像、壁面压力以及分离区尺度。Hawboldt^[9]和Chung^[10]等人研究了来流马赫数8时激波入射作用下的凸拐角层流边界层，发现当入射激波脚位于拐角下游2到4倍边界层厚度时，凸拐角对激波诱导的流动分离具有限制作用；同时也发现激波入射到拐角下游时，壁面拐角扩张引起的顺压梯度会降低激波与湍流边界层的干扰强度，相应的壁面压力脉动峰值更小。对于高超声速流动，需要提高雷诺数才能确保壁面流动为湍流边界层。White^[11]等人利用较大尺寸的模型对来流马赫数11.5时激波入射凸拐角问题做了大量壁面压力和热流的测量，发现当激波发生器达到10°时，边界层出现初始分离，当激波入射到拐角下游且距离拐角较近时，分离泡尺寸才会减小。以上研究表明，凸拐角附近的流动分离会受到激波入射位置和来流湍流条件等多种因素的影响。

事实上，合理的激波入射位置有助于实现对凸拐角附近激波与湍流边界层干扰诱导流动分离的有效控制，目前相关研究主要集中在数值模拟方面。为了考察更接近于进气道肩部流场速度，Zhang^[12]等人采用二维雷诺平均模拟方法研究了来流马赫数3.5、扩张拐角15°的激波入射凸拐角流动，对比了激波入射点位于拐角上游以及下游的两种情况，发现了四种激波与膨胀波的干扰模式，其中激波/激波/膨胀波、激波/膨胀波/激波两种干扰模式对壁面流动分离具有抑制效果。实际问题中，二维进气道前体壁面前缘通常会做钝化处理，这使得前体壁面边界层流动特性发生改变。为此，Li^[13]等人开展了二维雷诺平均模拟研究，针对不同的平板前缘钝化半径，研究了来流马赫数5、扩张拐角14°的激波入射凸拐角流动，发现前缘钝化导致边界层厚度增大，从而促使激波诱导流动分离；进一步，他们考虑了黏性效应下边界层的法向非均匀性，据此对激波入射位置和角度的预测进行修正。上述研究给出了时均意义上激波入射凸拐角壁面的流动特征。需要说明的是，激波/边界层干扰所产生的分离流动具有较强的非定常特征，而流经凸拐角的湍流边界层演化还需考虑流场的三维效应。因此，采用三维非定常数值模拟可获得更多流动细节，有利于开展对凸拐角附近激波与边界层干扰的精细研究。

本文采用大涡模拟方法，对凸拐角附近激波与边界层干扰问题进行三维非定常数值模拟研究。通过对比考察尖拐角和钝化拐角两种情况，分析了在凸拐角附近激波/边界层干扰的流场结构、非定常分离流以及湍流统计特性，从而深化对此类流动问题的机理认识，以期对工程设计起到指导作用。

1 计算描述和程序验证 1.1 计算描述

如图 1所示，本文所研究的模型中，凸拐角偏转角度为θ_c=14°。来流马赫数Ma=3.0, 总温T₀=810 K，总压P₀=1.27 MPa。根据以上来流条件的单位雷诺数Re=4.7×10⁶ m^-1。做为对比，本文分别考察尖拐角和钝拐角(倒圆半径为20 mm)两种模型。如图 1所示，本文将凸拐角肩点的流向坐标设为x= 0。

本文数值模拟求解的控制方程为Favre滤波后的三维可压缩N-S方程^[14-16]，其形式如下：

$ \frac{{\partial \bar \rho }}{{\partial t}} + \frac{{\partial (\bar \rho {{\tilde u}_i})}}{{\partial {x_i}}} = 0 $

(1)

$ \frac{{\partial (\bar \rho {{\tilde u}_i})}}{{\partial t}} + {\rm{ }}\frac{{\partial (\bar \rho {{\tilde u}_i}{{\tilde u}_j})}}{{\partial {x_j}}}{\rm{ }} + {\rm{ }}\frac{{\partial \bar \rho }}{{\partial {x_i}}} - {\rm{ }}\frac{{\partial ({{\tilde \tau }_{ij}} - {\rm{ }}\tau _{ij}^{{\rm{SGS}}} + T_{ij}^{{\rm{SGS}}})}}{{\partial {x_j}}} = 0 $

(2)

$ \begin{array}{l} \frac{{\partial \left( {\bar \rho \bar E} \right)}}{{\partial t}} + {\rm{ }}\frac{{\partial [\left( {\bar \rho \bar E + \bar \rho } \right){{\tilde u}_i}]}}{{\partial {x_i}}} - {\rm{ }}\\ \frac{{\partial ( - {{\tilde q}_i} + {{\tilde u}_j}\tilde \tau {_{ij}} + J_i^{{\rm{SGS}}} - Q_i^{{\rm{SGS}}} + \sigma _i^{{\rm{SGS}}} - H_i^{{\rm{SGS}}})}}{{\partial {x_i}}} = 0 \end{array} $

(3)

其中:|f和$ \tilde f$皆为滤波瞬时量，分别表示对变量f进行空间滤波和密度加权Favre滤波后得到的物理量；上标“SGS”表示亚格子(Subgrid Scale，SGS)项，具体形式参见作者前期发表文章^[14-16]。上述方程无量纲化采用来流参量，即密度ρ_∞、温度T_∞、声速a_∞、来流边界层厚度δ_i，由此可得基于来流边界层厚度的雷诺数Re_δi=23500。相应地，基于边界层动量厚度和壁面摩擦速度的雷诺数分别为Re_θ= 2070、Re_τ= 510。

计算中边界条件设置如下：出口设为等压出口条件，壁面设为无滑移绝热壁面。计算中采用Rankine-Hugoniot关系引入斜激波，斜激波角β₁为24°。斜激波相对于凸拐角下游扩张壁面入射角α约为10°。为了研究激波入射到拐角下游、流动分离受拐角限制情况下的非定常流动特性，本文根据Zhang^[12]、Li^[13]等人的理论，将激波无黏入射点设为x/δ_i = 2.5。此外，为模拟湍流边界层情况，入口来流条件参数设置如表 1所示。本文所给入口条件参数设置是基于Pirozzoli^[17]、Li^[18]、Sandham^[19]等人发展的人工合成湍流边层(synthetic inlet forcing)方法，该方法已被成功用于可压缩湍流边界层数值模拟研究^[20]。

参考前人经验^[21]，本文计算域在流向、法向和展向分别取为L_x×L_y×L_z=45δ_i×10δ_i×3.2δ_i。本文重点关注凸拐角附近流场，因此对该区域网格进行局部加密(如图 2所示)，并通过法向拉伸以提高壁面附近网格分辨率，展向采用均匀网格。为了有效捕捉流动结构和非定常特征，三维计算网格数取为N_x×N_y×N_z=1793× 201×271，凸拐角附近网格分辨率为Δx⁺×Δy_min⁺×Δz⁺=25.4×0.6×12.1，计算时间步长取为a_∞Δt/δ_i=0.005。为提高计算效率，将计算域分为336个子块进行并行计算。为获得足够的流场信息，计算经过流场样本采集时间为80δ_i/a_∞，样本数为1000。本文中，〈〉表示基于时间和空间展向的平均；对于Favre滤波变量，{}表示密度加权平均，即|{$ \tilde f$ }=〈ρ $ \tilde f$ 〉/〈ρ 〉|。相应地，瞬时脉动量为|f^′=f －〈$ \tilde f$〉，f^″=$ \tilde f$ －{$ \tilde f$ }。

1.2 结果验证

在本文研究中，采用人工合成湍流边界层方法以使来流边界层迅速地发展为湍流边界层。图 3给出了流场充分发展后位于凸拐角上游x/δ_i=-10处van Driest变换后的平均速度分布。Smits^[22]等人关于可压缩湍流边界层的研究表明，经过van Driest变换后的平均速度应在黏性底层(viscous sublayer)满足线性律(linear law)分布U_VD⁺=y⁺；在对数律区(logarithmic layer)则满足对数律(log-law)U_VD⁺=5.1+2.44lny⁺。图 3所示的平均流向速度很好地符合了线性律及对数律分布。

通过湍流边界层瞬时结构如不同法向平面的流向速度条带，可考察来流湍流边界层的发展程度。如图 4所示，充分发展的湍流边界层内存在沿流向延伸的长条带结构，这种条带结构直接反映了高速流体与低速流体沿展向交替出现的现象，以及湍流边界层流动的间歇特性。这与前人的研究结果是一致的^[20]。

此外，展向空间关联函数在边界层内的分布可以验证展向计算域尺度的合理性^[20-21]。图 5给出了展向空间关联函数在法向四个位置的分布。如图所示，随着展向空间尺度增加，密度、流向速度、法向速度、以及展向速度的展向空间关联函数在不到半个展向计算域的距离就已衰减趋为0。可见，本文所选取的展向计算域尺寸是足够大的，且在展向方向采用周期性边界条件是合理的。

2 结果与讨论 2.1 平均流场特征

凸拐角附近的激波与边界层干扰问题包含激波干扰、流动分离等多种流动现象，这些流动现象可以通过时间和展向平均得到的流场显示，其中激波结构可以通过压力梯度幅值辨识。如图 6中的蓝色背景云图所示，超声速流动经过拐角之后，形成膨胀波系(expansion waves)，而入射激波(incident shock)与厚度为δ_i的边界层干扰之后，产生一道反射激波(reflected shock)。图 6中的紫色流线显示，湍流边界层在受到入射激波干扰作用后发生流动分离，在拐角后缘出现分离泡。分离泡在壁面起始处为分离点，结束处为再附点。同时，通过壁面附近黑色的涡量幅值云图分布可以看出，激波脚与分离泡之间出现了清晰的剪切层。

图 6(a)显示，在尖拐角情况下，膨胀波系分布较为集中，剪切层从拐角后缘平滑脱泻，分离泡比较饱满。相对地，在钝拐角情况下，图 6(b)中膨胀波系更为分散，分离泡高度相对平滑，剪切层与壁面之间的距离较小。

流场结构与壁面力学特征密切相关。图 7(a)-(c)依次给出壁面压力〈p_w〉、流向压力梯度d〈p_w〉/dx以及壁面摩擦系数C_f沿流向的分布曲线，其中实线和虚线分别对应为尖拐角和钝拐角的情况。如图 7(a)所示，在尖拐角情况下，由于拐角附近分布有局部集中的膨胀波系，壁面压力在拐角出现短暂下降，但在随后的激波入射作用之下，压力沿流向迅速上升。相应地，图 7(b)中流向压力梯度在尖拐角的上游壁面为负值(即顺压梯度，d〈p_w〉/dx < 0)，并在拐角处达到最小值；而在尖拐角下游为正值(即逆压梯度，d〈p_w〉/dx>0)。对于钝拐角情况，图 7(a)中壁面压力在拐角上游就开始下降，随后则沿流向逐渐上升。相应地，如图 7(b)所示，壁面顺压梯度值也只是位于钝拐角上游，这与图 6(b)所示在钝拐角上游出现膨胀波是相符的。

在超声速流动中，剧烈的逆压梯度分布可导致边界层发生分离，分离点前后壁面摩擦系数符号发生改变^[23]。因此，通过计算壁面摩擦系数分布可确定边界层分离情况。如图 7(c)所示，在尖拐角情况下，C_f首先沿流向逐渐上升，并在拐角处达到局部极大值，然后陡降为负值。而在钝拐角情况下，C_f的局部极大值出现在拐角上游。拐角下游的两个逆压梯度峰值之间即为分离泡区域，据图 7(b)可得，尖拐角情况下的分离泡流向尺寸约为1.52δ_i，钝拐角情况下约为1.33δ_i。根据前人研究^[24]，对于可压缩绝热壁流动，在本文所给的参数下，壁面摩擦系数理论值为C_f= 1.98×10^-3，而图 7(c)中显示上游壁面摩擦系数C_f= 1.82×10^-3，两者符合良好，这也说明来流湍流边界层已经充分发展。

2.2 瞬时流场特征

图 6的平均流场显示入射激波清晰可见，反射激波则较为模糊，表明激波经凸拐角区域干扰作用后非定常特性明显增强。图 8(a)和图 8(b)分别显示了在尖拐角和钝拐角两种情况下，展向中心截面瞬时密度云图呈现“皱纹状”分布^[25]。边界层流动的间歇结构与壁面斜交并向下游运动，且经入射激波作用之后，凸拐角下游的“皱纹状”脉动结构能够一直维持。

图 9给出基于Q准则的瞬态三维涡结构，该准则定义如下：

$ Q = \frac{1}{2}{\rm{ }}\left( {\mathit{\tilde \Omega }{_{ij}}\mathit{\tilde \Omega }{_{ij}} - \tilde S{_{ij}}\tilde S{_{ij}}} \right) $

(4)

这里$\tilde S{_{ij}} = \frac{1}{2}\left( {\frac{{\partial \tilde u{_i}}}{{\partial {x_j}}} + {\rm{ }}\frac{{\partial {{\tilde u}_j}}}{{\partial {x_i}}}} \right) $表示速度梯度张量的对称部分，$ \mathit{\tilde \Omega }{_{ij}} = \frac{1}{2}\left( {\frac{{\partial \tilde u{_i}}}{{\partial {x_j}}} - \frac{{\partial {{\tilde u}_j}}}{{\partial {x_i}}}} \right)$表示速度梯度张量的反对称部分。对于Q>0的区域，流体的转动速率大于应变速率，此时该区域流动涡结构占主导。从图 9可以看出，在来流湍流边界层内区存在大量不规则的发卡涡(马蹄涡)结构，而在边界层外区出现了“茎状”结构，同时有较规则的发卡涡结构，这些流动结构混合在一起形成湍流边界层中被称为“发卡包裹”的大尺度涡结构^{[20, 26]}。当经过拐角附近时，位于边界层外区的大尺度涡结构未发生明显减弱，而位于边界层内区的小尺度涡结构则明显受到顺压梯度影响发生变化。其中，在尖拐角情况下，出现在分离点附近的小尺度涡结构较少，这与尖拐角附近出现相对较强的顺压梯度有关，而在再附点附近，由于经过入射激波诱导的逆压梯度作用，小尺度涡结构迅速增强。在钝拐角情况下，拐角附近的顺压梯度较小，小尺度涡结构受其影响相对较小，因而在分离区仍存在大量小尺度涡结构。

图 10和图 11所示的瞬时涡量幅值云图证实了凸拐角附近大、小尺度涡结构的上述变化情况。如图 10所示，流动绕过凸拐角后，边界层外区的大尺度涡结构并未明显减弱，这表明凸拐角产生的顺压梯度对大尺度涡结构的影响并不明显。如图 11所示，与钝拐角情况相比，在尖拐角情况下分离点附近的近壁小尺度涡结构明显更少，这是由于流动绕过尖拐角产生更强的顺压梯度，因而近壁小尺度涡结构受到更强的抑制作用。

2.3 流动非定常特性

激波与湍流边界层干扰流动问题中，涉及到入射激波、反射激波以及分离泡的复杂非定常运动^[27]。通过计算流场瞬时数值纹影图(密度梯度幅值| ∇ρ|)，图 12给出了尖拐角情况下典型时刻展向中心截面内的数值纹影图，时间间隔为1.25δ_i/a_∞。从图中可以看出，激波入射引起流动分离，使得拐角之后的流动偏折加剧，从而在膨胀波后形成压缩波，并最终形成反射激波。对于入射激波来说，来流的顺时针旋涡结构会使其偏折，导致激波角增大，这一现象与Pirozzoli^[21]等人研究的斜激波入射平板边界层流动类似。另一方面，尖拐角后缘的剪切层在遇到入射激波之前结构明显，这与图 9所示的在这一区域小尺度涡结构急剧减弱有关，而剪切层在受入射激波干扰之后迅速失稳，致使波后流场小尺度涡结构增强。

图 13给出了钝拐角情况下的展向中心截面内的数值纹影图。与尖拐角情况不同的是，在遇到入射激波之前，小尺度涡结构受钝拐角顺压梯度影响较小，剪切层能够持续失稳脱涡，从而进一步诱导入射激波脚发生小幅度的非定常运动。

图 14给出了尖拐角情况下的典型时刻瞬时壁面摩擦系数，时间间隔为1.25δ_i/a_∞。如图所示，再附点位置随时间前后移动，而分离点位置基本保持不变。尖拐角处壁面摩擦系数较高，此时逆压梯度的作用较小，不足以使流动分离点越过尖拐角。

钝拐角情况下的瞬时壁面摩擦系数如图 15所示。相对于尖拐角的情况，钝拐角处的壁面摩擦系数明显降低。特别地，流动分离点不再固定出现在拐角处，而是在拐角附近前后移动，同时再附点也发生前后移动。

壁面压力信号的频谱分析有助于进一步理解流场的非定常特性。这里表征频率的无量纲量为St=fa_∞/δ_i。由于功率谱密度函数G(f)描述了脉动压力均方根σ在频域区间的分布特性，我们可以采用fG(f)/σ²来分析壁面压力频谱的分布特性^[28]。图 16和图 17显示了尖拐角和钝拐角情况下在流动分离点与再附点附近的壁面压力频谱。从图中可以看出，分离泡区域流动呈现出明显的高频脉动特性，对应的频率与上游湍流边界层的特征频率相近^[26]，这表明流场中湍流仍保持充分发展状态，并未因凸拐角作用而发生再层流化。

在尖拐角情况下，图 16(a)中的压力频谱分布表明，分离点附近的压力脉动主要为高频脉动，而在再附点附近高频压力脉动相对减弱，低频压力脉动有所增强，在St ~ 0.5处出现局部离散峰值(见图 16(b))。以上压力频谱分布特性与下述事实有关：剪切层在分离点附近尚未明显失稳，而在再附点区域则已因激波作用而失稳。在钝拐角情况下，分离点和再附点附近St=0.15 ~ 0.8的压力低频脉动特征较为明显(见图 17(a，b))，其原因是剪切层在分离点区域已失稳形成涡结构。

2.4 流场统计特性

图 18给出了压力脉动均方根p_rms^′在壁面沿流向分布的情况。在尖拐角情况下，压力脉动在拐角处出现局部峰值，随后沿流向略有下降，但在再附点附近再次升高。Gerolymos^[29]等人的研究表明，在可压缩流动中，除了激波运动之外，流场压力梯度与速度脉动之间的相互作用也是产生压力脉动的原因之一。尽管在尖拐角处分离激波没有明显运动，但是该处存在较强的顺压梯度，并与来自上游湍流边界层的速度脉动相互作用。因此，在尖拐角附近出现了较为明显的脉动压力峰值。

在钝拐角情况下，壁面压力脉动沿流向出现二次抬升，即在分离泡前端抬升一次，随后在再附区域再次抬升。类似现象在斜激波入射平板边界层问题中也有发现。上述壁面压力脉动分布特性与常规的激波与边界层相互干扰问题类似^[23]。

湍动能可以用来表征流体脉动运动的强弱。图 19给出了以来流参数无量纲化的湍动能(Turbulent Kinetic Energy, TKE)分布。图 19(a)显示，流动经过尖拐角后立即产生分离，导致流体脉动运动增强。一方面，来流湍流边界层大尺度结构与激波相互作用通常会引起激波的“褶皱”^[30]，从而产生激波非定常运动。另一方面，拐角后缘的剪切层失稳强度由小尺度涡结构决定，且与拐角附近流动沿流向从顺压梯度转变为逆压梯度密切相关。从图 19(a，b)可以看出，与尖拐角情况相比，钝拐角情况下湍动能在分离泡周围的分布显著不同。在尖拐角情况下，流动经过拐角先受顺压梯度作用导致小尺度涡结构衰减，后受逆压梯度作用导致剪切层迅速失稳。相对而言，流动经历了更为急剧的压力梯度转变过程，剪切层脉动相对更为强烈，因而剪切层区域湍动能较强。而钝拐角情况下的顺压梯度较弱，对小尺度涡结构影响相对较小，因此流场的局部脉动运动主要出现在小尺度涡结构较强的激波脚和分离泡区域。

通常在凸拐角边界层问题中，流动经过拐角之后湍动能减弱。例如，Sun^[2]等人在研究扩张角为4°的凸拐角边界层流动发现，经过拐角之后湍动能强度减弱约50%，并在拐角下游约15倍边界层厚度处仍未恢复至其拐角上游的值。而在本文所研究的斜激波入射下的凸拐角流动中，如图 20所示，在经过拐角后湍动能仅在靠近壁面的区域有所下降，在距壁面较远处因激波作用反而升高。可见，入射激波作用是经拐角后湍流强度流向变化的主导因素。

3 结束语

本文采用大涡模拟方法研究了来流马赫数Ma_∞=3、雷诺数为Re_θ=2070的凸拐角附近湍流边界层与入射斜激波干扰问题。主要考察了尖拐角和钝拐角两种情况下的基本流场结构、非定常分离流以及湍流统计特性。本文主要结论如下：

1) 流场存在凸拐角引起的顺压梯度和入射激波引起的逆压梯度。这两种压力梯度的耦合作用使得流场出现激波、膨胀波、分离泡以及剪切层等复杂流动结构。在尖拐角情况下顺压梯度较为局部集中，导致壁面摩擦系数在拐角处出现峰值，分离点因而难以越过拐角向上游运动；而在钝拐角情况下，分离点则能够越过拐角出现在拐角上游。

2) 流动显示表明，流体经过尖拐角先受顺压梯度作用导致小尺度涡结构减弱，后受逆压梯度作用导致剪切层迅速失稳；而在钝拐角情况下流场受压力梯度影响较弱，上游生成的小尺度涡结构能够一直维持。瞬时壁面摩擦系数和频谱分析表明，在尖拐角情况下，分离泡的运动主要受再附点运动的影响；而在钝拐角情况下，分离点和再附点都会前后移动。

3) 壁面压力脉动均方根分析表明，尖拐角附近的压力脉动峰值与压力梯度与速度脉动之间的相互作用有关，而钝拐角压力脉动则在再附点附近达到峰值。流场脉动特性分析表明，激波与分离区干扰使湍流强度沿流向迅速增强，湍流脉动剧烈的区域主要集中在激波脚、剪切层以及再附点附近。