旋翼翼型非定常动态失速特性的CFD模拟及参数分析

赵国庆，招启军，王清

旋翼翼型非定常动态失速特性的CFD模拟及参数分析

赵国庆，招启军，王清

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

构建了一套基于运动嵌套网格技术和可压缩RANS方程的旋翼翼型非定常流动特性模拟的高效、高精度的CFD方法。首先，发展了基于Poisson方程求解的围绕翼型的粘性贴体正交网格生成方法，并提出了基于最小距离法(MDM)改进策略的运动嵌套网格生成方法，克服了弹簧法可能导致网格畸变的不足;其次，为准确模拟由湍流分离和气流再附引起的气动力的迟滞效应，基于RANS方程、双时间方法和高阶插值格式，建立了旋翼翼型非定常气动特性分析的高精度数值方法，并采用能够较好捕捉气流分离现象的S-A湍流模型;再次，针对旋翼后行桨叶动态失速时桨叶剖面来流速度较低、迎角较大的特点，为解决低来流速度时L-B半经验模型在旋翼翼型非定常动态失速计算中的局限性，并克服可压缩方程对低速流场计算收敛困难和精度低的问题，建立了基于Pletcher-Chen低速预处理方法、FAS多重网格法和隐式LU-SGS方法相结合的高效数值方法。应用发展的方法，分别针对NACA0012、SC1095旋翼翼型静态和轻度、深度动态失速进行计算，精确捕捉了气动力迟滞效应以及翼型前缘脱体涡的产生、对流和脱落过程，验证了本文方法的有效性;最后，着重针对NACA0012动态失速状态，开展了振荡参数对旋翼翼型非定常动态失速特性影响的分析，研究结果表明翼型迎角平均值、振幅及减缩频率的变化均能引起迟滞效应的改变并使得气动力峰值发生有规律的前、后移现象等。

旋翼;翼型;动态失速;N-S方程;运动嵌套网格;参数分析

0 引言

旋翼工作在严重非对称、非定常的涡流场中，旋翼桨叶的挥舞、周期变距以及畸变尾迹(诱导)形成的非均匀入流，导致桨叶剖面在不同方位角处的迎角有很大差别。当旋翼桨盘载荷很高时，后行桨叶工作在较大的迎角状态，容易产生气流分离进而出现复杂的非定常动态失速现象。旋翼动态失速现象虽然能够增大升力峰值，但同时造成阻力和力矩的突增，且翼型气动力中心不再稳定，这对旋翼的振动特性有重要影响，从而严重制约着直升机气动性能和飞行速度的提高［1］，因此其准确预测有助于直升机飞行性能包线的划定、振动水平的分析和设计等。与此同时，翼型作为旋翼的基本组成元素，其非定常动态失速特性反映了旋翼动态失速的基本特性。然而，旋翼及翼型的动态失速现象的机理仍在探索之中，对脱体涡的运动引起的气动力迟滞效应、逆压梯度的产生及大小等规律都未能彻底掌握，关于旋翼翼型非定常动态失速特性的研究一直是直升机空气动力学领域的研究焦点和难点。因此，进行旋翼翼型动态失速的研究具有重要的理论和实际应用价值。

目前，旋翼翼型动态失速特性的研究主要采用数值方法和试验方法。试验研究主要是针对风洞中振荡的翼型进行气动力、压力的测量［2-3］，因而可以直观地认知旋翼翼型动态失速中气动力的迟滞效应。但是旋翼翼型动态失速的风洞试验复杂、代价大、周期长并受试验测量技术的限制，并且很难捕捉动态失速过程中涡的流动细节、气流分离及再附等现象，因而其只能进行原理性或特定工况的研究。数值计算方法可以克服试验研究中存在的这些不足和难题，针对旋翼翼型动态失速可进行多种工况及流动细节的模拟，已成为旋翼翼型非定常动态失速特性研究的一条重要途径。

旋翼翼型动态失速模拟的数值方法包括基于试验数据的半经验模型和计算流体力学(CFD)方法。国外学者已发展了一系列的半经验动态失速理论模型，如Leishman-Beddoes(L-B)［4］、ONERA［5］、Johnson［6］等模型。其中，最为著名的是L-B二维翼型动态失速模型，该模型是由几个从翼型试验数据中获得经验参数，采用指数法建立的半经验模型。L-B模型在达到一定精度的要求下，大大缩短了计算时间，在直升机旋翼气动载荷、气动弹性以及飞行力学等方面有着广泛的应用［7-9］。但L-B等半经验模型并不是严格的预测方法［10］，仅对特定翼型适用，同时L-B模型对动态失速中气流再附时气动力的计算值与试验值相比有很大偏差，并且针对如旋翼后行桨叶的低来流速度、大迎角状态，其无法完全胜任翼型非定常气动力的准确模拟;此外，L-B模型只能针对气动力进行有效计算，无法模拟翼型非定常流场中复杂的气流分离及再附现象中的涡流动细节，从而很难更深入地探索动态失速的形成机理。

Leishman指出，只有通过基于Navier-Stokes方程的CFD方法才能完整地描述旋翼翼型的动态失速现象［11］。Jameson［12］较早地提出了非定常流动模拟的双时间方法，并采用Eluer方程对翼型的动态振荡过程进行数值模拟，有效地计算了翼型浅度失速时升力系数的迟滞效应，显示了CFD方法在模拟翼型非定常动态失速方面的优越性。由此，国内外学者相继开展了翼型动态失速现象的CFD数值分析研究［13-17］。现有CFD研究结果表明，动态失速的显著特征是具有能量的翼型前缘干扰涡的运动，涡运动引起翼型压力场、力和力矩等的瞬时变化，导致气动力表现出明显的非线性迟滞效应。然而，这些研究大多针对固定翼飞机的跨音速、小迎角的工作状态，对于旋翼后行桨叶翼型所处的低马赫数、大迎角状态的数值研究不多，并且针对翼型振荡参数对失速特性的影响分析也很少涉及，而上述这些对旋翼后行桨叶动态失速特性的研究却至关重要。

虽然通过CFD方法已经能够有效对旋翼翼型动态失速中气动力特性进行预测，但旋翼翼型非定常动态失速特性分析的CFD方法的发展及应用还远未成熟，对动态失速现象的模拟结果与试验值仍存在较大差异［12］。这主要是因为:(1)旋翼翼型非定常涡流场中存在较大的涡量非物理耗散，主要来源于数值方法的截断误差;(2)目前针对翼型动态失速的研究多采用变形网格方法展开，虽然通过弹簧方法等网格变形方法能够有效避免畸形网格的出现，但仍难保证网格正交性和光滑性，另外，变形网格技术需要满足几何守恒定律，这均影响计算精度。

针对这些问题，本文采用基于改进最小距离法和Hole Map方法生成绕旋翼翼型的运动嵌套网格;将Roe-MUSCL格式、隐式LU-SGS方法相结合，以基于S-A湍流模型的RANS方程模拟大范围气流分离现象以及翼型前缘脱体涡的流动过程;针对旋翼后行桨叶来流速度低、L-B模型在低速状态的计算失真等问题，并扩大可压缩N-S方程的应用范围，建立了FAS聚合式多重网格和Pletcher-Chen预处理方法结合的高效CFD方法。最后着重开展了旋翼翼型动态失速特性及深度失速状态下振荡参数影响分析，获得了一些有意义的结论。

图1 绕翼型的C型网格Fig.1C-type grid around airfoil

1 网格生成方法

1.1旋翼翼型贴体网格生成

旋翼翼型贴体网格是进行动态失速分析的基础，网格的正交性、光滑性以及翼型表面附近的网格布置直接影响数值模拟的精度。为生成壁面网格尺度、正交性满足湍流求解要求的光滑过渡的翼型网格，本文发展了基于Poisson方程求解的翼型网格生成方法。计算平面的Poisson控制方程为:

式中，α、β、γ为坐标变换系数，ψp、ψq为控制源项，其分别控制翼型网格的正交性和法向网格密度。

图1为生成的NACA0012翼型的整体C型网格示意图，其中翼型表面布置239个网格点，翼型法向60层网格，并设置第一层网格距离壁面1×10-6c(c为翼型弦长)。由图可以看出，网格正交性和贴体性良好，从而为N-S方程准确捕捉粘性及动态失速过程中气流分离提供基础。

图2 基于改进的MDM方法的贡献单元搜索过程Fig.2Procedure of donor cell searching by improved MDM method

1.2 运动嵌套网格生成方法

目前国内外针对翼型动态失速的模拟多采用变形网格方法，并通过几何守恒律计入网格变形对流动通量的影响。但当翼型运动幅度较大时，变形网格方法仍可能导致网格的畸变，从而影响计算精度。

基于此，本文采用运动嵌套网格方法进行旋翼翼型动态失速的非定常流场的数值模拟。首先，分别生成随翼型运动的贴体网格和固定不动笛卡尔背景网格;然后，采用Hole Map方法确定背景网格在翼型网格上的洞边界，在此基础上，提出了最小距离法(Minimum Distance Method，MDM)的改进策略，并进行背景网格人工内边界的贡献单元搜索，从而生成翼型与背景网格的运动嵌套网格。

MDM方法有一个不足，即搜索得到的与背景网格点P最近的翼型网格点可能是局部最近点。针对这一问题，文献［18］在翼型上下表面各选取一个网格点作为搜索起点，但这同时增加了贡献单元搜索的复杂性。本文提出了一种改进的MDM方法，采用该方法对背景网格上的点P进行贡献单元的搜索过程如图2所示，具体过程如下:

(1)初始化所有翼型网格点标记为0，并选取翼型网格点S(II，JJ)为起始点，其中II和JJ分别为沿翼型弦向和法向的网格点编号;

(2)以S点为中心，判断点S的四个相邻网格点到P的距离与|SP|的大小关系。跳过标记为1的点，若存在标记为0的网格点距离P点最近且小于|SP|，则令该点为新的S点，并将判断过的网格点标记为1，否则，将S点记为O，转到第(4)步;

(3)重复步骤(2)，直到找到O点;

(4)找出与O点关于尾迹奇异面对称的网格点，若|OP|小于|P|，则进入第(5)步;反之则将点记为S，释放O点，并返回第(2)步;

(5)通过矢量法判断点P在哪一包含O点的网格单元中，若不在这些单元中，则拓展网格范围判断，直到找到点P的贡献单元。

MDM方法搜索点是在一条曲线上进行的，类似于一维搜索，搜索次数为O［4×(IMAX+JMAX)］，可以有效地节省贡献单元搜索的时间。

翼型网格边界单元的贡献单元由两个一维搜索得到。嵌套网格之间的信息传递通过双线性插值完成。为减小插值过程中引入的数值误差，在挖洞时将背景网格洞边界远离流场梯度较大的翼型表面。图3给出了翼型网格运动、网格嵌套关系及贡献单元搜索结果示意图。

图3 围绕振荡翼型的运动嵌套网格边界单元及贡献单元示意图Fig.3Schematic diagram of boundary and donor cells of moving-embedded grids around oscillated airfoil

2 旋翼翼型非定常流场数值计算方法

2.1 流场控制方程求解

综合考虑旋翼翼型动态失速的非定常流场的复杂性，在进行动态翼型气动特性分析时采用雷诺平均N-S方程。基于面单元为dS的控制体Ω，积分形式的控制方程为:

其中，W、Fv分别为守恒变量和粘性通量。

基于翼型网格的运动，修正无粘通量为:

式中，Vt=nx为控制体表面的逆变速度，nx、ny表示控制面单位外法矢分量，则对流通量为:

其中u、v分别为x和y方向速度。

本文采用格心有限体积法在空间方向上进行离散。无粘通量采用Roe-MUSCL格式进行计算:

其中，A-Roei+1/2为Roe格式算子，i+1/2为左右单元L、R的交接面，WL、WR分别由MUSCL格式插值获得:

式中，Δ-、Δ+分别为向后、向前差分算子。κ为自由参数，本文取κ=1/3，即三阶迎风格式。

为模拟旋翼翼型流场的非定常特性，采用双时间方法进行时间推进。时间导数采用向后二阶差分进行离散，方程(2)可隐式离散为:

为准确模拟旋翼翼型动态失速流场的非定常特性，要求在进行时间推进时取很小的时间步长，这会降低数值模拟的效率。针对这一问题，本文采用隐式LU-SGS格式进行时间推进，从而有效加大时间步长，提高计算效率。时间推进过程为:

其中，L、D、U分别为下三角、对角和上三角矩阵。

2.2 S-A湍流模型

本文N-S方程湍流粘性系数计算采用S-A湍流模型。相对于代数湍流模型(B-L模型)，S-A模型能够更好地捕捉翼型振荡引起的大范围气流分离现象以及翼型前缘脱体涡运动过程，因而更适合于旋翼翼型动态失速过程的分析。S-A模型的控制方程为准运动涡粘性系数珓ν的输运方程，其积分形式为:

式中，Fc、Fν、QT分别为粘性系数输运方程的对流通量、粘性通量和单元源项，表达式如下:

对流项和粘性通量分别采用逆风、中心差分格式离散，时间推进方法与N-S方程一致。由珓ν最终求得运动涡粘性系数νt=珓νfv1，fv1为珓ν的函数。

2.3 多重网格方法

多重网格方法(Multigrid method)［19］能够将细网格上的低频误差转化为粗网格上的高频误差，从而达到在各自网格上消除其高频误差的目的，其工作量小、计算效率高，尤其适合旋翼及翼型复杂非定常流场的加速计算。本文采用FAS格式两重网格方法，首先，将细网格上求解得到流场守恒变量W+h和残值传递到粗网格上:

将修正量δW2h=Wn2h+1－W(2h0)叠加到细网格守恒变量中，得到新的守恒变量值:

为简化计算，本文在计算时只在网格数量大的翼型贴体网格上采用多重网格方法。

2.4 低速预处理方法

旋翼后行桨叶一般工作在低速状态，本文采用可压缩方程进行低来流速度状态翼型非定常动态失速计算时，会遇到“刚性”问题。这是因为当地流动速度和音速相差较大时，时间相关方程对流项特征矩阵的最大和最小特征值的比值也会很大，使计算收敛很慢。低速预处理方法能够通过改变方程的特征值来解决这一问题。预处理后的N-S方程为:

其中，WT=(p，u，v，T)T，Γ为预处理矩阵，本文采用Pletcher-chen［20］预处理方法。预处理后的N-S方程进行离散得到的残值为:

3 旋翼翼型非定常动态失速模拟

3.1 定常状态翼型流场的计算

为验证本文低速预处理、多重网格方法在流场模拟中的有效性，以NACA0012翼型为例，首先对可压状态Ma=0.55，α=8.34的流场进行了模拟，计算及试验［21］结果见图4。由图可知，本文预处理方法能够针对可压状态进行有效计算，并且结合多重网格方法能够更有效地提高流场求解的效率。

然后，考虑到旋翼后行桨叶的低速状态(来流马赫数一般小于0.35)，针对低来流马赫数Ma=0.1、α=3.59°的状态进行了计算。图5给出了收敛历程与预处理前后的翼型压强系数分布。由图可以看出，采用预处理方法能够消除可压方程在低速状态计算的压力分布的振荡，表明预处理方法能够有效提高求解精度;同时低速预处理、多重网格方法也能明显增强流场计算的收敛性。

图4 可压状态翼型流场收敛历程及压强系数分布Fig.4Convergence history and distribution of pressure coefficient for airfoil at compressible state

图5 低速状态翼型流场收敛历程与压强系数分布Fig.5Convergence history and pressure coefficient distribution of airfoil under low-speed condition

为验证本文方法对定常状态失速特性的计算能力，对NACA0012翼型在Ma=0.3，Re=4×106状态的静态失速特性进行了数值模拟，图6给出了翼型升力、力矩系数随迎角的变化。由图可知，本文计算获得的翼型升力、力矩系数与试验值［22］吻合较好，并且能够有效捕捉翼型的失速迎角，表明本文方法在翼型定常状态气动特性模拟方面的有效性。

为验证S-A模型对粘性计算的模拟精度，本文对RAE2822翼型的粘性流场进行了计算，计算状态为Ma=0.75，α=3.19°，Re=6.2×106，网格尺寸为369× 65，法向第一层网格距离壁面1×10-6c，以满足粘性计算的需要。图7给出了S-A模型和B-L模型计算的翼型上表面摩擦力系数分布及0.9c处速度型与试验值［21］的对比，从中可以看出，在捕捉分离区的能力上，S-A模型比B-L模型具有明显优势，也表明本文计算方法的有效性。

图6 升力、力矩系数随迎角变化曲线Fig.6Lift and momentum coefficients vary with AOA

图7RAE2822翼型上表面摩擦力系数分布及0.9c处速度型Fig.7Friction coefficient and velocity profile at 0.9c on the upper surface of RAE2822 airfoil

3.2 轻度失速算例分析

对于二维流动，翼型绕l/4弦长位置作正弦振荡，迎角变化规律为:

其中，α0、αm分别为迎角的平均值和振幅，无量纲时间t*=tU∞/c，缩减频率k=πfc/U∞。U∞为来流速度，c为翼型弦长。

为验证本文方法的有效性，选取NACA0012翼型的非定常流场进行求解。计算状态为典型前行桨叶工作状态:Ma=0.6，Re=4.8×106，其迎角变化规律为α=3.16°+4.59°sin(2×0.0811×t*)。

为准确地模拟每个时间步翼型的非定常气动力，将翼型一个振荡周期分为288个时间步，子迭代为200步，并进行了5个周期的计算。图8分别给出了采用B-L、S-A湍流模型计算收敛后的翼型升力、力矩系数的迟滞回线，由图可以看出，S-A模型的升力、力矩系数的计算值较B-L模型计算值与试验值［2］吻合更好，鉴于此，本文算例均采用基于S-A湍流模型的非定常N-S方程流场数值计算方法对振荡翼型的气动力变化进行有效模拟。

图9给出了S-A模型计算的翼型在一个周期内不同迎角的表面压强系数分布，图示箭头向上表示翼型的上仰过程，箭头向下表示翼型的下俯过程。由图可以看出，本文计算值与试验值吻合良好。

图8 翼型升力、力矩系数迟滞回线计算值与试验值的对比Fig.8Comparisons of lift and momentum coefficients of airfoil between calculated and experimental data

图9 不同迎角时翼型表面压强系数分布计算值与试验值的对比Fig.9Comparisons of pressure coefficient distributions of airfoils at different AOA

3.3 深度失速算例分析

为探索翼型深度动态失速时流场的分布特征，本文针对旋翼后行桨叶典型工况:Ma=0.2、α0=15°、αm=10°、k=0.05的状态进行了数值计算。图10为计算得到的翼型升力迟滞回线与L-B模型计算值及试验值的对比。从图中可以看出，马赫数低于0.3时，L-B模型计算得到的翼型升力系数与试验值差距很大，而本文在升力系数峰值、迟滞回线面积以及气流分离的再附的计算值与试验值［21］更接近，表明本文建立的方法能够较好地模拟较低来流速度下的翼型非定常动态失速特性。

图11给出了翼型半个周期的振荡过程中翼型表面涡量分布云图。从图中可以看出，在翼型上仰过程中，随着翼型迎角的增大，翼型前缘涡逐渐加剧，并且有后缘涡产生。当翼型迎角超过动态失速临界值时，前缘产生脱体涡，并进一步导致翼型前缘的气流分离。当翼型从最大迎角位置下俯时，由于相对来流速度的增大，前缘脱体涡以及气流分离现象逐渐加剧，当前缘涡流过翼型后缘并进入尾迹区时，翼型上表面气流完全分离，同时伴随着升力的突然减小。当翼型迎角恢复到足够小时，气流会在前缘再附。前缘涡的生成、沿翼型表面的对流直至脱落导致上仰和下俯过程翼型在迎角相同时的升力不对称，从而引起翼型升力的迟滞效应。

图10 翼型升力系数迟滞回线Fig.10Lift coefficient hysteresis loop of airfoil

图11 翼型上表面半个周期内不同迎角时的涡量分布云图Fig.11Vorticity contours on the upper surface of airfoil at different AOA in half a period

3.4SC1095旋翼翼型深度失速特性分析

开展了SC1095先进旋翼翼型深度动态失速时气动特性的计算。该翼型为美国先进的通用直升机黑鹰旋翼采用的主要翼型。针对后行桨叶典型的工况Ma=0.302，α0=9.78°，αm=9.9°，k=0.099的进行了数值计算。图12给出了计算得到的翼型升力迟滞回线的本文计算值、L-B模型计算值及试验值的对比。从图中可以看出，相对于L-B模型，本文计算值与试验值［3］吻合更好，表明旋翼翼型非定常流场模拟方法能够针对不同旋翼翼型在低马赫数状态的动态失速特性进行有效计算。

图12SC1095翼型升力系数迟滞回线计算值与L-B模型及试验值的对比Fig.12Comparisons of lift coefficient hysteresis loop of SC1095 airfoil among calculated，L-B model and experimental data

4 旋翼翼型运动参数影响分析

由旋翼翼型动态失速的运动方程(15)可知，翼型动态失速特性与俯仰运动的平均迎角、振幅及缩减频率均有关。本文通过对NACA0012翼型在不同参数组合下的气动力、力矩特性进行数值模拟，开展翼型动态失速特性的参数分析。

基准状态为Ma=0.3，α0=10°，αm=10°，k=0.1，进行参数分析时只改变所分析参数，其余参数与基准状态一致。

4.1 振幅的影响

针对旋翼后行桨叶的大迎角的动态失速特性，对三个不同振幅αm=5°、10°、15°分别进行了数值模拟。图13给出了在不同振幅时翼型升力、力矩系数随迎角的变化规律，并给出了定常状态升力系数的计算值。从图中可以看出，随着振幅的增大，迟滞效应增强。在翼型迎角增大时，不同振幅对应的升力系数在相同迎角时计算结果基本相同，并且发生失速时的迎角相对于定常状态均有所推迟。这是因为当迎角大于静态分离迎角时，分离点推迟及前缘脱体涡的影响使得升力不发生衰减，直到前缘涡从翼型后缘脱离。随着前缘涡的脱离，由于翼型上表面气流分离及涡诱导速度的损失，翼型升力系数会急剧衰减。翼型迎角减小时，同样存在分离气流再附着的延迟。随着振幅的增大，气流再附延迟更明显，即升力迟滞回线的面积更大。与此同时，随着翼型前缘涡的运动、脱落及气流分离，压力中心后移，低头力矩有一个陡增，并且振幅的增大使得低头力矩峰值在增大的同时发生后移;当迎角减小时，压力中心向翼型前缘移动，并且当迎角减小到一定值时，会产生一个抬头力矩峰值，随着振幅的增大，迟滞效应加剧，使得抬头力矩的峰值发生前移。

图13 不同迎角振幅时翼型气动力系数迟滞回线Fig.13Aerodynamic coefficient hysteresis loop of airfoil at different amplitudes

4.2 缩减频率的影响

直升机旋翼工作状态的缩减频率在0.05～0.15之间，因此，本文针对k=0.05、0.1、0.15三个状态分别进行了数值模拟。图14给出了在不同缩减频率时翼型升力、力矩系数随迎角的变化规律。由图可以看出，随着缩减频率的增大，气动力系数迟滞效应更加强烈，迟滞回线的面积增大，升力系数峰值增大，并且失速发生位置后移。与此同时，低头力矩的峰值随着缩减频率的增大而增大并且后移，抬头力矩峰值随缩减频率的增大发生前移。

缩减频率表征了旋翼翼型动态失速特性的强弱程度，频率越大，动态特性越强烈，反之，动态特性越弱。随着缩减频率的增加，一方面，在翼型迎角增大过程中，分离点的滞后作用更加显著，故动态情况下分离点的位置向后移动，即分离迎角增大;另一方面，在迎角减小的过程中，分离点再附着的延迟也更显著，故升力系数的减小量也更大。

图14 不同缩减频率时翼型气动力系数迟滞回线Fig.14Aerodynamic coefficient hysteresis loop of airfoil at different reduced frequencies

图15 不同平均迎角时翼型气动力系数迟滞回线Fig.15Aerodynamic coefficient hysteresis loop of airfoil at different average AOA

4.3 平均迎角的影响

针对翼型三个不同平均迎角α0=5°、10°、15°分别进行了数值模拟，图15给出了在不同平均迎角时翼型升力、力矩系数随迎角的变化规律。从图中可见，随着平均迎角的增大，升力系数峰值增大并且后移;此外，随着平均迎角的增大，气动力的迟滞效应明显增强，并且气流分离后的再附现象随着平均迎角的增大有逐渐消失的趋势。

平均迎角同迎角振幅一起决定了翼型动态过程的幅值范围，即动态特性的程度。最大迎角(平均迎角+正幅值)越大，动态特性越明显。当最大迎角略大于失速迎角时，翼型处于小分离状态，且前缘涡的强度也不是很大，故升力系数的变化不是很剧烈，只形成一个小的迟滞环(如图中平均迎角为5°的状态)，即轻度失速。但当最大迎角远大于静态失速迎角时，则翼型发生完全分离，前缘涡也得到了显著的加强，故迟滞环的面积更大(如图中平均迎角为10°以上的状态)，即深度失速。

5 结论

本文建立了基于运动嵌套网格技术的旋翼翼型动态失速数值模拟方法，并以NACA0012和SC1095旋翼翼型为研究对象，进行了翼型动态失速数值计算验证，并与试验值以及L-B模型计算结果进行了对比。在此基础上，着重开展了NACA0012旋翼翼型动态失速特性的参数分析。综合本文的计算分析结果，主要得出以下几点结论:

(1)本文建立的CFD方法能够精确模拟翼型定常、非定常状态的气动特性，适合旋翼后行桨叶低速、大迎角工作状态的模拟分析;

(2)采用多重网格法和低速预处理法不仅能够提高流场求解的效率，在低速状态更能够提高流场的求解精度;并且运用该方法在较低来流速度下能够取得比L-B模型更精确的翼型升力和力矩特性;

(3)翼型在上仰过程中产生的前缘脱体涡和气流分离会在翼型迎角减小时有所加剧，当迎角减小到一定值时产生的脱体涡和气流分离现象消失，气流发生再附，并且再附时翼型迎角比分离时有所滞后，进而引起翼型气动力的迟滞效应。

(4)一方面，随着翼型迎角振幅的增大，气流分离和再附延迟现象更为明显，因此气动力迟滞效应更加明显，迟滞回线面积增大;另一方面，随着振幅的增大，升力系数、低头力矩的峰值也会增大并发生后移，抬头力矩的峰值前移。

(5)随着缩减频率的增加，在翼型迎角增大过程中，分离点的滞后作用更加显著，故动态情况下分离点的位置向后移动，即分离迎角增大。相同的，在迎角减小的过程中，分离点再附着的延迟也更显著，故升力系数的减小量更大。

(6)最大迎角(平均迎角+正幅值)越大，动态特性越明显。轻度失速时，翼型处于小分离状态，前缘涡的强度不大，因此，升力系数的变化平缓，只形成一个小的迟滞环;深度失速时，翼型发生完全分离，前缘涡强度很大，进而产生大面积迟滞环。

［1］Geissler W，Raffel M，Dietz G，et al.Helicopter aerodynamics with emphasis placed on dynamic stall springer［M］.Berlin Heidelberg，2007.

［2］Landon R H.Compendium of unsteady aerodynamic measurements［R］.AGARD-R-702，1982.

［3］McAlister K W，Pucci S L，McCroskey W J，et al.Experimental study of dynamic stall on advanced airfoil sections［R］.NASA Technical Memorandum 84245，1982，2:330-380.

［4］Leishman J G.Validation of approximate indicial aerodynamic functions for two-dimensional subsonic flow［J］.Journal of Aircraft，1988，25(10):914-922.

［5］McAlister K W，Lambert O，Pitot D.Application of the ONERA model of dynamic stall［R］.NASA TP-2399，1984.

［6］Dindar M，Kaynak U.Effect of turbulence modeling on dynamic stall of a NACA0012 airfoil［R］.AIAA 1992-0027.

［7］Leishman J G，Beddoes T S.A semi-empirical model for dynamic stall［J］.Journal of the American Helicopter Society，1989，34:3-17.

［8］Sheng M，Galbraith R A McD，Coton F N.A modified dynamic stall model for low mach numbers［C］//45thAIAA Aerospace Sciences Meeting and Exhibit，Reno，Nevada，2007:1-17.

［9］Song C Y，Xu G H.Computations of unsteady dynamic stall responses on rotor airfoils［J］.Acta Aerodynamica Sinica，2007，25(4): 461-467.(in Chinese)宋辰瑶，徐国华.旋翼翼型非定常动态失速响应的计算［J］.空气动力学学报，2007，25(4):461-467.

［10］邵明玉，杨茂，陈凤明.旋翼翼型动态失速特性的数值仿真研究［J］.计算机仿真，2012，29(7):70-74.Shao M Y，Yang M，Chen F M.Numerical simulation study on dynamic stall characteristics of rotor airfoil［J］.Computer Simulation，2012，29(7):70-74.

［11］Leishman J G.Principle of helicopter aerodynamics［M］.2007，525-561.

［12］Jameson A.Time dependent calculation using multigrid with application to unsteady flows past airfoils and wings［R］.AIAA 1991-1596.

［13］Ko S，McCroskey W J.Computations of unsteady separating flows over an oscillating airfoil［J］.AIAA Journal，1997，35:1235-1238.

［14］Lee T，Gerontakos P.Investigation of flow over an oscillating airfoil［J］.J.Fluid Mech.，2004，512:313-341.

［15］Qian W Q，Leung R C K.Numerical simulation ofairfoil dynamic stall incorporating transition model［J］.Acta Aerodynamica Sinica，2008，26(1):50-55.(in Chinese)钱炜祺，Leung R C K.考虑转捩影响的翼型动态失速数值模拟［J］.空气动力学学报，2008，26(1):50-55.

［16］Dumlupinar E，Murthy V R.Investigation of dynamic stall of airfoils and wings by CFD［R］.AIAA 2011-3511，2011.

［17］Wernert B E，Schreck P，Raffel S.Investigation of threedimensional dynamic stall using computational fluid dynamics［J］.AIAA Journal，2005，43:1023-1033.

［18］Yang W Q，Song B F，Song W P.Distance decreasing method forconfirming corresponding cells of overset grids and its application［J］.Acta Aeronautica et Astronautica Sinica，2009，30(2):205-212.(in Chinese)

杨文青，宋笔锋，宋文萍.高效确定重叠网格对应关系的距离缩减法及其应用［J］.航空学报，2009，30(2):205-212.

［19］Mavriplis D J.Multigrid solution of the two dimensional Euler equations on unstructured triangular meshes［C］//AIAA 25th Aerospace Sciences Meeting，Reno，Nevada，1987:1-16.

［20］Pletcher R H，Chen K H.On solving the compressible Navier-Stokes equations for unsteady flows at very low Mach numbers［R］.AIAA-93-3368:765-775.

［21］Terry L H.Viscous transonic airfoil workshop compendium of results［R］.AIAA-87-1460，1987.

［22］Lee T，Gerontakos P.Investigation of flow over an oscillating airfoil［J］.Journal of fluid Mechanics，2004，512:313-341.

勘误表

Zhao Guoqing，Zhao Qijun，Wang Qing
(National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics，Nanjing University of Aeronautics and Astronautics，Nanjing210016，China)

A high-efficiency and high-precision CFD method for simulating the unsteady dynamic stall of rotor airfoil has been established based on moving-embedded grid and compressible RANS equations.Firstly，the generation method of viscous and orthogonal body-fitted grid around the rotor airfoil is developed by solving Poisson equations.Meanwhile，aiming at overcoming the shortcoming of spring simulation approach which may result in the distortion of grid，an improved Minimum Distance Method is proposed to generate the embedded grid around airfoil.Secondly，in order to simulate the hysteresis effect of aerodynamic forces caused by the turbulence separation and re-attachment of the flow，a high-precision method on the analysis of unsteady aerodynamic characteristics of rotor airfoil is developed by employing RANS equations and dual-time method.The S-A turbulence model is employed to capture the separation phenomenon of flow around airfoil.Thirdly，according to the conditions of low-speed inflow and high AOAs ofthe retreating blade，together with the limitation of L-B semi-empirical model on the calculation of unsteady dynamic stall of airfoil，a combination method of Pletcher-Chen preconditioning，FAS multigrid approaches and implicit LU-SGS scheme is established to overcome the problems of convergence difficulty and insufficient precision of compressible equations.The steady，mild and deep dynamic stall cases of NACA0012 and SC1095 rotor airfoils are calculated using this previously mentioned method，the hysteresis effect and the formation，convection，shedding of the vortical disturbance are well captured，the effectiveness of numerical simulation method on dynamic stall is verified.Finally，focus on the deep stall of NACA0012 airfoil，the influence analyses of parameters on the unsteady aerodynamic forces of rotor airfoil are carried out，and the results demonstrate that the exchanges of averaged AOA，amplitude and reduced frequency may cause a variational hysteresis effect and regularly changes of peak value of aerodynamic force.

rotor;airfoil;dynamic stall;N-S equations;moving-embedded grid;parameters analyses

序号年、卷、期、页码更正前更正后1《空气动力学学报》2014年32卷第6期第740页中文摘要第三行在此研究中，翼型做沉标和俯仰运动在此研究中，翼型做沉浮和俯仰运动2《空气动力学学报》2014年32卷第6期编辑部2014年纪要第五行……审稿和编辑出版工作，从2014第4期起连续刊载……审稿和编辑出版工作，从2014第5期起连续刊载3《空气动力学学报》2014年32卷第6期编辑部2014年纪要第十六行目前已经准备就绪并将在2014年实施的项目有目前已经准备就绪并将在2015年实施的项目有

V211.52;V211.3

Adoi:10.7638/kqdlxxb-2013.0010

Simulations and parametric analyses on unsteady dynamic stall characteristics of rotor airfoil based on CFD method

0258-1825(2015)01-0072-10

2013-01-23;

2013-02-25

国家自然科学基金项目(11272150);江苏高校优势学科建设工程资助项目

赵国庆(1984-)，男，山东五莲人，博士研究生，研究方向:旋翼流动控制、旋翼计算流体力学、直升机空气动力学.E-mail: zgq198495@nuaa.edu.cn

招启军(1977-)，男，教授、博士生导师，研究方向:直升机计算流体力学、直升机空气动力学及流动控制、气动噪声、总体设计等.E-mail:zhaoqijun@nuaa.edu.cn

赵国庆，招启军，王清.旋翼翼型非定常动态失速特性的CFD模拟及参数分析［J］.空气动力学学报，2015，33(1):72-81.

10.7638/kqdlxxb-2013.0010.Zhao G Q，Zhao Q J，Wang Q.Simulations and parametric analyses on unsteady dynamic stall characteristics of rotor airfoil based on CFD method［J］.Acta Aerodynamica Sinica，2015，33(1):72-81.