1 引言

现阶段随着机体/发动机一体化设计概念的提出，对超声速流道高效转弯变形的需求也日益增长，圆截面等流道设计再次成为国际上的研究热点^{[1, 2]}，这类流道通常需要实现三维空间复杂曲面的均匀过渡，因此也是空气动力学研究的重要方向，能够使流动更加均匀，避免压缩波汇聚成激波，壁面逆压梯度引起的分离。此外圆形截面流道具有结构、防热^{[3, 4]}等方面的优势，并广泛应用于风洞喷管、扩压器、发动机燃烧室等部件中。

在喷管、扩压器、进气道设计中，基于特征线法的设计理论较为成熟^[5~7]，其核心概念就是所谓的“消波”。消波过程通过处理马赫波的入射与反射比例，将部分或全部的马赫波消除^{[8, 9]}，达到改变马赫数、压力及速度的目的，实现流场的加速/膨胀、减速/压缩、转弯、变形等功能。一旦出现消波不完全，流场内的马赫数就会出现非预期的分布，致使边界层出现非预期的堆积、分离^[10]，例如REST进气道中的流向涡始终是反压前传和压力损失的重要因素^[11]。因此，分析马赫波在圆形流管的传播过程，具有重要的工程与理论价值。

采用给定边界参数分布的方式以达到控制流场波系转化的目的，这是典型的边界反问题。2013年南京航空航天大学的全志斌等设计了壁面压力分布可控的单边膨胀喷管^[12]。2015年国防科技大学的刘红阳根据轴线马赫数分布提出了特征线追踪方法^[13]。2016年厦门大学的李怡庆等设计了壁面压力分布可控的高超声速进气道^[14]。

2016年，蓝庆生等在三维特征线^[15~19]的基础上构造了基于双特征线的压力反问题求解器^{[20, 21]}，本文在此求解器的基础上进行圆截面超声速流动的消波过程研究。本文拟研究圆形截面流管因局部变形而导致的消波不完全问题^{[22, 23]}，该问题涉及到马赫波的三维相交、反射、聚焦、耗散^{[24, 25]}。若需要消除非预期的马赫波，必须在马赫波传播的路径上引入几何的变形，使其按照预期压力分布进行反射，这种几何变形过程也是十分复杂的。

2 三维消波算法 2.1 特征线消波技术

压力波的消波需要用几何变形来实现，这是一个几何与流动强烈耦合的复杂问题。超声速内流中，给定入口来流条件，通过调整流道壁面的压力分布，可以获得实现该压力分布的流道几何型面，这是典型的压力反问题。

为实现压力反问题的求解，作者提出三维超声速压力反问题的特征线求解技术（iMOC-3D），如图 1（a）所示，已知平面1（黑色）上各点的坐标和流动参数，以及平面2（红色）预期的壁面压力分布，通过如图 1（b）所示的壁面单元过程，可以求得平面2上所有点的坐标和流动参数，从而完成空间上沿x方向的一个步进。通过平面1和平面2的不断推进，逐渐生成与给定壁面压力相容的三维型面。

图 1（b）中，在5点建立基准面，该基准面由向量OP1和向量OP2共同构成，其中OP1与5点左右两点的连线相互平行，OP2为向量（1，0，0）。直线56为流线微元。密切面由矢量V和上游的β共同构成，V为5点和6点的平均速度矢量，β为流动相切条件中的壁面外法线。流线56为密切面和基准面旋转θ角的平面相贯获得，旋转角度的正方向满足旋转轴为OP1的右手定则。从6点向来流方向追踪半个马赫锥，三条双特征线从6点向后延伸，与初值面相交于1，2，3点，1~3点及5点参数由拟合了5点和8个相邻点的双变量内插多项式确定，6点坐标和参数可以通过求解沿双特征线及流线上的特征方程和相容方程获得。特征方程及相容方程如下

$ {\rm{d}}x/{\rm{d}}t = u, {\rm{d}}y/{\rm{d}}t = v, {\rm{d}}z/{\rm{d}}t = w $

(1)

$ \begin{array}{l} \left[{{u^2}-\left( {{\boldsymbol{V}^2}-{a^2}} \right)} \right]{\left( {{\rm{d}}x} \right)^2} + \left[{{v^2}-\left( {{\mathit{\boldsymbol{V}}^2}-{a^2}} \right)} \right]{\left( {{\rm{d}}y} \right)^2} + \left[{{w^2}-\left( {{\mathit{\boldsymbol{V}}^2}-{a^2}} \right)} \right]\cdot\\ {\left( {{\rm{d}}z} \right)^2} + 2uv\left( {{\rm{d}}x} \right)\left( {{\rm{d}}y} \right) + 2uw\left( {{\rm{d}}x} \right)\left( {{\rm{d}}z} \right) + 2vw\left( {{\rm{d}}y} \right)\left( {{\rm{d}}z} \right) = 0 \end{array} $

(2)

$ \rho u{u_{\rm{t}}} + \rho v{v_{\rm{t}}} + \rho w{w_{\rm{t}}} + {p_{\rm{t}}} = 0 $

(3)

$ {p_{\rm{t}}} - {a^2}{\rho _{\rm{t}}} = 0 $

(4)

$ \begin{array}{l} \rho a{n_x}{u_{\rm{t}}} + \rho a{n_y}{v_{\rm{t}}} + \rho a{n_{_z}}{w_{\rm{t}}} - {p_{\rm{t}}} + \rho {a^2} \times \\ \left[\begin{array}{l} \left( {n_x^2-1} \right){u_x} + \left( {n_y^2-1} \right){v_y} + \left( {n_z^2-1} \right){w_z} + \\ \left( {{u_u} + {v_x}} \right){n_x}{n_y} + \left( {{u_z} + {w_x}} \right){n_x}{n_z} + \left( {{v_z} + {w_y}} \right){n_y}{n_z} \end{array} \right] = 0 \end{array} $

(5)

因此存在一个单调关系

$ {P_{6{\rm{i}}}} = f\left( {{\theta _i}} \right) $

(6)

6点压力会随着参数θ的增大而减小，反之亦然。通过给定6点压力，反求对应的θ值，从而获得6点的三维坐标，实现三维超声速压力反问题的求解。

2.2 iMOC-3D精度验证

文献[20]中采用PM膨胀波的理论解对iMOC-3D求解器的精度进行验证，为了文章的完整性，本小节进行简单陈述。设计的算例模型如图 2所示，外折角为20°，来流马赫数为1.4，总压为101.325kPa，总温为288K，气体常数为287J/(kg·K)，比热比为1.4。下壁面角点处产生的膨胀波扇在上壁面反射。假设图中矩形区域尚未受到反射膨胀波的影响，故其精确流场可用普朗特-迈耶函数求得，如图 3中黑色马赫数等值线所示。从图 2矩形区域顶点（0.2449，0.25）处追踪一条流线，将其上的压力分布输入iMOC-3D求解器，从而获得与该压力分布相容的壁面坐标和流场。

初始截面所在位置如图 2所示，设Z方向垂直于纸面向外，则初始截面在YZ平面内，分别沿Y，Z方向设定40，20个计算节点。X为流动方向，沿流向步进80步，获得如图 3所示的红色马赫数等值线。从图中重合部分可以看出，iMOC-3D计算的马赫线和理论解几乎完全重合。流线上各参数与精确解的最大相对误差见表 1，可以看出各参数的最大相对误差在万分之几量级，意味着iMOC-3D继承了特征线高精度的特性。

3 圆形截面超声速消波流道设计

如图 4所示，在圆截面流道中，若存在边界层堆积、分离或者由加工误差导致流道型面发生变化（如向内凸起），形成一个稳定的局部高压区，这个高压扰动以三维压力波的形式向下游扩散。因此，设计具有局部高压区的圆截面流道，以此来探索压力波在超声速圆管中的传播、相交、聚焦及耗散作用。即，给定如图 4所示的压力分布，在x=0.15m，φ=π处设置局部高压区，其中x为直角坐标系下的X轴方向的流动距离，φ表示在YZ平面内与Y轴正方向的夹角，局部高压区的区域半径为r=0.1m，圆心处叠加的压力峰值为p_K=1%p，内部压力可表示为r的的四次多项式，在高压/低压区边缘以二阶导数连续的方式向均压区过渡，峰值处一阶导数为0，均压区的压力一直保持为来流静压。

图 5为流管入口的O-block网格，每个Block网格量为25×25。流管半径R₀=0.5m，入口为马赫2.45的平行均匀流，静压p=160kPa，静温T=3149.67K，密度ρ=0.177kg/m³，比热比γ=1.25。

将如图 4所示的压力分布作为iMOC-3D的输入项，当自由来流流经给定的局部高压区时，相应的壁面向内凸起产生稳定的局部高压，高压扰动在流动中传播时，后续壁面通过局部变形来保证壁面压力为设定值。由此获得与该压力分布对应的三维构型，如图 6所示。

由于iMOC-3D求解器所设计的型面为无粘型面，而三维粘性修正的技术尚不成熟，并且开展纯三维内流实验难度较大，考虑到采用数值模拟是气动设计的常见验证手段，因此采用Fluent商用软件来验证本文的设计模型。图 7将iMOC-3D计算的流场压力与Fluent解进行对比，模拟的控制方程为Euler方程，网格量约为200万，图中对比了z=0m截面y=0.4m的压力分布，可以看出，两种求解器所得流场基本相同。因此采用iMOC-3D求解压力波的三维传播、反射、消扰等物理过程具有一定的可靠性。

4 结果与讨论

为了探索马赫波在传播过程中的反射、相交、聚焦和耗散，本文设计了如表 2所示的算例，通过算例间的对比，获得几何和流动耦合的一些规律。

4.1 压力波的聚焦与三维传播

算例A1~A6分别在x=0.15m，φ=π处存在局部扰动，图 8（a），（b）分别为A1~A3和A4~A6算例中z=0m截面的压力云图，图中可以直接看出扰动的峰值直接影响压力波的强度；压力波每经过一次反射，反射点下游位置均存在“人”型波系。图 9为A3算例中不同流向横截面的压力云图，对比图 8（a）中的A3算例，可以看出：来流流经局部高压区后，出现压缩波与膨胀波的马赫锥相互嵌套的三维传播形式，以图 9（b）最为明显；圆弧壁面对压力波有聚焦作用，使得压力波在壁面反射之后存在一个聚焦点，如图 9（d），9（f），表现了扰动传播的三维特性。

4.2 消波壁面的几何变形

如图 10所示，其中横坐标表示壁面点与圆心的距离与R₀的比值，φ表示在YZ平面内与Y轴正方形的夹角。图中展示了算例A1的壁面几何变型，图的右侧以x=3m，5m，7m处的横截面为例，壁面宏观上的形变非常微小。将x=3~8m的横截面沿周向展开，可以看出，沿扰动源的流向位置，其壁面几何变形主要呈现两个波峰夹一个波谷的形状，且该形变的量级要远远大于其余位置的形变。

图 11展示了算例A1中φ=0，φ=π/2，φ=π流向壁面R/R₀值随流向的变化。对比图 11（a）和（b）R/R₀值的数量级，可以看出沿φ=π流向的形变要远远大于其他位置（如φ=0，φ=π/2），并且沿φ=π流向将趋于以一个稳定的速度向内变形；图 11（b）中可以看出，当扰动入射至φ=0和φ=π/2壁面上时，对应壁面开始发生形变，并且变形总体上呈现向外扩张的趋势。综上所述，给定图 4所示的压力分布，对应曲面的形变为：在扰动源的下游位置（φ=π），呈现两峰一谷的波形，且该波形的形变量级要远远大于其余位置；当扰动源为局部高压时，波谷形变大于波峰形变，随着流向，波形逐渐被拉长，即波峰愈加凸出，波谷愈加凹陷，且凹陷的速度将趋于定值，据此推测凹陷的距离可能将圆管掰成两瓣，而波圆管的其余位置呈现外扩的趋势。因此，图 4所示压力分布对应的超声速流动，扰动的存在使得壁面不能逐渐回归原型而呈现发散的趋势，从而说明存在局部高/低压扰动的压力反问题的解是一个发散解，展现了流动的非线性。

4.3 压力波与壁面变形的耦合

将A1算例的圆管壁面沿φ=0切开，展开成图 12所示的∂R/∂x云图，下方为z=0m截面的压力云图，由于不同位置的壁面形变相差较大，同时考虑压力波是造成壁面突然变化的主要因素，因此采用∂R/∂x来描述扰动在壁面传播的轨迹。

图 12中可以清晰地看出马赫锥和壁面的相贯线，结合z=0m截面的压力云图、图 9压力波的三维传播方式，以及图 10，11壁面变形规律，可以发现：由于在φ=π流向存在局部高压扰动，在高压内部的迎风面时压力升高，产生压缩波，对应壁面向内凹陷，∂R/∂x小于0，伴随着压力逐渐升高，∂R/∂x的绝对值逐渐增大；越过峰值流至背风面时压力降低，产生膨胀波，壁面向内凹陷速度减慢，∂R/∂x的绝对值逐渐减小，此处涉及曲率与压力曲线、波系的关系，限于篇幅就不进行展开讨论。当壁面恰好过渡至均压区时，内部流场正处于膨胀波后，为保证壁面压力不变，该膨胀波的入射点通过改变下游曲率，消除膨胀波，产生压缩波，如此往复，使得在扰动下游的流场存在交替出现并逐渐减弱的膨胀波、压缩波，使∂R/∂x逐渐趋于稳定，直至再次遇见反射的压力扰动。

4.4 多扰动源产生压力波的相交过程

为了进一步探索压力波的传播、相交过程，设计了如表3所示的C1~C4算例。C1~C3算例在（0.15，π/2），（0.15，3π/2）处均存在压力扰动，它们分别对应为局部低/低、低/低、高/低压区，而算例C4的两个局部低压区分别设置在（0.15，π/2），（1.05，3π/2）。图 13为算例C1和C4的∂R/∂x云图，图中可以看出由于扰动的存在，使得壁面产生了极为复杂的变化，而壁面曲率的变化也使得流场更加复发。

图 14对比了算例C1~C4中y=0截面的压力云图，可以看出不同强度、不同位置的压力波笔直地穿过了另一侧扰动源发出的压力波，现将单一扰动如A1~A6相应变换后进行数值叠加，并与组合扰动C1~C4进行对比，如图 15所示，叠加之后流场的压力均方根与组合扰动C1~C4吻合得较好，因此这些弱压力波的相交是一个线性叠加过程，没有体现流动的非线性。同时压力均方差随流向逐渐减小，可知压力波经过一些反射、相交、聚集之后，逐渐扩散，扰动的峰值减小，范围增大，避免汇聚成强扰动。

5 结论

本文对圆截面超声速流动的消波过程进行了研究，可以得到以下结论：

（1）在局部高压区下游，出现压缩波与膨胀波马赫锥相互嵌套的三维传播形式，同时圆弧壁面对压力波有聚焦作用。

（2）在消波过程中圆截面流道的形变主要呈现为两峰一谷/两谷一峰的波形，波形沿流向不断扩大而不能回归原型。扰动的存在使得壁面变形呈现发散的趋势，说明对应的压力反问题的解是一个发散解，显示了流动的非线性。

（3）压力波的消除问题与壁面曲率变化直接相关。

（4）弱压力波的相交未发现强非线性，可进行线性叠加。