“Click”运动高升力机制

夏 震 昊, 唐 斌, 于 紫 江

(大连理工大学 能源与动力学院, 辽宁 大连 116024 )

0 引 言

昆虫的飞行方式一直吸引着研究人员的注意．在研究的早期,按照定常流动分析昆虫悬停的结果是大黄蜂不会飞[1]．Weis-Fogh[2]的研究促使人们从流动的非定常方向思考昆虫产生高升力的问题．Ellington等[3]通过对鹰蛾翅膀的拍动模型进行可视化处理,观察到模型上存在展向流动,展向流动稳定了失速涡,使鹰蛾翅膀能够继续保持高升力．Liu等[4]通过计算流体力学方法,分析了展向流动稳定失速涡所形成的前缘涡,并解释了各阶段涡量的不同形态．后续研究发现,随着雷诺数Re的变化,前缘涡的结构和稳定性会发生变化[5]．Birch等[6]通过实验分别研究了移除涡流核心能量和通过尖涡与尾涡引起的向下流动降低有效攻角这两种假说,分析了它们对限制前缘涡增长、保持其稳定性的影响．Lentink等[7]通过可视化模型翼上的流动特性发现,在低Rossby数下,翼的扫掠产生了向心和科里奥利加速度,从而使前缘涡得到稳定．Jardin[8]的研究表明,当Re超过200时,科里奥利效应主导了流动稳定性,并且离心效应对流场没有显著影响．当Re低于200时,黏性效应的作用变得突出．近年来,Chen等[9]对旋翼的研究结果表明,当Re从1 000降低到10时,前缘涡形成周期的持续时间和稳态强度均显著下降．此外,有助于前缘涡稳定性的主要机制可以在不同Re水平上变化．随后,他们将昆虫翅膀的拍动简化为平板俯仰运动,通过研究发现,俯仰速率增加及俯仰轴向后缘移动都会终止前缘涡的增长,再通过强烈的向下对流使其向后缘移动．然而,在这两种情况下,展向流动和涡旋伸展(或压缩)的作用是不同的,从而导致不同的升力产生机制[10]．

在建立扑翼运动的运动学模型时,研究人员经常使用对称的上、下行程来简化翅膀运动学的建模[11-12]．在进行CCl4麻醉的双翅目昆虫解剖研究时,Boettiger等提出了“Click”机制假说[13]．该机制之所以被称为“Click”机制,是因为Boettiger等发现,当作用在双翅目昆虫翅膀结构上的激励力达到一个临界值时,昆虫的翅膀突然快速地从一个终点状态拍动到另一个终点状态,进而完成上(下)拍动行程．产生这种现象的原因是翅膀结构中的双稳态特性,也被国外生物学家使用拟声词称为“Click”机制．Thomson等[14]在前人工作基础上建立了“Click”机制的简化机构模型,使用突变理论进行了运动学分析．Brennan等[15]进一步建立了集中质量动力学模型,完成了动态特性研究．Tang等改进了Brennan等的动力学方程,在不同阻尼特性下进行了模型的比较分析．研究发现,在低于共振频率和二次阻尼条件下,“Click”机构较简谐机构拥有更大的优势[16]．Kok等[17]发现当拍动频率低于固有频率,即频率比为0.3时,使用不对称扑翼微飞行器的性能提高了75%;且在弹性力很大的柔性系统中,不对称性可以提高扑翼驱动系统的气动效率．Tang等[18]制造出一种采用“Click”运动的飞行器模型,并对模型结构进行了动力学研究．研究发现,当以一个激励周期内动能与输入能量比作为判据时,非对称结构的性能略好于对称结构．Bose等[19]改进了Brennan等给出的“Click”机构动力学模型,将其与考虑扑翼和非定常流之间非线性相互作用的势流求解器耦合,研究系统在扑翼响应、空气动力荷载和流场尾流模式方面的分叉行为．

扑翼飞行器的小型化和微型化是其研发的重要目标,而“Click”运动在机构的固有频率较大,即质量较小时,能够凸显双稳态结构的优势．双稳态结构具有在两种稳态之间的突弹跳变现象,较广泛地应用于软体机器人、仿生非线性工程结构等[20-22]．“Click”机构就是双稳态结构的一种,其主要优点是显著降低了体积和能源成本,有利于扑翼飞行器的小微型化．已有研究主要从动力学、升力与阻力及能量消耗等角度[14-19]对“Click”机构进行研究,并与简谐运动的仿生机构进行比较．然而,通过求解三维不可压缩非稳态Navier-Stokes方程,对具有“Click”运动的模型翼的空气动力学及流场特性研究却较少．本文在三维流场中,在模型翼扭转的基础上,分别在水平方向采用“Click”运动和简谐运动,对昆虫悬停状态下模型翼的空气动力学特性进行仿真研究,计算“Click”运动和简谐运动的升力与阻力,并分析拍动过程中翼表面的流线、压力和涡结构变化,探究“Click”运动产生高升力的非定常机制．

1 计算模型

1.1 “Click”机构的动力学模型

图1为“Click”机构示意图．BE和DF为悬臂梁,刚体ABC与刚性杆CD在C点铰接,翼固连于刚体ABC上．图1(a)为结构未变形的初始状态．当C点受到外界激励f(t)时,悬臂梁BE和DF在末端弯曲,储存势能．当C点越过图1(b)所示的不稳定平衡位置时,悬臂梁BE和DF储存的势能得到释放,使C点加速向下移动,最终达到图1(c)状态．

(a) 下拍始点

(b) 中间阶段

(c) 上拍始点

此模型为双稳态结构,存在3个平衡位置,图1(a)、(c)为两悬臂梁未变形时的稳态平衡位置,图1(b)为两悬臂梁处于最大弯曲状态的不稳定平衡位置．在图1(b)状态时,任何一点微小的扰动都会导致结构向另外两种状态快速转变,从而将悬臂梁BE和DF储存的能量在较短时间内释放出来,使固连于刚体ABC上的翼处于快速下拍或快速上拍的运动模式,即“Click”运动．

1.2 翅膀运动的坐标系

如图2所示,引入拍动角(φ)、抬升角(θ)和攻角(α)描述翅膀的运动．本文采用Kolomenskiy等[23]对实验室中大黄蜂悬停姿态研究中的测量数据,并在其基础上进行近似．假设保持拍动角和攻角随时间的变化规律不变,将抬升角近似为0°,即翅膀完全在拍动平面内运动．根据以上描述,建立大黄蜂悬停时翅膀的坐标系．橙色部分为翅膀,其在拍动平面内做周期运动．

图2 翅膀拍动示意图

1.3 坐标变换

如图3所示,描述翼相对于地球的运动涉及两个坐标系之间的关系：在空间建立固连于翼的牵连坐标系(oxyz)和惯性坐标系(OXYZ),其中牵连坐标系的y轴与翼展方向一致,x轴与翼弦方向一致,翼在拍动平面上运动．通过欧拉角变换,将相对于牵连坐标系的坐标转换为相对于惯性坐标系的坐标．

图3 牵连坐标系和惯性坐标系

根据理论力学刚体定点运动的解析描述,牵连坐标系绕着y轴旋转的变换矩阵为

(1)

绕着z轴旋转的变换矩阵为

(2)

因此,由牵连坐标系oxyz关于惯性坐标系OXYZ的欧拉变换矩阵可知,翼的牵连坐标与惯性坐标的关系为

(3)

1.4 简谐运动和“Click”运动

1.4.1 简谐运动 本文在研究中所采用的简谐运动形式为实验观测的近似拟合形式[23]．此时,翼绕z轴在拍动平面内运动的平动表达式为

(4)

翼绕y轴在拍动平面内运动的扭转表达式为

(5)

式中：f为翼的拍动频率,t为时间．使用式(3)～(5),可以将翼的牵连坐标转换为惯性坐标．

1.4.2 “Click”运动 本文重点研究在外部激励频率低于机构固有频率和二次阻尼条件下,“Click”机构的运动方程,并仍使用式(5)描述翼的扭转运动．采用Tang等研究中的运动模型及参数作为与简谐运动对比的运动方程[16]：

(6)

1.4.3 “Click”运动方程的解 使用数值方法求解式(6)所示的“Click”运动方程,在相同角振幅和运动周期条件下,比较简谐和“Click”运动的角位移,如图4所示．其中,时间t与运动周期T=2π/ω的量纲一比值τ为横轴,φ为纵轴．图5为两种运动在翼尖的速度vt的变化曲线,在角位移和运动周期相同的条件下,“Click”运动的最大速度约为简谐运动的2倍．“Click”运动的角位移和速度达到最大时,角位移和速度曲线都产生了波动,且速度曲线的波动更明显．

图4 简谐和“Click”运动的角位移

图5 简谐和“Click”运动的翼尖速度

1.5 求解设置

图6 计算模型及网格划分

对70 mm×70 mm×35 mm的区域采用四面体网格进行加密．对近翼范围进行3层加密,翼表面网格数为21 393,内部采用3层加密,加密部分网格数为4 467 162,未加密部分网格数为2 698 993,总网格数为7 269 405．加密的目的是使流场拥有较高的分辨率,更好地反映近翼范围的空气动力学特性．

采用计算流体力学的压力基数值方法,求解三维不可压缩Navier-Stokes方程．由于本研究模拟的是大黄蜂悬停时的飞行状态,因此远场设置为压力出口边界．流场环境采用Kolomenskiy等[23]研究中的仿真环境,空气密度ρ=1.197 kg/m3,运动黏度ν=1.53×10-5m2/s．

2 模型验证

为了验证模型计算的准确性,采用UDF定义翼的运动方式,计算出升力和阻力,并与Kolomenskiy等[23]的结果进行对比．同时,将翼运动产生的涡结构图与Ellington等[3]的鹰蛾扑翼机构的实验结果进行对比．

图7(a)、(b)分别是升力和阻力随着时间的变化曲线．纵轴rl和rd分别代表双翼升力和阻力与仿真昆虫所受重力比值的量纲一升力和阻力(下称升力和阻力)．本研究中采用量纲一力是为了更直观地比较不同运动方式所产生的力与昆虫所受重力之间的对比．

在前0.5个周期中,Kolomenskiy等的结果较快地出现短暂接近0.2的升力峰值,随后达到高于2的升力峰值,如图7(a)所示．而仿真结果在初期出现短暂向下的趋势之后便一直呈现上升状态,最终达到接近2.4的升力峰值．翼在下拍过程中产生升力峰值,两者走向基本相同．在最大升力峰值处,仿真结果大约比Kolomenskiy等的结果高15%．在后0.5个周期中,Kolomenskiy等的结果同样较早地出现了接近0.3的短暂升力峰值,随后一直增加至最大升力峰值．最大升力峰值处为上拍过程,仿真结果和Kolomenskiy等的结果的极值几乎相同．Kolomenskiy等的平均升力为0.87;仿真的平均升力为1.09,与实际的悬停状态更加接近．

(a) 升力曲线

图7 仿真结果与Kolomenskiy等的结果对比

如图7(b)所示,在仿真和Kolomenskiy等的结果的前0.5个周期中,仿真的阻力峰值高于Kolomenskiy等的阻力峰值约15%．在τ>0.4后,Kolomenskiy等的阻力高于仿真的．随后,在后0.5个周期中,Kolomenskiy等的结果出现向上的阻力峰值,而仿真结果虽然出现阻力向上的趋势,但是幅度较小,在τ接近0.8时两者阻力相同,约为-1.8．在τ>0.8时,仿真和Kolomenskiy等的结果趋势相同,但是前者仍然高于后者约20%．Kolomenskiy等的平均阻力为0.025,仿真的平均阻力为0.03．

Kolomenskiy等的结果和仿真的升力与阻力变化曲线大致相同,但是图7(b)在τ=0.6～0.8时有一定差别．这是由于Kolomenskiy等使用的昆虫翅膀运动的实测数据并不是完全的简谐运动,而仿真使用的式(4)、(5)为简谐运动．

图8(a)为Ellington等[3]的鹰蛾扑翼机构实验结果,图8(b)、(c)为简谐运动的仿真结果,其中图8(b)为正视翼前缘方向,图8(c)为俯视翼前缘方向．

(a) 鹰蛾扑翼机构实验结果

(b) 简谐运动正视翼前缘方向

(c) 简谐运动俯视翼前缘方向

在下拍过程的3个阶段中,仿真结果将涡量等值面与速度进行耦合,得到前缘涡的形态与实验结果基本吻合．图8(a)第一个图为实验下拍开始阶段,此时前缘涡尚未形成,表征前缘涡的烟雾从翼前缘的容器中释放出来,在图8(a)第二、三个图中可以清楚地观察到蓝色的前缘涡．

如图8第二行所示,在下拍中间阶段,仿真结果与实验结果具有相同的螺旋形前缘涡,且翼尖部分前缘涡有脱离迹象．如图8第三行所示,在下拍末期,实验结果显示前缘涡从翼展60%处脱离翼表面并破裂,仿真结果展示了与实验结果相同的前缘涡破裂形态．

将三维仿真结果分别与已有的仿真和实验结果进行对比,可以看出本研究中所采用的计算模型和仿真参数能够达到较高的准确性．

3 结果与讨论

本章将对比分析简谐运动与“Click”运动的升力与阻力,探究两者最大升力处(分别为τ=0.29和0.35)的前缘涡结构和展向的气动特性．

3.1 简谐运动与“Click”运动的升力与阻力

图9为简谐和“Click”运动的升力与阻力曲线．如图9(a)所示,升力曲线也呈现出双峰的特点,即上拍和下拍过程的最大升力峰值．在前0.5个周期中,当τ<0.3时,简谐运动的升力大部分高于“Click”运动的．当τ>0.3时,“Click”运动越过不稳定平衡位置,储存的势能得到快速释放,产生了接近10倍简谐运动的瞬时升力,之后的小幅振荡现象使得升力也产生了波动．这是由“Click”运动的双稳态特性决定的,将此运动作为平动形式应用于翼上,会在运动开始阶段与结束阶段产生小幅振荡现象．后0.5个周期中的升力与前0.5个周期相似．在整个周期中,“Click”运动的平均升力为2.5,简谐运动的平均升力为1.0．由此可见,“Click”运动比简谐运动能够产生更高的升力．如图9(b)所示,两种运动的阻力曲线均从0点开始,先增大到阻力极值后逐渐下降到负值(此时的阻力为推力)．当τ<0.3时,简谐运动的阻力高于“Click”运动的．当τ>0.3时,“Click”运动的阻力迅速增大至18,随后振荡下降至0以下．当τ在0.8左右时,阻力接近极小值-10．当运动的第一个周期结束后,两种运动的阻力回升至0．在整个周期中,“Click”运动的平均阻力为0.35,简谐运动的平均阻力为0.07,均较接近0．在悬停状态下,只有平均阻力接近0,昆虫才能保持水平方向的稳定．

(a) 升力曲线

b) 阻力曲线

3.2 下拍过程

3.2.1 下拍过程中的压力和涡量云图 选取简谐运动和“Click”运动的最大升力处进行定性分析,对比两种运动在不同截面处的压力云图和涡量云图．扑翼运动的升力主要来自于非定常力,其依赖于涡环量的变化及其涡结构的变形[24-25],因此对展向涡量和压力云图的研究是必要的．

涡量是流体速度矢量的旋度,是描述旋涡运动的重要物理量,单位为s-1．三维涡量计算公式为

ω=∇×v=(ωxωyωz)

(7)

图10对比了两种运动沿着翼展方向的压力云图和涡量云图,可以从图中研究不同位置处的前缘涡形态和由前缘涡带来的压力变化．图10(a)、(b)分别为简谐运动的压力云图和涡量云图,图10(c)、(d)分别为“Click”运动的压力云图和涡量云图．

图10 下拍过程中不同截面处的压力和涡量云图

对于简谐运动,随着翼展方向长度的增加,前缘涡从截面2 mm处一直增加到截面10 mm处,在截面10 mm处前缘涡已经逐渐呈现出脱离翼的状态,在截面13.5 mm处,前缘涡完全脱离了翼的前缘．随着翼展方向长度的增加,尾缘涡逐渐增大．由于翼的向前运动,尾缘涡相对于翼存在向后运动的趋势,但是未发生脱离翼的情况．

与简谐运动不同,“Click”运动沿着翼展方向的前缘涡并未脱离翼表面,而是紧紧地依附在前缘处,仅在翼尖附近有脱离翼表面的趋势．为了更好地对比两种运动的压力云图和涡量云图,本文将压力云图和涡量云图的标尺设为相同范围．对比两种运动的压力云图可以看出,由前缘涡的附着而产生的不失速机制导致其在升力最大时刻的低压区和高压区范围更大,翼上下表面的压力梯度也更大,因此使得“Click”运动产生更高的升力．

3.2.2 下拍过程的翼表面的压力云图与流线图 利用最大升力处的流线图可以探究简谐运动和“Click”运动在最大升力处前缘涡在翼表面呈现不同形态的原因,即“Click”运动高升力的涡量特性．

图11(a)、(b)分别是简谐运动和“Click”运动的压力云图与流线图．流线图较为清晰地展现了翼前缘涡量的形态．与图10给出的沿翼展方向的压力和涡量截面图的结论吻合,简谐运动的前缘涡在翼展的60%左右就已经出现脱离翼的趋势,而“Click”运动的前缘涡在翼展的80%之后才逐步脱离翼表面,与翼尖涡相融合．Liu等[4]在研究中发现,展向流动的存在稳定了前缘涡,使得鹰蛾在飞行过程中能维持不失速机制,从而产生高升力．在本研究中,无论简谐运动还是“Click”运动都存在一定范围内的展向流动,在翼尖前缘涡脱离部分未观察到展向流动的存在,由此可以推测展向流动和前缘涡的不脱落存在一定的联系．

(a) 简谐运动

(b) “Click”运动

图12表示两种运动在最大升力时刻及前后相同时间间隔的涡等值面形态变化图,图12(a)～(c)和(d)～(f)分别表示简谐运动和“Click”运动的涡结构变化,其中图12(b)、(e)分别表示简谐运动和“Click”运动在最大升力处的涡形态图．从图中可以看出,在两种运动的最大升力处,简谐运动的前缘涡在靠近翼尖部分已经逐渐破裂并脱离翼表面,而“Click”运动在最大升力处的前缘涡仍然保持完整．在图12(c)中,简谐运动的前缘涡已经完全破裂并脱离了翼表面;在图12(f)中,“Click”运动的前缘涡与翼尖涡相连,并逐渐脱离翼表面．在图12(a)中,简谐运动的前缘涡已经生长到完全形态;而在图12(d)中,“Click”运动的前缘涡还处于生长形态．图12的涡形态图与图11的流线图所展现的涡结构吻合较好．

图13(a)、(b)分别表示简谐运动和“Click”运动下拍过程中最大升力处的压力云图．两种运动在翼前缘中部的压力相较于翼根和翼尖低,可以推测展向流动的动力来自于展向的压力梯度．展向流动的存在能够将展向涡度的能量及时输运,从而限制了前缘涡的增长,保持了前缘涡的稳定[26-27]．如图9(a)所示,“Click”运动的最大升力落后于简谐运动,根据涡形态图分析,可能是展向

(a) τ=0.19(简谐运动)

(b) τ=0.29(简谐运动)

(c) τ=0.39(简谐运动)

(d) τ=0.25(“Click”运动)

(e) τ=0.35(“Click”运动)

(f) τ=0.45(“Click”运动)

(a) 简谐运动

(b) “Click”运动

压力梯度对展向涡度的能量输运限制了前缘涡的形成,从而延缓了最大升力到来的时刻．简谐运动的压力梯度和低压区范围明显小于“Click”运动,因此,前者的展向流动强度小于后者,则涡能量输运效率也小于后者．在此种条件下,简谐运动的前缘涡增长速率高于“Click”运动,前缘涡将会很快破裂并脱离翼表面．

3.3 上拍过程

本节将分析上拍过程中的不失速机制,即在上拍过程中前缘涡不脱落的现象,重点研究两种运动在最大升力处的流线图和沿着翼展方向截面的涡量云图,分析前缘涡在不同截面的形态．

图14(a)、(b)分别为简谐运动和“Click”运动的涡量云图．在上拍过程升力最大处,简谐运动和“Click”运动在涡量的耗散上区别较为明显．简谐运动下拍过程产生的涡量还未完全耗散,依然存在于翼前进方向的后部,而“Click”运动在下拍过程产生的涡量大部分已经耗散,其快速上拍的特点能够使翼快速地摆脱下拍过程所产生的尾涡干扰,保持较高的升力．

图15(a)、(b)分别为简谐运动和“Click”运动的压力云图与流线图．两种运动在上拍过程最大升力处的前缘涡的形态差异并不明显：在翼尖部分都产生了脱离翼尖前缘的趋势．从图15(a)可以看出,简谐运动仍然存在较多的尾涡,这也与图14(a)出现的较多尾涡现象相吻合．简谐运动的低压区范围明显大于“Click”运动的低压区范围,但是压力梯度却小于后者．

图14 上拍过程中不同截面处的涡量云图

Fig.14 Vortex nephograms at different cross sections during the upstroke

(a) 简谐运动

(b) “Click”运动

4 结 论

(1)相较于简谐运动,“Click”运动在拍动过程中能够保持更高的前缘涡稳定性．简谐运动在下拍升力最大处前缘涡从翼展的60%处开始破裂,而“Click”运动在翼展的80%处才出现破裂的状态．因此,“Click”运动能够产生更高的升力．

(2)扑翼运动能够保持前缘涡稳定的部分原因是,其沿着翼展方向的压力先减小后增大,翼中部的低压区使得沿着翼展方向存在着压力梯度．简谐运动的压力梯度较“Click”运动的压力梯度小,因此“Click”运动的展向流动更大,对展向涡度的能量输运效率就会更高,能够限制前缘涡的增长,保持更加完整的前缘涡形态,产生更高的非定常升力．

(3)在上拍过程中,两种运动的前缘涡形态相差不大,但“Click”运动的快速上拍能够让翼更快地脱离尾涡的干扰,以保持升力．