Resolución de la ecuaciones de Euler para un Sólido Rígido

Tenemos 5 cubos idénticos, cada uno de lado 𝑑, masa 𝑀. Se sueldan formando una T: 3 cubos en horizontal (centrados en la parte superior), 2 cubos en vertical desde el cubo central.

El centro de masas se encuentra sobre el eje vertical de simetría, a una altura $\frac{3d}{5}$ por debajo del centro del cubo central de la base.

Respecto al CDM, la triz de inercia es:

$ \mathbf{I}_{CM} = \begin{pmatrix} \frac{121}{30}Md^2 & 0 & 0 \\ 0 & \frac{17}{6}Md^2 & 0 \\ 0 & 0 & \frac{103}{15}Md^2 \end{pmatrix} $

donde los ejes principaleson el longitudinal paralelo al lado vertical, el transversal paralelo al lado horizontal, y el perpendicular al plano de la T (todos ellos pasando por su CDM).

Generamos ahora la matriz de inercia como objeto de SymPy.

Ecuaciones de Euler para un Sólido Rígido

Las ecuaciones de Euler describen el movimiento de rotación de un sólido rígido alrededor de su centro de masas, en un sistema de referencia ligado al cuerpo (es decir, que gira con él). Estas ecuaciones son fundamentales en la dinámica de cuerpos rígidos y se expresan en términos de los momentos principales de inercia y las componentes de la velocidad angular.

Sean $( I_1, I_2, I_3)$ los momentos principales de inercia del cuerpo respecto a sus ejes principales, y $ (\omega_1, \omega_2, \omega_3) $ las componentes de la velocidad angular $( \boldsymbol{\omega} )$ en dichos ejes. Las ecuaciones de Euler son:

\ \begin{aligned} I_1 \frac{d\omega_1}{dt} + (I_3 - I_2)\omega_2 \omega_3 &= M_1 \\ I_2 \frac{d\omega_2}{dt} + (I_1 - I_3)\omega_3 \omega_1 &= M_2 \\ I_3 \frac{d\omega_3}{dt} + (I_2 - I_1)\omega_1 \omega_2 &= M_3 \end{aligned}

donde $( M_1, M_2, M_3 )$ son las componentes del momento externo (torque) aplicado al cuerpo respecto a los mismos ejes.

Estas ecuaciones se derivan de la segunda ley de Newton para la rotación,

$ \mathbf{M} = \frac{d\mathbf{L}}{dt} $,

considerando que el sistema de referencia está en rotación con el cuerpo. La aparición de los productos cruzados entre las componentes de $\boldsymbol{\omega} $ se debe al efecto giroscópico, que es crucial en muchos sistemas físicos, como satélites, giróscopos o vehículos en rotación.

En el caso especial en que no hay torques externos $( \mathbf{M} = \mathbf{0} )$, el sistema describe la dinámica libre del sólido rígido, donde pueden aparecer fenómenos como la inestabilidad de Dzhanibekov.

Ecuaciones del movimiento

Escribos las ecuaciones de Euler utilizando SimPy.

import sympy as sp

# Definir variables
I1, I2, I3 = sp.symbols('I1 I2 I3')  # Momentos de inercia
M1, M2, M3 = sp.symbols('M1 M2 M3')  # Momentos externos
omega1, omega2, omega3 = sp.symbols('omega1 omega2 omega3', cls=sp.Function)
t = sp.Symbol('t')  # Variable de tiempo

# Definir funciones de velocidad angular
omega1 = omega1(t)
omega2 = omega2(t)
omega3 = omega3(t)

# Definir ecuaciones de Euler
omega1_dot = sp.diff(omega1, t)
omega2_dot = sp.diff(omega2, t)
omega3_dot = sp.diff(omega3, t)

Euler1 = sp.Eq(I1 * omega1_dot + (I3 - I2) * omega2 * omega3, M1)
Euler2 = sp.Eq(I2 * omega2_dot + (I1 - I3) * omega3 * omega1, M2)
Euler3 = sp.Eq(I3 * omega3_dot + (I2 - I1) * omega1 * omega2, M3)

# Mostrar ecuaciones
display(Euler1, Euler2, Euler3)

$\displaystyle I_{1} \frac{d}{d t} \omega_{1}{\left(t \right)} + \left(- I_{2} + I_{3}\right) \omega_{2}{\left(t \right)} \omega_{3}{\left(t \right)} = M_{1}$

$\displaystyle I_{2} \frac{d}{d t} \omega_{2}{\left(t \right)} + \left(I_{1} - I_{3}\right) \omega_{1}{\left(t \right)} \omega_{3}{\left(t \right)} = M_{2}$

$\displaystyle I_{3} \frac{d}{d t} \omega_{3}{\left(t \right)} + \left(- I_{1} + I_{2}\right) \omega_{1}{\left(t \right)} \omega_{2}{\left(t \right)} = M_{3}$

Las únicas reacciones son las que aparecen en el punto de apoyo, que no pueden general momento. Por lo tanto, imponemos las condiciones

$M1=M2=M3=0 $

Además, los momentos de inercia los hemos definido previamente. Escribimos las ecuaciones y les damos nombre, para poder resolverlas más adelante.

Euler1 = sp.Eq(I1 * omega1_dot + (I3 - I2) * omega2 * omega3, 0)
Euler2 = sp.Eq(I2 * omega2_dot + (I1 - I3) * omega3 * omega1, 0)
Euler3 = sp.Eq(I3 * omega3_dot + (I2 - I1) * omega3 * omega1, 0)

# Mostrar ecuaciones
display(Euler1, Euler2, Euler3)

$\displaystyle I_{1} \frac{d}{d t} \omega_{1}{\left(t \right)} + \left(- I_{2} + I_{3}\right) \omega_{2}{\left(t \right)} \omega_{3}{\left(t \right)} = 0$

$\displaystyle I_{2} \frac{d}{d t} \omega_{2}{\left(t \right)} + \left(I_{1} - I_{3}\right) \omega_{1}{\left(t \right)} \omega_{3}{\left(t \right)} = 0$

$\displaystyle I_{3} \frac{d}{d t} \omega_{3}{\left(t \right)} + \left(- I_{1} + I_{2}\right) \omega_{1}{\left(t \right)} \omega_{3}{\left(t \right)} = 0$

Para resolver el sistema de ecuaciones anterior de forma analítica, usamos la función "dsolve" de SimPy. Como argumentos de entrada, tenemos que meter las ecuaciones que hemos escrito anteriormente.

El sistema es integrable porque conserva la energía cinética y el módulo del momento angular. Esto nos ayudará a reducir la complejidad de las ecuaciones diferenciales.

• Energía cinética

$$ T = \frac{1}{2}(I_1 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2) $$

• Momento angular total

$$ L^2 = (I_1 \omega_1)^2 + (I_2 \omega_2)^2 + (I_3 \omega_3)^2 $$

# Definición en SymPy
T = sp.Rational(1, 2)*(I1*omega1**2 + I2*omega2**2 + I3*omega3**2)
L2 = (I1*omega1)**2 + (I2*omega2)**2 + (I3*omega3)**2

Usamos las constantes para eliminar variables.

omega2_sq, omega3_sq = sp.symbols('omega2_sq omega3_sq')

eq_energy = sp.Eq(I1*omega1**2 + I2*omega2_sq + I3*omega3_sq, 2*T)
eq_momentum = sp.Eq((I1*omega1)**2 + (I2**2)*omega2_sq + (I3**2)*omega3_sq, L2)

sol = sp.solve([eq_energy, eq_momentum], (omega2_sq, omega3_sq))
sol

{omega2_sq: omega2(t)**2, omega3_sq: omega3(t)**2}

Sustituyendo en la ecuación de $\dot{\omega}_1$, obtenemos una ecuación de la forma:

$$ \dot{\omega}_1^2 = f(\omega_1) $$

lo que implica que la dinámica se reduce a una ecuación diferencial de tipo elíptico. Para encontrar la solución analítica, habría que utilizar funciones elípticas.

# Sustituimos las soluciones de w2^2 y w3^2
omega2_sq_expr = sol[omega2_sq]
omega3_sq_expr = sol[omega3_sq]

# Ecuación de Euler para w1
domega1 = ((I2 - I3)/(I1)) * sp.sqrt(omega2_sq_expr * omega3_sq_expr)

# Elevamos al cuadrado para evitar raíces
eq_reduced = sp.simplify(domega1**2)

solucion = sp.Eq(omega1_dot**2, eq_reduced)

# Mostrar ecuaciones
display(solucion)

$\displaystyle \left(\frac{d}{d t} \omega_{1}{\left(t \right)}\right)^{2} = \frac{\left(I_{2} - I_{3}\right)^{2} \omega_{2}^{2}{\left(t \right)} \omega_{3}^{2}{\left(t \right)}}{I_{1}^{2}}$

Esta solución puede ser bastante complicada de interpretar. Por eso vamos a hacer una evaluación numérica.

Resolución numérica

Usamos numpy para evaluar numéricamente la evolución de w1, w2 y w3. Por simplicidad, tomamos $M=1$ kg, $d=1$ m. Sustituímos los valores que tenemos para los momentos principales:

import numpy as np
import matplotlib.pyplot as plt


# Valores arbitrarios (importante: que sean distintos)
I1 = 121/30   # eje intermedio
I2 = 17/6   # eje menor
I3 = 103/15   # eje mayor

Ecuaciones de Euler

Plantemaos las ecuaciones para resolverlas numéricamente. Las resolvemos mediante el método de Euler.

def euler_equations(w, t):
    w1, w2, w3 = w

    dw1 = ((I2 - I3)/I1) * w2 * w3
    dw2 = ((I3 - I1)/I2) * w3 * w1
    dw3 = ((I1 - I2)/I3) * w1 * w2

    return np.array([dw1, dw2, dw3])

# -----------------------------
# Integración (método simple)
# -----------------------------
dt = 0.01
T = 80
N = int(T/dt)

Valores para $\omega$

Probamos con un valor inicial, donde la velocidad en el eje $\vec{e}_1$:

$\omega = (\omega_1, \omega_1/100, \omega_1/100)$

Este eje es el intermedio de los momentos de inercia. Esperamos una oscilación inestable.

# Velocidad angular inicial:
# giro principalmente alrededor del eje intermedio + pequeña perturbación
w = np.array([1.0, 0.01, 0.01])

Resolvemos numéricamente las ecuaciones de Euler:

trajectory1 = []

for _ in range(N):
    trajectory1.append(w.copy())

    # método de Euler explícito
    w = w + dt * euler_equations(w, 0)

trajectory1 = np.array(trajectory1)

# -----------------------------
# 5. Visualización
# -----------------------------
t_vals = np.linspace(0, T, N)

plt.figure()
plt.plot(t_vals, trajectory1[:,0], label='ω1')
plt.plot(t_vals, trajectory1[:,1], label='ω2')
plt.plot(t_vals, trajectory1[:,2], label='ω3')

plt.xlabel("Tiempo")
plt.ylabel("Velocidad angular")
plt.legend()
plt.title("Efecto Dzhanibekov (inestabilidad del eje intermedio)")

plt.show()

Por comparar, probamos con una velocidad inicial dominante en alguno de los otros ejes:

# Velocidad angular inicial:
# giro principalmente alrededor del eje menor + pequeña perturbación
w = np.array([0.01, 1.0, 0.01])

trajectory2 = []

for _ in range(N):
    trajectory2.append(w.copy())

    # método de Euler explícito
    w = w + dt * euler_equations(w, 0)

trajectory2 = np.array(trajectory2)

# -----------------------------
# 5. Visualización
# -----------------------------
t_vals = np.linspace(0, T, N)

plt.figure()
plt.plot(t_vals, trajectory2[:,0], label='ω1')
plt.plot(t_vals, trajectory2[:,1], label='ω2')
plt.plot(t_vals, trajectory2[:,2], label='ω3')

plt.xlabel("Tiempo")
plt.ylabel("Velocidad angular")
plt.legend()
plt.title("Efecto Dzhanibekov (estabilidad del eje menor)")

plt.show()

# Velocidad angular inicial:
# giro principalmente alrededor del eje menor + pequeña perturbación
w = np.array([0.01, 1.0, 0.01])

trajectory3 = []

for _ in range(N):
    trajectory3.append(w.copy())

    # método de Euler explícito
    w = w + dt * euler_equations(w, 0)

trajectory3 = np.array(trajectory3)

# -----------------------------
# 5. Visualización
# -----------------------------
t_vals = np.linspace(0, T, N)

plt.figure()
plt.plot(t_vals, trajectory3[:,0], label='ω1')
plt.plot(t_vals, trajectory3[:,1], label='ω2')
plt.plot(t_vals, trajectory3[:,2], label='ω3')

plt.xlabel("Tiempo")
plt.ylabel("Velocidad angular")
plt.legend()
plt.title("Efecto Dzhanibekov (estabilidad del eje mayor)")

plt.show()

Como conclusión final, y a la vista de las gráficas, podemos decir que la rotación en torno a los ejes $\vec{e}_2$ y $\vec{e}_3$ (los de momentos de inercia mayor y menor) es una rotación estable, pero la rotación en torno al eje intermedio ($\vec{e}_1$) es inestable, lo que se traduce en una transferencia de la rotación hacia los otros ejes.

ANEXO: animación

Por probar, podemos intentar hacer una animación con matplotlib del objeto rotando.

# =========================================
# Animación 3D
# =========================================
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from matplotlib import animation
from IPython.display import HTML

d = 1.0

# Centros respecto al CM
centers = np.array([
    [-d, 3*d/5, 0],
    [0, 3*d/5, 0],
    [d, 3*d/5, 0],
    [0, -2*d/5, 0],
    [0, -7*d/5, 0]
])

def cube_vertices(center, size):
    cx, cy, cz = center
    s = size/2
    return np.array([
        [cx-s, cy-s, cz-s],
        [cx+s, cy-s, cz-s],
        [cx+s, cy+s, cz-s],
        [cx-s, cy+s, cz-s],
        [cx-s, cy-s, cz+s],
        [cx+s, cy-s, cz+s],
        [cx+s, cy+s, cz+s],
        [cx-s, cy+s, cz+s]
    ])

faces_idx = [
    [0,1,2,3],[4,5,6,7],
    [0,1,5,4],[2,3,7,6],
    [1,2,6,5],[4,7,3,0]
]

cubes = [cube_vertices(c, d) for c in centers]

# Momentos de inercia
I1, I2, I3 = 4.0, 2.0, 6.0

def euler(w):
    w1, w2, w3 = w
    return np.array([
        ((I2 - I3)/I1)*w2*w3,
        ((I3 - I1)/I2)*w3*w1,
        ((I1 - I2)/I3)*w1*w2
    ])

# Simulación
dt = 0.1
steps = 400
w = np.array([1.0, 0.01, 0.01])

ws = []
for _ in range(steps):
    ws.append(w.copy())
    w = w + dt * euler(w)

# Rotaciones
def rotation_matrix(axis, theta):
    norm = np.linalg.norm(axis)
    if norm == 0:
        return np.eye(3)
    axis = axis / norm
    K = np.array([
        [0, -axis[2], axis[1]],
        [axis[2], 0, -axis[0]],
        [-axis[1], axis[0], 0]
    ])
    return np.eye(3) + np.sin(theta)*K + (1-np.cos(theta))*(K@K)

R = np.eye(3)
Rs = []

for w in ws:
    theta = np.linalg.norm(w)*dt
    R = rotation_matrix(w, theta) @ R
    Rs.append(R.copy())

# Figura
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

def draw_frame(i):
    ax.cla()
    R = Rs[i]

    for cube in cubes:
        rotated = (R @ cube.T).T
        faces = [[rotated[j] for j in face] for face in faces_idx]
        ax.add_collection3d(Poly3DCollection(faces))

    ax.set_xlim(-2,2)
    ax.set_ylim(-2,2)
    ax.set_zlim(-2,2)

ani = animation.FuncAnimation(fig, draw_frame, frames=steps, interval=50)

HTML(ani.to_jshtml())

WARNING:matplotlib.animation:Animation size has reached 21010437 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.

Once Loop Reflect