机器鱼多鳍协同推进加速——滑行行为数值模拟

李浩宇,李宗刚,杜亚江,王小波

(1.兰州交通大学 机电工程学院, 兰州 730070; 2.兰州交通大学 机器人研究所， 兰州 730070)

1 引言

仿生水下机器人是一种运用真实鱼类的游动机理在水下进行快速游动、快速转弯等平面或空间运动的机器人系统，因其机动性好、推进性能强以及隐身性好的优势在水下勘探、生物监测、军事侦察与打击等方面具有重要的应用价值，是近年来国内外研究的热点。

海洋中鱼类的游动可分为身体和/或尾鳍(body and/or caudal fin,BCF)游动模式以及中央鳍和/或对鳍(median and/or paired fin,MPF)游动模式[1-3]，国内外科研人员通过深入研究鱼类的游动机理，研发了各类独具特色且拥有不同游动性能的仿生机器鱼[4]。为分析解决机器鱼运动与周围流场结构的相互作用问题，现阶段大多数研究利用数值模拟方法对流场进行处理以得到涡结构、速度矢量等物理量，进而阐明机器鱼在流场中的游动机理，优化机器鱼的行为。针对流场中胸鳍或尾鳍对流场的作用问题，李宁宇[5]将单胸鳍放入流场，分析研究胸鳍摆动时不同参数对流场结构的影响，发现单个胸鳍摆动时胸鳍的推力和效率达到最大时相位差应保持在90°；Zhang等[6]、Krishnadas等[7]等通过数值方法模拟了尾鳍形状对推进性能的影响，结果表明，在相同运动规律下，仿金枪鱼形状的尾鳍在摆动时的推进性能和效率最佳；冯亿坤[8]研究了机器鱼在C型启动过程中的水动力性能，分析了c型启动2个阶段产生的涡环和射流；Xu等[9]运用重叠网格法，将刚性胸/尾鳍放入流场中，仅靠胸鳍的摆动来推动鱼体直游和转弯；刘焕兴等[10]、Wu等[11]对仿生机器鱼的尾鳍柔性变形自主游动进行了研究。以上研究在分析单个胸鳍摆动或尾鳍波动的水动力性能中取得了一些成果，但并未考虑胸尾鳍同时运动时的协同作用。

胸尾鳍协同游动即在上述文献的基础上，将胸鳍的摆动和尾鳍的波动相结合，通过两者的共同作用使鱼体做出相应动作，它相较于单纯的胸鳍或尾鳍推动，在推进性能与效率上有明显优势，通过胸鳍涡与尾鳍涡的融合，机器鱼拥有更好的推进性能与机动性能；在胸尾鳍协同游动的研究中，Li等[12]分析了河豚的柔性多鳍动力学，将含柔性背鳍、臀鳍、尾鳍的仿河豚模型放入流场进行推进，分析了其水动力学性能和涡流结构。在多体动力学研究方面，Kanso等[13]分析了离散的鳗鱼模型在理想流体中的自推进游动；John等[14]研究了连续的鳗鱼模型在自主推进时的游动速度与尾涡的变化。

在其他仿生机器鱼相关的研究工作中，谢鸥等[15]针对针对仿生机器鱼作业时的环境识别问题，提出了一种基于人工侧线(ALL)的近壁面流场识别方法，采用CFD方法提取数据并运用多层前馈神经网络建立预测模型，通过验证，该方法具有良好的预测效果；Zheng等[16]建立了由压力传感器阵列组成的人工侧线系统来感受机器鱼游动时周围流场的变化，该研究验证了人工侧线系统的有效性和实用性；刘科显等[17]通过理论分析、CFD仿真以及估计的方法研究了二维多关节机器鱼远场流速与攻角预测；和岩辉等[18]基于中枢模式发生器(CPG)理论结合模糊控制器提出了一种可实现机器鱼定向游动的精确控制方法，通过仿真与实体实验验证了该算法的稳定性；梁旭[19]研制了一种高频驱动在线便变刚度柔性仿生机器鱼，并进行了该机器鱼的游动特性实验，验证了鱼类依靠改变刚度从而改变固有频率来匹配摆动频率的可行性；Yu等[20]研制了一种无线电控制的多关节机器鱼，且通过模糊强化学习算法实现了多个机器鱼之间的协调，并成功应用到国际水中机器人比赛中；此外，Yu等[21]还研究了多连杆机器鱼的转向控制问题，开发了一种使用方法控制机器鱼的转弯步态，并通过仿真和实体实验验证了该方法的适用性。

基于流体力学(CFD)方法，数值模拟了仿生机器鱼在水平面内的胸尾鳍协同自主推进直线游动，其中胸鳍采用3自由度刚性摆动，尾鳍采用柔性波动；胸鳍的摆动采用重新设计的3自由度摆动曲线，并设置了4组不同占空比，分析占空比对机器鱼游动速度以及水动力性能的影响。鱼体由静止开始，施加特定的运动规律进行直游，在经过一定时间的加速后进入稳定巡游阶段。对整个过程中鱼体与胸鳍的速度矢量、压力云图以及流场的三维涡结构进行分析，揭示鱼类直线游动的机理，为之后的胸尾鳍协同推进研究提供参考。

2 机器鱼的计算模型

2.1 胸鳍的计算模型

仿生机器鱼的三维几何模型与相关尺寸如图1所示，其中胸鳍与鱼体采用分离式结构，胸鳍进行3自由度摆动，两侧胸鳍初始位置贴近鱼体以模仿真实鱼类。该模型的运动部分有3个，分别是两侧的刚性胸鳍以及柔性尾鳍，2个胸鳍通过将上下拍动、前后拍动以及摇翼运动3个自由度运动耦合，以实现期望的运动轨迹，其中机器鱼整体长度为1 000 mm，宽度为205 mm，尾鳍高度为370 mm，胸鳍展长为170 mm，最大弦长为130 mm。鱼体分为PA和PB段，其中PA段为鱼头刚性部分，长度为400 mm，PB段为尾鳍柔性部分，长度为600 mm。胸鳍的运动学模型可由下式给出：

(1)

式中：φR、φFL、φF分别代表胸鳍的前后拍、上下拍、摇翼运动的欧拉角;φRC、φFC、φFLC分别代表前后拍、上下拍、摇翼运动角的平均值;φR0、φF0、φFL0分别代表3个运动的幅值; dφF、dφFL代表相位差;ω为运动角速度;t为时间。

图1 仿生机器鱼的三维模型及相关尺寸

图2为胸鳍侧视图中随体坐标系坐标轴与胸鳍外缘的交点随时间变化的运动轨迹，其中每个椭圆表示特定幅值下胸鳍前后拍翼角与上下拍翼角的相对关系，所有椭圆的交点即胸鳍在t=0时刻的位置，该曲线为胸鳍前后拍与上下拍2个自由度的约束条件，此外，由于该图为胸鳍XOZ平面的运动曲线，摇翼运动与另外2自由度的运动相对关系并不包含在该胸鳍摆动曲线内。

图2 胸鳍的摆动曲线

“加速-滑行”游动模式是鱼类等水生生物常见的游动方式，自然界中的鱼类在游动时胸鳍的摆动总呈现出一定的间歇运动，即胸鳍除了向后拍动，借助水的作用力推动鱼体游动的加速阶段之外，还包括回摆的过程以及2个摆动周期之间胸鳍相对静止的过程，由于胸鳍并未提供向前的推力，鱼体本身靠惯性作用向前滑行，因此以上2个阶段为机器鱼的滑行阶段。在“加速—滑行”行为的研究中，万宏[22]分析了鱼尾鳍波动时的“加速1滑行”游动模式，发现在2个运动周期之间加入空白的滑行阶段可提高机器鱼游动时的效率，降低能耗比。将该模式代入胸鳍的运动中，可防止胸鳍摆动时动力阶段产生的涡过早地被恢复阶段胸鳍的回摆破坏掉，充分利用胸鳍后摆所产生的推力。

占空比，即在胸鳍的完整摆动周期中摆动阶段需要的时间所占用的比例，其表达式如下：

(2)

式中：Tb为一个周期内胸鳍摆动阶段的时间;Tc为滑行阶段的时间。

为方便直观的显示占空比，现取2个周期内前后拍翼角的变化来表示胸鳍的运动，将占空比为1∶1、0.8∶1、0.6∶1、0.4∶1 时胸鳍的运动表示出来，如图3所示，在周期一定的情况下通过改变胸鳍摆动的时间来调整占空比，在相同的时间内，不同占空比拥有相同周期。

2.2 鱼体的计算模型

由于尾鳍采用了柔性波动，故运动规律的表示与胸鳍不同，设鱼体的顶点处为惯性坐标系的原点，鱼头向鱼尾的方向为X轴，垂直鱼体水平方向为Y轴，垂直鱼体竖直方向为Z轴，鱼体的运动模型定义如下：

图3 不同占空比下前后拍翼角的变化

(3)

式中：y为鱼体在X位置的横向位移，鱼体的中轴线在一个周期内的波动曲线如图4所示，由于鱼体采用柔性波动，在UDF中运用节点运动将每个时间步内鱼体中每一个节点的位移表示出来，并输出到鱼体模型上。

2.3 流场中的其他参数

机器鱼在游动过程中的推力系数与侧向力系数可分别由下式表达：

(4)

式(4)中：ρ为流体密度;Fx、Fy分别为鱼体受到的推力与侧向力;S为鱼体的投影面积。

鱼体的推进效率为平均输出功率与平均输入功率的比值，平均输出功率为鱼体运动速度与游动过程中平均推力的乘积，平均输入功率为

(5)

式(5)中:σ为鱼体的微元受到的力；V为该微元的速度，该乘积在鱼体的表面积分。由此可得鱼体的平均推进效率为

(6)

图4 尾鳍波动曲线

3 流场的控制方程与数值方法

该流场的控制方程为非定常、连续、不可压缩的三维Navier-Stokes方程，其表达式如下：

(7)

其中数值计算的网格如图5所示，由于模型需要自推进来向前游动，设机器鱼体长为L，计算域的尺寸为8L×1.5L×0.5L，且X方向的网格均匀分布，以保证在鱼体游动过程中计算保持相同精度；流场网格划分形式为四面体网格，模型的壁面采用三角形网格。

图5 计算域的网格划分

计算区域的边界条件为：入口边界与出口边界的速度和压力梯度均为0。由于在近壁面的流场中需要精度较高的结果，因此，本文实验中所采取的湍流模型为k-ωSST模型，且采用低Re数校准，采用SIMPLEC算法对连续方程中的压力和速度进行耦合，其中压力、动量、湍流动能均采用二阶迎风格式，扩散项采用二阶中心差分格式，瞬态方程采用二阶隐式。监测鱼体在3个方向上的水动力系数与X方向上的速度与位移。

本文采用有限体积法(finite volume method,FVM),并利用动网格技术，通过刚体运动与节点运动来实现胸鳍的摆动和尾鳍的波动，并使用弹簧近似网格光顺和局部网格重构方法以防止网格变形太大导致网格负体积的现象出现。该技术通过一套网格便实现了鱼体、胸鳍运动的同时整体向前推进，不失精度的同时计算量相对较小。

鱼体在计算域中采用牛顿定律以向前游动，鱼体在t时刻前进的速度可由下式表示：

(8)

式中:Ab和Af分别为鱼体与胸鳍的表面积;Fx(t)为机器鱼在t时刻受到的X方向的合力，包括鱼体收到的X方向的静压与左右胸鳍受到的力;m为机器鱼的整体质量。

4 数值计算及结果分析

4.1 网格收敛性的验证

在数值模拟中，为验证网格的收敛性，取3种网格尺寸，对每一种网格进行同一条件的数值计算，最终计算出每一种网格的水动力系数如图6所示，其中网格1、2、3分别对应尺寸粗、中、细，通过图6可知，3种网格下的计算结果差异较小，因此，网格设置与数值计算方法合理且收敛，以下仿真结果均通过第二套中尺寸网格计算得到。

图6 网格收敛性验证

4.2 占空比对推力系数与速度的影响

通过对不同占空比的仿真计算，达到稳定巡游状态下其中一个周期的机器鱼推力系数与升力系数、其中一侧胸鳍推力系数与所有周期的机器鱼推进速度如图7所示。

由图7(a)可知，单个胸鳍的最大推力系数幅值随占空比的减小而增大，占空比的改变导致了胸鳍摆动的频率的改变，频率越大，占空比越小，胸鳍推进效果越好，但是功率消耗也越大；由图7(b)可知，机器鱼的升力系数幅值同样随占空比减小而增大，且作用时间逐渐减少，在胸鳍滑行阶段，机器鱼的升力基本为0，可见机器鱼的升力主要是胸鳍的摆动产生的；由图7(c)可知，随着占空比的减小，机器鱼的最大推力系数逐渐增大，且在所有占空比下，推力系数正方向的数值略小于负方向，因此，在整体上机器鱼的推力略大于阻力，机器鱼的推力作用时间随占空比的减小而减小，到滑行阶段时，推力系数波动较小；由图7(d)可知，随着占空比的减小，稳定巡游速度逐渐增大，但是速度波动也逐渐增大，机器鱼达到稳定巡游状态的时间不随占空比的改变而改变。

图8为不同占空比下机器鱼整个游动过程的平均推力、平均升力与升阻比。由图8(a)可知，随着占空比的增大，平均升力系数逐渐增大，且增大的幅度逐渐减小，占空比为1∶1即没有中间滑行阶段时平均升力系数最大，占空比0.8∶1时机器鱼的平均升阻比最大；由图8(b)可知，平均推力系数随着占空比的增大而增大，结合图7，占空比越大，推力系数峰值越小，但是平均推力系数越大。

图7 不同占空比下推力系数和速度随时间的变化

图8 不同占空比下的平均推力系数、 升力系数与升阻比

4.3 流场速度矢量的分析

在胸鳍的运动规律中，取φRC=φFC=φFLC=0，取φR0=π/4、φF0=π/6、φFL0=π/12，取dφF=dφFL=π/6，尾鳍的波动方程中，取c0=0.02、c1=-0.08、c2=0.2，波数k=π、ω与胸鳍保持一致，在静水中进行自推进数值模拟，对胸鳍摆动占空比为1∶1时通过在静水中的数值模拟得到的速度矢量图进行分析。

取某个周期的6个典型时刻，每个时刻的速度矢量如图9所示。在t=T/6时，鱼体尾鳍快速向上弯曲，在尾部左侧产生了较大射流，从而产生一个较大的作用在鱼体上的侧向力，通过上个周期的运动，鱼体侧后方存在一个较大漩涡，该漩涡作用在鱼体上产生了一定推力。随后，尾鳍摆动幅度逐渐减小，尾鳍处的射流逐渐指向左后方，该射流的反作用力推动鱼体向前游动，鱼体受到的推力逐渐增大，漩涡脱落，开始远离鱼体；胸鳍开始展开，其周围出现侧向射流，到t=T/3时，胸鳍继续通过摆动产生较小漩涡。

t=T/2时，尾鳍向左弯曲产生射流，该射流为鱼体提供较大的侧向力，而在鱼尾后方，新的漩涡逐渐开始产生，增强鱼体推力；胸鳍继续向外侧展开，得益于摇翼运动，胸鳍两侧压差小，鱼体受到的阻力较小。

t=2T/3至5T/6，尾鳍伸直，并开始向右侧弯曲，侧向射流进一步增强，导致鱼体受到较大侧向力，后方漩涡成型并逐步远离鱼体，产生射流，推动鱼体游动；通过尾鳍作用新的漩涡产生。此时胸鳍向后划水，产生一定推力。

t=T时，尾鳍和身体逐渐回到初始位置，射流再次在产生在尾鳍左侧，后方的漩涡完全成型；胸鳍也回到初始位置，开始下一周期的运动，其后方存在的射流加强了推力。

图9 流场的速度矢量

以上为一个周期内鱼体周围流场的速度矢量分析，通过尾鳍和胸鳍的摆动周期性的产生射流与漩涡，鱼体间歇性的受到推力与侧向力，推力推动了鱼体向前游动，侧向力总体上相互抵消，鱼体不会发生明显横向运动。

4.4 流场压力云图的分析

图10为上述周期同时刻机器鱼周围流场的压力分布。t=T/6时鱼体两侧压差较大，且鱼体后方有前一周期从鱼体周围脱落的低压涡；t=T/3时，鱼体两侧压差进一步增大，且低压区域向后运动，为鱼体游动提供了推力，胸鳍两侧产生了不利于向前游动的压差；t=T/2时，低压区域进一步脱落，在鱼体右侧同样产生了低压区域；t=2T/3时，鱼体左侧的低压区域脱落到鱼体后方，鱼体两侧压差再次增大，为鱼体提供了较大推力；t=5T/6时，胸鳍的摆动使两侧压差变大，对鱼体的前进产生了积极影响，鱼体右侧低压区域开始脱落；t=T时，鱼体左侧再次产生低压区域。

4.5 鱼体周围流场的三维涡结构

取上述同一周期同时刻鱼体以及周围流场，生成其三维尾涡结构，如图11所示，三维涡结构的识别采用Q判据，涡结构的颜色反映了涡量的大小。

胸/尾鳍协同自推进机器鱼在游动过程产生了十分复杂的流场涡结构，下面就稳定巡游状态下一个周期的6个时刻做简要分析。胸鳍和尾鳍的拍动产生了大小、方向不同的垂直于XY平面的涡环，这些涡环导致了机器鱼的前游。

图10 流场的压力分布

t=T/6时，在尾鳍处，之前的运动产生了两列涡环，尾鳍附近的涡强度较高，该涡为该时刻鱼体收到的推力与侧向力的主要来源，上个周期脱落的涡环逐渐远离鱼体，对鱼体的作用逐渐减弱；在胸鳍的运动中，上个周期产生的涡环仍在产生作用，胸鳍的摆动正在产生新的涡环。由于尾鳍的摆动，在t=T/3时涡环逐渐增大，涡强度也逐渐增大，且主要集中在尾鳍的一侧，一定程度上增强了尾鳍的侧向力；在胸鳍处，上个周期产生的涡环完全脱落，新的涡环开始产生。

通过连续作用，t=T/2时，尾鳍处逐渐产生新的涡环VB1，该涡环逐渐增大并向机器鱼的左侧缓慢移动；两侧胸鳍处也产生了新的垂直于XY平面的涡环VF1与VF2，此时胸鳍与尾鳍产生的涡分别作用在机器鱼上，暂未发生涡的融合现象。随后VB1向后方移动，到t=2T/3时，尾鳍处的新涡VB2诞生，VB1的头部与VB2的尾部相连接，且VB1与VB2在空间上相互垂直，共同对鱼体产生作用；胸鳍处VF1与VF2同时开始从胸鳍前缘与已经脱落的涡环一样，慢慢向尾鳍方向移动，逐步靠近尾鳍处的涡。

至t=5T/6，VB2逐渐增大，胸鳍的第一个涡环逐渐与VB1融合，鱼体受到的推力被加强。随着VB2向侧后方移动，t=T时，右侧胸鳍的涡开始与VB2融合，VB1与VB2开始逐渐分离并脱离鱼体，并向后移动，尾鳍通过运动开始产生新的涡环，胸鳍处，随着胸鳍的摆动，VF1与VF2强度增大，新的周期开始。

图11 流场的三维涡结构

4.6 实体实验验证

为验证数值计算的准确性，通过实验室现有机器鱼进行了静水中的直游实验，并通过超声波传感器监测机器鱼的前进方向实时位移数据，实验过程共持续8 s，其中6个时刻机器鱼位置及运动轨迹如图12所示，机器鱼在流场中前进的同时，由于尾鳍的周期性波动，产生较小的横向位移，但总体上呈现直游行为。CFD仿真与实体实验中机器鱼前进位移随时间的变化如图13所示，相同时间机器鱼的实际位移均小于仿真结果，其主要原因为实验中机器鱼外壳的近亲水性材料、实验时的壁面效应以及外界因素的干扰增加了游动阻力；可在机器鱼外壳填充疏水性材料、增加流场尺寸以及优化实验流程来减阻，以改进实验方法。实验总体误差在期望范围内，且总体游动趋势一致。

图12 机器鱼游动过程

图13 机器鱼位移随时间的变化

5 结论

采用CFD方法数值模拟了不同占空比下仿生机器鱼在水平面内的自主推进运动，其中仿生机器鱼带有一对分离式的3自由度胸鳍，其中胸鳍采用刚性摆动，机器鱼体采用柔性波动。通过分析不同占空比机器鱼自主推进过程中的鱼体游动性能与机理、水动力性能、流场的速度矢量、压力云图与三维涡结构，得出了以下结论：

1) 在机器鱼的自主推进游动过程中，通过胸尾鳍的协同推进，鱼体胸鳍后方周期性的产生射流与漩涡，通过射流与漩涡的共同作用，为机器鱼提供了向前游动的推力。

2) 在机器鱼的整个摆动过程中，胸鳍和尾鳍两侧均有较大压差，高压区域脱落到鱼尾附近推动机器鱼向前游动。

3) 在每个周期的游动中，流场的涡结构都产生了复杂的变化，鱼尾的摆动产生了卡门涡街，并且胸鳍产生的涡向后脱落，逐渐与尾鳍的涡融合，增加了涡的强度。

4) 对比不同占空比下机器鱼整体和胸鳍产生的推力系数，发现随着占空比的减小，推力系数的振幅逐渐增大，但是持续时间相应的减少，到滑行阶段时，推力系数的波动较小。在速度的对比中，发现随着占空比的减小，机器鱼的稳定巡游速度增大，但是速度波动较大，占空比越小，鱼体的推进效果越好，但是占空比的减小也导致了功率消耗的增加。