一种基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真方法

摘要

本发明公开了一种基于ANCF(Absolute Node Coordinate Formulation)的介电弹性体柔性鱼尾的动力学响应仿真方法，在柔性多体系统动力学理论中ANCF法的基础上建立了介电弹性体层合梁结构动力学模型，在此模型下进行动力学响应的仿真。另外，本发明利用MATLAB APP DESIGNER开发可视化应用程序，使得本领域内技术人员可以便利的通过本程序计算不同参数下的柔性鱼尾动力学响应，预测该结构的运动轨迹。

说明书

技术领域

本发明属于多体系统动力学建模领域，具体是一种基于ANCF的介电弹性体(Dielectric Elastomer,DE)柔性鱼尾的动力学响应仿真方法。

背景技术

智能材料DE驱动的柔性仿生机器鱼是典型的柔性多体系统，从柔性多体系统动力学建模理论着手，对考虑真实物理场环境和具备大摆幅游动能力的柔性仿生机器鱼系统的多场耦合动力学特性及其在软质智能材料驱动下的大变形机理进行深入研究很有必要，此外水下的流固耦合效应也是DE驱动器应用于仿生机器鱼系统动力学建模中需要考虑的因素。因此建立DE柔性鱼尾大变形动力学模型有着非常重要的意义。

陈宏等人在《基于MATLAB的机器金枪鱼巡游的动力学仿真》中基于大摆幅细长体理论，结合能量守恒原理以及水动力学理论，研究了仿生金枪鱼游动的力学机理，通过浮动坐标系建立鱼体的运动学物理模型并推导出动力学方程和效率计算公式，但这种基于浮动坐标系的方法，采用的是小变形假设，在处理柔性体较大程度变形问题上尚有局限性。Tewary等人在《Dynamic analysis ofdielectric elastomermembrane for actuation insoft fishrobots》中基于超弹性本构模型和标准流变模型建立了以双压电介电弹性体驱动的柔性鱼尾的动力学模型，研究了频率、温度、粘弹性和锥度效应对双压电介电弹性体驱动性能的关系，但是并未考虑水下流固耦合的影响。

发明内容

本发明提出了一种基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真方法。

实现本发明目的的技术解决方案为：一种基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真方法，其特征在于，包括以下步骤：

步骤1，将介电弹性体驱动的柔性鱼尾简化为层合梁模型，该模型为三层结构，上下两层为由电压驱动的介电弹性体驱动层，中间层为柔性从动层，设定介电弹性体驱动层和柔性从动层的几何参数、材料参数、网格参数，设定环境参数以及施加在介电弹性体驱动层上的驱动参数；

步骤2，基于绝对节点坐标法，选用三维二节点梁单元对层合梁模型的每层进行建模离散，并使用耦合单元模型推导出上层介电弹性体驱动层和中间柔性从动层之间的坐标转换矩阵以及下层节点弹性体驱动层和中间柔性从动层之间的坐标转换矩阵，基于Neo-Hookean超弹性本构模型和介电弹性体的力电本构模型，确定施加电压时的上下两层介电弹性体驱动层和中间柔性从动层的单元质量阵、单元刚度阵和单元弹性力列阵，结合简化流体阻力公式，确定所有单元广义外力列阵；

步骤3，根据牛顿-欧拉方程，确定每个梁单元的动力学方程，将各梁单元的动力学方程组装得到施加电压时整个介电弹性体层合梁鱼尾系的动力学方程；

步骤4，使用广义-α法作为求解算法，对系统动力学方程进行迭代求解，得到整个介电弹性体层合梁模型的位移、速度、加速度数据；

步骤5，对得到的位移、速度、加速度数据进行可视化处理，获得整个层合梁模型的运动轨迹图、末端节点竖向位移图。

一种基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真系统，基于所述的介电弹性体柔性鱼尾的动力学响应仿真方法，实现基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真。

一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时，基于所述的介电弹性体柔性鱼尾的动力学响应仿真方法，实现基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，基于所述的介电弹性体柔性鱼尾的动力学响应仿真方法，实现基于ANCF的介电弹性体柔性鱼尾的动力学响应仿真。

本发明与现有技术相比，其显著优点：(1)应用ANCF建立介电弹性体柔性鱼尾的动力学模型，更适用于计算鱼尾大变形的动力学响应，用斜率矢量代替转角坐标，不存在小变形、小转动假设，建立的方程具有常数质量矩阵，不存在科氏力和离心力，能精确地描述鱼尾的大变形行为，为柔性仿生机器鱼的动力学研究提供了新的技术方法。(2)采用的层合梁单元建模方法为分层模型，在耦合条件充分，保持层合梁单元建模精度的同时，减少了层合梁单元节点坐标数量，提高了计算效率。(3)考虑了多物理场的影响，包括考虑了智能材料介电弹性体的力电关系以及流体阻力的影响，通过改变参数设定能准确计算不同电压驱动以及有无流体阻力下，该结构的位移、速度、加速度的动力学响应。

附图说明

图1为层合梁模型示意图。

图2为层合梁电致动示意图。

图3为耦合单元边界约束条件示意图。

图4(a)为整体求解流程；图4(b)为求解方程的广义-α方法迭代流程图。

图5为上下层介电弹性体驱动电压示意图。

图6(a)为实施例考虑流体阻力的末端节点位移图；图6(b)是实施例不同电压参数下考虑流体阻力的对比图；图6(c)是实施例层合梁模型整体运动可视化图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本发明基于MATLAB对介电弹性体柔性鱼尾进行动力学响应仿真，包括以下步骤：

步骤1，设定介电弹性体驱动层和柔性从动层的几何参数、材料参数，设定环境参数以及施加在介电弹性体驱动层上的驱动参数，具体为：

(1)几何参数

上下层介电弹性体驱动层的几何参数：长为L，宽为W，高度为h

中间柔性从动层的几何参数：长为L，宽为W，高度为h

(2)材料参数

上下层介电弹性体驱动层的材料参数：密度为ρ

中间柔性从动的材料参数：密度为ρ

(3)网格参数

层合梁模型的网格参数：每一层的网格数量为n＝n

(4)环境参数

环境参数：流体密度ρ

(5)介电弹性体驱动层的驱动参数

驱动参数：上层驱动电压U1，下层驱动电压U2，电压周期T。

步骤2，基于绝对节点坐标法，选用三维二节点梁单元对层合梁模型的每层进行建模离散，并使用耦合单元模型推导出上层介电弹性体驱动层和中间柔性从动层之间的坐标转换矩阵以及下层节点弹性体驱动层和中间柔性从动层之间的坐标转换矩阵，基于Neo-Hookean超弹性本构模型和介电弹性体的力电本构模型，确定施加电压时的上下两层介电弹性体驱动层和中间柔性从动层的单元质量阵、单元刚度阵和单元弹性力列阵，结合简化流体阻力公式，确定所有单元广义外力列阵，具体方法为：

(1)基于绝对节点坐标法，选用三维二节点梁单元对该层合梁模型的每层进行建模离散：

建立全局惯性坐标系，沿中间柔性从动层的中轴线长度方向为x轴，宽度方向为y轴，厚度方向为z轴，以中间柔性从动层中轴线长度方向、宽度中线方向和厚度中线方向的交点为原点O。将上下层介电弹性体驱动层和中间柔性从动层按照网格参数，沿x轴方向离散为n＝n

建立单元局部坐标系，沿单元的中轴线长度方向为x

每个单元有两个节点，每个节点的位置在全局惯性坐标系下表示为：

r

其中S为单元形函数，为：

式中，ξ＝x

(2)使用耦合单元模型推导出上层介电弹性体驱动层和中间柔性从动层之间的坐标转换矩阵以及下层节点弹性体驱动层和中间柔性从动层之间的坐标转换矩阵：

构建耦合单元模型，需满足假设：1)耦合单元变形前后，上层介电弹性体驱动层梁单元与中间柔性从动层梁单元端面的切向平面始终共面；2)上层介电弹性体驱动层梁单元与中间柔性从动层梁单元在接触面上没有相对滑动，二者变形保持一致；

设上层介电弹性体驱动层梁单元的左端面中点为A，上层介电弹性体驱动层梁单元与中间柔性从动层梁单元左端面重合线上中点为C，C点属于上层介电弹性体驱动层梁单元左端面内；中间柔性从动层梁单元的左端面中点为A'，中间柔性从动层梁单元与上层介电弹性体驱动层梁单元左端面重合线上中点为C'，C'点属于中间柔性从动层梁单元左端面内。上层介电弹性体驱动层梁单元的右端面中点为B，上层介电弹性体驱动层梁单元与中间柔性从动层梁单元右端面重合线上中点为D，D点属于上层介电弹性体驱动层梁单元右端面内；中间柔性从动层梁单元的右端面中点为B'，中间柔性从动层梁单元与上层介电弹性体驱动层梁单元右端面重合线上中点为D'，D'点属于中间柔性从动层梁单元右端面内。

根据耦合单元模型需满足的假设，应有A-C-C’-A’四点共线、B-D-D’-B’四点共线，由此得到C与C'两点间的坐标约束方程和D与D'两点间的坐标约束方程：

其中，r

C和C'点位置处的绝对节点坐标列阵分别通过上下两层梁单元的形函数与单元坐标列阵相乘得到，表示为：

式中，S

C与C'两点间的坐标约束方程：

D与D'两点间的坐标约束方程：

其中，T

综合公式(6)和(7)，得到：

式中，T

(3)基于Neo-Hookean超弹性本构模型和介电弹性体的力电本构模型，确定施加电压时的上下两层介电弹性体驱动层和中间柔性从动层的单元刚度阵、单元弹性力列阵、单元质量阵：

考虑不施加电压时上下层介电弹性体驱动层梁单元和中间柔性从动层梁单元的弹性力列阵和刚度矩阵，不施加电压时，三者的弹性力列阵和刚度矩阵都由Neo-Hookean超弹性本构模型推导得出。

Neo-Hookean超弹性本构模型采用应变能密度函数来描述材料的各向同性、非线性、超弹性特征，其多项式应变能函数为：

U

其中C

其中r

S

定义J为变形梯度张量，表示为：

J＝[r

则右柯西-格林应变张量C的三个不变量分别表示为：

其中tr(C)为张量C的迹，λ

式中，

上下层介电弹性体驱动层和中间柔性从动层在变形过程中体积不会发生改变，即拥有体积不可压缩的特性，采用罚函数法来保证材料的不可压缩性，故体积能量罚函数为：

其中k表示不可压缩常数，取值为1×10

由式(9)和(15)得到总应变能密度函数：

对式(17)积分得到Neo-Hookean不可压缩模型的总应变能为：

由虚功原理，对式(18)进行变分，则得到Neo-Hookean不可压缩模型的单元弹性力阵，即中间柔性从动层单元的单元弹性力列阵：

将第一不变量I

变形梯度张量的行列式J对单元坐标列阵q偏导数为：

式中

将式(21)对单元坐标列阵q求导得到中间柔性从动层的单元刚度阵：

式中，变形梯度张量的行列式J对单元坐标列阵q

式(20)和(23)为不施加电压时，梁单元的弹性力与刚度阵，即中间柔性从动层的弹性力列阵与刚度阵的表达式。

当施加电压时，上下层介电弹性体驱动层梁单元的弹性力列阵和刚度矩阵要在Neo-Hookean超弹性本构模型的基础上考虑介电弹性体材料的力电本构模型。

介电弹性体的本构模型用自由能密度函数来描述，包含拉伸(U)和极化(V)两部分的贡献，其表达式为：

式中，U(λ

介电弹性体的名义电位移和名义电场存在以下关系：

将式(28)代入式(27)有

式中，C

则综合考虑了Neo-Hookean超弹性本构模型和介电弹性体的力电本构模型后，上下层介电弹性体驱动层梁单元的总应变能密度函数为：

将式(30)在整个体积上进行积分，求出上下层介电弹性体驱动层梁单元的总应变能：

由虚功原理对式(31)进行变分得到上下层介电弹性体驱动层梁单元的弹性力列阵分别为：

其中S

式中，j＝1时表示为上层介电弹性体驱动层梁单元，j＝3时表示为下层介电弹性体驱动层梁单元。

式(32)对单元绝对节点坐标列阵q

式(32)与(34)为上下层介电弹性体驱动层梁单元的弹性力列阵与刚度矩阵表达式；

上下层介电弹性体驱动层梁单元的单元质量阵为：

中间柔性从动层梁单元的单元质量阵为：

式中，ρ

(4)使用简化流体阻力公式推导出单元外力列阵

流体阻力F

其中，ρ

F

其中j＝1,3，j＝1表示为上层介电弹性体驱动层梁单元，j＝3表示为下层介电弹性体驱动层梁单元。

中间柔性从动层梁单元的广义外力列阵为：

步骤3，根据牛顿-欧拉方程直接推得每个梁单元的动力学方程，将各单元的动力学方程组装得到施加电压时整个系统的动力学方程，具体方法为：

使用坐标转换矩阵将并列的上下层介电弹性体驱动层梁单元和中间柔性从动层梁单元耦合在一起，耦合单元刚度矩阵为：

式中，

单元广义外力矩阵为：

式中，

单元质量矩阵为：

式中，M

单元弹性力列阵为：

式中，

利用牛顿方程可推出该层合梁单元模型的动力学方程为：

将所有单元的质量阵、弹性力阵、广义外力阵按节点顺序组装为结构整体的质量矩阵、弹性力阵和广义外力阵，考虑系统所受约束，采用拉格朗日乘子法引入约束方程，得到整个系统的动力学方程组：

式中，M为整体量阵，

步骤4，使用广义-α法作为求解算法对系统动力学方程进行迭代求解，得到整个介电弹性体层合梁模型的位移、速度、加速度数据，具体方法为：

广义-α法中参数初始化：

其中ρ是谱半径，表征了数值耗散程度，等于0时数值耗散最大，等于1时无耗散；h是积分步长。a的初始值是系统的起始加速度向量，λ是拉格朗日乘子向量，每一增量步结束后二者会进行更新。

在迭代计算中，第n+1步和第n步的迭代过程：

迭代步之间的误差为：

动力学方程的雅可比矩阵：

则相邻增量步之间的增量为：

于是更新位移、速度、加速度列阵与参数λ：

判断||R

若不小于限定误差，则返回至式(47)再次迭代直至满足限定误差，再返回至式(46)计算下一增量步。

步骤5：将储存起来的位移列阵进行可视化处理，具体方法如下：

在MATLAB APP DESIGNER中对位移列阵进行可视化处理，得到该三层层合梁模型的不同时刻的位形截图和末端点轨迹图。

实施例

为了验证本发明方案的有效性，开展如下实施例以验证。

本实施例中采用如表1中的参数设定，设定驱动电压为正弦电压，幅值10Kv和频率5Hz。

表1介电弹性体三层层合梁模型参数

使用如上述参数，对上下层介电弹性体层合梁施加如图5所示的错开半个周期的正弦电压，求得考虑流体阻力时的末端节点位移图，从图6(a)中可看出，在0～2s区间内，随着电压的施加，梁末端节点的横向位移逐渐增大，直至2s到达稳定状态，其后开始稳定摆动，较好的模拟了鱼尾在水中稳定摆动前进的情况。随后只改变电压幅值参数，保持其他参数不变，分别模拟了该三层介电弹性体层合梁模型在5kV、10kV、15kV电压作用下的运动，绘制了三种情况下的梁末端节点位移图并作对比，如图6(b)所示，5kV时稳定摆动幅值为m、10kV时稳定摆动幅值为m、15kV时稳定摆动幅值为m，故而梁末端节点摆动的稳定幅值随着电压幅值的增大而增大；如图6(b)所示，5kV时2.5s之后开始稳定摆动、10kV时2s之后开始稳定摆动、15kV时1.5s之后开始稳定摆动，故而整个介电弹性体层合梁模型达到稳定摆动状态的时刻随着电压幅值的增大而减小。

本发明利用绝对节点坐标法，对介电弹性体柔性鱼尾进行了动力学建模，并根据动力学方程编写了算法程序，能准确地输出在一段时间内结构各单元节点包括速度、位移、加速度等一系列动力学响应，能直观地看出整体运动轨迹变化，对于介电弹性体驱动柔性机器鱼的设计具有重要价值。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本申请范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请的保护范围应以所附权利要求为准。