M atlab求解理論力學問題系列(二)典型機構的運動分析

高云峰

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

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

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

本文著重介紹2個典型機構的運動軌跡、速度和加速度分析，以及利用Matlab中的符號推導功能，從位移函數求導得到速度和加速度的表達式，并進行畫圖和動畫顯示。

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

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

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

圖1 四連桿機構

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

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

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

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

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

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

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

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

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

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

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

已知a1,a2,a3,a4和θ，估計值是計算雅可比矩陣，并代入當前估計值

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

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

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

圖2 無解或多解的情況

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

圖3 BC桿上不同點的運動軌跡

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

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

2 函數求導獲得角速度、角加速度

求出角度后，對約束方程(1)求導出現角速度，有

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

對式(5)進一步求導得到角加速度的方程

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

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

圖5 角速度與θ的關系

圖6 角加速度與θ的關系

3 M atlab中的符號求導

如果運動學問題中沒有具體數值，全是參數，用Matlab也可以處理。

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

圖7 一類接觸問題

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

圖8 速度分析

Mtlab處理本問題時，只需要把一般位置角度的表達式寫出來，然后利用符號求導的方法，就可以得到角速度和角加速度。角度以逆時針方向為正，設初始時刻AB桿垂直，根據圖9有

圖9 運動學關系

根據式(8)有

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

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

圖12 系統(tǒng)運動參量與時間的關系

圖13 系統(tǒng)運動參量與角度的關系

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

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

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

圖11 求特定位置的運動量

運行后屏幕顯示結果為

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

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

利用參數替換函數，還可以把式(1)和式(11)中的時間用角度替換(具體略)。當圓心運動Δx=v O t時，圓盤因為純滾動的轉角為θ=Δx/r(圖14)，且P點在動系中有

圖14 接觸點的軌跡分析

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

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

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

4 小結

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

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