翼型动态失速气动力二次峰值数值模拟研究

井思梦，赵国庆，招启军

（南京航空航天大学直升机旋翼动力学国家级重点实验室，南京 210016）

动态失速是一种非定常流动分离现象，通常发生在大拉力或机动飞行状态的直升机旋翼上。动态失速一旦发生，会引起升力突降、旋翼反扭矩和桨叶振动突增等问题，从而限制直升机的最大飞行速度和机动性能。因此，动态失速机理一直是旋翼空气动力学领域的研究热点之一。

迄今，许多学者针对翼型开展了动态失速机理的研究［1⁃3］。通过试验测量翼型的气动载荷以及流场特性，可以直观地了解翼型动态失速过程中气动力的迟滞效应以及动态失速涡的发展。McCros⁃key 等［4］测量了NACA0012 翼型及其前缘修型翼型的动态失速特性，发现涡脱落现象是不同类型动态失速的主要共同特征。Wang 等［5］采用PIV 技术测量了OA209 和SC1095 翼型前缘涡的输运速度，研究表明前缘涡的输运速度主要受到翼型振荡频率的影响，当振荡频率增加时，前缘涡的输运速度也会提升。Geissler 等［6］开展了OA312 翼型动态失速特性的试验和数值模拟研究，发现涡量的发展、脱落和积累对翼型的动态失速特性有重要影响。

由于翼型动态失速的风洞试验较为复杂，成本高，且受到测量设备、技术的限制，因此只能开展有限工况下的研究。随着计算流体力学的发展，数值模拟逐渐成为翼型动态失速特性研究的主要途径［7⁃10］。

为了研究翼型外形对动态失速特性的影响规律，一些学者开展了相应的研究。王清等［11］和邹锦华等［12］开展了翼型前缘外形对动态失速特性的影响规律研究，结果表明：上翼面前缘外形对翼型的动态失速特性有显著影响，并且合理的变形能够抑制翼型的动态失速。王友进等［13］针对翼型厚度对动态失速特性的影响开展了数值模拟分析，研究发现：薄翼型的动态失速是由前缘分离引起的，而厚翼型的动态失速则是由后缘分离引起的，并在分离区域向翼型前缘扩展的过程中发生失速。Shar⁃ma 等［14］模拟分析了NACA 系列不同厚度对称翼型在低雷诺数下的动态失速特性，结果表明：较薄的NACA0009 翼型的动态失速是由层流分离泡的破裂引起的，而最厚的NACA0018 翼型的动态失速则是由湍流边界层的分离引起的。

此外，气动参数以及翼型运动参数对翼型的动态失速特性也有显著影响。Kim 等［15］和孔卫红等［16］开展了压缩性效应对翼型动态失速特性的影响研究，分析表明：小马赫数下压缩性效应的影响较小，而在高马赫数下，压缩性效应对翼型的动态失速特性起主导作用。宋辰瑶等［17］、赵国庆等［18］和杨鹤森等［19］对翼型动态失速特性开展了参数影响研究，揭示了平均迎角、迎角振幅及缩减频率等振荡参数以及雷诺数和马赫数等气动参数对气动力迟滞效应以及气动力峰值的影响规律。

直升机旋翼翼型工作的缩减频率一般在0.03～0.15 之间，赵国庆等［18］针对缩减频率为0.05，0.10 和0.15 等3 个典型工况下的旋翼翼型动态失速特性开展了数值模拟研究。分析表明：缩减频率表征了旋翼翼型动态失速特性的强弱程度，对旋翼翼型动态失速的迟滞效应、气动力系数峰值和失速迎角具有重要影响。当缩减频率为0.05 时，除了动态失速涡引起的气动力峰值外，翼型的气动力系数曲线还出现了二次峰值，然而文中并未开展对气动力二次峰值的细致研究。Mcalister 等［20］针对翼型动态失速的试验研究也表明：在某些动态失速状态下，翼型的气动力会出现明显的二次峰值。Choudhry 等［21］对动态失速升力特征进行了系统的综述研究，分析了低雷诺数下翼型以恒定速度抬头至指定迎角过程中的动态失速现象和升力变化，发现了由二次涡结构引起的升力二次峰值现象。

目前，在翼型动态失速机理方面，相关学者已经做了大量的研究工作。然而，针对翼型动态失速气动力二次峰值的研究工作相对较少。因此，本文建立了翼型非定常流场数值模拟方法，在此基础上开展了翼型动态失速气动力二次峰值的机理研究，并着重分析了翼型厚度、弯度及弯度位置等参数对气动力二次峰值的影响，获得了一些有意义的结论。

1 翼型非定常流场数值模拟方法

1.1 网格生成方法

通过求解泊松方程来生成围绕翼型的C 型网格，将泊松方程离散化后可以得到

图1 给出了围绕SC1095 翼型的C 型网格。从图1 可以看出，翼型网格具有良好的正交性和贴体性。

图1 围绕SC1095 翼型的C 型网格Fig.1 C⁃type grid around SC1095 airfoil

为了模拟翼型俯仰振荡时的非定常气动特性，采用运动嵌套网格技术实现翼型网格的运动以及翼型网格和背景网格之间的信息交换。采用Hole⁃Map 方法确定背景网格在翼型网格上的洞边界，采用Inverse⁃Map 方法进行背景网格人工内边界的贡献单元搜索。图2 给出了运动嵌套网格系统示意图。

图2 运动嵌套网格系统Fig.2 Moving⁃embedded grid system

1.2 流动控制方程及求解方法

采用积分形式的RANS 方程作为翼型流场的控制方程

式中：ρ表示密度；p表示压强；E和H分别表示单位质量的总能和总焓；Vr和Vt分别表示相对运动速度和由网格运动引起的牵连速度，且Vr=VVt，其中V表示绝对运动速度；ix与iy分别为x与y方向的单位矢量；τij为黏性应力；Θi为描述黏性和热传导的作用项。

采用CLORNS 求解器［22］对流场控制方程进行求解。空间离散采用Roe 格式，网格面两侧的流场原始变量采用三阶MUSCL 格式进行重构，时间推进采用LU⁃SGS 隐式格式。为了较好地模拟大迎角下翼型的气流分离和失速特性，采用S⁃A 湍流模型模拟流场的黏性系数，实现控制方程的封闭。

1.3 方法验证

采用所建立的CFD 方法对NACA0012 翼型在深度失速状态下的气动特性进行了模拟。迎角变化规律为α=14.91°+9.88°sin（2kt），来流马赫数Ma为0.283，缩减频率k=0.151，雷诺数Re=3.4×106，t为时间。在该算例中，计算域设置如图3 所示，翼型表面采用无滑移和绝热条件，背景网格外边界采用无反射远场边界条件。在划分翼型网格时，设置壁面第一层网格单元高度为4×10-6c，保证y+≈1。图4 给出了3 套不同网格的气动力系数计算值与试验值［20］的对比，3 套网格的信息在表1中给出。从图4 中可以看出：计算值与试验值吻合较好，表明本文的数值模拟方法能够准确地模拟翼型的动态失速特性；并且，在翼型上仰过程中，3 套网格的计算结果一致，在下俯过程中，中等网格和密网格的计算结果基本重合。因此，在后文的研究中，均采用了与此算例中等网格相同的网格量。

表1 NACA0012 翼型网格信息Table 1 Information about the three sets of grids

图3 计算域示意图Fig.3 Schematic diagram of computational domain

图4 翼型气动力系数计算值与试验值[20]对比Fig.4 Comparison of aerodynamic force coefficients be⁃tween calculated results and test data[20]

2 动态失速气动力二次峰值分析

为了探究翼型动态失速状态下气动力二次峰值的发生机理，对NACA0012 翼型在状态Ma=0.293，α=14.82°+9.90°sin（2kt），k=0.051，Re=3.82×106下的动态失速特性进行了模拟分析。

图5 给出了NACA0012 翼型的升力、阻力和力矩系数曲线，图中A～F为6 个不同时刻。从图5中可以看出，NACA0012 翼型在此状态下发生了深度失速，并且失速后升力、阻力和力矩系数出现了二次峰值。

图5 NACA0012 翼型气动力系数Fig.5 Aerodynamic force coefficients of NACA0012 airfoil

图6 给出了图5 中6 个不同时刻（A～F）NA⁃CA0012 翼型的压强系数分布云图和流线图。不同时刻对应的翼型迎角如表2 所示。不同时刻翼型表面的压强系数分布如图7 所示。

图6 NACA0012 翼型压强分布云图和流线图Fig.6 Pressure coefficient cloud and freestream lines of NACA0012 airfoil

图7 NACA0012 翼型压强系数分布Fig.7 Pressure coefficient distribution of NACA0012 airfoil

表2 不同时刻的翼型迎角Table 2 Angles of attack of airfoil at different time

从图6 可以看出，在A时刻，翼型上表面已经形成了动态失速涡（Dynamic stall vortex，DSV），而且翼型上表面中部出现了由动态失速涡导致的局部负压区域。动态失速涡在翼型前缘附近形成后，会沿着翼型上表面向后缘移动。一方面，当动态失速涡附着于翼型表面时，涡致升力的存在会使得翼型升力系数持续增加；另一方面，负压中心也在向后缘移动，导致低头力矩逐渐增加。接着，动态失速涡会从翼型表面脱落，涡致升力开始减小，同时，前缘负压峰值减小，二者共同作用导致升力系数下降，阻力、低头力矩进一步增加。在A时刻，动态失速涡即将从翼型表面脱落，此时翼型的升力系数达到峰值，阻力和低头力矩系数接近峰值，如图5 所示。

在B时刻，动态失速涡诱导形成了后缘涡（Trailing⁃edge vortex，TEV），导致后缘附近局部负压峰值的出现，如图7（a）所示。在主流和后缘涡的作用下，动态失速涡向上抬升，脱离翼型表面，使得由动态失速涡引起的局部负压减小，进而导致升力系数的突降。从图7（a）可以看出，与A时刻相比，B时刻翼型前缘负压峰值减小，同时翼型中部动态失速涡引起的局部负压减小，并且后缘涡引起了新的局部负压，使得翼型的阻力系数和低头力矩系数减小。

在C时刻，动态失速涡的强度进一步减弱，在后缘涡的作用下，动态失速涡又逐渐向翼型上表面靠拢。后缘涡上洗并膨胀，强度也逐渐减弱。此时，翼型的升、阻力系数和低头力矩系数进一步减小。

在D时刻，后缘涡上洗，从翼型后缘脱落并向下游移动至尾迹区，由后缘涡引起的局部负压消失。由于后缘涡的诱导作用，动态失速涡继续向翼型上表面靠拢，并重新附着于翼型上表面。此时，动态失速涡较弱，在翼型表面并没有形成明显的局部负压。但是，如图7（b）所示，与C时刻相比，D时刻的前缘负压峰值和翼型上表面负压明显增加，导致升力、阻力和力矩系数回升。

由图7（b）可知，在E时刻，翼型前缘负压峰值和上表面负压较D时刻进一步增加。此时，翼型升、阻力系数和低头力矩系数达到二次峰值。同时，从图6（e）可以看出，由于动态失速涡的诱导作用，翼型后缘附近形成了新的较弱的后缘涡。

在F时刻，动态失速涡向上抬升，远离翼型表面，分离区域扩大，气流分离加剧，导致前缘负压峰值和翼型上表面负压的减小，从而导致升力、阻力和力矩系数的下降。

总体而言，动态失速涡的形成，导致第一个气动力峰值的出现；动态失速涡诱导形成后缘涡，在后缘涡的作用下，动态失速涡向翼型表面靠拢，并重新附着于翼型上表面，导致翼型气动力二次峰值的出现。

3 翼型外形对气动力二次峰值的影响分析

为了探究翼型外形对气动力二次峰值的影响，对NACA 系列不同厚度、弯度和弯度位置的翼型开展了气动特性模拟分析。计算状态同第2 节。

3.1 翼型厚度影响

为了探究翼型厚度对气动力二次峰值的影响，对NACA0009、NACA0012 和NACA0015 这3 个不同厚度翼型的动态失速特性进行了模拟分析。翼型外形对比如图8 所示。

图8 不同厚度翼型外形对比Fig.8 Comparison of airfoils with different thicknesses

图9 给出了NACA0009、NACA0012 和NA⁃CA0015 这3 个不同厚度翼型的气动力系数对比。从图中可以看出，翼型厚度对动态失速特性有显著的影响。翼型厚度增加，动态失速迎角增加，升力、阻力和低头力矩系数峰值均增加。在线性升力段，相同迎角下厚度较大的翼型产生的升力较小，且厚度对线性段升力线斜率的影响较小。此外，3 个不同厚度翼型在此动态失速状态下均出现了气动力二次峰值。本文称两个峰值之间的气动力最小值为谷值，由图9 可知，翼型厚度越大，二次峰值出现的越晚，且二次峰值相对谷值的增量（下文简称“二次峰值增量”）越大。

图9 不同厚度翼型的气动力系数对比Fig.9 Comparison of aerodynamic force coefficients of airfoils with different thicknesses

图10 给出了3 个翼型的压强系数分布云图和流线图。从图10 中可以看出，在翼型上仰过程中迎角为18.00°的状态下，NACA0009 翼型上表面已经形成了动态失速涡，且动态失速涡移动至后缘附近，诱导形成了后缘涡。此时，NACA0012 翼型前缘形成了动态失速涡，而NACA0015 翼型的流动仍附着在翼型表面，没有发生明显的气流分离。当迎角增加至20.23°，NACA0012 翼型上表面的动态失速涡移动至后缘附近，诱导形成了后缘涡。此时，NACA0015 翼型后缘发生了气流分离。当迎角继续增加至21.78°，NACA0015 翼型上表面形成了动态失速涡，NACA0012 翼型的后缘涡脱落至尾迹区。

图10 不同厚度翼型的压强分布云图和流线图对比Fig.10 Comparison of pressure coefficient cloud and streamlines of airfoils with different thicknesses

NACA0009、NACA0012 和NACA0015 翼型上表面先后形成了不同强度的动态失速涡和后缘涡，动态失速涡和后缘涡相互作用，导致翼型气动力出现二次峰值。翼型厚度增加，推迟了动态失速涡和后缘涡的形成，从而推迟了气动力二次峰值的出现。此外，3 个翼型在此状态下均发生了前缘失速，翼型厚度增加，形成的动态失速涡和后缘涡强度增加，导致气动力二次峰值增量增加。

3.2 翼型弯度影响

为了探究翼型弯度对翼型动态失速状态气动力二次峰值的影响，对NACA0012、NACA2212 和NACA4212 这3 个不同弯度翼型的动态失速特性进行了模拟分析。翼型外形对比如图11 所示。

图11 不同弯度翼型外形对比Fig.11 Comparison of airfoils with different cambers

图12 给出了NACA0012、NACA2212 和NA⁃CA4212 这3 个不同弯度翼型的气动力系数对比。从图12 中可以看出，弯度对翼型动态失速特性有明显影响。迎角相同时，弯度较大的翼型产生的升力系数越大，低头力矩系数也越大。随着翼型弯度的增加，动态失速迎角、阻力和力矩发散迎角增加，气动力系数峰值呈现先增加后减小的趋势。NA⁃CA0012 和NACA2212 翼型的气动力系数出现了明显的二次峰值，而NACA4212 翼型的二次峰值则较小。而且，翼型弯度越大，翼型升力系数二次峰值出现的越晚，二次峰值增量呈现先增加后减小的趋势，并且与NACA0012 和NACA2212 翼型相比，NACA4212 翼型的升、阻力和力矩系数二次峰值增量均最小。

图12 不同弯度翼型的气动力系数对比Fig.12 Comparison of aerodynamic force coefficients of airfoils with different cambers

图13 给出了3 个翼型在不同时刻的压强系数分布云图和流线图。从图13 中可以看出，NA⁃CA0012 和NACA2212 翼型先后发生了前缘失速，即动态失速涡在翼型前缘附近形成。而NA⁃CA4212 翼型则先出现了后缘气流分离，随着迎角的增加，分离点向上游移动，并在翼型中后部形成了动态失速涡。由于动态失速涡的诱导作用，3 个翼型后缘附近形成了不同强度的后缘涡。当翼型发生前缘动态失速时，翼型弯度的增加会推迟动态失速涡和后缘涡的形成，从而推迟二次峰值的出现；且弯度增加会增强动态失速涡和后缘涡的强度，从而导致二次峰值增量增加。而当弯度增加到一定程度，翼型上表面后缘附近的曲率变化较大，导致翼型失速类型转变为后缘失速，与前缘失速情况相比，后缘失速情况下形成的动态失速涡和后缘涡的强度较小，因此气动力系数二次峰值增量较小。

图13 不同弯度翼型的压强系数分布云图和流线图Fig.13 Comparison of pressure coefficient cloud and streamlines of airfoils with different cambers

3.3 翼型弯度位置影响

为了探究翼型弯度位置对动态失速气动力二次峰值的影响，对NACA2212、NACA2412 和NACA2612 这3 个不同弯度位置翼型的动态失速特性进行了模拟分析。翼型外形对比如图14所示。

图14 不同弯度位置翼型外形对比Fig.14 Comparison of airfoils with different positions of camber

图15 给出了3 个不同弯度位置翼型的气动力系数对比。在线性升力段相同迎角下，弯度位置靠后的翼型产生的升力系数较大，低头力矩系数也较大。随着弯度位置的增加，翼型的动态失速迎角、阻力和力矩发散迎角减小，气动力系数峰值也有所减小。弯度位置越靠后，二次峰值出现得越早，且二次峰值增量越小。

图15 不同弯度位置翼型的气动力系数对比Fig.15 Comparison of aerodynamic force coefficients of airfoils with different positions of camber

图16 给出了3 个翼型的压强系数分布云图和流线图。从图16 中可以看出，3 个翼型都发生了前缘失速，且弯度位置越靠后，动态失速发生的越早，后缘涡的形成也越早，导致二次峰值出现的越早。而且，弯度位置越靠后，动态失速涡和后缘涡的强度越小，从而导致翼型气动力系数二次峰值增量越小。

图16 不同弯度位置翼型的压强系数分布云图和流线图Fig.16 Comparison of pressure coefficient cloud and streamlines of airfoils with different positions of camber

4 结论

基于URANS 方程，结合运动嵌套网格技术、LU⁃SGS 隐式格式、Roe⁃MUSCL 数值格式和S⁃A湍流模型，建立了翼型非定常气动特性数值模拟方法，开展了翼型动态失速气动力二次峰值机理研究和参数影响分析。综合计算分析结果，得出结论如下：

（1）在动态失速过程中，翼型上表面会形成动态失速涡，导致第一个气动力峰值的出现；动态失速涡诱导形成后缘涡，在后缘涡的作用下，使得动态失速涡向翼型表面靠拢，并重新附着于翼型上表面，从而导致翼型气动力二次峰值的出现。

（2）翼型厚度增加、弯度位置前移以及一定范围内弯度的增加，会推迟动态失速涡和后缘涡的形成，并增加动态失速涡和后缘涡强度，从而使得气动力二次峰值推迟，二次峰值增量增加。而当翼型弯度较大时，翼型的失速方式从前缘失速变为后缘失速。与前缘失速情况相比，后缘失速情况下形成的动态失速涡和后缘涡的强度较小，因此气动力系数二次峰值增量较小。