微型飞行器的低雷诺数、非定常气动设计与分析

2019-01-11昂海松肖天航郑祥明

无人系统技术 2018年4期

昂海松，肖天航，郑祥明

（南京航空航天大学航空宇航学院，南京 210016）

1 引言

微型飞行器（Micro Air Vehicle，MAV）又称微型无人机，与常规航空器的气动特性区别主要是低雷诺数和非定常特征，包括固定翼无人机、旋翼无人机和扑翼无人机。低雷诺数和非定常气动特性极大地影响微型飞行器的飞行稳定性，因此给飞行器的飞行安全和自动控制带来很大的困难。

2 微型飞行器的低雷诺数空气动力学分析

2.1 飞行雷诺数的含义

为了表征流体流动的某些力学特性，常用一些相似参数，如马赫数、雷诺数、费劳德数、斯特劳哈尔数、普朗特数等来描述。其中雷诺数（Reynolds Number）是由英国著名的力学家、物理学家、工程师雷诺［1-2］提出的。雷诺在流体力学上最重要的贡献是发现了粘性流体流动的相似律，他引入了一个常数，表述为：

式中ρ、V、μ分别为流体的密度、流速、粘性系数，L为特征长度，ν=µ/ρ称为运动粘度。

雷诺数是表征惯性力和粘性力之间相互关系的无量纲数。实验表明，雷诺数小时，流体各质点间的粘性力占主要地位，流体各质点平行于管路内壁有规则地流动，呈层流流动状态。流态由层流转变为湍流时的Re值称为临界雷诺数。对于平行流体流过光滑平板的情况，边界层由层流转变为湍流的临界雷诺数约在1×105～3×106之间。目前常规飞行器的雷诺数属于高雷诺数（Re＞1×106）的范围，鸟与昆虫的飞行雷诺数属于低雷诺数（Re＜1×105）的范围。

现代大型飞机，以波音747飞机为例，其机翼最大弦长约14.5m，最大飞行速度约为274m/s，可以求得雷诺数为1.36×108。再以一小型无人机为例（机翼最大弦长约0.6m，最大飞行速度约为60m/s），可以计算得到其飞行雷诺数为2.2×106。目前所研制的微型飞行器，以尺寸15～35cm为例，其飞行雷诺数是介于大鸟（如鹰、信天翁）与小鸟（蜂鸟及昆虫）之间，部分研究中的微型飞行器，如哈佛大学研制的3cm的Flying robotic insects，其飞行雷诺数接近昆虫。

图1 不同飞行物的雷诺数特性Fig.1 Reynolds number characteristics of different flying objects

2.2 雷诺数对气动特性的影响

2.2.1 低雷诺数的附面层特性分析

低雷诺数对流体流动的影响主要是粘性力的作用，着重反映在附面层的特性上。而飞行器表面多为曲面，其附面层特性不但与雷诺数大小有关，而且与飞行器表面的形状有关。小鸟与微型飞行器的雷诺数范围大多为1000≤Re≤200000，此雷诺数范围内，翼表面附面层，在不同情况下，可能是层流、湍流或分离流，甚至会出现层流分离泡的特殊形态［3］。初步分析雷诺数影响如下：

（1）1000≤Re≤10000时，翼面附面层是层流的，而且很难转捩为湍流。大部分昆虫（如蜻蜓和苍蝇）飞行就是处于这样的雷诺数状态。蜻蜓翅膀表面有锯齿状突起，可以产生涡流使得表面气流不发生分离。目前尺寸很小的微型飞行器试验机雷诺数很低，试验表明它在这种层流附面层状态下能保持机翼有较好的气动性能。

（2）10000≤Re≤30000时，翼面附面层仍然是层流，甚至使用人工扰动也不能起作用。对小型无人机的试验表明，附面层一旦分离将不会再附着到翼面上。

图2 烟风洞中低雷诺数层流附面层发生分离的实验Fig.2 Experiment of separation of laminar boundary layer with low Reynolds number in smoke wind tunnel

（3）30000≤Re≤70000时，是目前微型飞行器最关心的雷诺数区域，也是微型飞行器设计师最感兴趣的范围。在这种状态下的翼型选择非常重要，这主要是因为厚翼型（最大厚度≥6%）在层流转捩为湍流后具有很明显的层流分离滞后。同样在雷诺数为50000时，层流分离后的自由剪切流通常不再转捩为湍流回到翼面。薄翼型（最大厚度≤6%）在这个雷诺数状态范围的上界可以表现出良好的性能。

（4）Re≥50000时，大部分翼型都会发生“层流分离泡”。“层流分离泡”的发生不但与雷诺数有关，也与翼型形状和来流迎角有关。分离泡会随着雷诺数的减小而变大，因而会引起翼型性能的恶化。层流分离泡实际是一个不稳定的分离现象，在分离泡后端再附着处存在低频脉动。随着迎角增大分离泡也会前移和发生形状变化。

图3 SD7003翼型（迎角4o，雷诺数Re =60000）的附面层流线图［4］Fig.3 Boundary layer streamline diagram for airfoil of SD7003 （Angle of attack is 4o，Reynolds number is Re =60000 ）

（5）70000≤Re≤200000时，尽管对某些特殊翼型会出现“层流分离泡”的问题，但是一般的翼型都可以获得比较宽的层流范围而改善性能。在这个雷诺数范围内，层流附面层的阻力对雷诺数的变化比较敏感。目前小型无人机在这个范围，飞行雷诺数在不同速度时附面层特性会有大的变化，从而造成设计的复杂性。

（6）对于Re≥200000的情况，翼型性能得到显著改善。翼型的设计可以从大鸟的滑翔和小型无人机等得到大量设计经验。

2.2.2 低雷诺数对气动力的影响

大多数常规飞行器的翼型在低雷诺数时性能很不好，如将NACA0015翼型用在雷诺数为1.5×105的情况时，将会得到很不好的升力曲线和极不规则的阻力曲线。迎角为0o时该翼型就出现比高雷诺数时大得多的阻力，这是因为低雷诺数、零迎角时，机翼上下翼面都发生了气流分离，而当有小迎角时，阻力会有所减少。

翼型的附面层分离和转捩对雷诺数、压力梯度和环境扰动的变化很敏感。分离和转捩在确定附面层的发展中扮演了很关键的角色，影响着翼型的全部气动性能。而机翼和其他部件的气动特性将影响整个飞行器的静态、动态气动特性。因此，掌握附面层对低雷诺数的敏感度，对于设计出适于低雷诺数的翼型对微型飞行器至关重要。

（1）升力特性

低雷诺数时的翼型升力曲线特征与高雷诺数时有较大差别。当雷诺数在1×104～1.5×105之间时，最大升力系数表现得对雷诺数的变化很敏感。在雷诺数小于1×104时，由于粘性占主导地位，这种状态下随着雷诺数的减小，翼型升力曲线斜率随着降低，翼型升力曲线的线性段也随着减少，因此，翼型的最大升力系数总是随着雷诺数的减小而降低。

（2）阻力特性

图4 不同雷诺数与迎角下的翼型升力系数［4］Fig.4 Lift coefficients of airfoils at different Reynolds numbers and angles of attack ［4］

无论是层流附面层还是湍流附面层，随着雷诺数的降低，由于粘性的增强，阻力系数都将不断增加，尤其当雷诺数小于7×104时，阻力系数迅速增大。当雷诺数大于1×105时，湍流附面层的阻力系数高于层流附面层的阻力系数，但是在雷诺数低于这个值时，情况发生了颠倒（如图5所示），而且在雷诺数为1×105附近时，层流附面层的最小阻力系数变化很剧烈，而湍流附面层的最小阻力系数变化相对比较平缓。其原因是雷诺数降到1×105附近及以下时，层流附面层发生分离或出现分离泡，导致阻力增大和剧烈变化。

图5 不同雷诺数和迎角下的翼型阻力系数［4］Fig.5 Airfoil drag coefficients at different Reynolds numbers and angles of attack［4］

（3）升阻比

最大升阻比(CL/CD)max是影响飞行器巡航气动效率的主要参数。综合升力曲线和阻力曲线，可以得出最大升阻比随雷诺数的变化趋势。从图6可以发现在雷诺数低于1×105时，层流翼型的气动性能急剧恶化，湍流翼型虽然气动性能也在恶化，但是情况没有层流那么剧烈。因此雷诺数为1×105是一个临界点，微型飞行器的设计飞行雷诺数范围应该尽量避开这个点。遗憾的是，微型飞行器在起飞和着陆阶段，或者较大顺风情况下，仍然会出现雷诺数等于或小于1×105的飞行状态，这为飞行稳定控制带来困难。不过，我们可以尽量设计微型飞行器的巡航速度避开雷诺数1×105附近值，以避免由于雷诺数的变化引起飞行器的升阻比在这个点附近剧烈变化，导致飞行性能的不稳定。

图6 不同雷诺数下的翼型最大升阻比［4］Fig.6 Maximum lift-drag ratio of airfoil under different Reynolds numbers ［4］

在大雷诺数时，翼型表面粗糙可以引起附面层从层流到湍流的转捩，使摩擦阻力有很大增加，因此这种状态的飞行器总是希望表面做得很光滑。但是雷诺数小于1×105时，表面比较粗糙的翼型比那些光滑的翼型的升阻比变化更加平稳。

从上述低雷诺数对气动性能的影响，可以认为导致微型飞行器设计困难的因素包括：最大升力系数降低，意味着携带载荷的能力降低；最小阻力系数增加，意味着需要更高的推进功率；最大升阻比降低，意味着相同条件下飞行器的续航能力降低。

这几个不利影响同样体现在处于低雷诺数的微型飞行器动力装置上，包括微型螺旋桨或微型旋翼本身的气动效率降低。微型飞行器的尺寸越小，雷诺数越低，上述不利因素的影响就越显著。

微型飞行器的雷诺数约为6.5×104～1.4×105。与常规飞行器的飞行雷诺数相比，这已经是相当低了。在此低雷诺数范围内，微型飞行器的空气动力特性存在其有别于常规飞行器的特点，主要是粘性力的影响强烈，表现为层流分离效应和非定常效应。流场和气动性能易受湍流度和表面粗糙度等因素影响。同时，微型飞行器飞行速度和风速在同一量级，风速的变化会造成雷诺数的剧烈变化，从而使按常规思想设计的飞行器气动性能、稳定性和操纵性急剧恶化。

2.3 低雷诺数MAV的空气动力学计算方法

2.3.1 基于边界层的计算方法

由于低雷诺数翼型呈现粘性力作用较强的特征，传统的边界层理论不失为一种可用解析公式表达、物理概念清晰的分析方法，当然还必须有能求解附面层之外的 “外流”的方法。通常可运用有粘/无粘相互作用的方法，将无粘外部势流与粘性边界层相互迭代计算，以求得飞行器的压力分布和气动力。

运用边界层（附面层）理论来分析低雷诺数翼型特性，空气动力学专家作了不少探索。我们分别运用Thwaites积分方法计算层流边界层，运用Cebeci准则、Michel准则和Wazzan准则进行层流边界层转捩的判断，采用Head积分方法计算湍流边界层，用速度梯度参数λ(x)≤-0.09来判断边界层的分离，经实验证实这是一套行之有效的方法。

然后，进行有粘/无粘相互迭代计算。即边界层仅仅在壁面附近很小的区域内存在，在边界层外，作用在流体微团上的粘性力可以忽略，为无粘性流动。当飞行器流动为小迎角绕流无气流分离时，可以采用一种有粘/无粘相互迭代的方法来进行分析。即先不考虑边界层的存在，而用无粘流理论计算表面压强分布和速度分布，然后用求出的表面速度作为边界层外边界上的主流速度进行粘性边界层计算，找出表面摩擦阻力等各种参数。边界层的计算结果再作为无粘流场计算的边界条件进行无粘计算，反复迭代，就可得到整个流动区域的精确结果。

2.3.2 低雷诺数粘性定常流N-S方程的数值方法

常规基于密度的时间推进法对中等亚声速到超声速的N-S方程的数值方法具有良好的稳定性和收敛性，但该方法用于不可压、低马赫数或低雷诺数流场时会面临收敛速度慢、不稳定和精度低的所谓“刚性”问题。引起这一问题的根本原因在于低速时控制方程系统矩阵特征值对应的特征波速相差太大。

用时间推进技术求解低速/不可压流问题，目前主要的方法是虚拟压缩方法和预处理方法。虚拟压缩方法的基本思想是，在不可压流体方程的连续方程中引入一个虚拟压缩因子，再附加一项压力的虚拟时间导数，使压力与速度显示地联系起来。

时间导数预处理方法是另一种可解决这一问题的有效方法。它在不影响最终定常解的前提下，通过对可压缩流控制方程时间导数项的预处理，使方程系数矩阵的特征值保持在同一量级而不至于相差太大，解决了低速时系数矩阵的刚性问题，使传统的可压缩流方法的计算能力拓展到了低速、低雷诺数不可压流场。我们主要运用预处理方法，特别是在处理具有运动边界的问题上，作了几方面的改进［5］，有效地求解低雷诺数粘性定常流N-S方程，建立了计算微型飞行器低雷诺数粘性定常流的数值方法。

3 微型飞行器的非定常运动特征与设计

除了低雷诺数以外，非定常空气动力也是微型飞行器特殊的基础问题。微型飞行器有两类非定常气动问题：一类是，由于固定翼微型飞行器质量小，在大气环境风（尤其是阵风）的作用下会产生较大幅度的振荡和波动；另一类是，在机翼主动运动（如旋翼、扑翼）的微型飞行器中，其非定常运动特征会对控制系统的设计造成很大的难题。

3.1 固定翼微型飞行器的非定常运动特征

首先分析一下大气扰流的类型，对微型飞行器影响的大气扰流主要有以下几种型式：

（1）阵风

阵风又称作离散突风，表现为确定性的风速变化。为了测试阵风对微型飞行器的影响，我们利用南京航空航天大学的非定常风洞，可以人工产生阵风干扰。风洞非定常机构通过对系统输入信号，可以产生不同形式的阵风变化，如正弦运动、三角波运动、锯齿波运动等变化。其中对典型正弦速度脉动变化曲线速度信号进行了无量纲化处理［6-7］。无量纲速度表示为：

式中U(t)为实时风速，U为平均风速，R为无量纲幅值（R＜1），f为风速脉动频率。通过伺服控制器对非定常机构纵向运动油缸进行控制，由动态压力传感器或热线风速仪测试，也可以得到正弦波形和锯齿波形的阵风变化。

（2）低空风切变

低空是微型飞行器主要飞行的区域。在低空，由于受地形变化和建筑物的影响，风切变经常发生。风切变的表示是用空间两点处的风速除以两点间的距离，也就是说反映空间单位距离上的风速变化。风切变类型比较复杂，如锋面风切变、低空急流、气旋、低空对流等。

（3）大气湍流

自然界中的风由于摩擦、漩涡等原因，往往伴随着湍流。低空飞行环境中的风切变、热交换、地形诱导等因素都会引起湍流的产生。湍流模型一般用平均风上叠加连续随机脉动来表示。

（4）非定常风场

在大气扰动作用下，微型飞行器空速、机体加速度以及飞行姿态都会在瞬间剧烈变化，严重时甚至引起失速［8］。图7记录了微型飞行器在一次真实飞行过程中受到一个4级阵风影响下的部分状态参数变化情况。

考虑微型飞行器飞行过程中穿越非定常风场时，其空速、迎角、侧滑角随时间变化的规律。设微型飞行器飞行中，特定位置（xg，yg，zg）上的大气速度为风速Vwg，风速变化可以表征为：

图7 阵风作用下MAV部分飞行状态变化图［9］Fig.7 Flight state diagram of MAV under gust ［9］

3.2 旋翼微型飞行器的非定常运动特征

旋翼由发动机驱动作旋转运动，因此其叶片及其气动力处于非定常状态。无人机的旋翼一般由两片或更多的桨叶组成。旋翼绕中心桨毂旋转运动，提供无人机飞行的拉力和控制力矩。根据结构形式的不同，无人直升机旋翼可以分为铰接式、无铰式和无轴承式。最常见的是铰接式，全铰接式旋翼结构有挥舞铰、摆振铰和变距铰。挥舞铰使得桨叶能够做垂直于桨盘平面的周期性上下挥舞运动。变距铰通过变距拉杆、摇臂和自动倾斜器实现桨叶的变距（即改变桨叶的迎角）。旋翼的旋转和一定幅度的挥舞运动，就构成了旋翼的非定常运动的主要特征。螺旋桨也是依靠旋转桨叶产生拉力，虽然没有挥舞运动，但有一定的振动，有的螺旋桨具有变距机构，因此，螺旋桨也具有非定常运动的特征。

图8 铰接式旋翼结构［10］Fig.8 Articulated Rotor Structure［10］

图9 桨叶的挥舞运动Fig.9 The swinging motion of the roller blade

3.3 扑翼微型飞行器的非定常运动特征

3.3.1 鸟与昆虫飞行特征

在研究仿生扑翼飞行器之前，我们需要研究鸟与昆虫的飞行特征。早在远古时代，人类就希望能像鸟一样在空中自由翱翔，虽然飞机也是一种仿鸟飞行的形式，但固定翼飞机最基本的升力产生只是利用了鸟的滑翔原理，而鸟与昆虫的灵活扑动的飞行方式，至今尚未在载人飞行器上真正实现。20世纪末前后，仿鸟扑翼飞行才在微型扑翼飞行器上实现。

扑动飞行（Flapping Powered Flight）是一种在翅膀扑动下可以同时产生升力和推力的飞行方式。这也是鸟类、昆虫和蝙蝠等飞行生物长期进化形成的独特的空中运动方式。扑翼飞行器与固定翼飞机的重要不同点，是不需要螺旋桨或喷气推动前进，而是仅仅依赖扑翼的运动既产生升力又产生推力。

鸟翼结构比较复杂，其骨骼由上臂、前臂、腕骨和手骨组成，羽毛也分为不同层次形式和形状。昆虫的翅膀主要由翅梁、翅脉和翅膜组成，昆虫的翅膀结构是自然界中最轻的结构，没有骨髂和肌肉，由角素构成。

通过实验和观察，鸟类飞行不仅是扑动，其飞行姿态也是变化的，并且鸟的飞行轨迹呈波浪型。鸟翼的弦线并不总是与鸟身轴线平行的。小鸟（如蜂鸟）翅膀会绕翼臂作90o以上的转动，大、中型鸟的翅膀在飞行中也能作一定角度的转动。昆虫的翅膀也很复杂，即使是高频扑动的蚊子也由上下扑动和转动所组成。蜜蜂的翼尖扑动呈8字形轨迹。鸟的翅膀除了前飞时保持一定的迎角，下扑时翼后缘部分会发生柔性变形。大、中型鸟上扑时外翼还会作展向折弯变形以减少上扑阻尼力，并且使翼转动一个较大的迎角，使得上扑时仍然有一定升力。经过研究，随着飞行动物体型的减小，其扑翼频率呈不断增大的趋势。飞鸟的扑动频率在2～100Hz之间，如天鹅和苍鹭扑动频率为2Hz，海鸥为5Hz，鸽子为10Hz，最小的蜂鸟可达50Hz。而昆虫的扑动频率通常大于100Hz，曾经有超过1000Hz的记录［11］。

3.3.2 仿生扑翼微型飞行器的非定常运动及其实现途径

图10 鸟翼的构造Fig.10 The structure of bird wings

图11 鸟在不同阶段翅膀的形状和姿势的变化Fig.11 Changes in the shape and posture of bird wings at different stages

我们设计仿生扑翼微型飞行器尚不能完全设计成像鸟或昆虫一样翅膀及其复杂运动状态，目前只能作部分简化的主要扑动形态的仿生设计。下面以我们设计的三个不同类型扑翼飞行器为例，说明不同扑动形态的实现。

（1）单自由度扑动形态实现机构

图12给出了两种曲柄摇杆扑动机构。图12（a）给出的是目前常见且最简单的一种单自由度扑动机构形式，由曲柄2带动连杆3、4，再实现左右摇杆5、6扑动。其特点是结构简单、质量轻（适用于微型飞行器要求）、传动可靠、效率高。为了提高左右摇杆扑动的对称性，我们又提出了一种变异设计方案，如图12（b）所示，即用四连杆代替原二连杆的扑动摇杆，增加了运动灵活性。但需要注意的是，设计连杆3和连杆5的适当长度及夹角∠57，是优化此方案的关键。

图12 曲柄摇杆扑动机构Fig.12 Crank rocker flutter mechanism

我们以图13为例，作运动学分析。机构各杆长度和质量如下：转动摇臂AO长l1，其转角为θ1；左右连杆长l2，其质心分别为S2、S2，转角分别为θ2、θ2；左右扑动杆扛杆短段BC和B'C'长l3，其转角分别为θ3、θ3；固定铰链点C、C'的位置参数分别为（xC，yC），（xC，yC）。

现在分析各杆的位置、速度与原动曲柄之间的关系。右侧连杆AB和摇杆BC的转角位置为：

右侧连杆AB和摇杆BC的角速度分别为：

图13 扑翼机构运动简图Fig.13 Motion of Flapping Wing Mechanism

需要注意的是，上述运动还仅仅是扑翼展向扑动的运动形态实现。作为扑翼飞行器，还必须要有弦向扑动形态，实际上是两自由度扑动。最简单的弦向扑动形态实现是利用柔性扑翼结构，在展向扑动作用下利用惯性与结构弹性来实现弦向扑动。当然，也有研究设计专门的弦向扑动机构，相对比较复杂。

其中，俯仰转动α(t)中，αm为平均俯仰角，α0为俯仰振幅幅度，角速度w=2πf，扑动β(t)中，β0为展向扑动的扑动幅度；ψ1、ψ2分别为弦向俯仰和展向扑动的相位。

（2）三自由度扑翼微型飞行器的设计

昆虫翼运动在一个扑动周期内不但有垂直身体的上下扑动和绕翼轴线的俯仰运动，通常还有平行身体的前后划动。典型的是蜜蜂翼运动，其翼尖运动轨迹呈8字形，同时蜂翼还发生俯仰转动。

朱保利等［12］在其博士论文中提出实现三自由度扑动翼机械机构的三种设计途径，其中“合路驱动”设计思想如图14所示。把三维扑翼运动简化描述为：

其中：f为扑动频率，z1为上下扑动幅度，x1为前后划动幅度，θ为前后滑动相位角，为使扑翼运动轨迹为8字形，θ取90o。研究表明，θ=90o时推进效率最大。

图14 可实现扑翼扭转的8字形轨迹机构Fig.14 The mechanism for flapping wing torsion in the 8-shaped trajectory

利用上述机构，设计出三自由度扑翼微型飞行器运动简图。短轴Q1Q2与CF杆固联，两翼与短轴分别在Q1、Q2组成球销副，可保证两翼随CF杆作俯仰运动（Pitching）；机翼与机架分别在R1、R2处组成滑球副，可将C点的平面8字形轨迹传至翼尖的空间8字形，实现上下扑动（Plunging）和前后划动（Sliding）两个运动。

进一步采用轨迹优化方法，根据扑翼所要求的运动规律来最优化机构的尺寸参数。优化设计变量为各杆杆长、EG的位置角φ、机翼与CF的夹角γ、AE的位置角λ、曲柄AB与DE的初始位置角θ1,0、θ4,0，目标函数设为机构铰链C点再现的s个轨迹点和角位移与给定的轨迹点和角位移之差值平方和为最小的方法来建立。

机构实现的机翼在位置i时的角位移为：

其中：

（3）多段变形扑翼微型飞行器的设计

观察大鸟飞行翅膀下扑时，前臂会同时沿着垂直于翅膀的轴线转动肩关节、肘关节和腕关节以伸展整个翅膀。在下扑结束时，快速地反向转动以上三个关节以收缩翅膀并且向后扭转前臂，在快要到达上扑顶点的时候，鸟会再次快速伸展翅膀为下一次的下扑做准备。这种复杂的伸缩和弯曲变化，大大提高了气动效率［12］。

2012—2015年，南京航空航天大学微型飞行器研究中心参考海鸥和游隼的翅膀结构设计出一种主动结构形状变形扑翼。翔鹰仿海鸥扑翼飞行器翼展2200mm，弦长400mm，空机质量440g。

为了更好的描述大型鸟类翅膀的运动规律，设计扑翼柔性变体函数变速-折叠扑翼扑动函数，设置初始扑动角β0、初始扭转角0α0、时间t、内外段翼扑动偏差角函数δ(t)，以得到多段变形扑翼飞行器的扑动模型：

相较于适用单段扑翼的传统模型，多段变形扑翼飞行器的扑动模型更全面地考虑了扑翼扑动速度和角度随时间的变化，区分内外段机翼的变化，能更好地模拟大型鸟类的实际扑动方式。

4 微型飞行器的非定常空气动力学分析

以往计算流体力学的研究主要集中在常规飞行器的定常流场计算，非定常流场的计算研究刚逐步兴起。微型飞行器非定常气动计算问题要复杂得多。有人想把研究固定翼的线性非定常流动理论运用于扑翼，该理论假设机翼表面和尾流区的涡无限薄，但这一假设对低雷诺数、粘性影响大和大幅度的微型飞行器非定常运动，并不适用。

图15 主动结构多段弯扭变形扑翼结构设计Fig.15 Flapping wing structure with multistage bending and twisting

图16 南京航空航天大学微型飞行器研究中心研制的仿海鸥多段变形扑翼飞行器Fig.16 Seagull-like multi-segment deformable flapping-wing aircraft by NUAA MAV Research Center

对于微型飞行器低雷诺数非定常流的研究仍然处于初始阶段。在非定常流体动力学方程、离散格式、求解精度、动态网格处理、湍流模型等方面，低雷诺数非定常空气动力学的计算方法研究仍面临不小的挑战。低雷诺数非定常空气动力学风洞试验也面临试验理论和方法的突破和创新，同时测量仪器和动态采集设备需具有更高的精度和响应速度。

4.1 旋翼气动力计算的简化方法

4.1.1 叶素理论

旋翼气动力计算常用的简化方法之一是叶素理论，也就是借用机翼翼型气动力系数积分的方法。旋翼桨叶的各个角度和相关速度之间的关系可由图17表示。

图17 旋翼桨叶速度与角度之间的关系Fig.17 Relationship between rotor blade velocity and angle

根据叶素理论，旋翼的拉力系数［14］：

其中，s为旋翼实度，a为二维升力斜率（数值常用5.7值），θ为桨叶桨距角，Φ为来流角，λ为流入比：

当用75%桨叶半径处的桨距角来近似为桨叶平均桨距角，则拉力系数可近似为：

其中，λ流入比对于悬停时旋翼可用下式计算：

旋翼的拉力为：

4.1.2 有限片桨叶固定尾迹分析方法

由于动量理论无法给出诱导速度与螺旋桨桨叶攻角之间的关系，因而无法将其应用于螺旋桨桨叶的气动设计中；叶素理论虽然从桨叶剖面受力情况来分析问题，把桨叶的几何形状与其所受的气动力联系起来，但不能很好解决沿半径的诱导速度分布。为此我们采用了一种有限片桨叶固定尾迹分析方法［15］，其涡系模型如图18所示，为显示清楚起见，图中只画出了一片桨叶的涡系结构，螺旋桨其他桨叶由同样的涡系结构组成。

图18 有限片桨叶固定尾迹模型Fig.18 Fixed wake model with finite blade

该涡系由三部分组成：沿桨叶分布的变强度的附着涡；每片桨叶后缘逸出的由大量自由涡形成的螺旋涡面，该螺旋涡面向后无限延伸，其延伸方向由螺旋桨在空气中的前进速度和桨盘处的等效诱导速度确定；位于螺旋桨旋转中心的中心涡束。

该固定尾迹分析方法的步骤为：

第一步：首先取一等效诱导速度初值，再加上来流速度确定初始尾迹形状，并由比奥-沙伐尔（Biot-Savart）定理确定桨叶处的诱导速度分布。假设在流场中有一根环量为Γ的涡线，涡线上某一微段ds→对其周围M点处的流体质点所激起的诱导速度可按比奥-沙伐尔定理来确定，该诱导速度在直角坐标系中的分量可由下式确定：

其中，dvx、dvy、dvz、lx、ly、lz、dsx、dsy和dsz分别为在三个坐标轴上的分量。

第二步：由上步确定的诱导速度分布计算出桨叶各个剖面的气动迎角，通过剖面二元翼型的气动特性数据计算出剖面的气动力，沿桨叶积分得到整个螺旋桨的气动力；

第三步：由动量理论求出新的等效诱导速度并以此值重新确定新的尾迹形状；

第四步：重复上述三个步骤经反复迭代最终得到稳定的尾迹形状，再计算气动力。

图19 旋翼垂直升降非定常涡系Fig.19 Unsteady vortex system for vertical flight rotor

图20 旋翼前飞非定常涡系Fig.20 Unsteady vortex system of forward flight rotor

4.2 非定常空气动力计算的N-S方程数值方法

4.2.1 非定常流控制方程

三维非定常可压缩流体的雷诺平均N-S方程为：

其中，W为守恒变量，F（W）和Fv分别为对流通量和粘性通量。对任意的可变形或可移动的控制体Ω（t）积分，得到：

4.2.2 非定常流的预处理双时间步推进技术

时间导数预处理法有效地解决了计算定常低雷诺数流场时面临的刚性问题，但它同时也破坏了控制方程的时间精确性。因此，对非定常流场，仅用预处理方法还不能获得时间精确解。在此，我们结合双时间步推进技术［16］，在方程（28）中引入预处理的伪时间导数项，即：

式中，τ和t分别为伪时间和物理时间，为预处理矩阵。双时间步推进包含两个循环，即伪时间的内循环和物理时间的外循环，每个物理时间步当作定常问题用伪时间推进求解。无需时间精确的伪时间内迭代可以采用预处理、当地时间步长等技术提高收敛速度，而物理时间步则可以不受控制方程系统刚性问题的限制。

4.2.3 方程离散与求解

（1）有限体积法空间离散

在任意控制体Vi内，方程离散为：

由于预处理矩阵的引入，在控制体Vi和Vj之间的交接面Sij上，无粘通量的计算有别于无预处理的情况。我们采用一种类似Roe格式的通量差分分裂格式。

（2）时间离散

时间步的离散，伪时间用一阶后向差分，物理时间用隐式k阶后向差分，方程变为：

m和n分别为伪时间步和物理时间步。RESi为残值，定义为：

4.2.4 非结构嵌套网格方法

对相对运动幅度较大的流场，常规的网格变形技术已无法适用。我们设计了适应于运动幅度大、有结构变形微型飞行器非定常形态的嵌套网格方法［5］。嵌套网格技术把流场适当地划分为多个具有重叠部分的区域，各个区域分别生成独立的网格并在其上求解，在重叠区上通过网格间插值进行区域间信息交换。我们发展了一种非结构动态嵌套网格方法，适合流场中多体复杂运动的非定常流计算。

图21 仿蜻蜓扑翼微型飞行器计算嵌套网格Fig.21 Nested grid of dragonfly-like aircraft

图22 前后双扑翼非定常流场计算结果Fig.22 Results of flapping wings with horizontal direction in the unsteady flow field

图23 固定翼微型飞行器计算网格Fig.23 Computational grid of fixedwing micro air vehicle

图24 旋翼非定常流场计算结果Fig.24 Results of rotor wing in the unsteady flow field

4.3 低雷诺数非定常气动设计软件与应用

根据上述非定常流N-S方程的求解方法和我们发展的动态网格生成技术，肖天航［5］用C++语言开发了三维定常和动态网格非定常N-S方程的求解程序3D2MUFS（3-Dimensional Dynamic Mesh Unsteady Flow Solver）。该求解器适用于任意非结构混合网格，可用于计算低雷诺数非定常大幅度运动的微型飞行器无粘/粘性流场。3D2MUFS的程序流程如图25所示。

南京航空航天大学微型飞行器研究中心首先利用该软件仿真揭示了鸟扑动飞行产生推力的反卡门涡街。1996年Jones［17］和1999年Lai，Platzer等［18］通过流动显示试验进一步观察到了机翼扑动产生的卡门涡街或反卡门涡街的存在。当扑翼扑动速度增大时，尾流所示呈反卡门涡街，两排涡的方向与卡门涡街相反，涡间尾流在来流方向上有与来流同向的流动，则扑翼受到与尾流反向的作用力，即为推力［16］。

为进一步研究柔性扑翼推力产生的机理，我们计算仿鱼游动时，发现鱼虽然没有做全身扑动，但身体后部（特别是尾鳍）作柔性摆动，鱼便会前进。

为考察扑翼柔性对气动特性的影响，我们数值研究了三种不同柔性程度的仿鸟扑翼流场，并考虑扑动频率对每种柔性翼的影响。结果表明，St数很小时，流场尾迹的涡街尚未成型，扑翼受到阻力作用；每种柔性翼，周期平均推力和气动功率都随着St数的增大而增大，St数在0.2～0.3附近时，推进效率出现最大值。我们的计算结果与美国实验的反卡门涡街结果很相似，从理论上揭示了鸟扑动产生推力的机理。

南京航空航天大学微型飞行器研究中心应用该软件先后进行了CN系列固定翼微型飞行器、翠鸟扑翼飞行器、差动扭转金鹰扑翼飞行器、仿大鸟雄雕扑翼飞行器、仿海鸥多段变形扑翼飞行器、黑鹫多旋翼飞行器等二十几种微型飞行器的气动设计和飞行验证机试验，验证了该中心建立的一整套微型飞行器气动设计方法的有效性。