基于Lagrange方程和二階影響系數矩陣的Delta并聯機構約束力分析及仿真研究*

張 勇,郭志飛*,甄希金,張敬芳

(1.河北機電職業技術學院 機械工程系,河北 邢臺 054001;2.上海船舶工藝研究所,上海 200032)

0 引 言

并聯機構具有剛度好、運動精度高、易在線控制等優點,在許多工業領域的應用廣泛[1]。對其進行動力學分析對于并聯機構的結構設計、動力學性能優化和實時控制都具有重要的意義[2]6-7。

對并聯機構的動力學分析需要系統的方法[2]119,并聯機構的力學分析也是并聯機構研究的基本問題[3-5]。目前,國內外在并聯機構動力學分析方面取得了較大的成果,常用的并聯機構動力學分析方法有Lagrange方法、Newton-Euler法和Kane方法[6-7]。YANG D C H[8]利用Newton-Euler法研究了Stewart平臺的逆動力學問題;但該研究忽略了關節的摩擦和桿的力矩慣量,并且其假設桿的中心位于桿的軸線上。JI Zhi-ming[9]在1994年建立了Stewart平臺的動力學方程,并且該方程考慮了腿部轉動慣量對Stewart平臺的影響。CODOUREY A[10]介紹了基于拉格朗日方程和虛功原理的并聯機構動力學建模方法,其中Lagrange方法是以能量為基礎建立的動力學的數學模型,該方法能夠得到形式簡潔的動力學方程,十分適用于并聯機構的動力學分析[11];但是,該方法的不足之處是其無法進行約束力分析,Kane方法亦然。在多數情況下,Newton-Euler法被用于對并聯機構的約束力進行分析研究,但其過程復雜、方程數量龐大、易出錯[12,13]。

綜上,本文開展基于逆動力學的并聯機構約束力的分析研究,通過建立三平移Delta并聯機構的Lagrange動力學模型,使用二階影響系數矩陣求解滑塊的加速度,對三平移Delta并聯機構驅動滑塊的約束力進行分析和求解,為三平移Delta并聯機構結構設計及控制系統設計提供理論依據。

1 Lagrange動力學模型建立

本文所研究的三平移動Delta并聯機構如圖1所示。

圖1 三平移Delta并聯機構

圖1中,該并聯機構的主要組成部分有:固定平臺、導軌、滑塊、連桿、動平臺等。當3個滑塊以不同規律運動的時候,動平臺就可以實現沿X、Y、Z方向的平動,以及各種耦合形式的運動。

為了便于后續分析,筆者將該并聯機構進行簡化,機構位置結構簡圖如圖2所示。

圖2 機構位置結構簡圖

圖2中,固定平臺簡化為等邊三角形A1A2A3,其外接圓半徑為R;動平臺簡化為等邊三角形B1B2B3,其外接圓半徑為r;在動平臺上建立動坐標系,動坐標系原點位于動平臺B1B2B3的中心,x軸垂直于B1B2,y軸平行于B1B2,z軸由右手坐標系得到;在固定平臺上建立靜坐標,原點O位于三角形A1A2A3的角平分線的交點,坐標軸的建立和動平臺相同。

該并聯機構的主要結構參數如表1所示。

表1 并聯機構結構參數

本文研究的并聯機構的材料為45#鋼,動平臺的質量ma=0.336 kg,滑塊的質量為m1=0.167 kg,連桿的質量為m2=0.028 kg。

通過對該機構分析,得到該機構的Lagrange動力學模型[14-16],表示廣義坐標x、y、z方向上廣義力:

(1)

記為:

(2)

通過力的雅各比矩陣得出在驅動滑塊所需的驅動力表達式為:

(3)

鑒于文章篇幅,此處不再詳細列出其推導過程。

2 三平移Delta并聯機構加速度分析

2.1 三平移Delta機構的二階影響系數矩陣

根據三平移Delta并聯機構的性質可知,當確定3個驅動滑塊的輸入位移d1、d2、d3,動平臺的位姿可以由其質心坐標x、y、z唯一確定。

動平臺的位姿可以表示為U={x,y,z}T,對應的3個驅動滑塊的位置為d={d1,d2,d3}T,如圖3所示。

圖3 三平移Delta并聯機構簡化圖

根據圖3可列出位移方程組:

(4)

輸入位移d1、d2、d3和輸出廣義坐標x、y、z均隨時間變化,將式(4)對時間進行隱函數求導,可以得到:

(5)

整理后可以得到其矩陣形式,即:

(6)

(7)

式中:矩陣[G]—該三平移Delta并聯機構的一階影響系數矩陣(Jacobian矩陣)[17]182-185。

式(7)再次對時間求導,經簡化整理后可得:

(8)

式(8)可以寫成矩陣形式:

(9)

式中:[H]—該并聯機構的二階段影響矩陣(Hessian矩陣)[17]186-190。

從式(8,9)可得到二階段影響矩陣的表達式為:

(10)

將二階段影響矩陣展開,有:

(11)

由式(11)可以看出,二階影響系數矩陣[H]是一個三維立體矩陣,即二階影響系數矩陣中的每個元素均為一個3×1的列向量。每個元素表達式為[17]189-190:

(12)

2.2 三平移Delta并聯機構加速度分析和計算

設執行機構的加速度表示為:a0={ax,ay,az}T,3個滑塊的加速度為:ad={a1,a2,a3}T,根據式(9)可以求出驅動滑塊加速度,即:

(13)

式中:v0—動平臺的速度，v0={vx,vy,vz}T;[J]—一階影響系數矩陣,式(7)中矩陣[G]的逆矩陣;[Hd]{d=1,2,3}—二階影響系數矩陣。

[Hd]展開后的表達式為:

(14)

(15)

(16)

在該并聯機構的動平臺添加空間圓錐螺旋線的運動軌跡,其曲線方程為:

(17)

在該運動軌跡下,使用MATLAB求解三平移Delta機構3個驅動滑塊的加速度,得到的加速度時間曲線如圖4所示。

圖4 三平移Delta機構加速度時間曲線(MATLAB)

圖4中,從左向右依次為驅動滑塊1、滑塊2和滑塊3的加速度時間曲線.

2.3 三平移Delta并聯機構加速度仿真分析

在ADAMS中建立該并聯機構的虛擬樣機,如圖5所示。

圖5 三平移Delta機構加速度時間曲線

在并聯機構虛擬樣機的動平臺上添加運動驅動,運動驅動軌跡方程如式(17)所示。圖5為Delta并聯機構的運動仿真,顯示出動平臺的運動軌跡為空間圓錐螺旋線[18-19]。

仿真完成后,可得3個驅動滑塊的加速度時間曲線,如圖6所示。

圖6 三平移Delta機構加速度時間曲線(ADAMS)

圖6中,從左至右依次為滑塊1、滑塊2、滑塊3的加速度時間曲線。

在圖4中,t=1 s時,滑塊一的加速度a1=0.079 9 m/s2;t=2 s時,滑塊二加速度a2=-0.045 m/s2;t=3 s時,滑塊三的加速度為a3=0.21 m/s2。

而在圖6中,即ADAMS的仿真結果中,t=1 s,滑塊一的加速度a1=0.079 9 m/s2;t=2 s時,滑塊二的加速度a2=0.045 m/s2;t=3 s時,滑塊三的加速度a3=0.21 m/s2。可以看出,ADAMS仿真得到加速度的值與在MATLAB數值計算的結果吻合,驗證了該加速度逆解數學模型的正確性,也證明了ADAMS仿真的正確性。

3 驅動滑塊約束力分析

3.1 三平移Delta并聯機構的約束力分析基礎

為方便三平移Delta并聯機構約束力的分析,需要做一些準備工作,具體如下:

(1)以驅動滑塊1為例,設動平臺中心坐標在靜坐標系中表示為Om(xm,ym,zm),與驅動滑塊1相連的連桿B1C1兩端的坐標點分別用B1和C1表示:

(18)

(19)

可以得到連桿B1C1的向量為:

(20)

根據并聯機構的運動學反解,在已知動平臺位姿時,3個驅動滑塊的Z方向位移可以表示為d1、d2和d3:

(21)

篇幅所限，這里不再給出其具體求解過程。

(2)該Delta并聯機構的3個分支均為4S閉環,并且為平行四邊形,如圖7所示。

圖7 Delta機構4S平行四邊形閉環

球鉸1、2處連接驅動滑塊,球鉸3、4處連接動平臺,根據螺旋理論可以知道,球鉸1、2處的約束力方向均是沿著分桿13和分桿24的軸向[20],本文將該并聯機構進行簡化。

在圖7中,筆者將4S平行四邊形機構簡化成一個連桿,即球鉸1和2合并,3和4合并。

驅動滑塊受力分析圖8所示。

圖8 驅動滑塊受力分析圖

(3)連桿B1C1約束力分析如圖9所示。

圖9 連桿B1C1約束力分析

圖9中,在連桿B1C1中,作用在滑塊1中球鉸的約束力方向沿向量連桿B1C1軸線方向。為了方便求解連桿兩端球鉸水平方向約束力,筆者在連桿B1C1兩端球鉸球心處添加兩個局部坐標系Ol-xlylzl和Op-xpypzp,兩個局部坐標系的X、Y、Z軸均和世界坐標系的X、Y、Z軸方向一致。

3.2 約束力豎直方向分力求解及仿真

根據前面的分析,可以得到驅動滑塊豎直方向受力情況,如圖10所示。

圖10 驅動滑塊受力分析圖

該并聯機構的材質為鐵,滑塊1上的驅動力為F1,加速度為ad1,球鉸處對驅動滑塊在豎直方向的約束力分力用F21表示。根據圖10,在忽略摩擦力的前體下,滑塊1豎直方向的受力為:

F1+Fg+F21=ma

(22)

筆者在MATLAB對式(22)進行求解,得到球鉸豎直方向的分力如圖11所示。

圖11 MATLAB中驅動滑塊1處球鉸約束豎直方向分力

筆者在ADAMS中對該并聯機構進行動力學仿真,給3個驅動滑塊添加驅動,動平臺運動軌跡如式(17)。

對驅動滑塊處球鉸約束力豎直分力進行求解,得到球鉸的約束力豎直方向受力如圖12所示。

圖12 ADAMS中驅動滑塊1處約束力的豎直方向分力

從圖11、圖12可以看出:兩者的結果吻合,驗證了Lagrange方程和二階影響系數矩陣對球鉸連接處豎直方向約束力分析結果的正確性。

3.3 約束力水平方向分力分析

對約束力水平方向分力的分析將會涉及到連桿,需要建立連桿的質心連體坐標系ol-xlylzl和質心平動坐標系op-xpypzp,對連桿的轉動慣量進行坐標間變換,這里研究對象為連桿B1C1。

圖13 連桿的連體坐標系和質心平動坐標系

質心連體坐標系的原點仍定義在連桿B1C1質心處,質心連體坐標系的Z軸正向單位向量定義如式(23)所示,根據該并聯機構的結構尺寸,得到其表達式為式(24),為了便于表達,Z軸正向單位向量在靜坐標系各軸分量用zl1、zl2、zl3表示,即:

(23)

(24)

坐標系ol-xlylzl的其他坐標軸的正向單位向量定義表達式為:

(25)

(26)

根據式(24～26),可以確定出質心連體坐標系o1-x1y1z1相對于質心平動坐標系的方向余弦矩陣為:

(27)

根據剛體的慣量張量在連心坐標系間的變換公式[21],已知連桿B1C1在質心連體坐標系下的慣性張量Il和連體坐標系相對于質心平動坐標系的方向余弦矩陣Cpl,可以求出連桿B1C1在質心平動坐標系的慣性張量Ip。連體坐標系O1-x1y1z1為連桿的慣性主軸坐標系,可知連桿在質心平動坐標系中的轉動慣性張量Ip,即:

Ip=[Cpi]*Il*[Cpi]T

(28)

經分解后得到的連桿所受約束力如圖14所示。

圖14 連桿的所受約束力分解

根據上述分析以及剛體一般運動微分方程(8),可以寫出連桿沿質心平動坐標系的xp軸方向的受力平衡方程和繞yp軸的力矩平衡方程,即:

(29)

(30)

通過MATLAB對式(29)進行解算,并在ADAMS中獲取對應球鉸的水平方向X軸約束力。

MATLAB中球鉸C1的X方向分力如圖15所示。

圖15 MATLAB中計算得出球鉸C1的X方向分力

ADAMS中球鉸C1的X方向分力如圖16所示。

圖16 ADAMS中計算得出球鉸C1的X方向分力

在該并聯機構的ADAMS仿真模型中,一個分支有4個球鉸,而文中的數學計算模型對分支球鉸進行了合并,故ADAMS仿真對一個分支中同一端的兩個球鉸在X方向的約束分力進行求和。

由圖(15，16)可看出:在MATLAB和ADAMS中,計算得到的球鉸約束力的X方向分力比較接近,存在的細微偏差是由于合理假設導致的,該結果驗證了本文對約束力分析的正確性。

球鉸約束力的Y、X方向分量分析過程類似,此處不再重復。

4 實驗驗證

為了驗證本文的結論,筆者使用現有的三平移Delta并聯機構運動平臺進行實驗,如圖17所示。

圖17 三平移Delta并聯機構

圖17中,該并聯機構的控制單元采用Arduino板,對執行機構進行軌跡規劃,使其運動軌跡如圖5中的螺旋線。

筆者使用MPU-6050加速度計采集驅動滑塊1的豎直方向的加速度數據,并使用均值濾波得到較平滑的數據;在MATLAB中打開該加速度文件,并繪制折線圖,如圖18所示。

圖18 MPU-6050加速度計采集的滑塊1的加速度數據

筆者根據式(22)和滑塊1的豎直方向加速度數據得到其豎直方向的合外力,如圖19所示。

圖19 實驗所得的滑塊1所受合外力

從圖19中可以看到:驅動滑塊1的合外力和圖11中通過公式計算所得的合外力數據趨勢吻合;其最大值、最小值也比較接近。該結果驗證了本文結論的正確性。

5 結束語

本研究使用Lagrange方程和二階影響系數矩陣,建立了三平移Delta并聯機構驅動滑塊球鉸處的約束力方程,運用MATLAB對約束力方程進行了求解,并用ADAMS方程對約束力方程進行了驗證;最后,通過實驗數據對仿真結果進行了驗證。

研究結果表明:

(1)使用Lagrange方程和二階影響系數矩陣建立的約束力方程,其表達式規范、簡潔,方程數量減少約1/2;

(2)將二階影響系數矩陣和Lagrange方程結合,分析并聯機構的約束力,簡化了約束力的推導過程,得到的約束力分析結果準確,有一定程度的方法創新。

下一階段,筆者將會把本研究成果推廣應用到其他構型的并聯機構,以提升該約束力算法的通用性、穩定性和運算速度;并且,筆者還將致力于開發具有自主知識產權的并聯機構專用三維動力學分析軟件系統。