基于全局稳定性分析的高超声速飞行器转捩预测方法

摘要

本发明涉及一种基于全局稳定性分析的高超声速飞行器转捩预测方法，通过以下步骤实现：步骤一，计算得到目标高超声速飞行器表面的层流基本流场，得到三维边界层；步骤二，对步骤一得到的基本流场的三维边界层，在不同流向位置，进行二维全局稳定性分析，得到多个不稳定模态的增长率云图；步骤三，计算步骤二得到的每个不稳定模态在每个流向位置的最大增长因子，用Nmax曲线表示出来，并对不同模态的Nmax曲线再次取包络线得到Nmax_all曲线，表示出保守估计的不稳定模态的最大增长因子；步骤四，根据转捩判据Nt预测目标高超声速飞行器的转捩位置。

说明书

技术领域：

本发明涉及飞行器转捩预测技术领域，具体涉及一种高超声速飞行器边界层的转捩预测方法，尤其涉及一种基于全局稳定性分析的高超声速飞行器转捩预测方法。

背景技术：

在高超声速飞行条件下，层流区和湍流区的摩阻和热流相差很多，如果可以准确预测转捩位置，在设计高超声速飞行器热防护系统时，就可以避免完全按照湍流状态进行保守设计而导致的飞行器重量过大的问题，从而大幅度提高飞行器的飞行距离和机动性，对于高超声速飞行器的研发具有重要意义。

目前常采用的转捩预测方法有：基于线性稳定性理论的半经验的传统eN方法、求解抛物化稳定性方程的方法(PSE)、直接数值模拟(DNS)和转捩模式，其中，在工程应用中，eN方法通常被认为是预测转捩位置的一种有效方法。

eN方法是一种基于线性稳定性理论的转捩预测方法，传统上针对的是一维剖面。它是通过计算不稳定波的线性增长倍数来预测转捩的，其转捩判据，即N的临界值，往往依赖于实验或者实际经验给出，因此这是一种半理论半经验的转捩预测方法。

某一频率扰动的幅值沿流向演化过程中的放大倍数可以用放大因子N的e指数，即e

因此放大因子N值的表达式为：

其中，ξ

传统eN方法针对的是一维基本流剖面，在稳定性分析时仅考虑基本流剖面沿法向的变化。然而高超声速飞行器表面基本流场多存在基本沿法向和展向均剧烈变化的区域，同一流向位置下是二维基本流剖面，此时传统eN方法无法准确预测该区域的转捩位置。因此亟需提出一种针对高超声速飞行器的同时考虑基本流剖面沿法向和展向变化，即基于二维全局稳定性分析的转捩预测方法。

发明内容：

本发明的目的在于提供一种基于全局稳定性分析的高超声速飞行器转捩预测方法，该方法基于二维全局稳定性分析方法，考虑了基本流剖面同时沿法向和展向剧烈变化的特点，弥补了传统的eN方法在此处只考虑基本流剖面沿法向变化的缺陷，同时将多个不稳定模态的N

一种基于全局稳定性分析的高超声速飞行器转捩预测方法，通过以下步骤实现：

步骤一，计算得到目标高超声速飞行器表面的层流基本流场，得到三维边界层；

步骤二，对步骤一得到的基本流场的三维边界层，在不同流向位置，进行二维全局稳定性分析，得到多个不稳定模态的增长率云图；

步骤三，计算步骤二得到的每个不稳定模态在每个流向位置的最大增长因子，用N

步骤四，根据转捩判据Nt预测目标高超声速飞行器的转捩位置。

其中，步骤二中，在不同流向位置，进行二维全局稳定性分析的方法为：采用适用于高超声速飞行器三维边界层的无矩阵的二维全局稳定性分析方法，建立矩阵形式的可压缩无量纲N-S方程；将瞬时量表示为基本流和扰动量之和，用瞬时量的方程减去基本流的方程，再略去高阶小量，得到线性扰动方程；再使用局部近似平行流假设，将扰动表示成流向波的形式，最终得到无矩阵二维全局稳定性分析所需要的无频率的二维全局稳定性方程。

进一步地，在不同流向位置，进行二维全局稳定性分析的方法，按照下列的步骤执行：

1)建立笛卡尔坐标系下矩阵形式的可压缩无量纲N-S方程：

其中，x、y、z分别表示笛卡尔坐标系下的流向、法向和展向，φ表示瞬时量矢量φ＝(ρ,u,v,w,T)

2)将瞬时量表示为基本流与扰动量之和，即

其中变量φ′表示扰动矢量(ρ′,u′,v′,w′,T')

Γ、A、B、C、D、V

3)线性扰动方程(2)通过坐标变换得到曲线坐标系下的线性扰动方程：

其中，ξ、η、ζ分别表示曲线坐标系下流向、法向和展向，φ′表示扰动矢量φ'＝(ρ′,u′,v′,w′,T')

V

V

V

V

其中J表示Jacobi矩阵行列式，表达式为

4)推导得到无矩阵的二维全局稳定性方法所用的控制方程，即无频率的二维全局稳定性方程，采用局部平行流假设，将扰动设成流向波形式

其中

其中Γ、

进一步地，步骤二中，得到多个不稳定模态的增长率云图的方法为：采用Arnoldi方法计算二维全局稳定性方程(7)的时间模式特征值λ，通过特征值λ和频率ω的关系得到扰动的时间模式增长率ω

进一步地，步骤二中，对于方程(7)，特征值λ和频率ω的关系为

其中，步骤三N

其中，步骤四转捩判据Nt的确定方法为：转捩判据Nt通过飞行实验或地面风洞实验数据获得和修正，设某频率的扰动在向流向传播的过程中，幅值达到初始幅值的e

本发明提出一种基于全局稳定性分析的高超声速飞行器转捩预测方法，与现有技术相比的优点在于：

(1)高超声速飞行器表面多存在基本流场沿法向和展向均剧烈变化的区域，本发明可以预测该区域的转捩位置。

(2)本发明采用了二维全局稳定性分析方法，同时考虑了基本流剖面沿法向和展向的变化，具有物理依据，转捩预测结果比较可靠。

附图说明：

图1为本发明步骤框图

图2为目标高超声速飞行器模型示意图

图3为目标高超声速飞行器模型基本流流向速度云图

图4为飞行器模型基本流场边界层内6个最不稳定模态增长率云图

图5为6个最不稳定模态的N

图6为6个最不稳定模态的N

具体实施方式:

下面结合具体实施方式及附图对本发明进行详细说明。在下面的描述中，出于解释而非限制性的目的，阐述了具体细节，以帮助全面地理解本发明。然而，对本领域技术人员来说显而易见的是，也可以在脱离了这些具体细节的其它实施例中实践本发明。

在此需要说明的是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的设备结构和/或处理步骤，而省略了与本发明关系不大的其他细节。

本发明实施例提供了一种基于全局稳定性分析的高超声速飞行器转捩预测方法，如图1所示，包括如下步骤：

步骤一，计算得到目标高超声速飞行器的层流基本流场；

该步骤中，本实施例的目标高超声速飞行器为长短轴之比为3:1的椭圆锥模型，钝头半径为0.95mm，短轴半径为41mm，飞行马赫数Ma＝6，攻角为0度，侧滑角为0度，采用直接数值模拟(DNS)计算基本流场，无量化时使用的特征长度为钝头半径，特征速度为自由来流速度，飞行器模型示意如图2所示；

该步骤中，边界层内需要保证了法向网格数量为100以上，以保证后续全局稳定性分析的需求；计算所得基本流流场的流向速度云图如图3所示，由于对称性，只计算了1/4模型的结果；

步骤二，对步骤一得到的基本流场的三维边界层，在不同流向位置，进行二维全局稳定性分析，得到多个不稳定模态的增长率云图；

该步骤中，全局稳定分析的步骤如下：

1)从笛卡尔坐标系下矩阵形式的可压缩无量纲N-S方程出发，方程如下所示：

其中，x、y、z分别表示笛卡尔坐标系下的流向、法向和展向，φ表示瞬时量矢量φ＝(ρ,u,v,w,T)

2)将瞬时量表示为基本流与扰动量之和，即

其中变量φ′表示扰动矢量

3)方程(2)通过坐标变换得到曲线坐标系下的线性扰动方程

其中，ξ、η、ζ分别表示曲线坐标系下流向、法向和展向，φ′表示扰动矢量φ'＝(ρ′,u′,v′,w′,T')

V

V

V

V

其中J表示Jacobi矩阵行列式，表达式为

4)推导得到无矩阵的二维全局稳定性方法所用的控制方程，即无频率的二维全局稳定性方程，解的是时间模式问题。采用局部平行流假设，将扰动设成流向波形式

其中

将式(6)代入曲线坐标系下的线性扰动方程(3)中，整理得到无矩阵的全局稳定性方法所用的控制方程，即无频率的二维全局稳定性方程

其中Γ、

5)采用Arnoldi方法计算二维全局稳定性方程(7)的时间模式特征值λ，通过特征值λ和频率ω的关系得到扰动的时间模式增长率ω

该步骤中，Arnoldi方法是本领域公知的方法，可通过调用开源的ARPACK程序包实现，具体参考：R.B.Lehoucq,D.C.Sorensen,C.Yang.ARPACK Users Guide:Solution ofLarge Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods.[J].Journal of Chemical Physics,2015,157(S14):044119。

该步骤中，特征值λ和频率ω的关系为

取固定时间间隔T

其中L(T

其中λ为特征值，

其中ω的实部ω

综上所述，特征值λ和频率ω的关系为

该步骤中，Gaster变换的形式为：

-α

其中c

该步骤中，对于目标高超声速飞行器模型，得到了6个最不稳定模态的增长率沿流向分布的云图，分别是Y-Mode1、Y-Mode2、Y-Mode3、Z-Mode1、Z-Mode2、Z-Mode3,如图4所示；

步骤三，计算步骤二得到的每个不稳定模态在每个流向位置的最大增长因子，用N

该步骤中，认为边界层内不同频率的扰动可以持续增长到转捩位置，给定频率小扰动的幅值从起始位置至流向某位置的放大倍数由式(14)给出：

其中ξ

该步骤中，基于多个不稳定模态的增长率云图，在不同频率下，将所述最不稳定模态的增长率沿流向进行积分，可以获得各频率扰动的幅值放大因子N值沿流向坐标的分布，用一条曲线将不同频率波N值沿流向坐标分布的最大值包络起来，包络线代表各种频率扰动在此处可能达到的最大幅值放大因子，即得到N

步骤四，根据转捩判据Nt预测目标高超声速飞行器表面的转捩位置；

该步骤中，认为只要某频率的扰动在向流向传播的过程中，幅值达到初始幅值的e

该步骤中，本实施例的转捩判据Nt为8.5，通过已知高超声速飞行器模型的风洞实验结果和数值计算结果标定得到。

该步骤中，本实施例中结合目标高超声速飞行器模型的N

如上针对一种实施例描述和/或示出的特征可以以相同或类似的方式在一个或更多个其它实施例中使用，和/或与其它实施例中的特征相结合或替代其它实施例中的特征使用。

应该强调，术语“包括/包含”在本文使用时指特征、整件、步骤或组件的存在，但并不排除一个或更多个其它特征、整件、步骤、组件或其组合的存在或附加。

本发明以上的装置和方法可以由硬件实现，也可以由硬件结合软件实现。本发明涉及这样的计算机可读程序，当该程序被逻辑部件所执行时，能够使该逻辑部件实现上文所述的装置或构成部件，或使该逻辑部件实现上文所述的各种方法或步骤。本发明还涉及用于存储以上程序的存储介质，如硬盘、磁盘、光盘、DVD、flash存储器等。

这些实施例的许多特征和优点根据该详细描述是清楚的，因此所附权利要求旨在覆盖这些实施例的落入其真实精神和范围内的所有这些特征和优点。此外，由于本领域的技术人员容易想到很多修改和改变，因此不是要将本发明的实施例限于所例示和描述的精确结构和操作，而是可以涵盖落入其范围内的所有合适修改和等同物。

本发明未详细说明部分为本领域技术人员公知技术。