Симуляція рідини методом Латтіс-Больцмана

const W = 200, H = 100; // розміри сітки const Q = 9; // вектори швидкості D2Q9: центр(0), осьові(4), діагональні(4) const CX = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; const CY = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; // Рівноважні ваги const WT = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; // Пари bounce-back (напрямок i ↔ OPPOSITE[i]) const OPPOSITE = [0, 3, 4, 1, 2, 7, 8, 5, 6]; // Масиви функцій розподілу: f[q * W * H + y * W + x] const f = new Float64Array(Q * W * H); // поточні розподіли const fNew = new Float64Array(Q * W * H); // після стрімінгу // Макроскопічні поля const rho = new Float64Array(W * H); // густина const ux = new Float64Array(W * H); // швидкість x const uy = new Float64Array(W * H); // швидкість y const wall = new Uint8Array(W * H); // 1 = тверда клітинка

// Рівноважний розподіл для заданих макроскопічних величин function feq(q, rho0, u0, v0) { const uc = CX[q] * u0 + CY[q] * v0; // dot(c, u) const u2 = u0 * u0 + v0 * v0; return WT[q] * rho0 * (1 + 3*uc + 4.5*uc*uc - 1.5*u2); } // Задаємо рівномірне поле з потоком праворуч і густиною 1 const U0 = 0.1; // швидкість на вході (ґраткові одиниці) for (let y = 0; y < H; y++) { for (let x = 0; x < W; x++) { const idx = y * W + x; rho[idx] = 1.0; ux[idx] = U0; uy[idx] = 0; for (let q = 0; q < Q; q++) { f[q * W * H + idx] = feq(q, 1.0, U0, 0); } } } // Кругла перешкода у центрі-ліворуч const cx = Math.floor(W * 0.25), cy = Math.floor(H * 0.5), r = 10; for (let y = 0; y < H; y++) { for (let x = 0; x < W; x++) { if ((x-cx)**2 + (y-cy)**2 <= r*r) wall[y*W+x] = 1; } }

const OMEGA = 1.7; // швидкість релаксації; вʼязкість ν = (1/omega - 0.5) / 3 // OMEGA=1.7 → ν ≈ 0.069 (помірне Re) function collide() { for (let y = 0; y < H; y++) { for (let x = 0; x < W; x++) { const idx = y * W + x; if (wall[idx]) continue; // Обчислюємо макроскопічну густину та швидкість let ρ = 0, vx = 0, vy = 0; for (let q = 0; q < Q; q++) { const fi = f[q * W * H + idx]; ρ += fi; vx += fi * CX[q]; vy += fi * CY[q]; } rho[idx] = ρ; ux[idx] = vx / ρ; uy[idx] = vy / ρ; // Зіткнення BGK for (let q = 0; q < Q; q++) { const fi = f[q * W * H + idx]; const feqi = feq(q, ρ, vx/ρ, vy/ρ); f[q * W * H + idx] = fi - OMEGA * (fi - feqi); } } } }

function stream() { for (let y = 0; y < H; y++) { for (let x = 0; x < W; x++) { const idx = y * W + x; for (let q = 0; q < Q; q++) { // Клітинка призначення після стрімінгу в напрямку q const nx = x + CX[q]; const ny = y + CY[q]; if (nx < 0 || nx >= W || ny < 0 || ny >= H) continue; const nidx = ny * W + nx; if (wall[nidx]) { // Bounce-back: обертаємо напрямок fNew[OPPOSITE[q] * W * H + idx] += f[q * W * H + idx]; } else { fNew[q * W * H + nidx] = f[q * W * H + idx]; } } } } fNew.copyWithin(0, 0); // ні — міняємо посилання f.set(fNew); fNew.fill(0); }

// Лівий вхід: гранична умова швидкості Зу-Хе (сталий приплив) function inletBC() { for (let y = 1; y < H - 1; y++) { const idx = y * W + 0; ux[idx] = U0; uy[idx] = 0; const ρ = (f[0*W*H+idx] + f[2*W*H+idx] + f[4*W*H+idx] + 2*(f[3*W*H+idx] + f[6*W*H+idx] + f[7*W*H+idx])) / (1 - U0); f[1*W*H+idx] = f[3*W*H+idx] + (2/3)*ρ*U0; f[5*W*H+idx] = f[7*W*H+idx] - 0.5*(f[2*W*H+idx]-f[4*W*H+idx]) + (1/6)*ρ*U0; f[8*W*H+idx] = f[6*W*H+idx] + 0.5*(f[2*W*H+idx]-f[4*W*H+idx]) + (1/6)*ρ*U0; } } // Правий вихід: нульовий градієнт (копіюємо f з передостаннього стовпця) function outletBC() { for (let y = 0; y < H; y++) { const idx = y * W + W - 1; for (let q = 0; q < Q; q++) { f[q * W * H + idx] = f[q * W * H + (y * W + W - 2)]; } } }

const canvas = document.getElementById('lbm'); const ctx2d = canvas.getContext('2d'); const imageData = ctx2d.createImageData(W, H); function render() { for (let y = 0; y < H; y++) { for (let x = 0; x < W; x++) { const idx = y * W + x; const p = (y * W + x) * 4; if (wall[idx]) { imageData.data[p] = 40; imageData.data[p+1] = 40; imageData.data[p+2] = 40; imageData.data[p+3] = 255; continue; } // Завихреність: ∂uy/∂x − ∂ux/∂y (скінченні різниці) const l = idx - 1, r = idx + 1; const b = idx - W, t = idx + W; const curl = (uy[r] - uy[l]) * 0.5 - (ux[t] - ux[b]) * 0.5; // Відображаємо завихреність у синій–білий–червоний const c = Math.max(-1, Math.min(1, curl * 10)); const R = c > 0 ? Math.floor(c * 255) : 0; const B = c < 0 ? Math.floor(-c * 255) : 0; const G = 0; imageData.data[p] = R; // R imageData.data[p+1] = G; // G imageData.data[p+2] = B; // B imageData.data[p+3] = 255; } } ctx2d.putImageData(imageData, 0, 0); } function loop() { requestAnimationFrame(loop); inletBC(); collide(); stream(); outletBC(); render(); }