孫馳宇 沈惠平, 王一熙 許正驍 袁軍堂

(1.南京理工大學(xué)機(jī)械工程學(xué)院, 南京 210094； 2.常州大學(xué)現(xiàn)代機(jī)構(gòu)學(xué)研究中心, 常州 213016)

0 引言

機(jī)構(gòu)剛度是指機(jī)構(gòu)動平臺在外部載荷作用下，彈性構(gòu)件發(fā)生形變而產(chǎn)生位移的度量。剛度分析對并聯(lián)操作手實(shí)際應(yīng)用于參數(shù)選擇和結(jié)構(gòu)設(shè)計都具有重要的參考價值。

目前，主要的剛度分析方法有虛擬關(guān)節(jié)分析法[1-5]、有限元分析法[6-8]以及矩陣結(jié)構(gòu)分析法[9-10]。虛擬關(guān)節(jié)分析法將桿件設(shè)為剛體，在關(guān)節(jié)處通過建立柔性關(guān)節(jié)來描述桿件和關(guān)節(jié)所累積的柔性量；有限元分析法對桿件和關(guān)節(jié)的建模都需建立準(zhǔn)確的物理模型，因此精度高，但其計算量較大；矩陣結(jié)構(gòu)分析法將桿件和關(guān)節(jié)看作單元，相比于其他兩種方法，其分析效率高，但無法直接得到笛卡爾系的剛度矩陣。

基于以上3種剛度分析方法，國內(nèi)外學(xué)者進(jìn)行了諸多剛度建模實(shí)例分析。胡波等[11]給出了考慮約束反力產(chǎn)生變形的剛度建模方法，基于虛功原理建立機(jī)構(gòu)靜力學(xué)方程，通過分析機(jī)構(gòu)在合力作用下的變形量給出機(jī)構(gòu)剛度；徐東濤等[12]基于傳統(tǒng)的機(jī)構(gòu)剛度映射矩陣建立機(jī)構(gòu)剛度模型，并提出以機(jī)構(gòu)彈性變形評價機(jī)構(gòu)剛性特性的方法；曲海波等[13]通過鎖定機(jī)構(gòu)驅(qū)動副，運(yùn)用互易積運(yùn)算求解動平臺反螺旋力系，再以此推導(dǎo)出機(jī)構(gòu)剛度及動平臺的廣義位移；周玉林等[14]基于機(jī)構(gòu)各構(gòu)件的彈性變形，利用小變形疊加原理導(dǎo)出機(jī)構(gòu)靜力學(xué)方程，再運(yùn)用正交變換得到機(jī)構(gòu)靜剛度；楊超等[15]運(yùn)用螺旋理論和應(yīng)變能方法研究了具有2R1T三自由度的2UPR-RPU過約束并聯(lián)機(jī)構(gòu)的靜彈性剛度性能，模型考慮了桿件和關(guān)節(jié)的柔度；CECCARELLI等[16]研究了力的傳遞及力作用下的機(jī)構(gòu)變形，以此推導(dǎo)出機(jī)構(gòu)的剛度矩陣；MAJOU等[3]通過參數(shù)化剛度分析建立了正交矩陣的相容模型，在各向同性結(jié)構(gòu)中計算剛度矩陣元素，從而得到機(jī)構(gòu)剛度；PASHKEVICH等[17]提出利用6自由度虛擬彈簧建立剛度模型的方法；YAN等[18]基于Castigliano第2定理，利用應(yīng)變能法推導(dǎo)出一般平行四邊形機(jī)構(gòu)剛度的代數(shù)表達(dá)式，得到機(jī)構(gòu)的整體剛度矩陣。

虛擬彈簧法基于虛擬關(guān)節(jié)分析法，是一種通過在彈性連桿的末端增加虛擬彈簧來描述連桿的線性/旋轉(zhuǎn)變形以及變形之間的耦合特性的建模方法[19-20]，可對并聯(lián)機(jī)構(gòu)的支鏈進(jìn)行單獨(dú)建模，相比于上述建模方法，該方法不需要推導(dǎo)各類復(fù)雜的雅可比矩陣或運(yùn)動映射關(guān)系，僅需利用機(jī)構(gòu)的運(yùn)動學(xué)逆解即可計算機(jī)構(gòu)剛度(含處于奇異位形時)。

三平移并聯(lián)操作手在定位、抓取、裝卸等各種產(chǎn)業(yè)工藝操作上具有較廣泛的應(yīng)用。筆者團(tuán)隊提出一種新型零耦合度的非對稱三平移并聯(lián)操作手，該機(jī)構(gòu)可得到位置正解的解析解，同時，還具有運(yùn)動部分解耦特性，其運(yùn)動控制及軌跡規(guī)劃較易，已完成其運(yùn)動學(xué)分析和計算[21]。本文基于虛擬彈簧法對此機(jī)構(gòu)進(jìn)行剛度建模與特性分析，給出機(jī)構(gòu)整體剛度在工作空間中的分布規(guī)律，以及該機(jī)構(gòu)在不同方向、不同高度截面上的扭轉(zhuǎn)、線性剛度特性，以期為該操作手的動力學(xué)分析、樣機(jī)設(shè)計及其實(shí)驗(yàn)研究提供理論基礎(chǔ)。

1 機(jī)構(gòu)描述

(RPa∥3R)2R+RPa三平移并聯(lián)機(jī)構(gòu)的拓?fù)浣Y(jié)構(gòu)如圖1所示，它由動平臺1、靜平臺0及2條混合支鏈(HSOC)組成，坐標(biāo)系如圖1所示[22]。

圖1 零耦合度且運(yùn)動解耦的非對稱三平移并聯(lián)機(jī)構(gòu)Fig.1 Asymmetric three-translational parallel mechanism with zero coupling degree and motion decoupling

混合支鏈Ⅰ由支鏈A、B、C組成，其中4個球副(Sa、Sb、Sc、Sd)組成的平行四邊形結(jié)構(gòu)(Pa(4S))，支鏈A由驅(qū)動副R11串聯(lián)Pa(4S)組成；R11軸線與平行四邊形SaSb邊平行；支鏈B為驅(qū)動副R21串聯(lián)R22與R23組成，它們軸線相互平行；支鏈C由R12串聯(lián)R13組成，它們軸線相互平行。因此，混合支鏈Ⅰ記為：(RPa(4S)3R)⊥2R?；旌现ф湤?yàn)轵?qū)動副R31串聯(lián)R32后，又串聯(lián)由4個R副(Re、Rf、Rg、Rh)組成的平行四邊形結(jié)構(gòu)(Pa(4R))，進(jìn)一步再串聯(lián)R33組成，R31、R32、R33的軸線相互平行，因子串ReR32Rf、RgR33Rh分別等效于2個球副(S)，因此，該平行四邊形相當(dāng)于4個球副組成的平行四邊形，但不存在繞其對角線的轉(zhuǎn)動，故混合支鏈Ⅱ可記為：RPa。

該機(jī)構(gòu)自由度為3，當(dāng)取靜平臺0上的3個轉(zhuǎn)動副R11、R21、R31為驅(qū)動副時，動平臺1可實(shí)現(xiàn)沿x、y、z軸的三維平移，且沿y、z軸向的位移僅由驅(qū)動副R11、R21確定，因而具有輸入-輸出部分運(yùn)動解耦性[23]。

2 剛度模型建立

2.1 單桿剛度矩陣

建立剛度模型時，可將單桿看作懸臂梁，通過分析其末端變形求解剛度矩陣。在單桿受到力/力矩作用時，根據(jù)歐拉-伯努利梁理論，可得連桿的撓曲線方程，在彎曲變形很小且材料服從胡克定律的情況下，撓曲線方程是線性的，因此，考慮到連桿受力/力矩的耦合情況，采用疊加法計算連桿在力/力矩作用下的柔度矩陣，再對柔度矩陣求逆，即可得到單桿剛度矩陣為[20，25]

(1)

式中G——楊氏模量

Ix、Iy、Iz——截面關(guān)于x、y、z軸的慣性矩

l——單桿桿長

E——彈性模量

A——桿件截面積

式(1)是計算機(jī)構(gòu)各支鏈及整體剛度矩陣的基礎(chǔ)。為便于后續(xù)表述和計算，將式中各元素的表達(dá)式分別記作Kij。

2.2 4S平行四邊形結(jié)構(gòu)剛度矩陣

平行四邊形結(jié)構(gòu)無法直接采用單桿的剛度矩陣進(jìn)行計算，需要將其作為獨(dú)立結(jié)構(gòu)進(jìn)行剛度建模。在4S平行四邊形結(jié)構(gòu)中，建立圖2所示的局部坐標(biāo)系，其中，平行四邊形的兩短桿可視為剛性構(gòu)件，兩長桿可視為懸臂梁。o1點(diǎn)為四邊形結(jié)構(gòu)末端基點(diǎn)；x1軸方向與兩短桿中點(diǎn)連線重合，指向o1點(diǎn)；y1軸方向垂直于四邊形平面向外，z1軸方向由右手法則確定。

圖2 4S平行四邊形結(jié)構(gòu)的剛度模型Fig.2 Stiffness model of 4S parallelogram structures

4S四邊形結(jié)構(gòu)的剛度矩陣為

(2)

其中，Sq表示sinq，下同。

2.3 4R平行四邊形結(jié)構(gòu)剛度矩陣

對于4R平行四邊形結(jié)構(gòu)，局部坐標(biāo)系建立方法同圖2，如圖3所示，平行四邊形的兩短桿可視為剛性構(gòu)件，兩長桿可視為懸臂梁。

圖3 4R平行四邊形結(jié)構(gòu)的剛度模型Fig.3 Stiffness model of 4R parallelogram structures

因此，4R四邊形結(jié)構(gòu)的剛度矩陣為

(3)

其中，Cq表示cosq。

因末端點(diǎn)o2沿z2方向的力受到被動副Re、Rf與被動副Rh、Rg的補(bǔ)償而消失，因此，該4R四邊形結(jié)構(gòu)的末端可用5-DOF的虛擬彈簧來描述其剛度特性。

2.4 一般支鏈的剛度建模

2.4.1運(yùn)動方程的建立

一般支鏈通常由驅(qū)動副(Ac)、主動桿、從動桿以及被動副組成，其虛擬彈簧模型如圖4所示。

圖4 4R一般支鏈的虛擬彈簧模型Fig.4 Virtual spring model of 4R branch chain

1-DOF的虛擬彈簧，表示驅(qū)動副Ac的伺服剛度，其變形可表示為Δθ0；6-DOF的虛擬彈簧，表示桿件在笛卡爾坐標(biāo)系中三自由度的旋轉(zhuǎn)變形和三自由度的拉伸變形，主動桿和被動桿上彈簧的變形量可分別表示為(Δθ1,Δθ2,…,Δθ6)和(Δθ7,Δθ8,…,Δθ12)。

基于圖4給出支鏈中彈簧變形和被動關(guān)節(jié)變形到末端變形之間的一般性運(yùn)動方程為

(4)

其中

式中 Δt——靜笛卡爾坐標(biāo)系中機(jī)構(gòu)末端變形

Δφ——旋轉(zhuǎn)變形 Δp——拉伸變形

將式(4)的支鏈運(yùn)動方程表示成螺旋運(yùn)動形式，即

(5)

其中

式中 $——機(jī)構(gòu)末端參考點(diǎn)變形的旋量

Δθi——虛擬彈簧變形量

Δψi——被動副運(yùn)動位移

2.4.2靜力學(xué)方程的建立

(6)

因此，外力fi所做虛功為

(7)

(8)

而被動關(guān)節(jié)受力后會發(fā)生被動運(yùn)動，所以在靜平衡狀態(tài)下，被動運(yùn)動不做功，即

(9)

由式(8)、(9)可得

(10)

(11)

于是，支鏈的靜力平衡方程可表示為

(12)

其中

令支鏈i的笛卡爾剛度矩陣為Ki，則有

fi=KiΔt

(13)

式中fi——支鏈i所受的外力

由式(12)、(13)可得

(14)

由于每條支鏈所受外力為機(jī)構(gòu)整體所受外力的分量，因此，一般整體的剛度矩陣為

(15)

2.5 機(jī)構(gòu)的剛度建模

該三平移并聯(lián)機(jī)構(gòu)由2條HSOC構(gòu)成，其支鏈拓?fù)浣Y(jié)構(gòu)簡化后如圖5所示。

圖5 機(jī)構(gòu)的支鏈拓?fù)浣Y(jié)構(gòu)簡圖Fig.5 Schematics of branch chain

2.5.1HSOCⅠ的剛度建模

HSOCⅠ由支鏈A、B并聯(lián)組成子并聯(lián)機(jī)構(gòu)后再串聯(lián)支鏈C組成，因此，先分別求出支鏈A、B的剛度矩陣，即可求得子并聯(lián)機(jī)構(gòu)的剛度矩陣；再將它看成一個整體與支鏈C串聯(lián)進(jìn)行剛度建模，即可得到HSOCⅠ的剛度矩陣。

(1)支鏈A的剛度矩陣

支鏈A的剛度模型如圖6所示，其中，驅(qū)動副(Ac)對應(yīng)笛卡爾坐標(biāo)系中1個旋轉(zhuǎn)自由度的虛擬彈簧變形；主動桿桿2對應(yīng)笛卡爾坐標(biāo)系中6個自由度的虛擬彈簧變形，包含3個旋轉(zhuǎn)變形和3個線性變形；4S平行四邊形結(jié)構(gòu)對應(yīng)笛卡爾坐標(biāo)系中2個自由度的虛擬彈簧(2.2節(jié))。在支鏈A中，所有被動副均在4S平行四邊形結(jié)構(gòu)內(nèi)，可將其看成整體結(jié)構(gòu)，因此，可認(rèn)為支鏈A不包含被動副。

圖6 支鏈A的虛擬彈簧模型Fig.6 Virtual spring model of branch chain A

由式(12)、(14)可得支鏈A的靜力方程為

其中

式中KA——支鏈A的剛度矩陣

KAc——驅(qū)動剛度

(2)支鏈B的剛度矩陣

支鏈B的剛度模型如圖7所示，其中，驅(qū)動副(Ac)對應(yīng)笛卡爾坐標(biāo)系中1個旋轉(zhuǎn)自由度的虛擬彈簧變形；主動桿桿6和從動桿桿7各對應(yīng)笛卡爾坐標(biāo)系中6個自由度的虛擬彈簧變形。

圖7 支鏈B的虛擬彈簧建模型Fig.7 Virtual spring model of branch chain B

由式(12)、(14)可得支鏈B的靜力方程為

其中

eB=[0 -1 0]T

式中KB——支鏈B剛度矩陣

在求得支鏈A和B的剛度矩陣后，由式(15)可得子并聯(lián)結(jié)構(gòu)的整體剛度矩陣為

Ksub=KA+KB

(3)混合支鏈的剛度矩陣

將支鏈A、B組成的子并聯(lián)機(jī)構(gòu)看成一個整體與支鏈C串聯(lián)進(jìn)行剛度建模，混合支鏈的剛度模型如圖8所示，其中，從動桿桿8對應(yīng)笛卡爾坐標(biāo)系中6個自由度的虛擬彈簧變形；桿C1E1、F1G1因桿長遠(yuǎn)短于其他桿件，因此，將它們看作剛性構(gòu)件，不考慮在其末端建立虛擬彈簧。

圖8 HSOCⅠ的虛擬彈簧模型Fig.8 Virtual spring modeling of HSOC Ⅰ

由式(12)、(14)可得HSOCⅠ的靜力方程為

其中

e1=[1 0 0]Te2=[0 1 0]T

e3=[0 0 1]T

式中K1——混合支鏈Ⅰ的剛度矩陣

2.5.2HSOCⅡ的剛度建模

HSOCⅡ的剛度模型如圖9所示，其中，驅(qū)動副(Ac)對應(yīng)笛卡爾坐標(biāo)系中1個旋轉(zhuǎn)自由度的虛擬彈簧變形；主動桿桿6對應(yīng)笛卡爾坐標(biāo)系中6個自由度的虛擬彈簧變形；4R平行四邊形結(jié)構(gòu)對應(yīng)笛卡爾坐標(biāo)系中5個自由度的虛擬彈簧(由2.3節(jié)說明)。被動副R32、R33各存在1個自由度的微小變形。

圖9 HSOCⅡ的虛擬彈簧建模Fig.9 Virtual spring modeling of HSOCⅡ

由式(12)、(14)可得HSOCⅡ的靜力方程為

其中

式中K2——混合支鏈Ⅱ的剛度矩陣

由式(15)可得機(jī)構(gòu)的整體剛度為

K=K1+K2

3 剛度分析

3.1 機(jī)構(gòu)參數(shù)

機(jī)構(gòu)各桿件的尺寸參數(shù)如表1所示。為了減輕機(jī)構(gòu)質(zhì)量，提高機(jī)構(gòu)強(qiáng)度，機(jī)構(gòu)桿件全部選擇碳纖維材料。機(jī)構(gòu)3個驅(qū)動輸入采用相同的驅(qū)動電機(jī)，其驅(qū)動剛度取為5×104N·m/rad。

表1 機(jī)構(gòu)桿件尺寸參數(shù)Tab.1 Dimension parameters of mechanism

3.2 數(shù)值算例

根據(jù)第2節(jié)中的剛度矩陣建模，可計算出機(jī)構(gòu)在某一姿態(tài)下的整體靜剛度矩陣?，F(xiàn)取動平臺基點(diǎn)p的坐標(biāo)為(0.02,0.03,0.8)，可得該姿態(tài)下機(jī)構(gòu)的整體剛度為

(16)

可知剛度矩陣K為6×6方陣，其中，對角線前3項(xiàng)為機(jī)構(gòu)在x、y、z軸方向的扭轉(zhuǎn)剛度(單位：N·m/rad)；后3項(xiàng)為機(jī)構(gòu)在x、y、z軸方向的線性剛度(單位：N/m)。

機(jī)構(gòu)在運(yùn)動過程中，剛度隨著機(jī)構(gòu)末端基點(diǎn)位置的變化而變化，為研究其分布規(guī)律，定義剛度矩陣中對角線上6項(xiàng)數(shù)值的平均值η[24]，即

(17)

利用Matlab軟件計算η在機(jī)構(gòu)工作空間中的分布，如圖10所示。

由圖10可知，當(dāng)動平臺在y軸方向上越靠近兩側(cè)時，機(jī)構(gòu)的整體剛度越大；為更加清晰地分析z軸方向?qū)C(jī)構(gòu)整體剛度的影響，利用Matlab軟件計算機(jī)構(gòu)工作空間中z=0.75 m和z=0.90 m兩截面的η分布，如圖11所示。

圖10 工作空間中機(jī)構(gòu)的剛度分布Fig.10 Stiffness distribution of mechanism

圖11 不同z值截面的機(jī)構(gòu)剛度分布Fig.11 Stiffness distribution of mechanism with different z sections

由圖11可知，動平臺在不同高度下，剛度變化趨勢相似；但隨著z的減小，機(jī)構(gòu)整體剛度變大。

為了更加具體地分析剛度矩陣中各方向的扭轉(zhuǎn)剛度和線性剛度的分布情況，取機(jī)構(gòu)工作空間中z=0.75 m和z=0.90 m兩截面進(jìn)行分析，如圖12所示。

圖12 不同z值截面的剛度矩陣主對角線元素分布情況Fig.12 Distributions of main diagonal elements of stiffness matrix with different z-sections

由圖12可知，機(jī)構(gòu)x軸方向的扭轉(zhuǎn)剛度最大，y、z軸方向的扭轉(zhuǎn)剛度相近，原因?yàn)橹ф淏不能繞x軸方向旋轉(zhuǎn)，因此，x軸方向的扭轉(zhuǎn)剛度高于其他兩軸方向。z軸方向的線性剛度最大，原因?yàn)樵擃悪C(jī)械手在工作時，其所受的外部線性載荷主要由z軸方向承擔(dān)，因此，需要z軸方向有較大的線性剛度。剛度特性均隨著z的減小而增大，說明機(jī)構(gòu)動平臺離定平臺越近，剛度越大，因此，只要機(jī)械手在剛抓取物件時不發(fā)生變形，則在其向上提升物件的運(yùn)送過程中，機(jī)構(gòu)也不會變形，這一特性符合機(jī)械臂的抓取性能需要。

3.3 數(shù)值驗(yàn)證

為了驗(yàn)證數(shù)值算例的正確性，現(xiàn)將該機(jī)構(gòu)的簡化模型導(dǎo)入ANSYS Workbench，進(jìn)行有限元分析(FEA)。當(dāng)機(jī)構(gòu)末端基點(diǎn)p的位置坐標(biāo)為(0.02,0.03,0.8) m時，在動平臺上施加單位力、力矩，將動平臺看作剛性無限大的柔性構(gòu)件，對機(jī)構(gòu)進(jìn)行剛?cè)狁詈戏治?，網(wǎng)格劃分如圖13a所示；考慮約束和受力，如圖13b所示。由此得到動平臺相對于靜坐標(biāo)系各個方向產(chǎn)生的微小位姿變化，其值分別對應(yīng)虛擬彈簧法(VSM)中所得剛度矩陣的逆矩陣中對角線上6個元素的數(shù)值。

對式(16)求逆，可得機(jī)構(gòu)的柔度矩陣為

圖13 機(jī)構(gòu)的剛?cè)狁詈戏治鯢ig.13 Rigid-flexible coupling analysis of mechanism

將有限元分析所得的機(jī)構(gòu)變形結(jié)果與虛擬彈簧法所得的機(jī)構(gòu)變形結(jié)果進(jìn)行對比，并得到其相對誤差，如表2所示。

由表2可知，理論變形與仿真變形的誤差在-6%內(nèi)。因此，虛擬彈簧法的建模精度滿足工程設(shè)計要求。

圖14 兩機(jī)構(gòu)的η值分布Fig.14 η distributions of two mechanisms

4 (RPa∥3R)2R+RPa與Delta機(jī)構(gòu)的剛度比較

將三平移(RPa∥3R)2R+RPa機(jī)構(gòu)與Delta機(jī)構(gòu)進(jìn)行剛度比較，兩者選取相同的機(jī)構(gòu)參數(shù)(3.1節(jié))，取工作空間中z1=0.85 m、z2=0.80 m、z3=0.75 m處的3個截面進(jìn)行研究。

根據(jù)文獻(xiàn)[25]中Delta機(jī)構(gòu)的運(yùn)動學(xué)反解，使用Matlab編程，得到3個截面中兩機(jī)構(gòu)的剛度參數(shù)η在x∈[-0.04,0.04] m，y∈[-0.04,0.04] m范圍內(nèi)的分布，及其在xoy平面的投影，如圖14所示。

由圖14可知，(RPa∥3R)2R+RPa機(jī)構(gòu)與Delta機(jī)構(gòu)，在y軸方向上越靠近兩側(cè)時，機(jī)構(gòu)的整體剛度越大。Delta機(jī)構(gòu)的剛度分布曲面更為規(guī)則且剛度變化更為平穩(wěn)。(RPa∥3R)2R+RPa機(jī)構(gòu)的剛度大于Delta機(jī)構(gòu)的部分在z1、z2、z3截面中所占的面積(圖中綠色部分)分別為60.54%、58.69%、56.37%(由劃分網(wǎng)格的面積估算)，均大于50%，因此，可得(RPa∥3R)2R+RPa機(jī)構(gòu)的剛度較大于Delta機(jī)構(gòu)。

5 結(jié)論

(1)運(yùn)用虛擬彈簧法對(RPa∥3R)2R+RPa三平移機(jī)構(gòu)進(jìn)行剛度建模，得到了該機(jī)構(gòu)笛卡爾空間的剛度矩陣。

(2)給出了該機(jī)構(gòu)在工作空間中的整體剛度分布，并分別對x、y、z軸方向的扭轉(zhuǎn)、線性剛度進(jìn)行分析。結(jié)果表明：機(jī)構(gòu)x軸方向的扭轉(zhuǎn)剛度最大，y、z軸方向的扭轉(zhuǎn)剛度相近，z軸方向的線性剛度最大；機(jī)構(gòu)動平臺與定平臺距離越近，剛度性能越穩(wěn)定。

(3)對機(jī)構(gòu)變形進(jìn)行有限元分析，并與虛擬彈簧法所得機(jī)構(gòu)變形結(jié)果進(jìn)行對比，結(jié)果表明，理論變形與仿真變形的誤差在-6%內(nèi)，驗(yàn)證了虛擬彈簧法所得剛度結(jié)果的正確性。

(4)(RPa∥3R)2R+RPa和Delta機(jī)構(gòu)的剛度特性對比表明，Delta機(jī)構(gòu)的剛度分布曲面更為規(guī)則，且剛度變化更為平穩(wěn)，(RPa∥3R)2R+RPa的剛度大于Delta機(jī)構(gòu)。