基于CEL法的柔性仿生机器鱼模型及巡游速度研究

张开升 邹庆彪 赵波 张保成

(中国海洋大学工程学院，山东青岛266100)

仿生机器鱼作为水下机器人的一种，具有高效、高速和灵活的特点，逐渐成为探索水下世界的重要手段。与传统的多刚体仿生机器鱼相比，柔性仿生机器鱼由柔性材料制造而成，具有更好的柔顺性、鲁棒性、高效性和低扰动性，能够更好地模拟鱼类的游动姿态和适应更加复杂多变的流场环境，逐渐成为国内外研究的热点。但当前国内外对柔性仿生机器鱼的研究还主要集中于采用实验手段进行柔性仿生机器鱼的功能实现研究，如Alvarado［1］制作了基于硅胶材料的柔性机器鱼，Marchese等［2］制作了由流体系统驱动的柔性机器鱼，Wang等［3］设计制作了四关节机器鱼，并致力于实现其最大游泳速度和更高的推进效率，而在鱼体设计参数对鱼体游动性能的影响规律方面研究还存在不足，因而也导致仿生机器鱼样机与真实鱼类的游动性能存在较大差异［4］。另外，通过实验的方法对柔性机器鱼进行研究存在研究周期长、研究成本高的缺点，且在样机制作完成后其设计参数不能改变，从而导致无法开展较为全面的研究，致使研究效率较低。因此，通过数值仿真的方法对柔性仿生机器鱼进行研究是一种可行且高效的研究手段。目前，许多学者对仿生机器鱼进行数值仿真研究，如Xia等［5］基于CFD(计算流体动力学)，采用UDF(用户自定义函数)和动网格技术对仿生机器鱼自主推进直线运动进行研究;Kern等［6］建立了鳗鱼的三维自主游动模型，并研究了鱼体周围受力的分布情况和尾迹涡形成与脱落的具体过程;杨亮［7］采用CFD的方法模拟仿金枪鱼尾鳍摆动时产生的尾涡特点，并分析其与游动性能之间的关系;胡亚南［8］以浸入边界法为理论基础，通过自主开发的算法程序研究仿生鱼的游动机理;刘鹏等［9］采用CFD的方法研究了柔性翼推力及效率;李宁宇等［10］基于浸入边界有限体积法研究拍动翼运动学、涡动力学和力的产生之间的关系;张开升等［11］采用CFD的方法研究了仿生机器鱼的多机体协同推进效率。但上述研究并没有考虑因为流固耦合作用而引起的鱼体变形，因此造成仿真结果存在较大误差。

1964年，Noh［12］提出的耦合欧拉 -拉格朗日(简称CEL)法能分析物体运动导致的流体大变形问题，在解决流体与结构之间的复杂相互作用问题时具有强大优势，如王永虎等［13］基于CEL法研究某型飞机撞入水面的流固耦合过程;王懿等［14］基于CEL法研究船舶抛锚入泥的流固耦合过程。因此，CEL法能够分析柔性仿生机器鱼在游动过程中鱼体的柔性变形以及鱼体与流体之间的流固耦合作用，能够为柔性仿生机器鱼的研究提供一种新的研究思路。本文采用CEL法建立柔性仿生机器鱼自主游动的流固耦合模型，研究鱼体刚度、摆动幅值、摆动频率对柔性仿生机器鱼巡游速度的影响规律，据此得到仿生机器鱼达到较优巡游速度的设计参数组合。

1 耦合欧拉-拉格朗日法

传统的分析方法，如拉格朗日法 (Lagrangian)，可以用来分析物体的受力及变形，但网格须与分析目标的变形保持一致，即材料与单元的变形必须保持一致。该方法可以简化求解，在小变形的情况下能够得到较为准确的节点及位移等信息，但变形较大时会导致网格畸形，从而无法收敛致被迫中止计算。而欧拉分析法 (Eulerian)中网格的位置属性、形状属性不会变化，包括流体在内的材料可在网格中流动，不会出现无法收敛的问题，因此该方法适于分析和处理变形较大的问题。但其网格的固定致使无法准确地捕捉物体的边界信息，且不能得到结构的准确受力情况。两种算法中的网格单元特征及连续介质变形情况如图1所示。

图1 网格特征及连续介质变形情况Fig.1 Mesh characteristics and continuum deformation

1.1 相互作用及控制方程

欧拉材料的质量守恒方程、动量守恒方程、能量守恒方程及连续性方程在欧拉坐标系下进行时间积分，且材料之间的相互作用使用黏性行为，拉伸应力可以在两个欧拉材料之间穿过界面进行传递，并且在这些界面上不会发生滑动。

CEL法的接触公式基于增强浸入边界法，为了模拟物体之间的相互作用，两套网格之间的接触在有限元软件Abaqus中使用基于罚函数方法的通用接触。仿真分析时，系统能自动地计算并且追踪拉格朗日结构和欧拉材料之间的界面，而不需要为欧拉域生成一个协调网格。该方法中接触面间的法向接触压力与两个接触面网格节点之间的法向间距h有关，如图2所示。

图2 CEL法的接触判定示意图Fig.2 Contact determination schematic diagram of CEL method

如果发生节点对物质界面的贯穿，界面力就会传递到欧拉材料的锚定点上，接触压力如下所示:

式中:Fp为接触压力;K为罚刚度系数，取决于拉格朗日和欧拉材料的特性;h为锚点与网格种子之间的距离。

1.2 显式动力分析

在有限元软件Abaqus中使用CEL法进行仿真计算，采用的是Abaqus/Explicit显式分析求解器。该求解器采用显式中心差分时间积分法则［12］，可用于分析柔性鱼在流场中的摆动行为。该分析方法的稳定性限制了Abaqus/Explicit求解器所能采用的最大时间步长，因此为了保证求解的精度和可靠性，该求解器默认使用自动时间增量，使时间增量尽可能地接近而且又不超过稳定性限制。时间增量与整个模型的固有频率相关。

1.3 材料属性的定义方式

流体的材料属性通过状态方程 (简称EOS)中的Us-Up状态方程，来描述模拟流体动力学材料或几乎不可压缩流体和描述材料的冲击力、体积和能量特性［15］。在CEL法中，流体的材料被描述为黏性可压缩牛顿流体，在黏性的作用下流场会产生涡流，可用下式描述:

式中:σ为Cauchy应力张量;p为压强;η为剪切黏度;˙e为应变率。

2 基于CEL法的流固耦合有限元模型

在Alvarado［1］提出的鲔科和鲹科鱼类的外形方程基础上，利用三维软件建立仿生机器鱼的三维模型，其长度为0.5 m;整个鱼体为拉格朗日体。由于本文只考虑柔性仿生机器鱼的直线巡游运动，鱼头摆动的幅度及变形较小，且通常头部的刚度较大，因此，为简化模型，将仿生机器鱼有限元模型头部设置为刚体并约束其只能沿X方向运动;由于柔性尾鳍的运动模式通常是随着鱼体的运动而被动摆动，因此将尾鳍设置为刚体并使尾鳍绕其与后颈部的连接处做摆动和平动，如图3所示。

图3 仿生机器鱼三维模型示意图Fig.3 Schematic diagram of three-dimensional model of bionicfish

鱼体材料始终在流场压力和鱼体运动的作用下变形，尾鳍始终随着尾部被动摆动，且被动摆动的幅度与鱼体材料的变形和流场有关，最终仿生鱼在流固耦合的作用下完成最终的游动表现。

流体域必须选择为欧拉体。根据仿生机器鱼的尺寸和可能的游动范围，将流体域尺寸设置为:长2.50m、宽0.35m、高0.25m。

柔性鱼体的材料属性主要考虑密度ρ、弹性模量E。根据文献 ［16］，流体域中水被定义为近似不可压缩的黏性流体，材料的基本属性参数如下:线性冲击波波速Us的系数c0=1485m/s，粒子速度Up的系数s=0，材料常数Γ0=0，密度ρ=1000kg/m3，黏度 μ =1 ×10-3Pa·s。

为尽可能地避免鱼体材料由于弹性模量较小而在流体中发生拉伸变形，在仿生机器鱼沿鱼体方向的中间截面上添加一层具有较大弹性模量的膜单元层，如图4所示。

图4 仿生机器鱼模型中的薄膜层Fig.4 Membrane layer in biomimetic fish model

使用波动方程来描述鱼类的波状运动，因此在有限元软件Abaqus中，选择仿生机器鱼体干上的若干点，根据波动方程获得每个点的侧向位移数据并输入软件中。波动方程建立的基础是仿生鱼做直线连续运动，尾部做被动摆动的运动，该运动模型有如下假设:①横向行进波在鱼体每一个截面上都是一维的，沿着头部向尾部传播，且没有沿行进方向的扭动;②在波动游动中，通过身体伸缩来实现这种运动［17］。根据波动方程可知，鱼身中部的摆动幅度远小于尾部，因此在尾部选取更多的点，每个点之间的运动幅度通过样条差值获得。

根据笔者前期的研究［11］，基于牛顿欧拉法和遗传算法获得了本文有限元模型中使用的基本波动方程:

式中，h(x，t)为从鱼头到鱼尾方向体干上坐标为x处t时刻的侧向位移。

由CEL法可知，拉格朗日体的表面自动和欧拉体进行接触，通过对拉格朗日体表面进行积分得到节点力 (欧拉体施加的压力)。拉格朗日单元相当于给欧拉单元施加了流动约束，部分拉格朗日体通过覆盖欧拉网格的体积和面，将流动约束施加给欧拉物质，因此欧拉单元的尺寸和拉格朗日单元的尺寸之间的比例对耦合效果具有一定影响，即欧拉单元的网格质量对计算结果存在影响。除此之外，罚刚度系数也会对计算结果产生影响。根据现有的硬件条件，网格收敛性分析对比组的设置及结果如表1所示。由表1可知，网格组5、组4的计算结果趋向收敛状态。为减少计算成本以提高效率，本文采用组2的方案，因此也导致最终计算的结果与收敛值相比略有偏差。

表1 网格收敛性分析结果对比Table 1 Comparisons of mesh convergence analysis results

有限元模型的网格数量及网格类型见表2。

表2 鱼体及流体域网格数量及类型Table 2 Number and type of mesh of fish body and fluid domain

柔性仿生机器鱼自主游动的流固耦合有限元模型如图5所示。

图5 柔性仿生机器鱼自主游动的流固耦合有限元模型Fig.5 Fluid-Structure-Interaction finite element model for selfpropulsion flexible bionic fish

3 设计参数对仿生机器鱼巡游速度的影响

3.1 设计参数组合的确定

本文采用正交试验来确定设计参数组合。正交试验能够以较少的试验成本得到最优的试验方案和结果，在式 (3)的基础上，本文的正交试验主要考虑鱼体刚度、摆动幅值和摆动频率这3个参数。通过密度、弹性模量和泊松比定义鱼体刚度，弹性模量定为10～50MPa。由鲹科鱼类的摆动幅值而将摆动幅值定为0.10l(l为体长，0.5 m)至0.20l。摆动频率定为1.72～5.16 Hz。将每个参数分为5个水平 (见表3)。

表3 因素水平表Table 3 Factor-level table

对三参数五水平选用正交表L25(56)确定不同参数组合和试验安排，共25组试验。

3.2 各个设计参数对巡游速度的影响

由正交表，利用有限元软件Abaqus进行仿真计算，使用极差法对仿真结果进行统计，并得到各设计参数对巡游速度的影响规律。

3.2.1 鱼体刚度对巡游速度的影响

将正交实验的结果拟合可得到各个设计参数对游动性能的影响规律。鱼体刚度对巡游速度的影响规律如图6所示。由图6可知，巡游速度呈现先增后减，然后趋于平稳的趋势，最大值出现在30MPa附近，在40至50MPa之间稳态巡游速度呈现略微上升趋势。

图6 仿真样本内鱼体刚度对巡游速度的影响Fig.6 Effect of fish stiffness on speed of cruise in simulated samples

3.2.2 鱼体摆动幅值对巡游速度的影响

鱼体摆动幅值对巡游速度的影响规律如图7所示。由图可知，巡游速度总体呈现先增后减的趋势，在0.13l附近具有较高的巡游速度。

图7 仿真样本内摆动幅值对巡游速度的影响Fig.7 Effect of amplitude on speed of cruise in simulated samples

3.2.3 鱼体摆动频率对巡游速度的影响

鱼体摆动频率对巡游速度的影响规律如图8所示。由图8可知，巡游速度总体呈现随摆动频率单调递增的趋势，在较高的摆动频率下可以获得较大的推进速度。

图8 仿真样本内摆动频率对巡游速度的影响Fig.8 Effect of swing frequency on speed of cruise in simulated samples

3.2.4 较优设计参数组合

由图6-8可知，当鱼体刚度为30 MPa、鱼体摆动幅值为0.13l、鱼体摆动频率为5.16 Hz时，柔性仿生机器鱼能获得较优的巡游速度，且通过仿真结果可知该设计参数组合下仿生机器鱼的巡游速度能够达到约1.4l/s。

3.3 柔性鱼体的变形

柔性鱼体发生弯曲导致尾鳍的运动与鱼体存在一定的滞后，从而影响两者之间的协同击水。对仿真进行观测，发现在鱼体刚度较小时，鱼体后颈部出现压缩弯曲现象，导致尾鳍运动的明显滞后，且在摆动换向时尾鳍未完成当前方向的摆动即被鱼体拖动向另一个方向运动，随机选取一段游动过程如图9所示 (t0为游动过程的某个时刻，T为摆动周期)。尾鳍与鱼体不能良好地协同击水，从而导致了性能指标较差。

图9 鱼体刚度较小时尾部出现的弯曲现象Fig.9 Bending phenomena in the tail with small stiffness

后颈部的弯曲程度会随着鱼体刚度的增加而逐渐减小，尾鳍运动滞后的程度也随之逐渐减小。鱼体刚度在30 MPa附近时，尾鳍与鱼体能较好地协同运动，因而仿生机器鱼能获得较高的巡游速度。

为更清楚说明不同刚度对鱼体弯曲的影响，给出不同鱼体刚度在摆动幅值及摆动频率分别为0.20l及5.16Hz时，相同时刻后颈部的弯曲程度，如图10所示，其中左图均为尾鳍摆动方向即将改变前的情况，右图均为摆动方向改变后的情况。

图10 尾部的弯曲情况Fig.10 Bending situation of fish tail

鱼体摆动幅值和摆动频率也会影响鱼体的弯曲情况。对样本内尾鳍摆动角度的最大值进行统计(见图11)可以发现:尾鳍最大摆动角度随着摆动幅值的增加呈现上升趋势，摆动幅值对尾鳍的最大摆动角度影响较为明显，在摆动幅值为0.20l时，尾鳍最大摆动角度达到最大值，大于1.570 8 rad(90°)，造成尾鳍的位置在前进方向上超过尾柄，因此形成了明显的负击水行为，从而导致游动性能的下降;摆动频率对尾鳍最大摆动角度的影响较小，随着摆动频率的增加，尾鳍最大摆动角度呈现下降的趋势，在摆动频率为1.72 Hz附近时达到最大值 (约1.38rad)，而在5.16Hz附近时达到最小值 (约1.32rad)。

图11 仿真样本内摆动幅值和摆动频率对尾鳍最大摆动角度的影响Fig.11 Effects of amplitude and frequency on maximum swing angle of caudal fin in simulated sample

结合图7、图8、图11可以看出，在最大摆动角度为1.32rad附近时，柔性仿生机器鱼能达到较高的游动速度。

3.4 流场变化情况

选取新月形尾鳍中部的一部分，对经过该部分边缘的流线进行跟踪，如图12所示。由图可知:在摆动过程中，尾鳍的摆动会使附近流体加速，从而逐渐形成闭合的环线，即逐渐形成涡，其旋转速度和尾鳍摆动的角度和速度有关。

图12 尾鳍中部形成的流线示意图Fig.12 Schematic diagram of streamlines formed in the middle of caudal fin

仿生鱼游动过程中，流体材料的某一瞬时压力情况如图13所示。由图可知，尾鳍的摆动使流场产生高压区和低压区，尾柄处的流体压力明显较大，说明鱼体通过摆动从流场获得了较大的反作用力;尾鳍处的涡团逐渐脱落并散落在流场中，与另一侧脱落的涡团相对应;鱼体表面的流体压力比尾鳍处的压力小，前后的压力差驱动仿生鱼前进。

图13 仿生鱼游动过程的某一时刻Fig.13 A moment in the swimming process of bionic fish

4 结论

本文基于CEL法研究了柔性仿生鱼的自主游动过程，首先建立了柔性仿生机器鱼自主游动的流固耦合模型，然后在此基础上，通过在有限元软件Abaqus上进行仿真计算得到鱼体设计参数 (鱼体刚度、鱼体摆动频率、鱼体摆动幅值)对仿生机器鱼巡游速度的影响规律，最后据此得到较优的设计参数组合。在样本范围内，当柔性仿生机器鱼具有一定的鱼体刚度时能获得较好的游动性能，对于本文的模型，当鱼体刚度约为30 MPa时具有较好的游动性能;鱼体刚度能影响鱼体的弯曲程度，在此基础上，鱼体摆动幅值和摆动频率能够进一步影响鱼体的弯曲程度和尾鳍的最大摆动角度，从而影响游动性能，对于本文的模型，鱼体摆动幅值约为0.13l、摆动频率约为5.16Hz时具有较好的游动性能。因此，对于本文的模型，基于CEL法获得了较优的设计参数组合，当鱼体刚度约为30MPa、鱼体摆动幅值约为0.13l、摆动频率约为5.16 Hz时具有较好的游动性能，巡游速度达到约1.4l/s。