Lattice Boltzmann Fluid Simulation

const W = 200, H = 100; // grid dimensions const Q = 9; // D2Q9 velocity vectors: centre(0), cardinal(4), diagonal(4) const CX = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; const CY = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; // Equilibrium weights const WT = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; // Bounce-back pairs (direction i ↔ OPPOSITE[i]) const OPPOSITE = [0, 3, 4, 1, 2, 7, 8, 5, 6]; // Distribution function arrays: f[q * W * H + y * W + x] const f = new Float64Array(Q * W * H); // current distributions const fNew = new Float64Array(Q * W * H); // post-stream // Macroscopic fields const rho = new Float64Array(W * H); // density const ux = new Float64Array(W * H); // velocity x const uy = new Float64Array(W * H); // velocity y const wall = new Uint8Array(W * H); // 1 = solid cell

// Equilibrium distribution for given macroscopic quantities 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); } // Set uniform right-flowing field with density 1 const U0 = 0.1; // inlet velocity (lattice units) 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); } } } // Circular obstacle in the centre-left 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; // relaxation rate; viscosity ν = (1/omega - 0.5) / 3 // OMEGA=1.7 → ν ≈ 0.069 (moderate 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; // Compute macroscopic density and velocity 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 collision 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++) { // Destination cell after streaming in direction 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: reverse direction 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); // no — swap references f.set(fNew); fNew.fill(0); }

// Left inlet: zou-he velocity boundary (constant inflow) 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; } } // Right outlet: zero-gradient (copy f from second-to-last column) 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; } // Curl: ∂uy/∂x − ∂ux/∂y (finite differences) 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; // Map curl to blue–white–red 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(); }