1. 引言
扑翼飞行器模仿自然界飞行生物的飞行方式,主要通过翅膀获得升力和推力[1]。与固定翼飞行器和旋翼飞行器相比,扑翼飞行器具有体积小、重量轻、隐蔽性强、机动性好、可垂直起降、悬停等诸多优点。该机未来可广泛应用于军事和民用领域。近年来扑翼飞行器受到越来越多研究者的关注[2]。
扑翼飞行器包含许多影响气动性能的参数,如翼型、翼面形状、扑动频率、下行程与运动周期之比。张焱等[3]使用Fluent与重叠网格技术研究了翼型在不同速度下与频率下的气动性能。结果表明,频率增大时,平均升力系数也变大平均阻力系数也增大。飞行生物能够有高的飞行效率与它们的翅膀结构有很大关系。刘旭博等人[4]构建了刚性翼与柔性翼模型,采用流固耦合技术,对扑翼流场进行了计算。结果表明,由于波动时域的存在,柔性翼产生了更大升力。张后伟等人[5]使用数值计算方法研究了单根羽毛弯角对羽毛气动性能的影响。研究发现流速不大时,在起始阶段,羽毛的升阻力波动较大。机翼上下摆动带动羽毛弯曲,从而产生升力。飞行生物通常做多自由度运动,通常包括上下扑动、俯仰和展向折,运动方式是影响气动性能的主要因素[6]。张兴等人[7]搭建了可以模拟来流状态和感知瞬时力矩的实验平台,研究了单段翼与双段翼的升力。结果表明,双段翼的升力约是单段翼的1.2~1.5倍。Lin等人[8]研究了NACA0012翼型在静水中的自主推进性能,分别为纯扑动运动、纯俯仰运动和扑动–俯仰运动。结果表明,扑动运动的自主推进速度和效率高于俯仰运动。
飞行生物的扑翼运动是不对称的,下行程的时间通常不等于上行程的时间。梁志宏等人[9]使用动网格技术分析了三维扑翼的气动力与力矩。结果表明,扑翼的升力随下扑速度增大而增加,调节上下扑动时间能够改变升力。谢鹏等人[10]通过变速因子使上下扑动时间不同,采用条带理论计算了非变速非折叠扑翼与变速折叠扑翼的气动性能。结果表明,扑动非对称的升力更大。Chen等人[11]对串列翼型模型进行了研究,串列翼型做上下扑动和不对称的阜阳运动。结果表明,串联翼应该保持同步俯仰,这样能获得更好的气动性能。这些研究分析了相位差固定时的不对称性,以及在不考虑对称性的情况下相位差对气动性能的影响。因此,有必要研究不对称和相位差共同作用时的气动性能。所以首先研究了扑翼在不同时间的下行程下的气动性能,又分析了相同下行程时间时相位差对气动性能的影响。
计算流体动力学(CFD)通过数值求解N-S方程来获得扑翼流场的物理量,并且还能考虑流场的粘性。张后伟等人[12]使用Fluent中线弹性实体光顺法研究了扑翼的气动性能,分析了升阻比、升力及阻力分别随流速、频率、攻角的变化规律。Zhu等人[13]使用计算流体力学方法研究了俯仰角对推力与升力的影响。Khanh Nguyen等人[14]通过计算流体力学研究了悬停机翼的运动方式对气动效率的影响,发现上下扑动与甩打运动结合时,能使升力增大5%。李佩等[15]对翼型周围流场进行了研究,发现扑翼运动使流场的速度发生变化,空气流动状态由层流变为湍流。重叠网格法将复杂的流动区域划分为具有简单几何边界的子区域,每个子区域中的计算网格独立生成。流场信息通过插值在重叠区域的边界处匹配和耦合。
扑翼下行程与运动周期之比为0.3~0.7。研究了时间不对称和相位差对扑翼气动性能的影响。研究对象为作扑动和俯仰运动的二维NACA0014翼型。根据曾锐等人[16]对相位差的研究,扑动和俯仰运动的
相位差分别为
和
。在曾的文章中,相位差为
时气动性能较好,相位差为
时气动性能较差。采用数值计算方法和重叠网格技术分析了扑翼的气动特性及扑翼周围的流场。
2. 翼型运动模型与气动参数
2.1. 运动模型
翅膀为飞行生物提供动力,运动是多维的,翅膀的运动可以用简谐函数来表示。扑翼运动不对称时,扑动的位移变化可以用分段函数表示,表达式如下:
(1)
其中
是最大扑动幅度,是一个无量纲数。
(2)
其中
是h的无量纲数,那么
代表扑动的幅度。L是特征长度并等于弦长c。T是运动周期,k是下行程与运动周期T的比值。通过改变k来调节扑翼的不对称运动。选择
。图1是扑动运动示意图,纵轴表示无量纲的位移,横轴表示出了弦长。图2显示了不同k时、在一个周期内扑动的位移,纵轴表示扑动运动的无量纲位移
。
俯仰运动的角度变化规律如下:
(3)
其中
是弧度制参数,
是俯仰运动的幅值,
是扑动和俯仰运动之间的相位差。当翼弦高于水平轴时,俯仰角定义为负值。研究了相位差为
和
的情况。图3的左侧表示相位差为
时俯仰角的变化,右侧表示相位差为
。
Figure 1. Plunging motion of two-dimensional NACA0014 airfoil
图1. NACA0014二维翼型的扑动运动
Figure 2. Changes of the dimensionless number h1 at different k
图2. 不同k时无量纲量h1的变化
Figure 3. Pitching motion of two-dimensional airfoil
图3. 二维翼型俯仰运动示意图
2.2. 气动参数
升力系数和阻力系数是反映气动性能的重要参数。扑翼的升力系数和阻力系数可定义如下:
(4)
(5)
其中L是升力,D是阻力,
是空气密度,S是翼型面积。因为数值计算中使用的翼型是二维翼型,所以S等于弦长。
和
分别是升力与阻力系数,都是与时间有关的物理量。
为了定量比较一个周期内的升力和阻力,平均升力系数和平均推力系数定义如下:
(6)
(7)
其中,
是平均升力系数,
是平均推力系数。均值可以定量地比较气动性能。
3. 数值计算方法
3.1. 控制方程与求解方法
扑翼飞行的流场属于低速不可压缩非定常粘性流场,雷诺数较小。在笛卡尔坐标系中,控制方程可由纳维-斯托克斯方程(N-S方程)表示,包含不可压缩连续方程和动量方程。求解算法为压力和速度耦合。采用Fluent软件对控制方程进行求解,结合重叠网格技术得到了瞬时流场的物理量。
3.2. 网格划分和边界条件
Figure 4. Grid and boundary conditions
图4. 网格与边界条件示意图
图4为在ICEM建立的包括翼型和背景网格的网格系统。翼型为NACA0014,忽略翼型变形。为了更好地捕捉翼型周围的流体流动,翼型周围的网格被加密。第一层网格的厚度为0.001 mm。计算区域的左侧边界为速度入口边界条件。速度值代表向前飞行的速度。右侧边界为压力出口边界条件,上下边界是对称边界条件,翼型是壁面边界。采用重叠网格技术和UDF处理翼型运动。翼型周围的网格是组分网格,它在背景网格平面内移动。数值计算时首先删除背景网格中组分网格对应的区域,然后通过相邻网格单元的插值实现组分网格与背景网格之间的数据传输。
当使用ICEM中的distance作为网格质量的标准时,最小网格质量为0.974;当the Aspect Ratio作为标准时,网格的最低质量为1。网格质量检验结果表明,网格精度较高。
3.3. 网格灵敏性和时间步长独立性验证
在减缩频率
、
、
、
、
的条件下,计算了三种不同网格系统的非定常流场,以检验网格灵敏性。减缩频率表示为:
(8)
其中
是来流速度,也是翼型向前飞行的速度。
测试网格灵敏性的网格系统为100 × 100节点(垂直翼弦方向100个,沿翼型弦向100个)、100 × 200节点、200 × 200节点,靠近翼型的第一层网格厚度为0.001 mm。每个网格系统模拟四个周期。每个周期中的时间步长数为500。机翼只做扑动运动。表1总结了每个网格系统的组分网格中包含的节点和网格的数量。B的网格在弦向的长度是A的两倍,在垂直弦向的长度相同。C的网格在垂直弦向的长度是B的两倍,弦向的长度相同。
Table 1. Number of grid nodes and total number of grids in grid system
表1. 网格系统的网格节点数与网格总数
网格系统 |
节点数 |
组分网格数量 |
A |
100 × 100 |
34,404 |
B |
100 × 200 |
70,101 |
C |
200 × 200 |
134,801 |
图5为不同网格系统在一个周期内的阻力系数变化,纵轴为阻力系数,横轴为时间。结果表明,网格数对翼型的阻力系数影响不大。在行程转换时,阻力系数略有不同。因为B具有较高的计算效率和精度,本研究采用网格B进行模拟。
图6为网格系统B在时间步长分别为
和
时的升力系数,纵轴表示升力系数,横轴表示时间。验证时间步长无关性的计算条件与验证网格灵敏性的一致。结果表明两个时间步长的升力系数基本相同。较小的时间步长将减少计算时间,因此后期研究采用时间步长
的方案。
3.4. 数值计算方法有效性验证
有必要验证数值计算方法在分析非定常流问题时的准确性。本文用Miao [17]中的结果验证了本文分析方法的准确性。在数值模拟中,计算条件与2.3节中的一致.迭代时间步长为
,网格系统为B。h从下方穿过0线,表示翼型从下形程起始位置开始运动。图7中展示了翼型位置的变化、Miao文中的结果和本文使用重叠网格法的数值结果,纵轴为阻力系数,横轴为时间。
通过比较发现,前四个周期的计算结果与苗的计算结果一致,数值也基本一致。因此,本文的数值计算方法可以有效地分析上下行程时间比和相位差对翼型气动特性的影响。
Figure 5. Drag coefficient of different grid systems
图5. 不同网格系统的阻力系数
Figure 6. Lift coefficient of different time step
图6. 网格系统B在不同时间步的升力系数
Figure 7. Accuracy verification of simulation results
图7. 数值方法准确性验证
4. 结果和讨论
4.1. 非对称运动对翼型气动性能的影响
在下行程与运动周期之比的研究中,扑动运动幅值
为0.4,弦长c为0.1 m,马赫数
为0.01,减缩频率
,俯仰幅值
,相位差
,雷诺数
。迭代时间步长为
,网格
系统为B。图8和图9为
时瞬时阻力系数和瞬时升力系数随时间的变化曲线。在第一个运动周期内,翼型的升阻系数的变化还不稳定,从第二个循环开始,升阻系数的变化基本稳定。升力系数和阻力系数分别在上行程和下行程时有峰值。在上一节定义的翼型运动中,翼型先进行下扑运动,然后进行上扑运动。负阻力系数表示产生了推力。通过阻力系数曲线发现,随着k的增加,下行程时翼型阻力系数的峰值变大,上行程时阻力系数的峰值变小。比较
的阻力系数与
的阻力系数发现,下行程中的推力峰值变小,但上行程中的推力峰值变大。当
时,一个周期中下行程时间长于上行程时间,因此上行程速度较大。将
的阻力系数与
的阻力系数进行比较,下行程期间的推力系数峰值变大,但上行程期间的推力系数峰值变小。当
时,一个周期中的下行程时间短,上行程时间长,则下行程速度大。
从图9可以发现,下行程的升力系数为正,上行程的升力系数为负。随着k的增加,下行程翼型升力系数的峰值变小。上行程升力系数的峰值变小,此时,翼型产生负升力。通过比较
和
的升力系数发现,前者在下行程期间具有较大的升力系数峰值,而在上行程期间升力系数的负峰值较小。
时下行程中升力系数的峰值略高于
的,但上行程中两个峰值相差不大。与
相比,下行程的
时的升力系数峰值变小,上行程的升力系数峰值也变小,这表明上行程的负升力变大。
Figure 8. The curve of instantaneous drag coefficient in three periods with different k,
图8. 不同k时在三个周期的瞬时阻力系数
Figure 9. The curve of instantaneous lift coefficient in three periods with different k,
图9. 不同k时在三个周期的瞬时升力系数
为了定量分析非对称运动对扑翼气动性能的影响,分别计算了一个运动周期内升力系数和阻力系数的平均值。选择的周期是数值计算的第三个周期。图10的左纵轴表示平均升力系数,右纵轴表示平均阻力系数。当下行程时间小于上行程时间时,
时的平均升力系数与
时的平均升力系数相比显著提高,但
时的平均升力系数小于
时的平均升力系数。当下行程时间大于上行程时间时,非对称运动的平均升力系数均小于对称运动的平均升力系数,此时平均升力系数为负。当
时,平均阻力系数接近于0,小于
时的平均阻力系数,表明翼型产生的阻力较小。当
时,平均阻力系数小于0,表明翼型在整个运动周期中都产生推力。然而,
的平均升力系数最低。
时的平均阻力系数大于对称运动时的阻力系数,而
时的平均推力系数与
时相比有显著提高。综合平均升力系数和平均阻力系数,当翼型做非对称运动且下行时间足够短时,平均阻力系数和平均升力系数都较好。当翼型做时间不对称运动时,上行程时间长可以增加平均推力系数,但平均升力系数为负。
Figure 10. Mean value of lift and drag coefficient at different values of k,
图10. 不同k时一个周期的平均阻升力系数,
4.2. 不同相位差下非对称运动扑翼的气动性能
研究翼型做相同非对称运动时相位差
和
对扑翼气动性能的影响。在模拟中,扑翼做扑动运动和俯仰运动,其他模拟条件与4.1节相同。
从图9和图11可以看出,当相位差
相同时,随着k值的增大,升力系数的变化规律基本相同。在这两种情况下,下行程和上行程的升力系数的峰值均减小,上行程的升力系数的峰值为负。当相位差
不同时,通过比较k相同时的峰值升力系数发现,
的峰值升力系数大于
的峰值升力系数。在翼型上行程期间,
的峰值升力系数小于
的峰值升力系数。比较图8和图12可以发现,在机翼下行程和上行程期间,
的阻力峰值大于
的阻力峰值。阻力系数的峰值都是负值,表明翼型产生推力。
为了定量比较
与
在相同k下的气动性能,图13给出了
在不同k下的平均升阻系数。计算周期为翼型运动的第三个周期。表2详细显示了一个周期的平均升力系数和平均阻力系数。当
和
时,翼型的平均升力系数低于
时的平均升力系数。但是平均推力系数更高,这意味着减少下行程时间会降低平均升力系数并增加平均推力系数。
和
时翼型的平均升力系数和平均推力系数高于
时的翼型。此外,
具有更长的下行程时间。平均推力系数增大,但平均升力系数减小。通过比较相同k下的升力系数发现,仅在
时,
的平均升力系数比
的平均升力系数有显著增加,而在k的其他值下,机翼的平均升力系数都降低了。从表2中还可以发现,当翼型做非对称运动时,推力系数显著增加,不对称度越高,推力系数越大。
Table 2. Mean value of lift and drag coefficient at different values of k,
表2. 不同k与
时一个周期内的平均升力系数与阻力系数
下行程与整个
行程比值k |
平均升力系数
,
相位差
|
平均阻力系数
,
相位差
|
平均升力系数
,
相位差
|
平均阻力系数
,
相位差
|
0.3 |
0.04636 |
0.04636 |
0.70557 |
−0.21718 |
0.4 |
0.10394 |
−0.04434 |
−0.37882 |
−0.17023 |
0.5 |
0.15672 |
−0.03771 |
0.07313 |
−0.18463 |
0.6 |
0.19599 |
−0.04316 |
−0.21477 |
−0.21477 |
0.7 |
0.15886 |
−0.06087 |
−0.57834 |
−0.22542 |
Figure 11. The curve of instantaneous lift coefficient in three periods with different k values,
图11. 不同k时在三个周期的瞬时升力系数,
Figure 12. The curve of instantaneous drag coefficient in three periods with different k values,
图12. 不同k时在三个周期的瞬时阻力系数,
Figure 13. Mean value of lift and drag coefficient at different values of k,
图13. 不同k时一个周期的平均升阻力系数,
4.3. 翼型的流场分析
图14展示了在
和
时,NACA0012翼型在一个运动周期内的压力场变化。翼型做扑动运
动和俯仰运动。在图中,h代表翼型扑动运动的无量纲值。当翼型的位置从
变化到0时,下翼型表面的负压中心区域逐渐脱落,下翼型表面存在一个较大的正压区域。在此过程中,下翼面压力大于上翼面压力。当
时,上翼面形成负压中心区,负压中心区逐渐向翼型后缘移动,
时即将脱落。从
到
,下翼面的压力总是大于上翼面。因此,在整个下行程期间,下翼型表面的压力大于上翼型表面的压力,这也表明下行程产生升力。上翼型负压中心区域出现的时刻对应于升力系数达到正峰值的时刻。
表2表明,
和
的平均升力系数大于
和
的平均升力系数,这表明下行程时间过长不利于增加平均升力系数,尽管它会延长上下翼面之间的压差时间。在上行程的
期间,上翼面的负压区脱落,下翼面的负压区逐渐扩大。在
到
期间,负压中心区总是位于下翼面,并随着上
行程向后缘移动。在上行程期间,上翼面压力大于下翼面压力,这由负升力表示。在
时,上下翼面之间的压差最大,对应于升力系数的负峰值。
Figure 14. Pressure contour of NACA0012 in one cycle of motion when
and
图14. NACA0012在一个周期的压力云图变化,
,
5. 结论
本文建立了时间非对称的扑动和俯仰运动模型,采用数值计算方法和重叠网格技术分析了上下行程时间比和相位差对气动性能的影响。主要结论如下:
(1) 当
时,延长下行程时间能提高平均升力系数。下行程时间过长会延长上下翼型之间的压差时间,但不利于增加一个运动周期的平均升力系数。
(2) 在相位差
时,平均升力系数在
时显著增加。
(3) 在整个下行程期间,下翼型表面的压力大于上翼型表面的压力,这也表明下行程产生升力。在上行程期间,上翼面压力大于下翼面压力,这表示产生负升力。
(4) 翼型做不对称运动时,推力系数显著增加,不对称程度越高,推力系数越大。