Cuando se usan directamente las fórmulas de colisión elástica, aparece un problema común: una caja que cae al suelo no se detiene, sino que continúa hundiéndose lentamente. La fórmula elástica funciona bien para rebotes, pero falla al manejar el contacto sostenido entre cuerpos.
Enfoque basado en restricciones
La solución más extendida en motores físicos actuales es organizar las ecuaciones mediante restricciones en lugar de aplicar fórmulas físicas directamente. Una restricción es una regla que obliga a dos cuerpos a cumplir una condición específica (ecuación de restricción).
Restricción durante el contacto
Consideremos dos cuerpos: A (inferior) y B (superior). Para evitar que B se hunda, debemos mantener constante la distancia de penetración actual.
- Función de penetración: \( C = (\mathbf{a} - \mathbf{b}) \cdot \mathbf{n} \), donde \(\mathbf{n}\) es la normal de colisión.
- Para que \(C\) no cambie, su derivada debe ser cero: \( \dot{C} = 0 \).
La derivada de \(C\) respecto al tiempo es:
\[ \dot{C} = (\mathbf{v}_{pa} - \mathbf{v}_{pb}) \cdot \mathbf{n} \] donde \(\mathbf{v}_{pa}\) y \(\mathbf{v}_{pb}\) son las velocidades en los puntos de contacto (combinación de velocidad lineal y angular):
\[ \mathbf{v}_{pa} = \mathbf{v}_a + \boldsymbol{\omega}_a \times \mathbf{r}_a,\quad \mathbf{v}_{pb} = \mathbf{v}_b + \boldsymbol{\omega}_b \times \mathbf{r}_b \] con \(\mathbf{r}_a = \mathbf{a} - \mathbf{O}_a\), \(\mathbf{r}_b = \mathbf{b} - \mathbf{O}_b\).
Para lograr \(\dot{C}=0\), debemos ajustar las velocidades para que la velocidad relativa proyectada en la normal sea cero. El método común es usar impulso para modificar la velocidad.
Cálculo del impulso
El impulso cambia el momento lineal y angular:
- Velocidad lineal tras impulso: \(\mathbf{v}_a' = \mathbf{v}_a + \frac{\mathbf{j}}{m_a}\)
- Velocidad angular tras impulso: \(\omega_a' = \omega_a + \frac{(\mathbf{r}_a \times \mathbf{j})}{I_a}\)
La velocidad relativa en el punto de contacto tras el impulso debe ser cero en la dirección normal:
\[ (\mathbf{v}_{pb}' - \mathbf{v}_{pa}') \cdot \mathbf{n} = 0 \] Resolviendo se obtiene la magnitud del impulso normal:
\[ j_n = -\frac{(\mathbf{v}_{pb} - \mathbf{v}_{pa}) \cdot \mathbf{n}}{ \frac{1}{m_a} + \frac{1}{m_b} + \frac{(\mathbf{r}_a \times \mathbf{n})^2}{I_a} + \frac{(\mathbf{r}_b \times \mathbf{n})^2}{I_b} } \] Este resultado es idéntico al de una colisión elástica con coeficiente de restitución \(e=0\). Surge la pregunta: ¿por qué sigue siendo inestable?
Causas de la inestabilidad
- Errores de punto flotante: pequeños errores en el impulso se acumulan con el tiempo, provocando un hundimiento gradual.
- Discretización de puntos de contacto: se usan dos puntos de contacto independientes, pero al aplicar correcciones a uno, el otro puede desajustarse.
La solución es iterar la corrección varias veces dentro de un mismo paso de simulación. Con 3 iteraciones el hundimiento se reduce; con 10 es casi imperceptible; con 50 parece perfecto, pero en juegos reales 50 iteraciones por frame son inviables.
Implementación en código
Se define una clase CuerpoRigido que almacena propiedades lineales y angulares, y aplica impulsos. La clase Fisica2D gestiona la detección de colisiones y la resolución iterativa.
Clase CuerpoRigido (fragmento):
public class CuerpoRigido : MonoBehaviour
{
private int id;
public Vector2 tamano; // forma fija: caja
public float masa;
private float invMasa;
public Vector2 velocidad;
public Vector2 posicion;
public float inercia;
private float invInercia;
public float velAngular;
public float rotacion;
public void AplicarImpulso(Vector2 impulso)
{
velocidad += impulso * invMasa;
}
public void AplicarImpulsoTorque(Vector2 r, Vector2 impulso)
{
float torque = r.x * impulso.y - r.y * impulso.x;
velAngular += torque * invInercia;
}
public Vector2 VelocidadEnPunto(Vector2 r)
{
Vector2 vAng = new Vector2(-r.y, r.x) * velAngular;
return velocidad + vAng;
}
public void PreUpdate(Vector2 gravedad, float dt)
{
velocidad += (fuerza * invMasa + gravedad) * dt;
velAngular += torque * invInercia * dt;
// aplicar impulsos instantáneos
velocidad += impulsoFuerza * invMasa;
velAngular += impulsoTorque * invInercia;
impulsoFuerza = Vector2.zero;
impulsoTorque = 0;
}
public void PostUpdate(float dt)
{
posicion += velocidad * dt;
rotacion += velAngular * dt;
}
}
Clase Fisica2D (fragmento):
public class Fisica2D : MonoBehaviour
{
public int maxIteraciones = 10;
private List<CuerpoRigido> cuerpos = new List<CuerpoRigido>();
private Dictionary<ParColisionKey, ParColision> paresColision = new Dictionary<ParColisionKey, ParColision>();
void Step(float dt)
{
// PreUpdate: fuerzas continuas e impulsos
foreach (var c in cuerpos) c.PreUpdate(gravedad, dt);
// Detectar colisiones y actualizar pares
DetectarColisiones();
// Resolver restricciones iterativamente
for (int i = 0; i < maxIteraciones; i++)
{
foreach (var par in paresColision.Values)
{
ResolverContacto(par, dt);
}
}
// PostUpdate: integrar posiciones
foreach (var c in cuerpos) c.PostUpdate(dt);
}
void ResolverContacto(ParColision par, float dt)
{
var A = par.cuerpoA;
var B = par.cuerpoB;
for (int i = 0; i < par.numContactos; i++)
{
var contacto = par.contactos[i];
Vector2 ra = contacto.punto - A.posicion;
Vector2 rb = contacto.punto - B.posicion;
Vector2 vRel = B.VelocidadEnPunto(rb) - A.VelocidadEnPunto(ra);
Vector2 n = contacto.normal;
float vn = Vector2.Dot(vRel, n);
float invMasaTotal = A.invMasa + B.invMasa;
float raN = Vector2.Dot(ra, n);
float rbN = Vector2.Dot(rb, n);
invMasaTotal += A.invInercia * (Vector2.Dot(ra,ra) - raN*raN)
+ B.invInercia * (Vector2.Dot(rb,rb) - rbN*rbN);
float impulso = -vn / invMasaTotal;
impulso = Mathf.Max(impulso, 0); // solo fuerzas normales positivas
Vector2 j = impulso * n;
A.AplicarImpulso(-j);
A.AplicarImpulsoTorque(ra, -j);
B.AplicarImpulso(j);
B.AplicarImpulsoTorque(rb, j);
}
}
}
Observación sobre productos vectoriales
En 2D, la velocidad angular es un escalar (componente z). Para el producto vectoiral \(\boldsymbol{\omega} \times \mathbf{r}\) en 2D se usa:
public static Vector2 ProductoZ(float z, Vector2 v)
{
return new Vector2(-z * v.y, z * v.x);
}