CN108170878B - 一种超音速飞行器音爆预测方法 - Google Patents

Abstract

本发明涉及一种超音速飞行器音爆预测方法，采用CFD数值模拟与远场外插相结合的方法实现超音速飞行器的音爆预测，包括绘制满足要求的结构或非结构网格，确定CFD求解器以及求解方法；进行超音速飞行器的CFD求解；提取飞行器近场静态压力特征；通过任意数量的线性线段来近似表示初始波形的形状；周期性求解描述波形参数变化的一阶耦合微分方程组，得到下一时刻波形的信息；重复执行直至得到指定高度处的过压值的步骤。本方法系统地模拟了飞行器外形、激波、膨胀波等因素的影响，计算精度大大高于传统线性方法，同时更加逼真地模拟了压力波在大气中的传播情况。

Description

一种超音速飞行器音爆预测方法

技术领域

本发明涉及到计算流体力学领域，具体而言，涉及一种超音速飞行器音爆预测方法。

背景技术

随着空气动力学、电子系统以及材料和制造工艺等各个学科技术的进步和不断积累，研制超音速飞机的条件和时机已逐渐成熟。可以预见，在未来民用航空市场上超音速大型客机必将逐渐成为主角，这是人类航空朝着更快、更经济环保、更安全舒适的一个必然方向之一。

音爆是一种非线性气动声学现象。当飞机做超音速飞行时，其机鼻、机翼、尾翼等均会产生激波；另一方面，由于机身外形变化，激波之间会存在膨胀波系。这两种波系相互作用在增加了飞机近场流场复杂性的同时，产生的压力扰动随飞机一起运动即成为音爆的声源。由于该声源发出的声波具有巨大幅值，在大气层里经过一定距离衰减之后仍会对飞行通过的区域带来严重的噪声污染，引起人焦躁不安的情绪，甚至造成建筑物的损坏。美国联邦航空局(FAA)制定的《联邦航空条例》(FAR)91部817条款，不允许商用飞机或私人飞机在美国的陆地上空进行超音速飞行，其他国家和地区的民航管理机构也制定了类似的法规，目的是避免音爆对居民的影响。

大多数音爆预测和优化方法都基于Whitham的修正线化理论，此理论基于弱激波理论，是对线性理论的修正并且考虑光滑旋成体的扰动合并成激波。为了满足弱激波理论要求，假设流体是定常、无粘、可压缩、轴对称、无旋和等熵的。一些试验和分析对马赫数接近3.0时线性理论的有效性提出了质疑。在高马赫数区域，强激波导致明显的高阶熵增现象，线性方法完全忽略了这些效应，并且升力体附近区域存在明显的横流效应。

理想的计算方法是整个传播过程都采用CFD方法。但是这么做有两个阻碍：一是所需网格量太大，需要从巡航高度(一般60000英尺)到地面生成足够密的网格，计算效率无法保证。随着计算机的计算能力逐年提高，但仍难满足需要；二是CFD本身的功能问题。一般而言，CFD不是用来进行远场分析的工具，其主要目的是计算飞行器近场的空气动力学特性。远距离计算时，由于精度和分辨率的损失，会造成典型流场参数和激波较大的计算误差。

从目前资料来看，音爆特征研究已成为新一代超音速飞机发展的关键技术，国外在这一领域已做过较为深入地研究，而目前国内对超音速飞机音爆特征的研究还处于初步探索阶段，且音爆的风洞试验技术未曾开展过。

发明内容

本发明的目的在于克服现有技术的不足，提出一种CFD数值模拟与远场外插相结合的超音速飞行器音爆预测方法。这种方法可以用计算机高级程序语言实现，并通过计算机相应程序软件运行程序来预测超声速飞行器音爆特性。采用CFD耦合远场外插方法，先利用CFD数值模拟超音速飞行器近场附近流场特征，并利用CFD近场压力特征作为远场外插的初值，通过求解描述波形参数变化的一阶耦合微分方程组得到压力波信息，最终得到指定高度处音爆特征。本发明中的预测方法在保证预测结果准确度的同时，又兼顾计算效率，最大程度上满足了超音速飞行器音爆预测的要求。

考虑到现有技术的上述问题，根据本发明公开的一个方面，本发明采用以下技术方案：

一种超音速飞行器音爆预测方法，采用CFD数值模拟与远场外插相结合的方法实现超音速飞行器的音爆预测，具体包括以下步骤：

1)CFD数值模拟前期准备工作：根据超音速飞行器的构型、来流条件、关注区域信息，绘制满足要求的结构或非结构网格，确定CFD求解器以及求解方法；

2)根据步骤1)所绘制的结构或非结构网格以及所确定的CFD求解器，进行超音速飞行器的CFD求解，得到飞行器的近场压力特征，如果飞行器是加速或者机动状态，则需根据时间依赖性来对近场压力特征进行修改；

3)提取飞行器近场静态压力特征：从步骤2)求解得到的飞行器的近场压力特征中提取飞行器下方1-3倍特征长度处压力分布作为近场静态压力特征的初始波形，近场静态压力特征提取位置要距离飞行器足够近，以保证CFD求解器本身存在的数值耗散不会影响结果的逼真度；与此同时提取位置还要距离飞行器足够远以保证可以忽略横流和升力效应的影响；

4)考虑斜率、压力增长和持续时间参数的影响，通过任意数量的线性线段来近似表示初始波形的形状，得到描述波形参数变化的一阶耦合微分方程组；

5)周期性求解描述波形参数变化的一阶耦合微分方程组，得到下一时刻波形的信息，所述下一时刻波形的信息包括传播方向，传播时间、高度和下一时刻波形的斜率、压力增长和持续时间；

6)重复执行步骤5)，直至计算高度降为指定高度，即得到指定高度处的过压值，作为音爆特征，评估音爆特征与地面的相互影响。

其特征在于，所述步骤5)具体为：

描述波形参数变化的一阶耦合微分方程组为：

其中m_i是线段i的斜率

Δp_i是第i段和第i-1段连接处穿过激波的压力增长，无激波时为0；λ_i是第i段的持续时间ΔT；

C₁,C₂为过程变量，其表达式为：

其中γ＝1.4为常数，a₀为所在高度音速；p₀为所在高度大气压；ρ₀为所在高度密度；c_n＝c+v·n为波形传播速度；v为风速，n是波阵面法向量；A为射线管面积。a₀,p₀,ρ₀,A,c_n是高度z的函数，所以沿着射线管变化；

不均匀大气有风条件下，对于任意波阵面形状传播的波来说，C₁和C₂沿着射线管是变化的，但是，如果假设这些量在很小的时间增长内是常数，则上述波形变形方程可积分得到解：

当：

或者

时

或

时

近似得到：

通过上述关系式可以看出，当给定

后，主要就是计算出C₁,C₂，由C₁,C₂的表达式可以看出，计算的关键是得到所在高度a₀,p₀,ρ₀,c_n和A随时间变化的微小变化率，因此首先需要采用合适的计算射线路径的方法得到射线路径，然后得到沿射线路经a₀,ρ₀,c_n的时间变化率

射线管面积可以通过射线管上多个点的坐标和方向计算得到；

针对射线管初始值和初始位置、射线方向，射线管面积以及空间推进方法描述射线管推进过程：上一高度的射线位置和方向确定之后，就可以在空间推进到下一高度，在压力传播过程中，射线方向取决于所在高度音速和风速，当高度发生变化时，大气属性发生变化的同时，音速和风速随之变化，所以传播方向随着高度变化，下一射线点的传播方向为：射线上某点从上一高度传播到下一高度的过程中存在空间偏移量，下一高度射线的传播方向为上一高度的传播方向加上传播过程中的偏移量，将相对偏移量投影到空间三个方向上，并将所得分量加到原来方向向量上就得到新的射线点的传播方向。

有益效果：

1.传统音爆预测方法中，采用经典的线性超音速空气动力学得到初始飞机压力扰动，但是在整个传播过程中都假设流场是线性的，这就使计算精度下降，导致激波和膨胀波的强度计算有误。本方法与传统音爆预测方法相比，使用CFD方法，通过求解Euler方程得到超音速飞行器近场流动特征。本方法系统地模拟了飞行器外形、激波、膨胀波等因素的影响，计算精度大大高于传统线性方法。

2.传统音爆预测方法中，采用修正线化理论来计算音爆传播过程，本方法与传统音爆预测方法相比，考虑超音速飞行器加速度，温度，压力和风速梯度对波形振幅和非线性畸变的影响，更加逼真地模拟了压力波在大气中的传播情况。

附图说明

图1示出了根据本发明的超音速飞行器音爆预测方法的使用流程图。

图2示出了根据本发明的超音速飞行器音爆预测方法的任意线性线段及参数示意图。

图3示出了本发明中具体实施方式中射线管推进示意图。

图4示出了本发明实施例中CFD数值模拟结果示意图。

图5示出了本发明实施例中地面音爆特征数值结果示意图。

具体实施方式

下面结合实施例对本发明作进一步地详细说明，但本发明的实施方式不限于此。

如图1所示为本发明的流程图，从图1可知，本发明提供的一种CFD数值模拟与远场外插相结合的超音速飞行器音爆预测方法。本方法需要将计算区域分为两部分，分别为近场区域和远场区域，近场区域需要采用CFD求解Euler方程，CFD求解Euler方程的方法为常规方法，而远场区域采用外插方法，为本发明的重点，远场外插方法具体描述如下：

如图2所示为本发明采用的远场外插方法的任意线性线段及参数示意图。在本方法中，通过任意数量的线性线段来近似压力波形状。对于任意线段，主要考虑的参数包括m_i,Δp_i和λ_i，其中m_i是线段i的斜率

Δp_i是第i段和第i-1段连接处穿过激波的压力增长，无激波时为0；λ_i是第i段的持续时间ΔT。

关于参数m_i,Δp_i,λ_i的描述压力波变化的一阶耦合微分方程组：

C₁,C₂为过程变量，其表达式为：

其中γ＝1.4为常数，a₀为所在高度音速；p₀为所在高度大气压；ρ₀为所在高度密度；c_n＝c+v·n为波形传播速度；v为风速，n是波阵面法向量；A为射线管面积。a₀,p₀,ρ₀,A,c_n是高度z的函数，所以沿着射线管变化；

不均匀大气有风条件下，对于任意波阵面形状传播的波来说，C₁和C₂沿着射线管是变化的，但是，如果假设这些量在很小的时间增长内是常数，则上述波形变形方程可积分得到解：

当：

或者

时

或

时

近似得到：

通过上述关系式可以看出，当给定

后，主要就是计算出C₁,C₂。由C₁,C₂的表达式可以看出，计算的关键是得到所在高度a₀,p₀,ρ₀,c_n和A随时间变化的微小变化率，因此首先需要采用合适的计算射线路径的方法得到射线路径，然后得到沿射线路经a₀,ρ₀,c_n的时间变化率

射线管面积可以通过射线管上四个点的坐标和方向计算得到；

如图3所示为本发明采用的远场外插方法的射线管推进示意图，进一步，针对射线管初始值和初始位置、射线方向，射线管面积以及空间推进方法进行描述。

射线管由4根射线组成，每根射线的方向不同。其中第一根射线的方向为实际波阵面的传播方向，其余三根定义方向是为了计算传播距离进而计算面积。

声学射线由超音速飞行器发出，传播方向垂直于波阵面。它们代表声学扰动在大气中的传播路径。射线的初始方向由飞机附近垂直于马赫锥的射线决定。为了计算射线管面积，选择四根射线，通过时间增长和方位增长区分。这些射线的初始方向由飞机飞行参数和方位角决定。

A.定义射线起点和初始方向

下述公式中R(m,n)代表射线上点的坐标；N(m,n)代表射线方向；m＝1～3代表三个方向；n＝1～4代表第n根射线；H代表飞行高度；vx是来流速度的水平分量；vy是来流速度的垂直分量；v0(1)是水平风速；v0(2)是垂直风速；α是攻角；β是侧滑角；μ是马赫角；

是方位角。

第一根射线位置(飞行高度处)：

R(1,1)＝0.0

R(2,1)＝0.0

R(3,1)＝H

第二根位置(与第一根起始点重合)：

R(1,1)＝R(1,2)

R(2,1)＝R(2,2)

R(3,1)＝R(3,2)

第三根位置(考虑小时间内飞行速度和风速作用)：

R(1,3)＝R(1,1)+(vx+v0(1))·dt

R(2,3)＝R(2,1)+(vy+v0(2))·dt

R(3,3)＝R(3,1)+mach·a0·sinβ·dt

第四根位置(与第三根起始点重合)：

R(1,4)＝R(1,3)

R(2,4)＝R(2,3)

R(3,4)＝R(3,3)

射线方向：

第一根射线就采用上式计算；第二根射线将方位角偏移小量

代入计算；第三根考虑加速之后，马赫数，攻角，方位角的变化，将相应值变化之后代入；第四根射线在第三根的基础上同样偏移方位角

代入。

上述公式可以得到初始射线的方向，而起始位置中需要将1、2分离，3、4分离。采用的方法即是在射线方向上移动一小段距离，计算公式为：

R(1,n)＝R(1,n)+(a0·N(1,n)+v0(1))·dt

R(2,n)＝R(2,n)+(a0·N(2,n)+v0(2))·dt

R(3,n)＝R(3,n)+a0·N(3,n)·dt

B.空间推进

上一高度的射线位置和方向确定之后，就可以在空间推进到下一高度。下一高度的射线点位置可以由上式计算得到(注：每根射线的音速，风速，方向均不同)。

在压力传播过程中，射线方向取决于所在高度音速和风速，当高度发生变化时，大气属性发生变化的同时，音速和风速随之变化，所以传播方向随着高度变化。下一射线点的传播方向的计算思想为：射线上某点从上一高度传播到下一高度的过程中存在空间偏移量，下一高度射线的传播方向为上一高度的传播方向加上传播过程中的偏移量，即：

偏移量的计算过程如下：假设上一高度为h₁，音速和风速分别为a₁,v₁；下一高度为h₂，音速和风速分别为a₂,v₂，则高度差为dz＝h₂-h₁，此方向上音速和风速的梯度为

波阵面法向速度梯度为

整个过程经历的时间为t，则相对于dz的偏移量为

将相对偏移量投影到空间三个方向上，并将所得分量加到原来方向向量上就得到新的方向传播方向。

实施例一

下面以一个具体应用实例进一步说明本发明提出的一种超音速飞行器音爆预测方法：

具体实施方案中，采用的超音速飞行器为F-5E飞机，其相关参数如表1所示：

表1

参数	数值
		巡航马赫数	1.4
巡航高度	32000ft
		攻角	0°

应用本发明中的方法进行超音速飞行器音爆特性预测主要由以下步骤：

1.根据超音速飞行器的构型以及表1所示来流条件，进行超音速飞行器的CFD求解，主要求解方程为Euler方程。

2.给定飞行器近场静态压力特征。从步骤1求解得到的超音速空间流场中提取飞行器下方2倍特征长度处压力分布作为初始波形,压力提取点方位角为0°。图4为提取的近场压力特征。

3.通过任意数量的线性线段来近似压力波形状，求解描述波形参数变化的一阶耦合微分方程组直至计算高度降为指定高度，即得到指定高度处音爆特征。图5为采用本发明所述方法预测得到的地面音爆特征。

本发明方法已经通过Fortran90/95计算机高级程序语言实现，并可通过CompaqVisual Fortran编译执行。但本发明不限定实现所用程序语言以及运行软件。上述实例是为了阐述本发明，不对本发明的保护范围构成限制。凡与本发明设计思路及工作原理相同的实施方案均在本发明的保护范围内。

Claims (1)

1.一种超音速飞行器音爆预测方法，采用CFD数值模拟与远场外插相结合的方法实现超音速飞行器的音爆预测，具体包括以下步骤：

1)CFD数值模拟前期准备工作：根据超音速飞行器的构型、来流条件、关注区域信息，绘制满足要求的结构或非结构网格，确定CFD求解器以及求解方法；

2)根据步骤1)所绘制的结构或非结构网格以及所确定的CFD求解器，进行超音速飞行器的CFD求解，得到飞行器的近场压力特征，如果飞行器是加速或者机动状态，则需根据时间依赖性来对近场压力特征进行修改；

3)提取飞行器近场静态压力特征：从步骤2)求解得到的飞行器的近场压力特征中提取飞行器下方1-3倍特征长度处压力分布作为近场静态压力特征的初始波形，近场静态压力特征提取位置要距离飞行器足够近，以保证CFD求解器本身存在的数值耗散不会影响结果的逼真度；与此同时提取位置还要距离飞行器足够远以保证可以忽略横流和升力效应的影响；

4)考虑斜率、压力增长和持续时间参数的影响，通过任意数量的线性线段来近似表示初始波形的形状，得到描述波形参数变化的一阶耦合微分方程组；

5)周期性求解描述波形参数变化的一阶耦合微分方程组，得到下一时刻波形的信息，所述下一时刻波形的信息包括传播方向，传播时间、高度和下一时刻波形的斜率、压力增长和持续时间，具体为：

描述波形参数变化的一阶耦合微分方程组为：

其中m_i是线段i的斜率

Δp_i是第i段和第i-1段连接处穿过激波的压力增长，无激波时为0；λ_i是第i段的持续时间ΔT；

C₁,C₂为过程变量，其表达式为：

其中γ＝1.4为常数，a₀为所在高度音速；p₀为所在高度大气压；ρ₀为所在高度密度；c_n＝c+v·n为波形传播速度；c为声速，v为风速，n是波阵面法向量；A为射线管面积； a₀,p₀,ρ₀,A,c_n是高度z的函数，所以沿着射线管变化；

不均匀大气有风条件下，对于任意波阵面形状传播的波来说，C₁和C₂沿着射线管是变化的，但是，如果假设这些量在很小的时间增长内是常数，则上述波形变形方程可积分得到解：

当：(1)

或者

时

(2)

或

时

近似得到：

通过上述关系式可以看出，当给定

后，关键就是计算出C₁,C₂，由C₁,C₂的表达式可以看出，计算的关键是得到所在高度a₀,p₀,ρ₀,c_n和A随时间变化的微小变化率，因此首先需要采用合适的计算射线路径的方法得到射线路径，然后得到沿射线路经a₀,ρ₀,c_n的时间变化率

射线管面积可以通过射线管上多个点的坐标和方向计算得到；

针对射线管初始值和初始位置、射线方向，射线管面积以及空间推进方法描述射线管推进过程：上一高度的射线位置和方向确定之后，就可以在空间推进到下一高度，在压力传播过程中，射线方向取决于所在高度音速和风速，当高度发生变化时，大气属性发生变化的同时，音速和风速随之变化，所以传播方向随着高度变化，下一射线点的传播方向为：射线上某点从上一高度传播到下一高度的过程中存在空间偏移量，下一高度射线的传播方向为上一高度的传播方向加上传播过程中的偏移量，将相对偏移量投影到空间三个方向上，并将所得分量加到原来方向向量上就得到新的射线点的传播方向；

6)重复执行步骤5)，直至计算高度降为指定高度，即得到指定高度处的过压值，作为音爆特征，评估音爆特征与地面的相互影响。

Images