鸽状扑翼飞行器气动特性研究

刘旭博，于纪言

(南京理工大学 智能弹药国防重点实验室, 南京 210094)

扑翼飞行器是一种模拟鸟类飞行的仿生飞行器，在无人机发展中占有重要位置且应用需求巨大[1]，而相比于固定翼飞行器和旋翼飞行器，扑翼飞行器具有更高的飞行效率与机动性，因此扑翼飞行器也是近年来的研究热点。

国外研究扑翼机起步较早，最早在1909年Knoller和Betz[2]发现二维翼型在扑动过程中可以同时产生升力和推力。1975年Weis-Fogh[3]解释了昆虫升力产生的拍合机制。1997年，Van Den Berg等[4]在对飞蛾飞行的研究中发现了“延时失速”机制，并在后来被Dickinso等[5]在果蝇模拟实验中证实。近十年来，国外渐渐入手将理论转化到实践，较为典型的是在珠海航展上露面的smart bird，该扑翼机是由德国Festo公司在研究银鸥的基础上仿生制作的，是现代人类较为成功的仿生扑翼机。

国内对扑翼机的研究起步较晚，集中于近20年，但也取得了不少成就。曾锐与昂海松[6]对扑翼的扑动模型做了相关研究并指出了扑动频率与扑动幅度对平均升力系数的影响。王姝歆[8]着重研究了仿生翅翼的中的尺度率问题。孙茂等[9]对蝴蝶与蜜蜂的飞行力学进行了相关的理论研究。另外，我国对扑翼的研究主要体现在实验分析上，2006年至今，多位学者都对针对扑翼飞行器进行了相关的风洞实验研究，并与不少关于扑翼机的理论研究相印证[10-13]。

尽管人类对扑翼机研究已久，然而现代扑翼机的飞行效率仍然远远小于鸟类的飞行效率，要想达到自然界鸟类的飞行效率，显然仍然需要更多的研究积累。本文在前人对鸟类扑翼结构研究的基础上，使用计算流体力学方法，采用三维鸽状有限翼展扑翅模型研究翅翼拍动时的气动特性，并对比不同攻角下的翅翼升阻力情况。

1 扑动函数与扑翼建模

1.1 扑动函数

在仿真之前需要先得到扑翼的扑动规律，扑翼的结构研究在近年来已趋于成熟，其中公认的易于实现、扑动过程对称且振动较小的扑翼机构为双曲柄连杆机构，如图1：

图1 双曲柄连杆机构

这里为了简化模型，设置上下扑动幅值在45°左右，且扑动对称。由于扑动对称，机构没有急回特性，因此机构的扑动角度曲线可以近似成一个简谐运动的函数曲线，另外根据文献[8]，鸟类的扑动角度变化规律确实类似于一个三角函数的周期性变化。

扑翼角度函数的格式如下：

θ=bsin(2πft)

(1)

其中:θ为翅翼与水平面的夹角，b幅值，即翅翼扑动时翅翼能达到的极限位置，f为扑翅的频率。t为扑动时间。

由于扑翼角度函数控制着翅翼的扑动规律，因此各参数对气动特性有着明显的影响。根据文献6和文献7，一个翼型的扑动幅值b越大，其周期内平均升力系数就越大，在0～60°范围内基本呈线性相关。扑动频率f对气动特性的影响和翅翼翼型有关，基本规律是：升力系数在扑动频率的某一点会达到峰值，在峰值前升力系数随着扑动频率的上升而增大，在峰值后升力系数随着扑动频率的上升而减小。

以鸽子为扑翼飞行器的仿生对象，通过观察鸽子飞行的高速摄影录像可知，鸽子的扑动频率在6～10 Hz左右，具体设置函数时，采用8 Hz的扑动频率。以θ表示翅翼与水平面的夹角，时间t为自变量，初始位置时翅翼水平，则扑翼角度函数如下：

(2)

求导得到其角速度函数：

(3)

1.2 扑翼建模

该扑翼机翅翼以鸽子为仿生对象，典型鸽子的飞行参数与翅翼特征参数数据如表1。

表1 典型鸽子参数

根据以上参数与具体鸽子的翅翼形状(如图2)，进行翅翼的建模：

图2 鸽翅

根据鸽翅的生物学构造，鸽翅的前缘部分主要由肱骨支撑，量得肱骨直径约为0.01 m。因此仿生翅的前缘半径选择为0.005 m。翅翼的后缘部分为羽毛尖端，厚度可以忽略不计，因此构建模型时后缘厚度为0。鸽翅的厚度最大处位于肱骨与副羽的相接部分，约在弦向20%的位置，厚度为0.016 m。翼展前缘长度约0.13 m，后缘长度为0.3 m。鸽翅的下表面展开时没有明显的弧度，为平面型下表面。鸽翅的翅翼模型如图3，仿生翅模型的具体参数如表2所示。

图3 翅翼模型

参数符号数值弦长c0.1 m展长b0.3 m前缘半径r0.005 m最大厚度t0.016 m最大厚度位置p20%前缘长度q0.13 m后缘长度h0.3 m展弦比λ3

对翅翼建模完成后，构建外流域，外流域大小与流域流体的流速有关，由于是低速飞行，雷诺数较小，流场变化范围不大，因此取外流域的长宽大约是翅翼的6倍即可。构建后的外流域如图4：

图4 计算域

2 计算前处理

1) 计算方法及网格划分策略

Fluent在进行流场计算时，主要采用了有限体积法。有限体积法的基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。

在CFD计算过程中，翅翼形状相对复杂，且翅翼是不断扑扇的，因此这里采用四面体网格，四面体网格有较好的复杂模型适应性，因此，本文采用四面体网格。在大尺度的变形下，计算迭代中会重构运动边界附近的网格，从而适应部件新位置的流场运算。使用商业软件ICEM CFD对：流动区域进行离散，最终生成的网格质量大于0.3，如图5所示。

图5 计算域网格

2) 计算模型

由于研究的是扑翼运动，气体流速较慢，气体密度基本无变化，采用常数值。流场以湍流为主，采用低雷诺数下的k-e方程进行计算。计算时主要使用的方程如下：

质量守恒方程：

(4)

式中：ρ为气体密度，x,y,z为坐标，u,v,w为速度矢量，t为时间。

动量方程：

(5)

低雷诺数下的k-ε方程：

Gk+Gb-ρε-YM+Sk

(6)

(7)

其中ρ为流体密度，k为湍动能，ε为湍流耗散率，Gk是层流速度梯度而产生的湍流动能，Gb是由浮力而产生的湍流动能，YM是在可压缩湍流中由于过渡的扩散而产生的波动。αk和αε是两个方程的普朗特数，Sk和Sε是在计算中自定义的参数。另外其中的数值C1、C2、C3具体如下：

(8)

C2=1.9

C3=1.0

3) 边界条件与动网格

为了模仿鸟类飞行时空气流向，设置外流域正面及周边四个面为速度入口边界，流速为11 m/s，方向为x方向。设置翅翼后方的一个面为压力出口流场，在有一个大气压的环境压力情况下，压力出口设置为0。翅翼部分默认为wall边界。

动网格是本文研究中的一个重要手段，将1.1节中得到的角速度函数写入UDF文件，而后在Fluent中编译。由于网格变形较大，则同时采用动网格中的Remeshing和Smoothing的方法对网格进行更新。在最终选择时间步长时应该考虑到网格特征尺寸与流速的影响，否则容易遇到更新网格失败的情况，从而导致计算中断。最终时间步长大小应该基本遵守以下公式：

(9)

其中:d为网格特征尺寸，v为流体流速。在以上公式计算的基础上也可以适当调整从而节省时间与资源。在本文的计算过程中，时间步长取0.000 05 s。

3 鸽状扑翅气动特性的数值计算

3.1 翅翼拍动升、阻力变化特点

在Fluent中设置监测时，一般得到的是升阻力系数Cl、Cd。无量纲量升阻力系数就可以直接表征升阻力的大小。升阻力系数的定义如下：

(10)

(11)

其中:L为升力，D为阻力，ρ为气体密度，v为气体特征流速。

Fluent中设置气体来流速度11 m/s，设置攻角为8°，可以监测到升阻力系数在周期内的变化情况如下：

由图6、图7可以看出，随着翅翼的周期性扑动，升力与阻力也随着时间出现周期性的变化，在翅膀下扑的过程中，当翅翼在水平位置时，翅翼的面积在竖直方向上的投影达到最大值，此时竖直方向上，即升力为最大值，当翅翼向上扑动复位时，翅翼与气流的受力是向下的，因此产生了向下的力，则升力为负值。同理当翅翼的扑扇速度达到最大值时，翅翼与气流的相互作用最强，此时在翅翼上的压力也达到最大值，因此即使翅翼在竖直面上的投影面积不变，但由于作用于翅翼上的压力的变化，阻力也因此发生相关的变化。绘制不同时刻下翅翼上下表面的压力云图，如图8所示

从图8可以看出，当翅翼扑翅速度在最大值附近时(t=0.005 s)，翅翼下翼面压力远远大于上翼面，因此产生较大的升力。当翅翼接近极限位置时(t=0.030 s)，下翼面压力仍然较上翼面大，但此时上翼面的压力值与下翼面压力值已经比较接近，从而产生的升力较小。当翅翼的扑动方向反转时(t=0.045 s、0.060 s)，则上翼面压力比下翼面大，但压力差并不显著，因此产生负的升力，但负升力相对产生的正升力而言并不大。为了研究扑翼的升力产生机理，对气动涡的研究意义非凡。翅翼在扑扇过程中展向某平面涡，如图9所示。

图6 升力系数曲线

图7 阻力系数曲线

图8 不同时刻的翅翼压力云图

图9 不同时刻x=10 cm处截面的速度矢量图

由图9可以观察出在平行于YZ轴的某个截面处涡的发展情况，当扑翅在水平位置附近时(t=0.005 s)，由于翅翼上下压差，造成了气流从下翼面绕流到上翼面，从而产生了翼尖涡，当翅翼继续下拍时(t=0.010 s、0.015 s)，上下翼面的压差逐渐减小，涡强度也逐渐减弱，且涡从翼尖向翼根扩散，当翅翼达到极限位置附近时(t=0.025 s)，翼尖涡逐渐破裂消失。下面从三维空间具体研究翅翼周围整个流场的涡的状态：

由图10可以看出，在翅翼刚开始扑扇时(t=0.005 s)，翅翼左侧边缘部分产生了前缘涡，翼尖产生一定强度的翼尖涡，并且在翅翼的边缘前缘涡与翼尖涡汇合，新汇合的涡朝着翼根出发展并渐渐减弱破裂。随着翅翼的继续运动(t=0.010 s、0.015 s)，前缘涡逐渐减弱并朝着弦向发展。当翅翼扑扇到极限位置时，前缘涡强度已经非常小(t=0.025 s)，前缘涡与翼尖涡形成的新的汇合涡也逐渐脱离翼展并减弱。

图10 不同时刻的翅翼三维涡核图

3.2 不同攻角对升力、阻力的影响

鸟类在飞行过程中有往往有不同的飞行姿态，比较典型的姿态变化便是翅翼与来流方向的夹角的变化。本文借助Fluent，对不同攻角下的翅翼的升力与阻力做了数值计算，但在表征翅翼的升力阻力能力时，一般采用周期内的平均升阻力系数作为衡量标准，具体计算方法如下：

(12)

(13)

使用以上公式，对结果进行近似计算，最终得到-8°、0°、8°、16°、攻角时的翅翼周期内平均升阻力系数，如图11和图12：

由图11、图12可以看出，扑翼产生的升力随着攻角的变大，不断增大，攻角为负值时，升力也为负值。扑翼产生的阻力在0°附近时最小，之后随着攻角的增大，阻力也增加。但不同的是，升力随着攻角的增加基本上是线性的，而阻力随着攻角的变化规律性并不强。

图11 不同攻角下的平均升力系数

图12 不同攻角下的平均阻力系数

以上的升阻力系数的变化与鸟类不同时刻的飞行姿态相印证，例如，当鸟类起飞时，身体总是与飞行方向有一个不小的角度差，这里是因为起飞时需要较大的升力，当鸟类在天空中快速飞行时，身体总是与飞行方向保持平行，则主要是为了获得较小的阻力。许多鸟类在雨天快速接近地面捕捉昆虫时，总是将身体下倾，此时升力为负，快速下降，将自身的势能转化为动能，通过调整翅翼角度快速掠过地面。

4 结论

1) 仿生扑翼飞行器飞行过程中的升阻力随时间进行周期性变化，一个周期内影响升阻力变化的原因主要有扑动速度与翅翼所处位置在升阻力方向的面积投影。

2) 当翅翼在水平位置附近时，前缘涡强度较大，随后朝着翼根处方向发展并扩散，当达到极限位置附近时，前缘涡与翼尖涡汇合，并逐渐脱离翼展。

3) 不同攻角下，在一定范围内，翅翼的平均升力随攻角的变化呈线性上升趋势，但阻力的变化较为复杂，在攻角为0°时，阻力最小。另外得出的规律也与鸟类不同的飞行姿态相印证，解释了鸟类起飞、掠飞时升阻力的具体情况。