基于广义Burgers 方程的声爆传播特性大气湍流影响

王迪，冷岩，杨龙，韩忠华，钱战森，*

1.西北工业大学 航空学院 翼型、叶栅空气动力学国家级重点实验室，西安 710072

2.中国航空工业空气动力研究院 高超声速气动力/热技术重点实验室，沈阳 110034

3.中国航空工业空气动力研究院 高速高雷诺数气动力航空科技重点实验室，沈阳 110034

声爆作为发展超声速民机不可回避的关键瓶颈问题，是未来新一代绿色超声速民机能否取得商业成功的重要羁绊［1-5］。当飞机进入超声速飞行时，其周围流场产生的激波与膨胀波系经过大气环境中的各种宏、微观效应的影响，到达地面后通常会演化成一道N 型波。然而，在地表附近的大气边界层内往往存在有大气湍流等非静止扰动，严重干扰声爆特征的演化过程，从而对声爆预测带来了一定的不确定性。对于未来中国超声速民机飞行轨迹规划，精准的声爆预测方法是其中的关键技术［6-8］。因此，开展大气湍流效应对声爆特性的影响规律研究具有重要意义。

大气边界层通常是指离地球表面约0.3～3 km高度的底层大气［9］，由于受地球表面粗糙度和热对流效应影响较为强烈，该区域的大气运动具有明显的湍流性质［10］，在短期内呈现出高度非定常性，导致声爆信号的幅值发生随机性增强或衰弱现象。图 1 展示了NASA 曾开展的飞行试验测量结果［11］，可以发现大气湍流对声爆特征的影响主要呈现出两个特点［12-14］：一是激波后可能存在“钉状细致结构”，导致声爆的过压峰值较大，最终演化成P 型波（Peaked Waveform）；二是激波峰值可能被抹平，导致声爆上升时间的随机性增加，最终演化成R 型波（Rounded Waveform）。

图1 不同大气条件下测量得到的地面声爆波形［11］Fig.1 Sonic boom waveforms measured under different atmospheric conditions on ground［11］

目前的远场声爆预测技术主要是以线性化和拟线性化为基础［15-17］发展而来的几何声学方法，该方法以Fermat 理论来确定声射线的轨迹，建立了声压与声射线管横截面积之间的关系。1968 年，Randall［18］考虑分 层大气 的影响，提出了 大气物性沿高度变化对声爆特性影响的处理方法。1972 年，Thomas［19］以射线法为框架发展的波形参数法成为最为经典的快速声爆预测方法，但是该方法对压力波传播过程中的幅值变化仍采用线性化方法处理，导致波形预测结果仍不理想。Robinson［20］、Cleveland［21］、Crow［22］考虑到真实大气环境中的热黏性吸收和分子弛豫效应的影响，建立了基于广义Burgers 方程的远场声爆传播模型，能够更全面地涵盖大气的宏观和微观效应，得到了国内外的广泛应用。近年来，国内在声爆预测方面也积累了一些成果。王刚［23］和冷岩［24］等采用波形参数法进行了典型标模的数值预测和影响因素分析。张绎典［25］和乔建领［26］等使用广义Burgers方程法开展了声爆远场计算和应用研究。王迪等［27］开展了广义Burgers方程的高阶格式离散方法研究。总体来说，静止大气条件下的声爆预测方法发展的相对比较成熟，但上述方法并未考虑到大气湍流效应的影响，无法解释声爆波形受大气湍流作用产生的幅值和相位的随机性变化。

由于声爆波和大气湍流相关作用的过程本质上是湍流与激波等复杂结构在近无限大空间范围的多尺度相互作用问题。为了简化求解过程，研究人员试图从声爆波和大气湍流的物理特征出发，通过建立合理高效的数学模型来描述大气湍流导致的散射和衍射等复杂效应。其中，最具代表的是由Zabolotskaya和Khokhlov［28］及Kuznetsov［29］提出的KZK 方程（Khokhlov-Zabolotskaya-Kuznetsov Eqution），Piacsek［30］、Locey［31］和Luquet［32］发展的NPE 方程（Nonlinear Progressive Equation）以及Dagrau［33］等推导 的HOWARD 方 程（Heterogeneous One-Way Approximation for the Resolution of Diffraction ）等形式的多维弱非线性方程，它们能够较为准确地模拟声爆信号受湍流效应作用后产生的形态变化，并与飞行试验测量结果吻合较好。但是，基于多维弱非线性方程的预测方法仍面临物理机理不甚明确，模型发展不够成熟等问题，且对于声爆从机体产生到传播至地面的这类大尺度空间问题，求解多维弱非线性方程在计算量上依旧是难以接受的。1993 年，Robinson［34］提出了通过在射线法中添加随机湍流速度场的方法模拟大气湍流对声爆传播的影响。该方法思路简单易行，可以快速预测声爆的传播过程，且能够在一定程度上反映声爆波形幅值受湍流效应影响的变化规律。2009 年，Yamashita 等［35］通过该方法进一步探究了均匀各向同性大气湍流对声爆传播的不确定性影响，并与XB-70 飞机的飞行实测结果进行比较，获得了良好的一致性规律。2020 年，冷岩等［36］受Yamashita 工作的启发，将此方法应用于航空工业空气动力研究院（以下简称航空工业气动院）研发的ARI-Boom 声爆预测程序中，详细分析了大气湍流对地面声爆过压峰值、声射线传播路径及地面到达点位置的影响。然而上述研究工作主要是基于波形参数法开展的，该方法由于没有考虑到大气分子对声能量的吸收和弛豫作用，导致远场声爆波形过压峰值相对偏大且无法计算声爆的上升时间，尤其是在大气湍流的扰动下无法对地面声爆波形做出准确判断。相比之下，基于广义Burgers 方程的远场声爆预测方法［20-22，25-27］由于考虑到真实大气环境中的热黏性吸收和分子弛豫等物理效应，能够更加准确地反映大气对声爆传播特性的影响，目前还未见到基于广义Burgers 方程开展大气湍流影响的研究工作。

本文在前期波形参数法的基础上，开展了基于广义Burgers 方程的均匀各向同性大气湍流效应对声爆传播特性的影响研究。通过对中国航空工业空气动力研究院发展的ARI-Boom 声爆预测程序［24，27，36-37］进行改进，采用基于广义Burgers 方程的远场声爆预测方法，结合射线法框架，在其中添加随机湍流速度场，建立了一套可考虑热黏性吸收、分子弛豫等物理效应的大气湍流声爆影响快速预测方法，并分析了湍流强度、大气边界层厚度等大气湍流参数对典型远程超声速民机的声爆特性影响规律。

1 计算模型与方法

考虑大气湍流效应影响的声爆快速预测方法的整体研究思路可以分为3 个步骤：

步骤 1通过求解三维Navier-Stokes 方程获得飞机绕流流场的近场空间压力分布。使用中国航空工业空气动力研究院开发的ARI-Overset数值模拟软件［37-39］，采用AUSMDV 通量差分格式［40］求解Navier-Stokes方程得到声爆的近场压力分布，作为远场计算的输入信号。根据经验［37，39］，为了准确且高效地捕捉激波，在近场区域的计算中采用混合型计算网格，详细的参数设置及验证信息这里不再给出。

步骤 2在大气边界层内，通过von Karman能量谱［41］，采用离散Fourier 模态有限和生成均匀各向同性脉动速度场，模拟大气湍流对声爆传播过程中的速度扰动。von Karman能量谱可定义为

式中：k=(k1，k2，k3)为波矢量；L0=0.032 m 为长度尺度；η=0.005 为Kolmogorov 尺度。

能量谱确定后，采用N阶离散Fourier 模态有限和生成的随机速度扰动可表示为

式中：x=(x，y，z)为湍流场中给定点的位置坐标为振幅，kn为每一个波数向量kn对应的模量；Ψn为相位的随机数；σn为单位法向量。由前期经验［36］发现，离散Fourier 模态阶数N取200 时较为合适，其他变量的取值可参考文献［35-36］。由于大气边界层的厚度通常介于0.3～3 km 的范围内，模拟大气边界层的区域在x（飞行方向）、y（飞行侧向）、（z高度方向）3 个方向的尺度分别为6.0 km×0.24 km×3.0 km，并确保湍流场网格的位置能够覆盖声爆信号的传播路径。每个方向的网格点等距分布，间距为40 m，总网格量为67 500。

步骤 3以近场声爆压力特征作为初始条件，通过沿声射线求解广义Burgers 方程，模拟声爆信号从飞机近场穿过包含湍流边界层在内的真实大气环境直至传播到地面的整个过程。整体研究思路如图 2 所示。

图2 整体研究思路示意图Fig.2 Schematic diagram of overall research

无量纲化后的广义Burgers 方程可表示为［19］

式中：P为无量纲声压；σ和τ分别为无量纲声射线长度和波形延迟时间；ρ0和c0分别为环境大气密度和声速；S为声线管面积；Γ为热黏性吸收系数；θυ和Cυ分别为无量纲分子弛豫时间和系数。

采用算子分裂法将式（3）分解成5 项独立方程依次进行求解［27］。其中，前3 项方程在波形方向上分别使用五阶WENO 格式、四阶Pade 格式和Crank-Nicolson 格式，沿射线传播方向使用四步四阶Runge-Kutta 格式；后2 项方程具有解析解。声爆的传播路径及射线管面积可由几何声学的射线法［42］技术获得，其表达式为

式中：R(i)为射线路径；N(i)为波阵面单位法向量。增量ΔR(i)和ΔN(i)可通过式（6）和式（7）计算得到：

式中：V0(i)表示声爆信号所在高度的平均风速。在无湍流条件下，V0(i)=(V0x，V0y，0)仅需考虑两个水平方向的风速分量；在大气湍流的作用下，V0(i)=(V0x+u'，V0y+v'，w')，其 中u'、v'和w'表示大气湍流的速度分量，即步骤2 中利用von Karman 能量谱生成的随机速度扰动。

图3 给出了声爆信号沿声射线穿过湍流速度场的示意图，本文生成的随机湍流场主要作用于声射线传播路径和声线管面积，继而改变声爆特征的演化过程。由于离散的声射线位置与生成的湍流速度场不一定刚好在同一网格点上，需要通过三线性插值的方法将湍流速度分量插值到射线路径点。另外，声爆的传播时间相对于湍流速度的响应时间尺度来说非常小，因此可以假设生成的湍流场是冻结的。

图3 声爆信号穿过湍流场示意图Fig.3 Schematic diagram of sonic boom passing through turbulent field

2 计算结果与分析

本文使用的计算模型是航空工业气动院提出的一类远程超声速民机［43］，其采用大长细比、大后掠角箭形翼和鸭式气动布局方案，可有效降低超声速飞行条件下的激波阻力和声爆强度，几何外形如图 4所示。该模型的特征长度为67 m，飞行高度为15 km，巡航马赫数为1.6，飞行攻角为4°，提取离机体机身长度1.5 倍距离处（h/L=1.5）的空间压力分布曲线作为远场的输入波形，图 5 给出了该算例的近场压力云图和提取的空间压力分布曲线［43］。

图4 一类典型远程超声速民机构型［43］Fig.4 Illustration of supersonic civil aircraft configuration［43］

图5 近场压力云图及空间压力分布曲线［43］Fig.5 Near-field pressure contour and pressure distribution［43］

图6 给出了使用广义Burgers 方程法与波形参数法在无湍流条件下计算得到的地面声爆波形对比结果，取地面反射因子为1.9。可以发现，基于广义Burgers 方程的远场声爆预测方法由于考虑到大气分子对于声能量的热黏性吸收和分子弛豫作用，得到了更符合物理特性的声爆信号，即波形中的激波不再是理想化的强间断，而是有限厚度的强压缩波，波形的最大过压峰值也相对有所下降。另外，与波形参数法相比，广义Burgers 方程模型能够较准确地计算声爆的上升时间，在本算例中约为1.80 ms。上升时间是计算声爆感觉噪声级的重要参数，在相同的过压峰值下，上升时间越短，声爆的感觉噪声级越高，人们感受到的声音越尖锐。

图6 无湍流条件下广义Burgers 方程法与波形参数法计算结果对比Fig.6 Comparison of calculation results of augmented Burgers equation and waveform parameter method without turbulence

图7 为湍流条件下，分别采用广义Burgers 方程法和波形参数法获得的地面声爆波形。每种方法分别选取100 组计算结果作为统计样本，这里假设湍流边界层厚度为1.5 km，生成的湍流随机速度场的均方根速度值Vrms约为2.5 m/s。其中，黑色曲线表示不受湍流影响的地面声爆特征；绿色和红色曲线分别表示声爆过压峰值受湍流效应影响后出现增强或减弱的现象。根据统计结果表明，采用广义Burgers方程法得到的地面声爆过压峰值范围介于19.44～41.85 Pa 之间；而波形参数法的过压峰值范围介于12.89～228.03 Pa 之间。显然，基于广义Burgers方程的预测方法由于能够表征热黏性吸收、分子弛豫等大气物理效应，其过压峰值的波动范围处于合理的区间内，与NASA 的飞行试验测量结果规律相一致［11］。由于波形参数法对声爆传播过程中的幅值变化采用了线性化假设，通过该方法获得的声爆过压峰值不够准确，尤其是经过大气湍流的干扰后，其变化趋势与实际的飞行试验测量结果有所偏差，一般情况下均有明显的过度放大。因此，本文建立的方法能够更加真实地反映大气湍流对声爆传播特性的影响。

图7 湍流条件下广义Burgers 方程法与波形参数法计算结果对比Fig.7 Comparison of calculation results of augmented Burgers equation and waveform parameter method with turbulence

另外，相比于前期的典型超声速公务机［36］，本文采用的远程超声速民机声爆波形更加复杂，该方法依然能给出复杂波系的大气湍流影响规律。虽然广义Burgers 方程仍是一维标量方程，但是其能够在一定程度上更加合理地模拟大气湍流对声爆特征幅值产生的随机性影响。同时，与KZK［28-29］、NPE［30-32］和HOWARD［33］等多维方程模型相比，一维广义Burgers 方程具有计算量小的优势，能够快速判断大气边界层内的湍流效应对声爆特性的影响规律，在未来超声速民机飞行轨迹规划的早期阶段，可方便快速地预估地面声爆强度受湍流影响的阈值。

2.1 湍流强度对声爆特性的影响

为了探究大气湍流强度对地面声爆特性的影响规律，假设大气边界层厚度为1.5 km，通过设置随机湍流速度场的均方根速度值Vrms分别为1.0、1.8、2.5 m/s，模拟低、中、高3 种湍流强度的影响。图 8 为3 种湍流条件下的地面声爆特征对比图，每组条件分别统计了100 组计算样本。可以发现，随着湍流强度的提升，声爆的过压峰值范围及最大过压峰值均显著增大，而最小过压峰值有变得更低的趋势，说明地面声爆特征发生的随机性变化会更加明显，具体的统计结果如表 1所示。

图8 不同湍流条件下的地面声爆特征Fig.8 Ground sonic boom signatures for different turbulence intensities

表1 不同湍流强度下地面声爆特征统计Table 1 Distribution statistics of ground sonic boom signatures for different turbulence intensities

图9 给出了不同湍流强度下地面声爆过压峰值的概率分布柱状图。对于低湍流强度条件，大部分统计结果与无湍流条件下声爆过压峰值的比值（Δpmax-tur/Δpmax-no-tur）主要集中在0.9～1.1 之间，说明大气湍流对声爆传播过程中的扰动较弱，因此地面声爆特征的幅值变化幅度较小。当湍流强度持续增大后，声爆的过压峰值不再稳定在较小的区间，根据统计，Δpmax-tur/Δpmax-no-tur介于0.9～1.1 之间的概率分别为93%、66%和38%，声爆过压峰值的变化范围有进一步扩大的趋势，可见高湍流强度会使声爆强度增强的概率显著提升。当湍流扰动速度达到2.5 m/s时，湍流条件下的最大声爆过压峰值达到无湍流条件下的1.54倍。

图9 不同湍流强度下地面声爆过压峰值的概率分布柱状图Fig.9 Ground peak overpressure probability distribution histogram for different turbulence intensities

图10 给出了所有条件下声爆到达点的位置分布。可以观察到，湍流条件下的声爆到达点位置（黑色方块）随机分布在无湍流条件下声爆到达点的周围（红色圆形，x=14.10 km，y=0 km处），且随着湍流强度的增大，声爆到达点围成的范围也呈现出扩大的趋势。其中，涵盖95%到达点围成的圆的半径分别为1.65 m、3.05 m 和3.55 m，说明高湍流强度不仅可能会增强对声爆特性产生的随机性变化，还有可能使声爆到达点的位置更加分散，继而影响声爆毯的面积。表 2列出了更详细的统计结果，声爆到达点沿飞行方向（x方向）的波动范围要明显大于沿飞行侧向（y方向）的波动范围，该规律与波形参数法［35-36］的预测结果相一致。但总体上达到点分布区域仍在比较小的范围，说明大气湍流对声爆到达点影响整体较小。

图10 不同湍流条件下声爆到达点位置分布Fig.10 Distribution of sonic boom arrival points on ground for different turbulence intensities

表2 不同湍流强度下声爆到达点位置分布统计Table 2 Distribution statistics of sonic boom arrival points for different turbulence intensities

2.2 湍流边界层厚度对声爆特性的影响

气象观测结果表明，地球表面的湍流边界层厚度通常介于0.3～3 km 的范围，会伴随地形和气象条件有所变化。为了探究湍流边界层厚度对声爆特性影响的一般性规律，本节假设随机湍流场的均方根速度值Vrms为2.5 m/s，设置边界层厚度z从0.5 km 依次增加到3 km，每个状态条件下同样统计了100 组数据样本，计算得到的地面声爆特征如图 11 所示。很明显，声爆特征的随机性变化范围随着湍流边界层厚度的增加有扩大的趋势。另外，根据表 3 的统计结果显示，声爆特征增强或减弱的概率与湍流边界层厚度没有明显的规律可循，而声爆过压峰值的标准差持续增大，说明声爆幅值的波动程度将更加剧烈。

表3 不同边界层厚度的地面声爆特征统计Table 3 Distribution statistics of ground sonic boom signatures for different turbulence boundary layer heights

图11 不同边界层厚度的地面声爆特征Fig.11 Ground sonic boom signatures for different turbulence boundary layer heights

图12 给出了不同湍流边界层厚度条件下，地面声爆过压峰值与无湍流条件的对比结果。对于4 种不同厚度的边界层，Δpmax-tur/Δpmax-no-tur介于0.9～1.1 之间的概率分别为81%、61%、44% 和38%，该变化规律与图 11 的观察结果相一致，进一步说明随着边界层增厚，大气湍流效应对声爆特性产生的随机性影响更加明显，其最大过压峰值将逐渐升高。当z=3 km 时，声爆到达地面的最大过压峰值达到了无湍流条件下的2.04 倍，这将对人们的正常生活产生严重的影响，因此在飞行轨迹规划时应尽量避免湍流边界层较厚的区域。

图12 不同湍流边界层厚度下地面声爆过压峰值的概率分布柱状图Fig.12 Ground peak overpressure probability distribution histogram for different turbulence boundary layer heights

图13 给出了4 种湍流边界层厚度对应的声爆地面到达点的位置分布。根据表4，当边界层厚度为0.5 km 时，有95%的声射线落点都在以无湍流条件下声爆到达点的位置为圆心、半径为2.30 m 的圆内。当湍流边界层厚度依次增加到1，2，3 km 时，其涵盖95%声爆到达点围成的圆的半径依次增加至3.61，5.09，6.70 m，说明随着湍流边界层厚度的增加，声爆到达地面的位置同样呈现出向周围分散的趋势，但达到点仍分散在比较小的范围。

表4 不同湍流边界层厚度下的声爆到达点统计Table 4 Distribution statistics of sonic boom arrival points on ground for different turbulence boundary layer heights

图13 不同湍流边界层厚度下的声爆到达点位置分布Fig.13 Distribution of sonic boom arrival points on ground for different turbulence boundary layer heights

3 结论

基于von Karman 能量谱，采用离散Fourier模态有限和生成均匀各向同性大区湍流随机速度场，结合广义Burgers 方程远场声爆预测方法，建立了可考虑热黏性吸收、分子弛豫等物理效应的大气湍流声爆影响快速预测方法，并开展了典型远程超声速民机声爆传播特性的大气湍流效应影响研究。得到的主要结论如下：

1）相比于前期基于波形参数法框架的声爆预测方法，建立的方法能够更合理地表征热黏性吸收和分子弛豫等大气物理效应，可以更真实地反映大气湍流对声爆传播特性的影响规律，该方法的计算结果也更符合真实飞行试验测量结果的变化趋势。

2）使用的算例相比于前期的典型超声速公务机生成的声爆波形更加复杂，但该预测方法仍能给出复杂波系的大气湍流影响规律；同时，该方法可以对地面声爆信号进行快速预测，在未来超声速民机的飞行轨迹规划中具有独到优势。

3）湍流强度和大气边界层厚度是影响声爆特性的重要因素，随着湍流强度和边界层厚度的增加，声爆过压峰值的波动范围有明显增大的趋势。在极端条件下，声爆信号受大气湍流影响后的最大过压峰值可能达到无湍流条件下的2 倍以上。因此，在规划飞行路线时，应该仔细评估当地湍流强度和大气湍流边界层厚度的影响。

4）受大气湍流效应的影响，声爆到达地面的位置呈现出随机性分布，且随着湍流强度和边界层高度的增加更加分散；同时，声爆到达点沿飞行方向的波动范围要明显大于沿飞行侧向的波动范围，该规律与波形参数法的预测结果相一致；但总体上达到点分布区域仍在比较小的范围，说明大气湍流对声爆到达点影响整体较小。