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 el análisis de interacción fluido-estructura (FSI) o en las simulaciones del Método de Elementos Discretos (DEM), es esencial describir con precisión el movimiento de los cuerpos rígidos. Un cuerpo rígido es un objeto sin deformación, donde la distancia entre cualquier par de puntos se mantiene siempre constante.

La velocidad de una partícula arbitraria $i$ dentro de un cuerpo rígido se descompone de la siguiente manera:

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

Donde $\mathbf{v}_{cm}$ es la velocidad del centro de masa, $\boldsymbol{\omega}$ es la velocidad angular y $\mathbf{r}_i$ es el vector de desplazamiento desde el centro de masa hasta la partícula $i$.

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}$

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
    ])

Advertencias Prácticas#

1. Olvidar la normalización del cuaternión A medida que se repite la integración numérica, puede ocurrir que $|\mathbf{q}| \neq 1$, haciendo que la matriz de rotación sea no ortogonal. Se debe realizar obligatoriamente q /= norm(q) en cada paso de tiempo.

2. Confusión del sistema de coordenadas del tensor de inercia $\mathbf{I}_0$ se define en el sistema de coordenadas fijo al cuerpo (body frame). Si el torque se calculó en el sistema de coordenadas del mundo (world frame), se debe aplicar tras convertirlo a $\mathbf{I}_{world} = \mathbf{R} \mathbf{I}_0 \mathbf{R}^T$.

3. Selección del paso de tiempo Para garantizar la estabilidad numérica, se debe cumplir la condición $\Delta t < 2 / \omega_{max}$ (similar a la condición CFL) para la frecuencia natural máxima $\omega_{max}$ del cuerpo rígido. Cuanto más rápido gire el objeto, menor será el $\Delta t$ necesario.

4. Interpolación de cuaterniones (SLERP) Al interpolar entre dos orientaciones, se debe utilizar la interpolación lineal esférica (SLERP) en lugar de la interpolación lineal (LERP) para asegurar una rotación de velocidad constante.

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

La combinación de la dinámica de cuerpos rígidos y FEM es la base para las simulaciones de FSI, DEM y dinámica multicuerpo. Una implementación correcta de la representación de rotación basada en cuaterniones permite simular de forma estable cualquier rotación tridimensional sin el bloqueo del cardán.

Cambia eje y ω para confirmar que la rotación cuaterniónica no sufre gimbal lock.