飞机超声速机动飞行条件声爆预测方法

冷 岩，张劲柏，钱战森，*

(1.中国航空工业空气动力研究院，沈阳 110034；2.高超声速气动力/热技术重点实验室，沈阳 110034；3.北京航空航天大学 航空科学与工程学院，北京 100191)

0 引言

自“协和号”首飞以来，声爆问题很大程度上决定了超声速客机能否进入商业运营及取得商业成功，因而成为了研制超声速民用飞机必须解决的关键难题[1-3]。特别地，针对起飞阶段不可避免的超声速加速过程和快速机动飞行时产生的聚焦声爆现象（如图1 所示[4]）的研究成为了研制新一代环保型超声速飞行器的重要课题。

图1 不同飞行状态下四种聚焦声爆示意图[4]Fig.1 Illustration of four types of caustics under supersonic flight operations[4]

只有在超声速巡航状态下地面声爆特征才是标准N 波。超声爆现象一般产生于飞机起飞后的超声速加速阶段或变马赫数状态下发生转弯、拉升等机动动作时，一般呈非对称U 型，其过压峰值通常是超声速巡航状态下产生的N 波的2～5 倍[5-6]。因此，量化焦散面上的聚焦声爆对超声速飞行器研制具有重要意义。经典理论在几何声学框架下开展声爆数值预测研究[7-8]，其中相位函数由射线路径决定，信号幅值由射线管面积决定，沿每条射线的非线性效应解释了压力信号从近场复杂激波/膨胀波系到地面简单N 波的变化。以焦散面为中心的超声爆区域射线管面积为零，导致几何声学方法在焦散面附近因奇异性而失效。

国外超声爆现象的研究历史悠久，试验测量主要是在20 世纪六七十年代。Haglund 和Kane 针对加速引起的聚焦声爆开展了详细研究[9-10]。其他研究团队针对超声速飞行器加速和转弯引起的聚焦声爆开展飞行测量研究，并成功捕捉到超声爆现象[11]。Sanai[12]等证明利用弹道靶实现超声爆测量的可能性。在开展试验研究的同时，也从数值建模和分析的角度对超声爆进行了研究。数值建模研究始于对具有解析解[12]的线性Tricomi 方程和非线性Tricomi方程[13-14]（nonlinear Tricomi equation,NTE）的推导，之后众多学者[15-19]用NTE 方程对超声爆进行建模，进而预测得到焦散线附近的声爆信号。在声爆预测及远场传播过程中，必须考虑包含损耗机制在内的真实大气效应影响，热黏吸收效应和分子弛豫效应是两个主要耗散机制。通过在NTE 方程中添加吸收/色散效应，Salamone[20]等提出广义非线性Tricomi 方程（lossy nonlinear Tricomi equation,LNTE）并数值求解，与飞行试验数据对比表明，在聚焦声爆的研究中添加耗散机制是非常必要的。非线性Tricomi 方程为混合型双曲/椭圆方程，Sescu[21]等通过添加虚拟时间变量以保证守恒形式的方程为双曲型，并采用非线性限制器来控制求解域的振荡。在此基础上，Kanamori[22]等开发了超声爆评估工具FFnoise。

国内在近场声爆CFD 计算与基于非线性Burgers方程的远场声爆预测方面取得了较大进展。西北工业大学[23-26]、中国空气动力研究与发展中心[27]、中国航空工业空气动力研究院[28-30]、中国航空研究院[31]、北京航空航天大学[32]等多家单位都开展了相关研究，发展了一系列声爆的近、远场预测方法。但是，针对超声速飞行器机动状态下的超声爆现象的研究，尚未见公开发表的文章。

鉴于此，受Salamone 等[20,33]工作的启发，本文基于广义Tricomi 方程建立表征焦散线附近压力波演化过程的空间模型，采用伪时间推进和分裂算法，在时域和频域分开求解以达到数值模拟超声爆的目的，并开发相应计算机代码ARI_Superboom，线性解析解分析和飞行试验数据验证表明了模型及算法的准确性。

1 超声爆及其模型方程

1.1 超声爆

当飞行器超声速飞行时，由近场复杂激波/膨胀波系引起的压力脉动沿声射线传播。声射线垂直于马赫锥，因此在固定参考系下，巡航状态生成的声射线间彼此平行，互无交叉。当飞行器加速飞行时，随着马赫数的增加，马赫角逐渐减小，声射线间的交叉导致焦散的形成，如图2 所示。进一步，在衍射和非线性效应的共同作用下会在焦散附近产生最大聚焦区域。当声射线进入聚焦区域，沿射线传播的波形就会发生畸变，进而产生聚焦声爆或非对称U 波。这种变形发生在聚焦区域的任何位置，包括焦散面和地面之间的交叉点周围，因此，如图2 所示，可在地面观测到由加速度引起的超声爆现象。

图2 超声速飞行器加速状态聚焦声爆形成示意图[20]Fig.2 Illustration of the focus boom formation during the acceleration phase of a supersonic flight vehicle[20]

以匀加速为例（图2），Ma从0.96 加速到1.35，马赫数随时间变化率为0.004/s，飞行高度为13.716 km。加速过程中，当马赫数大于截止马赫数后出现聚焦区域。标准大气无风状态下截止马赫数为1.153 2，根据真实飞行状态，初始计算状态设为Ma=1.156，每隔1 s 计算一次声射线，得到如图3 所示声射线射线和焦散线。

图3 ARI_Boom 预测所得声射线和焦散线Fig.3 Rays and caustic predicted by ARI_Boom

1.2 模型方程

在三维空间中，焦散是一个面；在二维空间中，焦散是一条线，如图3 和图4 中红色线所示，其中x轴与焦散线相切于O点，z轴与焦散线垂直于O点。焦散线上方为强干扰区（明亮区），反映了入射压力波在焦散线附近的强非线性相互作用，在此区域内的每个位置都有入射和出射两束射线通过。焦散线下方为消散波区（阴影区），反映由于衍射效应投射的波，一般强度较弱。

图4 射线在焦散附近变化示意图[33]Fig.4 Illustration of rays in the vicinity of the caustic[33]

广义Tricomi 方程表征的是焦散附近压力波的演化图像，控制焦散线附近声爆信号变化的无量纲形式的广义Tricomi 方程如下[20]：

方程（1）左端前两项表示衍射效应，属于线性项；第三项和第四项表示自然来流影响，也属于线性项，Max=u0x/c0、Maz=u0z/c0分别表示自然来流马赫数在x向和z向的分量，u0为飞行器所在高度的来流速度；第五项为非线性项；第六项为耗散项，包含了热黏性和分子弛豫效应。无量纲声学压力=(p−p0)/pac，声学马赫数pac为输入波形特征参考压力（最大过压峰值）；ρ0、p0、c0分别为大气环境密度、压力和声速，ε为与焦散线相切于O点的声射线的曲率半径。

为了避免求解多维空间问题，Tricomi 方程对焦散线附近的现象进行了简化，建立了一维空间模型。通过合理地选择计算区域，可以借助定义与2 个维度空间位置相关的无量纲时间参数来表征压力波从入射到出射的演化过程；这一无量纲参数既可以看成是沿着焦散线的空间发展，也可以看作是压力自身的时间发展。式（1）中无量纲时间变量的定义为：

其中Rcau为焦散线曲率半径。

方程（1）中 ε为小参数，反映特征波长λac=c0/fac与衍射层厚度 δ的比率：

定义 µ为与衍射相关的非线性效应的衡量值，对于超声爆情形其量值在0.1 左右。根据式（5）非线性项系数可改写为：

其中，下标 υ可取为1 和2，分别代表氧气和氮气分子弛豫效应；τυ为分子弛豫时间；c∞,υ为冻结声速，c∞,υ=c0+∆c∞,υ，∆c∞,υ为声速增量，其近似表达式及参数可见参考文献[23，36-37]。

2 广义Tricomi 方程的数值解法

2.1 数值求解

广义Tricomi 方程是关于z和的双曲/椭圆混合型方程，通常引入伪时间变量后采用时间推进求解，当伪时间收敛时，就可以得到压力波信号在焦散线附近的分布。数学上的求解区域及图像如图5 所示，其中，即计算区域在方向为2 倍衍射边界层厚度。

图5 Tricomi 方程求解区域示意图[22]Fig.5 Illustration of the computational domain of the Tricomi equation[22]

方程（1）的数值求解过程中，引入虚拟时间变量σ，则方程变为：

基于傅里叶变换，数值求解通过两步分裂算法分别在时域和频域中实现。针对每一个无量纲角频率ωn，在频域中求解衍射和吸收/色散效应项，非线性项在时域内求解。在分裂算法的每一步中，设置伪时间增量的大小一致并限制其量值以避免非线性求解中出现激波过冲现象[29]。

第一步：在频域内求解衍射项、耗散项、x向大气风项和z向大气风项；

在该步的离散求解时，对于伪时间项采用一阶前差，z方向的导数项取中心差分，具体计算格式为：

第二步：在时域内采用激波捕捉法求解非线性项，如下式：

该步的时间离散承接上步，右端项采用经典的二阶Roe 格式。简略形式如下：

2.2 边界条件

上述数值求解过程中每一步都需要设置边界条件。时域内，假设输入波形到达前和输出波形离开后声学压力为0，因此，

式（14）中，F表示经过pac无量纲化的输入波形。值得注意的是，上述辐射边界条件可允许在不知道输出波形的情况下使用。频域内，对应的辐射边界条件经傅里叶变换后形式为[34]：

2.3 焦散线附近压力波演化过程模拟方法

图6 给出超声速飞行器匀加速状态聚焦声爆预测过程及Tricomi 方程求解区域。整个过程包括以下三步：第一步，基于航空工业气动院自主研发的ARI_Overset 高精度数值模拟软件，采用CFD 手段数值求解三维Navier-Stokes（N-S）方程得到近场流场结构并提取模型下方设定位置处空间压力分布，压力分布提取位置一般为1～3 倍飞行器特征长度处；第二步，基于自主研发的ARI_Boom 声爆预测程序开展射线追踪和声爆远场传播计算，获得焦散线及所需声射线信息和距离地面Zδ处声爆信号，此时所得声爆信号类似于标准N 波；第三步，将第二步所得声爆信号作为聚焦计算初始输入波形，通过自主研发的基于广义Tricomi 方程的ARI_Superboom 超声爆预测程序获得地面或指定位置聚焦声爆信息，由于声衍射等效应的联合作用，声爆信号从幅值、形状等方面发生了明显的变化。第一、第二步涉及的近场CFD 预测和远场声爆在文献[28-29]中均有验证，本文的重点主要体现在第三步焦散线附近压力波演化过程模拟方法中，为了验证其准确性，聚焦区域输入波形采用文献数据。

图6 聚焦声爆预测过程示意图[33]Fig.6 Illustration of the focus boom prediction[33]

3 方法验证

3.1 线性Tricomi 方程分析

一种快速的机动飞行声爆远场传播预测方法是采用线性Tricomi 方程的解析解计算聚焦区域的声爆。线性Tricomi 方程相对于广义Tricomi 方程仅保留了声衍射项，方程变为：

其中，IFFT 表示傅里叶逆变换，sgn为符号函数，Ai为艾里函数。

图7 初始输入N 波波形Fig.7 N-shape wave for the initial input

图8 计算所得空间压力云图ⅠFig.8 Pressure contour Ⅰobtained by the numerical simulation

图9 预测结果与解析解对比Fig.9 Comparison between the numerical and analytical solutions

3.2 SCAMP 项目飞行试验

NASA Dryden 中心基于F-18 飞机开展的SCAMP（superboom caustic analysis and measurement project）飞行试验提供大量聚焦声爆数据以验证数值预测方法和代码。试验过程中由81 个麦克风组成的地面线性阵列测量超声速飞机产生的聚焦声爆现象。通过飞行轨迹记录、高空/地面大气条件监测、声学数据推算等手段获得射线追踪和CFD 模拟所需的飞行条件。

本文选取“状态C”[33]作为验证算例进行讨论。F-18 飞机的飞行马赫数为1.23，马赫数变化率为0.003 5/s，俯仰角速度为−0.25°/s，飞行高度为12.8 km，迎角为2.3°，飞机特征长度为46.3 m，近场特征提取位置为3 倍特征长度处，近场特征如图10 所示[28]。

图10 聚焦声爆预测初始输入波形Fig.10 Initial input wave for the focus boom prediction

图11 计算所得空间压力云图ⅡFig.11 Pressure contour Ⅱ obtained by the numerical simulation

图12 预测结果与飞行试验数据对比Fig.12 Comparison between the numerical solution and the flight test data

图13 非线性项对预测结果影响对比Fig.13 Comparison of the numerical results calculated with and without the nonlinear term

图14 耗散项对预测结果影响对比Fig.14 Comparison of the numerical results calculated with and without the viscous term

4 结论

本文基于广义Tricomi 方程建立真实大气效应的聚焦区域压力波演化过程模型开发自研程序ARI_Superboom，采用频域和时域相结合的方式离散模型并数值求解，提高了计算精度和效率；鉴于飞机超声速机动飞行条件声爆预测方法的复杂性，基于航空工业气动院自研软件建立超声速飞行器匀加速状态下聚焦声爆预测分析手段，进而获得聚焦区域超声爆信息，并通过线性Tricomi 方程分析和SCAMP 项目飞行试验两个算例验证了ARI_Superboom 程序的准确性和过程分析的可行性；基于SCAMP 项目飞行试验实测数据开展效应分析研究，对比数据表明，在超声爆数值预测中必须考虑非线性效应和耗散机制的影响，但是影响区域集中在明亮区，对于衍射效应主导的阴影区影响较小。