Euler's Equations for Rigid Body Rotation

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

L'esercizio simula la $\textbf{dinamica rotazionale di un corpo rigido}$ utilizzando le $\textbf{equazioni di Eulero}$ nel sistema solidale al corpo e ne integra numericamente l'evoluzione temporale.

Equazioni di Eulero

Per un corpo rigido con momenti di inerzia principali $I_1, I_2, I_3$, le equazioni di Eulero in presenza di momenti esterni $(N_1, N_2, N_3)$ sono: $$\large \begin{aligned} I_1 \dot{\omega}_1 &= (I_2 - I_3)\,\omega_2 \omega_3 + N_1, \\ I_2 \dot{\omega}_2 &= (I_3 - I_1)\,\omega_3 \omega_1 + N_2, \\ I_3 \dot{\omega}_3 &= (I_1 - I_2)\,\omega_1 \omega_2 + N_3. \end{aligned} $$

Nel codice queste equazioni vengono riscritte in forma esplicita come: $$\large \dot{\omega}_i = \frac{(I_j - I_k),\omega_j \omega_k}{I_i}

• \frac{N_i}{I_i},

\quad (i,j,k)\ \text{ciclica}. $$

Integrazione numerica

Le equazioni di Eulero vengono integrate numericamente tramite il $\textbf{metodo di Eulero esplicito}$, che approssima l'evoluzione temporale della velocità angolare come: $$\large \boldsymbol{\omega}(t + \Delta t) = \boldsymbol{\omega}(t)

\dot{\boldsymbol{\omega}}(t),\Delta t. $$

La simulazione utilizza:

• momenti di inerzia costanti $I_1 = 1$, $I_2 = 2$, $I_3 = 3$,
• una velocità angolare iniziale $\boldsymbol{\omega}(0) = (0.1,\,0.2,\,0.3)$,
• assenza di momenti esterni ($\mathbf{N} = \mathbf{0}$),
• un passo temporale $\Delta t = 0.01$,
• un tempo totale di simulazione $T = 10$.

Caso libero (assenza di coppie)

Poiché i momenti esterni sono nulli, il sistema descrive un $\textbf{corpo rigido libero}$. In questo caso:

• l'energia cinetica rotazionale è costante,
• il modulo del momento angolare è conservato,
• le componenti della velocità angolare possono oscillare a causa dell'anisotropia del tensore di inerzia.

Analisi dei risultati

Durante la simulazione viene calcolato il modulo della velocità angolare: $$\large |\boldsymbol{\omega}(t)| = \sqrt{\omega_1^2(t) + \omega_2^2(t) + \omega_3^2(t)}. $$

Il codice produce un grafico di $|\boldsymbol{\omega}|$ in funzione del tempo, che permette di verificare qualitativamente la conservazione dell'energia e la stabilità della rotazione nel caso libero.

Conclusione

L'esercizio mostra come:

• formulare le equazioni di Eulero per un corpo rigido,
• integrarle numericamente,
• analizzare l'evoluzione della velocità angolare,
• visualizzare il comportamento dinamico del sistema.

Si tratta di un esempio didattico di dinamica rotazionale non lineare di un corpo rigido tridimensionale.

import numpy as np
import matplotlib.pyplot as plt

def euler_equations_omega_dot(I1, I2, I3, omega, torque):
    """
    Calculate the derivatives of the angular velocities using
    Euler's equations.
        :param I1: Moment of inertia about axis 1.
        :param I2: Moment of inertia about axis 2.
        :param I3: Moment of inertia about axis 3.
        :param omega: Current angular velocities [omegal, omega2, omega3].
        :param torque: External torques [Nl, N2, N3].
    :return: 
        Derivatives of angular velocities [domegal/dt, domega2/dt, domega3/dt].
    """
    omega1, omega2, omega3 = omega
    N1, N2, N3 = torque
    domega1_dt = (I2 - I3) * omega2 * omega3 / I1 + N1 / I1
    domega2_dt = (I3 - I1) * omega3 * omega1 / I2 + N2 / I2
    domega3_dt = (I1 - I2) * omega1 * omega2 / I3 + N3 / I3

    return np.array([domega1_dt, domega2_dt, domega3_dt])

def integrate_euler_equations(I1, I2, I3, omega0, torques, dt, T):
    """
    Numerically integrate Euler's equations over time using Euler's
    method.
        :param I1: Moment of inertia about axis 1.
        :param I2: Moment of inertia about axis 2.
        :param I3: Moment of inertia about axis 3.
        :param omega0: Initial angular velocities [omegal, omega2, omega3].
    :param torques: External torques [N1, N2, N3].
    :param dt: Time step for integration.
    :param T: Total simulation time.
    :return: 
        Array with angular velocities over time.
    """
    n_steps = int(T / dt)
    omega = np.zeros((n_steps, 3))
    omega[0] = omega0
    for step in range(1, n_steps):
        omega_dot = euler_equations_omega_dot(I1, I2, I3, omega[step-1], torques)
        omega[step] = omega[step-1] + omega_dot * dt
    return omega

# Simulation parameters
I1, I2, I3 = 1.0, 2.0, 3.0 # Moments of inertia
omega0 = np.array([0.1, 0.2, 0.3]) # Initial angular velocities
torques = np.array([0.0, 0.0, 0.0]) # External torques
dt = 0.01 # Time step
T = 10.0 # Total simulation time 

# Compute angular velocities over time
omega_timeline = integrate_euler_equations(I1, I2, I3, omega0, torques, dt, T) 

# Print results
#for i, omega in enumerate(omega_timeline):
#    print(f"Time {i*dt:.2f} s: Omega = {omega}")

# Time array
time = np.linspace(0, T, len(omega_timeline))

# |omega|
omega_magnitude = np.linalg.norm(omega_timeline, axis=1)

# Plot |omega| vs time
plt.figure()
plt.figure(figsize=(4, 4))
plt.plot(time, omega_magnitude)
plt.xlabel("Time [s]")
plt.ylabel("|ω| [rad/s]")
plt.title("Angular velocity magnitude vs time")
plt.grid(True)
plt.show()

<Figure size 640x480 with 0 Axes>