34. Simular partículas transportadas por un flujo

34.1 Dos descripciones que se complementan

La descripción euleriana asigna velocidad, presión y otras propiedades a cada punto del espacio. La descripción lagrangiana sigue la posición y velocidad de entidades identificables.

Para simular partículas transportadas por un flujo combinamos ambas: evaluamos el campo euleriano en la posición instantánea de cada partícula y usamos ese dato para actualizar su estado lagrangiano.

34.2 Estado mínimo de una partícula

En dos dimensiones, el estado dinámico básico es:

Diámetro y densidad suelen ser parámetros. Posición y velocidad cambian. Mantener esta separación facilita reiniciar, guardar y comparar simulaciones.

34.3 Partícula trazadora ideal

Una trazadora sin inercia adopta instantáneamente la velocidad local del fluido:

No necesita una ecuación independiente para velocidad. Es útil para visualizar trayectorias y transporte, pero no representa sedimentación ni retraso inercial.

34.4 Partícula inercial con arrastre de Stokes

En un modelo diluido, para esfera pequeña y Re_p bajo:

El primer término relaja la velocidad de la partícula hacia la del fluido. El segundo combina gravedad y empuje bajo las hipótesis del modelo.

34.5 Velocidad de deslizamiento terminal

En un flujo uniforme y verticalmente constante, el modelo de Stokes predice:

Su magnitud coincide con |ρ_p − ρ_f|gd²/(18μ). Esta comprobación conecta la ecuación dinámica con el tema anterior.

34.6 Campo de velocidad como función

function campoVelocidad(x, y, tiempo) {
  const velocidadBase = 0.25;
  const omega = 1.2;
  const dx = x - 0.5;
  const dy = y - 0.5;
  const envolvente = Math.exp(-(dx * dx + dy * dy) / 0.12);
  return {
    x: velocidadBase - omega * dy * envolvente,
    y: omega * dx * envolvente
  };
}

console.log(campoVelocidad(0.4, 0.6, 0));

La función devuelve m/s para posiciones en metros. El parámetro tiempo queda disponible aunque este ejemplo sea estacionario.

34.7 Trayectoria y línea de corriente

En un campo estacionario, la trayectoria de una trazadora coincide con una línea de corriente. En un campo variable con el tiempo pueden ser diferentes.

Una captura instantánea de flechas no determina por sí sola por dónde pasó una partícula. La trayectoria requiere integrar en el tiempo.

34.8 Interpolación de campos discretos

Si el campo proviene de una malla, rara vez existe un nodo exactamente en la posición de la partícula. Debe interpolarse.

En dos dimensiones, la interpolación bilineal combina los cuatro nodos de la celda. También debe definirse qué ocurre fuera de la malla y cerca de obstáculos.

34.9 Euler explícito para trazadoras

Con paso Δt:

function pasoTrazadora(particula, campo, tiempo, dt) {
  const u = campo(particula.x, particula.y, tiempo);
  return {
    ...particula,
    x: particula.x + u.x * dt,
    y: particula.y + u.y * dt,
    vx: u.x,
    vy: u.y
  };
}

const campoUniforme = () => ({ x: 0.3, y: -0.1 });
const trazadoraInicial = { x: 1, y: 0.5, vx: 0, vy: 0 };
const trazadoraSiguiente = pasoTrazadora(
  trazadoraInicial,
  campoUniforme,
  0,
  0.2
);

console.log({ trazadoraInicial, trazadoraSiguiente });

34.10 Error de Euler

Euler tiene error global de orden Δt. Reducir el paso aproximadamente a la mitad debería reducir el error aproximadamente a la mitad en el régimen asintótico.

Un dibujo suave no demuestra precisión. La prueba correcta compara pasos sucesivamente menores o una solución conocida.

34.11 Punto medio o RK2

Para una trazadora, el método del punto medio evalúa primero una posición intermedia:

function pasoPuntoMedio(posicion, campo, tiempo, dt) {
  const k1 = campo(posicion.x, posicion.y, tiempo);
  const medio = {
    x: posicion.x + 0.5 * dt * k1.x,
    y: posicion.y + 0.5 * dt * k1.y
  };
  const k2 = campo(medio.x, medio.y, tiempo + 0.5 * dt);
  return {
    x: posicion.x + dt * k2.x,
    y: posicion.y + dt * k2.y
  };
}

const campoUniforme = () => ({ x: 0.2, y: -0.1 });
console.log(pasoPuntoMedio({ x: 0, y: 0 }, campoUniforme, 0, 0.5));

RK2 cuesta dos evaluaciones del campo por paso, pero suele mejorar mucho la precisión de trayectorias curvas.

34.12 El problema de una relajación rígida

Euler explícito aplicado a dv/dt = (u − v)/τ puede volverse inestable si Δt es grande respecto de τ. Partículas muy pequeñas poseen τ diminuto.

Reducir Δt hasta resolver el menor τ puede ser costoso. Para arrastre lineal con u aproximadamente constante durante el paso existe una actualización exponencial exacta.

34.13 Actualización exponencial del arrastre

Definimos la velocidad de equilibrio local:

Si u se considera constante en Δt:

Esta actualización es estable para cualquier Δt en la relajación lineal, aunque un paso grande todavía puede resolver mal las variaciones espaciales del campo.

34.14 Implementación de partícula inercial

function pasoInercial(particula, fluido, campo, tiempo, dt, g = 9.81) {
  const u = campo(particula.x, particula.y, tiempo);
  const tau = particula.densidad * particula.diametro ** 2 /
    (18 * fluido.viscosidadDinamica);
  const gravedadEfectiva = g * (1 - fluido.densidad / particula.densidad);
  const equilibrio = { x: u.x, y: u.y + tau * gravedadEfectiva };
  const decaimiento = Math.exp(-dt / tau);
  const vx = equilibrio.x + (particula.vx - equilibrio.x) * decaimiento;
  const vy = equilibrio.y + (particula.vy - equilibrio.y) * decaimiento;
  return {
    ...particula,
    x: particula.x + vx * dt,
    y: particula.y + vy * dt,
    vx,
    vy
  };
}

const agua = { densidad: 1000, viscosidadDinamica: 0.001 };
const esferaInicial = {
  x: 0,
  y: 0,
  vx: 0,
  vy: 0,
  densidad: 2500,
  diametro: 100e-6
};
const flujoUniforme = () => ({ x: 0.2, y: 0 });
const esferaSiguiente = pasoInercial(
  esferaInicial,
  agua,
  flujoUniforme,
  0,
  0.01
);

console.log({ esferaInicial, esferaSiguiente });

La posición usa la velocidad nueva, una actualización semiimplícita sencilla.

34.15 Número de Stokes

Compara el tiempo de respuesta de la partícula con una escala temporal del flujo T_f:

St ≪ 1 indica seguimiento cercano; St ≳ 1 indica inercia apreciable. En un flujo de escala L y velocidad U puede usarse T_f = L/U, de modo que St = τ_pU/L.

34.16 Reynolds de partícula durante la simulación

En cada instante:

Si deja de ser pequeño, el arrastre lineal de Stokes pierde validez. La simulación debe advertirlo o cambiar a una correlación no lineal.

34.17 Elegir el paso temporal

Δt debe resolver las escalas relevantes:

• tiempo de variación del campo;
• tiempo para cruzar una celda o detalle geométrico;
• tiempo de respuesta de partículas, según el integrador;
• colisiones y contacto con fronteras;
• precisión requerida de las trayectorias.

Una actualización estable puede seguir siendo inexacta.

34.18 Condición tipo CFL para advección

Si Δx es la escala espacial que debe resolverse:

Para seguimiento explícito suele buscarse que la partícula no salte muchas celdas por paso. El límite preciso depende del método de interpolación e integración.

34.19 Paso físico y frecuencia de pantalla

requestAnimationFrame se adapta a la pantalla y no ofrece un Δt físico constante. Usar directamente el intervalo entre cuadros vuelve el resultado dependiente del rendimiento.

Una estrategia es acumular tiempo real y ejecutar pasos físicos fijos hasta ponerse al día, con un límite para evitar una espiral de retraso.

34.20 Acumulador de tiempo fijo

function consumirTiempo(estadoReloj, deltaReal, dt, avanzar) {
  estadoReloj.acumulado += Math.min(deltaReal, 0.05);
  let pasos = 0;
  while (estadoReloj.acumulado >= dt && pasos < 100) {
    avanzar(dt);
    estadoReloj.acumulado -= dt;
    pasos++;
  }
  return pasos;
}

const reloj = { acumulado: 0 };
const pasos = consumirTiempo(reloj, 0.0167, 0.002, () => {});
console.log(`Pasos ejecutados: ${pasos}, resto: ${reloj.acumulado.toFixed(4)} s`);

34.21 Fronteras periódicas

En un dominio periódico, una partícula que sale por un borde reaparece por el opuesto:

function envolver(valor, minimo, maximo) {
  const longitud = maximo - minimo;
  return minimo + ((valor - minimo) % longitud + longitud) % longitud;
}

console.log(envolver(1.12, 0, 1));
console.log(envolver(-0.08, 0, 1));

El doble módulo maneja valores negativos de JavaScript.

34.22 Fronteras reflectantes y absorbentes

Una pared reflectante invierte la componente normal de velocidad, quizá multiplicada por un coeficiente de restitución. Una frontera absorbente elimina la partícula.

Detectar el cruce después de un paso grande puede colocar la partícula dentro de la pared. Para mayor precisión se calcula el instante de impacto o se subdivide el paso.

34.23 Inyección y eliminación

Un flujo continuo de partículas requiere una tasa de inyección. Si la tasa es fraccionaria por paso puede acumularse un residuo para no depender de redondeos.

Eliminar elementos mientras se recorre un arreglo puede saltar índices. Alternativas: recorrer hacia atrás, filtrar al final o mantener una lista de partículas activas.

34.24 Acoplamiento una vía y dos vías

En acoplamiento de una vía, el flujo mueve partículas, pero ellas no modifican el flujo. Es razonable para concentraciones pequeñas.

En dos vías, las fuerzas de reacción regresan al fluido. En cuatro vías también importan colisiones partícula–partícula. Aumentar el nivel de acoplamiento cambia las ecuaciones y el costo computacional.

34.25 Colisiones y volumen excluido

Las partículas puntuales pueden atravesarse. Si el problema exige contacto, deben detectarse solapamientos y calcularse impulsos o fuerzas de contacto.

La detección por todos los pares cuesta O(N²). Rejillas espaciales y árboles reducen candidatos cuando N es grande.

34.26 Movimiento browniano

Para partículas microscópicas puede añadirse un desplazamiento aleatorio consistente con difusión. Su amplitud escala con √Δt, no con Δt.

La fuente aleatoria debe ser reproducible mediante semilla si se desean pruebas. Un término de ruido exige interpretar la ecuación estocástica y sus unidades.

34.27 Diagnósticos durante la simulación

• cantidad de partículas activas;
• mínimos y máximos de posición y velocidad;
• Re_p máximo y promedio;
• número de cruces e impactos;
• tiempo simulado y pasos ejecutados;
• presencia de NaN o valores infinitos;
• sensibilidad al reducir Δt.

Una animación plausible puede ocultar una simulación numéricamente incorrecta.

34.28 Detectar estados inválidos

function validarParticula(p, indice) {
  for (const propiedad of ["x", "y", "vx", "vy"]) {
    if (!Number.isFinite(p[propiedad])) {
      throw new Error(`Partícula ${indice}: ${propiedad} no es finita`);
    }
  }
}

validarParticula({ x: 0.2, y: 0.4, vx: 0.1, vy: 0 }, 0);
console.log("Estado válido");

34.29 Reproducibilidad

Guardar semilla, parámetros, versión del modelo, Δt y condiciones iniciales permite repetir un resultado. La tasa de refresco de pantalla no debe cambiar la física.

Una prueba útil ejecuta dos veces sin dibujar y compara estados finales. Renderizar es una observación del modelo, no parte de su evolución.

34.30 Rendimiento

Antes de optimizar, se mide. Los costos frecuentes son evaluación del campo, interpolación, colisiones y dibujo.

Arreglos tipados, recorridos simples, menos asignaciones temporales y dibujo por lotes pueden ayudar. Reducir precisión física para ganar cuadros por segundo no debe hacerse sin cuantificar el impacto.

34.31 Accesibilidad y control de movimiento

Una animación debe ofrecer pausa y avance manual. Si el sistema indica prefers-reduced-motion, es apropiado iniciar pausada.

Los resultados importantes también deben aparecer como texto, porque el canvas por sí solo no comunica el estado a tecnologías de asistencia.

34.32 Laboratorio de transporte de partículas

Compará trazadoras e inerciales en un flujo uniforme con vórtice localizado. El borde horizontal es periódico y las paredes superior e inferior son reflectantes.

34.33 Ejercicio propuesto

Una esfera de d = 100 μm y ρ_p = 2500 kg/m³ se libera en reposo dentro de un flujo uniforme u = 0,20 m/s. El fluido tiene μ = 0,001 Pa·s. Ignorá gravedad y trabajá en una dimensión.

1. Calculá τ_p.
2. Implementá Euler explícito con Δt = 0,0001 s hasta t = 0,01 s.
3. Compará con v(t) = u[1 − exp(−t/τ_p)].
4. Compará la posición con x(t) = u[t − τ_p(1 − exp(−t/τ_p))].
5. Repetí con pasos menores y verificá convergencia.

const d = 100e-6;
const rhoP = 2500;
const mu = 0.001;
const u = 0.2;
const tau = rhoP * d ** 2 / (18 * mu);
const dt = 1e-4;
const tiempoFinal = 0.01;
let v = 0;
let x = 0;

for (let t = 0; t < tiempoFinal - 1e-12; t += dt) {
  v += (u - v) / tau * dt;
  x += v * dt;
}

const vExacta = u * (1 - Math.exp(-tiempoFinal / tau));
const xExacta = u * (tiempoFinal - tau * (1 - Math.exp(-tiempoFinal / tau)));

console.log(`tau = ${(tau * 1000).toFixed(4)} ms`);
console.log(`Euler: v = ${v.toFixed(6)} m/s, x = ${x.toFixed(8)} m`);
console.log(`Exacta: v = ${vExacta.toFixed(6)} m/s, x = ${xExacta.toFixed(8)} m`);
console.log(`Error relativo en x: ${(Math.abs(x - xExacta) / xExacta * 100).toFixed(3)} %`);

34.34 Ideas para recordar

• Una simulación lagrangiana evalúa el campo euleriano en cada partícula.
• Una trazadora sigue u; una partícula inercial responde en un tiempo τ_p.
• St = τ_p/T_f mide la capacidad de seguir cambios del flujo.
• Re_p debe comprobar la validez del arrastre de Stokes.
• La actualización exponencial estabiliza la relajación lineal, pero no elimina el error espacial.
• El paso físico debe ser independiente de la frecuencia de pantalla.
• Las condiciones de borde forman parte del modelo, no solo del dibujo.
• Diagnósticos, semillas y estudios de paso son esenciales para confiar en el resultado.
• La animación debe poder pausarse y comunicar resultados también como texto.

En el próximo tema integraremos los conceptos del curso en un laboratorio completo de fluidos y partículas.