Рідина Латтіс-Больцмана у 200 рядках

D2Q9 = 2 виміри, 9 векторів швидкості. Кожна клітинка ґратки зберігає 9 значень розподілу f[0..8] — густину ймовірності частинок, що рухаються в кожному напрямку:

// Напрямки швидкості D2Q9: [ex, ey] для i = 0..8
//   6  2  5
//   3  0  1    (0 = спокій, 1-4 = осьові, 5-8 = діагональні)
//   7  4  8
const EX = [ 0,  1,  0, -1,  0,  1, -1, -1,  1];
const EY = [ 0,  0,  1,  0, -1,  1,  1, -1, -1];
const W  = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; // ваги
const OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]; // індекс протилежного напрямку

const NX = 200, NY = 100;
const N  = NX * NY;

// 9 розподілів на клітинку, збережені послідовно: f[i*N + (y*NX + x)]
const f    = new Float64Array(9 * N);
const fNew = new Float64Array(9 * N);
const wall = new Uint8Array(N);   // 1 = тверда перешкода

function at(x, y) { return y * NX + x; }

Ініціалізуйте кожну клітинку до рівноважного розподілу для рівномірного горизонтального потоку зі швидкістю u0:

const u0 = 0.1; // початкова швидкість потоку (у ґраткових одиницях, тримайте < 0.3)

function feq(i, rho, ux, uy) {
  const eu = EX[i] * ux + EY[i] * uy;
  const u2 = ux*ux + uy*uy;
  return W[i] * rho * (1 + 3*eu + 4.5*eu*eu - 1.5*u2);
}

function initFlow() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      for (let i = 0; i < 9; i++) {
        f[i*N + cell] = feq(i, 1.0, u0, 0.0);
      }
    }
  }
}
initFlow();

Стрімінг: кожен розподіл f[i] переміщується до сусідньої клітинки в напрямку i. Це просто шаблон копіювання масиву:

function stream() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      for (let i = 0; i < 9; i++) {
        // Джерело: ЗВІДКИ прийшов цей розподіл?
        const sx = x - EX[i], sy = y - EY[i];
        if (sx < 0 || sx >= NX || sy < 0 || sy >= NY) {
          fNew[i*N + cell] = f[i*N + cell]; // край: залишаємо
        } else {
          fNew[i*N + cell] = f[i*N + at(sx, sy)];
        }
      }
    }
  }
  f.set(fNew);
}

Оператор зіткнень BGK релаксує f до його рівноваги feq. Час релаксації tau керує вʼязкістю:

// Кінематична вʼязкість: nu = (tau - 0.5) / 3
// tau = 1.0 → nu = 1/6 (помірна вʼязкість)
// tau = 0.6 → nu = 1/30 (менша вʼязкість, Re ≈ u0 * L / nu)
const tau = 0.6;
const omega = 1.0 / tau; // = 1/tau

function collide() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      if (wall[cell]) continue; // пропускаємо тверді клітинки

      // Обчислюємо макроскопічну густину та швидкість
      let rho = 0, ux = 0, uy = 0;
      for (let i = 0; i < 9; i++) {
        const fi = f[i*N + cell];
        rho += fi;
        ux  += EX[i] * fi;
        uy  += EY[i] * fi;
      }
      ux /= rho; uy /= rho;

      // Релаксуємо до рівноваги
      for (let i = 0; i < 9; i++) {
        const eq = feq(i, rho, ux, uy);
        f[i*N + cell] += omega * (eq - f[i*N + cell]);
      }
    }
  }
}

Для твердих стінок: частинки, що стрімляться у клітинку-стінку, відбиваються назад тим самим шляхом, яким прийшли (умова прилипання):

// Додаємо круглу перешкоду
function addCircle(cx, cy, radius) {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if ((x-cx)**2 + (y-cy)**2 < radius**2) wall[at(x,y)] = 1;
    }
  }
}
addCircle(NX * 0.3, NY * 0.5, NY * 0.1);

// Застосовуємо bounce-back ПІСЛЯ стрімінгу
function bounceBack() {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const cell = at(x, y);
      if (!wall[cell]) continue;
      for (let i = 0; i < 9; i++) {
        // Міняємо місцями з протилежним напрямком
        const opp = OPP[i];
        const tmp = f[i*N + cell];
        f[i*N + cell] = f[opp*N + cell];
        f[opp*N + cell] = tmp;
      }
    }
  }
}

Скидайте лівий стовпець до потоку на вході та задавайте нульовий градієнт для правого стовпця на кожному кроці:

function applyBoundaries() {
  // Лівий вхід: фіксована швидкість u0
  for (let y = 0; y < NY; y++) {
    const cell = at(0, y);
    for (let i = 0; i < 9; i++) f[i*N + cell] = feq(i, 1.0, u0, 0.0);
  }
  // Правий вихід: копіюємо з x = NX-2 (нульовий нормальний градієнт)
  for (let y = 0; y < NY; y++) {
    const c = at(NX-1, y), c2 = at(NX-2, y);
    for (let i = 0; i < 9; i++) f[i*N + c] = f[i*N + c2];
  }
}

// Головний цикл симуляції
function step() {
  stream();
  bounceBack();
  collide();
  applyBoundaries();
}

Завихреність (∂uy/∂x − ∂ux/∂y) гарно підкреслює вихори. Використовуйте ImageData для прямого запису пікселів — це набагато швидше за fillRect:

const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
canvas.width = NX; canvas.height = NY;
const img = ctx.createImageData(NX, NY);

function getVelocity(x, y) {
  const cell = at(x, y);
  let ux = 0, uy = 0, rho = 0;
  for (let i = 0; i < 9; i++) {
    const fi = f[i*N + cell];
    rho += fi; ux += EX[i]*fi; uy += EY[i]*fi;
  }
  return { ux: ux/rho, uy: uy/rho };
}

function render() {
  for (let y = 1; y < NY-1; y++) {
    for (let x = 1; x < NX-1; x++) {
      const p = (y * NX + x) * 4;

      if (wall[at(x, y)]) {
        img.data[p] = img.data[p+1] = img.data[p+2] = 60;
        img.data[p+3] = 255;
        continue;
      }

      // Завихреність методом скінченних різниць
      const r = getVelocity(x+1, y), l = getVelocity(x-1, y);
      const t = getVelocity(x, y+1), b = getVelocity(x, y-1);
      const curl = (r.uy - l.uy) - (t.ux - b.ux);

      // Відображаємо завихреність у колір: відʼємна = синій, додатна = червоний
      const v = Math.max(-0.1, Math.min(0.1, curl)) / 0.1; // [-1, 1]
      img.data[p]   = v > 0 ? Math.floor(255 * v) : 0;
      img.data[p+1] = 0;
      img.data[p+2] = v < 0 ? Math.floor(-255 * v) : 0;
      img.data[p+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}

// Анімація: виконуємо кілька фізичних кроків за кадр для швидкості
let running = true;
(function loop() {
  if (!running) return;
  for (let s = 0; s < 10; s++) step(); // 10 під-кроків за кадр
  render();
  requestAnimationFrame(loop);
})();