零耦合度部分運動解耦三平移并聯機構剛度建模與分析

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

(1.南京理工大學機械工程學院, 南京 210094； 2.常州大學現代機構學研究中心, 常州 213016)

0 引言

機構剛度是指機構動平臺在外部載荷作用下，彈性構件發生形變而產生位移的度量。剛度分析對并聯操作手實際應用于參數選擇和結構設計都具有重要的參考價值。

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

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

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

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

1 機構描述

(RPa∥3R)2R+RPa三平移并聯機構的拓撲結構如圖1所示，它由動平臺1、靜平臺0及2條混合支鏈(HSOC)組成，坐標系如圖1所示[22]。

圖1 零耦合度且運動解耦的非對稱三平移并聯機構Fig.1 Asymmetric three-translational parallel mechanism with zero coupling degree and motion decoupling

混合支鏈Ⅰ由支鏈A、B、C組成，其中4個球副(Sa、Sb、Sc、Sd)組成的平行四邊形結構(Pa(4S))，支鏈A由驅動副R11串聯Pa(4S)組成；R11軸線與平行四邊形SaSb邊平行；支鏈B為驅動副R21串聯R22與R23組成，它們軸線相互平行；支鏈C由R12串聯R13組成，它們軸線相互平行。因此，混合支鏈Ⅰ記為：(RPa(4S)3R)⊥2R。混合支鏈Ⅱ為驅動副R31串聯R32后，又串聯由4個R副(Re、Rf、Rg、Rh)組成的平行四邊形結構(Pa(4R))，進一步再串聯R33組成，R31、R32、R33的軸線相互平行，因子串ReR32Rf、RgR33Rh分別等效于2個球副(S)，因此，該平行四邊形相當于4個球副組成的平行四邊形，但不存在繞其對角線的轉動，故混合支鏈Ⅱ可記為：RPa。

該機構自由度為3，當取靜平臺0上的3個轉動副R11、R21、R31為驅動副時，動平臺1可實現沿x、y、z軸的三維平移，且沿y、z軸向的位移僅由驅動副R11、R21確定，因而具有輸入-輸出部分運動解耦性[23]。

2 剛度模型建立

2.1 單桿剛度矩陣

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

(1)

式中G——楊氏模量

Ix、Iy、Iz——截面關于x、y、z軸的慣性矩

l——單桿桿長

E——彈性模量

A——桿件截面積

式(1)是計算機構各支鏈及整體剛度矩陣的基礎。為便于后續表述和計算，將式中各元素的表達式分別記作Kij。

2.2 4S平行四邊形結構剛度矩陣

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

圖2 4S平行四邊形結構的剛度模型Fig.2 Stiffness model of 4S parallelogram structures

4S四邊形結構的剛度矩陣為

(2)

其中，Sq表示sinq，下同。

2.3 4R平行四邊形結構剛度矩陣

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

圖3 4R平行四邊形結構的剛度模型Fig.3 Stiffness model of 4R parallelogram structures

因此，4R四邊形結構的剛度矩陣為

(3)

其中，Cq表示cosq。

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

2.4 一般支鏈的剛度建模

2.4.1運動方程的建立

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

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

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

基于圖4給出支鏈中彈簧變形和被動關節變形到末端變形之間的一般性運動方程為

(4)

其中

式中 Δt——靜笛卡爾坐標系中機構末端變形

Δφ——旋轉變形 Δp——拉伸變形

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

(5)

其中

式中 $——機構末端參考點變形的旋量

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

Δψi——被動副運動位移

2.4.2靜力學方程的建立

(6)

因此，外力fi所做虛功為

(7)

(8)

而被動關節受力后會發生被動運動，所以在靜平衡狀態下，被動運動不做功，即

(9)

由式(8)、(9)可得

(10)

(11)

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

(12)

其中

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

fi=KiΔt

(13)

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

由式(12)、(13)可得

(14)

由于每條支鏈所受外力為機構整體所受外力的分量，因此，一般整體的剛度矩陣為

(15)

2.5 機構的剛度建模

該三平移并聯機構由2條HSOC構成，其支鏈拓撲結構簡化后如圖5所示。

圖5 機構的支鏈拓撲結構簡圖Fig.5 Schematics of branch chain

2.5.1HSOCⅠ的剛度建模

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

(1)支鏈A的剛度矩陣

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

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

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

其中

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

KAc——驅動剛度

(2)支鏈B的剛度矩陣

支鏈B的剛度模型如圖7所示，其中，驅動副(Ac)對應笛卡爾坐標系中1個旋轉自由度的虛擬彈簧變形；主動桿桿6和從動桿桿7各對應笛卡爾坐標系中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)可得子并聯結構的整體剛度矩陣為

Ksub=KA+KB

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

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

圖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所示，其中，驅動副(Ac)對應笛卡爾坐標系中1個旋轉自由度的虛擬彈簧變形；主動桿桿6對應笛卡爾坐標系中6個自由度的虛擬彈簧變形；4R平行四邊形結構對應笛卡爾坐標系中5個自由度的虛擬彈簧(由2.3節說明)。被動副R32、R33各存在1個自由度的微小變形。

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

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

其中

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

由式(15)可得機構的整體剛度為

K=K1+K2

3 剛度分析

3.1 機構參數

機構各桿件的尺寸參數如表1所示。為了減輕機構質量，提高機構強度，機構桿件全部選擇碳纖維材料。機構3個驅動輸入采用相同的驅動電機，其驅動剛度取為5×104N·m/rad。

表1 機構桿件尺寸參數Tab.1 Dimension parameters of mechanism

3.2 數值算例

根據第2節中的剛度矩陣建模，可計算出機構在某一姿態下的整體靜剛度矩陣。現取動平臺基點p的坐標為(0.02,0.03,0.8)，可得該姿態下機構的整體剛度為

(16)

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

機構在運動過程中，剛度隨著機構末端基點位置的變化而變化，為研究其分布規律，定義剛度矩陣中對角線上6項數值的平均值η[24]，即

(17)

利用Matlab軟件計算η在機構工作空間中的分布，如圖10所示。

由圖10可知，當動平臺在y軸方向上越靠近兩側時，機構的整體剛度越大；為更加清晰地分析z軸方向對機構整體剛度的影響，利用Matlab軟件計算機構工作空間中z=0.75 m和z=0.90 m兩截面的η分布，如圖11所示。

圖10 工作空間中機構的剛度分布Fig.10 Stiffness distribution of mechanism

圖11 不同z值截面的機構剛度分布Fig.11 Stiffness distribution of mechanism with different z sections

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

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

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

由圖12可知，機構x軸方向的扭轉剛度最大，y、z軸方向的扭轉剛度相近，原因為支鏈B不能繞x軸方向旋轉，因此，x軸方向的扭轉剛度高于其他兩軸方向。z軸方向的線性剛度最大，原因為該類機械手在工作時，其所受的外部線性載荷主要由z軸方向承擔，因此，需要z軸方向有較大的線性剛度。剛度特性均隨著z的減小而增大，說明機構動平臺離定平臺越近，剛度越大，因此，只要機械手在剛抓取物件時不發生變形，則在其向上提升物件的運送過程中，機構也不會變形，這一特性符合機械臂的抓取性能需要。

3.3 數值驗證

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

對式(16)求逆，可得機構的柔度矩陣為

圖13 機構的剛柔耦合分析Fig.13 Rigid-flexible coupling analysis of mechanism

將有限元分析所得的機構變形結果與虛擬彈簧法所得的機構變形結果進行對比，并得到其相對誤差，如表2所示。

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

圖14 兩機構的η值分布Fig.14 η distributions of two mechanisms

4 (RPa∥3R)2R+RPa與Delta機構的剛度比較

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

根據文獻[25]中Delta機構的運動學反解，使用Matlab編程，得到3個截面中兩機構的剛度參數η在x∈[-0.04,0.04] m，y∈[-0.04,0.04] m范圍內的分布，及其在xoy平面的投影，如圖14所示。

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

5 結論

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

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

(3)對機構變形進行有限元分析，并與虛擬彈簧法所得機構變形結果進行對比，結果表明，理論變形與仿真變形的誤差在-6%內，驗證了虛擬彈簧法所得剛度結果的正確性。

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