仿蝠鲼柔性潜水器翼型流场性能分析

邢 城，潘 光，黄桥高

（1. 西北工业大学 航海学院，陕西 西安 710072；2. 无人水下运载技术工信部重点实验室，陕西 西安 710072）

0 引言

“物竞天择，适者生存”蝠鲼经过长时间的进化与自然选择使得其本身具有优异的推进效能。影响蝠鲼推进效能的主要因素是蝠鲼的外形轮廓，而剖面轮廓作为蝠鲼外形轮廓的重要组成部分，其对蝠鲼运动时的流场特性起着决定性的作用。所以展开蝠鲼剖面轮廓形状对其运动特性的的影响是十分必要的。

国内外对蝠鲼已经进行了大量的研究，日本的研究人员通过试验研究了不同鳍的形状对仿生蝠鲼推进性能的影响[1-2]，试验结果表明不同形状的鳍对机器鱼的推进速度有着很大影响，并且通过对胸鳍上不同的行波参数下的推进性能比较，研究了行波参数和胸鳍推进速度间的关系。在国内方面，国防科技大学从仿生学层面较为系统地开展了牛鼻鲼胸鳍升力推进模式的仿生研究[3-4]，不仅基于流固耦合的方法分析了牛鼻鲼胸鳍推进模式的水动力性能，并设计和研制了仿牛鼻鲼潜水器，通过样机实航试验研究了不同的胸鳍摆动参数下潜水器的推进性能。北京航空航天大学机器人所在仿蝠鲼潜水器样机研制方面开展了大量的工作[5-7]。章永华等人[8]针对二维翼型的2种推进模式（振幅从前往后保持不变；振幅从前往后逐渐增加）进行了比较。研究表明在运动学参数（平均振幅、波速和频率）相同下，振幅逐渐增加的水下仿生鱼鳍推进器的游动速度要大于振幅保持不变的仿生鱼鳍推进器。

以上的实验研究以及数值分析大都主要集中在对于仿蝠鲼潜水器展长方向变形对流场特性的影响，很少涉及到剖面形状以及剖面运动状态对潜水器流场特性的影响。而剖面形状以及剖面运动状态对于整个潜水器的流场特性起着至关重要的作用。所以，本文针对剖面轮廓在大变形情况下对其流场特性的影响展开研究。

1 柔性翼计算模型的建立

根据图1的蝠鲼外所示的蝠鲼外形。在图2的轮廓中选取不同截面，如图 3应用不同的 NACA翼型作为该截面的翼型参数。通过三维软件得到如图4的仿生蝠鲼三维物理模型。

图1 仿生蝠鲼外形轮廓Fig. 1 Outline of bionic manta rays

图2 拟合后蝠鲼轮廓图/mmFig. 2 Outline of manta rays after fitting

选取其中 NACA0012翼型作为初始翼型建立柔性翼的运动模型，如图5，选取参考文献[9]中鳕鱼的振幅曲线，幅值函数为公式（2），其幅值曲线如图6所示。

图3 NACA翼型剖面Fig. 3 NACA airfoil section

图4 仿生蝠鲼三维模型Fig. 4 3D model of bionic manta rays

图5 NACA 0012的翼型曲线/mmFig. 5 Airfoil curve of NACA 0012/mm

由于通过观察发现蝠鲼的翼在摆动时与三角函数的波形近似，因此在幅值函数基础上引入正弦函数参量。考虑到变形是在翼型的基础上进行的，因此确定最终的运动方程为

式中：y0为在初始状态时翼型的初始纵坐标值；λ为波长（1、1.2、1.4、1.6、1.8）；f为频率（1、1.5、2、2.5、3、3.5、4）；t为时间（0～15 s）。

图6 柔性翼振幅幅值曲线/mmFig. 6 Amplitude curve of flexible wing/mm

根据得到的翼型波动方程，考虑到本文中翼型的变形比较大，运动规律也较复杂，利用动网格技术，采用弹性光顺与局部重构相结合的方案实现计算区域网格的更新。如图7所示为一个周期内不同时刻翼型运动与网络更新情况。

左侧边界条件设置为速度入口，速度大小为1 m/s，方向沿x轴正方向。右侧边界条件设定为自由出口，上下边界设定为滑移壁面。设定翼型壁面为无滑移壁面。采用有限体积法的离散方式，基于压力的一阶瞬态求解器，标准k-ε模型。设定时间步长为 0.001 s，待流场稳定后，对计算结果进行分析。

图7 一个周期内不同时刻翼型运动情况Fig. 7 Motion of airfoil with different time in one cycle

2 仿蝠鲼潜水器二维翼型流场性能计算结果分析

如图 8所示，柔性翼的运动在尾部产生了反卡门涡街。涡由柔性翼的前缘产生，随着柔性翼的运动沿着柔性翼表面向后移动，并最终随着柔性翼的周期性摆动从柔性翼的尾端脱落，涡脱离后的旋转方向与涡脱离柔性翼时翼梢的运动方向有关。

图8 柔性翼尾部涡量图Fig. 8 Vorticity diagram of flexible wing tail

2.1 频率对流场特性的影响

摆动频率是柔性翼运动时的一个关键参数，对运动过程中的流体动力特性的影响非常显著，取波长为1的情况为例，研究当摆动频率不同时，柔性翼的升力系数，阻力系数的变化规律。

如图 9所示，升力系数随着时间呈现周期性变化，变化周期与柔性翼运动的周期相同，通过对比不同频率的升力系数的变化范围可以发现，当波长一定时，摆动频率越大，柔性翼的升力系数幅值变化范围越大。升力系数的波形呈正弦形式，这是因为采用计算的翼型为NACA0012翼型，上下对称，在不考虑重力条件下，可以使得柔性翼在运动时的升力系数能够保证在一个稳定的区间内均匀波动。

从图10可知阻力系数的变化也是具有周期性的，阻力系数的变化周期与柔性翼摆动的周期是对应的，当波长一定时，阻力系数的变化幅值随摆动频率的增大而增大。图中的阻力系数小于0，表明柔性翼在摆动时所受到的力是延x轴负方向，即推力。这是因为柔性翼在运动时产生了反卡门涡街，使得柔性翼后方的流场速度加快，为柔性翼提供了向前的推力。

图9 波长为1时，升力系数随频率的变化规律Fig. 9 Variation rule of lift coefficient with frequency at 1 wavelength

图10 波长为1时，阻力系数随频率的变化规律Fig. 10 Variation rule of resistance coefficient with frequency at 1 wavelength

与升力系数曲线不同的是，阻力系数曲线在一个周期内有2个波峰，这是因为柔性翼在一个周期中会产生正反2个涡，所以可以得出柔性翼在一个周期内可以达到2次推力峰值。所以柔性翼的阻力系数在一个周期内达到2次波峰。

图11 波长L=1.8时，涡量云图随频率的变化情况Fig. 11 Variation of vorticity cloud pattern with frequency at 1.8 wavelength

以波长L=1.8为例，取15 s时的涡量图分析不同频率时流场的涡的形状。如图11所示，随着频率的增大，翼型波动产生的涡在离开翼梢后开始有向下偏转的趋势，且频率越大偏转的越明显。表明柔性翼在摆动时会产生沿y轴正向的升力。同时也可以较好地解释在图10中一个周期内阻力系数的两个峰值为什么会存在差异的问题。

2.2 波长对柔性翼流体动力特性的影响

在控制翼型运动的方程中，另外一个会对翼型的运动产生影响的就是波长，以频率等于1为例，研究其在相同频率下的升力系数，阻力系数随波长的变化关系。

由图12可知，升力系数随时间呈周期性变化，变化周期与翼型运动周期相同。频率一定时，随着波长的增加，升力系数变化的幅值也随之增加，这表明翼型收到的瞬时升力会随着波长的增大而增大。

由图13可知，阻力系数曲线的幅值随着波长的增大也在逐渐增大，在不同的波长情况下，柔性翼的运动存在相位差，这是由于在翼型的运动方程中存在x/λ项，致使波长不同时初相位也不同。柔性翼受到周期变化的推力，阻力系数的变化周期是柔性翼运动周期的2倍。

图12 频率为1时，升力系数随波长的变化规律Fig. 12 Variation rule of lift coefficient with wavelength at 1 frequency

图13 频率为1时，阻力系数随波长的变化规律Fig. 13 Variation rule of resistance coefficient with wavelength at 1 frequency

以频率为4 Hz为例，观察不同波长的涡量图在同一时刻的分布情况，如图 14。从图中不难看出，随着波长的增加涡的可见数量在减小，表明在波长较大的情况下，尾部涡量场的强度在降低，尾部涡更容易扩散消失发展成尾部湍流。

图14 f=4时涡量云图随波长的变化情况Fig. 14 Variation of vorticity cloud pattern with wavelength at frequency 4

3 结束语

通过分别对柔性翼的波动频率与波长进行分析，发现不同频率和波长与柔性翼的升力系数与阻力系数具有相关性。本文总结如下：1）在研究的频率范围内，波长一定时，升力系数的波动范围随频率的增大而增大，阻力系数的绝对值随频率的增大而增大。2）在研究的波长范围内，频率一定时，升力系数的波动范围随波长的增大而增大，阻力系数的绝对值随波长的增大而增大。