自行車系統的計算機模擬

姜星宇

(同濟大學 電子與信息工程學院，上海 200092)

從自行車發明以來，其穩定的原理一直是人們關注的內容.R.S.Hand在其著作中對自行車穩定性做出了定義:若自行車能夠經過垂直其運動方向的擾動后，傾斜角θ與操縱轉角β漸進地趨于豎直平衡運動，就是穩定的[1,2].David E. H. Jones在1970年發表的The Stability of the Bicycle[3]中通過制作無法騎行的自行車過程中指出，在自行車的穩定性上，陀螺的穩定效應并不是影響自行車平衡的主要因素，至少不是全部因素. 而前輪的幾何因素對中心高低的影響對于穩定性也是十分重要的. 由于兩個非完整約束使自由度從無窮小時間內的三個變成有限時間內的七個，加之前輪與車架之間具有十分復雜的的非線性幾何關系，因此解出嚴格的解幾乎是不可能的. 許多研究便對于模型進行了線性的簡化或者是只考慮其中某些廣義坐標，或者半定量地討論. 例如通過只考慮平衡位置附近的微小擾動將約束方程與拉格朗日函數簡化，最后得到線性微分方程[2]. 亦有將前叉設置為與輪心連線垂直的簡化方法[4].

對于自行車穩定性的分析也具有一定的實際意義. 關于自行車的自動化也就是無人操控保持平衡有不少研究[4,5]，這也是未來將交通工具智能化的趨勢. 而要做到這一點，首先就要分析其平衡的條件與原理.

在不能將其力學原理完全精確地計算出，直接實驗又限于成本與時間時，近似求解數值解對于研究自行車的穩定性顯得十分重要.本文以拉格朗日方程為基本原理，從分析力學的角度對自行車系統近似計算與分析.

1 基本原理

(1)

將第一項展開也就是

(2)

其中cαi為坐標的函數，λα為約束力變量. 為使約束方程中出現二階導，將約束方程求導:

(3)

要注意的是，初始條件一定要滿足約束，如果不滿足，那么約束方程右端就會是一個由初始條件決定的非零常數，而這往往不是我們所需要的.

2 方法驗證

根據上述討論，我們將近似計算方法運用到任意存在約束的力學系統當中，求解拉格朗日方程. 對于不同的系統，計算過程唯一的區別就是函數的形式，因此很方便利用計算機求解. 而對于參數，根據力學相似性可以進行歸一化，因此下面的討論均不考慮系數與單位的問題. 先以傅科擺驗證方法的可行性.

在旋轉非慣性系下的球面擺是一個完整約束. 以x軸為東，y軸為北，慣性系中速度vabs=v+Ω×R0+Ω×r，那么拉格朗日函數為

圖1 北極處傅科擺x-y圖及周期(Δt=0.01 s)

表1 傅科擺周期隨緯度的變化

3 自行車系統的分析

3.1 輪的約束與動能

為了研究自行車系統的運動，我們先來分析一下輪的純滾動. 輪的純滾動通常是一個非完整約束，它滿足接觸點速度為零，即drC+dφ×r=0. 其中rC為圓心位矢，φ為輪的角速度，r為輪上一點相對于圓心的位矢.

現將輪簡化為厚度遠小于半徑，對稱的平面剛體. 它可以用質心坐標(x,y,z)與輪主軸歐拉角(θ,φ,ψ)來描述.
圖2顯示，θ對應的是輪與水平面的夾角，φ對應的是輪平面與水平面交線與x軸的夾角，而ψ為自轉角. 根據垂直軸定理，I3=2I1=I. 對應的動能為

(4)

輪的約束為vC+Ω×r0=0 (vC為質心速度，r0為接觸點相對圓心的位矢) ，3個分量為：

也就是說無窮小運動自由度共有3個.

圖2 輪的位置

3.2 自行車模型

對于自行車的研究主要是基于用最少的廣義坐標描述系統的狀態，將完整約束盡可能體現在計算時自由度的減少上. 這樣最少需要7個廣義坐標(這里模擬時為簡單起見并沒有將z=Rsinθ兩個完整約束直接代換因此使用了9個廣義坐標). 而對于動能與勢能的表示，需要用這些坐標表示出前后輪的坐標與歐拉角，也就是R.S.Hand在研究自行車系統運動時提到的“auxiliary variables”輔助坐標[2].

圖3是一個簡化過后的自行車的模型. 后輪的輪心為C1，前輪的輪心為C2. 前輪的轉軸在圖中以虛線表示，前后輪在同一平面上時(即前輪未轉向)轉軸與C1C2的交點為C3，夾角為α. 于是C1C3作為X軸與后輪法向量作為Z軸構成車架參考系.

圖3 自行車的位形

一個自行車系統具有兩個純滾動的輪，也就是說兩個輪在接觸處共有6個約束. 而在兩輪之間仍有約束，它們都是完整約束. 為了確定自行車的位形，在不考慮純滾約束的情況下可以使用后輪質心位置C1(x1,y1,z1)，車架的歐拉角θ、φ、ψ，前后輪的自轉角ψ1,ψ2，前輪相對于軸轉過的角度β共9個廣義坐標描述. 參數為輪的半徑R，質量m，繞軸轉動慣量I，前輪轉軸角度α，車架的長度C1C3=d與C3C2=δ.

為求出系統的力學函數與約束，還需要用這9個坐標和完整約束對應的幾何關系表示出前后輪坐標與歐拉角共12個坐標. 后輪在笛卡爾坐標系下的位置為

r1=(x1,y1,z1);θ1=θ;φ1=φ;ψ1=ψ1

(5)

而前輪的位置要通過車架參考系變換來求解. 現需要將OXYZ車架參考系繞Z旋轉-ψ，再繞X旋轉-θ，最后繞Z軸旋轉-φ得到Oxyz地面系. 定義旋轉矩陣：

(6)

對應上述旋轉，3個矩陣自右向左變換就可以通過乘上矩陣D將OXYZ系中的坐標變換為Oxyz系中的坐標. 而在車架參考系中的前輪中心C2坐標為

(7)

那么前輪在地面參考系中的位矢便為r2=r1+DR2. 前輪法向量n=D(-sinαsinβ,-cosαsinβ,cosβ)，那么歐拉角：

于是我們可以得到兩組坐標的變換：

Q=(x1,y1,z1,θ1,φ1,ψ1,x2,y2,z2,θ2,φ2,ψ2)T，

q=(x1,y1,z1,θ,φ,ψ,ψ1,ψ2,β)T

4 模擬結果

在這里我們可以使用第3節中提供的約束矩陣C=C′J與拉格朗日函數：

直接利用非完整約束下的拉格朗日方程近似求解.函數中T1與T2對應前輪與后輪的動能，可用3.1中的動能表達式得到.U為勢能函數，即mgz1+mgz2.

而參數我們可以用國際單位制表示現實生活中較為接近的自行車參數. 如下表.

表2 模擬參數(國際單位制)

4.1 速度對穩定性的影響

圖4是不同初始速度下在滿足能量守恒的時間范圍內θ與φ隨時間的變化. Δt代表近似計算所取的間隔，T為模擬的總時間.θ反映車身的傾角，而φ反映車的轉角. 可以看見，速度太小時系統會迅速失去穩定性.

對于可以穩定運行的運動，考慮其坐標變化的總體趨勢. 速度大意味著輪的角動量大，也就是說提供同樣的力矩情況下，速度大的φ轉角變化越慢. 從圖中也可以很明顯地看出. 定性地來講，就是角動量越大，越難以改變其運動狀態，同時也更加穩定. 考慮一種極限情況，就是速度足夠大，相當于重力足夠小. 這樣的話自行車可以保持其初始狀態不變，φ與θ都是恒定值. 這樣就可以合理地推測速度越大，θ與φ狀態改變地越慢.

圖4 θ、φ與t關系圖(Δt=0.01 s,T=10 s,δ=0)

4.2 “前輪尾跡”對穩定性的影響

Jones在著作中提及“前輪尾跡”對自行車平衡的重要性. “前輪尾跡”描述轉軸-地面交點與前輪接觸點的相對位置. 如果轉軸交點在車輪交點之前，那么當前輪偏離中央的位置時，地面的摩擦力與車架對前輪的力形成的力矩將與偏轉方向相反，反之則相同. 這與腳輪(比如超市中的購物車)的運動具有相似的特征[6]. 現設置不同的δ=C3C2來驗證一下“前輪尾跡”對平衡的影響.

由圖5可得，轉軸與地面的交點位于前輪接觸點前的兩個運動δ=0,-0.4 m還是相當穩定的，并且距離越大振蕩平緩得越快. 而另一個就不是那么穩定了. 不過也并非前輪設置的越靠后越好，由于那樣的車過于穩定從而轉彎等操作顯得并不靈活. 可以看出，模擬結果和Jones在其著作中的論斷是一致的，前輪過于靠前是無法騎行的.

圖5 θ與t關系圖(Δt=0.01 s，T=10 s,v=4 m/s)

4.3 陀螺效應以及前輪的自由度對穩定性的影響

模擬如圖6，這兩種情況下輪都向一個方向倒下，并且θ的變化情況相同. 這是由于在這兩種情況下，都不會由于輪的旋轉產生x方向的力矩，θ坐標對應的廣義力只由重力產生，最后只能向一個方向傾斜.

再來重點討論一下陀螺效應的影響. Jones認為陀螺效應的影響是十分小的，但是此處的仿真和簡單的推理卻顯示沒有角動量的自行車是無法運行的. 再次回顧Jones的實驗，可以發現他是將前輪的角動量抵消，而后輪的角動量仍然是存在的. 于是經過模擬只改變前輪或者后輪某一個輪的動能函數，我們得到如圖7的結果.

圖6 θ與t關系圖(Δt=0.01 s，T=0.87 s,v=4 m/s)

圖7 θ與t關系圖(Δt=0.01 s，T=0.87 s,v=16 m/s)

Jones制作的自行車“URB I”[3]即圖中虛線所代表的只抵消前輪角動量的情況，也并不是完全不穩定的. 加之其實驗中并沒有完全使軸對齊且平面重合，并不能說是完全抵消，人還是可以操控這樣的自行車的. 不過顯然這樣會比普通的車更加地難以操控，Jones也提到這樣的自行車無法在沒有人的操控下穩定運行. 從圖中也可以清晰地看見，當后輪的陀螺效應被消除后，自行車仍能在前輪角動量的調節下保持相當穩定. 也就是說前輪對于穩定性的貢獻是遠遠大于后輪的.

以上模擬了自行車自由運動下的一個簡單模型，而在實際情況下，由于騎行者操控反饋的存在，這個系統成為一個非封閉的系統，情況變得更加復雜. 不過自行車自由和被操控時運動的穩定性一定是緊密聯系的. 因此分析力學的方法可以在近似計算自行車這樣自由度較高的系統上有著比較簡潔的應用.