基于绝对节点坐标法的压电驱动复合结构动力学特性

易灿明， 余海东， 王 皓,b

(上海交通大学 a. 上海市复杂薄板结构数字化制造重点实验室； b. 机械系统与振动国家重点实验室， 上海 200240)

软体机器人由于其动作灵活、形态多变和适用性强等特征，已成为目前机器人发展的研究热点.软体机器人的驱动主要采用复合结构，通过将硬度较大的驱动材料和柔性介质材料进行耦合，形成软硬结合的复合柔性结构.该结构具有的大变形特征使得传统的有限元建模方法难以对其运动和变形进行精确的描述[1]；同时考虑驱动器的存在，多单元耦合也给复合柔性结构的建模带来一定的难度.因此，建立精确的大柔性复合结构动力学模型对上述结构的运动控制及动力学特性分析具有重要意义.

在多单元耦合的复合结构建模方面，传统的基于浮动坐标或者转角坐标的有限元方法在建立多耦合单元时容易出现耦合条件不充分、在大变形情况下单元间出现穿透现象等问题[2].Carrera[3]总结了常用于多层板、壳结构的建模方法，其中单一板模型(ESL)假定各层板单元以单层板的形式变形，模型不够精确，且没有考虑柔性结构的大变形特征；在Layerwise模型中，通过在各层板之间插入高阶边界条件，提升了模型的精度，但同时也极大地增加了结构的自由度数量.Zhao等[4]基于不同单元之间的变形协调条件，提出了单元间的坐标转换矩阵，不仅提高了耦合建模的精度，同时降低了耦合单元的自由度数量.在处理复合结构大变形过程中刚体运动与柔性变形的耦合问题方面，Shabana等[5]提出绝对节点坐标法为复合结构建立高阶单元模型，并基于该方法建立了三维梁单元模型，进而运用连续介质力学方法[6]对单元的刚度矩阵进行了推导.绝对节点坐标法将单元节点坐标定义在同一个全局坐标系下，同时使用节点的斜率坐标代替传统的转角坐标对节点的方向进行描述，在解决结构大变形、大转动问题时具有更高的精度.随后，Mikkola等[7]通过该方法为三维板单元和壳单元建立了非增量有限元模型；赵春璋等[8]基于绝对节点坐标法建立了变截面梁单元模型，并针对截面特性对单元动力学性能的影响进行了研究；李彬等[9]通过大变形柔性梁系统验证了用绝对节点坐标法建模的正确性和精确性.在压电材料驱动器建模方面，Gilardi等[10]基于绝对节点坐标法建立了由压电驱动器驱动的大柔性梁模型，并研究其在高速旋转及振动控制下的静力学和动力学特性；Nada等[11]基于绝对节点坐标法对多层含压电薄片的薄板单元进行了有限元建模，并对其运动特性进行了分析.目前关于压电材料驱动器的研究大多通过传统的层合理论对压电驱动器和主体单元进行耦合，单元耦合建模精度有待提高.

鉴于绝对节点坐标法采用的高阶单元模型可对具有大变形特征的柔性复合结构运动过程中的刚柔耦合特性进行精确描述，同时考虑到耦合单元精确建模的需求，本文基于绝对节点坐标法，根据变形协调条件推导坐标关联矩阵，建立了大柔性梁板耦合单元模型，并对其动力学方程进行了求解.此外，引入压电材料本构方程，将梁单元作为复合结构的驱动器，用Matlab软件对压电驱动柔性板结构进行数值仿真分析，研究该结构在外加驱动下的静力学和动力学特性，并分析各项参数对其性能的影响，以期为软体机器人的结构及驱动设计提供参考.

1 基于绝对节点坐标法的柔性单元建模

1.1 单元的变形描述

根据绝对节点坐标法，柔性板单元的所有坐标均定义在全局坐标系下，三维四节点板单元模型如图1所示.图中：O-XYZ为全局坐标系；o-xyz为与板单元固连的浮动坐标系，描述板单元上任意一点的相对位置.

图1 三维四节点板单元模型Fig.1 Three-dimensional four-node plate element

该板单元模型包含A、B、C和D4个节点，每个节点由3个位置坐标和9个斜率坐标进行描述，板单元上节点i(i=A,B,C,D)的坐标可定义为

(1)

式中：ri为节点i的全局位置坐标；∂ri/∂k(k=x,y,z)为节点i处的浮动坐标系相对于全局坐标系在k轴上的斜率坐标.由于每个节点包含12个绝对节点坐标，所以整个板单元包含48个绝对节点坐标.

通过引入高阶插值函数，板单元上任意点P的位置在全局坐标系下可表示为

r=Se

(2)

式中：e为基于绝对节点坐标法建立的单元节点坐标；S为该单元对应的形函数矩阵，

S=[S1IS2I…S16I]

(3)

I为3×3的单位矩阵，且有

(4)

1.2 板单元质量矩阵

将式(2)对时间求导，可得到板单元上任意点的速度表达式：

(5)

因此板单元的动能可表示为

(6)

式中：ρ为材料密度；V为梁单元体积；M为单元常数质量矩阵，

(7)

1.3 板单元刚度矩阵

基于连续介质力学，板单元的变形梯度矩阵可表示为

(8)

(9)

式中：

(10)

由于该应变张量为对称的二阶张量，因此可将其改写成向量形式：

(11)

式中：

(12)

对于弹性体，根据广义胡克定律，其应变张量与应力张量存在如下对应关系：

σ=cE:ε

(13)

式中：σ为单元应力张量；cE为材料弹性系数张量.根据连续介质力学，对于各向同性的固体材料，弹性系数张量可以通过拉梅常数进行简化表达.简化后的弹性系数张量为

(14)

单元刚度矩阵可通过应变能函数求得.应变能函数表达式为

(15)

式中：V0为单元在初始构型下的体积.式(15)对单元节点坐标求偏导，可得单元广义弹性力为

(16)

由此可推导出刚度矩阵表达式为

K=(λ+2μ)K1+λK2+2μK3

(17)

式中：

(eTSce-1)Sc1]dV

(eTSbe-1)(Sa1+Sc1)]dV+

(eTSfe)Sf1]dV

(18)

1.4 板单元的动力学方程

基于Hamilton原理建立板单元动力学方程：

(19)

式中：Wa为广义外力做功.对方程中各项进行展开可得到简化后的柔性板单元动力学模型：

(20)

式中：Fa为系统广义外力阵，可由虚功原理推导得到，

(21)

F为系统所受广义外力.

2 压电驱动单元及复合柔性板单元建模

为了引入软体机器人的驱动方式，考虑在柔性板单元基础上增加压电驱动梁单元.

2.1 驱动梁单元变形描述

驱动梁单元采用三维两节点梁单元，其几何模型如图2所示.

图2 三维两节点梁单元模型Fig.2 Three-dimensional two-node beam element

梁单元包含Ax和Bx节点，根据绝对节点坐标法，其坐标定义在全局坐标系下，每个节点包含12个坐标，因此单元共包含24个绝对节点坐标：

(22)

通过引入高阶插值函数，可得到梁单元对应的单元形函数矩阵：

(23)

且有

(24)

2.2 耦合单元变形描述

复合柔性板结构由板单元和梁单元耦合而成，其几何模型如图3所示.

板单元与梁单元的耦合通过变形协调条件实现，满足2点基本假设：① 耦合单元变形前后，板单元与梁单元端面的切向平面始终共面；② 在运动过程中，板单元与梁单元的接触面不发生相对滑移.基于这2点基本假设，以图4所示的AD端面为例，A′、G、G′和Ax4点始终共线，并且G与G′点始终重合.

根据上述假设，由变形协调条件可得

图3 柔性梁板耦合单元模型Fig.3 Flexible composite beam-plate element

图4 耦合单元AD端截面示意图Fig.4 Cross section at side AD of the coupled element

(25)

G点的位置和斜率坐标可通过形函数矩阵表示为

(26)

同理，对于梁单元上的G′点有

(27)

综合式(25)～(27)可得

(28)

对于耦合单元BC端面，其示意图如图5所示.

图5 耦合单元BC端截面示意图Fig.5 Cross section at side BC of the coupled element

同样，由F和F′点处的变形协调条件可得

(29)

综合式(28)和(29)，有

(30)

(31)

式中：RBP为关联板单元与梁单元的坐标转换矩阵.通过该坐标转换矩阵，可用板单元的绝对节点坐标描述梁单元节点坐标，在保证耦合单元建模精度的同时，大大地缩减了耦合单元的坐标数量.

2.3 压电驱动材料本构关系

本文使用压电材料作为复合柔性单元的驱动器材料.压电材料本构关系可表示为

(32)

式中：eE为材料压电耦合系数；Q为施加的电场强度；D为电位移；ζ为相对介电常数.将本构方程展开可得

由材料的对称性有

c11=c22,c12=c21,c44=c55

c13=c23=c32=c31

e31=e32,ζ11=ζ22

由该本构方程可得驱动单元弹性势能表达式：

(33)

电场力做功为

(34)

2.4 复合板结构动力学方程

压电梁驱动器的质量矩阵由梁单元形函数矩阵推导得出:

(35)

基于变形协调条件得出的坐标转换矩阵，可得复合板单元的质量矩阵为

MBP=M+(RBP)TMBRBP

(36)

压电梁驱动器刚度矩阵可由式(33)对节点坐标2次求偏导得到，

(37)

电场力做功对节点坐标求偏导可得电场力表达式为

(38)

基于坐标转换矩阵可得复合板单元的刚度矩阵为

KBP=K+(RBP)TKBRBP+(RBP)TKERBP

(39)

复合板结构动力学方程的建立同样基于Hamil-ton原理，其表达式为

(40)

对式(40)进行展开并引入质量矩阵式(36)、刚度矩阵式(39)及广义外力表达式可得复合柔性板结构的动力学方程为

(41)

式(41)所示的代数微分方程可采用Runge-Kutta方法进行数值求解.本文采用商用软件Matlab中的ode45函数进行求解，该函数是一种自适应步长的Runge-Kutta数值解法，适用于求解具有非刚性特性的常微分方程.

3 压电驱动板结构数值仿真分析

3.1 静力学特性分析

数值仿真分析所选用的复合柔性板结构如图6所示.板的AD端完全固定，BC端自由，为悬臂结构.柔性板长度、宽度和厚度分别为 0.6、0.6 和 0.01 m.压电梁驱动器下表面与柔性板上表面贴合，Ax和Bx点分别位于柔性板AD端和BC端中点上方.压电梁驱动器长度、宽度和高度分别为 0.6、0.01 和 0.01 m.该结构纵向截面如图7所示.压电梁驱动器上下表面接通电源，从而产生压电驱动力.

图6 复合柔性悬臂板结构Fig.6 Cantilever plane with actuator element

图7 压电驱动板结构截面示意图Fig.7 Cross section of the plane with actuator element

本文的静力学计算案例中，柔性板的材料参数：E=100 MPa；ν=0.3；ρ=1 300 kg/m3.压电梁驱动器选用压电陶瓷材料(PZT)，其材料属性可由静态法试验测得，具体数值如表1所示.

表1 压电梁驱动器材料属性Tab.1 Material properties of piezoelectric beam

图8 不同驱动电压下柔性板的变形Fig.8 Deformations of the plane at different voltage

静力学计算中，在柔性板和驱动梁尺寸及材料参数不变的情况下，对驱动梁分别施加U=10，20，50，100，200 V的电压.柔性板在不同驱动电压值下的变形如图8所示.由图可知，驱动梁上施加不同的电压值时，柔性板产生不同程度的变形.这是由于压电材料具有特殊的逆压电效应，施加在压电驱动梁上下表面的电压产生的电场使得梁结构在厚度方向发生相应的结构变形，从而使得与其相连的柔性板随着压电驱动梁的结构变形而发生弯曲变形，且该变形量随驱动电压的增大而逐渐增大.

基于上述静力学计算结果，在E=100，50，10 MPa的3种不同条件下计算柔性板末端B点Z方向位移(dBZ)随驱动电压变化的情况，结果如图9所示.由图可见：柔性板在压电驱动作用下产生弯曲变形，柔性板末端位移量随驱动电压的增大呈增大趋势；柔性板末端位移增大幅度随弹性模量的减小而变大.

图9 不同驱动电压下柔性板末端B点Z方向位移Fig.9 Z-displacement of node B at different voltage

3.2 动力学特性分析

以图6所示的带驱动器的柔性板为例，分析不同电压参数及材料参数对板结构动力学特性的影响.柔性板长度、宽度和厚度分别取 0.6、0.6 和 0.01 m，驱动器长度、宽度和高度分别取 0.6、0.01 和 0.01 m.

压电材料参数见表1，驱动电压分别取50和100 V，柔性板的材料参数同前(E=100 MPa，ν=0.3，ρ=1 300 kg/m3).仿真时间为 0.1 s，输出柔性板末端B点X、Y、Z方向位移.由于仿真结果数据量庞大，所以只在仿真时间区间内均匀选取其中的30个点，且保证相邻点时间间隔相同，得到B点3个方向位移随时间(τ)变化的曲线如图10所示.由图可见：当对柔性板施加驱动电压时，板末端产生周期性的上下摆动；B点在X和Y方向上位移量较小，而Z方向的位移以近似余弦规律随时间变化；驱动电压由50 V增加至100 V时，末端位移明显增大，且末端运动规律相同，振动周期不变.

图10 驱动电压为50和100 V时B点位移随时间变化曲线Fig.10 Displacement of node B at 50 and 100 V

基于上述动力学计算结果，在E=100，50，10 MPa，U=100 V的条件下，对柔性板进行计算，以研究不同材料参数对板结构动力学特性的影响，计算得到的板末端B点Z方向位移随时间变化的曲线如图11所示.由图可见，随着柔性板弹性模量的减小，板末端在驱动电压作用下产生的位移量增大，且当弹性模量发生变化时，柔性板末端上下振动的幅度以及振动周期均发生改变.B点在不同弹性模量下的振动幅度(dBmax)及振动周期(T)如表2所示.由表可知，驱动电压不变时，柔性板末端B点Z方向位移随弹性模量的减小呈增大趋势，且振动周期相应增大.

图11 驱动电压为100 V时B点Z方向位移Fig.11 Displacement of node B at 100 V

表2 柔性板B点在不同弹性模量下的振动特性Tab.2 Kinetic statistics of node B under actuation

4 结论

本文基于绝对节点坐标法，对柔性板单元动力学方程进行了推导；同时引入压电材料本构方程，将梁单元作为驱动器，通过变形协调条件将压电梁驱动器和柔性板耦合，形成带驱动器的柔性复合板结构，并推导出该结构的动力学方程.通过对复合板结构进行静力学和动力学的仿真计算，分析了驱动电压参数及材料参数对复合板结构动力学特性的影响，得到的主要结论如下：

(1) 静力学计算中，当驱动梁上施加电压时，柔性板在驱动电压作用下产生弯曲变形.随着驱动电压的增大，柔性板弯曲变形程度增大，板末端位移量呈近似线性增大的趋势.

(2) 弹性模量对复合板结构的静力学特性具有一定影响.材料弹性模量越小，结构末端在驱动电压作用下产生的位移增幅随驱动电压的增大而增大.相同驱动电压条件下，柔性板末端位移随弹性模量的减小而非线性增大.

(3) 动力学计算中，当驱动梁上施加持续电压时，柔性板末端产生上下周期性的振动.随着驱动电压的增大，板末端振动的幅度相应增大，振动规律基本不变.当驱动电压一定而弹性模量减小时，板末端的振动幅度增大，同时振动周期也相应增长.