Rigid Body Dynamics¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
L'esercizio implementa un semplice modello di $\textbf{corpo rigido}$ soggetto a gravità, utilizzando concetti di meccanica classica, in particolare energia, equazioni di Eulero, formalismo lagrangiano e hamiltoniano.
Energia cinetica¶
L'energia cinetica totale del corpo rigido è composta da una parte traslazionale e una parte rotazionale: $$\large T = \frac{1}{2} m \, \mathbf{v}^2 + \frac{1}{2} \boldsymbol{\omega}^T \mathbf{I} \boldsymbol{\omega}, $$ dove:
- $m$ è la massa del corpo,
- $\mathbf{v}$ è la velocità del centro di massa,
- $\boldsymbol{\omega}$ è il vettore velocità angolare,
- $\mathbf{I}$ è il tensore di inerzia.
Energia potenziale¶
L'energia potenziale gravitazionale è definita come: $$\large V = - \mathbf{g} \cdot \mathbf{r}, $$ dove $\mathbf{g}$ è il vettore accelerazione di gravità e $\mathbf{r}$ la posizione del centro di massa.
Lagrangiana ed equazioni di Eulero--Lagrange¶
La Lagrangiana del sistema è: $$\large \mathcal{L} = T - V. $$
Il codice include una funzione che approssima numericamente le equazioni di Eulero--Lagrange: $$\large \frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial \dot{q}} \right) \frac{\partial \mathcal{L}}{\partial q} = 0,$$ utilizzando differenze finite.
Hamiltoniana
L'Hamiltoniana del sistema è calcolata come: $$\large H = \mathbf{p} \cdot \dot{\mathbf{q}} - \mathcal{L}, $$
e, nel caso del corpo rigido considerato, assume la forma: $$\large H = \frac{\mathbf{p}^2}{2m}\frac{1}{2} \boldsymbol{\omega}^T \mathbf{I} \boldsymbol{\omega}m , \mathbf{g} \cdot \mathbf{r}. $$
Dinamica rotazionale
La dinamica rotazionale del corpo rigido senza momenti esterni è descritta dalle \textbf{equazioni di Eulero}: $$\large \dot{\boldsymbol{\omega}} = \mathbf{I}^{-1} \left( \mathbf{I}\boldsymbol{\omega} \times \boldsymbol{\omega} \right). $$
Dinamica traslazionale¶
L'accelerazione del centro di massa è semplicemente: $$\large \ddot{\mathbf{r}} = \mathbf{g}, $$ coerente con il moto di un corpo soggetto esclusivamente alla forza di gravità.
Conclusione¶
Il codice calcola:
- energia cinetica e potenziale,
- Hamiltoniana del sistema,
- accelerazione del centro di massa,
- evoluzione della velocità angolare.
Si tratta quindi di un esempio didattico di meccanica analitica applicata a un corpo rigido tridimensionale.
import numpy as np
def kinetic_energy(mass, velocity, omega, inertia_tensor):
"""
Calculate the kinetic energy of a rigid body.
:param mass: Mass of the rigid body.
:param velocity: Velocity of the center of mass.
:param omega: Angular velocity vector.
:param inertia_tensor: Inertia tensor of the rigid body.
:return:
Kinetic energy.
"""
translational_ke = 0.5 * mass * np.dot(velocity, velocity)
rotational_ke = 0.5 * np.dot(omega, np.dot(inertia_tensor, omega))
return translational_ke + rotational_ke
def potential_energy(gravity, position):
"""
Calculate the potential energy of a rigid body.
:param gravity: Gravitational acceleration vector.
:param position: Position of the center of mass.
:return:
Potential energy.
"""
return -np.dot(gravity, position)
def euler_lagrange(L, q, q_dot, time):
"""
Evaluate the Euler-Lagrange equation for a given Lagrangian.
:param L: Lagrangian of the system.
:param q: Generalized coordinates.
:param q_dot: Derivatives of generalized coordinates.
:param time: Time vector.
:return:
Euler-Lagrange terms.
"""
dL_dq_dot = np.gradient(L, q_dot)
dL_dq = np.gradient(L, q)
d_dt_dL_dq_dot = np.gradient(dL_dq_dot, time)
return d_dt_dL_dq_dot - dL_dq
def Hamiltonian(L, q_dot, p):
"""
Calculate the Hamiltonian from the Lagrangian and conjugate
momenta.
:param L: Lagrangian of the system.
:param q_dot: Derivatives of generalized coordinates.
:param p: Conjugate momenta.
:return:
Hamiltonian of the system.
"""
return np.dot(q_dot, p) - L
def euler_rotation(omega, inertia_tensor):
"""
Equazioni di Eulero per corpo rigido senza momenti esterni
Solve Euler's equations for rotation.
:param omega: Angular velocity vector.
:param inertia_tensor: Inertia tensor of the rigid body.
:return:
Time derivatives of angular velocities.
"""
I = inertia_tensor
return np.linalg.inv(I) @ np.cross(I @ omega, omega)
def cm_acceleration(gravity):
"""Accelerazione del centro di massa"""
return gravity
def hamiltonian(mass, momentum, omega, inertia_tensor, gravity, position):
translational = np.dot(momentum, momentum) / (2 * mass)
rotational = 0.5 * np.dot(omega, inertia_tensor @ omega)
potential = -mass * np.dot(gravity, position)
return translational + rotational + potential
mass = 5.0
velocity = np.array([2.0, 0.0, 0.0])
momentum = mass * velocity
omega = np.array([0.0, 0.0, 2.0])
inertia_tensor = np.diag([1.0, 1.5, 2.0])
gravity = np.array([0.0, -9.81, 0.0])
position = np.array([0.0, 1.0, 0.0])
# Energie
T = kinetic_energy(mass, velocity, omega, inertia_tensor)
V = potential_energy(gravity, position)
# Hamiltoniana
H = hamiltonian(mass, momentum, omega, inertia_tensor, gravity, position)
# Dinamica
a_cm = cm_acceleration(gravity)
omega_dot = euler_rotation(omega, inertia_tensor)
print("Energia cinetica:", T)
print("Energia potenziale:", V)
print("Hamiltoniana:", H)
print("Accelerazione CM:", a_cm)
print("Accelerazione angolare:", omega_dot)
Energia cinetica: 14.0 Energia potenziale: 9.81 Hamiltoniana: 63.050000000000004 Accelerazione CM: [ 0. -9.81 0. ] Accelerazione angolare: [0. 0. 0.]
Cubo rotante rispetto ai tre assi¶
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from IPython.display import display, HTML
# ---------- Cube geometry ----------
L = 1.0
vertices = np.array([
[-L/2, -L/2, -L/2],
[ L/2, -L/2, -L/2],
[ L/2, L/2, -L/2],
[-L/2, L/2, -L/2],
[-L/2, -L/2, L/2],
[ L/2, -L/2, L/2],
[ L/2, L/2, L/2],
[-L/2, L/2, L/2]
])
faces = [
[0,1,2,3],
[4,5,6,7],
[0,1,5,4],
[2,3,7,6],
[1,2,6,5],
[0,3,7,4]
]
# ---------- Rotation parameters ----------
omega = np.array([1.0, 0.7, 1.5]) # [rad/s] rotation around x, y, z
dt = 0.05
theta = np.zeros(3)
def rotation_matrix_xyz(theta):
"""Return combined rotation matrix from Euler angles (x, y, z)"""
cx, cy, cz = np.cos(theta)
sx, sy, sz = np.sin(theta)
Rx = np.array([[1, 0, 0],
[0, cx, -sx],
[0, sx, cx]])
Ry = np.array([[cy, 0, sy],
[0, 1, 0],
[-sy, 0, cy]])
Rz = np.array([[cz, -sz, 0],
[sz, cz, 0],
[0, 0, 1]])
return Rz @ Ry @ Rx # ZYX convention
# ---------- Plot ----------
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection="3d")
ax.set_xlim(-1.5,1.5)
ax.set_ylim(-1.5,1.5)
ax.set_zlim(-1.5,1.5)
ax.set_box_aspect([1,1,1])
poly = Poly3DCollection(
[[vertices[i] for i in face] for face in faces],
facecolor="cyan",
edgecolor="k",
alpha=0.7
)
ax.add_collection3d(poly)
# ---------- Animation ----------
def update(frame):
global theta
theta += omega * dt
R = rotation_matrix_xyz(theta)
rotated = (R @ vertices.T).T
poly.set_verts([[rotated[i] for i in face] for face in faces])
return poly,
anim = FuncAnimation(fig, update, frames=200, blit=True)
# ---------- Save GIF ----------
anim.save("cubo_rotante.gif", writer="pillow", fps=20)
plt.close()
# ---------- Display in Jupyter ----------
display(HTML(f'<img src="cubo_rotante.gif" style="width:400px;height:400px;">'))
