陀螺的数值模拟（Matlab）

贡献者： addis

本文处于草稿阶段。

预备知识 1　刚体转动数值模拟（Matlab）

1. 进动和章动

图 1：运行结果，动画见这里。红色线段为角速度，浅蓝色线段为力矩。陀螺的上下摆动称为章动（nutation）。

　　该陀螺模型的惯性张量见习题 3 。

2. 一些解析分析

预备知识 2　刚体定点旋转的运动方程（欧拉角）

　　令陀螺底端为坐标原点，建立极坐标系，转轴的方向可以用极坐标 $\theta, \phi$ 表示，而绕轴转动的角度 $\psi$，那么 $\psi, \theta, \phi$ 可以看作 $z$-$y$-$z$ 欧拉角。角速度矢量为

\begin{equation} \boldsymbol{\mathbf{\omega}} = \dot\psi \hat{\boldsymbol{\mathbf{r}}} + \dot\theta \hat{\boldsymbol{\mathbf{\phi}}} + \dot\phi \hat{\boldsymbol{\mathbf{z}}} = (\dot\psi + \dot\phi\cos\theta) \hat{\boldsymbol{\mathbf{r}}} - \dot\phi\sin\theta\ \hat{\boldsymbol{\mathbf{\theta}}} + \dot\theta \hat{\boldsymbol{\mathbf{\phi}}} ~. \end{equation}

其中使用了 $ \hat{\boldsymbol{\mathbf{z}}} = \cos\theta\ \hat{\boldsymbol{\mathbf{r}}} - \sin\theta\ \hat{\boldsymbol{\mathbf{\theta}}} $。在 $ \hat{\boldsymbol{\mathbf{r}}} , \hat{\boldsymbol{\mathbf{\theta}}} , \hat{\boldsymbol{\mathbf{\phi}}} $ 的坐标系中，令

\begin{equation} \boldsymbol{\mathbf{I}} = \begin{pmatrix}I_1 & 0&0 \\ 0& I_2 &0 \\ 0& 0& I_2\end{pmatrix} , \qquad \boldsymbol{\mathbf{I}} ^{-1} = \begin{pmatrix}1/I_1 & 0&0 \\ 0& 1/I_2 &0 \\ 0& 0& 1/I_2\end{pmatrix} ~. \end{equation}

于是角动量为

\begin{equation} \boldsymbol{\mathbf{L}} = L_1(\dot\psi + \dot\phi\cos\theta) \hat{\boldsymbol{\mathbf{r}}} - L_2\dot\phi\sin\theta\ \hat{\boldsymbol{\mathbf{\theta}}} + L_2\dot\theta \hat{\boldsymbol{\mathbf{\phi}}} ~, \end{equation}

动能为（式 3 ）

\begin{equation} T = \frac{1}{2} \boldsymbol{\mathbf{\omega}} \boldsymbol\cdot \boldsymbol{\mathbf{L}} = \frac{L_1}{2}(\dot\psi + \dot\phi\cos\theta)^2 + \frac{L_2}{2}(\dot\phi^2\sin^2\theta + \dot\theta^2)~. \end{equation}

若质心离原点距离为 $a$，那么势能为

\begin{equation} V = Mga\cos\theta~, \end{equation}

这样就可以列出拉格朗日方程。

　　但与其列拉格朗日方程，倒不如直接用式 3 。重力力矩为

\begin{equation} \boldsymbol{\mathbf{\tau}} = Mga\sin\theta\ \hat{\boldsymbol{\mathbf{\phi}}} ~, \end{equation}

\begin{equation} \boldsymbol{\mathbf{\omega}} \boldsymbol\times \boldsymbol{\mathbf{L}} = -(L_2-L_1)(\dot\psi + \dot\phi\cos\theta)(\dot\theta\ \hat{\boldsymbol{\mathbf{\theta}}} + \dot\phi\sin\theta \hat{\boldsymbol{\mathbf{\phi}}} )~. \end{equation}

现在就可以代入转动方程

\begin{equation} \dot{ \boldsymbol{\mathbf{\omega}} } = \boldsymbol{\mathbf{I}} ^{-1}( \boldsymbol{\mathbf{\tau}} - \boldsymbol{\mathbf{\omega}} \boldsymbol\times \boldsymbol{\mathbf{L}} )~, \end{equation}

这样就可以解出 6 元常微分方程组。

3. 无章动的条件

图 2：第二套参数的运行结果，若适当给陀螺的初始角速度叠加一个绕 $z$ 轴的初始角速度，会发现章动消失，只剩下进动（precession）。

　　本节假设陀螺是轴对称的，在 $S'$ 系中相对转轴的角速度为 $ \boldsymbol{\mathbf{\omega}} _0$，而进动角速度，即 $S'$ 系相对于 $S$ 的角速度）为 $ \boldsymbol{\mathbf{\omega}} _a$（延 $z$ 方向)，那么总角速度为 $ \boldsymbol{\mathbf{\omega}} = \boldsymbol{\mathbf{\omega}} _0 + \boldsymbol{\mathbf{\omega}} _a$。显然此时 $ \boldsymbol{\mathbf{L}} = \boldsymbol{\mathbf{I}} \boldsymbol{\mathbf{\omega}} $ 的大小不变，处于转轴和 $z$ 轴所在的平面，且相对 $S'$ 系静止。所以角动量定理为

\begin{equation} M \boldsymbol{\mathbf{r}} _c \boldsymbol\times \boldsymbol{\mathbf{g}} = \boldsymbol{\mathbf{\tau}} = \dot{ \boldsymbol{\mathbf{L}} } = \boldsymbol{\mathbf{\omega}} _a \boldsymbol\times \boldsymbol{\mathbf{L}} = \boldsymbol{\mathbf{\omega}} _a \boldsymbol\times \boldsymbol{\mathbf{I}} ( \boldsymbol{\mathbf{\omega}} _0 + \boldsymbol{\mathbf{\omega}} _a)~. \end{equation}

对每个给定的 $ \boldsymbol{\mathbf{I}} $，$ \boldsymbol{\mathbf{\omega}} _0$ 和 $M, \boldsymbol{\mathbf{r}} _c$，该式存在唯一的解 $ \boldsymbol{\mathbf{\omega}} _a$。把初始条件中的角速度设为 $ \boldsymbol{\mathbf{\omega}} _0 + \boldsymbol{\mathbf{\omega}} _a$，则不会出现章动。

　　若陀螺转轴的极坐标的极角为 $\theta$，陀螺延转轴的转动惯量为 $I_1$，延另外两个垂直方向为 $I_2$，那么

\begin{equation} L_r = I_1(\omega_0+\omega_a\cos\theta), \qquad L_\theta = -I_2\omega_a\sin\theta, \qquad L = \sqrt{L_r^2 + L_\theta^2}~. \end{equation}

若 $ \boldsymbol{\mathbf{L}} $ 的极角为 $\theta_L$，则

\begin{equation} \theta_L - \theta = \arctan\frac{L_\theta}{L_r}~. \end{equation}

那么式 9 就是

\begin{equation} L\sin\theta_L \omega_a = Mga\sin\theta~, \end{equation}

解出 $\omega_a$ 即可。

未完成：在模拟中加入摩擦力（通过广义力实现），演示章动如何逐渐消失

未完成：使用四元数乘法 $qv\tilde q$ 旋转矢量而不是用旋转矩阵

未完成：使用体坐标系的欧拉方程而不是实验室坐标系

代码 1：top_sim.m

% 陀螺的数值计算
function top_sim
close all;
% === 参数设置 ===
M = 1; R0 = 1; a= 2; % 陀螺半径和长度
I0 = 2/5*M*R0^2*eye(3) + M*a^2*[1 0 0; 0 1 0; 0 0 0]; % 惯性张量
th0 = -pi/3;
q0 = [cos(th0/2); [1; 0; 0]*sin(th0/2)]; % 初始朝向（四元数）
w0 = quat2mat(q0)*[0; 0; 97.86]; % 初始角速度
g = 9.8; % 重力加速度
tau = @(q) cross(quat2mat(q)*[0;0;a], M*g*[0,0,-1]); % 力矩函数
tmin = 0; tmax = 12.2; Nt = 200; % 时间网格
RelTol = 1e-6; % 微分方程精度
% ===============

% 解微分方程
invI0 = inv(I0);
Y0 = [q0; w0];
f = @(t,Y)ode_fun(t, Y, I0, invI0, tau);
options = odeset('RelTol', RelTol);
[tt, Y] = ode45(f, [tmin,tmax], Y0, options);
t = linspace(tmin, tmax, Nt);
q = zeros(4, Nt); w = zeros(3, Nt);
for i = 1:4 % 四元数
    q(i,:) = interp1(tt, Y(:,i), t, 'spline');
end
for i = 1:3 % 角速度
    w(i,:) = interp1(tt, Y(:,i+4), t, 'spline');
end

% 验证角动量定理
verify(Y, tt, I0, tau);

% 播放动画
figure;
P1 = zeros(3, Nt);
for it = 1:Nt
    % 求位置
    R = quat2mat(q(:,it));
    clf; hold on; grid on; axis equal; view(33, 17);
    axis([-1, 1, -1, 1, -0.2, 1]*2*a);
    P1(:,it) = plot_top(R, R0, a);
    title(['t = ' num2str(t(it), '%.3f')]);
    % 画瞬时角速度
    plot3([0,w(1,it)],[0,w(2,it)],[0,w(3,it)], 'r');
    % 画力矩
    Tau = tau(q(:,it))/10;
    plot3([0,Tau(1)],[0,Tau(2)],[0,Tau(3)], 'c');
    % 画轨迹
    plot3(P1(1,1:it), P1(2,1:it), P1(3,1:it), 'm');
    xlabel x; ylabel y; zlabel z;
    drawnow();

    % 取消注释可将每一帧保存为 png 图片（当前目录下）
    saveas(gcf, [num2str(it) '.png']);
end
end

% 微分方程求导函数
function dY = ode_fun(~, Y, I0, invI0, tau)
q = Y(1:4); w = Y(5:7);
q = q / norm(q);
R = quat2mat(q);
RT = R';
Tau = tau(q); Tau = Tau(:);
dw = R*invI0*RT*(Tau - cross(w, R*I0*RT*w));
dq = 0.5*quat_mul([0; w], q);
dY = [dq; dw];
end

% 两个四元数相乘
function out = quat_mul(q1, q2)
s1 = q1(1); v1 = q1(2:4);
s2 = q2(1); v2 = q2(2:4);
out = [s1*s2 - dot(v1,v2); s1*v2 + s2*v1 + cross(v1, v2)];
end

% 由四元数 q 求旋转矩阵 R
function R = quat2mat(q)
s = q(1); vx = q(2); vy = q(3); vz = q(4);
R = [1 - 2*vy^2 - 2*vz^2, 2*vx*vy - 2*s*vz, 2*vx*vz + 2*s*vy;
    2*vx*vy + 2*s*vz, 1 - 2*vx^2 - 2*vz^2, 2*vy*vz - 2*s*vx;
    2*vx*vz - 2*s*vy, 2*vy*vz + 2*s*vx, 1 - 2*vx^2 - 2*vy^2];
end

% 画陀螺
function P1 = plot_top(R, R0, a)
[X,Y,Z] = sphere(10);
X = X*R0; Y = Y*R0; Z = Z*R0; Z = Z + a;
P = R * [X(:)'; Y(:)'; Z(:)'];
P1 = R * [0; 0; 2*a];
X = reshape(P(1,:), size(X));
Y = reshape(P(2,:), size(Y));
Z = reshape(P(3,:), size(Z));
surf(X, Y, Z, 'FaceColor', 'w', 'EdgeColor', 'b'); hold on;
plot3([0, P1(1)], [0, P1(2)], [0, P1(3)], 'LineWidth', 2, 'Color', 'k');
end

% 验证角动量定理
function verify(Y, t, I0, tau)
Nt = numel(t);
L = zeros(3, Nt); w = L; % 解的角动量
Ek = zeros(1, Nt); W = Ek;
q = Y(:, 1:4);
for it = 1:Nt
    w(:,it) = Y(it, 5:7);
    R = quat2mat(q(it,:));
    L(:, it) = R*I0*R'*w(:,it);
    Ek(it) = 0.5*dot(w(:,it), L(:, it));
    if (it > 1)
        W(it) = W(it-1) + dot(tau(q(it,:)), w(:,it))*(t(it)-t(it-1));
    else
        W(it) = Ek(it);
    end
end
L0 = zeros(3, Nt); % 冲量矩
L0(:, 1) = L(:, 1);
for it = 2:Nt
    L0(:, it) = L0(:, it-1) + tau(q(it,:)).' *(t(it)-t(it-1));
end

figure;
% 动能定理
subplot(5, 1, 1); ylabel E_k; hold on;
plot(t, Ek, 'r'); plot(t, W, 'b--');
% 角动量定理
for i = 1:3
    subplot(5, 1, 1+i);
    plot(t, L(i, :), 'r'); hold on;
    plot(t, L0(i, :), 'b--');
end
subplot(5, 1, 2); ylabel L_x;
subplot(5, 1, 3); ylabel L_y;
subplot(5, 1, 4); ylabel L_z;
% 角动量
subplot(5, 1, 5);
plot(t, w); hold on;
ylabel 'w_x, w_y, w_z';
legend({'w_x','w_y','w_z'});
xlabel t;
end

　　把以下参数替换到程序开始即可。

% === 参数设置 ===
M = 1; R0 = 1; a= 2; % 陀螺半径和长度
I0 = 2/5*M*R0^2*eye(3) + M*a^2*[1 0 0; 0 1 0; 0 0 0]; % 惯性张量
th0 = -pi/3;
q0 = [cos(th0/2); [1; 0; 0]*sin(th0/2)]; % 初始朝向（四元数）
w0 = quat2mat(q0)*[0; 0; 97.86] + [0; 0; 0.5]; % 初始角速度（矢量）
g = 9.8; % 重力加速度
tau = @(q) cross(quat2mat(q)*[0;0;a], M*g*[0,0,-1]); % 力矩函数
tmin = 0; tmax = 12.2; Nt = 200; % 时间网格
RelTol = 1e-6; % 微分方程精度
% ===============