Il corpo rigido a 6 gradi di libertà

Tre coordinate per dire dov'è, tre angoli per dire come è girato: sei numeri, e due leggi della meccanica per farli evolvere.

Sei gradi di libertà

Il simulatore tratta l'aereo come un corpo rigido: un oggetto che non si deforma. Il suo stato in ogni istante è dato da:

posizione del baricentro $\vec x$ e velocità $\vec v$ (assi terra);
orientamento $\mathbf q$ (un quaternione) e velocità angolare $\vec\omega$ (assi corpo).

Le forze cambiano la velocità, i momenti cambiano la velocità angolare. Tutto il resto (angolo d'attacco, velocità indicata, strumenti) si ricava da questi quattro vettori.

Il moto del baricentro: Newton

La seconda legge della dinamica, scritta negli assi terra (che si possono considerare inerziali: la rotazione della Terra è trascurabile per un aereo leggero):

Equazioni della traslazione

$$m\,\frac{d\vec v}{dt}=R(\mathbf q)\left(\vec F_{aero}+\vec T\right)+m\vec g+\vec F_{carrello},\qquad \frac{d\vec x}{dt}=\vec v$$

Le forze aerodinamiche e la spinta sono calcolate negli assi corpo e ruotate negli assi terra con $R(\mathbf q)$; peso e carrello sono già negli assi terra.

La rotazione: le equazioni di Eulero

Per le rotazioni la legge corrispondente è $d\vec H/dt=\vec M$, dove $\vec H=I\vec\omega$ è il momento angolare e $I$ il tensore d'inerzia. Conviene scriverla negli assi corpo, dove $I$ è costante; ma gli assi corpo ruotano, e la derivata di un vettore in un sistema rotante ha un termine in più (lo stesso che produce la forza di Coriolis):

$$\left.\frac{d\vec H}{dt}\right|_{terra}=\left.\frac{d\vec H}{dt}\right|_{corpo}+\vec\omega\times\vec H$$

Equazioni di Eulero

$$I\,\frac{d\vec\omega}{dt}+\vec\omega\times\left(I\,\vec\omega\right)=\vec M$$

Se gli assi corpo sono assi principali d'inerzia, $I=\operatorname{diag}(I_1,I_2,I_3)$ e per componenti:

$$I_1\dot\omega_1=M_1-(I_3-I_2)\,\omega_2\,\omega_3\qquad I_2\dot\omega_2=M_2-(I_1-I_3)\,\omega_3\,\omega_1\qquad I_3\dot\omega_3=M_3-(I_2-I_1)\,\omega_1\,\omega_2$$

Il termine $\vec\omega\times I\vec\omega$ è giroscopico: accoppia le rotazioni attorno ad assi diversi. È lui che fa precedere una trottola, e che rende instabile la rotazione attorno all'asse di inerzia intermedio (laboratorio qui sotto).

I momenti d'inerzia del Cessna

Asse (three.js)	Rotazione	Momento d'inerzia	Raggio giratore $\sqrt{I/m}$
x (ala)	beccheggio
y (verticale)	imbardata
z (fusoliera)	rollio

L'imbardata ha l'inerzia più grande (ruotano sia le ali sia la fusoliera), il rollio la più piccola (la massa è vicina all'asse della fusoliera, a parte le ali). Il raggio giratore è la distanza a cui bisognerebbe concentrare tutta la massa per avere la stessa inerzia. Il modello trascura il prodotto d'inerzia $I_{xz}$, che nell'aereo vero è piccolo ma non nullo.

Laboratorio — l'effetto racchetta da tennis (teorema dell'asse intermedio)

Asse di rotazione iniziale Disturbo

Si integrano le equazioni di Eulero senza momenti esterni, con le inerzie del Cessna, partendo da una rotazione di 1 rad/s attorno a un asse più un piccolo disturbo sugli altri due. Attorno all'asse minimo e massimo la rotazione è stabile; attorno all'intermedio si ribalta periodicamente (effetto Džanibekov: lo vide un cosmonauta sulla stazione Salyut 7 svitando un dado ad alette). Le curve sono le tre componenti di $\vec\omega$; la linea sottile è l'energia, che deve restare costante.

L'orientamento: perché i quaternioni

Gli angoli di Eulero (prua, beccheggio, rollio) hanno due difetti per il calcolo: in verticale si bloccano (gimbal lock, capitolo 1) e le loro derivate richiedono seni, coseni e divisioni per $\cos\theta$. Una matrice di rotazione 3×3 funziona sempre, ma ha 9 numeri per 3 gradi di libertà e, con gli errori di arrotondamento, smette lentamente di essere una rotazione.

Il quaternione unitario è il compromesso ideale: quattro numeri, nessuna singolarità, e si "ripara" con una semplice normalizzazione. Una rotazione di un angolo $\vartheta$ attorno all'asse unitario $\hat n$ è

$$\mathbf q=\left(\cos\frac\vartheta2,\;\hat n\,\sin\frac\vartheta2\right)=(w,\;x,\;y,\;z),\qquad w^2+x^2+y^2+z^2=1.$$

Comporre due rotazioni è il prodotto di quaternioni (non commutativo, come le rotazioni). L'orientamento evolve secondo

Cinematica del quaternione

$$\dot{\mathbf q}=\tfrac12\,\mathbf q\otimes(0,\vec\omega)\qquad\Longrightarrow\qquad \mathbf q_{n+1}\approx\frac{\mathbf q_n\otimes\left(1,\;\tfrac12\vec\omega\,\Delta t\right)}{\left\lVert\mathbf q_n\otimes\left(1,\;\tfrac12\vec\omega\,\Delta t\right)\right\rVert}$$

$\vec\omega$ negli assi corpo, quindi si moltiplica a destra. Per piccoli $\Delta t$, $(1,\tfrac12\vec\omega\Delta t)$ normalizzato è la rotazione di $\lVert\vec\omega\rVert\Delta t$ attorno a $\vec\omega$.

Nel simulatore

Tutto il "motore" del corpo rigido sta in 9 righe (js/fisica.js, righe 203–212):

st.vel.addScaledVector(Fw,dt/m);                 // Newton: v += F/m · dt
st.pos.addScaledVector(st.vel,dt);               //         x += v · dt
Iw.set(Ix*st.w.x,Iy*st.w.y,Iz*st.w.z);           // H = I·ω
tmp.crossVectors(st.w,Iw);                       // ω × (I·ω)
st.w.x+=(Tb.x-tmp.x)/Ix*dt;                      // Eulero: ω += I⁻¹(M − ω×Iω) · dt
st.w.y+=(Tb.y-tmp.y)/Iy*dt;
st.w.z+=(Tb.z-tmp.z)/Iz*dt;
dq.set(st.w.x*dt*0.5,st.w.y*dt*0.5,st.w.z*dt*0.5,1);   // (1, ω·dt/2)
st.q.multiply(dq).normalize();                   // q = q ⊗ dq, poi lunghezza 1

Prima di queste righe il programma ha sommato tutte le forze in Fw (assi terra) e tutti i momenti in Tb (assi corpo). Come si passa da queste equazioni a passi di calcolo discreti è l'argomento del capitolo 8.

Prova tu

Il Cessna rolla a 30°/s e intanto imbarda a 5°/s. Quanto vale il momento giroscopico $\vec\omega\times I\vec\omega$ sull'asse di beccheggio?
Scrivi il quaternione di una rotazione di 90° attorno all'asse verticale y. Verifica che ha lunghezza 1.
Perché nel codice si normalizza il quaternione a ogni passo? Che cosa succederebbe senza?

Soluzioni

$\omega_y=0{,}087$, $\omega_z=0{,}524$ rad/s (segni a parte). Componente x: $\omega_y I_z\omega_z-\omega_z I_y\omega_y=\omega_y\omega_z(I_z-I_y)=0{,}087\cdot0{,}524\cdot(1285-2667)\approx-63$ N·m: piccolo rispetto ai momenti aerodinamici (migliaia di N·m), ma c'è.
$\mathbf q=(\cos45^\circ,\,0,\,\sin45^\circ,\,0)=(0{,}707;\;0;\;0{,}707;\;0)$; $0{,}707^2+0{,}707^2=1$.
Il passo $(1,\tfrac12\vec\omega\Delta t)$ non ha lunghezza esattamente 1: senza normalizzare, il quaternione crescerebbe a ogni passo e la "rotazione" deformerebbe l'aereo (scala e taglio).