Метод ґраткового Больцмана — гідродинаміка з кінетичної теорії

1. Кінетична теорія — рівняння Больцмана

Класична гідродинаміка (Нав'є-Стокс) оперує макроскопічними полями: густиною ρ, швидкістю u та тиском p. Кінетична теорія дає мезоскопічний погляд: рідина описується функцією розподілу f(x, ξ, t) — густиною частинок у точці x з мікроскопічною швидкістю ξ у момент часу t.

Неперервне рівняння перенесення Больцмана керує еволюцією f:

Ключова ідея LBM — це дискретизація і простору, і швидкості. Простір стає регулярною ґраткою; неперервний простір швидкостей зводиться до малого скінченного набору дозволених швидкостей e_i. Член зіткнень апроксимується простою моделлю BGK. Попри ці радикальні спрощення, розклад Чепмена-Енскога показує, що отримані рівняння відтворюють рівняння Нав'є-Стокса для нестисливої рідини з точністю до другого порядку за числом Маха.

2. Ґратка D2Q9 — набір швидкостей і ваги

Назва «D2Q9» кодує вимірність і кількість швидкостей: 2 виміри, 9 напрямків швидкості. Кожен вузол ґратки поширює густину вздовж 9 напрямків (включно зі спокоєм).

Ваги w_i підібрані так, щоб ґратка коректно відтворювала перший і другий моменти (масу й імпульс) рівноважного розподілу. Умова Σ w_i = 1 забезпечує збереження маси.

3. Оператор зіткнень BGK

Наближення Батнагара-Ґросса-Крука (BGK) замінює справжній інтеграл зіткнень релаксацією до локальної рівноваги. Кожна функція розподілу f_i релаксує зі швидкістю 1/τ:

Час релаксації τ керує в'язкістю: τ = 1 дає ν = 1/6 (помірно в'язка). Стійкість вимагає τ > 0.5. Коли τ → 0.5, в'язкість → 0 і симуляція стає нестійкою для високошвидкісних потоків.

4. Крок поширення

Після зіткнення кожен пост-зіткненнєвий розподіл f_i* переходить до сусіднього вузла ґратки в напрямку e_i:

Отже, повний цикл оновлення LBM за один крок часу такий:

1. Обчислити макроскопічні ρ та u з поточної f
2. Обчислити рівноважні розподіли f^eq
3. Застосувати зіткнення BGK: f* = f − (1/τ)(f − f^eq)
4. Поширити: зсунути кожну f_i* на один крок у напрямку e_i
5. Застосувати граничні умови

5. Відновлення макроскопічних величин — ρ, u з f

Макроскопічні величини рідини — це моменти функції розподілу. Густина маси ρ — нульовий момент; імпульс ρu — перший:

Саме це локальне обчислення — що потребує лише 9 f_i в одному вузлі — робить LBM добре паралелізовним. Кожен вузол незалежний під час кроку зіткнення; лише крок поширення потребує даних від сусідів.

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

Стінки без проковзування — відбиття назад (Bounce-Back)

Найпростіша гранична умова без проковзування: будь-який розподіл, що прибуває у вузол твердої стінки, відбивається назад у протилежному напрямку — наче частинка відскочила від стінки.

Рухомі стінки — Зу-Хе

Для притоку/відтоку з відомою швидкістю u або відомим тиском ρ схема Зу-Хе аналітично екстраполює невідомі розподіли. Для лівого вхідного отвору зі швидкістю u = (u_x, 0):

7. Повна реалізація D2Q9 на JavaScript

// ── Константи D2Q9 ────────────────────────────────────────────────
const NX = 200, NY = 100;  // розмір сітки
const TAU = 0.6;            // час релаксації  (ν = (τ−0.5)/3 ≈ 0.033)
const OMEGA = 1 / TAU;      // швидкість релаксації
const N = NX * NY;

// Вектори швидкості: 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]; // індекс протилежного напрямку

// Функції розподілу: f[i * N + x + y*NX]
const f    = new Float32Array(9 * N).fill(0);
const fTmp = new Float32Array(9 * N);
const solid = new Uint8Array(N); // 1 = перешкода

function idx(x, y) { return x + y * NX; }
function fi(i, x, y) { return i * N + idx(x, y); }

// ── Ініціалізація: однорідний потік u_x = u0 ─────────────────────────────
function init(u0 = 0.1) {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      for (let i = 0; i < 9; i++) {
        f[fi(i, x, y)] = feq(i, 1.0, u0, 0.0);
      }
    }
  }
  // Додаємо круглу перешкоду ближче до лівого краю
  const cx = NX / 4, cy = NY / 2, r = NY / 8;
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if ((x-cx)**2 + (y-cy)**2 < r*r) solid[idx(x,y)] = 1;
    }
  }
}

// ── Рівноважний розподіл ──────────────────────────────────────
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);
}

// ── Один крок часу LBM ──────────────────────────────────────────
function step(u0 = 0.1) {
  // 1. Зіткнення + поширення (схема витягування)
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if (solid[idx(x, y)]) continue;

      // Обчислюємо макроскопічні величини
      let rho = 0, ux = 0, uy = 0;
      for (let i = 0; i < 9; i++) {
        const fi_val = f[fi(i, x, y)];
        rho += fi_val;
        ux  += fi_val * EX[i];
        uy  += fi_val * EY[i];
      }
      ux /= rho; uy /= rho;

      // Зіткнення BGK, поширення витягуванням
      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) {
          fTmp[fi(i, x, y)] = feq(i, 1.0, u0, 0.0); // приток
          continue;
        }
        if (solid[idx(sx, sy)]) {
          // Відбиття назад від перешкоди
          fTmp[fi(i, x, y)] = f[fi(OPP[i], x, y)];
          continue;
        }
        const f_src = f[fi(i, sx, sy)];
        fTmp[fi(i, x, y)] = f_src - OMEGA * (f_src - feq(i, rho, ux, uy));
      }
    }
  }
  fTmp.copyTo ? fTmp.copyTo(f) : f.set(fTmp);
}

// ── Рендер: колір за величиною швидкості ─────────────────────────
function render(ctx, canvas) {
  const img = ctx.createImageData(NX, NY);
  const d   = img.data;
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const p = idx(x, y);
      if (solid[p]) { d[p*4]   = 200; d[p*4+1] = 200; d[p*4+2] = 200; d[p*4+3] = 255; continue; }
      let ux = 0, uy = 0, rho = 0;
      for (let i = 0; i < 9; i++) {
        const v = f[fi(i, x, y)]; rho += v; ux += v*EX[i]; uy += v*EY[i];
      }
      ux /= rho; uy /= rho;
      const spd = Math.hypot(ux, uy) / 0.25; // нормалізуємо до [0,1]
      const t   = Math.min(spd, 1);
      d[p*4]   = (t * 56)  | 0;  // r
      d[p*4+1] = (t * 189) | 0;  // g (→ блакитний)
      d[p*4+2] = (t * 248) | 0;  // b
      d[p*4+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}

8. Застосування та розширення

• Доріжка вихорів Кармана: класичний еталонний тест обтікання циліндра — за Re ≈ 100 за перешкодою формується періодичний патерн зриву вихорів. LBM відтворює його природно.
• D3Q19 / D3Q27: 3D-варіанти ґратки з 19 або 27 швидкостями. Ідентична структура — лише більше напрямків швидкості й відповідно вищі витрати пам'яті.
• Багаточасова релаксація (MRT): замінює єдине τ релаксаційною матрицею у просторі моментів. Покращує стійкість за низької в'язкості (високих Re) і зменшує нефізичні зсувні артефакти.
• Багатофазний LBM: модель Шаня-Чена додає член міжчастинкової сили до кроку зіткнення, породжуючи фазове розділення, поверхневий натяг і динаміку крапель — повністю в межах підходу LBM.
• Тепловий LBM: друга функція розподілу відстежує енергію, уможливлюючи плавучу конвекцію (нестійкість Релея-Бенара).