球面兩自由度冗余驅動并聯機器人彈性動力學分析

李永泉 佘亞中 萬一心 張立杰

1.燕山大學河北省并聯機器人與機電系統實驗室，秦皇島，066004 2.燕山大學先進鍛壓成形技術與科學教育部重點實驗室，秦皇島，066004 3.燕山大學河北省重型機械流體動力傳輸與控制實驗室，秦皇島，066004

0 引言

隨著并聯機器人的發展，高速、輕巧、重載成為并聯機器人的發展方向。當輕質化的機器人在重載、高速的環境下工作時，由承受重載引起的各分支的彈性變形會影響機器人的運動學及動力學性能，并降低機器人的精度，因此，有必要對并聯機器人進行彈性動力學分析。一些學者已經對并聯機構的彈性動力學做了大量的研究工作。文獻［1］比較分析了6-PSS與8-PSS冗余驅動并聯平臺的彈性動力學特性。文獻［2-3］考慮連桿的變形，建立了空間3-RRS并聯機構的彈性動力學模型，并分析了其動力學特性。文獻［4］提出了一種簡單而精確的彈性動力學建模方法，建立了柔性并聯機器人動力學建模。文獻［5］提出了一種3-RPS并聯機構彈性動力學建模方法，并通過模態實驗驗證了該方法所得的固有頻率和振型。文獻［6］以3-PRS并聯機構為例，通過建立其彈性動力學模型，提出了分析其關節約束反力的方法。文獻［7］建立了5-UPS/PRPU冗余驅動并聯機床的彈性動力學模型，發現冗余驅動可以減小因彈性變形引起的誤差，但該種冗余驅動屬于關節式冗余驅動。上述并聯機構彈性動力學分析中均未涉及空間彎桿的彎曲、拉伸和扭轉變形問題，同時關于分支冗余驅動并聯機構的彈性動力學分析也鮮有報道。

1 彈性動力學建模

如圖1所示，本文研究的球面2-DOF冗余驅動并聯機構是在球面5R并聯機構的基礎上，通過添加一個驅動分支變成的冗余驅動并聯機構，它只包含轉動副，并且所有轉動副的軸線匯交于球心O。機構的輸出點為P(x,y,z)，它通過3個完全相同的3R串聯分支 PB1A1、PB2A2、PB3A3與正三角形定平臺A1A2A3相連接。驅動連桿的圓弧半徑為ri1(i=1,2,3)，圓弧角為α1，且驅動軸線與定平臺中垂線的夾角為 β。分支上連桿的圓弧半徑為ri2(i=1,2,3)，圓弧角為α2。實驗樣機如圖2所示，其結構參數為 α1=α2=90°，β=0°。

圖1 機構簡圖Fig.1 Mechanism sketch

在保證分析精度的情況下，為了便于分析計算，本文對模型進行以下簡化［8］：①機構的定平臺較其他構件剛度大，故視為剛體；②輸出軸OP假定為剛體；③機構在運動過程中，忽略剛體運動和彈性變形運動的耦合項，把真實運動視為剛體運動和彈性變形運動的疊加；④柔性梁單元的位移滿足三次多項式和線性分布；⑤所有的轉動副均為理想約束，即不考慮轉動副的變形與間隙。

圖2 實驗樣機Fig.2 Experimental prototype

1.1 單元動力學方程

將球面2-DOF冗余驅動并聯機器人中的各構件 AiBi和 BiP(i=1,2,3)當作柔性桿處理，輸出軸OP當作剛體處理。由于并聯機構中的系統支鏈一般較短，關節柔性變形引起的系統累計誤差較小，而且并聯機構為多閉環系統，對系統中的關節變形具有一定的約束作用，故關節柔性一般可以忽略不計。將機構的連桿AiBi和BiP分別設為單元i1和i2，整個機器人系統共包含6個柔性單元。

由于6個柔性單元模型都是結構一樣的空間曲梁單元，故本文對其中一個柔性單元BiP進行分析。如圖3所示，這里假定空間柔性曲梁單元發生軸向、橫向、扭轉變形，并用 δ =(δ1,δ2,…,δ18)Τ表示梁單元的廣義坐標向量，其中，δ1～ δ3與δ10～ δ12、δ4～δ6與 δ13～δ15、δ7～δ9與 δ16～δ18分別表示節點A、B處的彈性位移、彈性轉角和曲率。如圖4所示，依據有限元法，分別求出空間圓弧曲梁結構的柔性單元BiP的質量矩陣與剛度矩陣：

式中，N1、N2、N3、N4均為差值函數組成的向量，且向量中的各元素均為θ的函數。

將單元的動能與勢能代入拉格朗日方程，可以得到單元動力學方程：

就球面2-DOF并聯機器人系統而言，Fe是單元外加載荷的廣義力列陣，既包括與廣義坐標對應的集中節點力和力矩，也包括分布力和力矩的等效節點力，所有的力和力矩都是真實作用的外載荷；Qe是單元剛體慣性力列陣。

圖3 空間柔性梁單元廣義坐標圖Fig.3 Generalized coordinates of flexible space beam element

圖4 上連桿有限元圖Fig.4 Finite element of the connecting rod

1.2 支鏈動力學方程

圖5 所示為球面2-DOF冗余驅動并聯機器人的一個分支 AiBiP，對于驅動構件 AiBi(i=1,2,3)，可以看作懸臂梁，引入邊界條件可以得到Ai點處的彈性位移和轉角方向的節點變形均為零。Bi點處為轉動副，兩個相鄰的單元（單元i1和i2）分屬不同的構件，因而具有不同的轉角，在這兩個結點處，其中，繞Bi點處轉動副軸線方向的兩個曲率為零。構件BiP的節點P與剛性輸出軸連接，所以P點處轉動副軸線方向的兩個曲率為零。

圖5 單元坐標系下的廣義坐標Fig.5 Generalized coordinates in element coordinate system

如圖6所示，支鏈下連桿AiBi有11個非零廣義坐標，支鏈上連桿BiP有16個非零廣義坐標，分支AiBiP一共有27個廣義非零坐標。由于下連桿AiBi與上連桿BiP之間通過轉動副連接，轉動副的轉軸看作剛性軸連接，可知節點Bi處的3個彈性位移相同。此時，分支AiBiP的節點變形可用24個系統坐標表示為Ui=[ui1ui2… ui24]Τ，如圖5所示。由單元AiBi與BiP在單元坐標系下的動力學方程（式（3）），經過裝配可得分支AiBiP在定坐標系下的動力學方程：

式中，Ui為第i個支鏈 AiBiP上的單元廣義坐標，Ui∈R24×1；Mi為支鏈 AiBiP 的質量矩陣，Mi∈R24×24；Ki為支鏈 AiBiP 的剛度矩陣，Ki∈R24×24；Fi為支鏈AiBiP 的外加載荷廣義力矩陣，Fi∈R24×1；Qi為支鏈AiBiP 的剛體慣性力矩陣，Qi∈R24×1。

圖6 定坐標系下的廣義坐標Fig.6 Generalized coordinates in fixed coordinate system

1.3 系統彈性動力學模型

由于各支鏈桿件的變形引起的輸出軸的位置變化量U0=[uxuyuz]Τ，3個分支的P點與剛性輸出軸連接，根據3個分支與輸出軸的變形協調約束，有

式中，UPi為第i個分支P點的彈性位移；為第i個分支P點處的節點坐標系與定坐標系之間的變換矩陣。

取單個支鏈的廣義坐標為

球面2-DOF冗余驅動并聯機器人系統的運動約束方程為

其中，Ei×j為 i×j單位矩陣，0i×j為 i×j零矩陣。將式（7）代入式（4），可得

對式兩邊同乘可得

化簡整理式（9）可得

將式（10）中的分解為

采用系統廣義坐標U，并利用系統運動約束方程式，把廣義坐標表示的各個支鏈的動力學方程式裝配到一起，形成球面2-DOF并聯機器人系統的無阻尼彈性動力學方程：

式中，M 為系統總質量矩陣，M∈R66×66；K為系統總剛度 矩陣，K∈R66×66；U¨r為 系統剛體加速 度矩陣，U¨r∈R66×1；U 為系統廣義坐標矩陣；U¨為系統廣義坐標對時間的二次導數，即彈性加速度矩陣；F為系統廣義力矩陣，F∈R66×1。

2 固有頻率分析

根據建立的球面2-DOF冗余驅動并聯機構的彈性動力學模型，可得無阻尼自由振動方程：

由式（12）可得到求解2-DOF冗余驅動并聯機構的固有頻率方程：

由上式可求出第i階固有頻率ωi，而且可看出，機構的無阻尼固有頻率只與剛度矩陣和質量矩陣有關，而剛度矩陣和質量矩陣都與機構的結構尺寸和輸出軸在工作空間內的位姿有關，可以從改善剛度矩陣和質量矩陣的角度來提高固有頻率，從而優化系統的結構。

2.1 固有頻率的理論計算

當球面2-DOF冗余驅動并聯機構的輸出軸在初始位置時，由式（13）可以算出機構各階固有頻率的理論值，見圖7。

圖7 機構各階固有頻率Fig.7 Each order natural frequency of the mechanism

球面2-DOF并聯機器人的無阻尼固有頻率只與剛度矩陣和質量矩陣有關，而質量矩陣與剛度矩陣與機器人的姿態有關，因此，在工作空間內，機器人的各階固有頻率會隨著機器人輸出軸的姿態變化而變化。給定工作空間為φ∈(0°,20°)、φ∈(0°,360°)，并將工作空間轉換到直角坐標系下：

由式（13）可以算出工作空間內各階固有頻率，見圖8。由圖8可知，球面2-DOF冗余驅動并聯機器人的第一、二階固有頻率均在機器人處于初始位置時取得最大值，而第三階固有頻率在隨姿態角φ增大的工作空間邊緣取得最大值，但整體上前3階固有頻率在機器人的工作空間內的變化量較小。

圖8 機構在工作空間內的各階固有頻率Fig.8 Each order natural frequency of the mechanism in the workspace

2.2 固有頻率的仿真計算

在有限元分析軟件ANSYS Workbench中，利用其ANSYS Modal模塊功能可以進行樣機的模態分析。由于前文的理論分析部分不考慮基座對樣機模態分析的影響，故只保留三維模型中的桿件部分，在初始位置時，將簡化后的模型導入ANSYS，然后進行模態計算前處理：①定義材料與單元類型；②對桿件進行網格劃分（圖9）；③定義約束與轉動副連接。最后進行模態仿真計算，得到機構在初始位置時的各階固有頻率（圖7）的仿真值。通過圖7中初始位置的各階固有頻率的有限元仿真結果與理論計算結果的對比，可以看出，本文理論模型所得結果與有限元仿真結果的變化趨勢基本一致，可以在一定程度上證明本文模型的正確性與有效性。但是由于ANSYS仿真采用的是自由網格，而且網格劃分更加密集，使得初始位置時的各階固有頻率的有限元仿真結果比理論計算結果偏大。

圖9 機構連桿網格Fig.9 Grid chart of connecting rods

3 動態響應分析

3.1 彈性動力學方程的求解

為了方便計算，將式（11）改寫為

綜合考慮系統方程式求解的穩定性、精度和計算速度等方面因素，本文采用Newmark法求解球面2-DOF冗余驅動并聯機器人的動態響應，Newmark法是根據拉格朗日中值定理對t+?t時刻的速度向量作了近似假設，在參數滿足一定條件時，Newmark法是無條件穩定的。當已知系統方程在時刻t的解，推導t+?t時刻系統方程的解時，Newmark算法的積分形式為

其中，0≤μ1≤1和 0≤μ2≤0.5是根據積分精度和穩定性要求來確定的參數。由算法穩定性分析知，當 μ1≥0.5，μ2≥0.25(μ1+0.5)2時，Newmark積分是無條件穩定的，這時可以根據精度要求選擇時間步長?t。主要計算步驟［3］如下。

（1）通過對機器人系統進行建模分析，得到系統總質量矩陣M、總剛度矩陣K和阻尼矩陣C；

（2）根據系統的運動特性分析，獲得系統初始狀態向量 U0、U˙0、U¨0；

（3）選擇合適的時間步長?t以及參數 μ1和μ，并計算下列有關常數

2

（4）計算系統有效剛度矩陣：

（5）計算t+?t時刻的有效載荷向量：

（6）分別計算t+?t時刻的位移、速度、加速度：

3.2 數值計算

設置時間步長?t=0.01 s以及參數μ1=0.6和μ2=0.4，運動時間T=10 s。已知球面2-DOF并聯機器人樣機的材質為鋼，密度 ρ=7 800 kg/m3，拉壓彈性模量 E=210 GPa；剪切彈性模量G=80 GPa；各連桿都是矩形截面，長h=8 mm，寬b=20 mm。

球面2-DOF冗余驅動并聯機器人末端軌跡見圖10，其軌跡分為3段：第1段為 φ=0°，φ為0°～20°，用時2 s；第2段為 φ=20°，φ 為 0°～360°，用時6 s；第 3段為 φ=0°，φ 為 20°～0°，用時 2 s。為保證末端位移、速度、加速度平滑連續，每一段路徑均采用組合正弦函數軌跡規劃方法［9］對其進行軌跡規劃，得到機構輸出軸末端的理論位移、速度、加速度。由式（20）可以計算得到球面2-DOF冗余驅動并聯機器人輸出軸末端的彈性線位移動態響應曲線，見圖11。球面5R并聯機器人與球面2-DOF冗余驅動并聯機器人的彈性動力學模型的建模方法一致，只是少1個分支，同理可以得到球面5R并聯機器人輸出軸末端的彈性線位移動態響應曲線，見圖12。

圖10 機器人末端軌跡曲線Fig.10 Terminal trajectory curve of the robot

由圖11可以看出，系統中連桿的彈性變形會對球面2-DOF冗余驅動并聯機器人的輸出軸端點線位移產生影響，x方向線位移的最大誤差為0.069 mm，y方向線位移的最大誤差為0.042 mm，z方向線位移的最大誤差為0.058 mm。

圖11 球面2-DOF冗余驅動并聯機器人末端彈性位移曲線Fig.11 Elastic displacement curve of the spherical 2-DOF redundantly actuated parallel robot's end

圖12 球面5R并聯機器人末端彈性位移曲線Fig.12 Elastic displacement curve of the spherical 5R parallel robot's end

由圖12可以看出，連桿的彈性變形會使球面5R并聯機器人的輸出軸端點的線位移產生不可忽略的誤差，x方向線位移的最大誤差為0.12 mm，y方向線位移的最大誤差為0.032 mm，z方向線位移的最大誤差為0.062 mm。通過對比圖11與圖12中球面5R并聯機器人與球面2-DOF冗余驅動并聯機器人輸出端的彈性位移，可以看出輸出端彈性位移均沿x方向取得最大。由于本文所建的球面2-DOF冗余驅動并聯機器人彈性動力學模型考慮了空間彎桿的彎曲、拉伸和扭轉變形，情況比較復雜，而且與球面5R非冗余并聯機構比較，冗余分支的引入增加了1個分支在輸出端點處的彈性變形協調方程，使得3個分支在輸出端點處彈性位移之間在各方向不同程度互相抵消，導致本文機構與5R機構的動態響應位移偏差在不同方向的優劣不同。同時由于冗余分支引入了變形協調方程，即過約束，可以明顯提高機器人的輸出剛度，減小機器人輸出端的3個最大輸出彈性位移，從而減小了因彈性變形引起的誤差，即提高了球面2-DOF冗余驅動并聯機器人的精度。

4 結論

（1）基于有限元法建立了球面2-DOF冗余驅動并聯機器人的彈性動力學模型。

（2）比較了球面2-DOF非冗余與冗余驅動并聯機器人輸出端的動態響應特性，證明了冗余驅動的引入可以明顯減小機器人的輸出端的最大彈性位移，提高機器人的彈性動力學特性。

（3）采用本文方法可以對機器人系統結構進行優化，同時可以快速分析并聯機器人在整個工作空間下的動態特性，為包含空間彎桿的并聯機器人的動態特性設計提供了理論基礎。

［1］ ZHAO Yongjie，GAO Feng，DONG Xingjian，et al.Elastodynamic Characteristics Comparison of the 8-PSS Redundant Parallel Manipulator and Its Non-re?dundant Counterpart the 6-PSS Parallel Manipulator［J］.Mechanism and Machine Theory，2010，45（2）：291-303.

［2］ LIU Shanzeng，YU Yueqing，ZHU Zhencai，et al.Dy?namic Modeling and Analysis of 3-RRS Parallel Manip?ulator with Flexible Links［J］.Journal of Central South University，2010，17：323-331.

［3］ 劉善增，余躍慶，蘇麗穎，等.3-RRS并聯柔性機器人的動力學建模與頻率特性分析［J］.中國機械工程，2008，19（10）：1219-1224.LIU Shanzeng，YU Yueqing，SU Liying，et al.Dynam?ics Modeling and Frequency Analysis of a 3-RRS Flex?ible Parallel Manipulator［J］.China Mechanical Engi?neering，2008，19（10）：1219-1224.

［4］ 胡俊峰，張憲民，朱大昌，等.柔性并聯機器人動力學建模［J］.農業機械學報，2011，42（11）：209-211.HU Junfeng，ZHANG Xianmin，ZHU Dachang，et al.Elastodynamics Modeling of Flexible Parallel Robot［J］.Chinese Journal of Agricultural Machinery，2011，42（11）：209-211.

［5］ 落海偉，張俊，王輝，等.3-RPS并聯機構彈性動力學建模方法［J］.機器人，2014，36（6）：739-742.LUO Haiwei，ZHANG Jun，WANG Hui，et al.Elastody?namics Modeling Methode of 3-RPS Flexible Parallel Mechanism［J］.Robot，2014，36（6）：739-742.

［6］ ZHANG Jun，ZHAO Yanqin.Elastodynamic Modeling and Joint Reaction Prediction for 3-PRS PKM［J］.Jour?nal of Central South University，2015，22（8）：2971-2979.

［7］ 周鑫，許允斗，姚建濤，等.5-UPS/PRPU冗余驅動并聯機床彈性動力學分析［J］.中國機械工程，2016，27（9）：1139-1141.ZHOU Xin，XU Yundou，YAO Jiantao，et al.Elastody?namics Analysis of a 5-UPS/PRPU Redundantly Actu?ated Parallel Machine Tool［J］.China Mechanical En?gineering，2016，27（9）：1139-1141.

［8］ 盛東義，趙勇，王皓，等.Tricept并聯機構的彈性動力學性能［J］.機械設計與研究，2012，28（6）：24-27.SHEN Dongyi，ZHAO Yong，WANG Hao，et al.Kine?to-Elastic-Dynamics（KED）Analysis of Tricept Paral?lel Mechanism［J］.Machine Design and Research，2012，28（6）：24-27.

［9］ 李永泉.球面2-DOF并聯機構的理論分析與實驗研究［D］.秦皇島：燕山大學，2012.LI Yongquan.Theory Analysis and Experiment Re?search of Spherical 2-DOF Parallel Manipulator［D］.Qinhuangdao：Yanshan University，2012.