机械臂的刚柔耦合非线性动力学特性研究

空间RRRP机械臂的刚柔耦合非线性动力学特性研究王磊，王三民，牛治永(西北工业大学机电学院，西安710072)摘要机械臂的刚柔耦合非线性动力学特性分析是实现其运动精度控制的基础。针对某航天飞行器装载的RRRP型空间机械臂，采用模态组合函数描述各臂的弹性变形，将转动副角位移、移动幅线位移以及各臂的模态坐标作为广义自由度，利用Lagrange定理建立了空间RRRP机械臂的非线性动力学方程，采用4-5阶变步长Longgekuta法，对非线性微分方程组进行了数值求解。研究了机械臂刚柔耦合的非线性特性及结构参数变化对末端振动的影响，为进一步实现其结构优化和精度控制奠定基础。关键词：空间柔性臂，有限元，假设模态，拉格朗日方程，动力学中图分类号：TH113文献标识码：A0引言航天机械臂是航天飞行器上搭载机械手的系统，它可以用来带动机械手将卫星放入轨道或收回出故障的卫星，或将宇航员送至指定区域以完成各种复杂的空间作业。因此，现代飞行器大都配备有不同形式的机械臂[1]。这些机械臂通常具有质量轻、能耗低、机动性好等特点。然而，由于采用了轻质结构，其过大的变形对运动会造成很大影响。为了对其进行动态优化设计和运动精密控制，需要进行精确的动基金项目:陕西省自然科学基础研究计划项目(2005E228)第一作者王磊，硕士，1983年生力学分析[2]。近年来,考虑刚柔耦合的多柔体动力学研究已经成为机械动力学中最重要的研究方向之一[3]。然而，很多学者只是在平面系统内研究柔性机械臂的动力学特性，这就无法考虑系统在空间运动带来的向心力、科氏惯性力以及由此引起的陀螺力的影响，造成计算结果不能准确地应用到空间中去；另外，在已有柔性机械臂动力学分析中，基本上都采用了瞬时结构假设，将刚体运动与柔性变形分开考虑，因而忽略了刚柔耦合产生的非线性特性。本文对图1所示的某航天机械臂系统进行了刚柔耦合多柔体空间非线性动力学研究。图2为该机械臂的机构运动简图，它由3个转动副1R、2R、3R和一个移动幅组成。在进行系统的动力学分析时，不考虑飞行器姿态的调整对机械臂运动的影响。将飞行器基座视为固定点。由于RP1R2间的距离很短，所以将R1、R2、O三点视为重合。图1国际空间站上的机械臂Bm机械手集中质量P1111112将机械臂关节的驱动装置以及末端加持器和检测装置简化为集中质量。利用假设模态法和Lagrange方程建立该机械臂的动力学方程，采用Longgekuta法进行数值求解，得到系统广义坐标随时间的变化规律，为进一步实现柔性臂的结构优化和精度控制奠定基础。1动力学模型1.1坐标系和广义坐标的建立建立图3所示的总体坐标系，其原点O为机械臂的基座点，同时亦为旋转臂1的起点。整个机械臂绕Z轴旋转。局部坐标系OXYZΣ=1111oxyz1Σ=与旋转臂1固连，点与1o2R重合，1x轴为旋转臂1的理论轴线，其在OXY平面上的投影与X轴的夹角1θ（即绕Z轴的刚性转角）构成了第一个广义坐标。在垂直方向上，x1轴与平面OXY间图2机械臂机构运动简图1R3R2RElddρ旋转臂13333132Am关节集中质量BElddρ伸缩臂32222122Elddρ旋转臂2OA的夹角2θ（即旋转臂1在转矩1τ作用下的刚性转角）形成了第二个广义坐标。类似地可以建立旋转臂2的局部坐标系Σ=，将伸缩臂3与旋转臂2放在同一坐标系中，其在垂直方向上的转角2222oxyz23θ（即2x轴与水平面的夹角）为第三个广义坐标。伸缩臂3可以在旋转臂2内自由伸缩，并且提供了沿杆长方向上的广义坐标a。对柔性臂做如下假设和简化[1]：（1）材料为各向同性，其横向振动采用模态综合法表示为两个正交方向上的振动[2]11(,)()()(,)()()niyiyiiniziziiyxtqtxzxtqtxϕϕ==⎫=⎪⎪⎬⎪=⎪⎭∑∑(1)其中q、是与时间有关的模态广义坐标，()yit()ziqt()yixϕ、()zixϕ是选取的柔性臂关于局部坐标系iy轴和轴的模态基函数。iz（2）根据图2的结构，对旋转臂2和伸缩臂3的组合进行有限元模态分析。将伸缩臂3从完全收缩到完全伸出的过程分为9个状态，在marc中建立模型，并在原点处施加固定约束以定义悬臂梁，采用逆幂法分别求解其一阶固有振型，然后将振型曲线坐标值导入到MATLAB中。根据有限元分析得到的振型曲线形状，按照末端振型值为1的原则，设该系统的一阶固有振型函数为幂函数形式，即()2()/xxlaλϕ=+，式中λ的大小应与伸长量有关。按照最小二乘的原则逼近该曲线，采用黄金分割法寻找对应于每个模型的最优的aλ的值，得到伸出距离与近似基函数参数aλ的关系曲线。通过计算可以获得不同结构尺寸下的λ函数。表1给出了不同结构尺寸对应的λ函数。1232l5m5m6m21d80mm60mm80mm22d60mm40mm60mm3l4m4m4m31d60mm40mm60mm32d40mm0mm40mmZλ20.01670.02271.4952aa−+20.04020.07101.6537aa−+20.01270.02061.5095aa−+表1不同结构尺寸对应的λ函数而对于匀质悬臂梁，其一阶振型函数已知为coshcossinhsiniiiiixxxlllλλσλλ⎡⎤⎛⎞⎛⎞⎛⎞⎛⎞−−−⎜⎟⎜⎟⎜⎟⎜⎟⎢⎥⎝⎠⎝⎠⎝⎠⎝⎠⎣⎦xl，coshcossinhsiniiiiiλλσλλ+=+，同样可以按照上述方法拟合该函数曲线，得到其近似的幂函数形式。()1.5522()/xxlϕ=（3）在动力学建模过程中忽略轴向变形和剪切变形的影响，将每根柔性杆简化为Euler-Bernoulli梁。依据各柔性杆的支承与运动约束情况，并通过对一阶模态进行数值拟合，将旋转臂2和伸缩臂3的组合视为一根长度时变的非均匀截面臂，可得旋转臂1和组合臂的模态函数如下：1.552211()yxxlϕ⎛⎞=⎜⎟⎝⎠(2)11()sinzxxlπϕ=(3)222()()yzxxxlaλϕϕ⎛⎞==⎜+⎝⎠⎟(4)因此该系统的广义坐标列阵可以表示为()1231122Tyzyzqqqqaηθθθ=。而广义速度列阵可表示为()1231122Tyzyzqqqqaηθθθ=&&&&&&&&&。1.2系统动能图4为系统内任意点P的矢量表示。其中表示第i根臂上的任意点P在总体坐标系irΣ中的位置矢量，1iR−表示前一根臂的连接末端在总体坐标系Σ中的位置矢量。图3坐标系与广义坐标1τ2τ3τOAmBm()at1()tθ2()tθXY1x1y1z1()yt1()zt2()yt2()zt旋转臂1旋转臂2伸缩臂322z2yx3()tθ由图4的矢量关系知，在总体坐标系Σ中，旋转臂1上任意点位置矢量P101111()xyzrRAuuu=+++(5)式中()0000TR=，，，。()1100Txux=()11100Tyyyuqϕ=()1100Tzzuqϕ=1z矩阵为到的坐标转换矩阵，有1A1ΣΣ12112112112220ccscsAsccsssc−−⎡⎤⎢⎥=⎢⎢⎥⎣⎦−⎥(6)式中11coscθ=，22coscθ=，11sinsθ=，22sinsθ=。任意点P的速度矢量为()()()11111111122221111211222111221111211111111100yyyyyyyyyyyyTyyyyyyzzyyzzxxrAqAqqqxAAqqAqqxBCqqAqqϕϕϕϕθθϕθθϕϕϕθθϕϕϕϕ⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎛⎞∂∂⎢⎥=+⎜⎟⎢⎥∂∂⎝⎠⎢⎥⎣⎦+⎡⎤⎢⎥=+⎢⎥⎢⎥⎣⎦+&&&&&&&&&&&&&T1(7)式中12112111211000sccssA12Bccscsθ−−⎡⎤∂⎢⎥==−−⎢⎥∂⎢⎥⎣⎦(8)12121112222000cscsACsscsθ−12sc⎡⎤∂⎢⎥==−−⎢⎥∂⎢⎥−⎣⎦(9)令()11111111TyyzzDBxqqϕϕ=，()12111111TyyzzDCxqqϕϕ=，()11100TyyAAϕ=，，则()11100TzzAAϕ=1111122111yyzzrDDAqAqθθ=+++&&&&1&(10)因此旋转臂1的动能为11122111111011111122TlTyyzzTrrdxMqqqqθθθθρ⎛⎞⎛⎞⎜⎟⎜⎟⎜⎟⎜==⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠∫&&&&&&&&&&⎟(11)式中1M为旋转臂1的质量阵，表达式1111111111111121211211212110111111TTTTyzTTTlyzTTyyyzTzzDADADDDDDADADD1MdxAAAAAAρ⎛⎞⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟⎝⎠∫对称(12)将旋转臂2和伸缩臂3的组合视为一根长度和截面积时变的梁，在局部坐标系中，其两个方向上的一阶振型基函数均可表示为2Σ2()xxlaλϕ⎛⎞=⎜+⎝⎠⎟。按照旋转臂1求动能的方法，取旋转臂2的积分限为[]20,l,伸缩臂3的积分限为[]232,lalla+−+，分别推导旋转臂2、伸缩臂3、关节A以及机械手B的动能为：22222201122lTTTrrdx2Mρη==∫&&&η&(13)232233321122lalTTlaTrrdx3Mρη+−+==∫&&&η&(14)112111122TTAAAAAyyTmrrMqq2θθθθ⎛⎞⎛⎞⎜⎟⎜⎟==⎜⎟⎜⎜⎟⎜⎟⎝⎠⎝⎠&&&&&&&⎟&(15)图4系统内任意点的矢量表示ZixiyizXY1θ23()θθir1iR−yiyiqϕ⋅ziziqϕ⋅P1122331122221122TTBBBByyyzzTmrrMqqqqqaayqθθθθθθ⎛⎞⎛⎞⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟==⎜⎟⎜⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠&&&&&&&&&&&&&&&⎟&(16)式中，2M、3M、AM和BM为旋转臂2、伸缩臂3、关节A以及机械手B的质量阵。将各质量阵装配到总的质量阵M中，得到整个系统的总动能：12312TABTTTTTTMηη=++++=&&(17)1.3系统势能系统的弹性变形势能为1222322221111122011222222222220222222223322222121212lllalalyzVEIdxxxyzEIdxxxyzEIdxx++−⎡⎤⎛⎞⎛⎞∂∂⎢⎥=+⎜⎟⎜⎟∂∂⎢⎥⎝⎠⎝⎠⎣⎦⎡⎤⎛⎞⎛⎞∂∂⎢++⎜⎟⎜⎟∂∂⎢⎥⎝⎠⎝⎠⎣⎦⎡⎤⎛⎞⎛⎞∂∂⎢⎥++⎜⎟⎜⎟∂∂⎢⎥⎝⎠⎝⎠⎣⎦∫∫∫x⎥(18)式中1I、2I、3I分别为三根杆的截面惯性矩。将各模态函数代入上式并写成矩阵形式，可以得到用广义坐标表示的系统变形势能为14415511662277012012TyzzzyyTqkqkVqqkqqkKηη⎛⎞⎛⎞⎛⎜⎟⎜⎟⎜⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟⎜⎟⎜⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝⎠=11yzqq⎞⎟⎟(19)式中12211144111230117.037lyEIKEIdxxlϕ⎛⎞∂==⎜⎟⎜⎟∂⎝⎠∫(20)1224155111230112lz()()()()()()()()()()22232226622220222233222223222222333222233322123123123lylaylalKEIdxxEIdxxlEIlaEIlalalEIlaλλλλϕϕλλλλλλλλλ++−−3−⎛⎞∂=⎜⎟⎜⎟∂⎝⎠⎛⎞∂+⎜⎟⎜⎟∂⎝⎠−=−+−+−+−+−+−+∫∫(22)2223222772222022223326622lzlazlalKEIdxxEIdxxϕϕ++−⎛⎞∂=⎜⎟∂⎝⎠⎛⎞∂+=⎜⎟∂⎝⎠∫∫K(23)不计弹性变形对重力势能的影响，取XOY平面为零势能位置，则重力势能可以表示为()1121221223312233122331sinsin21sinsin21sinsin2sinsinsinABUmglmglmgllmgllalmgllaθθθθθθθθθ=++⎛⎞++⎜⎟⎝⎠⎡⎤⎛⎞++−+⎜⎟⎢⎥⎝⎠⎣⎦++(24)1.4动力学方程系统的拉格朗日方程为：idLLQdtηη⎛⎞