张兴伟，周超英，谢 鹏

(哈尔滨工业大学 深圳研究生院，518055广东 深圳，cyzhou@hit.edu.cn)

扑翼的气动特性研究是设计、优化和控制微型扑动飞行器(Micro Air Vehicle，MAV)的重要环节.近年来扑翼空气动力学特性的研究有了很大的发展，但对于MAV所处的低雷诺数下的非定常气动特性、扑翼柔性变形对气动特性的影响等一些突出问题还有待进一步研究.特别是随着材料、生物电池等技术的进步和发展，扑动飞行器进一步微型化、柔性化的发展趋势日趋明显，这些新型MAV的尺度、飞行雷诺数等参数都与小型昆虫的相似，因此，研究昆虫扑翅获取高升力和推力的飞行机制能为发展MAV提供一些参考和灵感，其中研究低雷诺数下柔性翅的气动特性尤为重要.近年来国内外学者对昆虫的扑翼飞行方式及其飞行机制进行了深入研究，取得了显著成果[1-5].

昆虫在飞行中通过扑动翅膀来调整升力和推力从而改变其飞行速度和高度.近年来人们通过对多种飞行昆虫的实验观测发现昆虫扑翼在扑动飞行中会产生不同程度的弯曲和扭转，不同昆虫这种弯曲和扭转的方式和程度也不同，这些变形对昆虫扑动飞行性能有显著的影响.此外，在MAV样机试验中也发现，采用仿生材料研制的扑动机翼在拍动过程中产生的柔性变形对扑翼的气动特性有着重要影响.早期的数值研究大都假设翅翼为刚性无变形，随着计算机技术和流体力学计算方法的进步和发展，越来越多的数值研究关注柔性变形对扑翼气动特性的影响.本文采用流体力学数值计算方法有限体积法求解Navier-Stokes方程，对柔性扑翼在悬停时的气动特性进行了研究，分析了扑翼的柔性变形对其悬停飞行性能的影响.一般而言昆虫的扑动飞行是一个三维问题，在考虑展向尺度足够长以及无展向扭转变形的前提下，可忽略展向流影响从而将翼型简化成二维模型来研究其扑动的非定常效应，对此Wang等[6]验证了二维模型数值计算结果与三维实验结果的可比性.为简化问题，在此选取的柔性扑翼模型也为二维模型.

1 可变形扑翼物理模型及运动模式

建立弦长为c，最大厚度为弦长2%，截面为椭圆的二维扑翼模型，翼的悬停飞行及其柔性变形如图1所示，其中:坐标xoy为固联于地面的惯性坐标系;x'ory'为固定在翼表面的移动坐标系;坐标原点or固定于距翼前缘c/4处.在悬停扑动过程中，x'ory'坐标系相对于惯性坐标系xoy做往复平动和转动.

图1 悬停中扑翼往复平动和转动示意图

翼的悬停扑动由往复平动和转动组成.在xoy坐标系下，平动位移H和转动转角α描述为

式中:A为扑动幅值;f为扑动频率;T为扑动周期(T=1/f);α0为转角幅值.通过实验和数值研究表明，当A/c=3～5时，扑翼前缘处产生的涡即前缘涡不会从翼表面完全分离，从而产生较大的升力[6－7].因此，在计算中 A 取值为 4.5c，α0取值为45°.α(t)与瞬时攻角有90°的相位差，即瞬时攻角 =90°－|α(t)|.扑动平面为水平面.

昆虫翅膀主要由位于前缘的翅梁、翅脉和翅膜组成.翅膀前缘的刚度比较大，在飞行中形变非常小可近似成刚性;后缘刚度较小，悬停时会产生相对来说较大的变形.对所采用的椭圆截面翼型，考虑到翼的弦向变形尺度相对于弦长很小，因此，假设变形前后翼的弦长和厚度均保持不变;翼的前缘部分不发生变形，后部分的变形起点位于距前缘1/4弦长处，即x'ory'坐标中的原点or点.在悬停过程中，从or点到翼前缘这部分不发生变形，仅做平动和转动;从or点到翼后缘这部分除了相对于坐标系xoy做平动和转动外，同时还产生在x'ory'坐标系下相对于or点的柔性变形.图2是翼型变形的示意图，d/c=0.25，翼型的后缘点为最大变形处.取弦向的最大弯曲尺度为a0c，其中无量纲系数a0表征了后缘端点相对于弦长的变形程度.

图2 悬停中翼变形的示意图

本文主要研究翼型弦向的动态弯曲变形对悬停气动特性的影响，因此计算中没有考虑由于展向扭转而引起的展向流影响.翼型在往复平动和转动复合运动中产生的动态柔性形变可用Miao＆Ho[8]的模型进行描述，即在 x'ory'坐标系下，从or点到翼后缘部分翼型的弦向变形表述为

式中 φ 为相位角.Read[9]，Heathcote[10－11]等发现，柔性翼的最优相位角为90°，因此，在此动态柔性模型的相位角φ取值为90°.b0为变形系数，当弯曲的起点or位于c/4处时，若翼型后缘点最大弯曲达0.1c，即取a0=0.1，则变形系数值为b0≈ 0.177 8.

利用坐标转换关系式为

即可实现翼型变形变量从x'ory'坐标系到xoy坐标系的转换.

2 计算方法

均匀不可压缩流体运动的控制方程为Navier-Stokes方程

式中:ui为流速，二维情况下 i=1，2，j=1，2;P为动压;雷诺数定义为Re=ρcuref/μ，其中ρ为流体密度，μ为流体粘性系数，uref为扑动平面内的平均平动速度.

采用有限体积法[12]，在控制体积中对上述控制方程进行积分，采用交错网格[13]实施积分方程在有限单元上的数值离散，应用PISO压力速度耦合法[14-15]迭代求解压力场和速度场.时间项采用一阶差分格式，对流项采用精度为二阶的迎风差分格式，粘性项用二阶中心差分格式.计算区域四周外边界取在距翼型壁面30 c处，边界条件给定为远场无来流条件;翼型表面边界条件为无滑移边界条件.计算中选取了三角型网格.采用动网格技术实现了对扑动和变形动边界的模拟，网格的运动变形满足空间守恒定律[16].计算中首先对网格密度以及时间步长对计算结果精度的影响进行了大量试算和验算，最终确定时间步长为Δt=0.001T，网格初始节点数为3.8 × 104.

翼型的升力系数CL和阻力系数CD分别定义为

式中:urms为平动速度的均方根值;FL和FD分别为瞬时升力和阻力.扑动周期内的平均升、阻力系数分别为

3 计算结果及分析

3.1 程序验证

首先，通过计算刚性翼悬停时的气动力及流场分布，检验了动网格计算的可靠性.图3给出了Re=75时悬停刚性翼的气动力计算结果，其翼型运动参数参照了文献[6]的模型参数.结果显示，所采用动网格方法的计算结果与文献[6]的结果符合较好.

3.2 弦向弯曲变化对气动力的影响

计算中假定了翼型不会因为流场变化而产生自适应的变形.为了研究柔性翼在悬停时不同变形尺度对气动力产生的影响，分别计算了6种不同变形幅度翼型在Re=140时的气动力及流场分布，其中 a0为 0.00，0.05，0.10，0.15，0.20 和0.30.

图4给出了a0=0.05时柔性翼在悬停运动状态下升、阻力系数的时程图，并与相同雷诺数下刚性翼的结果进行了比较.结果显示，柔性翼升、阻力系数的瞬时峰值分别由刚性翼的3.5和4.6提高到4.8和5.7，显然悬停飞行中柔性弯曲对瞬时气动力峰值的影响非常显著.另外，柔性弯曲变形没有改变瞬时气动力的分布，即在一个扑动周期内，刚性翼和柔性翼升、阻力峰值发生的时刻基本相同，显示旋转环流、尾迹捕捉等非定常效应与翼型的变形机制相比仍然是决定非定常气动力分布的主要因素.

图3 刚性翼悬停时升、阻力系数时程曲线(Re=75)

为研究柔性变形对气动特性的影响，对升、阻力系数的周期时均值及升阻比进行了比较分析.图5 为Re=140 时变形尺度系数a0为0.00，0.05，0.10，0.15，0.20和0.30时升、阻力系数时均值随a0的变化曲线.结果表明，当后缘产生一定程度弯曲时如a0=0.05时，虽然升、阻力系数时均值的变化幅度同瞬时值的变化幅度相比较小，但柔性变形还是使翼翅的时均气动性能有了一定的改变，升、阻力系数时均值略有增大，因此适当的柔性弯曲可以改善翼的悬停特性，这一结果与金晓怡等[17]试验观察的现象一致;但是，随着后缘点弯曲幅度a0继续增大，升、阻力系数时均值快速减小，Zhao等[18]通过样机试验发现了相同的现象，他们还发现当翼翅保持某一固定弯曲幅值时，阻力系数时均值随功角转动幅值增大而增大.图6为升阻比随变形尺度a0的变化规律，结果显示升阻比变化趋势与升、阻力系数时均值的变化趋势基本相同，柔性翼在a0=0.05，0.10，0.15 时得到了高过刚性翼的升阻比，但随着弯曲度的进一步增大，升阻比迅速减小，气动性能恶化.

图4 扑翼悬停时升、阻力系数时程曲线(Re=140)

图5 柔性翼时均升、阻力系数随变形幅度的变化曲线

图6 升阻比随变形幅度的变化曲线

3.3 弦向弯曲程度变化对涡量分布影响

图7和图8分别给出了刚性翼和不同变形幅度柔性翼在1/2扑动周期内的几个典型时刻的涡量分布图.当翼型位于最左端位置时(图7和图8中t/T=0.5)，上半个周期产生的前缘涡位于翼型右半面的前缘附近.随着扑翼向右运动，前缘涡向后缘发展，最后脱落产生较长的尾涡;同时，翼左半面的上缘处逐渐产生新的前缘涡(图7和图8 中 t/T=0.625).

图7 刚性翼在不同时刻的涡量分布图

对比图7，图8可知，翼的柔性变形使得扑动过程中产生的前缘涡涡强变弱(图7中t/T=0.875和图8(a)、图8(b)中 t/T=0.875).因此，相应的升、阻力时均值与刚性翼的相比较小.这一现象与文献[18]通过样机实验测量的研究结果相一致.柔性翼流场与刚性翼流场相比，后缘处的涡强变弱，且后缘处在半个周期内产生的尾迹在x方 向变长(图8(a)、(b)中t/T=1.0)，显示柔性弯曲不但改变了前缘、后缘涡的强度也改变了尾迹脱落涡的传输方向.

图8 柔性翼在不同时刻的涡量分布图

4 结论

1)扑翼的弦向弯曲变形能改变瞬时气动力的峰值.以后缘最大变形尺度为0.05倍弦长的柔性扑翼为例，升力和阻力最大峰值的增加分别可达37%和24%.

2)扑翼弦向弯曲变形对气动力时均值的影响相对瞬时值较小.适度的弯曲变形可以有效增大升力从而获得较高的升阻比，但过大的柔性弯曲会导致翼型的气动特性急剧恶化，升、阻力系数时均值及升阻比迅速减小.

3)扑翼扑动产生的非定常效应同翼的弯曲变形机制相比仍是决定气动力分布的主要因素.但柔性变形能在一定程度上改变前缘和后缘涡的涡强，以及脱落涡的传输方向.

[1]SANE S P.The aerodynamics of insect flight[J].Journal of Experimental Biology，2003，206(23):4191-4208.

[2]WANG Z J.Dissecting insect flight[J].Annual Review of Fluid Mechanics，2004，37(1):183-210.

[3]鲍麟，童秉纲.模型翼拍动中动态柔性变形效应的数值研究[J].中国科学院研究生院学报，2005，22(6):676-684.

[4]WALKER S M，THOMAS A L R，TAYLOR G K.Deformable wing kinematics in free-flying hoverflies[J].Journal of The Royal Society Interface，2010，7(42):131-143.

[5]周超英，林玉峰.扑翼气动力特性的数值研究[J].哈尔滨工业大学学报，2006，38(9):1403-1405.

[6]WANG Z J，BIRCH J M，DICKINSON M H.Unsteady forces and flows in low Reynolds number hovering flight:Two-dimensional computations vs robotic wing experiments[J].Journal of Experimental Biology，2004，207(3):449-460.

[7]DICKINSON M H，GÖTZ K G.Unsteady aerodynamic performance of model wings at low Reynolds numbers[J].Journal of Experimental Biology，1993，174:45-64.

[8]MIAO J M，HO M H.Effect of flexure on aerodynamic propulsive efficiency of flapping flexible airfoil[J].Journal of Fluids and Structures，2006，33(3):401-419.

[9]READ D A，HOVER F S，TRIANTAFYLLOU M S.Forces on oscillating foils for propulsion and maneuvering[J].Journal of Fluids and Structures，2003，17(1):163-183.

[10]HEATHCOTE S，MARTIN D，GURSUL I.Flexible flapping airfoil propulsion at zero freestream velocity[J].AIAA Journal，2004，42(11):2196-2204.

[11]HEATHCOTE S，GURSUL I.Flexible flapping airfoil propulsion at low Reynolds numbers[J].AIAA Journal，2007，45(5):1066-1079.

[12]VERSTEEG H K，MALALASEKERA W.An Introduction to Computational Fluid Dynamics:The Finite Volume Method[M].New York:Wiley，1995.

[13]王福军.计算流体动力学分析:CFD软件原理与应用[M].北京:清华大学出版社，2004.

[14]ISSA R I.Solution of the implicitly discretized fluid flow equations by operator-splitting[J].Journal of Computational Physics，1985，62(1):40-65.

[15]ISSA R I，AHMADI-BEFRUI B，BESHAY K R，et al.Solution of the implicitly discretized reacting flow equations by operator-splitting[J].Journal of Computational Physics，1991，93:388-410.

[16]LESOINNE M，FARHAT C.Geometric conservation laws for flow problems with moving boundaries and deformable meshes and their impact on aeroelastic computations[J].Computer Methods in Applied Mechanics and Engineering，1996，134(1/2):71-90.

[17]金晓怡，颜景平，周建华.仿生扑翼飞行机器人柔性翅试验研究及其理论解释[J].中国机械工程，2007，18(9):1028-1031.

[18]ZHAO L，HUANG Q，DENG X，et al.Aerodynamic effects of flexibility in flapping wings[J].Journal of the Royal Society Interface，2010，7(44):485-497.