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

郭明明，姚文莉

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

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

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

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

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

1 Udwadia-Kalaba方法的描述

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

(1)

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

(2)

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

(3)

(4)

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

(5)

(6)

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

2 修正Udwadia-Kalaba方法的描述

(7)

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

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

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

由庫(kù)侖摩擦模型可知，剛體受到的切向摩擦力的大小與它受到的法向約束力的大小|Fn|成正比，方向與兩個(gè)剛體的相對(duì)運(yùn)動(dòng)或者相對(duì)運(yùn)動(dòng)趨勢(shì)的方向相反。 切向摩擦力Ft可以表示為

(8)

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

Ft=dFn

(9)

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

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

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

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

圖1 計(jì)算流程

5 算例

圖2 含庫(kù)侖滑動(dòng)摩擦的滑塊擺桿機(jī)構(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，動(dòng)摩擦系數(shù)μ=0.56；初始條件，θ=-π/2，系統(tǒng)從靜止開始運(yùn)動(dòng)。

用修正Udwadia-Kalaba方法對(duì)受約束的多體系統(tǒng)進(jìn)行建?？煞譃槿剑菏紫?，用Lagrange方程建立不考慮約束的系統(tǒng)動(dòng)力學(xué)方程；然后，對(duì)約束方程求導(dǎo)，得到二階形式的約束方程；最后，將由約束產(chǎn)生的約束力施加到不考慮約束的動(dòng)力學(xué)方程上。

(10)

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

(11)

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

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

(12)

(13)

(14)

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

(15)

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

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

(16)

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

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

(17)

(18)

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

(19)

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

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