Туторіал · CFD · Алгоритми · JavaScript

📅 Липень 2026 ⏱ ≈ 30 хв 🎯 Середній – Просунутий

Метод Больцмана на ґратці приблизно за 200 рядків JavaScript

Вам не потрібен солвер Нав'є-Стокса, рівняння Пуассона для тиску чи генератор сітки, щоб симулювати реальний потік рідини. Метод ґраткового Больцмана дає обтікання циліндра — з повноцінною вулицею вихорів Кармана — за допомогою жменьки локальних, надзвичайно паралелізованих операцій над масивами.

1. Ґратка D2Q9

«D2Q9» означає 2 просторові виміри, 9 дискретних швидкостей. Популяція часток може або залишатись на місці, або переходити до одного з 8 сусідів (4 осьових + 4 діагональних) щокроку. Кожен напрямок має вагу w[i], що використовується у формулі рівноваги:

// Набір швидкостей D2Q9: індекс 0 = спокій, 1-4 = осі, 5-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]; // протилежний напрямок для bounce-back

Ґраткові одиниці: тут усе в «ґраткових одиницях»: Δx = Δt = 1, швидкість звуку ґратки c_s = 1/√3. Фізична в'язкість відповідає часу релаксації τ через ν = c_s²(τ − 0.5).

2. Структура даних

Дев'ять плоских Float64Array (по одному на напрямок) зберігають розподіли всієї сітки. Булева маска позначає комірки перешкоди — тут коло, розміщене на третині шляху в область:

const NX = 200, NY = 80;
const N  = NX * NY;
const f  = Array.from({length: 9}, () => new Float64Array(N));
const feq = Array.from({length: 9}, () => new Float64Array(N));
const rho = new Float64Array(N).fill(1);
const ux  = new Float64Array(N);
const uy  = new Float64Array(N);
const obstacle = new Uint8Array(N);

const cx = NX / 5, cy = NY / 2, R = NY / 9;
for (let y = 0; y < NY; y++)
  for (let x = 0; x < NX; x++)
    if ((x - cx)**2 + (y - cy)**2 < R*R)
      obstacle[y * NX + x] = 1;

3. Макроскопічні змінні

Густина і швидкість напряму випливають з розподілів як нульовий та перший моменти швидкості — жодного додаткового рівняння розв'язувати не потрібно:

ρ = Σᵢ fᵢ u = (1/ρ)\cdotΣᵢ fᵢ\cdoteᵢ

function computeMacroscopic() {
  for (let n = 0; n < N; n++) {
    let r = 0, vx = 0, vy = 0;
    for (let i = 0; i < 9; i++) {
      const fi = f[i][n];
      r += fi; vx += fi * ex[i]; vy += fi * ey[i];
    }
    rho[n] = r;
    ux[n]  = vx / r;
    uy[n]  = vy / r;
  }
}

4. Рівноважний розподіл

Рівновага BGK — це усічене до другого порядку розкладення розподілу Максвелла-Больцмана за локальною швидкістю u:

fᵢ^eq = wᵢ\cdotρ\cdot[1 + 3(eᵢ\cdotu) + 4.5(eᵢ\cdotu)² - 1.5|u|²]

function computeEquilibrium() {
  for (let n = 0; n < N; n++) {
    const r = rho[n], vx = ux[n], vy = uy[n];
    const usq = vx * vx + vy * vy;
    for (let i = 0; i < 9; i++) {
      const eu = ex[i] * vx + ey[i] * vy;
      feq[i][n] = w[i] * r * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * usq);
    }
  }
}

5. Зіткнення (релаксація BGK)

Кожен розподіл релаксує до своєї рівноваги зі швидкістю 1/τ. Саме в цьому єдиному відніманні живе в'язкість — більше τ означає повільнішу релаксацію та в'язкішу рідину:

fᵢ \leftarrow fᵢ - (1/τ)\cdot(fᵢ - fᵢ^eq)

const tau = 0.6;  // τ > 0.5 потрібно для стабільності; ν = (τ-0.5)/3
function collide() {
  for (let i = 0; i < 9; i++)
    for (let n = 0; n < N; n++)
      if (!obstacle[n])
        f[i][n] -= (f[i][n] - feq[i][n]) / tau;
}

6. Стрімінг

Стрімінг зсуває кожен розподіл на один крок ґратки вздовж його власного напрямку — чисте копіювання масиву, а обгортання пропускається на краях області (обробляється окремо граничними умовами):

function stream() {
  for (let i = 0; i < 9; i++) {
    const src = f[i], dst = new Float64Array(N);
    for (let y = 0; y < NY; y++)
      for (let x = 0; x < NX; x++) {
        const xs = x - ex[i], ys = y - ey[i]; // беремо з сусіда вгору за потоком
        if (xs < 0 || xs >= NX || ys < 0 || ys >= NY) continue;
        dst[y * NX + x] = src[ys * NX + xs];
      }
    f[i] = dst;
  }
}

Про продуктивність: промислові коди LBM стрімлять на місці, використовуючи хитрий трюк з індексами (або два буфери, що міняються місцями для кожного напрямку), щоб уникнути повного виділення масиву, показаного тут — але формулювання pull-based вище найлегше читати і правильно реалізувати з першого разу.

7. Граничні умови

Три різні правила замикають область: вхід з фіксованою швидкістю (спрощений Zou-He), вихід з нульовим градієнтом та bounce-back для твердих стінок і циліндра — на сьогодні найпростіша умова прилипання в LBM: просто реверсувати вхідні популяції на місці:

function applyBoundaries() {
  // Bounce-back: вузли перешкоди відбивають кожен вхідний напрямок
  for (let n = 0; n < N; n++) {
    if (!obstacle[n]) continue;
    const tmp = new Float64Array(9);
    for (let i = 0; i < 9; i++) tmp[i] = f[i][n];
    for (let i = 0; i < 9; i++) f[i][n] = tmp[opp[i]];
  }

  // Лівий край (вхід): фіксуємо однорідну горизонтальну швидкість u0
  const u0 = 0.08; // тримайте значно нижче c_s/√3 ≈ 0.577 для стабільності
  for (let y = 0; y < NY; y++) {
    const n = y * NX;
    ux[n] = u0; uy[n] = 0;
    let s = 0;
    for (let i = 0; i < 9; i++) if (ex[i] <= 0) s += f[i][n];
    rho[n] = s / (1 - u0);
    for (let i = 0; i < 9; i++) {
      const eu = ex[i] * u0;
      f[i][n] = w[i] * rho[n] * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * u0 * u0);
    }
  }

  // Правий край (вихід): нульовий градієнт — копіюємо стовпчик перед ним
  for (let y = 0; y < NY; y++) {
    const n = y * NX + (NX - 1), nPrev = n - 1;
    for (let i = 0; i < 9; i++) f[i][n] = f[i][nPrev];
  }

  // Верх / низ: bounce-back (тверді стінки каналу)
  for (let x = 0; x < NX; x++) {
    for (const y of [0, NY - 1]) {
      const n = y * NX + x;
      const tmp = new Float64Array(9);
      for (let i = 0; i < 9; i++) tmp[i] = f[i][n];
      for (let i = 0; i < 9; i++) f[i][n] = tmp[opp[i]];
    }
  }
}

8. Рендеринг завихреності

Сиру швидкість важко читати візуально; завихреність (ротор швидкості) робить почергові вихори Кармана видимими як червоно-сині смуги. Проста центральна різниця для ротора та розбіжна червоно-біло-синя колірна карта роблять свою справу:

function renderVorticity(ctx, img) {
  const data = img.data;
  for (let y = 1; y < NY - 1; y++)
    for (let x = 1; x < NX - 1; x++) {
      const n = y * NX + x;
      const dvdx = uy[n + 1] - uy[n - 1];
      const dudy = ux[n + NX] - ux[n - NX];
      const curl = (dvdx - dudy) * 40; // масштаб для видимості
      const p = (y * NX + x) * 4;
      if (obstacle[n]) { data[p] = data[p+1] = data[p+2] = 40; }
      else if (curl > 0) { data[p] = clamp(curl * 255); data[p+1] = 20; data[p+2] = 20; }
      else                { data[p] = 20; data[p+1] = 20; data[p+2] = clamp(-curl * 255); }
      data[p+3] = 255;
    }
  ctx.putImageData(img, 0, 0);
}
function clamp(v) { return Math.max(0, Math.min(255, v)); }

9. Складаємо разом — головний цикл

Кожен кадр анімації запускає чотири кроки по черзі — макроскопічні змінні, зіткнення, стрімінг, границі — а потім рендерить результат. Це весь алгоритм LBM, менше 200 рядків, включно з налаштуванням і рендерингом:

function step() {
  computeMacroscopic();
  computeEquilibrium();
  collide();
  stream();
  applyBoundaries();
}

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

// Ініціалізація у спокої з невеликою правою швидкістю всюди
computeEquilibriumAt(0.05);
function loop() {
  for (let s = 0; s < 4; s++) step(); // кілька кроків LBM на кожен кадр рендеру
  renderVorticity(ctx, img);
  requestAnimationFrame(loop);
}
loop();

Читання результату: при Re ≈ UD/ν ≈ 100–200 (підберіть τ та u0), слід за циліндром має почати зривати почергові червоно-сині пари вихорів — вулицю вихорів Кармана, описану у статті Стаціонарне та нестаціонарне обтікання циліндра. Повну теорію того, чому FDM, FVM та LBM дають ту саму фізику через зовсім різні дискретизації, дивіться у статті Чисельне рішення Нав'є-Стокса: FDM vs FVM vs LBM.

▶ Демо наживо

Часті запитання

Чого я навчуся в цьому уроці?

Побудуйте робочий CFD-солвер Lattice-Boltzmann (D2Q9) з нуля приблизно за 200 рядків JavaScript: стрімінг, BGK-зіткнення, bounce-back границі та рендеринг завихреності на canvas.

Які теми розглядаються в цьому уроці?

Цей урок охоплює такі теми: Ґратка D2Q9, Структура даних, Макроскопічні змінні, Рівноважний розподіл, Зіткнення, Стрімінг, Граничні умови, Рендеринг завихреності.

Скільки часу займає цей урок?

Цей урок займає приблизно 30 хв.

Які попередні знання потрібні?

Це урок рівня «Середній – Просунутий» — окрема попередня підготовка, крім базового JavaScript, не потрібна.