圆锥激波与平板边界层干扰的数值模拟研究

杨柳青，李 沁，尤延铖

（厦门大学航空航天学院，厦门 361102）

在超声速、高超声速飞行器上的内外流中，广泛存在因激波/边界层干扰（Shock-wave/boundary layer interaction，SWBLI）而引起的复杂分离流动现象。SWBLI 所产生的逆压梯度会引起流动分离和总压损失，增大流动的不稳定性，导致进气未启动等情况。这些流动现象将会给飞行器的飞行性能、表面热防护、结构设计等方面带来直接影响。

近年来，人们对三维激波/边界层相互干扰有了相当程度的了解，但目前的研究大多集中于后掠激波与边界层的相互干扰［1-5］。锥形激波/边界层的相互干扰（Concial shock-wave/boundary layer interaction，CSWBLI）在高速飞行器上也很常见，例如助推火箭顶部整流罩产生的锥形激波与火箭主体表面的边界层相互作用，但是这方面的研究相对较少。Migotsky 等［6］首先对CSWBLI 进行了数学分析；Panov［7］开展了锥形激波与平板边界层干扰的实验研究；Kussoy 等［8］在来流马赫数为2.2 的情况下进行了圆锥激波与轴对称边界层的实验研究；Gai 等［9］通过实验研究了来流马赫数为2 时，锥形激波与平面湍流边界层的相互作用；Hale［10］研究了来流马赫数为2.05 时锥形激波与平面湍流边界层的相互干扰情况；Zuo 等［11］模拟了Hale［10］的实验模型和流场条件，进行了高保真的直接数值模拟；之后，Zuo 等［12］又采用雷诺平均NS 方程（Reynolds averaged Navier-Stokes，RANS）的方法，数值模拟了在自由流马赫数4.0 条件下，由轴对称壁面发展而来的内锥形激波/边界层的相互干扰。这项工作旨在扩展先前Zuo 等［11］对外流CSBLI 的研究。利用直接数值模拟数据库［11］，Zuo［13］还进一步研究了CSWBLI 中分离泡的尺寸及其与压升的关系，建立了中等激波强度下壁面压升与分离泡几何形状的标度关系。

在已开展的针对CSWBLI 的研究工作中，干扰流场的实验结果缺乏足够的细节，更多雷诺数变化对干扰结果的影响有待研究；同时，缺乏三维计算和二维计算结果的对比分析以显示三维效应的影响。为了补充这方面的研究，本文采用两方程Menter-SST 模型，以Gai 等［9］所采用的实验模型和流场条件为基础，对圆锥激波/平板边界层干扰进行了数值模拟。通过定性与定量分析，进一步研究了入射激波与边界层交汇后干扰区的流动情况；同时还分别改变了激波发生器顶角和来流单位雷诺数以比较不同工况下的干扰结果，并探讨其流动分离的变化规律。此外还进行了二维情形的计算，将其与相应的三维结果进行了对比，分析两者间的差异，研究三维效应对数值模拟结果的影响。

1 研究对象及数值方法

1.1 数值方法

本文数值模拟使用风雷软件作为计算平台，控制方程为RANS 方程，其积分形式为

式中：Q为流动变量，n为控制体表面的外法向向量，FI和FV分别表示对流（无黏）通量和黏性通量。湍流模型选用两方程Menter-SST 模型，其具体表达式见文献［14］。时间离散格式采用Jameson［15］提出的双时间步隐式格式（Lower-upper symmetric-gauss-seidel，LU-SGS）。

1.2 研究对象

1.2.1 几何外形

本文的研究对象为Gai 等［9］所使用的由圆锥和平板组成的物理模型，其几何外形如图1 所示。

图1 几何外形示意图Fig.1 Geometrical sketch map

圆锥为激波发生器，圆锥半锥角θc可变，圆锥底面直径为20 mm；平板长279 mm，圆锥中心距平板高度h= 30 mm。为了方便比较，以平板前缘为原点，保证各计算条件下入射激波的入射点始终距离原点xc=169 mm。

1.2.2 网格及边界条件

为了方便生成网格并减少网格规模，对物面外形进行了简化处理：考虑到圆锥激波发生器上半部分对于激波/边界层干扰没有影响，因而生成的网格只包含圆锥下半部分；由于物面外形及其周围流场沿中心对称面z= 0 对称，因此仅生成对称面一侧的网格；激波发生器在风洞实验中都是有限长度的，而在简化网格中，将其底面一直延伸到了计算域的出口处，这样处理虽然会使激波发生器在原底面位置处形成的膨胀波系结束角度不同，但是由特征线理论［16］可知，不同转折角下产生的膨胀波在偏转同一角度的区域内的解相同，因此不会给激波/边界层干扰区带来额外影响。

计算域长宽高分别为190、55 和30 mm，起始位置（即来流入口）为平板前缘。根据圆锥高度和激波角，得到三维计算下圆锥顶点至来流入口的距离lc1，见表1。

表1 三维计算下的lc1Table 1 lc1 in three-dimensional cases

以θc=20°的圆锥（即算例4）为例，模型的计算域如图2 所示。

图2 算例4 的计算域示意图Fig.2 Sketch map of computational region for Case 4

计算网格采用结构网格，网格总体图见图3（a），网格总数约为540 万个。对变形较大的位置做了过渡处理且基本确保壁面位置的网格正交性，如图3（b，c）所示。平板上采用加密网格，平板的物面网格数为198×156（流向×展向），第一层网格厚度为0.001 5 mm（y+为2），网格增长比率为1.064 4。

图3 算例4 网格示意图Fig.3 Grid for Case 4

已知受三维效应的影响，入射激波角等效的二维流动与三维流动之间存在差异。为了进一步研究这些差异随流动参数改变的变化规律，本文还安排了与三维问题对应的二维算例，其入射激波角度与三维中心对称面上入射激波的角度相同。二维计算中，圆锥简化为斜楔。根据激波理论，对于同一激波角，三维计算和二维计算各自需要的顶角大小不同，激波发生器位置相应发生变化。此外，在实际计算过程中发现如果仍采用三维计算时的激波发生器高度，在二维计算时会出现壅塞现象，因此二维计算时需提高激波发生器的高度。根据后续测试，在二维情况下，斜楔高度设置为40 mm。根据斜楔高度和激波角，得到二维计算下相应的斜楔楔角θ2大小及其顶点与来流入口距离lc2，如表2所示。

表2 二维计算下的lc2Table 2 lc2 in two-dimensional cases

计算域入口采用超声速来流条件；出口为超声速出流，直接使用外推；平板壁面为无滑移、绝热边界条件。

1.2.3 计算工况

本文算例的入口来流状态如表3 所示。

表3 入口来流状态Table 3 States of incoming flow

2 计算结果与讨论

2.1 网格收敛性和数值验证

2.1.1 网格收敛性

使用算例4-2D，来流状态Run 2 来进行网格收敛性分析。参与比较的计算网格（流向网格数×纵向网格数）如表4 所示。

表4 计算网格Table 4 Computational grid

计算结果如图4 所示，图4 给出了3 种网格下的壁面压力分布和壁面摩阻系数分布，可以看到Grid 2 和Grid 3 计算得到的结果基本重合，而Grid 1 的结果与两者存在一些差距。根据网格收敛性分析，可以认为Grid 2 和Grid 3 的网格计算收敛，因此在后面计算中，中心对称面上网格采用Grid 2 密度的网格。

图4 网格收敛性分析Fig.4 Analysis of mesh convergence

2.1.2 数值验证

本文以文献［9］提供的实验数据作为数值验证的标准，图5 为流向截面壁面压力数值模拟结果与实验结果的比较，使用的是算例1，状态Run 2。

图5 中，实验结果在分离引起的压升位置上游存在明显的波动，推测该情况是实验测量过程中各种误差积累的结果。对比结果显示，使用风雷软件计算得到的数值模拟结果与实验测量值吻合得较好，壁面压力分布的最大误差控制在2%以内。随着逐渐远离中心对称面（z= 0 mm），各截面压升顶点逐渐降低，但每次降低的程度并不相等，越远离对称面，降低幅度越大；同时，各截面压升起始位置逐渐向下游移动，与壁面压升变化情况类似，不同展向位置截面向下游移动的程度并不一致，越远离对称面，压升起始位置移动幅度越大。压升峰值后的压降是因为圆锥底面产生的膨胀波传播到壁面导致的。

图5 算例1、状态Run 2 的流向截面壁面压力分布Fig.5 Wall pressure distribution in flow direction section for Case 1, Run 2

图6 为上游影响范围的数值模拟结果与实验结果的比较，图中“ref”指相应的实验数据，来自文 献［9］，参 与 比 较 的 算 例 有Case 1、Case 2 和Case 3，初始来流状态均为Run 2。上游影响范围定义为壁面压升的起始位置，可以表示壁面上的分离激波轨迹。从图中可以看出，除了θc= 14°时展向位置z= 24 mm 的数值模拟结果比实验结果有比较明显的上游偏移外，其余数值模拟结果与实验结果拟合得较好。由于边界层的存在，黏性计算下的激波轨迹相较于对应角度的无黏计算下的入射激波轨迹，其位置普遍向上游移动。由于远离对称面的激波效应减弱，所有黏性计算下的激波轨迹相较于无黏轨迹，越靠近出口，轨迹的曲率逐渐减小；对于4 个黏性计算下的算例，其激波轨迹向上游移动距离的变化幅度并不相等，14°、16°、18°的 变 化 幅 度 大 致 相 等，而20°相 较 于18°的结果，变化幅度并不大。此外，所有在中心对称面上的数值结果与实验结果相比，要略微向下游偏移。

图6 上游影响范围与无黏计算下的入射激波轨迹Fig.6 Upstream influence range and incident shock wave trajectory under inviscid calculation

2.2 不同半锥角对锥形激波边界层干扰的影响

计算不同半锥角时，来流状态均为Run 2。

2.2.1 中心对称面的壁面压力分布比较

图7 给出了不同θc角下二维、三维中心对称面（Z= 0）的壁面压力分布。

图7 不同θc 角下二维、三维中心对称面的壁面压力分布比较Fig.7 Comparison of wall pressure distribution between 2D cases and 3D cases in center-symmetric plane at different θc

由图7 中三维数值模拟结果间的对比表明，分离点随θc逐渐增大而向上游移动，原因是不同θc下产生的分离涡大小不同，导致各自的压升起始位置不同；压力比峰值也随着θc的逐渐提高而增大；此外，随着θc的增大，压升过程中逐渐显现出与分离流有关的压力平台，θc越大，压力平台越明显。二维数值模拟结果间的对比显示出与三维计算类似的规律。由于在二维计算的计算域内，圆锥底部产生的膨胀波还没有影响到壁面，因此图7 的二维计算结果没有表现出二次压升后的压降过程。

考虑到二维计算下没有横向流动，理论上二维计算下的壁面压升比应该会比对应的三维计算结果大，但是对比表明实际计算结果并不全是如此，需要将干扰区的压升分为一次压升和二次压升来分别讨论。对于一次压升幅度，除了θc=20°的计算结果外，其他3 个较小θc下的二维计算结果都接近对应的三维计算结果。一次压升是由分离激波引起的［17］，而分离激波的强度与分离涡尺寸密切相关，当θc较小时，三维计算受三维效应的影响并不显著，使得三维计算下中心对称面上的分离涡尺寸与对应的二维结果相差并不大（具体见2.2.2 节），因此两者的分离激波强度也相似，从而使得一次压升幅度大致相同；而当θc=20°时，由于受到的三维效应更加显著，二维计算的分离涡尺寸明显大于三维结果，因此其产生的分离激波更强，从而使得二维计算的一次压升幅度大于三维结果。综合以上分析，可以认为存在一个阈值，当θc大于等于该值时，中心对称面上会因为受到显著的三维效应影响而使得二维计算的一次压升幅度大于三维结果。压力平台后的二次压升则与再附过程有关（受到再附激波等因素的影响）［17］，从图7 可知，4 组计算在三维计算下的二次压升增长率都要高于对应的二维结果，而θc=20°时出现三维计算的二次压升幅度低于对应二维结果的情况，其原因是此时圆锥半锥角增大引起膨胀波向上游移动，使得再附过程更早地受到了膨胀波的影响。

2.2.2 中心对称面的旋涡结构和涡量分布比较

图8 给出了不同θc角下二维、三维中心对称面上的截面流线和展向涡量分布比较。以下从分离涡的尺寸、拓扑结构和涡量分布方面对图8 分别进行分析讨论。

图8 不同θc角下二维、三维中心对称面上的截面流线和展向涡量分布比较Fig.8 Comparison of streamline and vorticity between 2D cases and 3D cases in center-symmetric plane at different θc

（1）分离涡尺寸

三维和二维计算结果比较表明，由于保持不同θc下激波入射点位置不变，分离涡的涡核位置基本一致（在x= 166～167 mm），而分离点随θc的减小而逐渐向下游移动。结合表5 分析，分离涡的尺寸随θc增大而逐渐增大，其中宽度变化比较平均（二维计算中按2～3 倍增长；三维计算中每次增长1.5～2 mm），而分离涡的高度在θc=16°、18°、20°时是按照2～4 倍逐渐增大的，但是θc=16°下的分离涡高度相较于14°，增长了10 倍。文献［9］中，实验结果表明θc=14°时没有观察到分离现象，但数值模拟结果表明，分离还是存在的，只是尺度非常小。此外，图7 中的压力平台范围显然与分离涡的尺寸有直接关系，分离涡尺寸越大，压力平台就越明显。

表5 不同半锥角（楔角）的二维/三维分离涡尺寸对比Table 5 Comparison of separated vortex between 2D cases and 3D cases at different θc

考虑到气流的横向逸出，二维计算得到的分离涡应该要比对应的三维计算结果大，但是在目前的4 组算例中，只有θc=20°的结果符合预期。其他二维计算下的分离涡顶点大多位于相应三维计算下的分离流线和再附流线之间（或者接近再附流线顶点），而且分离点位置变化不大。

（2）分离涡内流线的拓扑结构

二维计算中，流线在分离涡内围绕涡核形成系列封闭流线，其边界是分离流线。而在三维流线中，分离涡不再是类似的封闭形状，而是开放形式；流线也不再是封闭的曲线，而是围绕涡核的螺旋形曲线；起于分离点的分离流线和滞止于再附点的再附流线是不同的流线，在图8 中可以看到两者间存在明显的间距。涡核处，流动向侧向逸出，这是一种典型的三维效应。

（3）涡量分布

二维结果和三维结果反映出同样的涡量大小分布规律，即分离涡内的涡量并不大，涡量较大的区域集中在边界层之中，边界层内涡量较大是因为边界层内受黏性作用形成剪切流。此外，由于边界层会受到分离涡的影响而发生抬升，因此对于不同θc的计算结果，分离涡附近的涡量分布与分离涡的尺寸有直接关系。而比较二维结果与相应的三维结果可以发现，两者间的差异与各自的分离涡大小有直接关系，因此可以认为涡量分布也会受到三维效应的影响。

2.2.3x=169 mm 位置的展向截面密度梯度分布比较

图9 给出了不同θc角下x=169 mm 的展向截面密度梯度云图。从图中可以分辨出较为清晰的来流边界层（黑框）、入射激波（蓝框）、反射激波（红框）等流动结构。截面上，密度梯度的大小随θc增大而增大；入射激波、反射激波等流动结构的尺寸也与θc密切相关，当θc增大时，截面上入射激波、反射激波、膨胀波系的尺寸均增大。

图9 不同θc 角下x = 169 mm 的展向截面密度梯度Fig.9 Density gradient of spanwise section (x = 169 mm)at different θc

2.3 不同雷诺数对锥形激波边界层干扰的影响

计算不同雷诺数时，使用的激波发生器外形均为Case 4。

2.3.1 中心对称面的壁面压力分布比较

图10 给出了不同Re数下二维、三维中心对称面的壁面压力分布。图10 显示三维计算下分离点随着Re增大而逐渐向下游移动，并且一次压升后的压力平台也变得越来越不明显，直至Re=3×109时，压力平台完全消失（说明该雷诺数下没有发生分离）。这表明入射激波与边界层相互作用产生的分离涡随Re增大而逐渐减小，一方面是因为边界层厚度减小，另一方面是因为高雷诺数下，边界层内速度亏损小，剖面含有更高的动能，对逆压梯度阻滞的抵抗能力更强。但是压升峰值随着Re增大而逐渐提高，且提高过程不完全是一个线性过程，Re= 3×106、3×107、3×108下的增长幅度大致相同，而Re= 3×109的压升峰值相较于Re= 3×108的结果，增长幅度很大。二维计算也显示出了与三维计算类似的规律，中心对称面的壁面压力变化情况总体上与分离涡的尺寸相匹配。

图10 不同Re 数下二维、三维中心对称面的壁面压力分布比较Fig.10 Comparison of wall pressure distribution between 2D cases and 3D cases in center-symmetric plane at different Re

比较三维计算结果和对应的二维计算结果发现，除了Re= 3×109时的计算外，其余3 个雷诺数下的二维计算结果相较于对应的三维计算结果，分离点都显著地向上游移动，且压力比基本上都大于三维结果。Re= 3×109时，不仅没有发生流动分离，而且二维计算下的壁面压升峰值也不及三维计算。Re= 3×106、3×107、3×108时三维计算与二维计算之间的差异可以用三维效应来解释，而Re= 3×109的计算结果则表明高雷诺数可能会减弱三维效应带来的影响。

2.3.2 中心对称面的旋涡结构和涡量分布比较

图11 给出了不同Re数下中心对称面的截面流线和展向涡量分布。以下从分离涡的尺寸、拓扑结构和涡量分布方面对图11 分别进行分析讨论。

（1）分离涡尺寸和拓扑结构

图11 的结果表明，随着Re增大，分离涡尺寸显著减小，分离点向下游移动，同时涡核位置基本保持不变（因为控制了激波的入射位置）；直至雷诺数增长至Re= 3×109时，入射激波与边界层的相互作用未发生分离。该规律对于二维计算和三维计算都成立，其原因与2.3.1 节中的分析相似。

图11 不同Re 数下中心对称面的截面流线和展向涡量分布比较Fig.11 Comparison of streamline and vorticity between 2D cases and 3D cases in center-symmetric plane at different Re

对比二维和三维的数值模拟结果，同样可以发现二维计算中分离涡内的流线是封闭的，而三维情况则是螺旋形曲线；同时结合表6，二维计算下的分离涡尺寸要大于相应的三维计算结果，这些都体现了三维效应的影响。但是进一步比较二维、三维分离涡尺寸的变化，注意到随着Re的增大，二维分离涡与三维分离涡之间的差距在逐渐缩小，直至Re= 3×109时都没有发生流动分离，流线几乎一致，这表明高雷诺数下中心对称面上二维计算与三维计算的流动结构差异减小。

表6 不同Re 数的二维/三维分离涡尺寸对比Table 6 Comparison of separated vortex between 2D cases and 3D cases at different Re

另外，Re= 3×109时，二维、三维计算都没有产生分离涡，考虑到高雷诺数对三维效应的抑制作用，是否存在一个雷诺数的阈值，可以使得三维计算时无分离涡，而二维计算时有分离涡，还需要进一步探索。

（2）涡量分布

与2.2.2 节的计算结果类似，二维、三维计算下涡量较大的区域均集中于边界层之中，而分离涡内的涡量并不大，且涡量分布与分离涡的尺寸（影响边界层状态）有直接关系，因此涡量分布也受到雷诺数变化和三维效应的影响。

2.3.3x=169 mm 位置的展向截面密度梯度分布比较

图12 给出了不同Re数下x= 169 mm 的展向截面密度梯度云图。从图12 可以看到较为清晰的来流边界层、入射激波、反射激波等流动结构。密度梯度大小随Re增大而增大；来流边界层、反射激波、膨胀波系的尺寸变化情况则与之相反，当Re增大时，这些流动结构的尺寸均明显减小。

图12 不同Re 数下x = 169 mm 的展向截面密度梯度Fig.12 Density gradient of spanwise section (x = 169 mm)at different Re

2.3.4 壁面极限流线比较

图13给出了不同Re数的壁面极限流线，极限流线图取第一层网格高度位置的水平截面进行绘制。

图13 不同Re 数的壁面极限流线Fig.13 Limiting streamlines of wall at different Re

除了Re= 3×109时的计算外，其他雷诺数下的计算都产生了反射激波诱导分离现象，相应的极限流线图上都存在相似的鞍点（红圈）、结点（蓝圈）、分离流线（红色虚线框）和再附流线（蓝色虚线框）等典型特征结构。Re= 3×109的计算结果则呈现典型的无分离弱干扰流动特征。对于Re=3×106、3×107、3×108的 结 果，随 着 雷 诺 数 增大，鞍点向下游移动，结点向上游移动，整体上分离流线和再附流线之间的间距在逐渐缩小，这一现象也反映出不同雷诺数下入射激波与边界层相互作用而产生的分离涡尺寸的变化。

2.3.5 壁面上游影响比较

图14 给出了不同Re数的上游影响范围和无黏计算下的入射激波轨迹。

图14 不同Re 数的上游影响范围与无黏计算下的入射激波轨迹Fig.14 Upstream influence and inviscid shock trace at different Re

与图6 显示的结果类似，黏性条件计算下的分离激波轨迹相较于无黏激波轨迹，轨迹位置均向上游移动，且越靠近出口位置，轨迹的曲率越小。随着Re增大，黏性计算下的激波轨迹向上游移动的距离逐渐减小，这与分离涡尺寸的变化情况相一致。

3 结 论

本文针对圆锥激波与平板边界层的相互干扰，使用两方程Menter-SST 模型进行了定常数值模拟。计算模拟了不同半锥角和雷诺数下的干扰流场，总结了不同参数下计算结果的变化规律，扩展了先前的实验和数值研究内容；此外，将三维计算结果与对应的二维计算结果进行了对比，分析了三维效应对计算结果的影响。

（1）随着半锥角的增大，入射激波与边界层干扰产生的分离涡尺寸增大，上游影响范围朝上游移动；壁面压升峰值提高，一次压升后的压力平台越发明显；干扰区的密度梯度大小增大，显示出的反射激波和膨胀波系等流动结构的尺寸也在扩大；中心对称面上的涡量分布则与分离涡尺寸直接相关。受三维效应的影响，中心对称面上的分离涡内流线的拓扑结构在二维计算中是封闭的，而在三维计算中则是开放的螺旋形曲线。此外，在分离涡尺寸和壁面一次压升方面，除了θc=20°的结果，其他3 个较小θc下的二维计算结果均接近对应的三维结果，因此可能存在一个阈值，当θc大于该值时，三维效应的影响才会明显显现出来。

（2）随着雷诺数的增大，入射激波与边界层干扰产生的分离涡尺寸减小，直至不发生流动分离，上游影响范围朝下游移动；壁面压升峰值提高，但是一次压升后的压力平台愈发不明显；干扰区的密度梯度大小增大，但是显示出的反射激波和膨胀波系等流动结构的尺寸减小；中心对称面上的涡量分布同样与分离涡尺寸直接相关。除了拓扑结构上的区别，受三维效应的影响，二维计算下的分离涡尺寸要大于相应的三维计算结果，但同时二维分离涡与三维分离涡之间的尺寸差距会随着Re的增大而逐渐缩小，这表明高雷诺数下中心对称面上两者流动结构的差异会减小。而在中心对称面的二次压升峰值方面，较小的3 个雷诺数下的二维结果大于三维结果，但是两者间的差距会随着Re的增大而逐渐缩小；当Re= 3×109时，三维情况下的二次压升峰值大于二维结果。