基于几何精确Euler-Bernoulli梁单元的柔顺机构动力学分析

张志刚，周翔，毛红生，王胜永

郑州轻工业大学 河南省新能源汽车轻量化设计与制造工程研究中心，河南郑州 450002

0 引言

柔顺机构利用柔性部件的大变形代替运动副来传递运动、力和能量[1]，具有部件数量少、易装配、质量轻、无摩擦磨损、振动噪声小、可实现微型化等优点，在微机电系统、微执行机构、精密定位等工程领域具有巨大优势和潜力[2].

柔顺机构的工作机理决定了其在运行过程中必然伴随柔顺杆件的大变形，这种强几何非线性特征使柔顺机构的分析设计较为困难.为了简化大变形柔顺杆的建模，L.L.Howell等[3]针对柔顺机构设计过程中如何快速计算柔顺杆端部位移的问题，提出了伪刚体模型方法，该方法将大变形柔顺杆件简化等效为两段由旋转副和扭转弹簧连接的刚性杆件，从而可以利用刚性机构分析理论和方法对柔顺机构进行分析.该方法在柔顺机构学研究领域被广泛关注，因等效模型中含有一个转动副(Revolute Joint)，因此又被称为1R伪刚体模型[4].虽然1R伪刚体模型能在一定范围内给出具有一定精度的杆端位移解，但精度较低.为了进一步提高1R伪刚体模型的精度，一些学者尝试通过增加等效刚体段和转动副个数或考虑轴向变形对弯曲挠度的影响对刚体模型进行了改进[5-9].

在满足静力学和运动学要求的基础上，如何建立柔顺机构精确数学模型进而改善和提高系统动力学特性，已经成为当前柔顺机构研究与应用领域中亟待解决的问题.为此，一些学者尝试将伪刚体模型方法推广到柔顺机构动力学问题的研究中[10-12].但伪刚体模型方法本质上仍是采用具有等效力-杆端位移关系的刚体构件模拟柔顺杆件，缺乏对柔顺杆大变形细节的描述能力，这也进一步制约了其在柔顺机构精确动力学仿真分析中的应用.

柔顺杆件在力学上可以抽象为大变形柔性梁，因此可将多体系统动力学研究领域中针对大变形梁所提出的建模方法应用于柔顺机构精确建模与动力学分析.其中最具有代表性的两类大变形柔性梁建模方法是绝对节点坐标 (Absolute Nodal Coordinate Formulation,ANCF) 方法和几何精确梁理论(Geometrically Exact Beam Theory,GEBT).ANCF梁单元[13-14]以总体惯性坐标系下节点的位置矢量坐标和梯度矢量坐标作为单元参数，具有单元质量矩阵为常数阵、系统方程不含向心力项和科氏力项等优点，被广泛应用于航空航天、车辆、机器人等工程领域.李鹏飞等[15]采用ANCF方法研究了平面固定-导向柔顺机构的变形与驱动力变化规律问题.张志刚等[16]考虑柔顺杆端部变形特征，提出了用于柔顺机构动力学建模的端部受约束ANCF梁单元.但ANCF梁单元节点自由度过多，这会导致系统动力学方程维数过高，从而影响计算效率.此外，ANCF梁单元还存在泊松闭锁、剪切闭锁等问题.

GEBT[17-18]在刚性截面假设基础上引入大转动参数描述截面转动，创造性地定义了不受大范围刚体运动影响的客观性应变度量，因而适用于包含大变形柔性梁刚柔耦合问题的建模与求解.目前该方法已经广泛应用于各类工程：廖明夫等[19]采用几何精确梁单元建立了风电组叶片非线性耦合模型，研究了非定常流场作用下叶片载荷波动问题；张健等[20]对螺旋弹簧变形进行整体插值，采用GEBT研究了弹簧刚度问题；陈进格等[21]基于几何精确梁模型建立了风力机叶片气弹模型，研究了风轮在正对来流条件下的气弹响应.胡建峰等[22]利用几何精确梁单元法对周边桁架天线的横杆、纵杆及索网进行建模，对不同传动效率和平均张力条件下的天线展开过程进行了动力学仿真.

目前基于GEBT的柔顺机构相关研究还鲜有报道.而大变形柔顺杆件属于细长梁，满足Euler-Bernoulli梁变形特征.鉴于此，本文拟基于几何精确Euler-Bernoulli梁单元[23]研究柔顺机构的动力学建模与仿真求解问题，采用严格满足Euler-Bernoulli梁变形假设的大变形几何精确梁单元，建立柔顺机构动力学模型并进行数值仿真计算，以期在保证计算精度的同时提高计算效率.

1 大变形Euler-Bernoulli梁单元的变形场插值

根据刚性截面假设，梁模型可以看作由无数个沿轴线分布刚性截面构成的多刚体系统，梁的运动和变形由截面形心运动和绕形心转动描述.因此，构造大变形梁单元的关键是能对梁形心线变形和截面转动提出合理的近似.为了满足传动需求，柔顺机构中的柔顺杆件一般被设计成大长细比，其变形特征符合Euler-Bernoulli 梁假设，即形心线切向始终保持与截面垂直.空间Euler-Bernoulli梁的运动学描述见图1.

图1 空间Euler-Bernoulli梁的运动学描述

由于梁主要承受横向载荷作用，变形也以弯曲为主，轴向应变通常很小，即ε<<1.因此，在构造Euler-Bernoulli梁单元时可以合理地忽略单元两端节点处的轴向应变ε1和ε2对形心线形状的影响，将整体惯性坐标系下截面形心矢径r近似为[23]

①

其中归一化坐标ξ=s/L，L为梁单元初始长度.

根据Euler-Bernoulli梁变形假设，梁截面单位法向矢量n可以由①式表示为n=r′/|r′|，其中r′=∂r/∂s代表对弧长s的导数.因此，要确定梁截面坐标系{n,t,b},只需确定另外两个基矢量t和b.为此，可以将梁截面坐标系{n,t,b}看作由左端梁截面坐标系{n1,t1,b1}经过两次连续转动得到，如图2所示.

图2 梁截面的转动

β=(β/sinβ)n1×n

然后将中间参考坐标系{n,t,b}绕着n方向转动α角得到最终截面坐标系{n,t,b}，其对应转动矢量为

α=αn

由上述梁截面坐标系的转动分解，梁截面坐标系方位矩阵可以表示为

Rθ=[ntb]=RαRβR1

②

由上述插值过程可知，上述单元变形场严格满足Euler-Bernoulli梁变形假设，因此不会出现剪切闭锁问题.

2 单元节点力

根据GEBT，大变形Euler-Bernoulli梁的内力虚功率δps可以表示为

③

其中，E为材料弹性模量；G为剪切模量；J,It和Ib分别为截面极惯性矩和两个截面惯性矩；ε为单元轴向应变；κn,κt和κb为曲率矢量κ在梁截面坐标系{n,t,b}下的分量，即κ=κnn+κtt+κbb.

④

单元轴向应变ε可由梁单元形心线插值函数①式表示为

ε=nTr′-1

⑤

将④式和⑤式代入③式，可得：

⑥

其中，单元广义节点力列阵Fs为

3 单元惯性力虚功率

几何精确梁模型保留了经典梁理论中的刚性截面假设，据此可以将梁看作由无数个沿弧长方向排列的刚性截面组成的特殊多刚体系统.因此，梁的惯性力虚功率可以表示为

⑦

⑧

⑨

⑩

其中

Wω=Tα(∂α/∂e)+RαTβ(∂β/∂e)+RαRβTφ1Dφ1

对⑩式关于时间求一次导数可以得到截面角加速度矢量：

其中，单元质量矩阵Mi和惯性力项Fi分别为

4 单元控制方程

作用于梁单元上的常见外力包括节点处所受的集中力外力F1，F2,集中力矩M1，M2,以及包括重力在内的分布力p(s)，则外力虚功率δpe可以表示为

其中，广义外力为

上述推导得到的大变形几何精确Euler-Bernoulli梁单元控制方程可以对柔顺杆件进行有限元离散，并按照自由度排序组装形成系统整体方程.由于描述柔顺机构系统运动学参量包含慢变的整体刚体运动学量和快变的弹性运动学量，因此最终形成的系统刚柔耦合方程需要采用刚性微分方程求解器仿真求解.

5 仿真结果与分析

为验证几何精确Euler-Bernoulli梁单元求解柔顺机构动力学问题的有效性，对典型算例进行数值仿真.在Matlab中编写了几何精确Euler-Bernoulli梁单元相关程序代码，并将得到的数值仿真结果记为GEBT/Euler.为了对比计算精度和计算效率，编写了三维ANCF梁单元[11]相关程序，并将仿真结果记为ANCF.此外，在机械系统仿真软件ADAMS中采用大变形梁类模块FE_Part搭建了相应算例模型，并将仿真结果记为ADAMS.

5.1 柔顺四连杆机构

图3 柔顺四连杆机构

分析可知，柔顺四连杆机构的运动和变形均发生在平面XY内.在搭建动力学方程时，刚性杆AB和BC采用刚体建模，柔顺杆CD均匀离散为5个单元，最终形成的系统方程采用Matlab刚性微分方程求解器ode15s求解，设定仿真时间为1 s.图4给出了采用不同方法时柔顺杆顶端C点位置坐标的时间变化曲线.

图4 柔顺杆顶端C点位置坐标的时间变化曲线

由图4可以看出，采用5个几何精确Euler-Bernoulli梁单元对柔顺杆进行离散便能得到与ADAMS十分吻合的数值结果，这初步表明了几何精确Euler-Bernoulli梁单元处理柔顺机构动力学问题的有效性；0～0.3 s时间段内，ANCF梁单元的数值结果与几何精确Euler-Bernoulli梁和ADAMS结果吻合较好，0.3 s后会出现一些偏差，究其原因,是由于ANCF梁单元的闭锁问题.ANCF梁单元闭锁问题是由于变形函数中存在弯曲与剪切的强耦合，当单元发生较大弯曲变形时会伴随产生过大虚假剪切变形，从而使单元表现过于刚硬.由图4还可以看出，前0.3 s柔顺杆弯曲变形较小，ANCF单元闭锁影响较小；随着柔顺杆弯曲变形的增大，在一个ANCF梁单元内产生了较大弯曲变形，对应的闭锁问题影响增大.因此在利用ANCF梁单元进行柔顺机构动力学分析时，必须考虑闭锁问题的影响.本文采用的几何精确Euler-Bernoulli梁单元可以避免剪切闭锁问题，适用于大变形柔顺机构的动力学仿真分析.

采用ANCF梁单元与几何精确Euler-Bernoulli梁单元进行柔顺四连杆机构动力学仿真计算的CPU时间分别为4 673.2 s和592.1 s，几何精确Euler-Bernoulli梁单元的计算效率明显高于ANCF梁单元.这是由于三维ANCF梁单元节点包含了节点位置矢量坐标和3个位置梯度矢量坐标，共12个参数，而几何精确Euler-Bernoulli梁单元的节点参数为节点位移矢量坐标和转动矢量坐标，共6个参数.相较于 ANCF 梁单元，采用几何精确Euler-Bernoulli梁单元可以明显降低系统参数维度.

5.2 空间圆弧导向柔顺机构

空间圆弧导向柔顺机构由柔顺杆OA和刚性杆AB构成(见图5)，其中柔顺杆两端分别与大地和刚性杆固结，刚性杆的另一端在B点与大地铰接.柔性杆OA与刚性杆AB之间的夹角为∠OAB=45°；刚性杆AB长为LAB=0.6 m，密度为ρAB=7860 kg/m3，截面面积为AAB=(0.06×0.06) m2；柔顺杆OA长LOA=1.2 m，截面面积为AOA=(0.04×0.04) m2，密度为ρOA=1000 kg/m3，弹性模量为E=1×107Pa，泊松比为v=0.27.该机构在铰接处B受到转矩M=15sin(πt/2) N·m的作用下，由图示初始位置开始运动.

图5 空间圆弧导向柔顺机构

在空间圆弧导向柔顺机构运行过程中，柔顺杆端A点的轨迹为一段圆弧.在这个算例中，刚性杆AB采用刚体建模，柔顺杆OA同样均匀地划分为5个单元，最后形成的系统动力学方程采用Matlab刚性微分方程求解器ode15s求解，仿真时间设定为1.2 s.图6给出了采用不同方法时柔顺杆端A点位置坐标随时间的变化曲线.

图6 柔顺杆端A点位置坐标随时间变化曲线

由图6可以看出，采用几何精确Euler-Bernoulli梁单元得到的结果与ADAMS几乎完全吻合；而由于存在单元闭锁问题，0.8 s以后，采用ANCF梁单元所得仿真结果与几何精确Euler-Bernoulli梁单元及ADAMS的数值结果相比，开始出现一定偏差.该算例进一步验证了几何精确Euler-Bernoulli梁单元对空间柔顺机构动力学问题具有良好的适用性和计算精度.

采用ANCF梁单元与几何精确Euler-Bernoulli梁单元对空间圆弧导向柔顺机构进行建模仿真计算所需CPU时间分别为3 390.1 s和198.6 s，几何精确Euler-Bernoulli梁单元对空间柔顺机构进行动力学仿真时的计算效率显著高于ANCF梁单元.

6 结论

应用GEBT研究了大变形柔顺机构的动力学建模及仿真求解问题，基于Euler-Bernoulli梁变形假设推导了适用于大变形细长梁建模的几何精确梁单元，有效避免了剪切闭锁问题的出现；以柔顺四连杆机构和空间圆弧导向柔顺机构作为典型算例，应用几何精确Euler-Bernoulli梁单元，建立动力学模型.与ANCF梁单元模型和ADAMS的数值结果进行了对比分析，验证了几何精确Euler-Bernoulli梁单元应用于大变形柔顺机构动力学问题时具有计算效率和计算精度上的优越性.