水下仿生贝壳电磁驱动器动力学分析

崔 磊，黄兴保，杨斌堂

（上海交通大学机械系统与振动国家重点实验室，上海 200240）

随着工业的飞速发展和人口的不断增长，陆地上可利用的资源日趋匮乏，因此合理开发和利用海洋中各类资源对人类的发展来说是必不可缺的。无人水下航行器和水下机器人在海洋探索、环境监测、资源开发、摄影记录等众多方面起着不可替代的作用，而水下推进方式直接决定了各类水下运动系统的动力特性[1-2]。海洋生物种类繁多，经过长久的进化，大多都具有适应水下生存的外形构造和水中运动的方式，根据已知海洋生物主要运动方式所产生的推进力进行分类，当前的水下仿生推进机制主要包括：模仿鱼类和海豚鲸鱼类，借助身体和鱼鳍波动或摆动时产生的涡流推进，并分为奇鳍/对鳍（Medial and Paired Fin,MPF）推进、身体/尾鳍（Body and Caudal Fin，BCF）推进[3-5]；模仿扇贝、水母和章鱼等将流体吸入后再快速喷出从而实现喷射推进[6]。

在上述的推进方式中，鳍/体的推进方式得到了广泛的关注，虽然喷射推进已经被证明有效，但相关的研究仍较少，特别是关于仿生贝壳推进机构的研究。目前已有的模仿章鱼或水母喷射运动的机器人，例如使用脉冲喷射推进技术开发的软体机器人最大能够达到4 cm/s 的运动速度，然而头足类的柔软体导致其难以生活在复杂而不稳定的水下环境[7]。研究发现，双贝类贝壳独特的外形结构能使其更好地存活在复杂的水下环境中并减小其运动过程中受到的流体阻力，此外，该贝壳的喷射推进方式容易实现且最高能达到每秒7 体长的运动速度[8]。但目前很少有使用双贝类贝壳作为研究对象的仿生推机构，有学者利用磁力驱动和电机驱动借助仿生贝壳研究其在低雷诺数下的往复运动[9]，但其外形与自然界中的贝壳差别较大，还有学者分别利用往复曲柄机构及旋转电机研制出了仿生贝壳机器人，最高能够达到每秒4.65体长的瞬时运动速度[10-11]。仿生贝壳在运动过程中始终保持往复转动，但之前学者的研究都以设计和实验分析为主，并未建立其动力学模型并进行进一步分析。

受双贝类贝壳结构和运动方式的启发，本文提出一种基于电磁感应的仿生贝壳水下驱动器，根据自然界中双贝类贝壳结构设计仿生贝壳执行机构，以单侧执行机构为研究对象，建立其动力学模型，提取模型中的关键参数之后构建驱动磁力矩以及永磁体之间磁力矩的数学模型并验证其准确性，根据所设计的仿生贝壳和动力学参数优化目标来确定永磁体的阵列排布位置，分析优化后的执行机构动力学模型及其动力响应。本文主要内容包括仿生贝壳执行机构动力学分析、驱动磁力矩和永磁体间磁力矩建模计算、动力学参数优化等。

1 仿生贝壳驱动器结构

双贝类贝壳通过往复的张开、闭合运动实现喷射流体推进，若将仿生贝壳放置于空心圆管内，它可以通过往复运动从而推动圆管中的流体运动实现驱动，还可以通过保持张开并与圆管内壁贴合从而阻断管道中的水流。受其启发，本文提出一种以仿生贝壳为执行机构，以电磁线圈、永磁体实现非接触式直接驱动的水下驱动器，如图1所示。

图1 仿生贝壳驱动器结构

通过定位销将仿生贝壳约束在圆管中使其仅能绕轴转动，再将均匀磁化的永磁体对称布置在仿生贝壳内侧或表面，向放置在圆管外侧的电磁线圈通入电流产生空间磁场，则永磁体在磁场中受到驱动磁力矩的作用带动执行机构运动。相比于电机驱动，电磁非接触式直接驱动具有结构简单、响应快的优点，改变通入电流的幅值和频率可以实现不同的驱动效果。此外，可以通过将多个驱动器同向串联来提高整体的驱动能力。

2 仿生贝壳执行机构动力学模型

2.1 执行机构一般动力学模型

对于仿生贝壳驱动器，执行机构始终在一定范围内进行往复转动，以单侧执行机构为研究对象，执行机构包括永磁体和仿生贝壳，如图2所示，通过定位销将执行机构约束在圆管内令其仅可绕轴转动，设其转动惯量为J，执行机构受到的力矩包括两侧永磁体之间的磁力矩T→pm和永磁体在通电线圈磁场中受到的磁力矩T→c，这两项均与永磁体体的形状尺寸、放置位置及磁化强度有关，此外，T→c还和通电线圈尺寸、位置及输入电流相关，由于定位销轴线垂直于地面，重力可忽略不计。

图2 执行机构动力学模型

当仿生贝壳模型、永磁体参数及位置确定后，两侧永磁体之间的磁力矩T→pm仅与执行机构转动角度θ有关，可以将T→pm等效为系统中弹性力产生的力矩，设其与转动角度的比值为弹性系数K，执行机构运动过程还受到水阻尼的作用，设阻尼系数为C，永磁体在通电线圈磁场中受到的驱动磁力矩为外力矩，因此单侧仿生贝壳的往复运动可以等效为执行机构的单自由度振动，忽略系统中可能出现的碰撞变形，根据刚体动力学方程可以建立其一般动力学方程：

对于质量连续分布的刚体，由转动惯量定义和平行移轴定理可得J的计算公式为：

易知与外部输入电流无关，且永磁体间磁力矩通常与它们之间距离的平方成反比，因此该模型中的弹性系数K一般是随转动角度改变，此外，由于永磁体的尺寸相对于线圈的尺寸不可忽略即动力学响应分析的外力矩幅值A随着执行机构转角变化而变化，因此该一般动力学模型为变刚度的非线性时变系统，系统复杂且难以分析其动力响应，之后将确定并优化各系统参数J、K、A从而简化动力学模型便于进一步分析。

2.2 磁力矩模型

确定系统动力学模型参数首先需要建立和的数学模型，即推导计算系统静磁力矩模型。

磁场是有旋无源场，计算永磁体在磁场中以及永磁体之间的磁场力和力矩时，常用的方法有等效磁荷法、等效电流法和有限元法等。本文分别根据等效磁荷法和等效电流法建立永磁体和电磁线圈的磁场模型。

等效磁荷模型将永磁体简化为等效磁荷的分布[12-13]，对于均匀磁化的永磁体，则可以等效为仅上下极面有磁荷分布的模型，σm为磁荷的面密度。假定始终大于，要使两侧的仿生贝壳同时张开或闭合，则两侧永磁体应对称且同极相对布置，如图3所示，当电磁线圈中通入图示方向的电流时，根据永磁体的极面位置易知此时两侧执行机构均有闭合的运动趋势，反之亦然。

图3 驱动磁力矩示意图

推导和的显式表达式需要先建立通电线圈和永磁体的磁场模型，分析其在空间任意点处产生的磁感应强度。

设有限长通电线圈长度为L，线圈内外径分别为R、R+rc，单位长度的匝数为n，电流大小为I，在图示直角坐标系中，由通电线圈的对称性，仅需计算YOZ平面内的磁场分布，其中Y、Z分量分别为磁感应强度的径向和轴向分量。根据毕奥-萨伐尔定律：

可以得到多层通电线圈在空间中任意点N（x0,y0,z0）处的磁感应强度沿y、z方向的分量为：

式中：μ0为真空磁导率，R1、R2分别为通电线圈的内外径。通电线圈在任意点N处的磁感应强度为：

根据等效磁荷模型，永磁体极面的磁荷在空间中任意点产生的磁感应强度为：

对于两个均匀充磁且同性极面平行相对的长方体形永磁体，几何尺寸分别为2a×2b×2c和2a′×2b′×2c′，磁铁的磁化强度矢量分别为和，且方向相反，在其中一个永磁体1 的几何中心建立全局坐标系OXYZ，在磁铁B的几何中心建立局部坐标系O′X′Y′Z′，永磁体2的几何中心在全局坐标系中的坐标为（x,y,z），研究永磁体2受到的磁场力，首先：

永磁体在点N处产生的磁感应强度在全局坐标系中可表示为：

由于执行机构仅能绕定位销轴线转动，本文仅需计算和的x方向分量。

考虑到磁荷极面之间的空间存在相对磁导率为μr的导磁介质，当导磁介质被外界磁场磁化后，极面上磁荷的面密度也会随之改变，改变后的磁化强度与磁场的磁感应强度及永磁体的剩磁的关系式可得：

在有限长通电线圈内，转轴位置处为（0,0,z′），则转轴处到永磁体两个极面上的任意磁荷的矢量分别为，则永磁体2受到的磁力矩为：

执行机构受到的驱动磁力矩的解析表达式不易获得，为了验证之前建立的磁场矩模型，采用有限元分析（FEA）的方法，借助Ansoft Maxwell 进行有限元仿真验算。利用Ansoft Maxwell软件可以计算出不同线圈匝数、通入电流、布置位置、初始永磁体转角的情况下永磁体在磁场中受到的磁力矩。通过上文所推导的磁力计算公式，通过MATLAB计算相同参数条件下的磁力矩的数值解，此为公式预测值。将公式数值解与有限元仿真解进行对比，就可以判断所建立模型以及推导公式的正确性和准确性。

仿真所选用的永磁体的材料为钕铁硼NdFe35B，相关计算参数如表1 所示。为了验证磁力矩公式，以永磁铁2中心坐标y、z以及夹角α为变量，计算出永磁体不同位置与布置角度下的力矩值，通过对比理论推导得到的数学模型的数值计算结果和有限元法Ansoft Maxwell仿真结果来进行验证。

表1 模型验证计算参数

由图4 可以看出，数学模型计算与有限元仿真得到的磁力矩在不同位置及角度都能够很好地吻合，相同参数对应的两条曲线上的数据点之间误差不超过5%。由于通电线圈磁场沿径向变化且永磁体的尺寸相对于线圈尺寸不可忽略，因此永磁体所处的空间磁场非匀强磁场，但为了便于进行永磁体的位置排布，仍近似认为永磁体在通电线圈磁场中受到的磁力矩与其中心面和线圈轴线夹角的余弦值为线性关系。

图4 数值计算与仿真计算结果对比

2.3 执行机构优化设计

执行机构的设计包括仿生贝壳的建模设计以及永磁体位置的排布，当外部输入电流一定时，执行的设计决定动力学参数，参数优化的主要目标是令弹性系数K和驱动磁力矩在执行机构转动过程为常数且执行机构的平衡位置为转动范围的中间位置，从而将系统简化为单自由度线性系统。

仿生贝壳的结构直接决定执行机构的转动惯量以及水阻力的大小。经过以往的大量研究发展，大多数双壳贝类都是以对数螺线的形式向外生长[14]，对数螺线具有良好的结构力学和流体力学性质，本文根据对数螺线极坐标方程：

设计仿生贝壳模型，如图5 所示，式中φ为螺线定角，ρ0为贝壳长度即转轴到仿生贝壳最远端的位置，再通过描点法和变比例扫描的方式构建仿生贝壳的三维模型。为了构造多组动力学参数进行对比分析，本文选取螺旋定角为40°，贝壳长度分别为70 mm、65 mm 建立两组仿生贝壳模型，对应圆管内径分别为65 mm、60 mm。

图5 仿生贝壳模型示意图

为了使仿生贝壳张开到极限位置与管壁碰撞冲击均匀分布从而避免其受到冲击导致部分区域破坏，设计时令其上表面在张开至最大角度时能够完全贴合圆管内壁，在其底面构造贝壳轮廓边线，轮廓线与中心面的交点为贝壳最前端的点，即该点到转轴距离为ρ0，设轮廓线上其他点到中心面距离为D，要使中心面变比例扫描后得到的执行机构能够在张开最大角度时与圆管内壁贴合，则轮廓线上的点到转轴处的距离为：

式中：R为圆管内径，通过描点法可以得到贝壳的轮廓线，这也是本文仿生贝壳的结构不同于真实贝壳的一点，即轮廓线均为椭圆但长轴和短轴位置不同。将中心面上由对数螺线围成的曲面沿着轮廓边线向两侧以变比例的方式扫描生成三维壳体，从而使得贝壳表面在转动一定角度后能够与圆管内壁贴合，构建的仿生贝壳转动范围为-12°到12°。

仿生贝壳形状确定后，执行机构在圆管内的转动范围也确定，首先优化驱动磁力矩。已知永磁体在匀强磁场中受到的磁力矩和其端面与磁场方向的为余弦函数关系，由于执行机构转动范围较小，为了便于分析，近似认为与转动角度的关系曲线为余弦函数图象的一段，因此要使近似为常数，采用阵列布置永磁体的方式，在仿生贝壳靠近转轴处到中间部分以不同的初始角度对称布置永磁体，令两部分永磁体受到的驱动磁力矩分别随着转角增大而逐渐减小和逐渐增大，从而使得在转动过程中尽量保持为常数且为较大值。

为了实现两侧仿生贝壳同时张开或闭合，两侧永磁体为同极相对布置，因此从在0到12°范围内为斥力矩，想要使永磁体间磁力矩从中间位置处反向并逐渐增大，令内侧永磁体z方向位置靠近甚至越过转轴位置，使其在运动过程中部分转到轴线的另一侧，根据静磁力矩数学模型可知，两侧永磁体越过转轴的部分之间由于距离较短会产生较大的反向磁力矩，从而使得永磁体间从磁力矩在中间位置处反向并增大。

本文通过先构建执行机构三维模型后进行电磁仿真计算，根据计算结果调整模型之后再重新计算直到达到优化目标，此时的执行机构示意图如图6所示。其中模型70-65-1、70-65-2为仿生贝壳长度为70 mm时两种永磁体分布情况，模型70-65-1中内侧永磁体布置角度为-24°，外侧永磁体水平布置，模型70-65-2则与之相反，模型65-60-1、65-60-2仿生贝壳长度为65 mm，永磁体布置与70-65-1、70-65-2 相同。永磁体和3D 打印的仿生贝壳材料分别为钕铁硼NdFe35B 和树脂，密度分别为7 750 kg/m3和1 200 kg/m3，由于执行机构形状特殊不易直接计算转动惯量，由Creo软件直接给出。

图6 优化后执行机构模型示意图

和与转角之间的关系分别如图7、图8 所示，优化后的四个模型与转角基本为线性关系，且平衡位置位于中间位置，驱动磁力矩基本不随转角改变，各角度位置对应的值与平均值误差不到5%，可以认为驱动器系统弹性系数K和驱动磁力矩在转动过程中不随转角改变，优化后的动力学参数如表2所示。

图7 两侧永磁体之间磁力矩与转角的关系曲线

图8 驱动力磁力矩与转角的关系曲线

表2 系统动力学参数

3 仿生贝壳驱动器动力响应分析

仿生贝壳驱动器在谐波电流输入下会呈现周期性张合运动，而执行机构最大转角和输入电流的幅值和频率存在很大关系。已知仿生贝壳驱动器动力学响应问题属于单自由度振动问题，其动力学方程描述如下：

式中：J为执行机构转动惯量，C为系统阻尼，K为执行机构综合刚度（结构刚度和磁性刚度），A为外力矩幅值，ω为驱动力矩变化频率，θ为执行机构转动角度。

为了得到执行机构的幅频响应特性，将其转动角度表示为：

式中：Y为执行机构的转角振动幅值。

将上式代入动力学方程中，可得仿生贝壳的转角振动幅值为：

其中：HR为实频部分，HI为虚频部分。

在相频特性中，相角可表示为：

仿生贝壳驱动器的外部激励是由阵列永磁铁在通电线圈磁场中受到的磁力矩提供，可以通过调控电流的幅值和频率来控制仿生贝壳的张合运动。驱动器通过模仿贝壳喷射流体进行驱动，其驱动性能与仿生贝壳转角振动幅值和振动频率有关，因此需要研究不同输入电流下的执行机构幅频特性，对不同贝壳模型和永磁体布局下的执行机构动力学特性进行评估，从而得到最优的仿生贝壳驱动器设计结构。

本文第二节给出了四个优化后的执行机构模型以及对应的动力学参数，应用该参数进行幅频特性分析。不同输入电流下的四个执行机构模型的幅频响应如图9所示，系统阻尼C为0.001 N·s/m，由图可知，在自然条件下不同执行机构的仿生贝壳驱动器均存在一个共振峰，模型70-65-1 的共振峰最大，且共振频率向左偏移；模型65-60-2的共振峰最小。随着电流的增大，四种结构形式的仿生贝壳驱动器的共振峰均变大，但共振频率位置则保持不变。因此，在最优频率的弱电流激励下，模型70-65-1对应的仿生贝壳驱动器更容易实现较大的偏转角，有利于实现小输入电流下的系统驱动。

图9 不同输入电流下的幅频响应

由于执行机构在圆管内运动转动存在转动范围限制，其在30度偏转角限位条件下的四个执行机构的幅频响应如图10所示，可以看出存在限位的情况下，四个执行机构均在低频激励下达到限定转角，且幅频响应差异减小，但随着频率的增加，模型65-60-1的最优频率宽度最大。

图10 I=2 A时存在限位条件下的幅频响应

仿生贝壳动力学响应特性，包括实频特性、相频特性和虚频特性，如图11所示。由图11（a）可知，模型70-65-1 的实频-虚频圆最大，从图11（b）可看出，四个执行机构在外部激励为10 Hz附近时均出现相位跳跃现象，模型70-65-1 结构相位跳跃频率提前，模型65-60-1相位跳跃频率滞后。

图11 I=2 A时执行机构动力学响应特性

考虑水的阻尼作用，本文还分析了不同系统阻尼对仿生贝壳驱动器动力学响应的影响，图12给出了不同系统阻尼下的模型70-65-1的幅频响应。

图12 I=6 A时不同系统阻尼下的模型70-65-1的幅频响应

由图可知，随着阻尼的增加，共振峰出现被削平现象，共振效应减弱，当阻尼为0.003 N·s/m时，则振动幅值出现了短暂的平台效应，而后随着频率的增加迅速衰减。

4 结语

基于电磁感应和双贝类贝壳的外形结构及喷射运动机制提出了仿生贝壳驱动方案及驱动器结构，建立了驱动器执行机构的一般动力学模型并提取出动力学参数，包括弹性系数K、转动惯量J和外力矩幅值A，构建了永磁体在通电线圈中受到的磁力矩和永磁体之间的磁力矩的数学模型并通过Ansoft Maxwell 电磁仿真计算永磁体在不同位置及夹角下受到的磁力矩，数值计算与仿真结果误差不到5%从而验证了模型的准确性，依据对数螺线极坐标方程并利用Creo 建立了仿生贝壳模型并通过永磁体阵列排布优化了动力学参数，优化后的系统弹性系数K和外力矩幅值均为常数，最后完成了动力学响应计算和分析。根据结果可知，优化后的各仿生贝壳驱动器动力学幅频响应均存在一个共振峰值，且在外部激励为10 Hz 附近时均出现相位跳跃现象，当系统阻尼为0.003 N·s/m时，执行机构振动幅值出现短暂平台效应而后迅速衰减。本文研究结果可为基于仿生贝壳电磁驱动器的流量控制系统和水下航行系统提供动力学分析参考。