高云峰

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

受力分析是理論力學(xué)教學(xué)中的重要內(nèi)容，而多點摩擦問題是其中的難點。在目前的教學(xué)安排中，靜力學(xué)有少量求解多點摩擦的習(xí)題，學(xué)生可以畫圖、列出方程并求解；但是在動力學(xué)中，幾乎沒有涉及多點摩擦的習(xí)題，而且也不要求解方程。

一個實際問題如何簡化為力學(xué)模型、如何分析列出方程、如何求解、如何驗證，是工科學(xué)生必需掌握的基本技能。當(dāng)然一門課程不可能教會學(xué)生處理實際問題的所有環(huán)節(jié)或流程，但在適當(dāng)?shù)那闆r下讓學(xué)生了解這一過程，讓部分感興趣的學(xué)生深入探究，還是有可能的。因此我在教學(xué)中有意識地設(shè)計了一些教具，讓學(xué)生有機會分析處理一個“完整問題”-- 問題的簡化、分析、計算、制作模型驗證。

1 實際問題：一個小巧的力學(xué)演示裝置

這是一個專門設(shè)計的裝置(圖1)，利用激光切割后，加上減速電動機、電池和開關(guān)，就可以拼裝出一個爬行裝置，且裝置尾巴長度可調(diào)(圖2)。

圖1 爬行裝置的設(shè)計圖

圖2 實際的爬行裝置

為了突出“完整問題” 的分析過程而不陷入復(fù)雜的計算，可以通過調(diào)整重心位置讓裝置在平面內(nèi)運動。觀察發(fā)現(xiàn)，裝置的運動狀況依賴于它的尾長：短尾時輪子一部分時間打滑，一部分時間不打滑，裝置前進得比較慢；長尾時輪子不打滑，前進得稍快些。

定性解釋很容易：尾長影響重心相對位置，導(dǎo)致尾部和輪子處的壓力不同。尾短時重心靠近尾部，輪子處壓力小會打滑；尾長時重心靠近輪子，輪子處壓力大不會打滑。

但是具體分析計算就復(fù)雜得多，輪子旋轉(zhuǎn)時，輪子和尾部受到的壓力、摩擦力與裝置的加速度和角加速度有關(guān)，不容易直觀判斷哪一點會打滑。為了進行具體的分析計算，首先就需要建立一個力學(xué)模型。

2 力學(xué)模型：如何合理簡化？

設(shè)OXY是地面參考坐標(biāo)系，Cxy是隨體坐標(biāo)系，C是裝置主體(不包輪子)的質(zhì)心，任意時刻，尾巴A和輪子邊緣B接觸地面。裝置的幾何尺寸參數(shù)為：尾巴長度a1；質(zhì)心距離尾部a2，距離底部h1，軸心O1距離尾部a2+a3，距離底部h2；輪長為l(圖3)。

圖3 裝置的力學(xué)模型

裝置分為兩部分，主體和輪子。設(shè)主體質(zhì)量為m1，輪子質(zhì)量為m2；主體質(zhì)心C 的坐標(biāo)為xC，yC；車身俯仰角度為α (繞z 軸)；觀察發(fā)現(xiàn)車輪以勻角速度ω0順時針轉(zhuǎn)動，設(shè)t 時刻θ =π/3-ω0t，因此θ ∈(-π/3，π/3) 為某輪接觸地面時的角度范圍。由于A、B 始終接觸地面，可以得到α 與輪子轉(zhuǎn)角θ的關(guān)系，根據(jù)軸心O1距離水平面的高度，有

從式(1) 解出α，求導(dǎo)后α、? 也都已知。類似可以得到y(tǒng)C， yO1與α 或θ 的關(guān)系。

為了保證裝置能正常運動，輪子長度需要滿足2 個條件：(1) l cos(π/3)＞h2，保證輪子總能接觸地面；(2)輪子不能太長，保證重心總落在A、B 之間，否則裝置會翻倒。

3 數(shù)學(xué)模型

裝置在平面內(nèi)運動時， 涉及的運動參數(shù)有xC，yC，α，力參數(shù)有FA，NA，F(xiàn)B， NB。注意到裝置中A、B 兩點接觸地面，哪點會打滑是多點摩擦問題中的難點。通常假設(shè)某點不打滑，然后得到計算結(jié)果，再判斷結(jié)果是否合理。在本問題中，存在三種模式。

3.1 模式一(假設(shè)A 不動，B 打滑)

設(shè)t 時刻A 點坐標(biāo)xA(t) = x*A，yA= 0 已知，根據(jù)圖3，其他各點的坐標(biāo)為

相應(yīng)各點速度、加速度都可以求出，從而各點的位置、速度、加速度都是已知函數(shù)。

根據(jù)受力圖(圖4)，利用動靜法，可以列出X和Y 方向力的平衡方程、B 點打滑的摩擦力補充方程、對A 點力矩的平衡方程，為

圖4 裝置的受力分析

這時方程(5) 中所有的運動學(xué)參數(shù)都是已知量，因此是代數(shù)方程，容易解出。

解出后要驗證是否滿足(1) 所有壓力、摩擦力均大于0；(2) FA≤μNA是否成立。如果條件都滿足，說明模式一假設(shè)成立(A 不動，B 打滑)，然后對t+Δt 時刻進行計算；否則說明模式一不成立，轉(zhuǎn)入模式二。

3.2 模式二(假設(shè)B 不動，A 打滑)

設(shè)t 時刻B 點坐標(biāo)xB(t) = x*B，yB= 0 已知，根據(jù)圖3，其他各點的坐標(biāo)為相應(yīng)各點的速度、加速度都可以求出，從而各點的位置、速度、加速度都是已知函數(shù)。

根據(jù)受力圖(圖4)，利用動靜法，系統(tǒng)方程為

方程(9)也是代數(shù)方程，解出后要驗證是否滿足(1) 所有壓力、摩擦力均大于0；(2) FB≤μNB是否成立。如果條件都滿足，說明假設(shè)模式二成立(B 不動，A 打滑)，然后對t+Δt 時刻進行計算；否則說明模式二不成立，轉(zhuǎn)入模式三。

3.3 模式三(假設(shè)A、B 均打滑)

方程(10) 中f1， f2， f3是α， α， ? 的函數(shù)，5 個方程5 個未知數(shù)，可以求解。求解結(jié)束后要驗證是否滿足(1)所有壓力、摩擦力均大于0；(2)FA≤μNA，F(xiàn)B≤μNB是否成立。如果條件都成立說明假設(shè)模式三成立(A 和B 都打滑)，然后對t+Δt 時刻進行計算；否則說明公式或程序中可能存在錯誤，要返回檢查。

但是如何求解微分代數(shù)方程，需要一些技巧，在下面數(shù)值計算中會介紹。

4 數(shù)值計算

有了方程，可以算出裝置運動和受力隨時間變化規(guī)律，然后把數(shù)據(jù)進行可視化處理。數(shù)值求解前通常先要畫出計算框圖，考慮各種可能情況(圖5)。下面采用Matlab 程序進行計算，可以在計算結(jié)束后直觀從屏幕上查看虛擬裝置如何運動，方便與實際裝置進行比較。

Matlab 程序采用矩陣和列陣計算很方便，可以把前面的方程統(tǒng)一改寫為AX =B 的形式，例如方程(10) 改寫為

圖5 求解方程的計算框圖

在模式一、模式二中，調(diào)用X = inv(A)*B 的命令格式就完成了求解。但模式三的微分代數(shù)方程處理起來要復(fù)雜一些：先把C當(dāng)作代數(shù)量，用X = inv(A)*B 求解出來，設(shè)解出C=*C(這時*C是具體數(shù)值)，再把C=*C當(dāng)作微分方程，化為標(biāo)準(zhǔn)的一階微分方程組后，直接調(diào)用榮格庫塔法(Runge-Kutta methods) 求解常微分方程。

具體過程是：設(shè)t 時刻系統(tǒng)各個參數(shù)(位置、速度) 都已經(jīng)算出，現(xiàn)在要通過微分方程C=*C計算出t+Δt 時刻的參數(shù)。設(shè)y1=xC，y2= ˙y1，這樣就把二階微分方程化為標(biāo)準(zhǔn)的一階微分方程組以及相應(yīng)的初始條件，為

積分時間段為[t，t+Δt]，求解常微分方程調(diào)用的格式為

程序式(13) 中參數(shù)按順序含義為：“t” 是存放積分時間；“iy” 是存放積分結(jié)果；“ode45” 是求解常微分方程的標(biāo)準(zhǔn)函數(shù)；“rg_kt” 是微分方程表達式的子函數(shù)；“[time， time+hh]” 是積分時間段；“y” 是初始值；“options” 是積分的誤差選項。

值得說明的是，考慮到數(shù)值計算存在誤差，因此在判斷摩擦關(guān)系是否違約時要稍微放松一點要求，假設(shè)積分允許誤差為ε=1.0×10-6，則把摩擦關(guān)系F ≤μN 放松為

數(shù)值計算得到的結(jié)果符合定性分析的結(jié)論，但是有更多細節(jié)。 下面是裝置在一個周期內(nèi) (θ ∈[-π/3，π/3]) 各參數(shù)變化的情況。

圖 6 表明了尾長與各點打滑時長的關(guān)系：(1) 尾巴越長，A 點打滑時間越多；存在一個臨界尾長a*1，a1≥a*1時A 點始終打滑。(2) 大部分情況下裝置不會出現(xiàn)A 和B 都打滑的情況，只在臨界尾長附近才會出現(xiàn)，且存在的時間都很短，在實際中不容易觀察到。

圖7 顯示了a1=5 mm 時主體質(zhì)心加速度的變化，一方面看出模式轉(zhuǎn)換時加速度會有跳變，另一方面可以看出質(zhì)心加速度比重力加速度小3 個數(shù)量級，因此可以認為裝置運動是“準(zhǔn)靜態(tài)”過程，從而可以從靜力學(xué)角度近似估計出a*1：如果輪子在x 方向距離質(zhì)心最遠時有NB≥NA，則A 點會一直要打滑。分別對圖4 中的A 和B 取矩，有

圖6 尾長與打滑時長的關(guān)系

考慮到實際上α 是小量，從式(15) 中近似得到

代入裝置的數(shù)據(jù)，得到a*1≈39 mm，接近圖6 中的40 mm。

圖7 尾長與質(zhì)心加速度的關(guān)系

為什么在臨界尾長a*1附近才會出現(xiàn)A 和B 均打滑的情況? 從圖8 和圖9 中可以看出，當(dāng)A 和B壓力曲線有交點時，會出現(xiàn)模式轉(zhuǎn)換。當(dāng)壓力曲線交點區(qū)域較窄時(斜率大)，只在模式一和模式二間進行轉(zhuǎn)換，不會出現(xiàn)模式三；當(dāng)壓力曲線交點區(qū)域較寬時(斜率小)，三種模式都會出現(xiàn)。這可以理解為：臨界尾長a*1的裝置在t →0 時，A 和B 壓力就很接近，此時質(zhì)心速度很小，數(shù)值計算時不能一步跨出這個區(qū)域，就會出現(xiàn)A 和B 都打滑這種模式。在非臨界尾長情況下，若A 和B 壓力接近，此時質(zhì)心速度較大，數(shù)值計算時可以一步跨出這個區(qū)域，不會出現(xiàn)A 和B 都打滑的狀況。

表1 中的數(shù)據(jù)支持了A 和B 都打滑的必要條件：剛開始的時刻，且兩者的壓力、摩擦力都很接近，以及在這種狀況下如果按模式一或模式二計算，都會出現(xiàn)摩擦關(guān)系違約的情況(a1= 36 mm，t =0.081 4 s，θ =53.939 4°)。

圖10 表明一個周期內(nèi)同一尾長A 點的壓力隨時間變小，因為時間增加時B 點更靠近重心；而同一時刻尾巴越長A 點壓力越小。圖11 表明通過判斷摩擦是否違約后再選擇適當(dāng)?shù)哪Ｊ剑煌查L時A 點摩擦關(guān)系都沒有出現(xiàn)違約情況，B 點也類似。

圖8 兩種模式的情況

圖9 三種模式的情況

表1 不同模式計算出的摩擦力與壓力

圖10 不同尾長時FA 和NA

圖11 不同尾長時FA/NA

5 小結(jié)

實驗、理論分析和數(shù)值計算都表明：爬行裝置的運動模式與尾長有密切關(guān)系。且理論分析、數(shù)值計算的結(jié)果都與實際情況吻合。

(1) 本案例是一個“完整問題”，包括設(shè)計、制作、建模、理論分析、數(shù)值計算?？梢宰寣W(xué)生更好地把多種知識和能力進行整合。

(2) 裝置尾短時輪子會打滑，前進距離較小；尾長時輪子不打滑，前進距離大。存在三種打滑的模式，根據(jù)摩擦關(guān)系是否違約來選擇某一種模式。

(3) 只有在特定尾長時才可能出現(xiàn)A 和B 都打滑的模式。這種模式的時長很短，在實際中不容易觀察到。

(4) 當(dāng)A 和B 中有一點靜止時，求解代數(shù)方程就可以得到位置和受力的結(jié)果；如果A 和B 都打滑，需求解微分代數(shù)方程。