彰化．山上．清水: 剛體轉動的數值模擬

Mar 15, 2014

剛體轉動的數值模擬

這邊介紹一個簡單的剛體轉動的數值模擬方法，這如果當教科書例題的話好像還不錯。我們以陀螺為例子，對陀螺的尤拉方程(Euler equation)的貼體角速度(body angular velocity)直接進行方向餘弦遞推數值積分。只要給定初始條件，即可追蹤模擬陀螺的完整運動。以此方法寫出的程式非常簡單，這邊以Python程式編寫。最後模擬了陀螺的四種經典運動：尖點，有環，無環，等周速運動。完整的物理理論證明也在我另一篇提供，因此這邊只寫下方法。

剛體轉動(Rigid body rotation)的尤拉方程如下，

重新整理成以下形式就可用簡單的四階Runge Kutta數值方法解出ω_b(t)，

一般認為解出貼體角速度ω_b(t)並無太大用處，因此一般書上在這裡會轉而求諸尤拉角(euler angle)並且用Lagrangian mechanics求出完整的解析解，不過完整的陀螺解析解包含橢圓積分(elliptic intergral)因此數值計算上相對我們這邊介紹得更為複雜。事實上若給定貼體角速度在空間座標上的初始值，我們可以以數值積分方法的初始值問題一步步求出剛體的貼體座標在空間座標中的運動，雖然我們沒有得到完整的解析解，我們還是可以完整模擬陀螺的三維運動。(此方法屬於方向餘弦矩陣遞推法，其物理證明在另一篇完整的證明中提供。)

我們考慮陀螺的body對稱單位軸z軸，假設z軸在空間中的初始位置為z(t₀)，我們以此為出發點，計算z(t₀)沿著ω_b(t₁)向量旋轉|ω_b(t₁)|*dt角度後，轉動後新的z軸為z(t₁)，接著再沿著ω_b(t₂)向量旋轉|ω_b(t₂)|*dt角度後轉到z(t₂)，以此可以算出z(t_1~N)。(這邊當然需要考慮微小轉動向量的特性infinitesimal rotation)

以下是python模擬出的四種經典運動的locus圖(body z軸跟單位圓的交叉點的軌跡圖)，參數也有寫在圖中。