基于修正Udwadia-Kalaba方法的剛體系統(tǒng)動力學(xué)模擬

郭明明，姚文莉

(青島理工大學(xué) 理學(xué)院，青島 266525)

Udwadia-Kalaba方法是南加州大學(xué)的兩位教授，UDWADIA和KALABA，在1992年提出的一種針對受約束復(fù)雜機械系統(tǒng)動力學(xué)建模的新方法[1]。 該方法與其他傳統(tǒng)的動力學(xué)建模方法有所不同，它在解決受約束系統(tǒng)動力學(xué)問題時，不需要引入輔助變量，可以得到解析形式的動力學(xué)方程和顯示形式的加速度表達。 可以說，該方法為動力學(xué)問題的求解提供了一種新思路。

隨著研究的不斷深入，Udwadia-Kalaba方法被推廣到了廣義坐標(biāo)表達的形式[2]，它的應(yīng)用范圍也被擴展到了非理想約束系統(tǒng)[3-6]，奇異質(zhì)量矩陣系統(tǒng)[7]等一般的機械系統(tǒng)[8-9]。 現(xiàn)在，Udwadia-Kalaba方法憑借其在表達上的簡潔性和一般性，已經(jīng)被廣泛應(yīng)用于多體系統(tǒng)的動力學(xué)建模[10-19]，同時系統(tǒng)理想約束力及加速度的顯式表達方法為基于逆動力學(xué)的控制問題提供了新的思路和理論基礎(chǔ)。

然而，Udwadia-Kalaba方法在擁有上述諸多優(yōu)勢的同時，也存在著一些問題。 當(dāng)使用Udwadia-Kalaba方法來處理含有庫侖滑動摩擦的多剛體系統(tǒng)時，求得的理想約束力和非理想約束力無關(guān)，這與實際情況是不符合的。 對此，姚文莉等[20]提出了修正方法，研究了滑動狀態(tài)下理想約束力的表達方法。

對于考慮干摩擦的剛體系統(tǒng)的運動過程，滑動-黏滯及摩擦自鎖等非光滑現(xiàn)象不容忽略。 本文借助于修正的Udwadia-Kalaba方法，研究了包含滑動和黏滯兩種狀態(tài)的非光滑剛體系統(tǒng)動力學(xué)方程，并建立了數(shù)值求解的計算框架。

1 Udwadia-Kalaba方法的描述

考慮無約束的機械系統(tǒng)，其動力學(xué)方程可以表示為

(1)

給此系統(tǒng)施加m個約束：

(2)

實際上，式(2)中包含了各種常見的約束：完整或非完整約束，理想或非理想約束，定常或非定常約束等。 由于這組約束的存在，系統(tǒng)不僅要受到主動力的作用，還要受到約束力的作用，在它們的共同作用下，系統(tǒng)的運動將發(fā)生改變，系統(tǒng)的動力學(xué)方程變?yōu)?/p>

(3)

(4)

Udwadia-Kalaba方法給出了解析形式的理想廣義約束力和非理想廣義約束力的表達式：

(5)

(6)

式(6)中：C為n維向量，它需要通過對某個具體的動力學(xué)系統(tǒng)進行試驗和對比等操作才可獲得。

2 修正Udwadia-Kalaba方法的描述

(7)

對于上述兩處修正的理論推導(dǎo)和證明，這里不再贅述。 本文將結(jié)合具體實例進行動力學(xué)仿真，來證明修正Udwadia-Kalaba方法的正確性。

3 剛體系統(tǒng)動力學(xué)的滑動-黏滯處理方法

對于含有庫侖滑動摩擦的剛體系統(tǒng)，其運動過程中可能會出現(xiàn)滑動-黏滯、摩擦自鎖等非連續(xù)現(xiàn)象。 由于這些非連續(xù)現(xiàn)象，會造成系統(tǒng)動力學(xué)方程的不連續(xù)或者分段連續(xù)，給方程的建立和數(shù)值計算帶來困難。

由庫侖摩擦模型可知，剛體受到的切向摩擦力的大小與它受到的法向約束力的大小|Fn|成正比，方向與兩個剛體的相對運動或者相對運動趨勢的方向相反。 切向摩擦力Ft可以表示為

(8)

為了便于計算，F(xiàn)t與Fn的關(guān)系可以寫成：

Ft=dFn

(9)

可以發(fā)現(xiàn)，對應(yīng)不同的運動狀態(tài)，F(xiàn)t有不同的表達式。

可以發(fā)現(xiàn)，對于含有雙邊約束的系統(tǒng)，由于兩個剛體之間接觸面的變換，法向約束力Fn會出現(xiàn)正值或者負值的情況，此時摩擦力的表達式就會出現(xiàn)絕對值問題，這會造成系統(tǒng)動力學(xué)方程求解的困難。 目前對于此問題，大多數(shù)文獻都是采用基于線性互補(LCP)的方法來進行解決[21-25]。

4 非光滑剛體系統(tǒng)動力學(xué)問題的計算流程

含單點接觸的剛體系統(tǒng)動力學(xué)的計算流程如圖1所示。

圖1 計算流程

5 算例

圖2 含庫侖滑動摩擦的滑塊擺桿機構(gòu)

系統(tǒng)參數(shù)：mA=1 kg，m=2 kg，l=0.5 m，F(xiàn)0=23 N，w=π/5，c= -0.8 N·m·s，靜摩擦系數(shù)μ0=0.78，動摩擦系數(shù)μ=0.56；初始條件，θ=-π/2，系統(tǒng)從靜止開始運動。

用修正Udwadia-Kalaba方法對受約束的多體系統(tǒng)進行建模可分為三步：首先，用Lagrange方程建立不考慮約束的系統(tǒng)動力學(xué)方程；然后，對約束方程求導(dǎo)，得到二階形式的約束方程；最后，將由約束產(chǎn)生的約束力施加到不考慮約束的動力學(xué)方程上。

(10)

不考慮約束時系統(tǒng)的動力學(xué)方程可以寫成：

(11)

3)將由上述約束產(chǎn)生的約束力施加到不考慮約束的動力學(xué)方程上，

應(yīng)用修正后的Udwadia-Kalaba方法：

(12)

(13)

(14)

綜合式(8)—(14)，可以得到滑塊受到的法向約束力的表達式：

(15)

滑塊受到的切向摩擦力可由Ft=dFn求得。

應(yīng)用未修正的Udwadia-Kalaba方法，綜合式(5)(10)(11)(14)，可以得到滑塊受到的法向約束力的表達式：

(16)

滑塊受到的切向摩擦力可由Ft=dFn求得。

同時，再應(yīng)用經(jīng)典的剛體動力學(xué)方法，列出此系統(tǒng)的動力學(xué)方程：

(17)

(18)

綜合式(17)(18)，可以得到滑塊受到的法向約束力的表達式：

(19)

滑塊受到的切向摩擦力可由Ft=dFn求得。

分別用三種方法建立了含摩擦滑塊擺桿機構(gòu)的動力學(xué)方程，用Matlab進行數(shù)值仿真，結(jié)果如圖3—8所示。

由圖3可以看出，由于庫侖滑動摩擦的存在，系統(tǒng)出現(xiàn)了滑動-黏滯現(xiàn)象。

初始時，滑塊靜止，水平主動力F(t)從0開始增加，在F(t)的值達到滑塊的最大靜滑動摩擦力之前，滑塊受到的切向摩擦力大小等于F(t)的值，方向與F(t)的方向相反，滑塊處于黏滯狀態(tài)，擺桿沒有發(fā)生擺動。

隨后，當(dāng)F(t)的值達到滑塊的最大靜滑動摩擦力時，滑塊受到的切向摩擦力大小等于滑塊的最大靜滑動摩擦力的值，方向與滑塊的運動趨勢方向相反，滑塊處于臨界狀態(tài)，滑塊即將開始運動，擺桿也即將開始擺動。

然后，滑塊在水平主動力F(t)和切向摩擦力以及滑塊和擺桿的慣性力作用下運動，滑塊處于滑動狀態(tài)。 由于F(t)是一個周期變化的力，所以它會導(dǎo)致系統(tǒng)的受力情況發(fā)生周期變化。

將圖3與圖4對比，可以發(fā)現(xiàn)，圖3在滑塊處于滑動狀態(tài)時，法向約束力Fn的值在某個數(shù)值附近發(fā)生平穩(wěn)的波動，這是由于擺桿的輕微擺動引起的，而此時摩擦力Ft也發(fā)生平穩(wěn)的波動，這符合滑動狀態(tài)時摩擦力的表達式。 圖4在滑塊處于滑動狀態(tài)時，摩擦力Ft沒有在某個數(shù)值附近平穩(wěn)的波動，這不符合動滑動摩擦力的特征。 圖5是用經(jīng)典的剛體動力學(xué)方法仿真后的結(jié)果，和圖3一致。 這是必然的結(jié)果，因為將式(15)中矩陣的元素代入，進行進一步的計算，可以得到和式(19)一致的表達式，這個已有推導(dǎo)證明[20]。 圖3—5的仿真結(jié)果表明，對于Udwadia-Kalaba方法的修正是正確的。

圖6和圖7是應(yīng)用修正Udwadia-Kalaba方法仿真的滑塊水平方向速度、加速度和滑塊位移。 在水平主動力F(t)的值達到滑塊的最大靜滑動摩擦力的值之前，滑塊水平方向速度、加速度和水平方向位移均為0；當(dāng)水平主動力F(t)的值達到滑塊的最大靜滑動摩擦力的值時，系統(tǒng)處于臨界狀態(tài)，在下一時刻，滑塊水平方向速度、加速度和水平方向位移均將大于0；隨后，滑塊的運動狀態(tài)將隨著F(t)的變化，發(fā)生周期變化。 圖7中滑塊的豎向位移始終為0，表明在此算例中滑塊下表面始終與滑道保持接觸。

圖8是應(yīng)用修正Udwadia-Kalaba方法仿真的擺桿位移和擺動角度。 初始時，擺桿豎直向下，與水平方向夾角為-π/2，隨著滑塊的運動，擺桿在慣性力和阻力矩作用下發(fā)生擺動。從仿真的數(shù)據(jù)可以看出，擺桿擺動的角度范圍是[-2.084，-1.058]，表明擺桿始終在水平面下方發(fā)生輕微的擺動，不足以引起滑塊下表面離開滑道。 由于擺桿的擺動，擺桿的水平方向位移和豎向位移都發(fā)生波動現(xiàn)象。

文中選擇的算例，在給出的系統(tǒng)參數(shù)和初始條件下，法向約束力始終為正值。 在改變滑塊和擺桿的質(zhì)量比值或者改變初始條件時，滑塊的下表面可能會離開滑道，法向約束力可能會出現(xiàn)可正可負的情況，這將是作者下一步的研究內(nèi)容。

6 結(jié)論

在解決受約束系統(tǒng)動力學(xué)問題時，Udwadia-Kalaba方法與傳統(tǒng)的動力學(xué)方法相比，可以不引入輔助變量，獲得解析的動力學(xué)方程和顯示的加速度表達，操作更加方便，可以為動力學(xué)問題的求解提供一種新思路。

對于考慮干摩擦的剛體系統(tǒng)，系統(tǒng)中的非理想約束力和理想約束力是相互耦合的，本文應(yīng)用修正的Udwadia-Kalaba方法，對包含滑動和黏滯兩種狀態(tài)的非光滑剛體系統(tǒng)動力學(xué)方程進行了數(shù)值模擬，證明了修正Udwadia-Kalaba方法的正確性和有效性。 修正后的Udwadia-Kalaba方法可以更可靠地應(yīng)用于解決非光滑剛體系統(tǒng)動力學(xué)問題。