基于DMD方法的翼型大迎角失速流动稳定性研究

2018-06-29 11:06:32叶正寅

空气动力学学报 2018年3期

叶坤, 叶正寅, 武洁, 屈展

(西北工业大学翼型叶栅空气动力学国家级重点实验室, 陕西西安 710072)

0 引言

失速是航空航天领域中一种普遍存在的流动现象，失速研究具有重要的学术和应用价值。传统认为静态失速是当迎角增大到某一定值时，在翼型上表面会出现大规模分离流动，从而导致升力突然下降和阻力突然增加及非定常压力脉动放大等负面效应，这对飞行器的整体气动特性和操纵性能产生严重影响[1]。随着航空航天新技术的不断发展，对现代飞行器的性能要求越来越高，如大迎角机动性，苏27过失速“眼镜蛇”机动时迎角超过了120°，用于改善飞机起降性能的增升襟翼也经常在大迎角条件下工作，特别是利用迎角40°左右存在的二次升力峰值和阻力是着陆时增升增阻的一种理想选择[2]。因此，只有深入理解大迎角复杂分离流动演化规律，才能更好地控制失速现象，使飞行器可用迎角逾越静态失速这一传统禁区。

深入理解大迎角复杂分离流动演化规律，首先需要获取较准确的流场。近年来，LES方法和DES方法被广泛地应用于分离流动的模拟中[3-7]。在获取较准确的流场信息后，通常采用非线性理论、模态分解等方法对流动结构进行分析。从非线性动力学角度讲，静态失速中存在滞后、突跳及分岔现象。Mittal等[8]通过数值模拟对翼型静态失速迎角附近的迟滞现象进行了研究，并发现迟滞现象的发生是由于流动具有对历史行为的记忆能力。张家忠等[9]对翼型失速的非线性动力学特性进行研究，发现静态失速是一种“鞍—结”分岔现象。另一方面，在获得流动结构的众多方法中，本征正交分解方法(Proper Orthogonal Decomposition， POD)被广泛使用。自Lumley[10]首先将POD方法引入湍流研究以来， Sirovich[11]引入了snapshot 方法来研究波动流的动力学问题。康伟等[12]基于POD方法，采用非线性动力学理论对翼型绕流的多模态耦合机制进行分析，研究了模态耦合作用与流动稳定性的关系。

翼型大迎角分离涡的演化包含了涡量生成、积聚和耗散，其中的流动机理非常复杂。本文采用DMD方法对分离涡的演化过程进行稳定性分析。首先基于SST-DES方法对NACA0012进行数值模拟得到非定常流场数据。实际翼型大迎角分离流动会存在明显的三维效应，且通常基于DES方法进行计算时，对物面展向长度和网格数量有一定的要求，文献[22]中对此进行了较详细的分析。但受计算资源的限制，本文数值模拟中网格展向仅拉伸一层。事实上，对于一些经典的三维分离流动，例如圆柱、凹腔等流动，大部分学者及文献也均是从二维的角度进行计算和分析。

1 计算方法及验证算例

三维非定常N-S方程在直角坐标系中的积分守恒形式为:

(1)

其中Ω为控制体，∂Ω为控制体单元边界，Q为守恒变量，Q=[ρ,ρu·ρv,ρw,ρE]T，F(Q)为无黏通量，FV(Q)为黏性通量。其中ρ、P、T、E分别为密度、压强、温度和单位质量流体的总能，u、v、w分别为x、y、z方向的速度分量。

采用基于剪切应力输运SSTk-w两方程湍流模型的SST-DDES混合方法对流场进行模拟，采用有限体积方法进行空间离散，空间格式采用Roe，时间推进采用LU-SGS。

目前公开的翼型大迎角实验结果较少，本文采用文献[23]中NACA0012的实验数据。实验条件为Ma=0.15、Re=6×106。通常采用DES方法进行计算时，对物面展向长度和网格数量有一定的要求，文献[22]中对此进行了较详细的分析。但受计算资源的限制，数值模拟中网格展向仅拉伸1层。图1为生成的结构网格，翼型弦长为1 m。为了避免使用拉伸比很大的“C”型网格，采用420×180的“O”型结构网格。其中翼型上表面300个点，下表面120个点，物面网格第一层高度为5×10-6m，满足y+<1，附面层径向采用40层网格，增长率为1.25。为准确模拟分离，对上翼面附近区域的网格进行加密，网格层数为100层，增长率为1，之后采用40层过渡到远场，远场距离物面约为100倍弦长。非定常计算中的时间步长为5×10-4s。

图2为升力系数(时间平均下的升力系数)随迎角的变化曲线。实验数据中，迎角为17.13°时，升力系数达到最大，当迎角为18.21°时，升力系数突然大幅下降。本文计算结果中，当迎角为17°时，升力系数达到最大。为进一步细致了解17°至18°之间的流动，在迎角17.5°至18°之间，每隔0.1°进行一次非定常计算。当迎角由17.7°进一步增加到17.8°时，升力系数突然大幅下降，且开始呈现周期性变化。图3为17.7°和17.8°的升力系数随时间的变化，图4为流场流线和涡量云图。上述结果表明本文的计算结果与实验数据吻合较好，说明本文计算结果是可信的。

2 DMD方法

DMD方法只基于流动快照而不受模型限制，通过提取流动中的模态，准确地描述流动结构。对于线化流动而言，被提取的模态可以用来表征全局稳定模态。对于非线性流动而言，DMD的模态描述了在数据序列中起支配作用的流动结构。以下为其数学推导过程。

(2)

第一步，假设线性映射A是快照vi到vi+1的映射，即:

vi+1=Avi

(3)

且在全部区域和整个时间段内的采样都满足上述映射关系。对于缓慢变化的系统，乘法准则是上述假设的基础。

={Av1,Av2,Av3,…,AvN-1}

(4)

A≈UFU*

(5)

(6)

其中，Σ是一个r×r的对角矩阵，有非0的对角元素{σ1,…,σr}，并且，U和V为酉矩阵。

U∈CM×r，并且U*U=I

V∈Cr×N，并且V*V=I

将式(5)和式(6)代入式(4)中，可得：

(7)

F即为A的最优低维估计矩阵，由此通过求解F的特征值和特征向量即可得到DMD分析结果。

假设F是一个r维子空间内的精确映射矩阵而不是对A的估计，那么

xk+1=Fxk

(8)

则A的POD模态U可以看作从xk到vk的近似映射，则

vk≈Ux

(9)

如果F的特征向量个数与秩相等且都线性独立，那么F可以写做如下形式

(10)

(11)

(12)

(13)

由此可得：

[v1,v2,…,vN-1]≈[φ1,φ2,…，φr]·

(14)

动态模态分解获得的最重要的结果是模态φi和模态对应的特征值μi,lg(μi/Δt)的实部为对应模态的

放大率，其虚部表示对应模态的频率。在进行稳定性判断时，放大率为正，则对应的模态发散；放大率为负，则对应的模态收敛；若放大率为0，则对应的模态为稳定极限环模态。同样也可以将μi置于复平面的单位圆上进行判断，若在单位圆内，则收敛；单位圆外，则发散；单位圆上，则为稳定周期性模态。

3 计算结果分析

本文在验证算例中实验条件的马赫数和雷诺数下对NACA0012大迎角失速翼型流场的涡稳定性演化规律进行研究。首先，对迎角在17.8°至40°范围内的状态进行非定常计算，其中22°至24°的计算中，迎角间隔为1°。图5为Strouhal数(St)随迎角的变化曲线，

St=fC/U∞

(15)

式中f为频率，C为弦长，U∞为来流速度。

可以看出，St随迎角的增加并不是呈现单调减小的规律。当迎角在17.8°至23°范围内，随迎角的增加，St基本呈线性减小；当迎角由23°增加到24°时，St大幅减小；当迎角在24°至29°范围内，St随迎角增加缓慢减小，并且迎角为29°时，St达到最小，为0.166168；迎角由29°增加到30°时，St有微幅的增加；当迎角进一步由30°增加到31°时，St大幅增加；迎角从31°增加到40°时，St再次基本呈线性减小。因此，采用DMD方法分别对以下3种迎角状态下的分离流场进行稳定性分析：(1) 迎角为17.8°时，流动刚刚进入失速状态，称为失速初始状态；(2) 迎角增加到24°，流动进一步分离，St突然大幅减小，称为浅失速状态；(3) 迎角进一步增加到29°时，St达到最小，称为深失速状态。

由于在上述三种迎角下，流场最后均表现为稳定极限环状态，升力系数均呈现周期性变化。因此，采用DMD方法进行稳定性分析时，DMD采样中仅提取一个周期的瞬态流场数据，三种迎角下的样本数量分别为65、182、234。

3.1 DMD稳定性分析

图6为升力系数随时间的变化。三种状态下，升力系数在一个周期内的变化相差较大。失速初始状态下，升力系数一个周期内呈现类似正弦曲线的变化，周期性变化规律简单，一个周期中仅有一个极大值和一个极小值，最大值与最小值相对均值的振荡幅值分别为28.6%和32.2%。浅失速状态中，升力系数一个周期内的变化较复杂，一个周期中存在四个极大值和四个极小值，在5/8T时刻附近达到的极大值略小于1/8T时刻的最大值，最大值与最小值相对均值的振荡幅值分别为24.4%和18.6%。深失速状态中，升力系数在一个周期中存在两个极大值和两个极小值，且在5/8T时刻附近达到的极大值小于1/8T时刻的最大值。最大值与最小值相对均值的振荡幅值分别为38.3%和17.6%。

图7为升力系数FFT分析结果。可见FFT可提取其主频及高阶倍频，图中仅标出前五阶频率，按照St对应的峰值对St进行排序，三种状态的顺序各不相同。失速初始状态，St的顺序为：St1，St2，St3，St4，St5；浅失速状态，St的顺序为：St2，St4，St1，St3，St5；深失速状态，St的顺序为：St2，St1，St4，St3，St5。

表1为DMD方法提取频率与FFT分析结果的对比，两者之间频率最大差异为0.15037%，说明DMD方法可以非常准确地提取流动中的频率。