FEM para la dinámica de cuerpos rígidos: Manejo de rotaciones basado en cuaterniones y tensores de inercia

Ecuaciones Básicas de la Dinámica de Cuerpos Rígidos#

En un cálculo de FSI o en un montón de DEM, un cuerpo rígido solo tiene que mantener fijas las distancias entre puntos. No hay deformación. Lo que suele romperse primero es la representación de la rotación.

La velocidad de la partícula $i$ es traslación del centro de masa más una rotación.

$\mathbf{v}_i = \mathbf{v}_{cm} + \boldsymbol{\omega} \times \mathbf{r}_i$

$\mathbf{v}_{cm}$ es la velocidad del centro de masa, $\boldsymbol{\omega}$ la velocidad angular y $\mathbf{r}_i$ el desplazamiento desde el centro de masa.

El control de eje y el de ω cambian el eje de giro y el tiempo de una vuelta, aunque la magnitud de la velocidad angular se quede igual.

Al subir ω el cuerpo cierra la vuelta antes. Con el eje sobre una diagonal del cuerpo, se observa si algún ángulo queda trabado en un solo eje de coordenadas.

Momento Lineal y Momento Angular#

El momento lineal total a partir de la segunda ley de Newton es:

$\mathbf{P} = M \mathbf{v}_{cm}, \qquad \frac{d\mathbf{P}}{dt} = \mathbf{F}_{ext}$

Donde $M = \sum_i m_i$ es la masa total. La aceleración del centro de masa es:

$M \dot{\mathbf{v}}_{cm} = \mathbf{F}_{ext}$

El momento angular total es:

$\mathbf{L} = \mathbf{I} \boldsymbol{\omega}, \qquad \frac{d\mathbf{L}}{dt} = \boldsymbol{\tau}$

Donde $\mathbf{I}$ es el tensor de inercia (inertia tensor) y $\boldsymbol{\tau}$ es el par motor o torque total.

Tensor de Inercia#

El tensor de inercia con respecto al centro de masa es:

$I_{jk} = \sum_i m_i \left( |\mathbf{r}_i|^2 \delta_{jk} - r_{i,j} r_{i,k} \right)$

Extendiéndolo a un continuo:

$\mathbf{I} = \int_V \rho \left( |\mathbf{r}|^2 \mathbf{E} - \mathbf{r} \otimes \mathbf{r} \right) dV$

Donde $\mathbf{E}$ es el tensor identidad y $\mathbf{r} \otimes \mathbf{r}$ es el producto externo (outer product).

Cuando un objeto gira, el tensor de inercia se transforma mediante la matriz de orientación $\mathbf{R}$.

$\mathbf{I}(t) = \mathbf{R}(t) \, \mathbf{I}_0 \, \mathbf{R}^T(t)$

Donde $\mathbf{I}_0$ es el tensor de inercia inicial en el sistema de referencia del cuerpo (body frame), el cual solo necesita calcularse una vez al inicio de la simulación.

Representación de la Rotación mediante Cuaterniones#

Representar la rotación con ángulos de Euler provoca el problema del bloqueo del cardán (gimbal lock). Para evitar esto, se utilizan cuaterniones (quaternions).

Un cuaternión $\mathbf{q} = (q_0, q_1, q_2, q_3)$ consta de una parte escalar y una parte vectorial.

$\mathbf{q} = q_0 + q_1 \mathbf{i} + q_2 \mathbf{j} + q_3 \mathbf{k}$

Una rotación de un ángulo $\theta$ sobre un eje unitario $\hat{\mathbf{n}}$ es:

$\mathbf{q} = \cos\frac{\theta}{2} + \sin\frac{\theta}{2}\left(n_x \mathbf{i} + n_y \mathbf{j} + n_z \mathbf{k}\right)$

La rotación compuesta de dos cuaterniones $\mathbf{q}_1$ y $\mathbf{q}_2$ se expresa como un producto.

$\mathbf{q}_{12} = \mathbf{q}_2 \mathbf{q}_1$

Utilizando esto, la matriz de orientación $\mathbf{R}$ se expresa explícitamente en términos de los componentes del cuaternión.

La derivada temporal de un cuaternión está vinculada a la velocidad angular $\boldsymbol{\omega}$.

$\dot{\mathbf{q}} = \frac{1}{2} \begin{pmatrix} 0 \\ \boldsymbol{\omega} \end{pmatrix} \otimes \mathbf{q}$

La diagonal de $I_{jk}$ es $m(|\mathbf{r}|^2-r_j^2)$. Fuera de la diagonal queda $-m r_j r_k$. Si toda la masa está sobre el eje $x$, entonces $r_y=r_z=0$, $I_{xx}=0$ y $I_{yy}=I_{zz}=m x^2$. Alrededor de ese eje no hay brazo de palanca. $\mathbf{I}_0$ es esa suma, guardada una vez en el marco del cuerpo. En el marco del mundo solo hay que rehacer $\mathbf{R}\mathbf{I}_0\mathbf{R}^T$. Multiplicar un tensor del cuerpo por una velocidad angular del mundo inclina el eje del torque. La misma congruencia aparece al cambiar de coordenadas un tensor de rigidez, que es lo que sigue la transformación del tensor constitutivo dentro del elemento.

En $\theta=0$, $\cos(\theta/2)=1$ y $\sin(\theta/2)=0$, así que $\mathbf{q}=(1,0,0,0)$. Es el cuaternión con el que arranca el código. En $\theta=\pi$, $\cos(\pi/2)=0$ y $\sin(\pi/2)=1$: el escalar se apaga y la parte vectorial es solo el eje $\hat{\mathbf{n}}$. En $\theta=2\pi$, $\cos\pi=-1$ y $\mathbf{q}=(-1,0,0,0)$. La rotación coincide con la identidad. Solo cambió el signo. Un cuaternión cubre una vuelta completa dos veces. La entrada $(1,1)$ de $\mathbf{R}$ es $1-2(q_2^2+q_3^2)$, que vale 1 para esa identidad, y las entradas fuera de la diagonal valen 0. El $1/2$ delante de $\dot{\mathbf{q}}$ es el mismo factor que $\theta/2$ en la fórmula del ángulo. Integrarlo tal cual aleja $|\mathbf{q}|$ de 1. Un tetraedro lineal tiene 4 nodos y 3 desplazamientos por nodo, 12 grados de libertad. La restricción rígida los deja en 3 traslaciones más 3 rotaciones, 6 en total. $\sum_i L_i=1$ solo mantiene cerradas las coordenadas. No fabrica esos 6 por su cuenta.

Discretización FEM y Funciones de Forma Tetraédricas#

Al discretizar el campo de desplazamiento de un cuerpo rígido en FEM, las funciones de forma (shape functions) de un elemento tetraédrico son:

$N_i(\xi, \eta, \zeta) = L_i, \quad i = 1, 2, 3, 4$

Donde $L_i$ son coordenadas baricéntricas adimensionales y satisfacen la condición $\sum_i L_i = 1$.

En elementos tetraédricos lineales, la interpolación del desplazamiento es:

$\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{4} N_i(\mathbf{x}) \, \mathbf{u}_i$

Al aplicar las restricciones de cuerpo rígido, los grados de libertad se reducen a 6 (3 lineales + 3 rotacionales).

Ejemplo de Implementación en Python#

import numpy as np
 
class RigidBody:
    def __init__(self, mass, I_body):
        self.M = mass
        self.I0 = I_body          # Tensor de inercia en el body frame (3x3)
        self.x_cm = np.zeros(3)   # Posición del centro de masa
        self.v_cm = np.zeros(3)   # Velocidad del centro de masa
        self.q = np.array([1., 0., 0., 0.])  # Cuaternión (w, x, y, z)
        self.omega = np.zeros(3)  # Velocidad angular
 
    def quaternion_to_rotation(self):
        w, x, y, z = self.q
        return np.array([
            [1-2*(y**2+z**2),  2*(x*y - w*z),  2*(x*z + w*y)],
            [2*(x*y + w*z),  1-2*(x**2+z**2),  2*(y*z - w*x)],
            [2*(x*z - w*y),    2*(y*z + w*x),  1-2*(x**2+y**2)]
        ])
 
    def inertia_world(self):
        R = self.quaternion_to_rotation()
        return R @ self.I0 @ R.T
 
    def step(self, F_ext, tau_ext, dt):
        # Integración del movimiento lineal
        a_cm = F_ext / self.M
        self.v_cm += a_cm * dt
        self.x_cm += self.v_cm * dt
 
        # Integración del movimiento rotacional
        I_world = self.inertia_world()
        alpha = np.linalg.solve(I_world, tau_ext - np.cross(self.omega, I_world @ self.omega))
        self.omega += alpha * dt
 
        # Actualización del cuaternión
        omega_q = np.concatenate([[0.], self.omega])
        q_dot = 0.5 * quaternion_multiply(omega_q, self.q)
        self.q += q_dot * dt
        self.q /= np.linalg.norm(self.q)  # Normalización del cuaternión unitario
 
def quaternion_multiply(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2
    ])

Sin la normalización, R deja de ser una rotación#

Si $|\mathbf{q}|$ no vuelve a 1 en cada paso, las columnas de $\mathbf{R}$ pierden longitud unidad. La línea q /= norm(q) es esa división. Un torque calculado en el marco del mundo y multiplicado por $\mathbf{I}_0$ tal cual manda $\boldsymbol{\alpha}$ por otro eje. El tensor del mundo es $\mathbf{I}_{world} = \mathbf{R} \mathbf{I}_0 \mathbf{R}^T$ antes de que $\mathbf{L}=\mathbf{I}\boldsymbol{\omega}$ cuadre. El techo explícito es $\Delta t < 2 / \omega_{max}$. Si $\omega_{max}$ crece, el paso se acorta primero.

Una mezcla recta entre dos orientaciones se hunde dentro de la esfera unidad. SLERP es la curva que vale $\mathbf{q}_0$ en $t=0$ y $\mathbf{q}_1$ en $t=1$.

$\mathbf{q}(t) = \mathbf{q}_0 \left(\mathbf{q}_0^{-1} \mathbf{q}_1\right)^t$

La interpolación lineal no se queda en esa esfera. Contar la libertad de rotación junto con la deformación del elemento hace que la restricción se coma el movimiento, del mismo modo que el bloqueo por cortante en una cáscara delgada endurece una placa que debería flexionarse. Un cuerpo rígido tiene 6 grados de libertad, no 12.