機械臂系統的動力學建模與分析

蘇華勇，李海艷

（廣東工業大學機電工程學院，廣州 510006）

0 引言

機械臂具有控制范圍廣及其操作靈活性高等優點，在工業各領域中得到廣泛應用，并且機械臂系統在當下的發展需求材料更輕質、操作范圍更廣及運行精度更高[1]。因此，需考慮機械臂在運動過程中各節臂發生彈性變形對位置控制的影響，這就是典型的柔性多體動力學問題，其主要特點是機械臂系統在運動的過程中不但經歷著大的剛性轉動，也伴隨著變形運動，并且這兩種運動是高度耦合的[2]。機械臂的動力學建模是機械臂運動分析的基礎，典型的動力學分析方法主要有拉格朗日法、牛頓-歐拉法和Kane 方法等等，目前較常用的是拉格朗日法，它在分析過程中不必求機械臂的內作用力而只需求其對應的速度[3]。由于彈性變形有無限多自由度，動力學模型的精確解是無法得到的，通常將彈性變形離散成有限自由度來近似[4-6]。典型的離散方法有模態分析法[7]和瑞利-里茲法[8]等等，目前常用的是通過瑞利-里茲法離散彈性變形。為了方便模型的求解，常取彈性變形離散形式前兩階近似，但沒有分析不同階數對模型求解準確度的影響。此外，動力學模型屬于高度復雜的非線性方程組，從而忽略了模型中一些小變形項以簡化分析的難度，但這同樣影響了動力學模型求解的準確度。

本文以三節臂機械臂為研究對象，如圖1 所示。通過拉格朗日法和瑞利-里茲法推導出準確的動力學模型，分析彈性變形對機械臂運動的影響及彈性變形離散不同階數對模型求解準確度的影響。

圖1 三節臂機械臂

1 機械臂動力學建模

在拉格朗日法中，只要求出機械臂系統的總動能Ek和總勢能Ep，得到它們的差值L＝Ek-Ep，最后建立拉格朗日動力學方程：

式中：Qi為系統的廣義坐標；Q˙i為對應的廣義速度；Fi為系統的廣義力；n為完整約束方程的個數。

1.1 參考坐標系

如圖2 所示，分析機械臂系統首先需要建立一系列的坐標系來表示各機械臂的位置，其中xoy 為整體坐標系，x1o1y1、x2o2y2、x3o3y3為局部坐標系，局部坐標系分別以各臂桿產生彈性變形的首末端的連線作為坐標系的橫軸，這樣建立坐標系的好處是各臂桿自身產生的彈性變形不會影響到其他臂桿的位置描述，從而降低建立動力模型的復雜度。

圖2 機械臂坐標系示意圖

在機械臂系統中，把各臂桿均考慮為均質桿，驅動電機及機械手夾持裝置均考慮為集中質量，并設臂桿i 的質量為mi，長度為li，電機i的質量為mi-1，機械手夾持裝置的質量為M3，臂桿i 在整體坐標系中的偏轉角為θi，將臂桿的位置向量從局部坐標系過渡到整體坐標系需要借助旋轉矩陣Ai：

其中i＝1，2，3。

1.2 機械臂的總動能

假設機械臂系統各臂桿在相對坐標系上一點的位置向量為pi[ xi, 0 ]T，則產生彈性變形wi( xi, t )后的位置為pi′＝[xi,wi(xi,t)]T，末端位置表示為 pli＝[li,0]T，借助旋轉矩陣可將在局部坐標系上的位置轉換成在整體坐標系上的位置，轉換坐標系后的位置向量分別表示如下：

通過對位置向量關于時間t 求導后便可得到對應的速度。因此，臂桿的總動能表示為：

式（6）的前一項表示臂桿i 在運動時產生的動能，后一項表示表示相應臂桿末端負載的動能。

1.3 機械臂的總勢能

機械臂在運動過程中，總勢能的變化除了各臂桿運動時產生的重力勢能Eg1，還有臂桿產生彈性變形產生的重力勢能Eg2和變形產生的應變能Es，各勢能可以表示為：

式中：EIi為臂桿i的剛度；E為彈性模量；I為慣性矩。

因此，機械臂總勢能可以表示為：

1.4 瑞利-里茲法

由于機械臂臂桿的彈性變形是與距離和時間有關的量，通常需要將彈性變形離散成有限自由度來近似，典型的離散方法有模態分析法和瑞利-里茲法。目前常用的是瑞利-里茲法，它也是彈性連續體力學中求近似解的最基本方法。該方法是將彈性變形項離散成里茲基函數與對應時間系數乘積的級數展開式，以往的方法一般都是取級數展開式的前2階進行近似。為了進行對比驗證，此次選取了級數展開式前4階來近似彈性變形項，各臂桿離散后的彈性變形可表示如下：

1.5 建立動力學方程

將彈性變形的離散形式分別代入到動能和勢能的求解式中，最后得到的總動能和總勢能表達式都是僅與時間t有關的多項式。

拉格朗日函數可表示為L＝Ek-Ep，方程的廣義坐標表示為Q，機械臂系統的廣義坐標除了臂桿的偏轉角外，還有考慮各臂桿發生彈性變形造成的變化，它也是受時間t影響的量，在彈性變形的離散形式中以各階基函數對應的時間系數表示。方程的廣義力表示為F，在旋轉運動中的廣義力為驅動轉矩，并設電機的轉矩為T，臂桿1 受到電機1 的驅動轉矩和電機2反作用的驅動轉矩，臂桿2受到電機2的驅動轉矩和電機3的反作用驅動轉矩，臂桿3 受到電機3 的驅動轉矩。因此臂桿i受到的驅動轉矩τi為：

最后，機械臂的拉格朗日動力學方程表示如下：

式中：Q＝[θ,q]T， θ＝[θ1,θ2,θ3]， q ＝ [q11,q12,q13,q14,q21,，。

從式（13）可以看到，機械臂的動力學模型包含了15個方程，并且該模型屬于復雜非線性強耦合的方程組。它既包含了機械臂各臂桿間剛性運動的耦合，還包含了剛性運動與臂桿彈性變形的耦合。前者表示機械臂的大范圍剛性運動，而后者會使機械臂系統運動產生擾動。該模型考慮了彈性變形的影響，因此對機械臂系統描述更加全面，但是也很大程度地增加模型的求解難度。

圖3 動力學模型求解流程圖

1.6 動力學方程求解

以往的方法對方程組的求解過程中，為了方便求解就一定程度地簡化了模型，把方程中的小變形項都省略了，但是這無疑對方程的準確度有一定影響。因此，動力學方程的求解不應該舍棄小變形項對模型的貢獻。

通過Mathematica對公式推導和整理，拉格朗日方程組最后可以整合成My¨＝z 的形式，其中M15×15是個滿秩的對稱矩陣，是廣義坐標的二階導數項，將其他剩余的項整合到z15×1里面，且方程組存在唯一解。

由于該方程組屬于非線性的常微分方程組，可以借助MATLAB 軟件的ode 函數求解，但是方程一般需要滿足格式y˙＝f (t,y) 的形式，因此要對動力學方程作進一步處理并求解，流程如圖3所示。

2 數值仿真

現在通過具體的參數進行數值仿真，三節臂機械臂的長度分別為 l1＝l2＝l3＝1.5 m，各臂桿的質量分別為 m1＝m2＝3 kg，m3＝2.5 kg，電機的質量分別為M0＝M1＝M2＝5.5 kg，末端夾持裝置的質量 M3＝4.5 kg，重力加速度 g＝9.8 m/s2，剛度 EI1＝EI2＝EI3＝406 N·m2，電機驅動力矩T1＝1 150 sin(2πt)3+850(N·m)；T2＝452sin(2πt)3+450(N·m)；T3＝155 sin(2πt)3+150(N·m)。仿真時間為1.5 s，并設初始條件Q ＝Q˙＝ 0 。

為了分析產生變形前后對位置的影響，考慮將機械臂的末端位置進行對比，可以表示為：

如圖4 所示，機械臂末端位置在變形前后的曲線變化存在一定的差異，在變形前的曲線是很光滑的，機械臂彈性變形后的曲線在局部放大處更是可以看到很明顯的振動，這正是彈性變形產生的抖動現象導致的。另外，機械臂各臂桿的偏轉角在變形前變化后角度也產生了一定的偏差。在變形前的曲線也是很光滑的，在機械臂發生彈性變形后在局部放大處也是可以看到很明顯的鋸齒狀。這些現象都可以說明機械臂系統發生的彈性變形對位置控制有一定程度影響。

圖4 機械臂變形前后對比

由于彈性變形有無限多自由度，動力學模型的精確解是無法得到的，通常將彈性變形離散成有限自由度來近似。為了分析不同階數近似的彈性變形對求解的影響，用彈性變形近似的前2階和前4階進行對比，其誤差分析如圖5所示。

可以看出，選取前2階和前4階近似彈性變形對模型求解準確度的影響并不大。但是，選取展開式近似的階數越大，動力學模型規模就會更大更復雜，求解的難度也會陡增，因此選取階數多少需要根據實際情況進行綜合權衡。

圖5 彈性變形2階與4階近似誤差

3 結束語

本文通過以三節臂的機械臂作為研究對象，考慮彈性變形對機械臂系統的影響，建立動力學模型并進行數值仿真求解與分析，可以得到如下的結論。

（1）機械臂系統在運動過程中產生的彈性變形對臂桿的位置控制產生的擾動較大。因此，機械臂系統的控制有必要考慮彈性變形的影響，彈性變形的擾動特性需要進一步分析。

（2）機械臂的彈性變形離散近似取前2階和前4階對模型的求解準確度影響不大，并且階數增加會導致模型更加復雜，階數的選取也需要根據實際情況綜合權衡。