文章快速检索     高级检索
  空气动力学学报  2019, Vol. 37 Issue (4): 663-674  DOI: 10.7638/kqdlxxb-2018.0267

引用本文  

乔建领, 韩忠华, 丁玉临, 等. 基于广义Burgers方程的超声速客机远场声爆高精度预测方法[J]. 空气动力学学报, 2019, 37(4): 663-674.
QIAO J L, HAN Z H, DING Y L, et al. Sonic boom prediction method for supersonic transports based on augmented Burgers equation[J]. Acta Aerodynamica Sinica, 2019, 37(4): 663-674.

基金项目

国家自然科学基金(11772261);航空科学基金(2016ZA53011);CFD前沿技术项目(2015-F-016)

作者简介

乔建领(1994-), 男, 山东菏泽人, 博士研究生, 研究方向:超声速客机声爆预测与低声爆优化设计.E-mail:qiaojianling@163.com

文章历史

收稿日期:2018-11-12
修订日期:2019-05-09
基于广义Burgers方程的超声速客机远场声爆高精度预测方法
乔建领1,2 , 韩忠华1,2 , 丁玉临1,2 , 池江波1,2 , 孟冠宇1,2 , 宋笔锋1 , 宋文萍1,2     
1. 西北工业大学 航空学院 超声速客机研究中心, 西安 710072;
2. 西北工业大学 航空学院 气动与多学科优化设计研究所, 西安 710072
摘要:声爆高精度预测技术是新一代环保型超声速客机设计的核心关键技术之一。基于广义Burgers方程发展了可考虑"大气风"效应的远场高精度预测方法,开发了声爆预测程序"bBoom",并研究了近场声爆信号提取位置和"大气风"对远场声爆计算结果的影响。首先,给出了广义Burgers方程和声爆传播射线的计算方法,重点讨论了方程中计算分子弛豫效应和热黏吸收效应的关键参数。其次,通过简单轴对称构型标模算例、NASA的C25D构型和洛马的LM1021构型等复杂超声速客机算例对所发展的方法进行了验证,表明本文发展的方法在预测远场声爆时具有较高可信度。最后,基于所发展的方法对比了由不同近场提取位置传播到地面的波形,研究了飞机向不同方向飞行时,"大气风"对地面声爆强度和地面影响域的影响。结果表明:对于类C25D标模构型,为了确保远场声爆预测结果具有较高精度,应取机身下方约3倍机身长度位置处的近场压强信号作为传播方程的输入;另外,"大气风"会影响地面声爆强度及地面影响域,在预测时有必要加以考虑。
关键词超声速客机    声爆预测    广义Burgers方程    声射线方法    大气风    
Sonic boom prediction method for supersonic transports based on augmented Burgers equation
QIAO Jianling1,2 , HAN Zhonghua1,2 , DING Yulin1,2 , CHI Jiangbo1,2 , MENG Guanyu1,2 , SONG Bifeng1 , SONG Wenping1,2     
1. Supersonic Transport Research Center, School of Aeronautics, Northwestern Polytechnical University, Xi'an 710071, China;
2. Institute of Aerodynamic and Multidisciplinary Design Optimization, School of Aeronautics, Northwestern Polytechnical University, Xi'an 710071, China
Abstract: The high-fidelity sonic-boom prediction is crucial for the design of low-boom supersonic transport aircraft. In this article, a far-field sonic-boom prediction method solving the augmented Burgers equation is studied, and a prediction code called "bBoom" is developed. First, the numerical method for solving augmented Burgers equation and a ray tracing method in propagation are studied. Some significant parameters of molecular relaxation processes and thermoviscous absorption in the Burgers equation are presented. Second, some benchmark cases from the 2th AIAA Sonic Boom Prediction Workshop, including a simple axi-symmetric body, a NASA low-boom configuration, and LM1021 configuration, are employed to validate the code. The results show that the developed code is accurate and reliable for predicting the far-field waveforms. At last, far-field sonic boom waveforms propagated from different near-field pressure signals at different locations are compared, and the distinction of waveforms in different flight azimuthal angles with atmospheric winds is investigated. It is shown that, in order to improve the accuracy of predicting the sonic boom for configurations like a C25D model, the near-field signal should be extracted at a distance of three body lengths off the aircraft. It is also shown that atmospheric wind profiles have a significant impact on the sonic-boom waveform and sonic-boom carpet on the ground.
Keywords: supersonic transport    sonic boom prediction    augmented Burgers equation    ray tracing    atmospheric wind profiles    
0 引言

超声速客机是未来民机发展的重点方向之一。早在20世纪中后期,英法和前苏联就研制了以“协和”和“图-144”为代表的第一代超声速客机。然而其商业运营以失败告终,究其原因,严重的声爆问题是造成其最终停止商业运营的重要原因之一。在吸取第一代超声速客机失败教训的基础上,近年来,美国、欧洲和日本等掀起了新一代环保型超声速客机的研究热潮,并制定了一系列研究计划[1-3](如图 1为美国波音公司的“N+2”代超声速客机方案),重点突破以声爆预测及其抑制为核心的一系列关键科学与技术问题。


图 1 波音公司拟研发的下一代超声速客机[4] Fig.1 The next generation of supersonic transport proposed by Boeing[4]

声爆问题是制约新一代超声速客机研制的核心关键问题。远场声爆的高精度预测是新一代超声速客机设计中声爆评估和优化设计的基础和前提,是其中核心的关键技术之一。目前,国内外发展的声爆预测理论和方法主要有修正线化理论[5-7]、简化声爆预测方法[8-9]、CFD全流场求解方法[10-12]、近场CFD计算耦合远场传播方法[13-16]。由于受到计算量和计算精度的限制,近场CFD与远场传播相结合的方法,成为了最主要的声爆高精度预测方法。其中,远场传播方法主要有波形参数法和广义Burgers方程法。基于广义Burgers方程的方法能够较准确地计算声爆上升时间,因而目前已经成为远场传播的主要方法。

广义Burgers方程在声爆预测领域的应用可追溯到20世纪80年代。1981年,Pierce[17]推导了考虑分子弛豫、热黏效应的广义Burgers方程,并在频域中进行求解。1991年,Kang[18]和Robinson[19]分别在博士论文中通过求解广义Burgers方程,研究了分子弛豫效应、大气分层及风梯度对远场波形的影响,并且求解过程只在时域中进行,不需要频域和时域之间的转换,减少了计算量和数值误差。1993年,Robinson[20]基于考虑分子弛豫、大气分层及风梯度影响的广义Burgers方程,发展了第一个非线性声爆预测程序“ZEPHYRUS”,该程序能够较好地计算波形结构和上升时间。1995年,Cleveland[21-22]采用算子分裂法对大气中声爆的传播过程进行模拟,研究了分层大气和几何声学扩散对远场声爆上升时间的影响。目前,广义Burgers方程已经考虑分子弛豫、热黏效应、大气分层、几何声学和经典非线性的影响。国际上大部分高精度的远场声爆预测程序都是基于该方程进行求解的。例如,2011年NASA研究者Rallabhandi[23]开发的“sBOOM”程序,2012年日本研究者Yamamoto等[24-25]开发的“Xnoise”程序等,均采用了Cleveland的方法。

国内从2009年开始,西北工业大学宋笔锋教授团队开展了修正线化理论、简化声爆预测方法、波形参数法的声爆预测研究[26-30]。近年来,随着声爆高精度预测方法的发展,国内开展了一些近场声爆CFD计算方面的研究[31-35],基于非线性波动方程(nonlinear wave propagation equation)[36]和广义Burgers方程的远场声爆预测方法研究[37],其中基于广义Burgers方程的方法已经可以考虑分子驰豫、热黏性、分层大气和几何声学因素的影响。但是,为了获得高精度的远场声爆预测结果,需要进一步研究近场数据提取位置对远场波形的影响和“大气风”对远场声爆强度及地面影响域的影响。一方面,近场数据提取位置关系到远场声爆的计算结果和近场计算域的选择,涉及到近场计算量的大小,进而影响声爆的高精度评估速度。另一方面,在真实大气中,“风”是普遍存在的,它会改变声爆的原始传播路径,对远场声爆强度和地面影响域有重要影响,因此在更接近真实大气条件的声爆预测中,还需要进一步研究“风”的影响。

本文开展了针对广义Burgers方程的离散数值求解研究,发展了一套能够考虑“大气风”影响的远场声爆计算方法,并开发了声爆高精度预测程序“bBoom”。采用美国航空航天学会(AIAA)第二届声爆预测研讨会(SBPW-2)提供的标模算例,验证了方法和程序的正确性和精度。在此基础上,开展了近场压强提取位置和“大气风”对远场波形及影响域的影响研究。本文研究工作可作为进一步开展考虑真实大气效应的高精度远场声爆预测、以及低声爆高升阻比超声速客机设计与优化研究的基础。

1 广义Burgers方程及其数值求解 1.1 广义Burgers方程

考虑几何声学、大气分层、热黏吸收和分子弛豫的影响,无量纲化的广义Burgers方程[23]可表示为:

$ \begin{array}{l} \frac{{\partial P}}{{\partial \sigma }} = - \frac{1}{{2A}}\frac{{\partial A}}{{\partial \sigma }}P + \frac{1}{{2{\rho _0}{c_0}}}\frac{{\partial \left( {{\rho _0}{c_0}} \right)}}{{\partial \sigma }}P\\ + P\frac{{\partial P}}{{\partial \tau }} + \frac{1}{\mathit{\Gamma }}\frac{{{\partial ^2}P}}{{\partial {\tau ^2}}} + \sum\limits_j {\frac{{{C_j}}}{{\left( {1 + {\theta _j}\frac{\partial }{{\partial \tau }}} \right)}}\frac{{{\partial ^2}P}}{{\partial {\tau ^2}}}} \end{array} $ (1)

式中各个变量的说明如表 1所示。右边各项的物理含义分别为:第一项是几何声学中能量的重新排布对声爆传播的影响;第二项是大气分层对声爆传播的影响,它和第一项合称为几何声学效应,通常根据Blokhintzev声学不变量进行求解,在考虑“大气风”时,其表达式为式(2);第三项是经典Burgers方程的非线性项;第四项是热黏吸收效应的影响;第五项是各分子弛豫过程对声爆传播的影响。

$ P\sqrt {\frac{{c_n^2A}}{{{\rho _0}c_0^3}}} = {\rm{const}} $ (2)
表 1 方程(1)中各变量的说明 Table 1 Explanation of variables in equation (1)

其中,cn=c0+ W · NW为当地大气风速, N为波面单位法向量。

分子弛豫效应和热黏吸收效应在声爆传播计算中具有重要的作用,下面将详细介绍这两个效应涉及的关键参数。

(1) 分子弛豫效应

式(1)分子弛豫效应项中,需要确定小信号下的声速增量(Δc)j[38]和相应的弛豫时间τj[39]

声速增量(Δc)j的近似表达式如下:

$ {\left( {\Delta c} \right)_j} = \frac{{{c_0}}}{2}\frac{{{{\left( {\gamma - 1} \right)}^2}}}{\gamma }\frac{{{n_j}}}{n}{\left( {\frac{{T_j^ * }}{{{T_j}}}} \right)^2}{{\rm{e}}^{ - T_j^ * /{T_j}}} $ (3)

式中:nj/n为第j个弛豫分子占分子总量的比例;Tj*为分子振动的特征温度;Tj为当前分子振动的温度。在只考虑氧气O2和氮气N2分子弛豫效应的情况下,nO2/n=0.21,nN2/n=0.78,T*O2=2239 K,T*N2=3352 K。

不同分子的弛豫时间τj不同,根据半经验公式,有:

$ \left\{ \begin{array}{l} {f_{r,{{\rm{O}}_2}}} = \frac{{{p_0}}}{{{p_s}}} \times \left( {24 + 4.04 \times {{10}^4}h\frac{{0.02 + h}}{{0.391 + h}}} \right)\\ {\tau _{{{\rm{O}}_2}}} = \frac{1}{{2{\rm{ \mathsf{ π} }}{f_{r,{{\rm{O}}_2}}}}} \end{array} \right. $ (4)
$ \left\{ \begin{array}{l} {f_{r,{{\rm{N}}_2}}} = \frac{{{p_0}}}{{{p_s}}}{\left( {\frac{{{T_{{\rm{ref}}}}}}{{{T_0}}}} \right)^{1/2}} \times \left( {9 + 280h \times {{\rm{e}}^{ - 4.17\left[ {{{\left( {{T_{{\rm{ref}}}}/{T_0}} \right)}^{1/3}} - 1} \right]}}} \right)\\ {\tau _{{{\rm{N}}_2}}} = \frac{1}{{2{\rm{ \mathsf{ π} }}{f_{r,{{\rm{N}}_2}}}}} \end{array} \right. $ (5)

式中:p0T0分别为环境大气的压强和温度;Trefps分别为参考温度和相应的参考压强,Tref=293.15 K,ps=101325 Pa; h为大气中水分子绝对湿度(单位:%),它和相对湿度hr的关系为:

$ h = 100\frac{{{p_w}}}{{{p_0}}} = {h_r}\frac{{{p_{{\rm{sat}}}}}}{{{p_0}}} = {h_r}\frac{{{p_{{\rm{sat}}}}}}{{{p_s}}}\frac{{{p_s}}}{{{p_0}}} $ (6)

其中:psat为水的饱和蒸汽压。其与参考压强ps的关系为:

$ \lg \frac{{{p_{{\rm{sat}}}}}}{{{p_s}}} = - 6.8346{\left( {\frac{{{T_s}}}{{{T_0}}}} \right)^{1.261}} + 4.6151 $ (7)

其中:Ts=273.16 K。

(2) 热黏吸收效应

式(1)热黏吸收效应中,需要确定扩散系数δ[38],其计算式为:

$ \delta {\rm{ = }}\frac{\mu }{{{\rho _0}}}\left( {\frac{4}{3} + \frac{{{\mu _B}}}{\mu } + \frac{{{{\left( {\gamma - 1} \right)}^2}\kappa }}{{\gamma R\mu }}} \right) $ (8)

式中:μμB分别为剪切黏性系数和体积黏性系数,其比值近似为0.6;κ为热传导系数;R为气体常数。其中,μκ是温度的函数,其表达式为:

$ \left\{ \begin{array}{l} \frac{\mu }{{{\mu _r}}} = {\left( {\frac{{{T_0}}}{{{T_r}}}} \right)^{3/2}}\frac{{{T_r} + {T_\mu }}}{{{T_0} + {T_\mu }}}\\ \frac{\kappa }{{{\kappa _r}}} = {\left( {\frac{{{T_0}}}{{{T_r}}}} \right)^{3/2}}\frac{{{T_r} + {T_A}{{\rm{e}}^{ - {T_B}/{T_r}}}}}{{{T_0} + {T_A}{{\rm{e}}^{ - {T_B}/{T_0}}}}} \end{array} \right. $ (9)

式中:Tr=300 K,相应地,μr=1.846×10-5kg/(m·s), κr=2.624×10-2 W/(m·K); Tμ=110.4 K; TA=245.4 K; TB=27.6 K。

1.2 声爆传播射线及声线管面积

声爆传播射线是广义Burgers方程的传播域,即方程沿声爆传播射线进行求解。在大气分层效应和风梯度存在的情况下,声爆的传播路径通常不是直线和规则圆弧。根据几何声学理论,有限小信号的扰动波在大气中的传播路径遵循Snell准则[40],即在有风情况下,声爆传播射线的微分方程为:

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}\mathit{\boldsymbol{R}}}}{{{\rm{d}}t}} = {c_0}\mathit{\boldsymbol{N}} + \mathit{\boldsymbol{W}}\\ \frac{{{\rm{d}}\mathit{\boldsymbol{N}}}}{{{\rm{d}}t}} = \left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{N}} \otimes {\mathit{\boldsymbol{N}}^{\rm{T}}}} \right)\nabla \left( {{c_0} + \mathit{\boldsymbol{W}} \cdot \mathit{\boldsymbol{N}}} \right) \end{array} \right. $ (10)

式中:R为声线路径矢量;N为波面单位法向量;W为当地风速;c0为环境大气声速;⊗为克罗内克积;I为单位矩阵。

对(10)式进行离散,可得:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{R}}\left( {i + 1} \right) = \mathit{\boldsymbol{R}}\left( i \right) + \Delta \mathit{\boldsymbol{R}}\left( i \right)\\ \mathit{\boldsymbol{N}}\left( {i + 1} \right) = \mathit{\boldsymbol{N}}\left( i \right) + \Delta \mathit{\boldsymbol{N}}\left( i \right) \end{array} \right. $ (11)

其中,

$ \left\{ \begin{array}{l} \Delta \mathit{\boldsymbol{R}}\left( i \right) = \left( {{c_0}\mathit{\boldsymbol{N}} + \mathit{\boldsymbol{W}}} \right)\left| {_i} \right.\Delta t\\ \Delta \mathit{\boldsymbol{N}}\left( i \right) = \left[ {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{N}} \otimes {\mathit{\boldsymbol{N}}^{\rm{T}}}} \right)\nabla \left( {{c_0}\mathit{\boldsymbol{W}} + \mathit{\boldsymbol{N}}} \right)} \right]\left| {_i} \right.\Delta t \end{array} \right. $ (12)

式中:Δt为声爆在空间传播过程中第i步到第i+1步的时间间隔。

在计算声线管面积时,通常的做法是计算由四条声线围成的“管”[41],如图 2所示。图中Δφ为声线在飞机周向上的角度增量。则声线管面积的计算式为:

$ A = \frac{1}{2}\left\{ {\left( {{\mathit{\boldsymbol{R}}_4} - {\mathit{\boldsymbol{R}}_1}} \right) \times \left( {{\mathit{\boldsymbol{R}}_3} - {\mathit{\boldsymbol{R}}_2}} \right)} \right\} \cdot {\mathit{\boldsymbol{N}}_1} $ (13)

图 2 声线管面积求解示意图 Fig.2 Sketch of ray tube area calculation
1.3 数值求解方法

用数值方法求解方程(1)时,直接对分子弛豫过程进行离散比较困难,通常的求解策略是采用“算子分裂法”[22]。在很小的空间推进步Δσ内,分别单独计算式(1)中各个效应对声学压强的影响,即依次单独求解式(14)-(18),并将前一方程的解作为后一方程的输入,从而达到各个效应解耦的目的。已经证明,当Δσ足够小时,采用算子分裂法的计算结果收敛于式(1)的解[42-43]

$ \frac{{\partial P}}{{\partial \sigma }} = \sum\limits_j {\frac{{{C_j}}}{{\left( {1 + {\theta _j}\frac{\partial }{{\partial \tau }}} \right)}}\frac{{{\partial ^2}P}}{{\partial {\tau ^2}}}} $ (14)
$ \frac{{\partial P}}{{\partial \sigma }} = \frac{1}{\mathit{\Gamma }}\frac{{{\partial ^2}P}}{{\partial {\tau ^2}}} $ (15)
$ \frac{{\partial P}}{{\partial \sigma }} = - \frac{1}{{2\mathit{A}}}\frac{{\partial A}}{{\partial \sigma }}P $ (16)
$ \frac{{\partial P}}{{\partial \sigma }} = \frac{1}{{2{\rho _0}{c_0}}}\frac{{\partial \left( {{\rho _0}{c_0}} \right)}}{{\partial \sigma }}P $ (17)
$ \frac{{\partial P}}{{\partial \sigma }} = P\frac{{\partial P}}{{\partial \tau }} $ (18)

采用有限差分方法对上述5式进行离散求解的过程,在Cleveland的博士论文[14]第四章有详细论述,本文不作过多介绍。但需要注意的是,在计算考虑“大气风”的几何声学效应时,Blokhintzev不变量的离散形式如下:

$ p_i^{k + 1}{\left( {\sqrt {\frac{{c_n^2A}}{{{\rho _0}c_0^3}}} } \right)_{k + 1}} = p_i^k{\left( {\sqrt {\frac{{c_n^2A}}{{{\rho _0}c_0^3}}} } \right)_k} $ (19)

式中:k表示空间推进过程中的第k步,i为第i个离散波形点。

根据算子分裂法和1.2节提到的射线追踪技术,本文开发了一套可考虑“大气风”效应的远场声爆高精度预测程序“bBoom”。图 3为该程序的流程。首先,将输入的近场声爆信号插值到均匀网格上;其次,运用1.2节的射线追踪技术求解声爆在大气中的传播路径,以便方程在空间方向上推进;然后,采用算子分裂法将式(1)右边各项效应依次进行单独求解,每一效应求解过程中均采用有限差分方法;最后,沿传播路径重复上一过程,直至传播到地面,即获得地面声爆波形。


图 3 所开发的bBoom程序求解广义Burgers方程的框架 Fig.3 Framework of "bBoom" code solving augmented Burgers equation using an operator splitting method

图 4为所开发的“bBoom”程序中,坐标系定义及声爆传播周向角ϕ的定义。坐标原点O定义为声爆由机体开始向外传播时的飞机位置,x轴为沿飞机轴线指向飞行方向,y轴为垂直于飞机轴线指向飞行员的右侧,z轴根据右手坐标系定义,其垂直于xOy平面向下。声爆传播周向角定义为由z轴向y偏转时为正,ϕ=0°时代表声爆向飞机正下方传播(Under-track)。另外,定义x轴正向与正北向的夹角β为飞机飞行方位角。针对有“大气风”的大气剖面,该角度对地面声爆计算具有重要作用。


图 4 “bBoom”程序所采用的坐标系统及传播周向角定义 Fig.4 Definition of coordinate system and roll angle of sonic-boom propagation in the code "bBoom"
2 算例验证与讨论 2.1 标模算例验证

以AIAA第二届声爆预测研讨会(SBPW-2)[44-46]提供的轴对称标模、低声爆概念机C25D构型和LM1021构型为例,验证所发展的方法对远场声爆预测结果的正确性。其中,轴对称标模是基于Cart3D流场求解器设计的旋成体,其正下方3倍模型长度处采用Euler方程计算的近场声爆信号与C25D构型相同。

2.1.1 轴对称标模算例

该算例主要以轴对称外形验证由不同近场提取位置传播到远场波形的计算结果。图 5为声爆预测研讨会提供的轴对称标模及近场压强提取位置示意图,模型长度为32.92 m,图 6为不同提取位置(H/L=1.0, H/L=3.0, H/L=5.0,L表示模型长度)的近场压强信号(由研讨会提供)。远场传播条件:Ma=1.6, α=0°,巡航高度为15.760 km,标准大气条件。图 7为标准大气的压强、密度、温度和相对湿度剖面。


图 5 轴对称标模近场压强提取位置示意图 Fig.5 Locations of near-field pressure extraction for an axi-symmetric body model


图 6 不同近场位置提取的压强信号对比 Fig.6 Comparison of near-field pressure signal extracted at different locations


图 7 标准大气剖面图 Fig.7 The standard atmosphere from U. S.

取地面反射因子为1.9,将计算的远场波形与NASA的sBOOM程序计算结果(基于广义Burgers方程的传播结果)、波形参数法及添加经验上升时间的结果进行了对比,如图 8~图 11所示。波形参数法的结果是基于Thomas公开的源码[13]计算的,并采用Plotkin的“3/p”经验方法[47]添加声爆上升时间,根据激波处扰动压强的双曲正切分布添加到原始波形参数法获得的波形上[48]。波形参数法及上升时间的添加方法已经集成到本团队开发的“FL-BOOM”[27]程序中。


图 8 由近场H/L=1.0处的压强信号传播至地面的波形对比 Fig.8 Comparison of ground waveforms propagated from the one body length below the axis-symmetric body model


图 9 由近场H/L=3.0处的压强信号传播至地面的波形对比 Fig.9 Comparison of ground waveforms propagated from the three body lengths below the axis-symmetric body model


图 10 由近场H/L=5.0处的压强信号传播至地面的波形对比 Fig.10 Comparison of ground waveforms propagated from the five body lengths below the axis-symmetric body model


图 11 基于广义Burgers方程的传播结果与波形参数法及添加上升时间的结果对比 Fig.11 Comparison of predicted sonic-boom waveforms using the approaches based on nonlinear Burgers equation, waveform parameter method and shock rise time added method

由于缺乏实验数据,我们将从不同近场压强提取位置计算的远场波形与NASA开发的sBOOM程序计算结果进行比较,结果表明两者几乎一致(如图 8图 9图 10所示)。与波形参数法结果相比(如图 11),基于广义Burgers方程的结果能够从本质上解决声爆上升时间计算问题。

2.1.2 低声爆概念机C25D构型标模算例

该算例用于验证不同周向角上远场声爆的计算结果。图 12为研讨会提供的C25D标模及近场压强提取位置示意图,模型长度仍为32.92 m,近场压强提取位置为距离飞机轴线3倍机身长度处,选取的两个周向角分别为0°和30°。图 13为相应两个周向角下的近场压强信号(由研讨会提供)。远场传播时,该模型以Ma =1.6、α =3.375°在高度15.760 km处水平飞行,传播大气为标准大气条件。


图 12 低声爆构型C25D及近场压强信号提取位置 Fig.12 Locations of near-field pressure extraction for the low-boom configuration C25D


图 13 不同周向角下C25D构型相应近场压强信号 Fig.13 Comparison of near-field pressure signals at the different azimuthal angles around the configuration C25D

取地面反射因子为1.9,同样将计算的远场波形与sBOOM程序计算的结果进行对比,如图 14图 15所示。由对比结果可知,在0°周向角(飞机正下方)与30°周向角时,计算结果与sBOOM计算结果符合很好,进一步表明所开发的程序针对复杂飞机外形的远场声爆计算仍具有较高可信度。


图 14 飞机正下方地面波形对比 Fig.14 Comparison of sonic-boom waveforms on the ground below the aircraft


图 15 周向角为30°时地面波形对比S Fig.15 Comparison of sonic-boom waveforms on the ground at a roll angle=30°
2.1.3 LM1021标模算例

该标模算例用于验证大气中存在“风”时,bBoom程序对远场声爆波形的计算结果。LM1021标模构型示意图及近场声爆信号提取位置如图 16所示。模型长度为71.12 m,巡航马赫数为1.6,巡航高度为16.764 km,在计算近场声爆信号或风洞试验时模型绕机头有向下2.1°的偏转(等同于来流迎角为2.1°)。在声爆传播周向角为0°、30°和-30°时,距离飞行路径3.1299倍机身长度位置处,提取的近场声爆信号如图 17所示。


图 16 LM1021构型及近场声爆提取位置和传播方向 Fig.16 LM1021 configuration, near-field pressure signals extraction location and the direction of sonic boom


图 17 LM1021构型3.1299倍机身长度处的近场声爆信号 Fig.17 Near-filed pressure signals at 3.1299 body lengths of LM1021

研讨会共提供了四种大气剖面用于计算该标模的远场声爆,其中两种存在“大气风”(大气季风)。本文选取有风大气剖面中的一种(大气剖面profile1)对bBoom程序的远场计算结果进行验证,图 18为大气剖面profile1示意图,详细数据读者可以参见[44]。需要说明的是,风剖面中“N_wind”和“E_wind”的数值为正时,分别表示风自南向北吹和自西向东吹。


图 18 SBPW-2提供的大气风剖面及相对湿度和压强剖面 Fig.18 Profiles of wind, relative humidity and atmospheric pressure at the atmosphere 1# provided by SBPW-2

远场声爆传播过程中,取地面反射因子为1.9,当标模飞行方向为正东方向(即β=90°)时,本文计算的远场声爆结果与sBOOM结果对比如图 19所示。由图可知,在周向角ϕ=0°和30°时,两者的计算结果几乎一致。另外,尽管在周向角ϕ=-30°时,本文计算结果与sBOOM结果存在差别,但与日本JAXA的Xnoise结果几乎一致。对比结果一定程度上能够说明本文发展的bBoom程序能够有效模拟“大气风”存在时的声爆传播过程。


图 19 在SBPW-2提供的大气剖面profile1下,远场声爆波形计算结果对比 Fig.19 Comparison of sonic-boom waveforms on the ground at the atmosphere 1# profile by SBPW-2
2.2 近场提取位置对远场声爆计算的影响

近场压强信号的提取位置关系到近场计算域的选取。当计算域选取较小时,空间中强激波膨胀波系来不及合并,压强扰动较强,如果将其作为广义Burgers方程的信号输入,那么由于小扰动假设不成立,就可能会得到错误的远场解。相反地,当计算域选择过大时,为保证能够很好地捕捉空间中激波膨胀波系,势必会增加网格量,无形之中就增大了近场的计算量,造成计算资源的浪费,影响声爆的高精度评估效率。

本节以C25D标模为例,将周向角ϕ=0°方向上1L、3L、5L(L为机身长度)位置处的近场声爆传播到远场,研究近场提取位置对远场声爆计算的影响。其中,近场声爆信号由NASA研究人员Aftosmis采用Cart3D求解器计算中等网格获得[49]。在标准大气下,取地面反射因子为1.9,则远场对比如图 20所示。对比可知,由近场H/L=3.0和H/L=5.0得到的远场波形除在后体波形存在微小差别外,其他位置都符合很好,而由近场H/L=1.0得到的远场波形与前两者相比存在很大差异。因此,针对该标模,在用CFD计算近场信号时,3倍左右机身长度处提取的近场信号作为传播方程的输入时具有较高可信度。


图 20 不同近场位置传播到远场的波形对比 Fig.20 Comparison of far-field waveforms propagated from different near-field locations
2.3 “大气风”对远场声爆计算的影响

本节以C25D标模为例,研究“大气风”风向、风速对远场声爆计算的影响。在下面两种讨论情况中,对C25D标模近场声爆进行传播的条件都为:巡航高度15.760 km,巡航马赫数1.6,巡航迎角3.375°,近场提取位置为3倍机身长度处,声爆传播周向角ϕ=0°,地面反射因子为1.9,大气条件为在标准大气下添加相应的风剖面。

2.3.1 风向对远场声爆计算的影响

为了研究“大气风”风向对远场声爆计算的影响,给定任意单向风剖面(风向沿海拔不变),如表 2所示,记为“wind_1”。表中描述的风向为自西向东,风速大小参照SBPW-2提供的1#大气剖面中的风速(如图 18左图)给定。

表 2 给定任意单向风剖面“wind_1” Table 2 Arbitrarily wind profile "wind_1"

为了不失风向的一般性,将飞机飞行方向的方位角分为β={0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°},即飞机分别向八个方向飞行(可参见图 4)。在该风剖面情况下,地面声爆波形对比如图 21所示。由于β=0°和180°时,风对周向角为ϕ=0°的声爆传播影响效果相同,因此这两种声爆波形相同。同样地,β=45°和135°、β=225°和315°时情况也相同。由对比图可知,在八种飞行方向下,地面波形形状相似,但飞机向正东方向飞行时,地面声爆超压值最大,而向正西方向飞行时,地面声爆超压值最小。


图 21 在风剖面“wind_1”下,飞机向不同方向飞行时,地面声爆波形对比 Fig.21 Comparison of far-field waveforms on different flight directions in the standard atmosphere with wind profile "wind_1"
2.3.2 风速对远场声爆计算的影响

为研究整个剖面风速对地面声爆计算的影响,在2.3.1节的风剖面“wind_1”基础上,将剖面风速增大1倍,记为“wind_2”。当飞机向正东和正西方向飞行时,在“wind_1”和“wind_2”两种风剖面下,地面声爆波形对比如图 22所示。由图可知,针对该单向风剖面,当风速增大时,飞机向正东方向飞行引起的地面声爆超压值增加;而飞机向正西飞行时,地面声爆超压值减小。


图 22 在wind_1和wind_2风速下计算地面声爆波形对比 Fig.22 Comparison of sonic-boom waveforms on the ground in the standard atmosphere with "wind_1" and "wind_2"

为了研究不同高度的风速对地面声爆计算的影响,在“wind_1”风剖面基础上构造两个风剖面。将“wind_1”风剖面5 km高度处的风速增加到40 m/s,0 km至5 km高度范围内风速线性变化,其他高度风速不变,记为“wind_3”。将“wind_1”风剖面10 km至20 km高度范围内的风速增加到40 m/s,其他高度风速不变,记为“wind_4”。

当飞机向正东和正西方向飞行时,在无风、“wind_1”、“wind_3”和“wind_4”四种情况下,地面声爆波形对比如图 23所示。由图可知,对于飞机向正东飞行的情形,声爆超压值由大到小的顺序为“wind_4”、“wind_1”、“wind_3”和无风大气。而对于飞机向正西飞行的情形,声爆超压值的顺序正好相反。由当前的风剖面对比可知,较低海拔高度的速度增加可能会削弱风对地面波形影响。


图 23 在无风、有风(wind_1、wind_3和wind_4)情况下计算的地面声爆波形对比 Fig.23 Comparison of sonic-boom waveforms on the ground in the standard atmosphere with "no wind", "wind_1", "wind_3", and "wind_4

由于“大气风”的影响效果是一个积分效应,从目前的结果来看,剖面风速对地面声爆波形的影响较为复杂。为了给出“大气风”对声爆超压值影响的定性结论,本文给出上述四种风剖面情况下,地面声爆超压值、传播时间和传播路径长度对比,如附表 1所示。前文提到,对于该风剖面,β=0°和180°、β=45°和135°、β=225°和315°时波形相同,因此,表中只给出了β=0°、45°、90°、225°、270°时的传播情况。根据表中数据,给出了在这4中风剖面情况下,声爆在大气中的传播时间与地面声爆超压值的关系,如图 24所示。除了“wind_3”风剖面的情况,传播时间越长,声爆超压值越小。

附表 1 地面声爆超压值、传播时间和传播路径长度对比 Appendix table 1 Comparison of sonic-boom overpressure, propagation time and path length


图 24 声爆超压值与传播时间关系 Fig.24 Relationship between far-field overpressure and propagation time
2.4 “大气风”对声爆地面影响域计算的影响

选取与2.3节中给定的“wind_1”和“wind_2”风剖面,以C25D标模向北飞行(β=0°)为例,考察“大气风”对声爆影响域的影响。C25D标模的计算状态为:巡航高度15.760 km,巡航马赫数1.6,巡航迎角3.375°,近场提取位置为3倍机身长度处。在三种大气环境下,声爆由飞机巡航高度传播到地面的射线路径及地面影响域的形成如图 25所示,地面声爆影响域对比如图 26所示。图 26中未标角度的点其周向角两两相差10°。


图 25 声爆传播射线及声爆地面影响域的形成示意图 Fig.25 Sketch propagation path of sonic boom and the sonic boom carpet formation


图 26 声爆地面影响域对比 Fig.26 Comparison of sonic-boom carpets

由图可知,在无风标准大气下,声爆由飞机向地面传播的射线关于飞行轨迹所在的铅垂面对称,其所形成的地面影响域关于飞行轨迹在地面上的投影线对称,射线与地面相交的最大周向角为50.34°,最小周向角为-50.34°。然而,在西风作用的情况下,最大周向角、地面声爆影响域的位置和范围发生变化,地面声爆影响域向东方向扩张。这是水平分层的“大气风”对声爆传播射线综合作用的结果。在水平分层大气环境下,声爆传播射线满足Snell准则,其会有向传播速度较低区域折射的趋势。

在此类剖面的侧风作用下,当传播射线的周向角大于某一角度时,相同周向角下传播射线与地面的交点向逆风方向移动,且射线与地面相交的最大周向角发生变化。针对该算例,在飞机飞行方向的迎风侧,当周向角小于-40°时,尽管侧风作用会使射线与地面相交位置向偏离飞行轨迹的方向移动,有增大该侧声爆影响范围的趋势,但是由于侧风作用下最大周向角减小,综合效应使该侧声爆影响域范围减小。而在飞行方向的顺风侧,情况则恰好相反,综合作用结果使该侧声爆影响域范围增大。由于顺风侧声爆影响范围的增加量大于迎风侧的减小量,造成该算例中计算的声爆影响范围向东扩张。

3 结论

本文基于广义Burgers方程的数值求解方法,发展了一套适用于预测超声速客机远场声爆的高精度方法,并自主开发了计算程序“bBoom”。采用AIAA第二届声爆预测研讨会提供的三个标模算例验证了本文开发程序的正确性。在此基础上,开展了近场压强提取位置对远场计算结果的影响研究,以及“大气风”对远场波形与声爆影响域计算结果的影响研究。结果表明:

1) 在标准无风大气和有风大气环境下,“bBoom”程序计算的远场声爆结果具有较高可信度。

2) 对于类C25D标模构型,为了确保远场声爆预测结果具有较高精度,近场压强信号提取位置应在机身下方约3倍机身长度处。

3) “大气风”对声爆传播的影响较为复杂,在文中给定的风剖面情况下,当Ma=1.6且顺风飞行时,由周向角0°传播到地面的声爆超压值增大,而逆风飞行时减小。

4) 针对文中给定单向风剖面,飞机在垂直于风向飞行的情形,当周向角大于某一角度时,相同周向角的声爆传播射线与地面的交点相比于无风情况向逆风侧移动,且最大周向角发生变化,它们的综合效应会改变地面声爆影响域。

由于近场声爆提取位置还受马赫数等因素的影响,同时飞行高度和风剖面也会对声爆产生影响,因此下一步拟开展不同马赫数下,近场提取位置对远场声爆计算结果的影响研究,以及不同飞行高度和不同风剖面对声爆的影响研究。

参考文献
[1]
PLOTKIN J K, MAGLIERI J D. Sonic boom research: history and future[R]. AIAA-2003-3575, 2003.
[2]
SAKATA K. Supersonic experimental airplane(NEXST) for next generation SST technology[R]. AIAA-2002-0527, 2002.
[3]
MORGENSTERN M J, STELMACK M, JHA D P. Advanced concept studies for supersonic commercial transports entering service in 2030-35 (N+3)[R]. AIAA-2010-5114, 2010.
[4]
WELGE H R, NELSON C, BONET J. Supersonic vehicle systems for the 2020 to 2035 timeframe[R]. AIAA-2010-4930, 2010.
[5]
WHITHAM G B. The flow pattern of a supersonic projectile[J]. Communications on Pure & Applied Mathematics, 1952, 5(3): 301-348.
[6]
WALKDEN F. The shock pattern of a wing-body combination, far from the flight path[J]. Aeronautical Quarterly, 1958, 9(2): 164-194. DOI:10.1017/S0001925900001372
[7]
HAYES D W, HAEFELI C R, KULSRCD E H. Sonic boom propagation in a stratified atmosphere, with computer program[R]. NASA CR-1969-1299, 1969.
[8]
CARLSON HW. Simplified sonic-boom prediction[R]. NASA TP-1978-1122, 1978.
[9]
PLOTKINJ K. A rapid method for the computation of sonic booms[R]. AIAA-1993-4433, 1993.
[10]
YAMASHITA R, SUZUKI K. Full-field sonic boom simulation in real atmosphere[R]. AIAA-2014-2269, 2014.
[11]
YAMASHITA R, SUZUKI K. Rise time prediction of sonic boom by full-field simulation with thermal nonequilibrium[R]. AIAA-2015-2583, 2015.
[12]
YAMASHITA R, SUZUKI K. Full-field sonic boom simulation in stratified atmosphere[J]. AIAA Journal, 2016, 54(10): 3223-3231. DOI:10.2514/1.J054581
[13]
THOMAS C L. Extrapolation of sonic boom pressure signatures by the waveform parameter method[R]. NASA TN D-1972-6832, 1972.
[14]
CLEVELAND R O. Propagation of sonic booms through a real, stratified atmosphere[D]. Austin: The University of Texas at Austin, 1995.
[15]
RALLABHANDI S K. Advanced sonic boom prediction using the augmented Burgers equation[J]. Journal of Aircraft, 2011, 48(4): 1245-1253. DOI:10.2514/1.C031248
[16]
TAKENO J, MISAKA T, SHIMOYAMA K, et al. Analysis of sonic boom propagation based on the KZK equation[R]. AIAA-2015-0745, 2015.
[17]
PIERCE A D. Acoustics: an introductionto its physical principles and applications[M].Acoustical Soc. of America, 1981.
[18]
KANG J. Nonlinear acoustic propagation of shock waves through the atmosphere with molecular relaxation[D]. PA: The Pennsylvania State University, 1991.
[19]
ROBINSON L D.Sonic boom propagation through an inhomogeneous, windy atmosphere[D]. Austin: University of Texas at Austin, 1991.
[20]
ROBINSON L D. Anumerical model for sonic boom propagation through an inhomogeneous, windy atmoshpere[R]. NASA CP-1992-1372, 1992.
[21]
CLEVELAND R O. Effect of stratification and geometrical spreading on sonic boom rise time[R]. NASA N95-14880, 1995.
[22]
CLEVELAND R O. Propagation of sonic booms through a real, stratified atmosphere[D]. Austin: The University of Texas at Austin, 1995.
[23]
RALLABHANDI S K. Advanced sonic boom prediction using the augmented Burgers equation[J]. Journal of Aircraft, 2011, 48(4): 1245-1253. DOI:10.2514/1.C031248
[24]
YAMAMOTO M, et al. Long-range sonic boom prediction considering atmospheric effects[C]//40th International Congress and Exposition on Noise Control Engineering 2011, INTER-NOISE, 2011.
[25]
YAMAMOTO M, et al. Numerical simulation for sonic boom propagation through an inhomogeneous atmosphere with winds[J]. AIP Conference Proceedings, 2012, 1474(1): 339-342.
[26]
冯晓强.声爆计算方法研究及在超声速客机设计的应用[D].硕士学位论文.西安: 西北工业大学, 2012: 9-15.
FENG X Q.The research of sonic boom prediction method and application in supersonic aircraft design[D]. Master Thesis, Xi'an: Northwestern Polytechnical University, 2012. (in Chinese)
[27]
冯晓强.超声速客机低声爆机理及设计方法研究[D].博士学位论文.西安: 西北工业大学, 2014: 15-47.
FENG X Q. The research of low sonic boom mechanism and design method of supersonic aircraft[D]. Xi'an: Northwestern Polytechnical University, 2014. (in Chinese)
[28]
冯晓强, 李占科, 宋笔锋. 超音速客机音爆问题初步研究[J]. 飞行力学, 2010, 28(6): 21-23.
FENG X Q, LI Z K, SONG B F. Preliminary analysis on the sonic boom of supersonic aircraft[J]. Flight Dynamics, 2010, 28(6): 21-23. (in Chinese)
[29]
冯晓强, 宋笔锋, 李占科. 低声爆静音锥设计方法研究[J]. 航空学报, 2013, 34(5): 1009-1017.
FENG X Q, SONG B F, LI Z K. Research of low sonic boom quiet spike design method[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(5): 1009-1017. (in Chinese)
[30]
FENG X Q, LI Z K, SONG B F. Research of low boom and low drag supersonic aircraft design[J]. Chinese Journal of Aeronautics, 2014, 27(3): 531-541. DOI:10.1016/j.cja.2014.04.004
[31]
徐悦, 宋万强. 典型低音爆构型的近场音爆计算研究[J]. 航空科学技术, 2016, 27(7): 12-16.
XU Y, SONG W Q. Near field sonic boom calculation on typical LSB configurations[J]. Aeronautical Science & Technology, 2016, 27(7): 12-16. (in Chinese)
[32]
徐悦.超音速飞行器地面音爆强度预测方法对比研究[C]//全国气动声学学术会议, 2016.
[33]
王刚, 马博平, 雷知锦, 等. 典型标模音爆的数值预测与分析[J]. 航空学报, 2018, 39(1): 164-176.
WANG G, MA B P, LEI Z J, et al. Simulation and analysis for sonic boom on several benchmark cases[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(1): 164-176. (in Chinese)
[34]
马博平, 王刚, 雷知锦, 等. 网格对声爆近场声爆预测的数值研究[J]. 西北工业大学学报, 2018, 36(5): 865-874.
MA B P, WANG G, LEI Z J, et al. Numerical investigation of influence of mesh property in nearfield sonic boom prediction[J]. Journal of Northwestern Polytechnical University, 2018, 36(5): 865-874. DOI:10.3969/j.issn.1000-2758.2018.05.008 (in Chinese)
[35]
霍满, 谭廉华, 林大楷, 等.声爆标模近场计算对远场预测影响研究[C]//首届空气动力学大会, 绵阳: 中国空气动力学会, CARS-04-2018-122, 2018.
[36]
钱战森, 刘忠臣, 冷岩, 等. OS-X0试验飞行器声爆特性测量与数值模拟分析[C]//首届空气动力学大会, 绵阳: 中国空气动力学会, CARS-04-2018-116, 2018.
[37]
张绎典, 黄江涛, 高正红. 基于增广Burgers方程的音爆远场计算及应用[J]. 航空学报, 2018, 39(7): 122039.
ZHANG Y D, HUANG J T, GAO Z H. Far field simulation and applications of sonic boom based on augmented Burgers equation[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(7): 122039. (in Chinese)
[38]
PIERCE D A. Acoustics-an introduction to its physical principles and applications[M]. New York: Acoustical Society of America, 1989: 587-588.
[39]
BASS E H, SUTHERLAND C L, ZUCKERWAR J A, et al. Atmospheric absorption of sound:further developments[J]. Journal of the Acoustical Society of America, 1995, 97(1): 680-683. DOI:10.1121/1.412989
[40]
ONYEONWU R O. The effects of wind and temperature gradients on sonic boom corridors[R]. Univ. of Toronto Inst. for Aerospace Studies TN-168, 1971.
[41]
YAMAMOTO M, SAKIA T. A unified approach to an augmented Burgers equationfor propagation of sonic booms[J]. Acoustical Society of America, 2015, 137(4): 1857-1866. DOI:10.1121/1.4916833
[42]
LEE Y-S. Numerical solution of KZK equation for pulsed finite amplitude sound beams in thermoviscous fluids[D]. Austin: University of Texas at Austin, 1993.
[43]
LEE Y-S, HAMILTON F M. Time-domain modeling of pulsed finite-amplitude sound beams[J]. Journal of the Acoustical Society of America, 1995, 97(2): 906-917. DOI:10.1121/1.412135
[44]
PARK M. 2nd AIAA sonic boom prediction workshop[EB/OL]. (2017-1-17)[2018-7-20] https://lbpw.larc.nasa.gov/
[45]
PARK M A, NEMEC M. Nearfield summary and statistical analysis of the second AIAA sonic boom prediction workshop[R].AIAA-2017-3256, 2017.
[46]
RALLABHANDI S K. Summary of propagation cases of the second AIAA sonic boom prediction workshop[R].AIAA-2017-3257, 2017.
[47]
PLOTKIN J K. State of the art of sonic boom modeling[J]. Journal of the Acoustical Society of America, 2002, 111(1): 530-536. DOI:10.1121/1.1379075
[48]
SALAMONE J. Sonic boom simulation using conventional audio equipment[R].Baltimore, Maryland, NOISE-CON-2004, 2004.
[49]
ANDERSON G R, AFTOSMIS M J, NEMEC M. Cart3D simulations for the second AIAA sonic boom prediction workshop[J]. Journal of Aircraft, 2019, 56(3): 896-911. DOI:10.2514/1.C034842