0 引言

自然界的鸟类经过亿万年的进化，形成了出色的飞行能力，其飞行效率、机动性都远远超出了目前所有的人造飞行器。具有较低能耗、噪声，同时具有优越机动性能的小型扑翼飞行器是未来航空器发展的一个重要方向，为开展此类研制工作，首先需要对鸟类扑翼运动过程中的升力、推力产生机制以及流动控制机理进行深入细致的研究。

鸟类翅膀的构造非常精细，由肌肉、骨骼、羽毛和多个关节组成，扑翼过程中伴随着翅膀的弦向转动、展向的收缩-伸展以及扭转^[1]，其飞羽在上拍、下拍过程也会规律性的收缩-合拢。鸟类的升力主要源于翼型本身以及翼型的非对称拍动，例如海鸥的翼型剖面和S1223翼型非常相似，在低速情况下具有很高的升阻比，而鸟翼在下拍过程中会完全打开，上拍过程中则会尽量收缩，且其上拍、下拍的时间也存在很大差异，这种非对称拍动的目的也是为了增强升力。

目前国内外常见的关于鸟类扑翼飞行的研究工作大概可以分为两类：一类通过观测，采用数据统计和量纲分析的手段，分析鸟类形态和各飞行参数之间的关系，并得到一些经验性的公式，如Greenewalt^[2]、Tennekes^[3]、Pennycuick^[4]等人的工作。另一类则是采用试验、理论或近似的数值分析方法，针对具体的扑翼过程开展研究，如Jones^[5]、Angela^[6]、Pennycuick^[7]、Rayner^[8]、Phlips^[9]等人的工作。其中“面元法”作为一种简单有效的数值模拟手段，在鸟类扑翼的研究中使用较多，如Smith^[10]采用该方法并结合有限元模型对柔性扑翼运动进行了研究，国内的昂海松研究团队^[11]、余永亮^[12]等也分别采用该方法对昆虫及蝙蝠的扑翼进行了研究。

虽然这些针对鸟类扑翼运动的研究已经初步揭示了鸟类的升力、推力产生机制，然而这些定性的认识仍不能满足扑翼飞行器设计上的需求。鸟类扑翼运动的许多细节问题如翼型、迎角、展向变形/扭转、上/下拍动时间比、拍动频率等对其力学特性有非常重要的影响，因此非常有必要针对这些细节问题进行深入研究。

计算流体力学的迅速发展为扑翼运动的精细化研究提供了条件，且已经在昆虫的扑翼研究中得到了广泛应用^[13-16]。然而由于鸟类扑翼过程的复杂性，三维情况下复杂扑翼过程的精细化数值研究还不多见，很多数值研究工作仍基于二维或者简单的三维扑翼开展^[17-18]。

本文以海鸥的扑翼飞行为研究背景，选择二维的S1223翼型作为研究对象，采用动态混合网格技术及非定常数值模拟方法，对其刚性俯仰-沉浮运动进行了研究，得到了下拍时间比、拍动迎角等参数的影响规律。这些工作为后续的复杂三维问题的数值研究奠定了基础。

1 数值计算方法

本文的数值模拟基于HyperFLOW软件平台^[19]开展。该软件平台是中国空气动力研究与发展中心研发的具有完全自主知识产权的大型CFD多学科通用求解平台，具有优越的体系架构和生态系统，并已集成了结构/非结构NS方程流场解算器、动态混合网格生成技术、飞行力学/流体动力学一体化算法等，可进行完全气体和化学非平衡气体的定常/非定常计算。以下对其中的动态混合网格技术以及非定常算法进行简要介绍。

1.1 动态混合网格生成技术

常用的网格变形技术有弹簧松弛法、求解偏微分方程的方法以及插值方法等。在之前的研究工作中，作者所在的课题组建立了弹簧松弛法和插值法相结合的混合网格变形技术，具有较好的变形能力和动网格生成效率，并通过结合局部网格重构技术，提高了针对大变形、大位移、相对运动等复杂动边界问题的适应能力^[20-22]。

在最近的研究工作中，我们进一步发展了基于径向基函数(RBF) 的网格变形技术，该技术具有优越的网格变形能力，是目前国内外CFD工作者争相研究的一个热点。而标准的RBF方法在处理大规模网格时效率极差，为了提高RBF网格变形技术的适用性，我们参考文献[23]的做法，通过选择有限的参考点来减少RBF算法中矩阵的规模，以提高计算效率。

1.2 非定常数值模拟方法

基于动态混合网格的NS方程解算器采用了格心型的有限体积格式，时间离散采用二阶的欧拉后插方法，为提高非定常计算效率，采用了双时间步算法和BLU-SGS隐式计算方法。算法具体细节请参见文献[22]。

动网格下的NS方程离散算法必须满足几何守恒律。运动网格下的几何守恒律(GCL) 特指体积守恒，对于控制方程，假设流场均匀，则可得到GCL方程如下：

(1)

其中V为控制体Ω的体积，∂Ω为控制体边界，w表示控制体外边界的运动速度，n为控制体边界的外法向。导数符号“d”表示跟随控制体的随体导数，而不是跟随流体的物质导数。

式(1) 表示控制体单元的体积变化率等于边界面法向运动速度的面积分，其在微分意义下是恒成立的，而计算所采用的离散格式也需要使式(1) 得到满足，才不至于在流场内引入额外的误差。

满足式(1) 的数值格式很多，详见Farhat^[24]、Kallinderis^[25]、刘君^[26]、Wang^[27]等的工作。这里我们采用了一种较为便捷的方法：直接对边界面的法向速度进行约束。以一阶欧拉隐式方法为例，式(1) 的离散形式为：

(2)

式中v=w · n表示控制体边界面的法向速度。对于单元的每一个离散的边界面，限定其法向速度的求解方法：

(3)

其中ΔV_jⁿ表示单元第j个面扫过的面积。则式(3) 是满足式(2) 的充分条件。同理，对于二阶欧拉隐式格式，应当选用如下的法向速度求解方法：

(4)

在随后的研究中，我们进一步从误差分析的角度对各种运动网格几何守恒算法进行了理论分析，发现可以将现有的几何守恒算法归为两类，即限制整体积分误差的“体限制方法”和限制每个边数值误差的“面限制方法”，并分别可以针对不同的时间格式写成统一的形式。进一步的对比发现，各种几何守恒算法的计算结果大体一致，但是我们发展的几何守恒算法在三维情况下更为简便。具体的理论分析详见文献[28]。

2 翼型及俯仰-沉浮运动

根据文献描述^[29]，本文选用二维S1223翼型作为研究对象，翼型弦长c=0.2m。将整个拍动周期分为下拍和上拍两个阶段，其沉浮运动规律如下：

(5)

其中dy为翼型的纵向位移，A为振幅，T为拍动周期，k表示下拍时间和整个拍动周期T的比值。本文令振幅A=0.2m。

翼型绕其1/4弦长处进行俯仰，俯仰角(低头为正) 变化规律如下：

(6)

其中θ₀表示上拍和下拍过程的起始拍动角，φ表示俯仰角变化与上拍-下拍运动之间的相位差。本文令φ=0、θ₀=0，因此俯仰角振幅θ_down、θ_up分别代表最大下拍角和最大上拍角，且最大拍动角出现在上拍或下拍过程的中间时刻。

本文只针对二维问题进行数值模拟，其相当于三维问题的一个展向截面(如图 1所示)。

取来流速度U=10m/s，翼型拍动频率。采用海平面的空气密度ρ=1.225kg/m³为参考密度，以弦长c为参考长度，以来流速度U为参考速度，得到雷诺数Re=1.35×10⁵，减缩频率为：

(7)

本文需要研究俯仰-沉浮运动的能耗特性，能耗系数C_P的定义如下：

(8)

其中f表示流体作用在物体边界上的力，w表示边界元的运动速度，∂Ω表示物体表面。

3 数值结果分析

采用混合网格离散计算域，物面附近采用各向异性的四边形单元。翼型外围采用三角形网格单元，并在翼面附近以及尾迹内进行了加密处理，网格单元总数8.6万(图 2)。计算采用SA湍流模型，不考虑转捩因素的影响，采用二阶空间和时间格式，一个拍动周期采用400个非定常时间步离散。

3.1 下拍时间比对升力及能耗的影响

观察发现^[30]，海鸥等鸟类巡航飞行时其下拍时间占到整个拍动周期的60%~80%，说明下拍时间比k与其力学特性关系密切。本节不考虑翼型的俯仰，仅考虑其上下沉浮运动，针对k=0.40~0.75等几种情况进行数值模拟。

图 3(a)所示为不同k值情况下升力系数在一个周期内的变化曲线，横坐标为拍动的相位角。下拍过程中(Phase < π)，在翼型的随体坐标系下来流的等效迎角为正，因此产生正升力，但是随着k值的增加，翼型下拍的速度减小，等效迎角减小，因此升力峰也随之减弱；上拍过程中，翼型等效迎角为负，因此翼型产生负升力，且随着k的增加，上拍速度增加，负升力峰增大。但是由于翼型本身具有较好的升阻比特性，在上拍和下拍的起始阶段(Phase=0、π附近) 仍然可以维持较大的升力系数，这对于增大整个拍动周期的平均升力是有益的。俯仰运动的瞬时能耗系数与翼型拍动速度的大小关系密切，在拍动速度最快的时刻(Phase=0.5π、1.5π) 会出现能耗的峰值(如图 3b所示)，随着k值的增加，下拍过程的能耗峰值减弱，而上拍过程能耗峰值则会增加。

图 4所示为拍动一个周期内几个典型相位的涡量云图，在Phase=0.5π、1.5π两个相位下，不同的k值对应的流场结构基本类似，只是由于k值的变化影响了其下拍和上拍过程的等效迎角，因此背风区内的旋涡大小、形态等存在一些差异；而Phase=0、π两个相位下，翼型的法向运动速度为0，但是加速度处在整个拍动周期的峰值，此时不同k值条件下流场的旋涡结构存在较大差异。以Phase=0为例，此时处于上拍结束的急剧减速以及下拍起始的急剧加速阶段，翼型带动周围流体急剧加速/减速，受“附加质量效应”影响较大，且由于本文设计的运动规律在此时的加速度不连续，因此导致升力系数存在“不规律”的波动(图 3c)。然而，由于翼型拍动的减缩频率较小，来流速度相对于拍动速度是一个较大的量，因此这些“不规则”的旋涡结构会很快消失在尾迹，升力系数也会很快进入单调增加的趋势。

上述的数值模拟结果表明，翼型的下拍过程对升力的贡献较大，而上拍过程则产生负升力。实际上，鸟类在扑翼过程中正是通过翅膀展向的折转来突出下拍的“有益”作用(下拍时翼面全部展开、增大有效升力面)，减弱上拍的“不利”影响(上拍时适当折叠扑翼，减小上拍迎风阻力面)。为了近似反应这种“三维效应”，本文对升力、阻力、能耗等引入一个近似的加权系数：

(9)

λ在0.5和1.0之间周期性变化，φ表示加权系数和扑翼运动之间的相位差，本文令φ=0。将计算得到的瞬时气动力系数乘上加权系数λ，则本文的二维问题可近似看作一个三维的翼型在做俯仰-沉浮运动的同时改变其展向长度。

考虑加权前后，扑翼运动的时均升力系数、时均能耗系数以及它们之间的比值随k值的变化如图 5所示：与原始模式(Original Mode) 相比，加权模式(Weighted Mode) 下时均升力系数整体上有所增加，时均能耗系数则在整体上有明显的降低。加权模式下，随着k的增加，升力系数呈增加的趋势，而时均升力/能耗比先增加后减少，并在k=0.65~0.7附近出现最大值。这个结论和自然界一些鸟类的观测数据保持一致^[29-30]。

3.2 拍动角的影响

本节取k=0.7，令θ_down=θ_up，分别针对-5°~45°等若干拍动角的情况进行了数值模拟，分析拍动角对扑翼过程气动力特性的影响。

一个拍动周期内升力系数、阻力系数的变化情况如图 6所示，图中横坐标为扑翼运动的相位。分别针对下拍过程和上拍过程进行单独分析。下拍过程中(Phase=0~π)，随着拍动角θ_down的增加，相对于翼型而言来流迎角减小，因此升力系数的峰值减弱，且当θ_down>25°时甚至会出现负的升力峰；就阻力系数而言，在一定的θ_down范围内为负值，说明产生了“推力”作用，但是随着θ_down的增加，“推力”的峰值逐渐减小，并最终导致阻力的产生。在上拍过程中(Phase=π~2π)，相对于翼型而言来流的等效迎角为负，因此当上拍角θ_up较小时，会出现负的升力峰值，但是随着θ_up的增加，负升力峰逐渐减弱，并在一定条件下变为正值。

采用和3.1节相同的方法，为升力、阻力系数引入公式(9) 所示的加权系数λ，得到拍动过程的时均力学系数如图 7所示，图中横坐标为最大拍动角。下拍过程的时均升力系数在θ_down=-2°时出现最大值，之后随着θ_down的增加而单调减小，在θ_down=30°左右变为负值；时均阻力系数的最小值出现在θ_down=0°附近，之后随着θ_down的增加而单调增加，在θ_down=21°左右变为正值。上拍过程的时均升力系数随上拍角的增加而单调增加，并在θ_up=25°附近变为正值，说明上拍过程仍可能产生升力；就阻力而言，在各个θ_up条件下其时均值均为正，只是在一定的条件下具有较小的量，从图 7可以确定出这个阻力最小的拍动角约为θ_up=20°。

这些数值结果表明，翼型的拍动角度对其升阻力特性影响极大，因此采用合适的拍动角是其提高升力、增加推力的关键。下拍过程对升力、推力的产生都是有利的，且θ_down=0°时具有较好的升阻力特性。由于本文所采用模型的局限性，在各个上拍角条件下，上拍过程中均无法产生推力，但是在合适的角度下能够保证较小的阻力，且当上拍角大于一定数值时仍能够产生升力。

3.3 优化的拍动方式

根据3.2节的分析，选择θ_down=0°，针对θ_up=20°的情况进行了数值模拟。

拍动一周期内升力系数、阻力系数的变化曲线如图 8所示，横坐标表示拍动的相位。拍动过程中几个典型状态下的瞬时流线及压力云图如图 9所示(为体现来流的“等效迎角”，速度场减去了翼型在y方向运动速度)。

下拍过程中(Phase=0~π)，翼型背风区内的流动分离较弱，仅在个别时刻在尾缘附近出现了较小的流动分离，因此整个下拍过程中翼型的升力系数均为正，且存在较大的升力峰值。从压力云图也可以看出，下拍过程中在翼型头部存在较强的负压区，可诱导出较大的升力以及正推力。

上拍过程的起始及结束阶段，翼型的升力系数仍为正值，仅在下拍过程的中间阶段出现负升力。这是由于上拍过程中翼型存在较大的上拍角度，相对于来流，翼型的等效“负”迎角减小，因此上拍的不利影响会得到减弱。

下拍-上拍的转换阶段(Phase=0、π)，由于本文所采用的运动规律在此时的俯仰角速度不连续，导致力学系数出现跳跃，但是这些瞬时的强非定常现象对整个拍动周期的影响相对较弱。

采用与3.1节相同的做法，引入公式(9) 所示的加权系数，将力学系数进行加权，得到的整个拍动过程的时均升力系数、阻力系数和能耗系数分别为1.10、-0.14、0.24。定义效率：

(10)

则扑翼的推进效率约为58.3%。

4 结论

本文以海鸥等大型鸟类的巡航飞行状态为研究背景，选择S1223翼型以及相应的扑翼运动参数作为研究对象，希望通过这些初步的研究，探讨翼型俯仰-沉浮运动过程的升力、推力的产生机理。

本文选用的扑翼方式具有较小的减缩频率，和昆虫等的扑翼方式不同，其非定常效应较弱，翼型的静态因素对气动力特性影响则较大。翼型通过沉浮运动改变来流的相对迎角，在下拍过程中，来流等效迎角为正，且由于非定常流动的迟滞作用，翼型背风区内的流动分离现象很弱，鉴于S1223翼型优秀的升阻比特性，气动力合力指向翼型前上方，同时产生“升力”和“推力”(图 10所示)。

文献[18]的数值模拟结果表明，三维柔性翼型在上拍过程仍有较大的推力产生。而在本文所选用的模型及运动方式条件下，上拍过程没有能够产生推力，其原因可能在于：(1) 本文选用的翼型在负迎角条件下具有较差的力学特性；(2) 本文采用刚性的俯仰-沉浮运动方式，没有考虑到实际情况下翼型的柔性变形；(3) 所选用的扑翼运动参数如频率、拍动角等存在差异。

下拍时间、拍动角等扑翼参数对翼型的气动力特性影响较大。在本文的数值模拟条件下，当下拍时间占整个拍动周期的65~70%左右时具有较好的时均升力/能耗比，而下拍角取为0°，上拍角取为20°时具有较好的升力、推力产生性能。

本文仅针对刚性的二维扑翼模型进行了数值模拟，初步分析了下拍时间、拍动角等的影响，对其推力、升力产生机理有了定性的认识。而三维情况下的扑翼运动具有更多的自由度，并且存在柔性变形如翼沿展向的弯曲、扭转等。这些三维效应对升力、推力的影响很大，下一步我们将重点开展这方面的研究工作。