基于重叠网格的翼型动态失速数值模拟

张团元,夏润泽,刘书岩,员海玮

(1. 海军航空大学,山东烟台264001;2.南京航空航天大学机械结构力学及控制国家重点实验室,江苏南京210016)

1 引言

动态失速是指在流场中进行周期振荡或临界机动动作等非定常运动的机翼/旋翼周围发生气流分离,从而造成动态迟滞的现象[1]。该现象在歼击机的过失速机动、直升机的桨叶旋转、风力机叶片的偏转等流动中非常普遍。其主要特征是翼型周围发生非定常分离流动,使得翼型发生失速直至气流重新恢复附着流动。动态失速现象的流动相较于静态失速而言,周围流场特征表现的更加复杂;升力、阻力、力矩特性的变化更为突出。在静态失速问题中,翼型迎角达到临界值时,由于气流分离气动力发生突变。但由于非定常分离涡及大尺度漩涡结构的存在,动态失速会表现出明显的气动力迟滞特性,气动力系数也极大的偏离了静态值,所以进行动态失速问题的相关研究变得尤为重要。

当前主要通过试验及仿真的方法研究动态失速问题。试验方法是在风洞中对缩比模型进行吹风试验得到气动力、压力等,进而直观认识到动态过程中的流动延迟现象。但其试验的前期准备和具体风洞试验过程非常复杂、成本高昂以及周期时间较长,并只能针对特定工况以及原理性的研究,很难捕捉到模型在振荡过程中脱体涡的运动形式。数值模拟计算可以克服以上试验中存在的不足和难题,针对模型进行多种工况下的动态失速模拟以及脱体涡整个运动过程的捕捉,该方法在动态失速问题的研究中已广泛应用和发展。

半经验模型方法和计算流体力学(CFD)方法是数值方法中研究非定常动态失速问题的两条重要途径。国外学者已经发展了一系列的半经验动态失速理论模型,如Gangwani模型[2]、ONERA模型[3]、Beddoes-Leishman(B-L)模型[4]。其中Gangwani模型是非常早期的半经验模型,而ONERA模型和Leishman-Beddoes(B-L)模型在二维气动力预测上较为成熟且应用广泛。以上两种方法都可以达到工程应用精度,并且大大缩短了计算时间。但半经验模型是基于特定模型试验数据实现的,不具有普适性。此外半经验模型方法仅在气动力模拟上具有较高的精度和效率,这种方法无法模拟出流场中复杂的分离现象,因此对于复杂流动的机理研究便无法更加深入。

计算流体力学(CFD)方法能够捕捉到动态失速过程中涡具体运动形式的细节,是动态失速问题研究的重要手段。国外学者L.Dubuc[5]等人采用隐式双时间步长方法求解Euler方程研究了NACA0012翼型跨音速小幅振荡的非定常运动,准确模拟了翼型浅失速状态下的迟滞效应。之后通过TFI方法开发一种网格变形算法,并通过襟翼的强迫振荡数值模拟验证了该方法的有效性[6]。国内学者王军利[7]采用改进的线性弹簧方法,同样是求解Euler方程计算了二维翼型以及三维机翼的非定常气动力问题。张兵[8]等跟据几何变形关系和网格的分块处理以及通过弹簧类比法和修正的TFI方法分别计算角点位移和子网格内部结点位移,形成了具有弹簧-TFI混合特征的动网格新方法。满洪海[9]以及雷延生[10]舍弃网格正交性,采用非结构网格通过弹簧近似法和局部网格重构方法实现网格变形。目前国内外在针对翼型非定常运动的研究中使用的动网格处理方法主要为网格变形法,如代数法、无限插值方法(TFI)、弹簧法等等[11]。其中TFI方法是一种基于多变量的代数插值方法,并且为保证网格的正交性只适用于结构网格下的小变形情况;弹簧法即将组成单元的线看作弹簧,通过与边长相关的公式计算得到弹簧刚度系数,从而将整个求解域构成一个弹簧网格。在边界节点运动后,求解网格节点处的静力平衡方程来计算新节点的位置,但其未考虑剪切力的作用,对于大变形来说容易出现网格畸变,从而产生负体积[12]。计算网格的正交性、光滑性以及近翼面贴体网格的布置是评价网格质量的重要因素,此外变形网格的方法需要满足几何守恒定律,这些都是影响计算精度的要素[13]。

针对网格变形及动态失速特性等问题,为避免网格变形带来的正交性下降和负体积的产生,本文采用重叠网格的方法及稳定的Roe格式进行动态失速数值模拟。验证了网格及计算方法的有效性;最后,根据数值模拟结果,分析了大迎角动态失速过程中脱体涡产生、分离以及再附着的运动过程以及翼型的气动特性。

2 模型及计算网格

2.1 物理模型

本文选用NACA0012对称翼型作为研究对象,其外形及俯仰运动形式如图1所示。

图1 NACA0012翼型运动形式示意图

翼型的运动形式如图1所示,绕翼型1/4弦线进行俯仰振荡,其迎角运动变化的规律为

α=α0+α1sin(kτ)

(1)

式中,α0是平均迎角,α1是幅值迎角。k是减缩频率,其值定义为

k=πfc/u∞

(2)

τ是无量纲化时间——减缩时间:

τ=2u∞t/c

(3)

式中,u∞——远前方来流速度;t——时间;f——振荡频率;c——翼型参考弦长。基于弦长的雷诺数由下式定义:

Re=ρu∞c/μ

(4)

式中,ρ——远前方来流空气密度;μ——动力黏度。

式(1)中,平均迎角α0是远前方来流的角度,幅值迎角α1为翼型的俯仰振荡角度,两项相加形成总迎角α的变化。

2.2 重叠网格

目前应用比较广泛的弹簧网格方法以及TFI插值方法虽然可以有效的避免负体积和畸形网格的出现,其对于小幅度的振荡运动是具有一定的准确度的,但当大变形和大幅度振荡的翼型运动是很难保证网格正交性和光滑性,产生网格畸变引起负体积导致无法计算。重叠网格在不改变网格形状的前提下,允许不同区域的网格进行独立并行计算,并通过插值的方式进行耦合。所以,本文采用重叠网格技术并通过对翼型动态失速问题进行研究与分析。

正因如此，GE Digital的“三步走”发展策略，即GE For GE、GE For Customers、GE For World是很务实的。在逻辑上，前两步都正确，但第三步踩了“知识壁垒”的红线，“For World”哪是那么容易？这中间恐怕省略了“For Other Industries”“For USA”等重要步骤，以及为了验证和优化这些重要步骤所必需的长期沉淀与反复打磨。GE Digital多跨出去的那一步，就是导致其业务停摆的短板——知识壁垒。

如图2所示,重叠网格由前景网格和背景网格两部分组成,通过数据插值将两部分网格合并后进行数值模拟分析。翼型周围的前景网格采用C-H型结构网格,以保证其良好的贴体性和正交性。背景网格采用正交性非常好的H型结构网格。计算网格形式及边界条件处理如图2所示。对翼型物面周围网格进行了加密处理,同时通过设置第一层网格高度保证y+≈1,为求解雷诺平均N-S方程提供捕捉粘性和大迎角气流分离提供良好的条件。

图2 NACA0012翼型重叠网格及局部放大示意图

本文不涉及翼型结构的弹性变形,所以动网格方面通过Fluent的用户自定义函数(User-defined functions)来控制重叠网格中前景区域的运动,从而实现控制翼型的简谐振荡运动。

3 离散格式

3.1 控制方程

为较好的模拟动态失速中的流动分离现象,控制方程采用对附着边界层湍流和分离湍流模拟具有良好精度的SST k-ω模型[14]。该模型借鉴了Johnson-King模型的思想,在计算中考虑到流动过程中的不平衡作用,得到的翼型动态失速的气动性能更加准确,适合于边界层流动、分离和转捩,因而采用该模型进行翼型动态失速的模拟分析。

在该模型中,涡粘性定义为

(5)

式中,Ω是涡量绝对值,Ω=|∂U/∂y|。F2是混合函数,表示为

(6)

其中,湍动能k运输方程和湍流比耗散率ω方程的具体形式是

τtijSij-β*ρωk

(7)

(8)

式(8)中,右端最后一项为交错扩散项。式中的生成项Pω≈γρΩ2。同时,识别函数F1是能够实现k-ε模型与k-ω模型之间切换,其定义为

F1=tanh{min[max(Γ1,Γ3),Γ2]}

(9)

Γ1,Γ2,Γ3可表示为

(10)

该模型下的参数β,γ,σk,σω可以表述为

φ=F1φ1+(1-F1)φ2(φ={σk,σω,β,γ})

(11)

其中,φ1和φ2分别适用于处理边界层内外流动,给定的系数值为

σk1=0.85,σω1=0.5,β1=0.075

σk2=1.0,σω2=0.856,β2=0.0828

(12)

γ可表示为

(13)

其它常数为

a1=0.31,β*=0.09,κ=0.41

(14)

SST k-ω湍流方程因其在预测动态失速迎角及气动力计算方面有较好的精度和效率[15],所以该方程目前应用较为广泛。本文在Fluent软件中,采用密度基求解器利用该湍流方程求解基于雷诺平均的N-S方程,从而得到俯仰振荡过程中的气动力及流场特性。

3.2 空间离散格式

目前,迎风(FDS、FVS)格式、中心格式、TVD格式和ENO格式是CFD计算中常用的格式。本文所采用的格式就是Godunov类求解器的FDS-ROE格式。

该方法是通过构造一个合理的常矩阵将非线性的Riemann间断分解问题转化为线性问题,其转化后的线性方程为

(15)

(16)

(17)

其中,k为自由参数,代表迎风格式的精度。Δ+、Δ-分别为向前、向后差分算子。

Roe格式一般被称为通量差分分裂(Flux Difference Splitting)格式,相对于其它计算格式,这种方法对激波和接触间断都具有较高的分辨率,同时也具有较好的稳定性。该格式已成为目前应用最广、评价最高的CFD空间离散格式之一。

4 非定常动态失速数值模拟

采用上述计算方法,通过小幅振荡及大迎角失速开展对动态失速特性的模拟研究,并将数值模拟结果与试验值[5],[17]进行对比分析,验证了计算方法和动网格的有效性。

计算所涉及的初始条件如下表所示。

4.1 小幅振荡非定常气动力计算

翼型做小幅振荡运动的计算条件为上表中的Case1所示。该状态的试验值在文献[5]给出。其它计算条件为:基于密度基求解器,采用压力远场边界,重叠区域边界设置为overset,壁面条件为无滑移边界条件,翼型绕其1/4弦线做正弦运动。计算过程为首先计算翼型在平均迎角下的定常流场,将定常流场的计算值作为非定常计算的初始值,以保证瞬态计算结果的准确性。计算结果如图3所示。

图3 NACA0012翼型非定常气动力系数曲线

图3分别给出了翼型运动过程中升力系数和力矩系数的变化规律。从图中可以看出,动态过程中升力与力矩系数与静态时明显不一致,随着迎角的变化均表现出非定常迟滞环现象。数值模拟的计算结果中除力矩系数中个别值以外,其它同试验值基本一致,验证了该方法对翼型小幅振荡下非定常气动力计算结果的精度。

4.2 大迎角动态失速特性分析

采用以上所述的数值方法模拟NACA0012翼型大迎角动态失速过程,试验条件和结果在文献[17]给出。算例的来流条件在表1的Case2中给出。雷诺数Re=3.45×106,静温T=288.15K,对应的黏度μ=1.7894×10-5Pa·s。其它条件设置和4.1节计算参数保持一致。

表1 NACA0012翼型数值模拟的计算条件

图4为计算结果中升力系数、阻力系数以及俯仰力矩系数随迎角变化的曲线,为保证计算结果的准确性同样与试验值进行了对比。在图中可以看出,在峰值处计算结果的绝对值略高于试验结果的绝对值,这是由于在动态失速问题中求解雷诺平均N-S方程时,对脱落涡所诱导的升力有较高的估算[18]。除峰值外,气动力等系数计算值与试验值吻合较好,能够较好的模拟低马赫数下的翼型大迎角动态失速状态下的气动特性。

图4 NACA0012翼型气动力系数迟滞回线

图5为NACA0012翼型在深失速俯仰振荡中翼型表面处不同迎角的流线图,较好的反映了翼型动态失速过程中前缘涡的产生、发展和脱落。对于NACA0012翼型而言,α在15°左右即会发生静态失速[19]。而动态失速的显著特点为失速迎角明显后移。

图5 NACA0012翼型的动态流场流线

通过图5可以看到,在α=15.47°↑时(↑和↓分别代表翼型上仰和下俯状态),翼型周围的流动为完全附着流动。翼型在上仰过程中,迎角逐渐增大,法向力系数也随之增大,同时翼型周围的流动开始出现变化,当迎角上仰角度大于16°时,如图5(b)所示开始在翼型前缘开始出现脱体涡,并且伴随着后缘涡的产生。随着迎角继续增大,流动在翼型前缘处发生分离,即前缘产生脱体涡,并向下游发展逐步演化成向下脱落的涡。当翼型运动到25°并开始作下俯运动时,翼型周围相对来流速度变大,前缘脱体涡和气流分离现象逐渐加剧,前缘涡沿翼型逐渐向后发展直至脱落,翼型的上下表面气流完全分离,导致升力系数骤然下降。与此同时,由于前缘涡往后缘发展并发生分离,导致压力中心后移,所以低头力矩系数会出现急剧上升现象。随着翼型下俯运动的继续,前缘再次产生脱体涡,所以图4(a)中升力系数会呈现一段上升现象。当运动到小迎角时涡消失,翼型表面恢复为附着流动。俯仰运动过程中正是由于前缘涡的产生、运动、脱落及再附着现象导致翼型在俯仰运动过程中同一迎角时而气动参数不同,表现出明显的迟滞现象,这是动态失速区别于静态失速的一个显著特征。

5 结论

本文基于重叠网格,采用雷诺平均N-S方程对NACA0012翼型进行了动态失速数值模拟,并通过动态流场分析了动态失速过程中脱体涡的运动特点。通过以上计算分析,主要的出以下结论:

1)本文所采用的动网格数值计算方法,由于不涉及网格变形,所以其具有良好的网格正交性,具有一定的精确度,可以有效的模拟翼型动态失速。该方法除了可以模拟小幅振荡,同时适合大迎角动态失速的计算分析。

2)动态失速过程中气动特性与静态失速中的不同,翼型的气动力会表现出明显的迟滞环现象。通过Case1和Case2可以看出迎角和减缩频率影响着迟滞环的形状和大小。

3)大迎角动态失速过程中,随着迎角的增大在前缘产生脱体涡,并随着上仰运动逐渐发展直至脱落,气流完全分离,导致升力骤降,低头力矩急剧上升。在迎角下俯到一定值时,气流会在前缘再附。气流分离与再附着的非定常效应导致翼型振荡过程中升力不对称,产生迟滞效应。