含非线性大变形构件的柔顺机构建模与分析

李鹏飞， 曹博宇， 汪振宇， 赵 晨， 王立鹏

(西安理工大学 机械与精密仪器工程学院， 西安 710048)

传统刚性机构借助于刚性构件和运动副实现运动和功能。柔顺机构则是一种利用构件自身弹性变形完成运动和力传递及转换的新型机构[1],具有减少构件数量和装配时间、简化加工工序、无摩擦磨损和传动间隙、降低振动和噪声等优点[2]，可应用于精密测量与定位、微机电系统(MEMS)、仿生机械、医疗器械等领域[3]。柔顺机构运动往往伴随着柔性构件的非线性大变形过程，须用非线性方程描述构件大变形造成的几何非线性，使得柔顺机构建模变得困难。

目前柔顺机构的建模方法主要有伪刚体模型法[4-5]、椭圆积分法[6-7]和有限元法[8]等，其中伪刚体模型法和椭圆积分法较为常用，计算效率较高。伪刚体模型法将柔顺杆件离散成为用扭簧铰接的刚性杆件，柔顺机构建模问题转化为刚性机构建模，这种转化是一种近似的转化，误差较大。Zhang等通过引入梁上拐点个数和梁末端弯矩方向等参数，运用椭圆积分法求解了梁在各种不同载荷条件下的大变形问题。王雯静等基于有限元法建立柔顺机构的动力学方程，但其梁单元基于小变形假设，对于描述大变形柔顺机构具有一定的局限性。Berzeri等[9]应用ANCF法分析了具有大变形构件的平面四杆机构。

本文简要介绍了椭圆积分法和ANCF法的基本原理，在比较两种方法建模过程基础上，对一种固定—导向柔顺梁的变形过程和受力进行仿真，设计实验装置，验证梁变形仿真结果的正确性。应用ANCF法对一种复杂柔顺双稳态机构进行建模分析，对比分析两种建模方法的特点。

1 柔顺机构建模

1.1 椭圆积分法基本原理

在现有的柔顺机构建模方法中，常用椭圆积分法计算梁的大挠度变形。图1所示为长度为L的悬臂梁在初始水平(虚线位置)、末端同时受到水平方向力Fx、竖直方向力Fy和弯矩MO共同作用时的变形示意图，坐标系xOy如图所示。根据Euler-Bernoulli梁理论，梁上任意一点A(x,y)处的弯矩M为

(1)

根据式(1)，可以推导得到最基本的椭圆积分法大挠度梁的支配方程为

(2)

式(2)中忽略了梁的轴向变形以及截面的剪切作用。通过椭圆积分求解式(2)时，需先假定梁上拐点的个数和梁末端受到的弯矩方向，会出现多解，且对未知数的求解次序有一定的要求，不便于计算求解。

1.2 ANCF法基本原理

1.2.1 ANCF有限单元

ANCF法推导建立的多体系统微分代数方程具有常数质量矩阵、不存在科氏力和离心力项等特点，适合于描述大转动大位移、大变形柔性多体系统。

图1 末端受力和弯矩作用的悬臂梁变形示意图

如图2所示，将柔顺机构构件离散为三维二节点梁单元，该梁单元上任意一点的全局位置表示为[11-12]

r=S(x,y,z)e(t)

(3)

式中:S为单元形函数矩阵;e为单元节点坐标向量。

图2 三维二节点梁单元

每个梁单元含有2个节点i、j，每个节点有12个绝对坐标，节点坐标e(t)定义为

(4)

1.2.2 单元动力学方程

将式(3)对时间求导，由单元的动能表达式可写出单元质量矩阵Me为

(5)

式中:ρ为材料密度;V为单元体积。

根据连续介质力学中非线性应变-位移关系可计算单元的广义弹性力。变形梯度J写为

(6)

式中:J0=∂X/∂x;X=Se0;x=[x，y，z]T，e0为单元初始位置时的节点坐标向量。则Green-Lagrange应变张量εm写为

(7)

式中:张量εm包含梁单元的轴向拉压、扭转、弯曲应变以及截面的剪切作用;I表示单位矩阵。考虑到应变张量εm的对称性，将其各元素写为向量形式ε。由应变向量ε和应力向量σ间关系可以求出单元的应变能为

(8)

将应变能U对节点坐标e求偏导，可获得弹性力向量

(9)

由式(5)、(9)，无约束条件下梁单元的动力学方程为

(10)

式中，Qe为单元外力向量。

将单元动力学方程可进一步装配成系统动力学方程。

1.2.3 弹性力的物理意义

由式(9)，单元弹性力Qs向量可写为

(11)

式中:对于单元的节点i，fi为作用在节点i上的力;fix、fiy和fiz与施加在截面上的力矩有关。当施加在截面上的力矩为T=[TxTyTz]T时，有下式成立[14]

(12)

将式(12)直观表达在单元上，如图3所示。

图3 单元弹性力分量

图4所示为柔顺机构的建模求解过程，可以清晰地得到本文所介绍的两种建模方法的应用流程。

图4 柔顺机构建模过程

1.2.4 基于OpenMP并行编程

以上计算过程中，质量矩阵是常数矩阵，弹性力矩阵是时变的，且多次采用高斯-勒让德积分，降低了计算效率。本文选用Intel visual fortran语言编制仿真程序，使用8核CPU的 DELL服务器硬件，在编制程序时采用基于OpenMP共享存储并行编程的并行计算思路，使得计算效率极大提高。

2 柔顺机构算例

2.1 固定-导向柔顺梁

如图5(a)所示的全柔顺双稳态机构两侧对称，共包含4根柔顺梁。为了简化分析，取其中如图5(b)所示的固定-导向柔顺梁作为分析对象，图中A端固定，另外一端B垂直向下运动，且B端截面转角保持不变。梁AB初始倾角为β=30°，长度L=30 mm，矩形截面尺寸b×h=1×0.2 mm，弹性模量E=1.6×1011Pa，泊松比ν=0.3。设B端位移按式(13)函数规律加载。

(a)全柔顺双稳态机构(b)固定-导向柔顺梁

图5 固定-导向柔顺机构

Fig.5 Fix-guided compliant mechanism

(13)

式中:Vs为B端的最终速度，Ts为速度从0加载到Vs所用的时间，t为加载时间。

图6所示为分别采用椭圆积分法和ANCF法计算固定-导向柔顺梁运动时梁轴线的变形结果。从图6中看出：两种方法得到的变形基本吻合。梁AB从起始位置开始运动，并发生屈曲。初始阶段，梁上有两个拐点，随着B端位移增大，靠近固定端A的拐点逐渐趋近A端直至消失，此后梁上只剩一个拐点。拐点数量的变化揭示了此结构两种稳态的跳转变化过程。需强调的是，在求解时，相较于椭圆积分法，ANCF法无须假设梁上拐点数目和末端弯矩方向。

采用ANCF法将柔顺梁划分为10个单元，如图7所示为距梁B端3 mm处点P(如图5(b))在x方向的位置随时间变化曲线。图中可以看出：正常加载时，0.03 s之后点P在x方向产生了剧烈的振动，这是因为0.03 s之前梁主要在轴向产生压缩，在大约0.03 s时发生了屈曲，突然弯曲产生了剧烈的振动。在0.22 s时梁上两个拐点变为一个拐点，梁的形态发生了跳转。

应用ANCF法分析该柔顺梁时，发现梁在屈曲的瞬间产生了较大的振动，而采用椭圆积分法则无法获得该振动现象。为了验证图7中的结果，在如图5(b)的Q点，梁屈曲前施加一个微小竖直向上的初始干扰，得到如图7所示初始加干扰的P点x方向的位置变化曲线。可以看出，由于施加的初始干扰与梁屈曲变形方向一致，使得屈曲变形变得容易，则梁的振动明显减小。当施加向下的干扰时，柔顺梁轴线变形变成如图8所示的另一种形态。

图6 固定-导向柔顺梁轴线变形

图7 点P在x方向上的位置变化

为验证图6和图8所示柔顺梁仿真计算结果，设计了如图9所示实验装置。图9中采用XYZ平台作为加载装置，弹性梁3通过夹具分别与Z轴下端1、固定底板2固接。

(a)实验台(b)实验台实物

图9 固定-导向梁变形实验台

Fig.9 Experiment facility of the fixed-guided beam

XYZ平台的Z轴由控制器控制，按式(13)规律垂直运动，实现柔顺梁一端运动、一端固定，同时自身发生变形。使用摄像机记录施力过程中柔性梁的变形。ANCF法的数值仿真结果和实验结果对比如图10所示，其中虚线为数值结果，实线为实验结果。可以看出理论和实验结果基本完全吻合，验证了理论计算结果的正确性。

根据式(12)，采用ANCF法求得图5(b)中B端运动过程中的力F如图11所示。图中，力F同样伴有剧烈震荡；若施加初始向上/下的干扰时，振动减小，与椭圆积分法计算得到的力结果不吻合。在力F的作用下，B端从初始位置开始向下运动，在梁屈曲前发生轴向压缩，力F迅速变大，当载荷大于屈曲临界值时，梁发生屈曲；此后力F则随位移的增大而减小，直至力F变为负值，梁发生了状态的跳转，拐点由两个变为了一个。图6和图8所示的两种变形形态得到的力F是相同的。需要指出的是：ANCF法中，梁屈曲前力F上升是由于梁发生了轴向压缩，而椭圆积分法中梁屈曲前力F上升段没有任何工程意义，只要给定足够小的位移，并且能够通过数值方法得到式(2)的解就能够使梁屈曲时的最大力F无限接近y轴，与实际情况不符。

由于整个梁产生了大范围的位移运动，因此为了方便观察梁截面产生的剪切变形，将梁截面在梁单元的单元局部坐标系中绘出，如图12给出了在梁AB上l=12 mm处梁截面的剪切变形情况。由图12可以看出在梁的变形过程中，截面产生了显著的剪切变形，截面尺寸拉伸和转角变化明显。

2.2 柔顺双稳态机构

为了进一步验证ANCF法的有效性，对如图13(a)所示柔顺双稳态复杂机构进行仿真。该机构包括一个滑块、两个柔性悬臂梁BO、B′O′和两个铰接-铰接型柔性屈曲梁AB、A′B′，BO和B′O′对称布置在滑块两侧。该机构具有两个稳态位置，滑块对应处于顶端或底端。如图13(b)，将相同拓扑特性的多组悬臂梁和铰接-铰接型屈曲梁串联，可实现滑块的精密微型直线运动。

图13(a)中，杆OB和BA长度均为1.25 m，两杆夹角为135°，截面为边长0.03 m的正方形，弹性模量E=2.07×1010Pa。为求得使该机构越过平衡位置(即杆BA水平时)施加的最大外力，同样在滑块上按式(13)加载位移载荷，直到机构到达第二个稳态位置结束。由于结构是对称的，只取其左半部分进行分析。图14所示给出了机构从第一个稳态位置开始运动到第二个稳态位置的变形过程。图15给出了机构运动过程中，施加的外力的变化规律。从图14、15可以看出，当t=1 s时，杆BA到达水平位置，施加的竖直力F变为0，当t=1.8 s时机构达到第二个平衡位置，施加的力F再次变为0。轨道对滑块的水平约束力如图16所示，可以看到水平约束力最大位置在杆BA越过平衡位置之后出现，并不是在杆BA水平时水平约束力最大。由变形可以进一步计算出O处的弯矩，如图17所示，其量值决定了机构的寿命。

t=0.06 s

t=0.08 s

t=0.10 s

t=0.12 s

t=0.14 s

t=0.16 s

t=0.18 s

t=0.20 s

(a) 与图6对应的变形(虚线为仿真结果，实线为实验结果)

t=0.06 s

t=0.08 s

t=0.10 s

t=0.12 s

t=0.14 s

t=0.18 s

t=0.20 s

t=0.22 s

(b) 与图8对应的变形(虚线为仿真结果，实线为实验结果)

图10 实验结果与仿真结果比较

Fig.10 Deformation comparison of the experimental results with that of the simulation

图11 力F和末端B位移的关系

图12 梁上l=12 mm处截面的剪切变形

(a)(b)

图13 柔顺双稳态机构

Fig.13 Compliant bistable mechanism

图14 柔顺双稳态机构的变形过程

经过以上两个算例的分析和比较，可以看到椭圆积分法基于欧拉伯努利梁理论，忽略了梁截面的剪切和轴向变形，在解支配方程时需要预先假定拐点的个数和梁末端的弯矩方向，对位置数的求解顺序也有一定的要求，计算用时少，但不确定性增加。而ANCF法可以计及梁的轴向变形、弯曲变形和截面的剪切变形，能更精确地描述大变形大位移大转动问题，通用性强，非常适合于柔顺机构建模与其动力学特性分析，采用多线程并行算法可以提高计算效率。

图15 竖直力F随时间t的变化规律

图16 水平力随时间t的变化规律

图17 点O处弯矩随时间t的变化规律

Fig.17 The relationship between the bending moment ofOand time

3 结 论

柔顺机构由于含非线性变形大运动构件，其建模变得困难。本文首先简述了基于Euler-Bernoulli梁理论的椭圆积分法和基于Timoshenko梁理论的绝对节点坐标法的基本原理，以固定—导向柔顺梁为例，仿真与实验结果证明了两种方法对简单柔顺机构建模与分析的有效性；其次采用ANCF法对柔顺双稳态机构进行建模与分析，验证了ANCF法对复杂柔顺机构动力学特性分析的有效性。通过对比两种方法的原理和仿真便利性，显示ANCF法对柔顺机构建模与分析更具有优越性，为柔顺机构这类现代机构学产品的工程设计与分析，提供了一种精确而高效的方法。