直升机的飞行性能好坏在很大程度上取决于旋翼的气动特性，而旋翼的气动特性又与旋翼翼型密切相关。相对于固定翼飞行器的机翼，直升机旋翼通常工作在严重的非定常气动环境中，特别是在前飞情况下，旋翼翼型一直处于动态失速状态。与定常状态下翼型的气动特性不同，动态失速状态下的翼型气动特性呈现一个明显的迟滞回线，由此带来的翼型气动力的变化将给直升机旋翼气动特性带来很多不利影响，例如失速颤振、振动载荷激增、噪声增强^[1-2]等，因此旋翼翼型的动态失速特性一直是直升机非定常空气动力学研究领域的难点和重点，开展旋翼翼型的动态失速特性的研究对于认识和改造旋翼气动特性有重要的实际意义和学术价值。

国外对于翼型动态失速特性的研究开展得较早，许多机构在20世纪中期就利用热线、烟流流场显示等手段对翼型动态失速特性做了细致的试验测量研究。其中，McCroskey^[3]对翼型的动态失速做了测力、测压的试验研究，分别对比了不同缩减频率、迎角和马赫数影响下的翼型升力系数和力矩系数，并给出了动态失速下不同时刻的压力系数分布。McAlister等^[4]对VR-7及NACA0012等翼型的动态失速特性进行了试验研究，给出了不同状态下的翼型气动力参数。Carr^[5]在1988年进一步从理论上描述了翼型动态失速发生、演化的过程，并分析对比了翼型气动参数，如速度、迎角和缩减频率等对动态失速特性的影响，并分析总结了三维效应对翼型动态失速特性的影响。VanDyken等^[6]研究了压缩性对二维失速影响的试验,马赫数直到0.45，雷诺数直到400万。Martin等^[7]在阿姆斯研究中心流体力学实验室25 cm×35 cm可压缩动态失速实验设施上进行了VR-12翼型试验，马赫数为0.2～0.45，缩减频率为0～0.1。近年来，Rbee^[8]对翼型前缘分离涡的特征和影响做了较系统的研究，测量了前缘分离涡对翼型压力系数的影响以及动态失速下翼型气动力参数的特征。中国国内翼型动态失速的研究相对于国外起步较晚，主要有孙茂等^[9]试验研究了翼型在等速上仰过程中的气动特性，此外，一些研究人员采用CFD方法对翼型在相应动态失速状态进行了数值模拟研究^[10]。

目前，关于旋翼翼型外形对动态失速特性影响的研究不多，因此有必要进一步研究和分析这一内在关系，以便于更合理便捷地设计出新的直升机旋翼翼型。试验及数值模拟表明^[11]，翼型前缘外形对翼型的动态失速特性的影响十分重要，不同的前缘外形能够显著地改变翼型在动态失速过程中的压力分布，进而影响到翼型的气流分离及前缘涡的生成。因此，本文重点研究和分析了翼型前缘外形对翼型动态失速特性的影响。基于NACA0012翼型设计了3类(每类2种，修改翼型1~6) 修改翼型，其中每种修改翼型各有两个不同变形量的前缘外形，并通过对比分析了不同变形翼型与NACA0012翼型的动态失速特性。

1 旋翼翼型流场求解方法 1.1 网格生成方法

围绕旋翼翼型的网格生成采用基于Poisson方程^[12]的椭圆方程网格生成方法，相应的控制方程为

$\left\{ \begin{align} & {{\xi }_{xx}}+{{\xi }_{yy}}=P\left( \xi ,\eta \right)\text{ } \\ & {{\eta }_{xx}}+{{\eta }_{yy}}=Q\left( \xi ,\eta \right)~ \\ \end{align} \right.$

(1)

式中：P(ξ,η)和Q(ξ,η)为控制函数。通过改变P(ξ,η)可以控制网格的正交性，改变Q(ξ,η)可以控制网格的疏密程度。在翼型动态失速计算中，采用嵌套网格的方法来满足翼型的俯仰要求。嵌套网格如图 1所示，背景网格为笛卡尔网格，大小为119×119，翼型网格为C形网格，大小为379×60(弦向×展向)。

1.2 旋翼翼型流场求解方法

为捕捉旋翼流场中边界层流动特性，提高旋翼翼型优化设计的精度，本文的流场求解算法中采用积分形式的雷诺平均N-S方程作为旋翼流场求解控制方程

$\partial \partial t\iiint\limits_{\Omega }{Wd\Omega }+\iint\limits_{\partial \Omega }{({{F}_{c}}-{{F}_{v}})}dS=0$

(2)

式中:W为守恒变量;Ω为网格体积；F_c,F_v分别为对流、黏性通量，表达式如下

$W = \left[ \matrix{ \rho \hfill \cr \rho u \hfill \cr \rho v \hfill \cr \rho E \hfill \cr} \right]$

(3)

${F_c} = \left[ \matrix{ \rho {V_r} \hfill \cr \rho u{V_r} + {n_x}p \hfill \cr \rho v{V_r} + {n_y}p \hfill \cr \rho H{V_r} + {V_t}p \hfill \cr} \right]$

(4)

${F_v} = \left[ \matrix{ 0 \hfill \cr {n_x}{\tau _{xx}} + {n_y}{\tau _{xy}} \hfill \cr {n_x}{\tau _{yx}} + {n_y}{\tau _{yy}} \hfill \cr {n_x}{\Theta _x} + {n_y}{\Theta _y} \hfill \cr} \right]$

(5)

式中：V_r=(u,v)^T为流场的相对速度；ρ，E，p，n=(n_x,n_y)以及H分别为空气密度、总能、压强、单位法矢和总焓；τ_ij，Θ_i为黏性相关项。在流场求解中，通过翼型网格在背景网格坐标系进行俯仰运动，故在通量求解中引入网格单元的移动速度V_t，从而对流通量中的速度变化描述为V_r=V-V_t,V表示流场绝对速度。

为计算动态失速过程中翼型流场的非定常特性，本文采用双时间法来进行时间离散。物理时间采用具有二阶精度的二阶后向差分，因此式(5) 可以离散为

$\frac{3{{W}^{n+1}}-4{{W}^{n}}+{{W}^{n-1}}}{2\Delta t}+R({{W}^{n+1}})=0$

(6)

式中：Δt为物理时间步长，在此基础上叠加伪时间推进

$\frac{d({{\Omega }^{n+1}}{{W}^{n+1}})}{d\tau }+{{R}^{*}}({{W}^{n+1}})~=0$

(7)

${{R}^{*}}({{W}^{n+1}})=\frac{3{{W}^{n+1}}-4{{W}^{n}}+{{W}^{n-1}}}{2\Delta t}+R({{W}^{n+1}})~$

(8)

为了提高计算效率，本文虚拟时间推进采用隐式LU-SGS方法进行时间离散^[13]。同时，为了更好地捕捉翼型振荡引起的大范围气流分离现象以及翼型前缘分离涡的运动过程，本文N-S方程湍流黏性系数的计算采用S-A湍流模型^[14]。

为了验证本文建立的CFD方法的准确性，图 2给出了一个典型的NACA0012翼型动态失速模拟结果。模拟的状态为缩减频率k为0.1，迎角变化为 α=12.0°+8.5°sin(ωt) ，来流速度马赫数为0.3。从图中可以看出，本文的模拟结果与试验值基本上吻合，很好地模拟了气流在动态失速状态下的再附着过程。

为了进一步验证本文方法对失速状态下的模拟准确性，图 3给出了NACA0012翼型在深度动态失速状态下的CFD模拟结果与试验值对比。缩减频率k为0.151，迎角变化为 α=14.91°+9.88°sin(ωt) ，来流速度马赫数为0.283。整体而言，CFD模拟结果与试验值^[15]相比吻合得很好，说明本文采用的动态失速数值模拟方法能够很好地模拟不同翼型的二维动态失速特性。

2 算例分析

为了研究和分析翼型前缘外形对翼型动态失速气动特性的影响，分别基于NACA0012翼型设计了3类不同前缘外形的新翼型，每一类包括2种不同变形量值的翼型，其外形同NACA0012翼型的对比结果如图 4所示。从图中可以看出，第一类修改翼型修改了翼型前缘附近(0～0.4c，c表示弦长)上表面的外形，从而增加了翼型的前缘半径，而且修改翼型2具有更大的变形量，即意味着有较大的前缘半径。第二类修改翼型旨在修改翼型的前缘附近下表面的外形，同NACA0012翼型相比，修改翼型的的前缘半径更大。第三类修改翼型重点考虑了翼型前缘附近弯度的变化，修改翼型在保持最大厚度不变的状态下，增加了翼型的弯度，而且修改翼型6的弯度与修改翼型5相比更大一些。

2.1 第一类修改翼型动态失速分析

图 5,6给出了第一类修改翼型在来流速度马赫数为0.3，缩减频率为0.075，迎角变化为10°±8°动状态下的动态失速特性对比图。从图 5中可以看出，修改翼型前缘上表面外形可以显著影响到翼型的动态失速特性。修改翼型的动态失速特性同NACA0012翼型的动态失速特性相比，具有更好的再附着特性，即在翼型下俯过程中，升力系数恢复得更快，同时阻力系数及力矩系数的峰值也更小。同时可以看出，修改翼型的变形越大，对动态失速特性的抑制也越明显。图 6给出了NACA0012翼型及修改翼型的流线图及涡量云图。图中“”表示迎角增大，“”表示迎角减小。从图中可以看出，同NACA0012翼型相比，修改翼型1和修改翼型2的前缘分离涡更弱，同时，前缘涡脱落时的翼型迎角也更大，而且变形越大，前缘涡也越弱。

分析表明，翼型前缘附近上表面凸起能够增加翼型前缘上表面附近的曲率，从而有效地降低了翼型前缘附近的逆压梯度，削弱了动态失速过程中产生的前缘分离涡的强度，进而有效地改善翼型的非定常气动特性。从修改翼型1同修改翼型2的对比中可以看出，上表面凸起变形越大，对动态失速特性的抑制也越大。

2.2 第二类修改翼型动态失速分析

图 7，8给出了第二类修改翼型在来流速度马赫数为0.3，缩减频率为0.075，迎角变化为10°±8°动状态下的动态失速特性对比图。从图 7的气动力曲线中可以看出，第二类修改翼型的动态失速特性同NACA0012翼型的动态失速特性相比没有明显的变化，这说明翼型下表面的变化对翼型动态失速特性的影响不敏感。从图 8的流线图可以看出，不同翼型的前缘分离涡的特点也基本相似，由此可以得出翼型下表面外形的改变并没有影响到翼型上表面前缘附近的逆压梯度的变化，因此没能对前缘涡产生明显的影响。

2.3 第三类修改翼型动态失速分析

图 9，10给出了第二类修改翼型在马赫数为0.3，缩减频率为0.075，迎角变化为10°±8°动状态下的动态失速特性对比图。从图 9的升力系数C_l曲线中可以看出，修改翼型的最大升力系数有所降低，但在大迎角下的气流再附着特性有所改善。与此对应，阻力系数C_d及力矩系数C_m的峰值有所减小，而且修改翼型6的减小量比修改翼型5的减小量更大一些。从图 10的流线图中可以看出，修改翼型的前缘分离涡的强度更小一些，而且前缘涡脱落的迎角也更大一些。因此，通过修改翼型前缘附近弯度，可以有效抑制翼型的动态失速，且变形越大，抑制作用越明显。

3 结论

本文采用运动嵌套网格方法，构建了一套适合于翼型动态失速特性模拟的数值方法。通过改变NACA0012翼型前缘外形，分析和对比了3种类型翼型的动态失速特性，可以得出如下结论：

(1) 翼型前缘附近上表面外形对翼型的动态失速特性具有重要的影响，通过改变翼型上表面外形，增加翼型的前缘半径，能够有效地抑制翼型动态失速涡的生成，从而改善翼型的动态失速特性。

(2) 翼型前缘附近下表面外形对翼型的动态失速特性不是很敏感，通过修改翼型下表面来增加翼型前缘半径对动态失速特性的影响是有限的。

(3) 通过修改翼型前缘附近的中弧线曲率，可以影响到翼型的动态失速特性，增加翼型前缘附近的中弧线曲率，能够抑制翼型的动态失速特性。