轴对称内锥流动中马赫盘的形成与演化

姬隽泽，李祝飞，*，司东现，张 涛，施崇广，尤延铖，杨基明

(1. 中国科学技术大学，合肥 230027；2. 厦门大学，厦门 361102)

0 引 言

近年来，内转式进气道[1-4]凭借诸多性能优势，在吸气式高超声速飞行器的研发中展示了诱人的潜力。然而，由于内转式进气道普遍是基于轴对称内收缩基准流场进行设计的，其内部流动存在汇聚效应[5-8]，与二元进气道相比更为复杂，甚至存在本质性差异。一般而言，二元进气道内部的激波相互作用通常表现为平面激波的入射和反射[9-10]，并存在规则反射和马赫反射两种反射类型[11]。然而，在轴对称内收缩流场中[12-14]，内锥形的入射激波会持续向中心汇聚增强，最终必然在轴线上发生马赫反射[15-16]。在工程实际中，马赫反射区域具有强压缩损失和弱抗反压能力等特性[12,17]，会给进气道性能带来严重的负面影响，因此应极力避免。

长久以来，通过在轴对称内收缩基准流场中设置中心体将轴线区域掩盖，进而规避马赫反射的措施[12,18]，在内转式进气道设计中被广泛采纳。实际上，这种措施在给进气道设计带来便利的同时，并未促进流动机理认识的提升，甚至制约了设计理论及方法的持续创新。因此，随着对内转式进气道流场品质的要求日益严苛，越来越多的研究者开始重视轴对称内锥流动中的激波反射问题，并尝试探索理论求解方法。谭廉华等[19-20]从轴对称内锥形激波和马赫盘入手，较早地给出了确定它们形状的理论关系式。近年来，弯曲激波理论[21-25]（Curved shock theory，简称CST）逐步得到广泛应用。该理论可以在已知内锥流动中入射激波和马赫盘形状的前提下，计算出反射激波的局部形状以及波后流动的高阶参数。将上述两种理论相结合，似乎可以逐次确定轴对称内锥流动中原本难以解析的激波形状，给理论求解轴对称内锥形激波反射问题带来了新的曙光。然而，应用这两种理论均涉及到马赫盘位置这一关键参数。由于目前缺乏对轴对称内锥流场中马赫盘形成机理及马赫盘下游流动的深入认识，从理论上确定马赫盘位置的研究仍鲜见报道。

根据平面激波反射理论可知[11]，当入射激波的角度β大于von Neumann准则规定的激波角βN时，有可能发生马赫反射。对于平面激波的马赫反射结构，入射激波、反射激波和马赫杆的波后气流偏转角关系（即三激波理论）是唯一确定的。由于平面激波的β保持恒定，三激波理论可以在平面激波的任意位置成立。于是，马赫杆的形成位置仅由其下游流动决定[11,26]。

与平面激波相比，轴对称内锥形激波随着向中心汇聚，自身的激波角β具有持续增大的特点，使得影响马赫盘位置的因素变得更为复杂。在轴对称内锥形激波的入射过程中，当β尚未增大至βN时，不能形成马赫反射结构。换言之，轴对称内锥形入射激波自身的汇聚增强过程会限制马赫盘可能存在的空间范围。进一步地，虽然三激波理论在β超出βN之后成立，但是三激波理论中的反射激波后的气流偏转角不是唯一确定的，而是与当地入射激波的β紧密相关。由此可见，在确定轴对称内锥形激波反射流场中马赫盘位置时，入射激波自身的强度变化需要予以考虑。然而，关于马赫盘位置与激波汇聚增强过程之间的关系，目前鲜见报道。

此外，与平面激波马赫反射类似的是，在轴对称内锥形激波反射流场中马赫盘下游流动也会影响马赫盘位置。Gounko[27]通过改变轴对称内锥模型的壁面长度，获得了不同的马赫盘下游流动条件及马赫盘位置，但并没有阐明马赫盘下游流动对其形成位置的影响机理。在下游流动影响下，马赫盘的形成及演化规律仍有待系统地研究。

综上所述，本文综合考虑轴对称内锥形激波汇聚过程及马赫盘下游流动两方面的因素，采用无黏数值模拟，考察内收缩直锥的壁面前缘角度和壁面长度对马赫盘位置的影响，揭示形成马赫盘的内在机制及其演化规律，以期丰富对轴对称基准流场的认识并支撑内转式进气道设计潜力的挖掘。

1 模型和数值方法

如图1所示，利用轴对称内收缩直锥生成内锥形激波，采用原点位于模型入口圆心处的轴对称坐标系，其中x为轴线方向，r为半径方向。该模型的壁面前缘半径R= 100 mm，前缘压缩角度为θw，无量纲的壁面长度为w/R。参考内转式进气道基准流场设计参数[28]，选取的壁面前缘角度范围为θw= 8°～15°，通过改变θw得到具有不同初始强度的激波，继而调节内锥形激波的汇聚增强过程。在每个给定的θw工况下，通过选取不同的壁面长度w/R来改变壁面尾缘膨胀波的位置，进而产生不同的马赫盘下游流动条件。当壁面太短时，尾缘膨胀波靠近上游，会对入射的内锥形激波产生干扰；而壁面太长时，内收缩比过大，超声速流场难以建立。因此，w/R取值约束在一定范围内[11,26]。表1列出了本文所有工况中θw和w/R的取值。表1中，w/R最大值和最小值为w/R可取的极限值，在本文中不作为主要讨论的工况。

图1 模型及计算域Fig. 1 Model and computational domain

表1 数值计算工况Table 1 Calculation conditions

鉴于本文主要关注内收缩直锥轴线附近的波系结构，黏性和边界层效应对远离壁面的主要激波及其反射结构的影响可以忽略[5,24]，因而通过求解理想气体（比热比γ= 1.4）的欧拉方程，对轴对称内收缩直锥流场进行CFD （Computational Fluid Dynamics）定常模拟。空间差分采用总变差减小（TVD）的迎风格式[29]，数值通量采用HLLC（Harten-Lax-van Leer-contact）近似黎曼求解器[30-31]计算。在所有计算工况下，均采用结构化网格将计算域离散。计算域如图1所示，除壁面和轴线边界外，其他边界均采用无反射的远场条件[19,24]。远场来流条件与文献[6]中保持一致，马赫数Ma∞= 6，静压p∞= 0.9 kPa。在计算过程中监测各方程残差，待残差下降4～5个数量级且波系结构稳定后，认为流场计算收敛。

为考核网格无关性，分别采用单元尺度为Δx/R≈Δr/R≈ 0.002和0.001的两套网格进行计算，其中Δx和Δr分别为单个网格沿x和r方向的长度。将获得的流场进行对比，发现马赫盘位置沿x方向的变化量仅约为0.001R，两套网格的主要波系结构几乎重合。因此，本文所有工况的数值模拟都采用单元尺度为Δx/R≈ Δr/R≈ 0.001的网格。

为验证本文数值模拟方法的可靠性，以θw= 10°、w/R= 1的工况为例，与文献[6]中的实验结果进行了对比。如图2所示，将数值模拟得到的马赫数等值线（红色实线）叠加在风洞实验[6]纹影照片上可以看出，数值模拟给出的马赫盘位置稍靠下游，与实验结果相距仅约为0.0008R，两者的主要激波反射结构基本重合。这也表明本文的数值模拟结果是可靠的。

2 结果与讨论

2.1 轴对称内锥流动基本结构

鉴于w/R不影响轴线附近的激波反射类型，本节统一选取w/R= 1.5，并以θw= 10°工况为例，简要介绍轴对称内锥流动中的基本激波结构。

图3（a）展示了θw= 10°、w/R= 1.5工况的流场马赫数（Ma）云图及等值线，图3（b）将激波反射结构所在区域（虚线框Ⅰ）进行了局部放大显示。结合图3（a）和图3（b）可以看出，壁面前缘产生的轴对称内锥形入射激波IS在向下游发展的过程中持续地汇聚增强，尤其在接近轴线（r/R= 0）时出现较为剧烈的弯曲。最终，IS在轴线附近发生马赫反射，汇聚过程被马赫盘MD终结，同时在三波点T处形成反射激波RS和剪切层SL（在本文讨论中不考虑SL的厚度，将其视为滑移线[11]）。壁面尾缘产生的膨胀波EW会透射过RS并入射在SL上，对SL产生干扰。图3（c）将MD下游的局部流场（图3（a）中虚线框Ⅱ）进行了放大显示，可以看出，在紧邻MD下游的区域，SL下方与轴线包裹的流管内是亚声速气流；随着SL向下游延伸，流管内的气流逐渐加速至超声速。换言之，在SL下方与轴线包裹的流管内，会形成声速喉道。在上、下游流动的影响下，SL的形态在声速喉道的形成过程中，扮演着重要的角色。

在来流Ma∞= 6条件下，由平面激波反射理论可知[11]，对于气流偏转角同样为θw（见表1）的平面激波而言，理论上均不存在马赫反射。对于轴对称内锥形激波而言，虽然其初始强度等同于气流偏转角为θw的平面激波（壁面前缘处的波后气流偏转角完全由壁面约束），但是在汇聚增强的过程中，内锥形入射激波IS会超出von Neumann准则对应的激波强度[11]（称为von Neumann强度[32]），使得马赫反射对应的入射激波强度条件得以满足（见图3）。鉴于马赫盘位置应当在入射激波强度达到von Neumann强度之后，允许马赫盘存在的空间范围取决于轴对称内锥形激波的汇聚增强过程。此外，在该空间范围内，三波点需与马赫盘下游的流动条件相匹配，因此马赫盘的位置受其下游流动的影响。

为深入探究影响马赫盘位置的机制，需要从两个方面入手：首先，通过改变壁面前缘角度，考察激波汇聚增强过程对马赫盘位置的影响；然后，通过改变壁面长度，进一步考察马赫盘下游流动对马赫盘位置的影响。

2.2 壁面前缘角度对马赫盘位置的影响

由2.1节可知，壁面前缘角度θw决定了轴对称内锥形入射激波的初始强度，并通过影响激波汇聚过程来限制允许马赫盘存在的空间范围。因此，阐明θw对内锥形激波汇聚过程的影响是十分必要的。

根据谭廉华等[19-20]提出的理论，轴对称内锥形入射激波形状r=r(x)可以表达为二阶非线性常微分方程式（1）和式（2）：

其中，rx和rxx分别表示r关于x的一阶导数和二阶导数，入射激波后的气流偏转角θi、压力pi以及马赫数Mai均是入射激波角βi的函数。对于本文而言，式（1）存在形如式（3）和式（4）的两个边界条件，

其中，壁面前缘处的激波角βw由Ma∞和θw确定。可见，θw是影响激波汇聚增强过程的重要参数。

以w/R= 1.5工况为例，采用Runge-Kutta法[33]对式（1）～式（4）进行积分，得到了不同θw下的入射激波形状，如图4所示。同时，在图4中还叠加了本文数值模拟获得的入射激波及马赫盘形状。对比图4中的结果可知，在不考虑发生马赫反射的情况下，谭廉华等[19-20]的理论能够较为准确地给出入射激波的形状。由于式（1）～式（4）是基于入射激波波后流线为直线这一假设[19-20]推导而来的，而实际流场中入射激波波后流线存在曲率[21]，因此，图4中的理论和数值模拟结果存在一定的偏差。然而，这种偏差很小，式（1）～式（4）仍不失为一种描述入射激波形状的有效理论。需要特别注意的是，式（1）～式（4）在轴线处（r/R= 0）存在奇性，在入射激波到达轴线之前，必然会发生马赫反射并形成马赫盘[15-16]。由于本节主要关注远离轴线处的入射激波，不影响对内锥形激波汇聚过程的讨论。

图4 不同θw下轴对称内锥形入射激波形状(w/R = 1.5)Fig. 4 Shapes of axisymmetric internal conical shocks at different θw (w/R = 1.5)

鉴于式（1）～式（4）对内锥形激波汇聚增强机理的反映不够直观，因而结合CST对图4中的结果进一步分析。针对任意曲面激波，CST方程[21]给出了波后当地沿流线的无量纲压力梯度的表达式。对于均匀来流下的轴对称内锥形入射激波而言，其波后沿流线的无量纲压力梯度P2i可表示为式（5）[22]：

式（5）中a1和a2均为与Ma∞和βi有关的系数，Sai和Sbi分别为流平面（激波当地来流和波后气流偏转方向共同决定的平面）和法平面（与流平面内激波当地切向垂直的平面）内的激波曲率[22]。在壁面前缘处（x/R= 0），可以将式（5）化简为式（6）[8]：

式（6）中系数c与初始激波强度正相关，Sci为激波的初始横向曲率。由于在壁面前缘处R是确定的（Sci确定），根据式（6）可知，c越大（θw越大），初始时激波汇聚增强越快（P2i越大）。

在向下游的汇聚过程中，由于激波波后流动方向不完全受壁面约束，式（6）将不再精确适用。因此，根据式（1）～式（4）计算得到的激波形状r=r(x)（见图4），再结合式（5）可以理论给出P2i的沿程分布。如图5所示，在壁面前缘处（x/R= 0），P2i随着θw的增加而增大，这与式（6）的描述相一致。前文图4已经表明，在下游同一x/R位置处，当θw越大时，入射激波角增大，并且入射激波到轴线的距离更近。这也意味着激波面的曲率Sai和Sbi均增大。相应地，根据式（5）可知P2i也随着增大（见图5）。进一步表明了随着θw的增大，内锥形激波汇聚增强得更快。换言之，在激波强度和激波面曲率的共同作用下，初始较强的激波（θw大）汇聚增强始终快于初始较弱的激波（θw小）。

图5 不同θw下轴对称内锥形激波波后沿流线的无量纲压力梯度分布(w/R = 1.5)Fig. 5 Distributions of the normalized pressure gradient behind the axisymmetric internal conical shocks at different θw (w/R = 1.5)

为了认识内锥形激波汇聚过程中激波自身强度的变化，如图6所示，在数值模拟流场中沿着入射激波前缘到三波点提取了激波前后的压比pi/p∞，并利用谭廉华等[19-20]的理论，由式（1）～式（4）给出了pi/p∞沿x/R的分布。图6中理论和数值模拟获得的pi/p∞分布整体符合较好，在接近三波点时，由于激波波后流线曲率较大[21]，理论结果存在一定的偏差。随着下游激波强度的持续增大，在内锥形激波的汇聚过程中，必然会达到von Neumann强度，图6中用虚线给出了本文Ma∞= 6条件下von Neumann准则对应的激波压比pN/p∞≈ 9.67。在不同的θw下，内锥形激波汇聚增强的快慢不同（见图5），导致入射激波上对应于von Neumann强度（pi/p∞=pN/p∞）的位置存在差异。在给定的θw下，根据式（1）～式（4）描述的入射激波形状，可以得到入射激波上von Neumann强度的对应点（xN/R，rN/R）。

图6 不同θw下轴对称内锥形激波压比的沿程分布(w/R = 1.5)Fig. 6 Pressure ratio distributions of the axisymmetric internal conical shocks at different θw (w/R = 1.5)

如图7（a）所示，在θw= 7°～16°范围内，利用式（1）～式（4）分别计算出了xN/R和rN/R随θw的变化曲线，并且标记出了本文数值模拟流场中各θw工况对应的xN/R和rN/R位置。从图7（a）可以看出，随着θw的增大，xN/R减小，而rN/R迅速增大。进一步地，将这些由式（1）～式（4）计算出的（xN/R，rN/R），在（x/R，r/R）平面内依次展示，可以得到如图7（b）所示的曲线。鉴于当入射激波强度达到von Neumann强度时，马赫反射才有可能发生，因而图7（b）中的（xN/R，rN/R）曲线和轴线（r/R= 0）共同限定了不同θw下允许马赫盘存在的空间范围。换言之，式（1）～式（4）为从理论上确定马赫盘的位置提供了一种可能的途径。

图7 不同θw下von Neumann强度对应点和三波点(w/R = 1.5)Fig. 7 Point of von Neumann strength and triple point at different θw (w/R = 1.5)

为对比三波点与（xN/R，rN/R）之间的位置关系，从数值模拟流场中提取了不同θw下的三波点位置和von Neumann强度对应点的坐标，叠加在图7（b）中。通过数值模拟结果可以看出，在不同的θw工况下，三波点均被约束在von Neumann强度对应点（xN/R，rN/R）的下游，两者相距很近，但并不重合。随着θw增大，轴对称内锥形激波更快地达到von Neumann强度，但马赫盘的位置并非完全依赖于激波汇聚增强过程。由此可见，仅通过式（1）～式（4）结合CST方程难以完全确定马赫盘的位置，还需要考虑下游流动条件。

2.3 马赫盘下游流动对其位置的影响

针对马赫盘下游的流动问题，一般可以将滑移线SL包裹的流管（见图3）视为准一维流动[11,34-36]。根据马赫盘面积Am和流管内声速喉道面积As，可以建立形如式（7）的等熵关系[11,26]：

式（7）中为紧邻马赫盘下游的气流平均马赫数。马赫盘位置与其下游流动之间的匹配，往往可以通过声速喉道来实现。

对于平面激波反射问题，已经发展出相应的几何模型[11,26]来描述马赫杆下游的流动。如图8所示，起初从三波点T′发出的滑移线SL′具有恒定的角度（即三波点后气流偏转角）。紧接着，壁面尾缘产生的膨胀波穿过反射激波RS′后，作用在SL′上，使得SL′的角度连续减小（D′A′段），进而形成了声速喉道。将这一几何模型及相应的二维激波/膨胀波关系与式（7）相结合，即可确定马赫杆MS′的位置[11,26]。这种方法为轴对称内锥形激波反射问题中马赫盘位置的确定，提供了极大的启示。然而，在轴对称内锥形激波反射流场中，马赫盘下游的流动在受到壁面尾缘膨胀波干扰之前，先需要匹配反射激波后的非均匀压力，导致滑移线经历的变化比二维情况要复杂得多。因此，需要对轴对称内锥形激波反射流场中马赫盘下游流动进行深入考察。

图8 平面激波反射示意图Fig. 8 Schematic of planar shock reflection

2.3.1 马赫盘下游流动特征

首先关注壁面尾缘膨胀波干扰之前，马赫盘下游的流动特征。对于不同的θw，当w/R充分大时（见表1），下游壁面尾缘产生的膨胀波离马赫盘较远，使得滑移线SL包裹的流管充分发展。根据流管自身能否演化形成声速喉道，将马赫盘下游流动划分为两类。

第一类工况：在壁面尾缘膨胀波干扰之前，马赫盘下游滑移线SL包裹的流管能够形成声速喉道。θw= 8°或10°且w/R充分大时，均属于此类工况。以θw= 10°、w/R= 3.5工况为例进行详细分析。图9展示了该工况下的流场局部无量纲压力（p/p∞）云图及其等值线，可以看到，入射激波IS和反射激波RS的波后压力都是不均匀的。以穿过三波点T1的流线表征马赫盘MD下游的滑移线SL，并以MD下游轴线（r/R= 0）上的马赫数表征流管内的马赫数分布，量化分析马赫盘下游的流动特征。

图9 θw = 10°、w/R = 3.5工况数值模拟流场的局部无量纲压力云图Fig. 9 Local normalized pressure contours of the numerical flow field at θw = 10° and w/R = 3.5

图10给出了SL包裹的流管面积A= πr2与马赫盘面积Am之比的沿程变化情况，其中，xm为马赫盘MD在轴线处的x坐标。同时，图10中还叠加了流管内的马赫数（Ma）分布，并标明了流管内首个声速喉道的位置。结合图9和图10可知，SL包裹的流管从T1开始快速收缩，使得紧邻马赫盘下游的亚声速气流加速。在A1点处，流管内形成了声速喉道。值得注意的是，在该工况下，壁面尾缘膨胀波与SL干扰的起始点D1不仅远离声速喉道A1，而且两者之间还有超声速区相阻隔。可见，声速喉道A1的形成原因，是独立于壁面尾缘膨胀波的。这与图8中的平面激波反射流场存在本质差异。

图10 θw = 10°、w/R = 3.5工况流管面积变化及其内马赫数分布Fig. 10 Area ratio of the stream tube and distribution of Mach number in the stream tube at θw = 10° and w/R = 3.5

进一步地，图11给出了反射激波RS（从T1点到F1点）波前、波后的压力以及滑移线SL上（从 T1点到D1点）的压力分布，以便于通过分析马赫盘下游流动如何与反射激波后的非均匀压力场相匹配，进而揭示马赫盘下游滑移线SL包裹的流管内声速喉道的形成原因。对比图11和图6可以看出，RS的波前压力与入射激波IS波后压力的变化趋势并不一致。这是因为在轴对称内收缩流场中，IS波后不同位置的气流在到达RS之前，经历了不同程度的等熵压缩过程（见图9）。这些气流经过RS再次压缩后，压力呈先降低（见图11的T1E1段）再升高（见图11的E1F1段）的变化趋势，而RS波后非均匀分布的压力将直接影响马赫盘下游的滑移线SL及流管形态。

图11 θw = 10°、w/R = 3.5工况反射激波前后及滑移线上无量纲压力分布Fig. 11 Pre-shock and post-shock normalized pressure distribution of the reflected shock, and normalized pressure distribution along the slip line at θw = 10° and w/R = 3.5

从图11可以看出，滑移线SL上T1B1段的压力连续下降，以匹配RS在T1E1段的波后压力；而滑移线SL上B1D1段的压力连续上升，以匹配RS在E1F1段的波后压力。结合图10和图11可知，在滑移线SL上T1B1段的降压过程中，SL包裹的流管内原本亚声速的气流连续加速至超声速。换言之，流管内声速喉道的形成是匹配RS波后非均匀压力的结果。RS和SL之间（见图9中T1D1F1区域）通过等熵膨胀或压缩来实现压力变化，尽管这种变化关系不如二维平面激波反射流场那样清晰，但从理论上建立分析方法仍然是将来值得尝试的研究课题。

第二类工况：在壁面尾缘膨胀波干扰之前，马赫盘下游滑移线SL包裹的流管无法形成声速喉道。当θw增加至12°或更大（θw= 15°）且w/R充分大时，均属于此类工况。图12以θw= 12°、w/R= 3.0工况为例，展示了流场的无量纲压力云图。图13给出了滑移线SL包裹的流管面积的沿程变化情况(A/Am)以及流管内的马赫数分布。图14展示了反射激波RS（从T2点到F2点）波前、波后的压力以及滑移线SL上（从 T2点到A2点）的压力分布。

图12 θw = 12°、w/R = 3.0工况数值模拟流场的局部无量纲压力云图Fig. 12 Local normalized pressure contours of the numerical flow field at θw = 12° and w/R = 3.0

图13 θw = 12°、w/R = 3.0工况流管面积变化及其内马赫数分布Fig. 13 Area ratio of the stream tube and distribution of Mach number in the stream tube at θw = 12° and w/R = 3.0

与前文θw= 10°工况类似，马赫盘下游流管内的压力也需要匹配RS的波后压力。有所不同的是，θw= 12°时RS在T2E2段的波后压力（见图14）不像θw= 10°工况那样迅速降低（见图11）。因而θw= 12°时马赫盘下游流管内气流的加速更为缓慢，以至于在压力谷值点B2处，气流仍未达到声速（见图13）。在B2点下游，由于要匹配RS在E2F2段波后压力的升高，流管内的亚声速气流连续减速增压（B2G2段），直至壁面尾缘膨胀波在D2点处对滑移线SL产生干扰。在马赫盘下游流管内的亚声速流动中，由于下游扰动能够传到上游，在D2点上游不远处的G2点已经能够感受到壁面尾缘膨胀波对SL的影响。由膨胀波带来的压力降低（见图14），使得流管内的气流急剧加速，直至在A2点处形成声速喉道（见图13）。可见，θw= 12°时马赫盘下游流管内的气流无法通过匹配RS的波后压力达到声速，而声速喉道的形成依赖于壁面尾缘膨胀波对SL的干扰。

图14 θw = 12°、w/R = 3.0工况反射激波前后及滑移线上无量纲压力分布Fig. 14 Pre-shock and post-shock normalized pressure distribution of the reflected shock, and normalized pressure distribution along the slip line at θw = 12° and w/R = 3.0

综上所述，在两类马赫盘下游流动中，反射激波RS波后压力下降的快慢程度不同，导致马赫盘下游流管内声速喉道的形成过程存在是否依赖于壁面尾缘膨胀波的差异性。由于马赫盘与声速喉道之间具有如式（7）所示的平衡关系，马赫盘位置对壁面尾缘膨胀波的依赖性也将不同。因此，需进一步考察壁面尾缘膨胀波对马赫盘位置的影响。

2.3.2 膨胀波对马赫盘位置的影响

在给定壁面前缘角度θw和入口半径R的条件下，w/R决定了流场的内收缩比以及壁面尾缘膨胀波的存在区域。在本文关注的w/R范围内，壁面尾缘膨胀波对马赫盘下游流动产生干扰的位置随着w/R的增加而向下游移动。

以θw= 10°为例，不同w/R工况下的入射激波及马赫盘形状如图15所示。从图15中马赫盘局部区域（虚线框）的放大图可以看出，不同w/R工况下，入射激波的形状均重合，但马赫盘位置存在差异。这种差异与2.3.1节讨论的声速喉道的形成机制紧密相关。

图15 不同w/R下轴对称内锥形入射激波及马赫盘形状(θw = 10°)Fig. 15 Shapes of axisymmetric internal conical shocks and Mach disks at different w/R (θw = 10°)

图16提取了θw= 10°不同w/R工况中马赫盘下游流管内的马赫数分布，其中T1标明了三波点的位置。以壁面长度充分大的w/R= 3.5工况为参照，A1点标明了其声速喉道的位置，而其他w/R工况下流管内的马赫数分布从T1下游不同位置开始与w/R= 3.5工况产生差异。如2.3.1节所述，这种差异是由壁面尾缘膨胀波对滑移线产生干扰而导致的。在w/R= 1.1～2.0范围内，壁面尾缘膨胀波对滑移线产生干扰的位置（图16中三角符号）在A1点上游，并辅助马赫盘下游流管迅速形成声速喉道（图16中矩形符号）。在w/R= 2.5～3.7范围内，壁面尾缘膨胀波对滑移线产生干扰的位置在A1点下游，此时马赫盘下游流管中声速喉道形成于A1点，不依赖于壁面尾缘膨胀波。

图16 θw = 10°工况不同w/R下的流管内马赫数分布Fig. 16 Distributions of Mach number in the stream tubes with different w/R at θw = 10°

进一步地，在θw= 10°工况下图17（a）分别以w/R= 3.5时的马赫盘位置xm*/R和声速喉道位置xs*/R为参考值，给出了不同w/R工况的马赫盘位置和声速喉道位置的变化情况。可以看出，在w/R=1.1～2.0范围内，随着w/R增大，声速喉道明显向下游移动；同时，马赫盘向上游移动。需要注意的是，马赫盘位置的变化幅度远小于声速喉道位置的变化幅度。马赫盘的位置仅在w/R＜ 1.5时变化明显，在w/R＞ 1.5之后变化微小。结合前文2.2节可知，由于在轴线附近入射激波的形状将越来越陡峭，马赫盘位置的明显移动，必然引起马赫盘面积的大幅变化。

图17 θw = 10°工况不同w/R下的马赫盘、声速喉道的位置和面积Fig. 17 Positions and areas of Mach disks and sonic throats with different w/R at θw = 10°

图17（b）以w/R= 3.5时的马赫盘面积Am*和声速喉道面积As*为参考值，给出了不同w/R工况下的马赫盘面积和声速喉道面积的变化情况。可以看出，由于两者遵循式（7）的关系，两者的变化曲线几乎重合。在w/R＜ 1.5时，马赫盘面积和声速喉道面积显著减小；在w/R= 1.5～2.5范围内，尽管声速喉道的位置随着w/R增大而向下游大幅移动（见图17（a）），其面积变化却很小（见图17（b））。这主要是由于在w/R= 1.5～2.0范围内，当壁面尾缘膨胀波入射滑移线时，马赫盘下游流管内的气流已经接近声速，流管的微小面积变化也会使得流管内快速形成声速喉道（见图16）。在w/R= 2.5～3.7范围内，声速喉道的形成不依赖壁面尾缘膨胀波，因此，从图17（a）和图17（b）中可以看出，声速喉道及马赫盘的面积和位置，均基本保持恒定。类似地，在θw= 8°、w/R= 1.5～4.3工况下，壁面尾缘膨胀波不影响声速喉道的形成。

对于θw= 12°和15°工况，在本文研究的w/R范围内（见表1），马赫盘下游流管内声速喉道的形成均依赖于壁面尾缘膨胀波。图18以θw= 12°工况为例，给出了不同w/R工况马赫盘下游流管内的马赫数分布。图19（a）和图19（b）分别展示了不同w/R工况马赫盘、声速喉道的位置和面积，其中，参考值xm*、Am*和As*均选自w/R= 3.0工况。结合图18和图19（a）可以看出，随着w/R增加，壁面尾缘膨胀波入射位置向下游移动，声速喉道的位置也随之向下游移动，几乎呈线性变化。然而，马赫盘位置的变化幅度远小于声速喉道位置的变化幅度，马赫盘向上游移动的趋势逐渐减缓，在w/R＞ 1.5之后变化微小。从图19（b）可以看出，马赫盘和声速喉道的面积同样在w/R＞ 1.5之后变化很小。其原因与θw= 10°、w/R= 1.5～2.0工况下马赫盘、声速喉道面积趋于稳定的原因类似。

图18 θw = 12°工况不同w/R下的流管内马赫数分布Fig. 18 Distributions of Mach number in the stream tubes with different w/R at θw = 12°

图19 θw = 12°工况不同w/R下的马赫盘、声速喉道的位置和面积Fig. 19 Positions and areas of Mach disks and sonic throats with different w/R at θw = 12°

在本文研究的轴对称内锥流场中，马赫盘下游声速喉道的形成机制存在是否依赖于壁面尾缘膨胀波的两类情况。对于不依赖壁面尾缘膨胀波的情况（如θw= 8°、w/R= 1.5～4.3工况以及θw= 10°、w/R=2.5～3.7工况），马赫盘下游流管通过与反射激波后的非均匀压力相匹配，进而形成声速喉道；对于依赖壁面尾缘膨胀波的情况（如θw= 8°、w/R= 1.2工况，θw= 10°、w/R= 1.1～2.0工 况 以及θw= 12°和15°工况），马赫盘下游流管还需要进一步匹配壁面尾缘膨胀波引起的压力降低，才能够形成声速喉道。这种匹配关系，为理论求解马赫盘位置提供了启示。

2.4 马赫盘位置的理论求解

结合2.2节和2.3节的分析可知，马赫盘位于入射激波上von Neumann强度对应位置所限定的范围内，而马赫盘下游的流动还需要同时满足反射激波下游的压力变化以及准一维等熵关系式（7）。换言之，确定马赫盘的位置需要求解反射激波下游的超声速流场。最近，Shi等[25]发展的弯曲激波特征线方法（Method of curved-shock characteristics，MOCC）具有高效且精确求解超声速流场的优势。本节在给定来流和壁面前缘入口条件下，借助MOCC求解马赫盘位置。

首先，假定马赫盘处于入射激波上的von Neumann强度对应点（见2.2节）。利用MOCC，可以准确求解入射激波从壁面前缘到三波点的形状，而且可以进一步获得入射激波下游的超声速流场，包括入射激波与反射激波之间的流场、反射激波形状及其波后流场。滑移线作为反射激波波后边界上的流线，其形状在利用MOCC的求解过程中可以确定[37]。根据滑移线形状，计算出滑移线包裹的流管面积，并找到流管面积的第一个最小值Asc。

上述滑移线形状仅是基于反射激波下游的流场得到的，而反射激波下游压力变化还需要与马赫盘下游流管内的压力相匹配。因此，利用准一维流动的等熵关系式（7），验证马赫盘下游流动是否符合2.3.2节中的描述。由于事先假定了马赫盘位置，式（7）中的马赫盘面积Am及流管内初始马赫数均为已知量[19]，利用式（7）可以得到声速喉道面积的理论值Ast。将Asc与Ast进行对比，若两者之差在设定的误差范围内，则认为Asc是准确的，表明马赫盘下游流动既能匹配反射激波波后的非均匀压力变化，又满足准一维流动的等熵关系，即假定的马赫盘位置正确；若在设定的误差范围之外，则需要将假定的马赫盘位置沿x方向调整后，重新计算一遍滑移线包裹的流管面积变化，并再次对比新获得的Asc与Ast的值。如此迭代求解，直至最终得到马赫盘的位置。

图20以θw= 12°工况为例，利用上述理论方法求解了马赫盘位置随w/R的变化情况。作为对比，图20中还给出了CFD流场中的马赫盘位置。可以看出，理论与CFD结果吻合良好，两者之间的最大差异仅为0.002。这表明，本节的理论方法，能够准确地预测出马赫盘随w/R的变化。鉴于MOCC在理论求解效率和精度方面的独特优势，可望在更宽的来流条件和内锥几何参数范围内，快速预测马赫盘的位置。

图20 θw = 12°工况不同w/R下马赫盘位置的理论和数值对比Fig. 20 Theoretical and numerical positions of Mach disks with different w/R at θw = 12°

3 结 论

在来流Ma∞= 6条件下，采用无黏数值模拟研究了轴对称直内锥流动中影响马赫盘位置的因素，主要得到了以下结论：

1）轴对称内锥形入射激波自身的汇聚增强过程影响马赫盘的位置。结合入射激波形状理论与弯曲激波理论，阐明了随着壁面前缘角度θw增大，入射激波的汇聚增强越快，并且能够更快地达到von Neumann强度，该强度的对应位置限定了允许马赫盘存在的空间范围。

2）马赫盘下游声速喉道的形成机制存在是否依赖于壁面尾缘膨胀波的差异性，进而对马赫盘位置产生不同的影响。当声速喉道的形成不依赖于壁面尾缘膨胀波时，声速喉道及马赫盘的位置和面积均不随w/R增大而变化。当声速喉道的形成依赖于壁面尾缘膨胀波时，随着w/R的增加，声速喉道的位置向下游大幅移动，但声速喉道及马赫盘的面积逐渐趋于稳定。

3）将弯曲激波特征线方法与马赫盘下游流管的准一维流动模型相结合，可以快速求解马赫盘的位置。