0 引言

扑翼飞行器一直是国内外很多研究机构的一个研究热点，其设计灵感来源于自然界的昆虫、鸟类等飞行生物。这些飞行生物经过了亿万年的进化，形成了非常出色的飞行能力。不同于固定翼飞行器，鸟翼即是升力产生机构，也是推力产生机构，其飞行效率更高、噪音更低，并且由于其气动力产生方式更为灵活，因此机动性能更加出色。和固定翼相比，扑翼飞行的气动力产生机制以及流动控制机制更加复杂，通过对自然界的鸟类、昆虫等生物扑翼飞行进行观察和研究，将有助于推动扑翼飞行器研制工作的进展。

中等体型鸟类的扑翼和昆虫的扑翼存在较大差异。首先，两者扑翼的减缩频率不一致。昆虫扑翼的减缩频率很高，雷诺数较低，在翼的拍动过程中相对于翼面的来流迎角很大，因此会产生非常明显的前缘涡以及后缘涡，其中前缘涡在翼的运动过程中保持不脱落，从而可以维持一个较高的升力^[1]。而中等体型鸟类飞行速度较快，扑翼的减缩频率较低，雷诺数较高，翼型本身的作用更加明显^[2-3]。其次，两者扑翼的结构和运动方式差异较大。昆虫翼一般由少量翅脉和大面积的柔性薄膜组成(如蜻蜓翼)，翼面作往复式旋转，拍动过程中往往伴随大角度的“8”字型扭转，翼面相对于翅脉的柔性变形较小，且主要是被动弹性变形，因此早期的研究中将其视为刚性薄板进行简化处理^[4-5]；鸟翼的结构以及变形机制较为复杂，其由肌肉、骨骼、羽毛和多个关节组成，构造非常精细。其扑动过程具有较多的自由度，如翼的拍动、折转、收缩、扭转、飞羽的收缩-合拢等，主动以及被动变形非常剧烈。

由于鸟类扑翼的减缩频率相对较低，雷诺数相对较高，翼型本身对升力的贡献作用更为明显，因此很多经典的研究工作正是基于翼型理论开展的，如Jones^[6]、Angela^[7]、Pennycuick^[8]、Rayner^[9]、Phlips^[10]等人的工作。这些近似的分析工作定性的给出了鸟翼扑动过程中非定常升力、推力的产生机制。其中“面元法”作为一种简单有效的数值模拟手段，在鸟类扑翼的研究中使用较多，如Smith^[11]采用该方法并结合有限元模型对柔性扑翼运动进行了研究，国内的昂海松研究团队^[12]、余永亮^[13]等也分别采用该方法开展了扑翼问题的研究。

虽然这些针对鸟类扑翼运动的研究已经初步揭示了鸟类的升力、推力产生机制，然而这些定性的认识还不能完全满足扑翼飞行器设计上的需求。鸟类扑翼运动的许多细节问题如折转收缩变形、展向扭转变形、相位等对其力学特性有非常重要的影响，其中必定蕴含着重要的增升、减阻、减能耗机制需要流体力学工作者去不断研究、挖掘。

计算流体力学的迅速发展为扑翼运动的精细化研究提供了条件，且已经在昆虫的扑翼研究中得到了广泛应用^[14-18]。然而由于鸟翼扑动过程的复杂性，三维情况下复杂扑动过程的精细化数值研究还不多见，很多数值研究工作仍基于二维或者简单的三维扑翼开展^[19-20]。本文基于文献资料构建了海鸥翼几何模型并设计了简化的扑动运动模型，采用动态混合网格技术及非定常数值模拟方法，对扑动过程进行了数值模拟，对拍动角、折转角的影响进行了对比分析。

1 数值计算方法

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

1.1 动态混合网格生成技术

在之前的研究工作中，作者所在的研究团队建立了弹簧松弛法和背景网格映射法相结合的混合网格变形技术，具有较好的变形能力和动网格生成效率，并通过采用局部网格重构技术，提高了针对大变形、大位移、相对运动等复杂动边界问题的适应能力^[22-24]。在最近的研究工作中，进一步耦合了基于径向基函数^[25](RBF)的网格变形技术，其具有优越的网格变形能力。标准的RBF方法在处理大规模网格时效率极差，为了提高其适用性，我们通过选择有限的参考点来减少RBF算法中矩阵的规模，以提高计算效率。图 1所示为本文的RBF及相应的参考点选择算法流程：计算准备阶段，根据物面的初始位移，通过贪婪算法确定出参与矩阵求解的参考点；动网格生成的每一步，先对物面点进行插值，并进行误差检验，如果物面误差不满足计算需求，则重新建立物面参考点序列，之后再对空间点进行插值。

根据文献中的观测数据^[31]，我们建立了海鸥翼的三维模型(如图 2所示)。其弦向截面为S1223翼型，根弦长c=0.2m，展长L=0.5m。本文采用混合网格离散计算域，翼前后缘附近采用非结构的四面体和三棱柱网格单元，其它区域采用六面体单元，并在尾流以及翼面上下等局部区域进行了网格加密处理，网格单元总数282万(由于外形比较简单，这里不再显示具体网格分布)。图 2所示的红点为通过贪婪算法选择的RBF参考点，物面网格点数4万，通过筛选只保留了1000个控制点，因此可以极大的提升动态网格生成效率。图 3所示为物面点的最大以及平均误差随参考点数目的变化情况，当选择1000个控制点时，物面的最大误差可保持在0.01mm以下。

1.2 非定常数值模拟方法

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

在之前的研究中，我们基于误差分析以及数值测试，对运动网格下的几何守恒律问题进行了细致的研究，将现有的几何守恒律算法^[26-29]归纳为两类：限制整体积分误差的“体限制方法”和限制每个单元边界面数值误差的“面限制方法”。在此基础上提出了一种简便的满足几何守恒律的算法，其只需要约束单元边界面的法向速度即可。以二阶欧拉隐式格式为例，采用如下的法向速度求解方法：

$ v_j^{n + 1} = \frac{{1.5\Delta V_j^n - 0.5\Delta V_j^{n - 1}}}{{S_j^{n + 1}\Delta t}} $

(1)

其中ΔV_jⁿ表示单元第j个面在n~n+1时刻扫过的体积，S_jⁿ⁺¹表示该面在n+1时刻的面积。利用该算法可以很好地保证动网格非定常计算的几何守恒，能得到较好的非定常计算结果。具体的理论分析及数值测试结果详见文献[30]。

2 翼型及运动方式

本文采用如下的坐标系定义：坐标原点位于翼根部翼型前缘顶点；x轴指向来流流动方向，y轴指向翼型上方，z轴指向展向方向。

真实海鸥的扑翼运动过程较为复杂，包括绕根部的拍动、沿展向的折转、扭转、收缩等。除此之外，还包括一些其它的提高力学特性的运动机制如飞羽的打开/合拢、翼的被动柔性变形等。因此完全考虑真实的扑翼过程是非常困难的。

本文将海鸥翼分为两段，靠近根部的第一段长度0.2m，剩下的为第二段，长度0.3m。只考虑翼型的折叠和拍动，因此将扑翼运动简化为如下两种简单运动的叠加：

1) 翼型的变形：第二段翼的折转，折转角定义为θ。

2) 翼型的刚性拍动：绕x轴的拍动(拍动角γ，下拍为正)，旋转中心设置在翼根部。

图 4给出了翼型扑翼运动的示意图。折转角及拍动角的变化规律如下：

$ \begin{array}{l} \theta = 0.5{\theta _{{\rm{max}}}}\left[{1-{\rm{sin}}\left( {\varphi + \Delta \varphi } \right)} \right]\\ \gamma = - {\gamma _{{\rm{max}}}}{\rm{cos}}\varphi \end{array} $

(2)

其中θ_max、γ_max分别为最大折转角、最大拍动角。φ表示扑翼运动的相位，Δφ表示翼的折转运动和拍动运动之间的相位差。相位φ的表达式如下：

$ \varphi = \left\{ {\begin{array}{*{20}{c}} {{\omega _{{\rm{down}}}}t, 0 < t < kT}\\ {\pi + {\omega _{{\rm{up}}}}\left( {t - kT} \right), kT \le t < T} \end{array}} \right. $

(3)

式中k表示下拍时间和整个拍动周期的比值，ω_down、ω_up分别为下拍、上拍过程的角频率：

$ {\omega _{{\rm{down}}}} = \frac{\pi }{{kT}}, \;\;\;{\omega _{{\rm{up}}}} = \frac{\pi }{{\left( {1 - k} \right)T}} $

(4)

k值取为0.7。本文的参数设定以海鸥的巡航飞行为研究背景，来流速度U固定在10m/s，0°迎角，扑翼频率f=10/π(Hz)。以翼根弦长c为参考长度，来流雷诺数Re=1.35×10⁵，减缩频率为：$\frac{{\pi fc}}{U} = 0.2$。在拍动角30°的情况下，前进比为3(前进比定义为来流速度和翼尖平均运动速度之比：$\frac{U}{{4f{\gamma _{{\rm{max}}}}L}}$)。

3 数值结果分析 3.1 拍动角的影响

以下分别针对拍动角γ_max=5°、10°、20°、30°等四个状态进行数值模拟，此时折转角为0。该状态下翼没有变形，仅做单自由度的拍动。

图 5(a)所示为四个状态下一个拍动周期内的升力系数、阻力系数以及能耗系数变化曲线。能耗系数定义为：

$ {C_P} = \frac{P}{{\frac{1}{2}{\rho _\infty }{U_\infty }L_{{\rm{ref}}}^3}}{\rm{, }}\;\;P = \oint\limits_S {\mathit{\boldsymbol{f}}\cdot\mathit{\boldsymbol{v}}{\rm{d}}\mathit{s}} $

(5)

其中S为翼型表面，f为固壁上的气动力，v表示固壁的运动速度，ρ_∞为来流密度，U_∞为来流速度，L_ref为参考长度，本文取L_ref=1.0m。

一个周期内升力系数、阻力系数均出现了一个波峰和波谷，而能耗系数出现两个波峰。升力系数峰值、阻力系数最小值均出现在下拍速度最快的时刻(t/T=0.35)；升力系数最小值、阻力系数最大值出现在上拍速度最快时刻(t/T=0.85)；而能耗系数峰值则分别出现在拍动速度最快的两个时刻。随着最大拍动角γ_max的增加，升力系数、阻力系数以及能耗系数的峰值均单调增大。图 6给出了下拍最快时刻几个典型状态下翼上下表面的压力云图，可见第二段翼的压力分布受γ_max影响最为明显：随着γ_max的增加，翼运动速度加快，第二段翼上表面的前缘出现较大的负压区，下表面出现较强的高压区，对升力、推力产生积极影响。

时均力学系数随最大拍动角γ_max的变化如图 5(b)所示。随着γ_max的增加，下拍过程(t/T=0~0.7)的时均升力系数单调增大，上拍过程的时均升力系数则单调减小，并在一定的范围内变为负值，而整个拍动周期的时均升力系数的变化则相对较小。下拍过程的时均阻力系数随γ_max的增加单调减小，并在一定γ_max范围内为负，说明起到了“推力”作用，而上拍过程则正相反，时均阻力系数随γ_max的增加而增大。由于下拍过程所占整个拍动周期的比重较大，因此总的时均阻力随着γ_max的增加而减小。上拍、下拍以及总的时均能耗系数均随着γ_max的增加而单调增大。

对于上述单自由度刚性拍动的翼型，增加最大拍动角能够较为明显的改善下拍过程的升力特性，并能够在下拍过程产生推力。这实际上是容易理解的，在相同的扑动频率情况下，增加最大拍动角相当于增大扑动速度，由此导致扑翼各截面“等效迎角”变化加大，进而导致升阻力特性的变化。然而，由于上述计算由于没有考虑翼型的弯曲、扭转等变形机制，上拍过程中气动力的不利影响较大，因此导致扑翼总的力学性能较差。

3.2 第二段翼折转角的影响

翼面沿展向的折转变形是鸟类扑翼过程中的一个显著特征，通过展向的折转，改变了翼型的有效迎风面积，因此可以改进拍动过程的升阻力特性。本节在3.1节的基础上，考虑第二段翼的折转，分析最大折叠角度、折转角相位等对拍动过程升阻力特性的影响。

首先令公式(2)中的Δφ为0，表明下拍最快的时刻折转角为0，上拍最快的时刻折转角最大。最大拍动角γ_max固定为30°，针对θ_max=0°~50°等若干状态进行数值模拟。

拍动一周期内的升力系数、阻力系数以及能耗系数变化曲线如图 7(a)所示。由于折转角按照正弦曲线变化，下拍过程中完全展开，折转角为0，因此下拍过程的气动力受最大折转角的影响较弱；而上拍过程中折转角最大，对气动特性的影响相对较大：随着θ_max的增加，上拍过程的负升力峰减弱，阻力、能耗的峰值均减小。图 8所示为上拍最快时刻(t/T=0.85)翼型的压力云图及涡量Q等值面，此时翼型折转角最大，随着θ_max的增加，第二段翼前缘的压力峰值逐渐减弱，翼梢前缘诱导的分离涡也呈减弱趋势。从时均力学特性的变化情况可以看出(图 7b)，上拍过程受折转角影响较大，随着θ_max增加，上拍过程的气动性能得到较为明显改善：升力增加，阻力、能耗均降低。而下拍过程所受影响相对较弱。图 9给出了不同折转角情况下整个拍动周期内的升阻力系数和能耗系数的比值，可见随着折转角的增大，升力能耗比增加，阻力能耗比减少，说明增加折转角对整个拍动周期是有益的。然而本文计算得到的能耗为翼型对流体所做的功，并不等同于翼型扑翼运动的总能耗。实际扑翼过程中存在一个从肌肉化学能到机械能再到有用功的转换效率问题，随着折转角度增加，翼型的运动幅度增加，结构自身加速、减速运动必然会造成额外的能量消耗，可能会导致能量转换效率的降低。此外，该定量的结论也仅针对本文简化的扑翼过程。在海鸥的实际飞行过程中，除了翼的折转运动之外，还存在明显的展向扭转运动和沿流向的收缩运动。在这些复合运动的作用下，上拍过程仍有可能产生升力，因此在实际的扑翼过程中折转角并不是越大越好。因此，在下文关于折转相位的分析中，根据文献观测数据^[31]，折转角度统一设定为θ_max=30°。

折转角相位差Δφ影响拍动过程中翼型展开-折叠的时机，其不仅影响到翼型有效迎风面积的变化，而且会影响到第二段翼的拍动速度。为分析其对翼型动态气动力特性的影响，针对Δφ=-60°~90°等若干状态进行了对比计算，计算中最大拍动角γ_max固定为30°。

Δφ＞0时，翼型展开-折叠的相位提前，在到达下拍中点之前翼型已经完全打开，而在到达上拍中点之前已经完全折叠，Δφ＜0时则相反。

图 10、图 11分别给出了拍动过程中的动态气动力系数和时均气动力系数。下拍过程中，随着Δφ从-60°~90°变化，升力峰、阻力峰、能耗峰值均提前，其原因在于翼完全打开的时间提前。

翼型绕根部的拍动角速度为：

$ \dot \gamma = {\gamma _{{\rm{max}}}}\omega {\rm{sin}}\varphi $

(6)

第二段翼的折转角速度为：

$ \dot \theta = - 0.5\omega {\theta _{{\rm{max}}}}{\rm{cos}}\left( {\varphi + \Delta \varphi } \right) $

(7)

当Δφ改变时，第二段翼的最大下拍速度会发生变化。例如当Δφ=0°时，在下拍过程的中间相位(φ=0.5π)时${\dot \gamma }$最大，${\dot \theta }$=0；而当Δφ=90°时，在下拍过程的中间相位时${\dot \gamma }$、${\dot \theta }$均为最大值。因此Δφ的变化必然会影响到气动力系数峰值的大小。下拍过程的时均阻力系数随Δφ的增加而单调减小，时均升力系数在Δφ≈30°时出现最大值。

上拍过程中，随着Δφ由负变为正，负的升力系数、阻力系数、能耗系数的峰值均单调增加，原因在于Δφ改变了上拍过程的最快拍动速度。从时均力学系数的变化曲线(图 11b)来看，相对“滞后”的折转角相位对上拍过程是有利的，可以减小阻力和“负升力”，并减少能耗。根据公式(7)可得，上拍过程中(φ=π~2π)，在φ=π~1.5π-Δφ区间内，第二段翼的折转角速度为正，因此可以减弱第二段翼的上拍速度。而随着Δφ的减小，这一区间是增大的，因此上拍的“不利”影响会得到减弱。

在本文所用模型和计算参数前提下，Δφ对于下拍和上拍过程的影响恰好相反，因此并没有对全周期的力学性能(图 11c)产生明显的有益影响。通过上述分析，折转运动最理想的情况应该是：下拍过程中，第二段翼向下折叠运动，以增益下拍速度，且在下拍速度最快的时刻保证较小的折叠角，以增加有效迎风面积；上拍过程中，第二段翼向下折叠运动，以减弱上拍速度，且保证上拍速度最快时折叠角较大，以减小有效迎风面积。通过文献[31]的观测数据可以看出，翼的折转角的变化规律不是简单的正弦曲线(图 12所示，psi2为折转角)，其折叠所用的时间远大于展开所用时间：在整个下拍过程，以及上拍过程的前半段，翼都处在折叠运动的过程中，在上拍过程快结束时再迅速展开，这一观测结果和本文的分析结果是相吻合的。

4 结论

本文通过对简化的海鸥扑翼过程进行数值模拟，分析了拍动角、折转角对非定常气动力以及流场结构的影响。拍动是扑翼的一个主要特征，随着拍动角度的增加，平均升力、推力均单调增加，同时能耗也明显增大；折转角则体现了翼的主动柔性变形，通过折转，改变上拍-下拍过程中翼的有效迎风面积以及第二段翼面的运动速度，对气动力影响较大。本文的数值模拟结果表明，折转角度主要影响上拍过程，尤其对能耗影响较大，而折转角相位对上拍及下拍过程的平均气动力均会产生明显影响。为了增益扑动过程的平均气动力特性，应当延长翼的折叠时间，减少翼的展开时间，这和实验观测的结果一致。

然而本文所得到的结论仅针对本文简化的扑翼过程，对于实际的扑翼过程，各种运动机制如拍动、扭转、收缩、折转是协同运作的，相互之间肯定存在明显的影响，单纯的将各个因素独立出来分析必然存在较大的近似。希望本文的工作能够为下一步针对更为真实扑翼过程的研究提供一些指引。