低雷诺数翼型局部振动非定常气动特性

李冠雄，马东立，杨穆清，郭阳

北京航空航天大学 航空科学与工程学院，北京 100083

临近空间飞行器由于在民用及军事领域有广泛的应用前景而成为国内外研究的热点[1]。临近空间飞行器飞行高度高，飞行速度低，因此处于低雷诺数流动范围，机翼弦向雷诺数一般在106以下[2]。低雷诺数条件下飞行器气动特性会显著下降并出现一些特殊流动现象，机翼升阻比急剧下降，翼型在小迎角时出现升力系数非线性效应，在大迎角时出现滞回现象[3]。其次，为减轻结构重量该类飞行器多采用轻质结构甚至薄膜结构，在飞行过程中轻质蒙皮不可避免地会出现振动现象，蒙皮振动会对机翼气动特性及流场结构产生显著影响[4]。由于在某些情况下蒙皮振动可提高机翼气动性能，因此蒙皮振动逐渐成为机翼主动流动控制的一个研究方向[5]。

为改善低雷诺数机翼气动特性，可采用多种流动控制措施，目前国内外研究最多的主动流动控制措施为抽吸气流动控制[6]，抽吸气可有效抑制机翼表面层流分离提高机翼气动特性[7]，但抽吸气流动控制需消耗较大能量，抽吸系统较为复杂且机翼表面开孔可能降低机翼结构强度[8]。随着智能材料与大柔性结构的发展，局部蒙皮振动为改善低雷诺数飞行器气动特性提供了新的途径[9]。目前，国内外学者已对机翼蒙皮振动开展了一系列研究，康伟等研究了低雷诺数下翼面局部振动的增升机理，研究表明非定常流动分离中旋涡之间的距离及其演化频率与振动频率的关系是影响翼型翼面局部振动增升效果的重要因素[10]。Gordnier 在低雷诺数条件下研究了薄膜机翼振动对气动特性的影响，研究结果显示，在大迎角状态下，薄膜机翼表面动态位移对机翼非定常流场有显著的影响，可延迟机翼失速，提高大迎角升力系数[11]。Song和Breuer利用风洞试验对哺乳动物仿生薄膜翼面的气动特性进行了研究[12]；Rojratsirikul等采用粒子图像测速(PIV)技术对低雷诺数下二维薄膜翼型的非定常气动特性进行了研究[13]。实验结果证实，弹性结构与流场的耦合作用有助于翼型气动特性的改善。

目前专家学者对机翼振动的研究多集中在固定位置局部蒙皮振动或蒙皮整体振动对机翼气动特性的影响，所研究的气动特性多为机翼大迎角升力特性[14-15]。对于局部振动位置的影响及振动对巡航迎角附近机翼阻力特性影响的研究还不深入。本文在前人研究的基础上，建立具有不同弦向位置、不同振动频率和振幅的局部振动力学模型，研究振动位置、振动频率及振幅对巡航迎角附近机翼升阻特性和流场结构的影响规律，并对蒙皮局部振动的增升减阻机理进行深入研究，得到合理的蒙皮振动流动控制方式，有效提高机翼气动性能。

1 研究方法

1.1 数值计算方法

通过求解雷诺平均Navier-Stokes方程，采用非定常方法对流场进行数值计算，控制方程的积分形式为

(1)

式中：V为控制体体积；t为时间变量；Q为守恒变量向量；S为控制体表面积；F为通过表面S的净通量向量(包含黏性项和无黏项)；n为表面S的单位外法向向量。

采用有限体积方法求解式(1)所示的控制方程，无黏项离散采用二阶迎风格式，黏性项离散采用二阶中心差分格式，时间推进方式采用双时间步法[16]。

间歇因子γ的输运方程为

(2)

式中：ρ为流体密度；Uj为j方向速度；xj为j方向坐标；μ为动力黏性系数；μt为涡黏系数；σγ为常数项；转捩源项定义为

Eγ1=ce1Pγ1γ

再层流化源项定义为

{Pγ 2=ca2ρΩγFturb

Eγ 2=ce2Pγ 2γ

其中：Flength为经验修正系数，其控制了转捩区的长度；Fturb用于避免层流边界层外部或黏性底层内部再次形成层流，Fturb=e-(RT/4)4,RT为相对黏度；s为应变率的大小；Fonset用于触发湍流脉动的间歇性；Ω为涡强度；间歇因子方程中的ce1、ce2、ca1、ca2为常数项。

(3)

式中：uj为j方向时均速度；σθ t为扩散项系数；源项定义为

其中：cθ t为源项系数；U为当地速度；Reθ t为转捩临界雷诺数；δ为当地附面层厚度；Fwake为尾迹区开关函数。

Menter等指出，为准确地捕捉层流和转捩边界层，在低雷诺数流动计算中，需使y+接近于1[18]。本文采用结构化的H型二维网格，翼型弦长为c，计算域边界距翼型前缘为20c，距翼型后缘为25c，上下边界距翼型弦线为15c。整个计算域及壁面附近网格如图1所示。

为验证本文非定常CFD计算方法及网格对低雷诺数翼型升阻特性时均值计算的准确性，选取Eppler387翼型进行数值计算并与风洞试验结果进行对比。为研究网格参数对计算结果的影响，选取3种具有不同参数的网格分别计算，进行网格敏感性分析，3种网格参数如表1所示。

在雷诺数为3.0×105条件下，用不同参数的网格计算得到得翼型升力系数CL和阻力系数CD随迎角α变化的时均化计算结果及美国Langley研究中心的低湍流度增压风洞(LTPT)试验数据[19]对比如图2所示。

由图2可知网格A计算得到的CL在各迎角下均小于试验值，CD均高于试验值。网格B和网格C计算得到的CL基本重合且与试验值吻合较好，CD在迎角小于8°范围内与试验值吻合较好，由于本文主要研究翼型小迎角气动特性，因此网格B和网格C均能满足本文计算要求，综合考虑计算精度及计算量，选取网格量较小的网格B作为计算网格。

图1 翼型局部振动的计算网格Fig.1 Calculation grid of airfoil with local oscillation

表1 3种不同网格的参数Table 1 Parameters of three types of grids

网格网格数量y+网格A1434510网格B409401网格C743850.1

图2 Eppler387翼型试验[19]与计算得到的气动参数的对比Fig.2 Comparison between CFD results and experimental results[19]of Eppler387 airfoil

为验证计算方法对非定常流动计算的准确性，数值模拟了NACA0015翼型俯仰振荡的非定常流场，并将计算结果与美国NASA Ames研究中心Piziali的试验结果[20]进行比较。通过非定常数值计算得到了升力系数的时滞曲线。计算时流动参数与试验参数保持一致，来流马赫数Ma=0.29，雷诺数Re=1.95×106。参考状态翼型迎角为4°，振幅为4.2°，振动频率为10 Hz。如图3所示，计算得到的升力系数总体上与试验数据符合得较好，说明本研究所采用的非定常数值计算方法具有足够的准确性。

图3 NACA0015翼型升力系数时滞曲线Fig.3 Hysteresis cycle of lift coefficient of NACA0015 airfoil

1.2 网格变形方法

本文采用基于弹簧光顺(Spring Smoothing)的动网格方法实现翼型表面局部变形[21]。弹簧光顺假设中，每两个节点之间的网格线被看做理想的弹簧。边界运动之前各个网格线的初始距离构成了网格的平衡状态。给定边界节点的位移会产生与该节点连接的所有弹簧位移成正比的力。由胡克定律得，该网格节点上的力为

(4)

式中：Δxi和Δxj为节点i和邻近节点j的位移；ni为与节点i相连接的节点数；kij为节点i和节点j之间弹簧的弹性系数。

平衡状态时所有弹簧所连接的节点上的合力必须为零，由此得到迭代方程为

(5)

式中：m为迭代次数。

由于边界节点位置更新后边界位移已知，利用Jacobi迭代方法可解得式(5)。求解收敛后，位置更新为

(6)

式中：n+1和n分别表示下一个时间步和当前时间步。

由于翼型蒙皮变形量较小，采用弹簧变形法足以保证本文的变形需求，并且较其他方法可减少计算工作量，提高计算效率。

2 力学模型

本文选取一种针对高空无人机优化设计的层流翼型SFYT7315作为研究对象，该翼型相对厚度为15.3%，最大厚度位置在30.7%弦长处，相对弯度为7.0%，翼型SFYT7315局部振动模型如图 4所示。本文研究中翼型弦长为c，振动表面弦向长度为0.1c，局部振动表面中心距翼型前缘距离分别为0.1c、0.2c、0.3c、0.4c、0.5c、0.6c、0.7c、0.8c。

在翼型局部振动表面施加周期性简谐运动，振动位移量表达式为

(7)

式中：Δy为垂直于弦向的振动位移；L为振动表面弦向长度；xc为振动点距振动起始位置的弦向距离；A为振动幅值；f为周期振动频率；t为时间。

假设机翼局部表面振动周期为T，在一个周期内的不同时刻，机翼局部表面变形情况如图5所示。

图4 翼型局部振动示意图Fig.4 Diagram of localized oscillation of airfoil

图5 一周期内翼型局部表面变形情况Fig.5 Local deformation of airfoil during one period

在翼型表面变形过程中，网格疏密会发生变化，可能对计算结果产生影响。为研究网格分布变化对计算结果的影响，以0.1c位置变形为例，在振幅A=0.001c条件下，局部蒙皮变形至最低点时，网格变得最为稀疏，对该情况下的网格在垂直于壁面方向进行加密，保证网格分布情况与变形前相同。加密前后振动表面附近网格分布如图6所示。

在翼型雷诺数Re=3×105、迎角为3°、表面固定的条件下分别对加密网格和未加密的初始网格进行计算，结果如表 2所示，表中：K=CL/CD为升阻比。

表2中给出的计算结果显示，与加密后网格的计算结果相比，未加密的初始网格计算得到的气动力系数误差很小，均在1%以内。由此可知，翼型表面振动过程中由于振幅较小，因此网格分布变化很小，由网格分布变化带来的计算误差可以忽略。

图6 振动表面附近网格分布Fig.6 Grids distribution near oscillation surface

表2 加密网格和初始网格的计算结果

Table 2 Calculation results of denser grid and original grid

网格CLCDK加密网格1.23470.0173771.0823初始网格1.23420.0175070.5257相对误差/%-0.04050.748-0.783

3 计算结果分析

本文分别计算了局部振动区域位于翼型上表面不同弦向位置，以不同振动频率和不同振幅振动时对翼型气动特性的影响。翼型雷诺数Re=3×105，迎角为3°，局部蒙皮振幅A=0.000 2c～ 0.001c。研究过程中翼型气动力参数主要有升力系数CL、阻力系数CD和升阻比K，下标base表示刚性翼型时均气动力。

3.1 振动位置对翼型气动特性的影响

为了研究振动弦向位置对低雷诺数翼型气动特性的影响，在振动频率f=5 Hz、振幅A=0.001c的条件下，对比不同弦向位置局部振动时翼型的气动特性，翼型时均气动力随振动表面弦向位置的变化如图7所示。图7中：LS为刚性翼型时均层流分离位置，RT为刚性翼型时均再附位置，LS与RT之间即翼型上表面层流分离泡。

由图7可知，振动表面弦向位置为0.1c和0.6c时，翼型升力系数和升阻比明显升高，阻力系数明显降低。在振动频率f=5 Hz条件下，0.1c位置振动时，翼型升力系数提高了4.87%，阻力系数降低了15.29%，升阻比提高了23.79%；0.6c位置振动时，翼型升力系数提高了4.94%，阻力系数降低了15.15%，升阻比提高了23.68%。振动位置在0.1c和0.6c时翼型所获得的气动力变化量相似，即0.1c和0.6c位置表面振动产生的增升减阻效果相同。

由图7可知，刚性翼型时均层流分离位置在0.506c处，时均再附位置在0.721c处，0.6c振动位置在翼型上表面层流分离泡中心附近，由以上研究可知翼型表面振动位置位于翼型前缘附近及层流分离泡中心附近时具有最好的增升减阻作用。在振动频率f=5 Hz条件下，不同位置振动时翼型上表面流线如图8所示。

由图8可知，在雷诺数Re=3×105条件下，刚性翼型上表面0.506c～0.721c范围内形成多个层流分离泡，层流分离范围较大，达到0.215c。0.5c位置振动时，层流分离泡形态变化不明显，仍有较大的层流分离范围。当0.1c和0.6c位置振动时，翼型上表面层流分离区域由多个小分离泡变为一个分离泡，分离范围缩小。这是由于0.1c位于机翼前缘附近，振动可以有效地向附面层内注入能量，并可以对下游较远范围产生影响，对层流分离泡造成有利干扰；0.6c位于层流分离泡内部，振动可直接为分离泡内部低速气流注入能量从而减弱分离。

图7 f=5 Hz时振动位置对翼型气动特性的影响Fig.7 Influence of oscillation location on aerodynamic characteristics of airfoil for f=5 Hz

图8 不同振动位置翼型流线图Fig.8 Streamlines of airfoil with oscillation at different locations

图9 不同振动频率下翼型气动力随时间变化曲线Fig.9 Variation of aerodynamic force of airfoil with time under different oscillation frequencies

在振动频率f=5 Hz条件下，0.1c和0.6c位置振动时翼型振动最大位移量Δym及气动力随时间变化曲线如图9所示。由图9可知，随着局部蒙皮振动，翼型升力系数及阻力系数呈现出周期性变化。0.1c和0.6c位置振动时，一个周期内翼型升力系数均高于刚性翼型，阻力系数均低于刚性翼型，局部振动可有效提高低雷诺数下翼型的气动特性。0.6c位置振动引起翼型气动力系数波动幅度大于0.1c位置振动，这是由于0.6c位置振动对翼型上表面层流分离泡的影响更明显，因此0.6c位置振动使翼型气动力非定常效应更显著。

通过气动特性曲线和振动位移曲线对比可知，翼型升力系数及阻力系数的变化与表面振动位移具有相同的变化趋势，气动力曲线与位移曲线之间存在一定的时差效应。与0.6c位置振动相比,0.1c位置振动的时差效应更显著，这是由于0.1c位置振动对附面层产生的影响传播至下游需要一定时间，而0.6c位置振动直接作用于层流分离泡，时差较小。

0.6c位置振动时，一个周期内翼型表面层流分离泡附近的流场结构如图10所示。由图10可知，随着局部蒙皮振动，翼型上表面层流分离泡在不同时刻呈现出不同形态，正是由于层流分离泡的周期性变化，才造成了翼型气动力系数的周期性波动。

3.2 振动频率对翼型气动特性的影响

为研究局部蒙皮振动频率对翼型气动特性的影响，以0.1c和0.6c位置振动为例，在振幅A=0.001c条件下翼型时均气动力系数随振动频率的变化曲线如图11和图12所示。

图10 一个周期内翼型流场结构Fig.10 Flow structures of airfoil during one period

图11 翼型气动力随振动频率变化曲线(振动位置：0.1c)Fig.11 Variation curves of aerodynamic force of airfoil with oscillation frequency (oscillation location:0.1c)

图12 翼型气动力随振动频率变化曲线(振动位置：0.6c)Fig.12 Variation curves of aerodynamic force of airfoil with oscillation frequency (oscillation location:0.6c)

从图11和图12中可看出，不同位置振动时，振动频率对翼型气动特性均有显著影响。由图11可知，0.1c位置振动时，随着振动频率的增加，翼型升力系数具有持续增加的趋势，阻力系数先减小后增加，升阻比先增加后减小。当振动频率为f=6 Hz时翼型升阻比达到最大值，此时与刚性翼型相比，升阻比提高23.8%，升力系数提高4.9%，阻力系数降低15.3%。由图12可知，0.6c位置振动时，随着振动频率的增加，翼型升力系数和升阻比出现一致的变化规律，即先增加后减小，阻力系数先减小后增加。当振动频率f=4 Hz时翼型升阻比达到最大值，此时与刚性翼型相比，升阻比提高23.7%，升力系数提高4.9%，阻力系数降低15.2%。

由以上分析可知，0.1c和0.6c位置局部振动对翼型有显著的增升减阻作用，其中减阻作用更加明显。为进一步研究蒙皮振动对翼型阻力特性的影响，图13给出了0.6c位置振动时翼型时均压差阻力及摩擦阻力随振动频率的变化曲线，图中：CDp和CDv为振动状态下翼型的时均压差阻力系数及时均摩擦阻力系数；CDpbase和CDvbase为刚性翼型的时均压差阻力系数及时均摩擦阻力系数。由图13可知0.6c位置振动时，翼型压差阻力减小，摩擦阻力增加。随着振动频率的提高，翼型压差阻力先减小后增加，而摩擦阻力先增加后减小。当振动频率f=4 Hz时，翼型压差阻力达到最小值，比刚性翼型降低23.3%，而摩擦阻力达到最大值，比刚性翼型提高11.6%。由于振动状态下翼型压差阻力的减小量大于摩擦阻力的增加量，因此翼型的总阻力降低15.2%。

0.6c位置振动时不同振动频率下翼型上表面湍动能分布如图14所示，上表面剪切力分布如图15所示。

图14中湍动能分布表征了翼型上表面转捩位置，由图14可知，局部振动使翼型上表面转捩位置提前，当振动频率由f=1 Hz增加至f=4 Hz时，转捩位置向前移动，与摩擦阻力随频率的变化规律一致。刚性翼型上表面时均转捩位置在0.656c，振动频率f=4 Hz时，转捩位置在0.608c，转捩位置向前移动了0.048c，导致翼型摩擦阻力系数增加。由图15可知局部振动对翼型振动区域上游的剪切力分布没有明显影响，而使翼型振动区域下游的表面摩擦力分布出现显著变化，随着振动频率由f=1 Hz增加至f=4 Hz，翼型表面剪切力显著增加且均高于刚性翼型，这是翼型摩擦阻力变化的重要原因。

图13 0.6c位置振动时压差阻力系数和摩擦阻力 系数随振动频率的变化Fig.13 Variation of pressure drag coefficient and viscous drag coefficients of airfoil with oscillation frequency when oscillation location is 0.6c

翼型上表面层流分离泡对翼型气动特性有显著影响，因此为进一步研究局部振动对翼型气动特性的影响机理，需研究局部振动对翼型层流分离特性的影响。0.6c位置振动时，不同振动频率下翼型上表面时均流线如图16所示。

图15和图16综合反映了翼型上表面层流分离泡随振动频率的变化情况，图 15中表面剪切力为负值的区域即为分离区，图15和图16表明，0.6c位置振动时，当振动频率由f=1 Hz增加至f=4 Hz的过程中，翼型上表面层流分离泡逐渐缩小，翼型等效外形发生变化，导致翼型气动特性

图14 0.6c位置振动时不同振动频率下翼型上表面 湍动能分布Fig.14 Turbulence kinetic energy on upper surface of airfoil under different frequency when oscillation location is 0.6c

图15 0.6c位置振动时不同振动频率下翼型上表面 剪切力分布 Fig.15 Wall shear on upper surface of airfoil under different frequency when oscillation location is 0.6c

图16 0.6c位置振动时在不同振动频率下翼型上表面流线Fig.16 Streamlines of upper surface of airfoil under oscillation with different frequency when oscillation location is 0.6c

发生变化。刚性翼型及0.6c位置振动频率f=4 Hz时的等效翼型如图17所示。

注： 黑色实线为原翼型；红色实线近似为层流分离泡和附面层共同作用下的等效翼型。图17 分离泡作用下的等效翼型Fig.17 Equivalent airfoils with separation bubble

等效翼型是指低雷诺数条件下，由原翼型与层流分离泡和附面层构成的假想翼型边界。等效翼型表面流体无分离流动，且与压力分布相比，黏性作用影响不大。等效翼型的无黏流动与原翼型的低雷诺数流动近似相同，因此可通过等效翼型的线性化分析，定性解释和分析原翼型的气动特性变化[22]。与刚性翼型的等效翼型相比，局部振动条件下，翼型上表面层流分离泡变小，等效翼型相对厚度减小，这是由于局部振动为上翼面附面层注入能量，使层流分离减弱。由翼型线化理论可知翼型厚度对翼型压差阻力有重要的影响，局部振动使翼型上表面层流分离范围缩小，从而等效翼型相对厚度减小，因此翼型压差阻力降低。

3.3 振幅对翼型气动特性的影响

在振动位置为0.1c和0.6c、振动频率f=5 Hz的条件下研究蒙皮振幅对翼型气动特性的影响，翼型时均气动力系数随振幅的变化曲线如图18所示。当局部蒙皮振幅在A=0.000 2c～0.001c范围时，随着振幅的增加，翼型时均升力系数及升阻比均增加。在0.1c位置振动条件下，当振幅由A=0.000 2c增加至0.001c时，升阻比提高0.996%；在0.6c位置振动条件下，当振幅由A=0.000 2c增加至0.001c时，升阻比提高0.916%。

图19给出了0.6c位置振动时在不同振幅条件下，翼型气动力随时间的变化曲线。从图中可以看出，不同振幅条件下，翼型气动力均呈现出周期性变化，随着振幅的增加翼型气动力波动幅度增大，流场非定常特性更加显著。当振幅由A=0.000 2c增加至A=0.001c时，局部振动为附面层注入更多能量，对层流分离的控制效果更好，翼型升力系数提高，阻力系数降低。

图18 翼型气动力随振幅的变化曲线Fig.18 Variation of aerodynamic force of airfoil with amplitude

图19 不同振幅下翼型气动力随时间变化曲线Fig.19 Variation curves of aerodynamic force of airfoil under different amplitude with time

4 结 论

1) 表面局部振动可有效提高低雷诺数翼型气动特性，与刚性翼型相比，表面局部振动可使机翼升力系数提高4.9%，阻力系数降低15.3%，升阻比提高23.8%。

2) 振动位置对翼型气动特性及流场结构有显著的影响，振动表面位于翼型前缘附近或位于层流分离泡中心时可有效控制翼型层流分离，从而提高翼型升力系数，降低阻力系数。

3) 振动频率对翼型表面层流分离及转捩位置均有显著的影响，随着振动频率增加，翼型气动特性出现最优值。与刚性翼型相比，表面振动使翼型转捩位置略向上游移动，摩擦阻力增加，且随频率的增加转捩位置越靠前，但振动使等效翼型相对厚度减小，压差阻力明显减小。

4) 当局部蒙皮振幅在A=0.000 2c～0.001c范围时，随着振幅的增加，流场非定常特性更加显著，翼型升力系数及升阻比增加。

[1] COLOZZA A, DOLCE J L. High-altitude, long-endurance airships for coastal surveillance: NASA/TM-2005-213427[R]. Washington, D.C.: NASA, 2005.

[2] DONGLI M, YANPING Z, YUHANG Q. Effects of relative thickness on aerodynamic characteristics of airfoil at a low Reynolds number[J]. Chinese Journal of Aeronautics, 2015, 28(4): 1003-1015.

[3] SELIG M S, GUGLIELMO J J, BROERN A P, et al. Experiments on airfoils at low Reynolds numbers[C]∥34th Aerospace Sciences Meeting and Exhibit. Reston, VA: AIAA, 1996.

[4] 刘强， 刘周， 白鹏. 低雷诺数翼型蒙皮主动振动气动特性及流场结构数值研究[J]. 力学学报, 2016, 48(2): 269-277.

LIU Q, LIU Z, BAI P. Numerical study about aerodynamic characteristics and flow field structures for a skin of airfoil with active oscillation at low Reynolds number[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 269-277 (in Chinese).

[5] GREENBLATT D, WYGNANSKI I J. The control of flow separation by periodic excitation[J]. Progress in Aerospace Sciences, 2000, 36(7): 487-545.

[6] REDHA W, DAVID B. Suction effects on the transition and reattachment of a transitional bubble[C]∥51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Reston, VA: AIAA,2013.

[7] WRIGHT M C M, NELSON P A. Wind tunnel experiments on the optimization of distributed suction for laminar flow control[J]. Journal of Aerospace Engineering, 2001, 215(6): 343-354.

[8] REDHA W, DAVID B. Effects of distributed suction on an airfoil at low Reynolds number[C]∥40th Fluid Dynamics Conference and Exhibit. Reston, VA: AIAA, 2010.

[9] LIAN Y S, WEI S. Membrane wing aerodynamics for micro air vehicles[J]. Progress in Aerospace Sciences, 2003, 39(6): 425-465.

[10] 康伟， 刘磊， 徐敏. 低雷诺数下翼面局部振动增升机理研究[J]. 航空学报, 2015, 36(11): 3557-3566.

KANG W, LIU L, XU M. Lift enhancement mechanism for local oscillation of airfoil surface at low Reynolds number[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(11): 3557-3566 (in Chinese).

[11] GORDNIER R E. High fidelity computational simulation of a membrane wing airfoil[J]. Journal of Fluids & Structures, 2009, 25(5): 897-917.

[12] SONG A, BREUER K. Dynamics of a compliant membrane as related to mammalian flight[C]∥AIAA Aerospace Sciences Meeting & Exhibit. Reston, VA: AIAA, 2007.

[13] ROJRATSIRIKUL P, WANG Z, GURSUL I. Unsteady fluid-structure interactions of membrane airfoils at low Reynolds numbers[M]. Berlin: Springer, 2010: 297-310.

[14] 裘进浩， 李大伟， 聂瑞. 增加翼型升力的局部振动流动控制技术[J]. 南京航空航天大学学报, 2012, 44(5): 638-644.

QIU J H, LI D W, NIE R. Vibration flow control technology for increasing lift of airfoil[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2012, 44(5): 638-644 (in Chinese).

[15] 叶正寅， 谢飞. 弹性振动对翼型失速迎角附近流场的影响[J]. 航空学报, 2006, 27(6): 1028-1032.

YE Z Y, XIE F. The effects of elastic vibration on the flow field near stall-incidence of the airfoil[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(6): 1028-1032 (in Chinese).

[16] QIAO Y H, MA D L, LI Z. Unsteady numerical simulation of the mutual disturbance between propeller and wing[J]. Transactions of the Japan Society for Aeronautical & Space Sciences, 2015, 58(5): 298-306.

[17] LANGTRY R, MENTER F. Transition modeling for general CFD applications in aeronautics[C]∥43rd AIAA Aerospace Sciences Meeting and Exhibit. Reston, VA: AIAA, 2005.

[18] MENTER F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[19] ROBERT J M, BETTY S W, BETTY F M. Experimental results for the Eppler 387 airfoil at low Reynolds numbers in the Langley low-turbulence pressure tunnel: NASA/TM-4062[R]. Washington, D.C: NASA Langley Research Center, 1988.

[20] PIZIALI R A. 2-D and 3-D oscillating wing aerodynamics for a range of angles of attack including stall[R]. Washington, D. C.: NASA, 1994.

[21] BLOM F J. Considerations on the spring analogy[J]. International Journal for Numerical Methods in Fluids, 2000, 32(6): 647-668.

[22] 白鹏， 崔尔杰， 李锋. 对称翼型低雷诺数小攻角升力系数非线性现象研究[J]. 力学学报, 2006, 38(1): 1-8.

BAI P, CUI E J, LI F. Study of the non-linear lift coefficient of the symmetric airfoil at low Reynolds number near the 0° angle of attack[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 1-8 (in Chinese).