高云峰

(清華大學(xué)航天航空學(xué)院,北京100084)

在傳統(tǒng)的運(yùn)動學(xué)教學(xué)中，為了避免復(fù)雜的求導(dǎo)而陷入純數(shù)學(xué)演算，通常強(qiáng)調(diào)基點(diǎn)法和點(diǎn)的復(fù)合運(yùn)動等具有物理含義的方法，求解系統(tǒng)在特定時(shí)刻或位置的運(yùn)動信息[1]。這樣處理的優(yōu)點(diǎn)是強(qiáng)調(diào)物理概念，缺點(diǎn)是對系統(tǒng)的整體運(yùn)動不容易了解，某些點(diǎn)的運(yùn)動軌跡有時(shí)難以想象。這樣訓(xùn)練出來的學(xué)生，對于機(jī)械運(yùn)動缺乏了解，也難以勝任未來的設(shè)計(jì)工作。

利用運(yùn)動學(xué)計(jì)算機(jī)輔助分析方法和相關(guān)軟件，可以很容易演示系統(tǒng)的整個(gè)運(yùn)動過程、求出任意點(diǎn)或剛體在任意時(shí)刻的速度、加速度、角速度、角加速度等運(yùn)動量。Matlab是實(shí)現(xiàn)計(jì)算機(jī)輔助分析的工具之一，可以幫助學(xué)生加深理解機(jī)械運(yùn)動的特點(diǎn)，為將來學(xué)生自己設(shè)計(jì)裝置提供了方便。

本文著重介紹2個(gè)典型機(jī)構(gòu)的運(yùn)動軌跡、速度和加速度分析，以及利用Matlab中的符號推導(dǎo)功能，從位移函數(shù)求導(dǎo)得到速度和加速度的表達(dá)式，并進(jìn)行畫圖和動畫顯示。

在傳統(tǒng)點(diǎn)的復(fù)合運(yùn)動中，要特別注意動點(diǎn)、動系的選擇，如果選擇不當(dāng)，就分析不下去，因此對很多同學(xué)而言有一定的難度。而采用Matlab，可以采用統(tǒng)一的方法，只要列出一般位置的表達(dá)式就可以利用符號推導(dǎo)功能，得到速度、加速度等運(yùn)動量的表達(dá)式。

1 M atlab中非線性方程的求解及動畫演示

案例1：如圖1，已知四連桿機(jī)構(gòu)ABCD，AB桿長為a1，BC桿長為a2，CD桿長為a3，AD距離為a4。若AB桿以勻角速度ω0轉(zhuǎn)動，初始θ0=0。求BC和CD桿的角度、角速度變化規(guī)律。

圖1 四連桿機(jī)構(gòu)

這個(gè)問題中尺寸不是特別給定的值，一般學(xué)生分析起來會有困難。而用Matlab可以進(jìn)行規(guī)范化處理：建立坐標(biāo)Axy，θ=ω0t已知，廣義坐標(biāo)為φ1和φ2。系統(tǒng)的約束方程為[2]

方程(1)是關(guān)于轉(zhuǎn)角φ1和φ2的非線性方程組，通常沒有解析解，下面給出一般的處理方法。

設(shè)有非線性方程組f i(x1,x2,···,x n),i=1,2,···,n，沒有解析解，可用多種數(shù)值計(jì)算的方法求解，這里以牛頓?拉普森法為例，其主要步驟如下：

步驟(1)：給定計(jì)算中的允許誤差ε，給出一組粗略的估計(jì)值

步驟(2)：計(jì)算f i(x?)，判斷|f i(x?)|<ε是否成立，若成立，當(dāng)前估計(jì)值即為一組數(shù)值解，求解結(jié)束。若不成立，到步驟(3)。

步驟(3)：計(jì)算非線性方程組的雅可比矩陣J(x)，代入當(dāng)前估計(jì)值得到J(x?)。多元函數(shù)f i(x1,x2,···,x n)的雅可比矩陣定義為

步驟(4)：類似一元函數(shù)的泰勒展開式，f(x)=f(x0)+f′(x0)(x?x0)+o(x?x0)，多元函數(shù)為

略去高階小量，從而得到關(guān)于修正值dx=[dx1,dx2,···,dx n]T的一組線性方程組

步驟(5)：解線性方程組，得到修正值dx。

步驟(6)：得到新的估計(jì)值x?=x?+dx，返回到步驟(2)。

具體回到案例1，把方程(1)改寫為

已知a1,a2,a3,a4和θ，估計(jì)值是計(jì)算雅可比矩陣，并代入當(dāng)前估計(jì)值

得到一組關(guān)于修正值[dφ1,dφ2]T的線性方程組

類似代碼X=inv(A)*B，可解出[dφ1,dφ2]T，總之按上述步驟1～6，可以解出任意時(shí)刻的角度。

在某些情況下，雅可比矩陣的行列式可能趨近于零或等于零，意味著方程(4)多解或無解。多解可能的原因是：機(jī)械系統(tǒng)由于設(shè)計(jì)導(dǎo)致在某個(gè)特定位置本身就是多解的，如圖2中的曲柄?滑塊機(jī)構(gòu)，OA運(yùn)動帶動B運(yùn)動通常不會有問題；但反過來B運(yùn)動帶動OA運(yùn)動(如蒸汽機(jī))，當(dāng)OAB共線時(shí)，OA下一時(shí)刻可以向上也可以向下(考慮動力學(xué)時(shí)OA桿由于慣性繼續(xù)運(yùn)動，解是確定的，但僅考慮運(yùn)動學(xué)時(shí)OA桿位置是多解的)。無解的可能是：機(jī)械系統(tǒng)由于尺寸關(guān)系導(dǎo)致運(yùn)動發(fā)生沖突，例如圖2中l(wèi)

圖2 無解或多解的情況

編程計(jì)算得到角度的變化關(guān)系后，可以算出任意時(shí)刻各鉸的位置，以及BC桿上不同點(diǎn)的運(yùn)動軌跡(圖3)：很明顯B點(diǎn)軌跡是圓，C點(diǎn)軌跡是圓的一部分(AB桿大范圍運(yùn)動時(shí)，CD桿只在小范圍運(yùn)動)，而在BC桿上不同的點(diǎn)軌跡就很復(fù)雜了。

圖3 BC桿上不同點(diǎn)的運(yùn)動軌跡

知道任意時(shí)刻系統(tǒng)的各鉸位置，在屏幕上畫出各鉸點(diǎn)的位置并連接起來，就得到了四連桿機(jī)構(gòu)在某一時(shí)刻的圖象，延遲一定的時(shí)間后再畫出下一時(shí)刻的圖象，就形成了動畫。本問題中動畫的源代碼見圖4，其中plot函數(shù)表示畫線段；h1是句柄，定義動畫中每一幀圖的特性，如顏色、線型、寬度等；set函數(shù)是讀取數(shù)據(jù)；drawnow函數(shù)是畫出當(dāng)前幀的圖案；pause函數(shù)表示暫停，讓每一幀有一定的視覺殘留以便更好地形成動畫(類似放電影，1秒鐘放24格圖像)[3]。

圖4 動畫部分的源代碼截圖

2 函數(shù)求導(dǎo)獲得角速度、角加速度

求出角度后，對約束方程(1)求導(dǎo)出現(xiàn)角速度，有

由于θ,φ1,φ2已在前面求出，因此得到關(guān)于 ˙φ1,˙φ2的一組線性方程組。類似X=inv(A)*B可解出角速度，從而可以獲得角速度隨時(shí)間或隨θ變化的關(guān)系(圖5)。

對式(5)進(jìn)一步求導(dǎo)得到角加速度的方程

此時(shí)θ,φ1,φ2,˙φ1,˙φ2都是已知值，因此得到關(guān)于¨φ1,¨φ2的一組線性方程組。類似X=inv(A)*B可解出角加速度，從而可以獲得角加速度隨時(shí)間或隨θ變化的關(guān)系(圖6)。

這只是一類典型機(jī)構(gòu)的例子，采用本文介紹的方法，任意機(jī)構(gòu)的運(yùn)動軌跡、剛體的角速度和角加速度、某點(diǎn)的速度和加速度都可以類似求出。而這些結(jié)果(圖3、圖5、圖6)都是傳統(tǒng)方法難以獲得的。

圖5 角速度與θ的關(guān)系

圖6 角加速度與θ的關(guān)系

3 M atlab中的符號求導(dǎo)

如果運(yùn)動學(xué)問題中沒有具體數(shù)值，全是參數(shù)，用Matlab也可以處理。

案例2：半徑為r的圓輪在水平桌面上作直線純滾動，輪心速度v O的大小為常數(shù)。搖桿AB與桌面鉸接，并靠在圓輪上(圖7)。當(dāng)搖桿與桌面夾角φ=60°時(shí)，試求搖桿的角速度和角加速度。

圖7 一類接觸問題

傳統(tǒng)處理方法：取AB桿為動系，O為動點(diǎn)(圖8)。通過速度和加速度分析(具體過程略)，可以得到

圖8 速度分析

Mtlab處理本問題時(shí)，只需要把一般位置角度的表達(dá)式寫出來，然后利用符號求導(dǎo)的方法，就可以得到角速度和角加速度。角度以逆時(shí)針方向?yàn)檎?，設(shè)初始時(shí)刻AB桿垂直，根據(jù)圖9有

圖9 運(yùn)動學(xué)關(guān)系

根據(jù)式(8)有

對式(9)求一階導(dǎo)數(shù)獲得角速度，求二階導(dǎo)數(shù)獲得角加速度。程序的源代碼見圖10。

因此可以畫出系統(tǒng)運(yùn)動全程的角度、角速度、角加速度變化曲線，取歸一化參數(shù)r=1，v O=1，得到的曲線關(guān)系見圖12和圖13。

圖12 系統(tǒng)運(yùn)動參量與時(shí)間的關(guān)系

圖13 系統(tǒng)運(yùn)動參量與角度的關(guān)系

源代碼中具體涉及幾個(gè)函數(shù)：(1)符號定義是“syms”。(2)函數(shù)求導(dǎo)的格式是diff(function,’variable’)，表示函數(shù)對某個(gè)變量求導(dǎo)。運(yùn)行后得到的結(jié)果為

可以看到，傳統(tǒng)方法的分析要很大的篇幅，而Matlab只需要幾行就得到了結(jié)果(順便說一下，前面式(5)和式(6)也可以利用符號推導(dǎo)從式(1)得到)。

如果要求特定位置的值，只需要把具體值代入，格式為subs(function,old,new)，表示用新的變量替換老的變量，這部分程序的源代碼見圖11。

圖11 求特定位置的運(yùn)動量

運(yùn)行后屏幕顯示結(jié)果為

借助Matlab還可以對題目進(jìn)行深入研究。例如，用傳統(tǒng)方法分析時(shí)，如果動點(diǎn)選為圓盤上的接觸點(diǎn)，速度還可能會正確，但是加速度沒有辦法分析了，這是為什么呢？可以看看接觸點(diǎn)相對AB桿的軌跡，這借助計(jì)算機(jī)很容易實(shí)現(xiàn)。

以φ=60°為系統(tǒng)初始位置，設(shè)圓盤上P點(diǎn)與AB桿接觸，且注意負(fù)號表示順時(shí)針方向，因此計(jì)算機(jī)推導(dǎo)出的結(jié)果與傳統(tǒng)方法相同，見式(7)。

利用參數(shù)替換函數(shù)，還可以把式(1)和式(11)中的時(shí)間用角度替換(具體略)。當(dāng)圓心運(yùn)動Δx=v O t時(shí)，圓盤因?yàn)榧儩L動的轉(zhuǎn)角為θ=Δx/r(圖14)，且P點(diǎn)在動系中有

圖14 接觸點(diǎn)的軌跡分析

在動系中畫出的軌跡如圖15所示，很像是旋輪線。如果P點(diǎn)的軌跡是旋輪線，則曲線的尖點(diǎn)應(yīng)該接觸AB桿；而實(shí)際上P點(diǎn)接觸桿時(shí)，相對速度不為零，但軌跡的切線方向與桿平行，不過該曲線的曲率在詳細(xì)分析前還不清楚。

這就回答了前面的問題：如果選P點(diǎn)為動點(diǎn)，在分析速度時(shí)認(rèn)為相對速度沿桿方向(相信這樣做的學(xué)生自己也不能說得很清楚，只是憑感覺)，碰巧是正確的；但是在加速度分析時(shí)，由于相對運(yùn)動的軌跡不是太清楚，相對切向加速度未知，相對法向加速度也未知，再加上AB桿的角加速度未知，未知數(shù)太多求解不了。

注意在地面上看P點(diǎn)的軌跡是旋輪線，通過反轉(zhuǎn)和平移，旋輪線與圖15中的軌跡完全相同，這個(gè)結(jié)論是否出乎意料？

4 小結(jié)

在運(yùn)動學(xué)中引入Matlab，可以方便了解系統(tǒng)整個(gè)運(yùn)動過程，獲得系統(tǒng)中任意點(diǎn)的運(yùn)動軌跡，以及速度和加速度隨時(shí)間或某角度的變化規(guī)律，直觀明了，讓學(xué)生更容易理解系統(tǒng)的運(yùn)動。

本篇著重介紹了非線性方程組的求解方法，以及畫圖和動畫的格式；另外介紹了參數(shù)方程的求導(dǎo)和替換。畫圖涉及到的Matlab函數(shù)是plot(畫直線)，而句柄函數(shù)、set函數(shù)和drawon函數(shù)的組合就可以實(shí)現(xiàn)動畫；另一個(gè)重要函數(shù)是diff(函數(shù)求導(dǎo))和subs(變量替換)。掌握這幾個(gè)命令后，就可以對一般的運(yùn)動學(xué)問題進(jìn)行全程的計(jì)算和動畫演示，獲得豐富的運(yùn)動學(xué)信息。更進(jìn)一步，借助Matlab可以更加深入研究運(yùn)動學(xué)的問題，解決傳統(tǒng)分析方法無法處理的很多問題。