Demo desarrollada en Processing en la que se simula el comportamiento de partículas cuando colisionan entre ellas y contra planos. Se han utilizado dos métodos para calcular estas colisiones.

Modelo de velocidades

Para este modelo se ha descompuesto la velocidad de las partículas para calcular su velocidad tangencial y normal.

Pero primero necesitamos comprobar cuanto hay que desplazar la partícula. Para ello usaremos la siguiente expresión:

$$ \Delta s = \vec{n_{plano}} · distancia $$ Donde n es la normal del Plano.

Luego se calcula la velocidad normal de la partícula:

$$ \vec{V_{n}} = (\vec{n_{plano}} · \vec{V_{particle}}) · \vec{n_{plano}} $$

$$ \vec{V_{t}} = \vec{V_{particle}} - \vec{V_{n}} $$

$$ \vec{V_{salida}} = \vec{V_{t}} - (\vec{V_{particle}} · \vec{n_{plano}}) · \vec{n_{plano}} $$

void computePlanesCollisionsVelocity( Particle par, PlaneSection p){
  float dist = p.getDistance( par.getPos());

  // Step 2 : Restitution
  float dvn = par.getVel().dot( p.getNormal().normalize());
  float as = par.getRadius() - dist;
  PVector delta = PVector.mult( p.getNormal().normalize(), as);
  par.setPos( par.getPos().add( delta));

  // Step 3 : Final Velocity
  PVector vn = PVector.mult( p.getNormal().normalize(), dvn);
  PVector vt = PVector.sub( par.getVel(), vn);
  PVector vs = PVector.sub( vt, PVector.mult( vn, Cr1));
  par.setVel( vs);
}

Para la colisión entre Partículas se trata de la misma forma, pero en lugar de hacer los calculos sobre la normal del plano, se utilizará el vector distancia entre ellas.

void computeParticleCollisionsVelocity( Particle p1, Particle p2, float timeStep){
  PVector p1p2 = PVector.sub( p2.getPos(), p1.getPos());
  float dist = PVector.dist( p1.getPos(), p2.getPos());

  // Descomposition
  PVector normal1 = PVector.mult( p1p2.normalize( null), p1.getVel().dot( p1p2) / dist);
  PVector tan1 = PVector.sub( p1.getVel(), normal1);
  PVector normal2 = PVector.mult( p1p2.normalize( null), p2.getVel().dot( p1p2) / dist);
  PVector tan2 = PVector.sub( p2.getVel(), normal2);

  // Step 2 : Restitution
  float l = p1.getRadius() + p2.getRadius() - dist;
  float vrel = PVector.sub( normal1, normal2).mag();

  if( vrel < 15){
    vrel = 15;
  }

  p1.getPos().add( PVector.mult( normal1, -l/vrel));
  p2.getPos().add( PVector.mult( normal2, -l/vrel));

  // Step 3 : Exit Velocity
  float u1 = normal1.dot( p1p2) / dist;
  float u2 = normal2.dot( p1p2) / dist;
  float v1_ = ((( p1.getMass() - p2.getMass()) * u1) + 2*p2.getMass()*u2) / (p1.getMass() + p2.getMass());
  float v2_ = ((( p2.getMass() - p1.getMass()) * u2) + 2*p1.getMass()*u1) / (p2.getMass() + p1.getMass());

  PVector normal1_ = PVector.mult( p1p2.normalize( null), v1_);
  PVector normal2_ = PVector.mult( p1p2.normalize( null), v2_);

  p1.setVel( PVector.add( normal1_.mult( Cr2), tan1));
  p2.setVel( PVector.add( normal2_.mult( Cr2), tan2));
}

Modelo de Muelles

En este modelo las colisiones se calculan situando un muelle entre los objetos que colisionan. No es necesario hacer restitución de la partícula por que de ello se encargará el muelle.

Se utilizará la fórmula de la fuerza del muelle:

$$ F_{muelle} = -k (L_{actual} - L_{reposo}) $$

Donde según nos interese iremos utilizando los parámetros necesarios.

• L actual será la distancía entre nuestro muelle y el plano o la otra partícula.
• L reposo será el radio de la partícula si colisiona contra el plano, o la suma de radios si colisiona contra otra partícula.

La implementación quedaría de la siguiente forma en colisiones contra Planos:

void computePlanesCollisionsSpring( Particle p, PlaneSection plane){
  float dist = plane.getDistance( p.getPos());

  float lRep = p.getRadius();
  PVector normal = plane.getNormal().normalize( null);

  float Fs = -_ks * ( dist - lRep);

  PVector Fs1 = PVector.mult( normal, Fs);
  PVector v1 = PVector.add( p.getVel(), Fs1.div( p.getMass()).mult( Cr2).mult( _timeStep));

  p.setVel( v1);
}

La implementación quedaría de la siguiente forma en colisiones contra otras partículas:

void computeParticleCollisionsSpring( Particle p1, Particle p2, float timeStep){
  PVector p1p2 = PVector.sub( p2.getPos(), p1.getPos());
  float dist = PVector.dist( p1.getPos(), p2.getPos());

  float lRep = p1.getRadius() + p2.getRadius();
  p1p2.normalize();

  float Fs = _ks * ( dist - lRep);

  PVector Fs1 = p1p2.mult( Fs);
  PVector v1 = PVector.add( p1.getVel(), Fs1.div( p1.getMass()).mult( Cr2).mult( timeStep));
  p1.setVel( v1);
}