//最近在看理力，隨便寫(xiě)點(diǎn)筆記。

//一個(gè)挺有意思的剛體現(xiàn)象，順手寫(xiě)了個(gè)數(shù)值模擬。

//和封面視頻配合觀看效果更佳。

01 起源

Dzhanibekov效應(yīng)是蘇聯(lián)宇航員Dzhanibekov于1985年在太空中執(zhí)行任務(wù)是偶然發(fā)現(xiàn)的：形狀不規(guī)則的螺母在自由轉(zhuǎn)動(dòng)時(shí)，會(huì)出現(xiàn)周期性的翻轉(zhuǎn)：

出現(xiàn)這種現(xiàn)象的原因是不規(guī)則物體的轉(zhuǎn)動(dòng)角速度和角動(dòng)量的方向有微小差別，導(dǎo)致角速度出現(xiàn)復(fù)雜的變化。

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

篇幅所限，我們只簡(jiǎn)述該問(wèn)題的理論模型及相關(guān)微分方程。詳見(jiàn)本文的封面視頻。想深入了解相關(guān)理論可參閱相關(guān)教材。[1]

為了準(zhǔn)確描述剛體角動(dòng)量和角速度的關(guān)系，我們不妨取固定在剛體上的坐標(biāo)系，并合理選取方向使之為剛體主軸，這意味著剛體的角動(dòng)量和角速度的關(guān)系可以簡(jiǎn)單表示為：

注意到剛體時(shí)刻在轉(zhuǎn)動(dòng)，固連在剛體上的坐標(biāo)系也在轉(zhuǎn)動(dòng)，故

而自由轉(zhuǎn)動(dòng)的剛體滿(mǎn)足

展開(kāi)以上方程，得到

由該方程組和初始角速度得到任意時(shí)刻剛體繞固連自身軸的角速度。

接下來(lái)用歐拉角描述固連剛體上的坐標(biāo)系和剛體質(zhì)心慣性系的位置關(guān)系，如圖所示：

有歐拉運(yùn)動(dòng)學(xué)方程聯(lián)系角速度與歐拉角：

由此給出三個(gè)歐拉角與時(shí)間的關(guān)系，然后即可導(dǎo)出慣性系與剛體系之間的變換矩陣。

03 代碼

我們使用mathematica對(duì)該問(wèn)題進(jìn)行數(shù)值模擬與動(dòng)畫(huà)繪制。

首先了解一下長(zhǎng)方體的繪制方法。我們將用長(zhǎng)寬高4:3:1的長(zhǎng)方體進(jìn)行模擬。

Graphics3D[Cuboid[{-4, -3, -1}, {4, 3, 1}], 
 PlotRange -> {{-7, 7}, {-7, 7}, {-7, 7}}]

注意到角速度的微分方程對(duì)轉(zhuǎn)動(dòng)慣量是齊次的，我們不妨設(shè)轉(zhuǎn)動(dòng)慣量為（不難證明長(zhǎng)方體繞對(duì)稱(chēng)軸的轉(zhuǎn)動(dòng)慣量正比于另兩個(gè)軸尺寸的平方和）：

Ix = 10;
Iy = 17;
Iz = 25;

角速度微分方程及相應(yīng)的初值：

Equations = {Ix omegax'[t] + (Iz - Iy) omegay[t] omegaz[t] == 0,
 ? Iy omegay'[t] + (Ix - Iz) omegaz[t] omegax[t] == 0,
 ? Iz omegaz'[t] + (Iy - Ix) omegax[t] omegay[t] == 0};
Initial = {omegax[0] == 0.001, omegay[0] == 2, omegaz[0] == 0};

Dzhanibekov定理指出，剛體繞轉(zhuǎn)動(dòng)慣量最小/最大的主軸轉(zhuǎn)動(dòng)較為穩(wěn)定，而繞轉(zhuǎn)動(dòng)慣量中等的主軸轉(zhuǎn)動(dòng)是不穩(wěn)定的。為了觀察到翻轉(zhuǎn)現(xiàn)象，這里指定剛體繞y軸轉(zhuǎn)動(dòng)，omegax作為微擾。

數(shù)值求解微分方程并定義函數(shù)存儲(chǔ)結(jié)果(以下只計(jì)算了前20秒)：

result = NDSolve[{Equations, Initial}, {omegax[t], omegay[t], 
 ? omegaz[t]}, {t, 0, 40}]
\[Omega]1[t_] = First[omegax[t] /. result];
\[Omega]2[t_] = First[omegay[t] /. result];
\[Omega]3[t_] = First[omegaz[t] /. result];

運(yùn)行界面如下。

同樣的方法求解歐拉角與時(shí)間的關(guān)系，并定義函數(shù)存儲(chǔ)結(jié)果：

eqns2 = {\[Theta]'[t] Cos[\[Psi][t]] + \[Phi]'[
 ? ? ? t] Sin[\[Theta][t]] Sin[\[Psi][t]] == \[Omega]1[t],
 ? -\[Theta]'[t] Sin[\[Psi][t]] + \[Phi]'[
 ? ? ? t] Sin[\[Theta][t]] Cos[\[Psi][t]] == \[Omega]2[t],
 ? \[Psi]'[t] + \[Phi]'[t] Cos[\[Theta][t]] == \[Omega]3[t]};
initial2 = {\[Theta][0] == 0.001, \[Phi][0] == 0.001, \[Psi][0] == 
 ? ?0.001};
result2 = 
 NDSolve[{eqns2, initial2}, {\[Theta][t], \[Phi][t], \[Psi][t]}, {t, 
 ? 0, 40}]
\[Phi]1[t_] = First[\[Phi][t] /. result2];
\[Theta]1[t_] = First[\[Theta][t] /. result2];
\[Psi]1[t_] = First[\[Psi][t] /. result2];

(b站專(zhuān)欄系統(tǒng)里的代碼塊功能沒(méi)有包含mathematica語(yǔ)法，導(dǎo)致上面這塊代碼高亮標(biāo)得亂七八糟，別管他就是了)

這里有一個(gè)小細(xì)節(jié)，initial2也即初始?xì)W拉角設(shè)置的時(shí)候，本來(lái)是該設(shè)為(0,0,0)的，但由于歐拉運(yùn)動(dòng)學(xué)方程中含有歐拉角正弦項(xiàng)，若在該初始值下開(kāi)始數(shù)值求解，則會(huì)在t=0時(shí)出現(xiàn)0=0，導(dǎo)致Mathematica無(wú)法求解。將角度設(shè)為0.0001，就解決了這個(gè)問(wèn)題。

接下來(lái)全部結(jié)果都有了，可以繪制動(dòng)畫(huà)了。我們利用幾何變換函數(shù)和歐拉旋轉(zhuǎn)矩陣函數(shù)將一開(kāi)始的長(zhǎng)方體旋轉(zhuǎn)到正確位置，繪制動(dòng)畫(huà)的代碼如下：

Animate[Graphics3D[
 ?GeometricTransformation[
 ? Cuboid[{-4, -3, -1}, {4, 3, 1}],
 ? EulerMatrix[{\[Phi]1[t], \[Theta]1[t], \[Psi]1[t]}, {3,1,3}]
 ? ]
 ?, PlotRange -> {{-7, 7}, {-7, 7}, {-7, 7}}],
 {t, 0, 40}]

轉(zhuǎn)起來(lái)了！

動(dòng)畫(huà)的具體效果見(jiàn)封面視頻。

參考文獻(xiàn)

[1] 周衍柏. 理論力學(xué)教程（第四版）[M]. 北京：高等教育出版社，2018.8：112~150.