两自由度微扑翼飞行器升阻力系数分析

李德宝， 舒宁， 王洪洋

（合肥工业大学机械与汽车工程学院，合肥230009）

0 引言

微型扑翼飞行器由于运动方式独特、灵活性强等特点，越来越受到关注。在数值理论研究方面，研究人员［1］分析了二维机翼在一定频率范围内推力产生的机理，研究发现在0.7 Hz附近推进效果最佳；研究人员［2］分析研究了低雷诺数下，昆虫翅翼做非定常运动的气动性能，发现当翅翼大攻角扑动时，对大升力的产生十分有利；研究人员［3］通过数值模拟，对昆虫悬停时高升力机理有了更为深入的解释。计算流体力学（简称CFD）分析软件Fluent是通过数值模拟计算，对研究流体运动、热传递等现象具有良好的优势。研究人员［4］利用Fluent软件分析比较了翅翼在三种不同拍动模式下的升力和升力系数。研究人员［5］运用Fluent软件模拟分析了仿蜻蜓翅翼扑动时的流场特性。研究结果表明，对比以前风洞实验等研究方法，利用Fluent软件使空气动力的计算变得便捷和准确，且成本低、周期短，能分析模拟不同状态并获得完整数据。

利用Fluent软件，根据所设计两自由度扑翼结构的扑动运动规律，对两种不同扑动形式进行了二维仿真模拟，分析研究两种拍动形式下流场特性，获得翅翼在一个拍动周期内升阻力系数的变化情况，并进行相互间分析比较。

1 基本控制方程

在对昆虫和鸟类的长期观察中，研究人员发现翅膀上下扑动的过程中，不仅有上下平动，而且同时伴随着翅翼绕某一轴线的规律转动，可以认为是低雷诺数下的湍流流动。黏性不可压缩流体（二维）的基本控制方程由如下两个方程组成。

连续性方程为

式中：p 为压力，Pa；v为流体的速度矢量，m/s；Re 为雷诺数，常规飞行器流场雷诺数Re一般都大于106，由于扑翼飞行器在低空飞行，尺寸和飞行速度都较小，故流场雷诺数Re要比常规飞行器Re小得多，一般都在104～105以下。扑翼飞行器雷诺数按照Ellington［6］的定义，可利用翅膀的平均弦长D和平均翼尖速度V进行计算：

式中：ψ为展翅角；λ为展弦比；f为频率；l为翅长；μ为运动黏度。

2 扑翼结构分析

采用双曲柄双摇杆机构，故对曲柄摇杆机构进行简要分析。如图1所示，在机构中，构件L1为原动件，构件L3在构件L2带动下做往复运动，翅翼固连在构件L3上，从而实现翅翼的上下拍打运动。四杆尺寸参数如表1所示。

图1 铰链四杆机构

根据机构各杆所构成的封闭矢量多边形，可以得出以下矢量关系式：

求解可得角位移分量形式为：

其中 θ1、θ2、θ3分别为各构件的方位角，对于给定任意一个曲柄位置θ1，均可以根据上式求出连杆和摇杆的方位角 θ2、θ3，相应的角速度和角加速度可通过求一阶和二阶导数得到。

在微型电机驱动下，通过减速机构，构件1以14 Hz的频率匀速转动时，翅翼的上下扑动角度及角速度随时间变化规律如图2所示，翅翼的最大扑动角和最小扑动角分别为 124°、50.37°，整个拍动幅度约为 73°。

表1 四杆机构参数mm

图2 角度和角速度变化曲线

考虑到研究人员［7］所设计的转动机构过于复杂，故在凸轮机构的基础上引入了球形铰链，有效地简化了扑翼结构，减轻了整机重量。如图3所示，由微电机驱动，经过减速机构带动曲柄转动，通过连杆传递实现翅翼上下往复运动，同时由于翅根处凸轮机构的转向作用，从而使得翅翼从单一沉浮运动演变为沉浮加俯仰的两自由度运动模式，图4为凸轮扭转角度变化规律，如图所示，翅翼可以实现约为40°的扭转。

图3 三维实体模型

图4 凸轮扭转角度变化曲线

3 扑动过程中的升阻力系数分析

3.1 网格划分及动网格设置

网格划分之前，选取合理计算区域十分关键，合理的求解区域将直接影响整个计算结果的准确性以及计算时间。采用椭圆面为翅翼的弦向截面，从而建立计算模型。将翼型置于800 mm×400 mm的矩形区域内，为提高网格适应性、计算精度和节约计算时间，网格采用三角形非结构化网格。利用Gambit生成网格，如图5所示，共有30 396个节点数个，网格数为60 122个。

图5 全局网格图

利用CFD软件Fluent中的用户自定义函数（UDF），使用void DEFINE_CG_MOTION函数编译动网格UDF文件控制翅翼的上下拍动的运动规律，从而建立所需的动态网格模型。网格更新主要有三种动网格方式：弹簧近似光滑模型、动态分层模型和局部重构模型。弹簧近似光滑模型计算精度高，但在计算区域变形较大时，会影响计算精度；动态分层模型要求与运动边界相邻的网格必须为楔形或者六面体网格；局部重构模型仅能用于四面体网格和三角形网格。根据翅翼模型实际运动情况综合考虑，确定采用弹簧近似光滑模型和局部重构模型相结合的方法进行动网格的生成。

3.2 沉浮运动的升阻力系数

在FLUENT中有很多湍流模型可供选择，在扑动过程中，翅翼的上下往复运动使得翅翼附近的空气流动复杂多变，因此拟采用RNG k-ε模型。针对扑翼模型在流场的扑动飞行，设定无穷远处自由来流速度U∞=3 m/s，在控制体上利用simple算法［8］对其进行数值离散。翅翼在流场中快速拍打，翼背（上表面）和翼面（下表面）产生压力差，从而产生升阻力。升阻力系数计算公式为

式中：ρ为空气密度；c为翼型的弦长；U为周期内翼型表面来流平均速度。

在Fluent中，可以直接对升阻力系数进行监测，图6为Fluent软件求解流程图。

图6 Fluent求解流程图

设定拍动频率为14 Hz，时间步长为0.001 s，计算300个时间步长，迭代了6000次后，升阻力系数变化规律曲线如图7所示。

图7 升阻力系数变化曲线

从图中可以看出在下拍过程中，升力系数急剧增加，约在0.25T时达到了峰值，此时在上表面形成的漩涡最大，上下表面的空气流动速度的差异使得上下表面形成了压力差，上表面压力小于下表面，便形成了升力。同时由于攻角存在，后缘涡强度大于前缘涡强度，表现为阻力系数的增加。上拍过程中升阻力系数变为负值。在整个拍动周期，净升阻力系数为正值。下拍时流场速度及压力分布如图8～图9所示。

图8 下拍时流场速度分布

图9 下拍时流场压力分布

3.3 沉浮加俯仰运动的升阻力系数比较

在凸轮机构的转向作用下，可以实现翅翼约为40°的扭转。翅翼在下拍过程中与来流方向呈10°夹角，到最低位置时经过凸轮机构的作用，使得翅翼在上拍过程与来流方向夹角为50°。可根据仿真数据将扭转变化规律简化为如下方程：

运用与上述相同的研究方法，将转动的运动模式加入翅翼的运动中，可以得到升阻力系数的变化曲线，如图10所示。图中，A为沉浮运动时的升阻力系数，B为沉浮加俯仰运动时的升阻力系数，通过曲线对比可知，沉浮加俯仰运动时，上拍过程中，升力系数明显升高，同时阻力系数也有所降低，整个拍动周期翅翼获得的净升力系数要大于单一的沉浮运动时的净升力系数，净阻力系数要小于单一沉浮运动时的净阻力系数。上拍时流场速度与压力分布如图11～图12所示。

图10 升阻力系数对比图

图11 上拍时流场速度分布

图12 上拍时压力流场分布

4 结论

CFD软件Fluent可有效模拟模型翼在流场中的运动情况，所提供的自定义函数（UDF）窗口把翅翼相对复杂的运动模式转化成可编译的外部C程序，有效解决了数值模拟中流场N-S方程所存在的“动边界问题”，结果表明运用Fluent软件模拟分析流场相关特性是可靠的。

研究结果表明，翅翼在不同的运动模式下的升阻力系数存在着差异，由于凸轮转向机构的作用，两自由度的机构在扑动中获得了更大的净升力系数，同时也降低了净阻力系数，有效增大了升力，减小了阻力。

研究局限在机构有无扭转情况下的升阻力系数，且均采用二维简化模型，对于具体翼型产生的影响及更复杂的扑翼运动规律还需更深入的研究。

［1］ Wang J.Vortex shedding and frequency selection in flapping flight［J］.JFluid Meeh，2000，410：323-341.

［2］ 兰世隆，孙茂.模型昆虫翼非定常运动时的空气动力特性［J］.力学学报，2001，33（2）：173-182.

［3］ Sun M，Hossein H.Force and flow structures of an airfoil performing some unsteady motions at small Reynolds number［J］.Aeta Aerodyamica Siniea，2000，18（Suppl）：96-102.

［4］ 朱洪波，谢进，陈永.用Fluent分析仿生翅翼运动产生的升力和升力系数［J］.机械设计与制造，2011(5)：68-70.

［5］ 张孝松，基于蜻蜓翅翼的仿生微扑翼飞行器机翼的有限元分析［D］：哈尔滨：哈尔滨工业大学，2006.

［6］ E1lington C.The aerodynamics of hovering insect flight∶I.the quasi-steady analysis ［J］.Philosophical Transactions of the Royal Society of London B，1984,305∶1-15.

［7］ 朱宝，王姝歆.两自由度扑翼机构及其运动仿真研究［J］.中国制造业信息化，2009,38（21）：24-28.

［8］ 周建华.微型拍翅式飞行机器人翅运动及控制系统研究［D］.南京：东南大学，2005.