1 引言
机翼扑动飞行的仿生飞行器具有高度仿生性和良好隐蔽性能[1]。近年来，许多受飞行生物启发的仿生扑翼机器人平台被研制出来[2]，昆虫级微型扑翼机有苍蝇[3-4]，delfly[2,5]等。鸟类作为常见的飞行动物也受到人们关注[6-7]。机翼可变形的机器鸟是仿生扑翼飞行器研究的一个活跃领域[8]。同时为提高扑翼机适应复杂场景稳定飞行的能力，对仿生扑翼飞行器的气动性能研究也显得尤为重要。
对于仿鸟扑翼飞行器研究，德国Festo公司仿照海鸥飞行姿态研制出的Smartbird[9-10]成功飞行，区别于以往的扑翼飞行器的机翼骨架结构，Smartbird将机翼分离成内段翼和外段翼，使得机翼运动过程中可做两段翅膀的运动形态不同变化。提出了内段翼提供主要飞行升力，外段翼提供主要飞行推力理论。Lentink团队研发出一种带有真实鸟类羽毛并成功飞行的飞行器PigeonBot[8],通过风洞实验和外场空中飞行，发现翅膀不对称变形以及羽毛位置不同可以控制飞行方向，实现仿生飞行器在复杂环境中灵活飞行。昂海松团队[11-12]研究与设计多段柔性变体扑翼飞行器，探讨了扑翼流场涡结构，总结升推力产生机理。杨永刚等[13]采用近似柔段翼思想建立扑翼机模型，并对此模型进行气动仿真实验，为柔性翼多自由度扑翼机研制提供理论基础。陈亮等[14]通过优化的叶素理论，对小型仿鸟飞行器在非常定流场条件下的气动特性进行研究。
本文分析鸟类基本参数以及飞行姿态特征，设计具有“拍打-折叠”扑翼机构，建立三维折叠翼扑翼机模型，研制可折叠“铰链式机翼”，实现机翼翼尖“8”字形运动轨迹[15]，接近真实鸟类翅膀扑打状态。利用CFD软件对折叠翼扑翼机进行了低雷诺数非常定流场情况下气动仿真模拟，分析折叠翼飞行器初始迎角，机翼展弦比和扑动频率等对气动特性作用。最终研制“羽毛版”可折叠机翼，此飞行器以扑动羽毛机翼产生动力，并成功起飞并持续飞行90多秒后安全着陆。
2 扑翼机机构建模仿真
2.1 鸟类扑翼运动机理
文献[12]通过对绿头鸭飞行动作进行研究，将鸟类飞行姿态分为以下阶段：
下扑阶段：翅膀外翼翼尖从最高点运动到最低点的过程，全过程翅膀处于全展开状态。
折弯阶段：鸟类翅膀下降至最低处，外段翼相对内段翼向下弯折。
上扑阶段：鸟类翅膀外翼翼尖从最低点向上运动过程，整个过程鸟类整体成“M”字型，此时翅膀处于收缩状态。
展平阶段：鸟类翅膀外翼相对内翼做舒展动作，外翼翼尖运动至最高点，整个翅膀全展开，准备进入接下来的下扑阶段。
2.2 扑翼机建模
折叠翼扑翼机运动模型参考海鸥外形，并在其翅膀复杂运动的基础上进行简化，在Solidworks平台上进行建模，设计如图1所示骨架模型。该骨架模型包括主机身，尾翼，内段翼和外段翼等，如图1(b)。鸟身(包括尾翼)长1 045 mm，单边翼展长约900 mm，其中内段翼长320 mm，外段翼长约 580 mm，机翼平均几何弦长280 mm。
图1 骨架模型示意图
Fig.1 The diagram of skeleton model
2.3 扑翼机传动机构设计
扑翼机两翼间扑打频率，幅度，变化相位在保持飞行器飞行稳定方面起到至关重要作用。在两翼对称拍打过程中，忽略环境因素的影响，受到的空气对机翼产生的力矩是相当等的。反之，两翼间不对称扑动，同时刻空气对两翼做功不同，产生滚转力矩，导致扑翼机无法平稳飞行。
采用双曲柄双摇杆机构，平面对称传动模式，如图2所示。优势在于确保左右机翼运动保持一致，实现机翼对称扑打。单边曲柄使用4个M2沉头螺钉将其固定在主齿轮，机翼通过曲杆与曲柄相连。电机转动，通过减速比i=1∶48的减速齿轮组，实现2个主齿轮相对圆周运动，带动两边机翼的同时上下扑打运动。
图2 双曲柄双摇杆机构简图
Fig.2 The diagram of double crank double rocker mechanism
2.4 机翼运动方程
分解鸟类翅膀复杂的运动形态，建立机翼模型，可实现鸟类飞行过程中扑打-折叠动作。针对机翼运动模式，建立如图3所示模型简图，点G为齿轮圆心，GA表示曲柄，AC为曲杆，长度已知；四边形BCEF为机翼的内段翼，其中点D为固定支撑点；EH表示外段翼，通过铰链结构与内段翼相连。机翼内段翼与外段翼运动轨迹方程，主要体现为E,H两点。本文模型中所有杆子长度单位均为毫米。
图3 机翼模型简图
Fig.3 The diagram of wing model
取x轴设定为水平地面。设转动角速度为ω0rad/s，求得在任意时刻，曲柄GA转动角度：
θ0=ω0t
(1)
其中θ0∈[0,2π]；θ1为图3所示辅助线lDA与水平地面的夹角，θ1∈[0,π]；θ2为内段翼运动过程中相对水平地面夹角，
为外段翼相对水平地面折叠角，
为曲杆BA运动相对水平地面夹角。 β为曲杆BA运动相对水平地面夹角。
扑翼机内段翼运动，求解对E点运动轨迹：
(2)
(3)
几何关系得出：
β=θ1+∠CADZ
(4)
根据几何关系可知，D为骨架上固定点，坐标为(xD,yD)；G点坐标(lGA,0)，则点A坐标表示为
(5)
根据三角形余弦定理：
(6)
通过固定点D与动点A之间几何关系：
(7)
联立等式(2)—式(7)，最终可求解出E点的具体运动轨迹，确定内段翼运动模式。同理，求解出点H的坐标，将周期内H的坐标相连，确立外段翼运动轨迹。
四边形BCEF为平行四边形，对H点：
(8)
θ3=∠FEI=90°-β
(9)
(10)
联立式(8)—式(10)，可求得H点运动轨迹。其中：lDE为固定点D至内段翼最外侧点间距离，lCD为内段翼最内侧点至定点D间距离，lAC为曲杆长度，lGA为曲柄长度，lEH为外段翼长度。
2.5 扑翼机运动仿真
根据2.4节推导出的扑翼机内段翼与外段翼运动轨迹方程，对设计的机构进行运动学求解。根据给定折叠翼机构设计参数：翼展长度，内、外段翼长度，弦长参数。在机械运动约束条件下：曲柄长度：25 mm≤lGA≤30 mm；曲杆长度：45 mm≤lAC≤50 mm；摇杆长度：lCD≤lCE；固定支撑点D与曲杆C点距离：40 mm≤lCD≤50 mm。实际飞行实验期望机翼扑动幅度50°≤θ2≤70°；扑翼机扑打频率f ≥ 3Hz。在Solidworks平台上对设计的扑翼机模型进行基本运动仿真及motion 分析。最终导出机翼运动轨迹，如图4所示。
图4 机翼运动轨迹图
Fig.4 Wing movement trajectories
机翼运动仿真路径，内段翼运动做相对简单上下扑打运动，外段翼相对内段翼具有折叠动作，轨迹曲线接近“8”字形运动，成功模拟真实鸟类扑翼运动的规律。
3 扑翼机气动特性研究
3.1 三维气动仿真平台搭建
为了达到在运动状态下，研究扑翼机机翼初始迎角，展弦比及机翼拍打频率对气动特性的影响，本文对23组不同机翼初始迎角，7组不同展弦比及5组不同拍打频率对折叠翼飞行器进行三维气动仿真实验和数据分析。选取机翼初始迎角变化范围在-10°～12°；模拟扑翼机展弦比为：1、1.5、2、2.5、3、3.2、4；模拟扑翼机机翼扑动频率：2 Hz、3 Hz、4 Hz、5 Hz、6 Hz。在软件平台XFlow中选用虚拟风洞模式对折叠翼飞行进行气动仿真，虚拟风洞基本设置参数如表1所示。
表1 虚拟风洞基本设置参数
Table 1 Parameter of virtual wing tunnel
名称说明参数计算区域虚拟风洞长方体4 m×2.5 m×3 m流体空气密度ρ=1.225 kg/m3温度T=288.15 K动力粘度系数μ=1.789×10-5 Pa·s入口来流速度VX=5 m/s出口压力边界相对入口压力Pout=0 Pa机身模型沿X轴正向沿各轴方向速度为V=0 m/s
3.2 扑翼机机翼初始迎角对气动特性影响
本文以NACA7412翼型为计算模型，探讨机翼初始迎角对气动特性影响。以地面坐标系为基准，将机翼未开始运动，处于平放状态下具有的迎角设定为初始迎角。通常情况采用无量纲数升力系数Cl，阻力系数Cd衡量飞行器运动中受到的升力及阻力。在XFlow平台下模拟不同机翼初始迎角下的Cl，Cd数据，并进行数据处理，得到如图5所示的曲线。选取初始迎角在-10°～12°，变化步长Δ=1°，约束条件：来流风速为5 m/s，雷诺数Re设定为1.6×105。
图5 不同初始迎角下升力、阻力系数(a)与
升阻比(b)变化曲线
Fig.5 Curve of lift coefficient,drag coefficient (a) and lift-drag ratio at different initial wing angle(b)
由图5(a)可以分析出，升力系数与机翼初始迎角成正比例关系。在机翼初始迎角为负值时，升力系数为负值，此时机翼扑动主要产生负升力，随着机翼初始角度增加，升力系数变为正值，并增大。初始迎角到达-7°时，升力系数变为正值，此时机翼扑动开始产生正向升力；阻力系数呈现随初始角度增加先逐渐减小后增大的趋势，图5(b)中升阻比变化曲线存在多个拐点，体现为随初始角度增大先增大后减小，再增大最后逐渐减小过程，同时在-5°后升阻比急剧增大，5°时升阻比达到峰值，随后曲线变得平缓。有关数据如表2所示，取-10°、-5°、0°、1°、5°、10°不同初始迎角做比较。在1°时受到阻力最小，通过改变飞行器平飞阶段翼面迎角，在该角度下受到最小空气阻力；当扑翼机初始迎角为5°飞行时，获得最大升阻比，在起飞阶段，有助于快速升空。
表2 不同初始迎角下平均气动特性参数
Table 2 Average aerodynamic characteristic parameters at different initial wing angle
初始迎角α / (°)平均升力系Cl平均阻力系Cd升阻比K-10-0.185 80.127 0-1.463 6-50.173 70.053 973.218 500.785 60.015 3851.079 310.891 30.013 7964.633 851.313 20.018 4871.060 6101.743 70.025 4168.622 6
3.3 机翼展弦比对扑翼飞行器气动特性影响
简化仿真模型，将内段翼与外段翼简化成一段机翼，取扑翼机左机翼为参照，改变机翼翼展，研究扑翼机展弦比为1、1.5、2、2.5、3、3.2、4对扑翼运动的影响。仿真时间为2 s，约束条件：来流速度为5 m/s，机翼扑打角度θ=60°，机翼扑动频率为4 Hz，即一个运动周期T=0.25 s，初始迎角为5°；机翼扑打满足急回特性规律，即机翼一个运动周期内，向下扑动时间大于上扑时间，并根据文献[12]指出扑翼运动近似正弦或者余弦变化规律。在XFlow中设定产生时间差的分段函数。扑翼机机翼一个周期内扑打角度变化表达式如下：
(11)
其中时间t′可通过求余运算，如式(12)所示，实现机翼连续周期性运动。
t′=t-⎣t/T」×T
(12)
式中： t为仿真实验运行时刻，T为运动周期。
φ为仿真实验未开始时，机翼初始位置与水平面夹角，调整φ值可改变机翼初始位置，本文取φ=0°；α为机翼翼尖扑动至最高(低)点，机翼平面与机翼初始位置平面夹角，本文取α=-30°。联立式(11)、式(12)，求解出机翼扑动角度θ变化，如图6所示。一个周期内，当角度变化斜率为正值，机翼为下扑，时间为0.15 s，角度变化斜率为负值，机翼上扑，时间为0.1 s。
图6 机翼扑动角变化曲线
Fig.6 Curve of flapping wing angle
一个运动周期内，取0.5～0.75 s进行气动分析，不同展弦比下扑翼机升力系数和阻力系数变化曲线如图7。
图7 不同展弦比下升力系数(a)与阻力系数(b)变化曲线
Fig.7 Curve of lift coefficient (a) and drag coefficient (b) at different aspect ratio
由图7(a)升力系数变化曲线结果可知，不同展弦比机翼的升力系数变化曲线大体相似。0.5～0.65 s为下扑阶段，升力系数曲线先增大后减小，在机翼运动值中点位置时，升力系数达到峰值，由于急回特性规律约束，此时机翼速度最快，受到空气反作用力最大，即产生升力最大，此后机翼减速下扑，升力逐渐减小；0.65～0.75 s上扑阶段，负升力逐渐增大后减小，上扑至中点位置达到负升力峰值。
图7(b)阻力系数变化曲线说明，不同展弦比的阻力系数在下扑阶段先增大后减小，大展弦比阻力系数的峰值在更早时间出现，在下扑阶段展弦比AR≥2.5的机翼阻力系数峰值在0.57 s，当展弦比AR≤2时，阻力系数峰值出现在0.61 s左右。上扑阶段机翼阻力系数先减小后增大，展弦比较小机翼变化不明显，随着展弦比增大，阻力系数最小峰值逐渐增大，此时正推力逐渐增大。
不同展弦比机翼对扑翼机气动特性影响如表3所示。在来流速度5 m/s，扑打频率为4 Hz的情况下，机翼的升力系数与阻力系数随着翼展增大而增大，升阻比先减小后增大。扑翼机升阻比在展弦比AR=2.5时最低，随后增大。
表3 不同展弦比下平均气动特性参数
Table 3 Average aerodynamic characteristic parameters at different aspect ratio
展弦比AR平均升力系数Cl平均阻力系数Cd升阻比K10.167 10.062 32.6821.50.165 50.081 42.03320.212 70.154 51.3772.50.285 10.238 51.19530.448 50.318 31.4093.20.634 00.357 11.77541.440 60.528 02.728
在实际飞行中，应当考虑机翼负载能力以及机翼结构强度问题。其中负载能力以飞行器净升力表示[16]，净升力变化曲线如图8所示。
图8 扑翼机净升力随展弦比变化曲线
Fig.8 Curve of net lift changing with aspect ratio
仿真结果分析考虑折叠翼扑翼机整体受到重力影响。由图8可知，在展弦比AR≤2时，扑翼机整体净升力为负值，此翼展下机翼无法带动整机升空，随着翼展增大，展弦比在2到3.2之间净升力从0 N提升至3 N，随后由于机翼增长，重量随之增加，净升力提升维度减缓。
图9(a)等所示为不同展弦比的扑翼机下扑阶段机翼压强分布图，取展弦比AR=1，3，4做比较。左侧为机翼上表面压强分布情况，右侧为同时刻机翼下表面受到压强，高压区域主要集中机翼最外侧部分，说明机翼扑动过程之中，最外侧受到空气反作用力最大，相应的，机翼结构越易受损，同时长翼展机翼运动消耗的能量更多。综合考虑材料、能源消耗、扑翼机整体重量以及气动效率，机翼展弦比AR=3左右可提供良好飞行性能，同时实现更为经济。
图9 不同展弦比下机翼压强分布云图
Fig.9 Static pressure diagram of wing at different aspect ratio
3.4 折叠翼扑动频率对飞行器气动特性影响
图10(a)和图10(b)分别表示不同的机翼扑动频率下升力系数和阻力系数曲线。实验条件为来流速度为5 m/s，机翼初始迎角为5°，机翼展弦比为3，选取仿真时间1 s。
随着扑动频率增大，机翼上下表面压强差值增大，扑翼机升力系数和阻力系数的波峰正向增大，波谷负向增加。由图10(b)分析得出：在机翼扑动频率4～6 Hz时，阻力系数负向峰值出现负数，即正向推力系数。不同扑动频率的平均升力系数与平均阻力系数如表4所示，结合图10和表4数据，扑翼机扑动频率在4～5 Hz获得较优气动效率。
图11(a)-图11(e)表示扑翼机机翼扑动频率逐渐从2 Hz增大至6 Hz，机翼附近空气流动变化。左侧为机翼下扑至最低点，右侧为机翼上扑，上扑过程中，外段翼相对内段翼做折叠运动。随着频率的增加，机翼附近的空气流速增大；下扑时，机翼上表面空气流速明显大于下表面，产生压差，使得飞行器获得升力，机翼折叠上扑时，由速度云图变化来看，机翼下表面空气流速大于上表面；折叠翼不对称运动，一个扑打周期内，空气对其做功不均等，使其飞行过程中获得更大的净升力。
图10 不同扑动频率下升力系数(a)与
阻力系数(b)变化曲线
Fig.10 Curve of lift coefficient (a) and drag coefficient (b)
at different flutter frequency
表4 不同扑动频率下平均气动特性参数
Table 4 Average aerodynamic characteristic parameters at different flutter frequency
扑动频率f /Hz平均升力系数 Cl平均阻力系数 Cd升阻比K20.419 30.280 11.49730.638 80.326 41.95740.797 60.373 02.13850.892 10.408 62.18360.907 70.430 72.107
图11 不同扑动频率下速度云图
Fig.11 Velocity contours at different flutter frequency
4 羽翼动力飞行器飞行实验
鸟类翅膀不仅实现拍打-折叠2个自由度运动，还可实现翅膀合拢-展开，非对称伸展等多自由度运动。文献[8]指出：鸟类控制翅膀羽毛不对称变形从而快速转向，实现复杂环境高效飞行。传统扑翼机采用膜状翼结构，如图12(a)所示，机翼覆膜区域连为一体，无法实现飞行器机翼变形。本文仿照真实鸟类翅膀分层结构，研制出羽毛机翼，如图12(b)所示，实现拍打-折叠2个自由度运动。
图12 机翼结构模型示意图
Fig.12 Model diagram of wing structure
羽毛机翼优势在于羽毛分层结构及羽毛独立性。在扑翼机未来研制中，控制每根羽毛的独立运动，实现机翼多自由度变形，使得人造机器鸟具有更强仿生性能，提高复杂环境飞行适应性。
本文研制羽毛机翼实现拍打-折叠2个自由度运动。扑翼机不同阶段姿态如图13(a)等所示。首次使用羽毛扑打产生动力实现人造机器鸟飞行[17]。扑翼飞行器扑打频率约4 Hz，展弦比为3.2，飞行时长90 s，证明仿真结果准确性，同时为研制机翼多自由度变形仿生机器鸟提供理论可行性和有效飞行平台。
图13 全羽翼折叠翼扑翼机外场飞行实验姿态图
Fig.13 Field flight experiment of full-wing folding wing flapper
5 结论
本文仿照真实鸟类飞行机翼运动机理，建立“铰链式”折叠翼扑翼机模型，在传统扑翼飞行器机翼基础上，通过铰链键将机翼分为内段翼和外段翼，并通过两段机翼间运动相位差，实现模仿鸟类扑打-折叠动作。通过模型运动仿真以及气动仿真实验，得出以下结论：
1) 建立折叠翼飞行器模型，实现内段翼往复运动，外段翼翼尖“8”字运动模式，成功模拟真实鸟类翅膀运动轨迹。证明了模型的可行性。
2) 仿真实验表明：扑翼机机翼初始迎角在5°时获得最大升阻比系数；扑翼飞行器机翼展弦比为3左右气动性能较优；提高机翼扑动频率可提升扑翼机的气动效率，并在4～5 Hz时取得最佳气动效率。
3) 羽翼动力扑翼飞行器外场成功飞行，验证了飞行器模型和仿真实验结果合理性、准确性，为后续仿生机器鸟复杂运动提供有效的实验数据。
[1] 昂海松.微型飞行器设计导论[M].西安.西北工业大学出版社,2012.
[2] Matěj,Karásek,Florian T,et al.A tailless aerial robotic flapper reveals that flies use torque coupling in rapid banked turns[J].Science,2018,361(6407):1089-1094.
[3] Ma K Y,Chirarattananon P,Fuller S B,et al.Controlled flight of a biologically inspired,insect-scale robot[J].Science,2013,340(6132):603-607.
[4] Wood R J.The First Takeoff of a Biologically Inspired At-Scale Robotic Insect[J].IEEE Transactions on Robotics,2008,24(2):341-347.
[5] De Croon G C H E,De Clercq K M E,Ruijsink,et al.Design,aerodynamics,and vision-based control of the DelFly[J].International Journal of Micro Air Vehicles,2009,1(2):71-97.
[6] Pennycuick C J.Modelling the Flying Bird[M].Academic Press,2013.
[7] Weimerskirch H,Bishop C,Jeanniard-Du-Dot T,et al.Frigate birds track atmospheric conditions over months-long transoceanic flights[J].Science,2016,353(6294):74.
[8] Chang E,Matloff L Y,Stowers A K,et al.Soft biohybrid morphing wings with feathers underactuated by wrist and finger motion[J].Science Robotics,2020,5(38):doi:10.1126/scirobotics.aay1246.
[9] Send W,Fischer M,Jebens K,et al.Artificial hinged-wing bird with active torsion and partially linear kinematics[C]//Proceeding of 28th Congress of the International Council of the Aeronautical Sciences.Brisbane,2012.
[10] Mackenzie D.A Flapping of Wings.[J].Science,2012,335(6075):1430-1433.
[11] 肖天航,段文博,昂海松.仿鸟柔性扑翼气动特性与能耗的数值研究[J].空气动力学学报,2011,29(6):709-718.
Xiao T H,Duan W B,Ang H S.Numerical investigation on aerodynamics and power requirement of a bird-likeflexible flapping wing[J].Acta Aerodynamica Sinica,2011,29(6):709-718.
[12] 曾锐,昂海松.仿鸟复合振动的扑翼气动分析[J].南京航空航天大学学报,2003,35(1):6-12.
Zeng R,Ang H S.Aerodynamic Computation of Flapping-Wing Simulating Bird Wings[J].Journal of Nanjing University of Aeronautics & Astronautics,2003,35(1):6-12.
[13] 杨永刚,顾新冬.多段翼多自由度扑翼飞行器气动特性仿真研究[J].系统仿真学报,2019,31(09):1827-1834.
Yang Y G,Gu X D.Simulation Study on Aerodynamic Characteristics of Flapping-wing Aircraft with Multi-stage Wing and Multiple Degree of Freedom[J].Journal of System Simulation,2019,31(09):1827-1834.
[14] 陈亮,管贻生,张宪民.仿鸟扑翼机器人气动力建模与分析[J].华南理工大学学报(自然科学版),2011,39(06):53-57,70.
Chen L,Guan Y S,Zhang X M.Aerodynamic Modeling and Analysis of Biomimetic Flapping-Wing Robot[J].Journal of South China University of Technology(Natural Science Edition), 2011,39(06):53-57,70.
[15] 姜森,郝永平,李伦,等.基于空间连杆机构实现“8”字形运动的扑翼机的设计[J].机床与液压,2019,47(13):76-80.
Jiang S,Hao Y P,Li L,et al.Design of Flapping Machine for the“8”Zigzag Movement Based on Spatial Linkage Mechanism[J].Machine Tool & Hydraulics,2019,47(13):76-80.
[16] 张志君,陈默,杨贺捷,等.基于XFlow的仿生扑翼飞行器机翼气动特性分析[J].东北大学学报(自然科学版),2021,42(06):821-828.
Zhang Z J,Chen M,Yang H J,et al.Analysis of Aerodynamic Characteristics of Bionic Flapping Wing Aircraft Based on XFlow[J].Journal of Northeastern University(Natural Science Edition),2021,42(06):821-828.
[17] Mattehw S.Eagle-inspired robot flies by flapping its feather-covered wings[EB/OL].(2021-03-03)[2021-07-12].https://www.newscientist.com/article/2269774-eagle-inspired-robot-flies-by-flapping-its-feather-covered-wings/.