0 引言

声爆是超声速民机实现商业运营需要解决的核心问题之一。近些年来，国内外科研工作者对声爆预测与低声爆设计的研究逐渐增多。美国航空航天协会(AIAA)先后组织了两届声爆预测研讨会^[1]，对现有的声爆预测技术进行了评估。Haas等提出了低声爆超声速客机多激波反设计方法^[2]。朱自强等对现有的降低声爆的优化设计方法进行了概述，提出以声爆最小为目标的高效、鲁棒的全机数值优化软件将是CFD研究者努力的目标^[3]。Kirz运用代理优化方法完成了AIAA声爆预测大会标模的低声爆优化设计^[4]。冯晓强等对混合网格的近场预测以及基于波形参数法的远场预测相结合的声爆预测方法进行了研究，建立了基于声爆预测方法、响应面模型和Pareto遗传算法的声爆/气动一体化设计方法^[5]。乔建领等将最新发展的代理优化算法与声爆预测方法相结合, 发展了一种具有全局优化能力的高效低声爆优化设计方法^[6]。这些工作都是着眼于优化设计方法，而鲜有考虑设计目标对于低声爆优化的影响。郝璇等开展了某公务机低声爆气动布局优化设计研究，得到同时降低地面声爆过压值与A计权声级的优化结果，两者存在一定的差异^[7]。与A计权声级相比，本文选择的史蒂文斯(Stevens)响度方法对地面人员声爆感受的描述更为详细具体，可以在声爆预测中更好地体现声爆的噪声水平。

人对于声音的第一感觉是响度，它与人的听觉器官和听觉神经对声音强度和频率响应有关。声音强度可以用声压级L_p体现，即人耳的听觉相应正比于声能量的对数值, 单位是分贝(dB)，声压级的计算方法见公式(1)，其中p₀=2×10^-5 Pa。

$ {L_p} = 20{\rm{lg}}(p/{p_0}) $

(1)

观察者在地面上接受到的远处传来的声爆，由于对不同频率上发出的相同强度的声爆没有统一的响应，因而用声压级作为声爆大小的评价尺度是不合适的。由此引入声爆响度级概念，采用总响度级作为衡量地面声爆强度的参数指标。

常用的噪声响度级指标主要包括感知噪声级(PNL)，A加权(SEL(A))、C加权(SEL(C))和线性(无加权)(SEL(U))声音曝光级，史蒂文斯响度级(PL)，茨维克(Zwicker)响度级(LLZd和LLZf)等^[8]。在提出的多种响度级评价量中，不同指标的差异主要是加权方式的不同，其中史蒂文斯法^[9]与茨维克法^[10]在国内外研究中较为常用。两种方法计算结果相近，但与后者相比，史蒂文斯提出的Mark Ⅶ模型由于计算方法更为简洁，因此选择史蒂文斯响度法作为声爆总响度级的计算方法。

本文将史蒂文斯响度方法集成至超声速飞机低声爆优化设计平台，对Seeb-ALR模型进行了低声爆优化设计，与常规近场最大过压值方法的优化结果进行对比，分析两者在模型轮廓与远近场过压值分布形态的差异，验证史蒂文斯响度方法在该超声速民机低声爆设计平台中的实用性。

1 史蒂文斯响度计算方法

Mark Ⅶ模型考虑到不同噪声频率之间的遮蔽效应，结合声爆特性，得到一组等响度曲线(不同频率具有同等响度声音的声音能量水平的曲线)如图 1所示。该方法认为响度指数最大的频带贡献最大，而其他频带由于前者的遮蔽而贡献减少，Mark Ⅶ模型设定标准参考信号为1/3倍频带3150 Hz的中心频率，取该中心频率上32 dB所产生响度为基本度量单位宋(Sone)(单位宋，参考声压级32 dB的响度或噪度)，之后利用等响度指数曲线将三分之一倍频程声压级转换为各频带响度，求出响度级(PLdB)。信号声压级每提高9 dB，响度将增加一倍。

在计算总响度时，乘上一个小于1的修正系数F(带宽因子)得到结果，F与响度有关，带宽因子随响度变化曲线如图 2所示。

史蒂文斯响度级方法的计算流程如下：首先将已经得到的地面声爆时域信号进行快速傅里叶变换(FFT)，将时域结果转换到频域上，然后得到每个频带上的声压级(倍频程或三分之一倍频程)，根据各频带中心频率的声压级，采用响度指数曲线确定各频带响度S_i，由公式(2)计算总响度S_t。

$ {S_t} = {S_m} + F(\sum {S_i} - {S_m}) $

(2)

其中:S_m为响度指数中最大值；S_i为第i个频带的响度指数。确定各频带响度及总响度后，最后通过公式(3)求出总响度级P。

$ P = 32 + 9{\rm{lo}}{{\rm{g}}_2}{S_t} $

(3)

2 低声爆优化设计平台

将史蒂文斯响度方法引入超声速民机低声爆优化设计平台，基于Mark Ⅶ模型构建远场声爆响度计算模块。集成后考虑响度的低声爆优化设计平台框架如图 3所示。

该平台基于Matlab语言开发，可以适用于超声民机低声爆设计，具有包含远近场声爆计算、模型参数化、远场声爆响度计算等功能模块。

2.1 远近场混合的声爆计算方法

声爆从飞机传播到地面距离较远，在单一计算域下进行全场CFD求解获得地面声爆信号难度较大。Hayes等提出了基于近场流场精确模拟和远场传播模型相结合的混合预测方法^[11]，其核心思想是将流场划分为近场与远场两个计算域(见图 4)，首先通过风洞试验或CFD数值模拟，获得飞机近场压力分布，声爆的远场传播过程使用远场传播模型来计算。

近场流场计算采用二阶迎风格式的有限体积法求解Navier-Stokes控制方程。在广义坐标下，三维非定常可压缩Navier-Stokes方程积分形式如下：

$ \frac{\partial }{{\partial t}}\int\limits_\mathit{\Omega } {\mathit{\boldsymbol{Q}}{\rm{d}}\mathit{\Omega }} + \oint\limits_{\partial \mathit{\Omega }} {({\mathit{\boldsymbol{F}}_c} - {\mathit{\boldsymbol{F}}_v})} \cdot \mathit{\boldsymbol{n}}{\rm{d}}\mathit{\boldsymbol{S}} = \int\limits_\mathit{\Omega } {{G_s}{\rm{d}}\mathit{\Omega }} $

(4)

其中$\mathit{\boldsymbol{Q}} = {\left[{\rho \; \rho u\; \rho v\; \rho w\; \rho E\; \rho \overline v } \right]^{\rm{T}}}$为守恒变量矢量，ρ为流体单元密度，u, v, w为速度分量，E为单元总内能和，${\overline v }$为湍流模型变量，n为控制体表面法向量，F_c和F_v分别为对流通量与黏性通量，G_s为湍流模型源项。根据不同计算精度与计算时间的需求，人工选取合适的求解器与求解方法对流场进行数值模拟。

远场传播计算采用波形参数方法得到地面声爆波形^[13]，该方法将风洞试验或CFD数值模拟得到的近场声爆波形作为输入，将近场声爆波形离散为若干线性信号段，每个信号段采用一组参数描述其在大气中传播，通过几何声学方法建立声爆信号大气传播波动方程并进行求解，在给定大气条件的情况下，得到整个传播路径上任意高度的声爆特征参数。波形参数法基本原理如图 5所示。

声爆信号的传播波动方程如下：

$ \frac{{{\rm{d}}{m_i}}}{{{\rm{d}}t}} = {C_1}{m_i}^2 + {C_2}{m_i} $

(5)

$ \frac{{{\rm{d}}\Delta {p_i}}}{{{\rm{d}}t}} = \frac{1}{2}{C_1}\Delta {p_i}({m_i} + {m_{i - 1}}) + {C_2}\Delta {p_i} $

(6)

$ \frac{{{\rm{d}}{\lambda _i}}}{{{\rm{d}}t}} = - \frac{1}{2}{C_1}(\Delta {p_i} + \Delta {p_{i + 1}}) - {C_1}{m_i}{\lambda _i} $

(7)

$ \begin{array}{l} \Delta {p_i} = {p_i} - {p_{i - 1}}\\ {\lambda _i} = {T_{i + 1}} - {T_i} \end{array} $

(8)

其中:m_i为第i段信号压强随时间的变化率，Δp_i为第i段信号通过激波的压强增量，λ_i为第i段信号的持续时间。

2.2 参数化方法与优化方法

参数化方法是采用数学函数定义一组可以控制几何外形的特征参数，对设计对象进行数学建模，从而通过调整这些特征参数修改几何外形。简单的曲线参数化可通过样条参数化方法，即Bezier曲线等直接描述模型外形。对于复杂外形，可采用CST(Class Shape function Transformation)方法，该方法采用类函数(Class function)与形函数(Shape function)相结合的方式进行几何外形描述，数学表达式如下：

$ \xi \left( \psi \right) = C_{N1}^{N2}\left( \psi \right)S\left( \psi \right) $

(9)

其中:C_N1^N2(ψ)为类函数，用来表达参数化对象几何形状分布解析函数；S(ψ)为型函数，形式为线性加权的Bernstein多项式，用来对外型进行调整。

优化方法部分采用遗传算法(GA)与粒子群优化算法(PSO)开展超声速民机的低声爆优化设计。与梯度优化算法相比，两种优化方法在超声速民机低声爆设计中均可以通过计算适用度来同时比对全局多个点，对于初始值依赖小^[14-15]。两者在寻优路径上各有优劣，可以根据实际问题需求选取合适的优化方法。

3 方法应用

基于改进后的超声速民机低声爆优化设计平台，分别选取史蒂文斯响度法计算得到的地面声爆总响度级与近场最大过压值为目标进行优化设计。

3.1 模型描述

本文计算采用的Seeb-ALR模型是一个轴对称、非升力体模型^[16]，主要用于声爆风洞试验的测量标定，具有较为详细的声爆风洞试验数据，因此被选作第一届AIAA声爆预测会议的研究模型之一^[17]。其模型几何构型如图 6所示。

3.2 数值验证

Seeb-ALR模型几何尺寸参考长度L=449.1 mm，初始马赫角为38.68°，计算迎角α=0°，来流马赫数Ma=1.6，计算雷诺数Re=6.42×10⁶。

计算网格采用六面体网格，网格总单元数为80万，如图 7所示。采用SA湍流模型，空间离散格式采用二阶迎风格式，进行定常计算迭代2000步至流场收敛。

图 8为对称面上的过压分布云图，可以看出锥体前缘头激波，膨胀波，尾激波等主要的流场分布特征。

在沿锥体轴线垂直距离h=1L处提取近场过压信号，与试验结果^[17]进行对比, 如图 9所示。可以看出，计算得到的前缘头激波上升位置、压力膨胀段、尾激波位置与试验结果均吻合良好，尾激波强度预测略大于试验结果，主要原因是由于试验支架造成测量结果有所偏差。结果表明，模型采用的CFD流场计算方法是可靠的，计算结果满足优化设计所需要求。

Seeb-ALR模型作为风洞试验模型, 在计算远场传播时将其等比放大到一个真实飞机尺寸l=51.816 m, 并给定传播距离H=15743 m，进而计算远场的声爆响度级。

3.3 几何参数化

Seeb-ALR模型几何构型比较简单，因此采用Bezier曲线进行参数化建模。选取6个控制点坐标值如表 1所示。

控制点分布如图 10所示，其中点1、点5与点6为固定点，点2、点3与点4在上下限内随机取值，构成计算样本区域。对得到的控制点构成Bezier曲线得到锥形轮廓线，旋成新的锥模型进行自动网格生成并完成流场计算。

3.4 设计方案 3.4.1 常规近场最大过压值优化设计

常规优化设计采用两轮优化，采用遗传优化算法，以近场最大过压值为优化目标，以控制点横纵坐标值为优化变量，约束值如表 1所示，得到最优方案之后，结合第一轮最优方案与最初控制点范围，确定第二轮控制点变化范围进行进一步优化。

两轮优化过程如图 11与图 12所示。经过两轮优化，近场最大过压值显著减少。表 2显示了经过前后两轮优化设计，最大过压值相比原模型降低了40.2%。

3.4.2 远场响度级优化设计

与常规近场优化设计相比，远场响度级优化采用同样的优化算法、优化变量与约束，以远场响度级为优化目标，其优化过程如图 13所示。对远场响度级优化显著降低了远场声爆响度级的大小。表 3展示了地面声爆响度级相比原模型下降了2.2PLdB。

3.5 结果对比

两种低声爆设计方案在各自的优化目标下都有效地降低了目标值。选取两种方案得到的最优结果进行对比，分析两种不同设计目标下结果的差异。

两种方案得到最优结果的锥形轮廓如图 14所示。可以看出，两种设计方案结果与原模型相比都在一定程度上削减了前缘厚度，地面声爆响度级方案获得的锥体轮廓线更为尖锐。

两种方案得到的近场过压分布曲线与地面过压分布曲线分别如图 15与图 16所示。可以看出：两组设计方案结果都有效降低了初始模型的最大过压值，使最大过压值相比原模型分别降低了40.6%和18.4%；前缘形状也同时影响了远场传播地面声爆信号形状，使得不同的设计目标得到的地面声爆信号分布差别很大；优化结果都在一定程度上推迟了尾激波的位置。

两种方案地面响度级对比如图 17所示，其使地面声爆响度级相比原模型分别降低了1.4 PLdB和2.2 PLdB。以地面史蒂文斯总响度级为优化目标得到的优化结果在显著降低声爆响度级的同时，也降低了最大近场过压值。

4 结论

本文将史蒂文斯的Mark Ⅶ模型集成进超声速民机低声爆优化设计平台，对Seeb-ALR模型进行了低声爆优化设计，并与常规近场最大过压值方法得到的优化结果进行对比，得到以下结论：

1) 该超声速民机低声爆设计平台实现了史蒂文斯响度级方法的应用，可以通过响度级反映对地面人员的影响，并将其应用于低声爆设计。

2) 两种设计结果都有效地减少了近场最大过压值与地面响度级，不同的设计目标得到的锥体轮廓线与地面声爆信号分布差别很大。

3) 常规近场过压值优化方案，虽然也降低了地面声爆声压级，但是与直接对地面声爆声压级优化方案得到的结果仍有差距。

4) 地面声爆声压级优化方案，对近场过压最大值影响较小，说明降低近场过压最大值并不是有效降低地面声爆的方式，通过修型改变近场过压值分布可以得到更低的地面声爆声压级。