拍动迎角对扑翼悬停飞行气动特性影响的数值模拟研究

Numerical Simulation on the Effect of Angle of Attack on Aerodynamic Forces about Flapping Wing in Hovering Motion

摘要: 为了研究悬停状态下扑翼飞行的运动参数对其所产生的非定常气动力影响，本文通过数值模拟方法模拟了不同拍动迎角条件下以提前翻转模式往复拍动二维平板的流场演化特性。通过比较不同迎角算例的瞬时与平均气动力，得到了不同迎角下的气动力变化规律。在拍动迎角增加的过程中，总共有四个因素受到迎角的影响，即尾迹捕捉机制，翻转机制，失速延迟机制，拍动迎角对气动力的升力分量的影响。在拍动行程开始的加速阶段，气动力主要受平板–涡相互作用的影响。而在匀速平动阶段，大迎角条件下的拍动由于受附着前缘涡的作用，可产生较大的升力，但若继续增大迎角会导致气动力的升力分量减小。在最后的减速–翻转阶段，由翻转机制提供的气动力随着迎角的增加而减小，这主要是因为大迎角拍动的平板对应着更小的转动的角速度，导致附加环减小，从而降低升力分量。综合上述各种机制的影响，对于所研究的拍动模式，扑翼的时均升力系数在拍动迎角为55度附近存在一个极大值，进一步增加迎角会导致升力的减小。

Abstract: In order to study the effect of kinematic parameters on aerodynamic force of flapping wing in hovering motion, the evolution of flow fields about a plate with advanced rotation under different angle of attack is numerical investigated. The result shows the variation of instantaneous lift coefficient and time averaged lift in different angle of attack. There are four mechanism changes by angle of attack, wake capture, rapid pitch rotation, stall delay, and the influence of angle of attack on the component of aerodynamic force on lift. In the beginning of the stroke, the aerodynamic force is mainly affected by plate-vortex interaction. In the constant velocity phase, the plate at high angle of attack can generate high lift by the leading edge vortex. However, a higher angle of attack can also reduce the component of aerodynamic force on lift. At the rotation-declaration phase, as the rotational velocity in the case with high angle of attack is lower than that of the case with small angle of attack, the added circulation and lift are also reduced. Considering the influence of the above mechanisms in the flapping motion present studied, there is a maximum value of the time-averaged lift coefficient of flapping wing near the angle of attack of 55 degrees, and further increase of the angle of attack will lead to the decrease of lift.

1. 引言

近年来微型飞行器(MAV)以其尺寸小，成本低，隐蔽性好等特点在信息搜集，环境监测和军事侦查领域有着广泛的应用，成为飞行器设计领域的新热点。但传统的固定翼飞行模式在如此低雷诺数条件下已不再适用。受到与其尺寸相当的动物飞行(昆虫、小尺寸的鸟或蝙蝠等)的启发，人们对扑翼飞行模式在MAV上的应用产生了浓厚的兴趣。扑翼飞行的广泛地存在于自然界中，比如大量的鸟类和昆虫都采用了扑翼飞行的方式 [1] [2]。昆虫可以通过扑翼扑动平衡自身的重力，蜂鸟的扑翼飞行中可以实现悬停和倒飞。这种具有出色气动性能的飞行模式引起了大量研究者的关注，针对昆虫和飞鸟的扑翼高升力机理，国内外开展了不少研究 [2]。

目前发现的升力机制可以分为非定常升力机制(失速延迟机制 [3] )和准定常升力机制(尾迹捕捉机制 [4]，翻转机制 [4]，附加质量机制 [5] )。其中，非定常升力机制比准定常升力机制更加重要。因为非定常特性是扑翼飞行的主要特性之一，并且扑翼行程结束的转动阶段产生的非定常升力比扑翼行程中间阶段的准定常升力更高 [4]。而在所有非定常升力机制中，尾迹捕捉机制是目前仍然有待探索机制之一。在尾迹捕捉机制提出的早期甚至存在一些争议。在Dickinson等人 [4] 提出尾迹捕捉机制之后，Sun等人 [6] 的数值模拟研究结果表明尾迹捕捉反而会降低气动力。直到Lua等人 [7] 通过二维实验观察了不同迎角下往复拍动的平板，发现了图1中的两类不同涡和平板相互作用。当尾迹涡对的诱导速度垂直于平板时，尾迹捕捉才会提高升力。Lua等人 [7] 的结果揭示了尾迹捕捉机制会在不同情况下带来不同的作用，初步解决了早期关于尾迹捕捉机制的争议。然而Lua等人 [7] 的研究中平板以固定迎角往复拍动，并没有以真实扑翼飞行的平动加转动的方式进行拍动。所以，在真实的扑翼拍动状况下，怎样的运动参数能够产生第一类尾迹涡–平板相互作用还有待进一步研究。

Figure 1. Schematic of two kinds of wake-wing interaction [7]. (a) The first kind of wake-wing interaction; (b) The second kind of wake-wing interaction

图1. 两类尾迹涡–平板相互作用示意图 [7]。(a) 第一类尾迹涡–平板相互作用；(b)第二类尾迹涡–平板相互作用

本文将通过二维数值模拟的方法，以梯形分段函数作为扑翼实际悬停飞行时的运动模型，对不同拍动迎角条件下平板的气动力特性与尾迹流场结构进行研究。通过进一步分析尾迹涡量场演化与平板表面压强分布的关系，得出产生两类不同尾迹涡–平板相互作用的运动参数条件，并解释其对气动特性的影响规律。

2. 计算模型描述与运动参数

在本文的工作中也主要以提前转动 [4] 作为扑翼的运动模式，考虑二维的悬停飞行，并且拍动平面在二维条件下简化为水平直线。如图2，主要考虑平板沿着扑动方向平动和绕着转轴的转动。地面坐标系的坐标原点固定在开始运动前的转轴处。主要运动参数包括转轴的位移和平板绕转轴转动的角度。扑翼转动转轴位置在距离前缘0.25弦长的位置。在本文中，扑翼被简化为一个有厚度的刚性平板，弦长为0.06 m，平板厚度为弦长的1/20。

Figure 2. Schematic of two-dimensional motion of wing. The direction of the velocity is given by red arrows (translation) and black arrows (rotation). The leading edge is noted by yellow point. The “Stroke L” and “Stroke Rare” is the left stroke and the right stroke, respectively

图2. 二维扑翼运动示意图。速度的方向由红色箭头(平动)和黑色箭头(转动)给出。前缘被黄色圆点标出。“Stroke L”和“Stroke R”分别是向左行程和向右行程

在本文中假设平板的平动速度和转动速度函数满足提前转动的梯形分段模型 [4]：

$u\left(t\right)=\{\begin{array}{l}-U_0\sin \left(\frac{\pi t}{2t_a}\right)t\in \left[0,t_a\right)\hfill \\ -U_{0}t\in \left[t_a,\frac{T}{2}-t_a\right)\hfill \\ U_0\sin \left[\frac{\pi \left(t-\frac{T}{2}\right)}{2t_a}\right]t\in \left[\frac{T}{2}-t_a,\frac{T}{2}+t_a\right)\hfill \\ U_{0}t\in \left[\frac{T}{2}+t_a,T-t_a\right)\hfill \\ -U_0\sin \left[\frac{\pi \left(t-T\right)}{2t_a}\right]t\in \left[T-t_a,T\right]\hfill \end{array}$ (1)

$\omega \left(t\right)=\{\begin{array}{l}0t\in \left[0,\frac{T}{2}-t_r\right)\hfill \\ \frac{{\omega }_0}{2}\sin \left[\frac{\pi }{t_r}\left(t-\frac{T}{2}+t_r\right)\right]-\frac{{\omega }_0}{2}t\in \left[\frac{T}{2}-t_r,\frac{T}{2}\right)\hfill \\ 0t\in \left[\frac{T}{2},T-t_r\right)\hfill \\ -\frac{{\omega }_0}{2}\sin \left[\frac{\pi }{t_r}\left(t-T+t_r\right)\right]+\frac{{\omega }_0}{2}t\in \left[T-t_r,T\right)\hfill \end{array}$ (2)

其中u是平动速度， $\omega$ 是转动速度， $U_0=0.011\text{m}/\text{s}$ 是最大平动速度， ${\omega }_0$ 是最大转动速度， $T=0.591\text{s}$ 是扑翼拍动的周期， $t_a=0.5\text{s}$ 是加速或者减速段的时间段， $t_r=1\text{s}$ 是转动的时间。在本文中角速度 ${\omega }_0$ 将随着转动幅度的增加而等比例增加。在图3中以迎角为60˚的算例为例子，给出了平动和转动速度随时间的变化的示意图。

Figure 3. The translational velocity and rotational velocity of the flapping wing

图3. 扑翼的运动速度随时间变化示意图

本文中的流体参数与Dickinson等人 [4] 的实验中相同运动粘度为 $\upsilon =1.15\times {10}^{-4}{\text{m}}^2/\text{s}$，密度为 $\rho =880\text{kg}/{\text{m}}^{\text{3}}$，使得雷诺数为 $\mathrm{Re}=U_{0}c/\upsilon =57.4$。这一设置在一般扑翼飞行的雷诺数范围 $O\left({10}^1\right)\text{-}O\left({10}^4\right)$ 之内 [8]。

在本文中设置了迎角从15˚增加到75˚的共6个算例，这6个算例的参数在表1中列出。在真实扑翼扑动飞行的过程中，迎角的范围是25˚~60˚ [8] [9]，本文所设置的算例参数能够覆盖此范围。

Table 1. Angle of attack and rotational velocity of different cases

表1. 不同算例的迎角和转动速度

3. 数值计算方法

3.1. 控制方程与求解

考虑到本文中算例的雷诺数为57.4，因此计算模型选择非定常、二维、层流、不可压的N-S方程，方程的描述如下：

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$ (3)

$\frac{\partial \left(\rho u_i\right)}{\partial t}+\nabla \cdot \left(\rho u{u}_i\right)=\frac{\partial {\tau }_{ij}}{\partial x_j}+F_i$ (4)

其中 $\rho$ 是流体密度， $u_i$ 是流体速度的分量， $u$ 是流体速度矢量， $F_i$ 为体积力， ${\tau }_{ij}$ 是包括静压和粘性摩擦力在内的表面力。主要计算方法采用有限体积法。空间离散化方案中，梯度使用基于单元的最小二乘方法，压力采用二阶方法，动量采用二阶迎风方法。瞬态离散化方案采用一阶显式方法。

3.2. 网格划分与计算设置

本文的计算基于方法重叠网格技术。在重叠网格的计算过程中，将计算域分为前景网格和背景网格。其中静止不动的部分为背景网格，移动的平板和平板附近的一定区域被划分为背景网格。图4给出了前景网格和背景网格的示意图。图中平板位于背景网格中间，且初始时刻的转轴位置为坐标原点。前景网格中，平板周围区域的边界被称为重叠网格边界。重叠网格边界可以任意选取，因此可以通过选取合适的重叠网格边界以保证较高的网格质量。在本文中，为了保证前景网格的正交性，将重叠网格边界设置为平板边界向外偏移R所形成的形状。背景网格则设置为L × W的长方形。在本文中 $R=8.37c,L=21.5c,W=38c$。前景网格和背景网格分别使用了分块结构化网格。在背景网格当中，L和W两个方向的网格节点数量157 × 385，前景网格的壁面上和垂直于壁面上的网格数量为596 × 394。网格总数一共为268,449。

背景网格边界处的边界条件为压力出口，且压力参考点在背景网格边界上。前景网格的平板壁面设置为无滑移壁面平板周围的边界设置为重叠网格边界。在所有算例中，时间步长设置为常数 $\Delta t=0.01\text{s}$。为了使得平板拍动运动达到稳定的周期性状态，总计算时长设置为8T，即平板往复拍动8个周期的时间。

Figure 4. Schematic of the computational domain, fore-ground domain and background domain

图4. 计算域及前景区域和背景区域示意图

3.3. 数值方法验证

为了验证计算设置的合理性，本文将对实验设置进行验证性模拟。本文的验证模拟是和Wang等人 [10] 的结果进行对比。为了便于对比，验证模拟的拍动函数、流体参数等设置与Wang等人 [10] 的研究中的对称转动一致。

Figure 5. Comparison between C_L curves of the validation case and the result of Wang et al.

图5. 验证算例与Wang等人结果的升力系数曲线对比

Figure 6. Comparison between vorticity contour of the validation case and the result of Wang et al. A and C are computational result of Wang et al. B and D are validation case’s vorticity contours

图6. 验证算例和Wang等人结果中的涡量云图对比。A和C为Wang等人的计算结果。B和D为验证性算例计算结果

图5和图6中分别表现了本文的验证模拟结果与Wang等人 [10] 的结果中升力系数和涡结构结果的对比。结果表明本中的计算设置能够还原真实流场中的涡结构并准确的计算出平板的气动力。

4. 计算结果及讨论

本文的研究目的是得到在不同迎角下的尾迹涡和平板相互作用的规律。因此将主要讨论3.1中给出的6个算例，讨论迎角的变化对尾迹涡和气动力的影响。在4.1中首先得到时均气动力随迎角的变化规律，分析了单次拍动行程中的气动力特征。在4.2中分析了不同迎角下的涡结构随时间变化的规律，通过压强分布运动证明了涡结构对气动力的影响，从而得到了迎角对尾迹捕捉的影响。

4.1. 平板的瞬时与平均气动力特性

在本文中升力系数的定义如下：

$C_L=\frac{F_L}{\frac{1}{2}\rho U_0^{2}c}$ (5)

其中 $F_L$ 是平板所受到的升力，升力通过计算平板上压强和粘性切应力在升力方向上的积分得到。时均升力则是升力系数在第8个周期的时均值：

$\overline{C_L}={\displaystyle {\int }_7^8{C}_L\text{d}{t}^{\text{*}}}$ (6)

其中 $t^{\text{*}}=t/T$，是用往复拍动周期进行无量纲化处理的时间。不同拍动迎角的时均升力系数在图7中给出。在图7中可以看出在迎角不超过55˚时，时均升力系数随着迎角的增加而增加。在迎角为15˚时，时均升力系数最低，为0.17，在迎角为55˚时，时均升力系数最高，达到了1.39。而迎角从55˚增加到75˚时，时均升力系数从1.39下降到了0.80。这说明迎角为55˚时能产生最大的时均升力。

Figure 7. Time-averaged C_L of cases with different angles of attack

图7. 不同迎角算例的时均升力系数

为了进一步研究升力系数随迎角变化的规律，图8中给出了第8次往复拍动时的瞬时升力系数曲线。在图8中将迎角大于55˚的算例和迎角小于55˚的算例分别分组比较。并且标注了10个不同的特征时刻。在加速阶段，迎角大于或者等于55˚的算例会出现一个峰值，这个峰值随着迎角的减小而降低。而迎角小于55˚的算例，在这一阶段升力系数不会出现峰值。在迎角45˚的算例中，加速段的升力系数曲线近似为直线。在更小的迎角中，加速阶段的升力系数曲线出现下凹。在加速阶段中，气动力机制主要是平板加速导致的附加质量效应和尾迹捕捉机制。因为对所有的算例，加速阶段的速度随时间的变化都保持相同的变化规律，所以导致加速阶段的气动力出现差异的原因可能是尾迹捕捉。这一假设在4.2中将对涡结构的演化是进行分析时进行检验。在匀速阶段，升力系数曲线表现出准定常特性，即随时间变化幅度很小。这一阶段的升力系数主要受到失速延迟机制 [2] [3] 和迎角对气动力的升力分量影响。根据Li等人 [11] 的对二维大迎角平板在无粘自由来流中的涡动力学分析，气动力近似保持不变的原因是前缘涡提供的气动力降低和后缘涡提供的气动力增加，因此总的气动力保持了平衡。当拍动迎角从75˚减小到55˚时，准定常升力系数随着迎角的减小而增加，这是迎角减小导致的升力系数分量增加。但是当迎角从55˚进一步减小到15˚时，准定常升力系数随着迎角的减小而降低。在转动阶段，升力系数曲线会首先出现一个峰值。这个峰值时平板快速转动导致的升力。迎角越小的算例在相同的时间内需要转过的角度越大，因此转动速度也就越大，从而能够产生更大的峰值。当平板转动到竖直位置时，平板的法向没有竖直方向的分量，因此升力系数几乎为0。当平板进一步转动时，平板开始减速，气动力回升到行程开始时的大小。

(a)(b)

Figure 8. C_L curves of cases with different angle of attack. (a) Cases with angle of attack within 15˚~55˚; (b) Cases with angle of attack within 55˚~75˚. Ten representative instants t₁~t₁₀ are depicted in this figure

图8. 不同迎角算例的升力系数曲线。(a) 15˚~55˚拍动迎角算例的瞬时升力系数；(b) 55˚~75˚拍动迎角算例的瞬时升力系数。图中标注了t₁~t₁₀的10个特征时刻

通过以上气动力分析，分别得出了三个运动阶段内气动力曲线随迎角的变化规律。在行程开始的加速阶段的峰值随着拍动迎角的减小而减小。在匀速运动阶段，迎角大于55˚时，迎角减小会使得气动力的升力分量增加，导致这一阶段升力系数增加。迎角小于55˚，迎角减小会导致匀速升力的减小。在转动阶段，拍动迎角减小会引起转动速度增加，进一步使得升力系数增加。但是以上规律中，行程开始时的尾迹捕捉，行程中匀速阶段的前缘涡引起的失速延迟等现象仍然需要通过涡结构的观察以进一步确认。

4.2. 尾迹涡结构对表面压强分布的影响

尾迹捕捉，失速延迟等机制都和涡结构有密切关系。Lua等人 [7] 指出的两类尾迹涡–平板相互作用也是由尾迹涡对和平板的相对位置决定的。因此对涡结构进行观察是非常必要的。在图9中展示了三个不同拍动迎角的涡量和压强云图，确认了在行程的不同阶段涡结构的演化和对气动力的作用。

(a) (b) (c)

Figure 9. Pressure (left) and vorticity (right) contours at ten instants during the eighth stroke for the cases with different angle of attack. (a) The case with α = 60˚; (b) The case with α = 55˚; (c) The case with α = 45˚. LEV is leading edge vortex and TEV is trailing edge vortex. The subscript L (or R) here indicates that the vortex is formed during the left stroke (or the right stroke). The index number in the subscript corresponds to the order of the vortex formation in each stroke

图9. 不同迎角的算例在第8个拍动行程内，10个特征时刻的压强(左侧)和涡量云图(右侧)。(a) α = 60˚的算例。(b) α = 55˚的算例。(c) α = 45˚的算例。在云图中负值的等值线由虚线标记。LEV代表前缘涡，TEV代表后缘涡，下标中的L(或R)代表向左的行程(或向右的行程)中产生的涡。数字下标对应每个行程中涡的产生的序号

在扑动行程开始的时刻，平板的气动力主要受到尾迹捕捉的影响。在对于迎角始终保持为90˚的平板，流动结构具有对称性。两个尾迹涡的位置关于平板的中垂线对称，尾迹涡诱导的速度只有垂直于平板的分量，此时发生的尾迹捕捉对应第一类尾迹涡–平板相互作用。当拍动迎角减小到(55˚，90˚)的范围内，对称性被破坏。但是来自于上个扑动行程的前缘涡LEV_L1和后缘涡TEV_L2的诱导速度在垂直于平板的方向上仍然有较大的分量，从而诱导流体撞击平板产生较大的气动力。在这一时刻的压强分布云图中，平板中央会出现高压区。在上一节的升力系数曲线上，这一时刻的和升力系数曲线出现了第一个峰值。这种尾迹捕捉过程依然对应第一类尾迹涡–平板相互作用。当拍动迎角进一步减小到(0˚，55˚)的范围内，涡对的诱导速度在垂直于平板的分量进一步降低，这对应第二类尾迹涡–平板相互作用。并且此时尾迹涡对和平板的相对位置也发生了变化，上个扑动行程的前缘涡会LEV_L1距离平板近，而后缘涡TEV_L2距离平板较远。上个扑动行程的前缘涡LEV_L1会在平板的迎风面产生低压区，降低了气动力。这种现象在Lua等人的研究被称为前缘涡的滑移 [12]。

在匀速平动阶段，平板的气动力和前缘涡引起的失速延迟机制有关。通过比较涡量云图可以看出迎角越小的算例，前缘涡形成的越慢。大迎角的平板能够更快地卷起前缘涡。而迎角为45˚的小迎角的平板直到开始转动前，前缘涡都没有形成。在匀速平动阶段，大迎角的前缘涡在上表面形成的低压区，导致这一阶段的气动力大于小迎角的算例。因此迎角越小，前缘涡形成的越晚，失速延迟机制越弱，升力越低。但是对于拍动迎角超过55˚的算例，迎角过大也导致气动力的升力分量减小。

在转动阶段，气动力主要受到翻转机制的影响。在转动阶段初期，t₆前后，较小迎角的算例能够产生较大的转动速度，从而产生更高的气动力峰值。在转动阶段t₇和t₈，涡量和压强分布云图表明在转动到垂直位置之后的一段时间之内，前缘涡LEV_R1和后缘涡TEV_R2产生的低压区使得平板出现负升力。这个负升力是t₈时刻的升力系数曲线出现谷值的原因。之后平板逐渐减速，前缘涡不再吸收涡量。随着涡量的耗散，涡的低压区产生的负升力逐渐减小，因此气动力出现回升。在此之后，大迎角算例中的尾迹涡已经在t₉时产生了第一类尾迹涡–平板相互作用，这进一步增强了升力。而低迎角的算例，随着迎角不断减小，平板受到的气动力在竖直方向的分量不断增加，因此负升力也逐渐增加。由于迎角较小，此时不会出现第一类尾迹涡–平板相互作用，因此升力系数会逐渐降低。

涡结构在平板拍动过程的演化验证了尾迹捕捉机制，失速延迟机制随着拍动迎角变化的规律。在迎角90˚时会发生第一类尾迹涡–平板相互作用。随着迎角减小到55˚附近的临界迎角，尾迹涡对和平板的相对位置发生了改变，使得尾迹涡对的诱导速度平行于平板。这使得尾迹捕捉转化为第二类尾迹涡–平板相互作用，进一步导致在行程开始的加速阶段气动力降低。当迎角减小时，导致前缘涡形成缓慢，从而在匀速运动阶段失速延迟机制产生的气动力较弱。当迎角过大时，迎角的减小也会导致气动力的升力分量增加，从而引起加速和匀速阶段的气动力的增加。

5. 结论

通过对不同迎角的二维提前转动平板了数值模拟研究。在拍动行程开始阶段，发现迎角是影响尾迹捕捉类型的因素之一。拍动迎角超过55˚附近的临界迎角，可以使得尾迹涡产生第一类尾迹涡–平板相互作用，从而增加行程开始阶段的升力系数。然而迎角过大时，迎角的减小也能导致气动力的升力分量增加，从而提高升力。在匀速阶段，迎角的减小会导致前缘涡的行程变慢，从而减弱失速延迟机制，进一步导致气动力降低。通过比较升力系数曲线验证了翻转机制，转动速度较快的算例能够在转动阶段产生较大的升力。

但是运动参数的影响较为复杂，除了迎角这一变量之外，还有拍动行程，拍动平面与水平面的夹角等运动参数。当其他运动参数发生变化时，尾迹涡和平板的相互作用类型会如何变化还需要进一步研究。

NOTES

^*通讯作者。