This interactive simulation demonstrates fluid dynamics by solving the incompressible Navier-Stokes equations in real-time. It uses Jos Stam's "Stable Fluids" algorithm, which is the industry standard for real-time smoke and fluid effects in games and visualizations.

The Navier-Stokes Equations

The Navier-Stokes equations describe how fluids (liquids and gases) move. For an incompressible fluid, they are:

Momentum Equation:

∂u/∂t = −(u·∇)u + ν∇²u − ∇p + f

Incompressibility Constraint:

∇·u = 0

Where:

• u = velocity field (how fast and in what direction the fluid moves at each point)
• ν = viscosity (how "thick" the fluid is — honey vs water)
• p = pressure field
• f = external forces (like your mouse input!)
• ∇² = Laplacian operator (measures how different a point is from its neighbors)
• ∇· = Divergence (measures if fluid is being created/destroyed)

The Stable Fluids Algorithm (3 Core Steps)

Jos Stam's brilliant insight was to split the equation into three sequential steps:

Step	Physics	What You See
1. DIFFUSE	Solve: ∂u/∂t = ν∇²u Viscosity spreads velocity; diffusion spreads density	Smoke "blurs" outward; thick fluids slow down
2. PROJECT	Enforce: ∇·u = 0 Solve Poisson equation for pressure, subtract gradient	Fluid doesn't compress; pushes around obstacles
3. ADVECT	Solve: ∂u/∂t = −(u·∇)u Backtrace particles along velocity field	Smoke moves with the flow; swirls form

Why "Stable"?

Traditional explicit methods "explode" when the time step is too large (the simulation becomes unstable). Stam's method uses implicit integration solved via Gauss-Seidel relaxation, which remains stable regardless of time step size. This is crucial for real-time applications where we can't afford tiny time steps.

Interactive Controls

• Click & Drag: Add density (smoke) and velocity (force) at cursor location. Fast movements = stronger forces.
• Run/Pause: Start or stop the simulation. Default state is paused.
• Reset: Clear all density and velocity.
• Step ◀ / ▶: Step backward/forward through simulation history frame by frame.
• Material: Quick presets for Smoke (low viscosity), Water (medium), Honey (high viscosity).
• Viscosity (ν): How "thick" the fluid is. High viscosity = velocity decays quickly (like honey). Low = flows freely (like water).
• Diffusion: How fast the smoke spreads out. High = smoke disperses quickly. Low = smoke stays coherent.
• Time Step (dt): Simulation speed. Higher = faster but may lose detail.
• Show Velocity: Toggle velocity field arrows. Color indicates speed: green (slow) → yellow → red (fast).
• Color Schemes: Grayscale (smoke), Thermal (heat map), Ocean (water), Plasma (scientific).

Demo Scenarios

The simulation includes several pre-configured scenarios demonstrating real-world applications:

Scenario	Description	Real-World Application
Free Draw	Manual interaction — draw smoke with mouse	General fluid exploration, education
Wind Tunnel	Constant horizontal airflow from left edge. Add obstacles to see flow separation, vortex shedding, and wake turbulence.	CFD / Aircraft Design: Engineers use wind tunnels to analyze drag, lift, and aerodynamic properties.
Pipe Flow	Channel with walls, parabolic velocity profile (Poiseuille flow). Demonstrates laminar flow through confined spaces.	Medical / Blood Flow: Simulates blood flow through arteries and veins. Used in cardiovascular research.
Fire	Upward buoyancy force based on density. Heat source at bottom creates rising flames with turbulent behavior.	Games / Movies: Real-time fire, smoke, and explosion effects in video games and visual effects.

Obstacle Types

Solid obstacles block fluid flow (velocity = 0 inside). Available shapes:

Obstacle	Shape	Flow Phenomena
Circle	Circular cylinder (radius 15 cells)	Classic von Kármán vortex street — alternating vortices shed behind the obstacle.
Rectangle	Rectangular block (20×16 cells)	Sharp corners cause strong flow separation and recirculation zones.
Airfoil (NACA 2412)	Realistic wing profile — 64 cells long, 12° angle of attack. NACA 2412: 2% camber at 40% chord, 12% thickness (Cessna 172 style)	Demonstrates lift generation — flow accelerates over curved upper surface (low pressure), decelerates under flat lower surface (high pressure).

Wind Speed Control

The Wind slider (0.5 - 5.0) controls inflow velocity for Wind Tunnel and Pipe scenarios. Higher speeds create more dramatic vortices but may require lower time steps for stability.

Understanding the Physics

Viscosity (ν): This is the "internal friction" of the fluid. At the molecular level, it's how much momentum is transferred between adjacent fluid layers. High viscosity fluids (honey, molasses) resist flow. Low viscosity fluids (water, air) flow easily. In the simulation, viscosity causes the velocity field to "blur" over time.

Diffusion: This represents how quickly "stuff" in the fluid spreads out due to random molecular motion. For smoke, high diffusion means it disperses into the air quickly. For dye in water, low diffusion means it stays concentrated longer.

Incompressibility (∇·u = 0): This constraint says that fluid cannot be created or destroyed at any point — what flows in must flow out. The PROJECT step enforces this by computing a pressure field that "pushes back" against any compression. This is why fluid flows around obstacles rather than through them.

Advection: This is the transport of quantities (density, velocity itself) along the flow. The semi-Lagrangian method traces particles backwards in time to find where they came from, then samples the field there. This "backtracking" approach is unconditionally stable — it won't explode even with large time steps.

The Gauss-Seidel Solver

Both the diffusion and pressure steps require solving large systems of linear equations. We use the Gauss-Seidel relaxation method:

x[i,j] = (x₀[i,j] + a·(x[i-1,j] + x[i+1,j] + x[i,j-1] + x[i,j+1])) / c

We iterate this formula 16 times per step. More iterations = more accurate but slower. The beauty of implicit methods is that they're stable even with few iterations.

Grid-Based (Eulerian) Approach

This simulation uses a fixed 128×128 grid where each cell stores:

• Density: How much "smoke" is at this location
• Velocity (Vx, Vy): How fast and in what direction fluid is moving

This is called the Eulerian approach (named after mathematician Leonhard Euler). The alternative is Lagrangian (tracking individual particles), which is used in SPH simulations. Grid-based methods are simpler to implement and understand, which is why they're standard for educational purposes.

Boundary Conditions

The simulation uses reflective boundaries:

• Velocity: Perpendicular component is negated at walls (fluid bounces off)
• Density: Simply copied to boundary cells (smoke doesn't disappear at walls)
• Corners: Average of neighboring boundary cells

Detailed FDM Calculation Process

The Finite Difference Method (FDM) approximates continuous derivatives using discrete grid values. Here's exactly how each step works:

Grid Setup

We use a 128×128 grid with 2 extra boundary cells on each side (130×130 total). Each cell is indexed as:

IX(i, j) = i + (N + 2) × j     where N = 128

Cell (i,j) has neighbors at (i±1, j) and (i, j±1). The grid stores:

• density[IX(i,j)] — smoke amount at cell (i,j)
• Vx[IX(i,j)] — x-component of velocity at cell (i,j)
• Vy[IX(i,j)] — y-component of velocity at cell (i,j)

STEP 1: DIFFUSE — Solving ∂u/∂t = ν∇²u

The Continuous Equation:

∂u/∂t = ν × ∇²u = ν × (∂²u/∂x² + ∂²u/∂y²)

FDM Discretization of Laplacian (∇²):

The second derivative is approximated using the 5-point stencil:

∇²u[i,j] ≈ (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] − 4×u[i,j]) / h²

where h = 1/N is the grid spacing.

Implicit Time Integration:

Instead of explicit (u_new = u_old + dt × ν × ∇²u_old), we use implicit:

u_new = u_old + dt × ν × ∇²u_new

Rearranging with a = dt × ν × N²:

(1 + 4a) × x[i,j] = x₀[i,j] + a × (x[i+1,j] + x[i-1,j] + x[i,j+1] + x[i,j-1])

Gauss-Seidel Iteration (16 passes):

for (k = 0; k < 16; k++) {           // 16 iterations
    for (j = 1; j <= N; j++) {
        for (i = 1; i <= N; i++) {
            x[IX(i,j)] = (x0[IX(i,j)] + a × (
                x[IX(i+1,j)] + x[IX(i-1,j)] +
                x[IX(i,j+1)] + x[IX(i,j-1)]
            )) / (1 + 4×a);
        }
    }
    set_bnd(b, x);                    // Apply boundary conditions
}

STEP 2: PROJECT — Enforcing ∇·u = 0

This step ensures mass conservation (incompressibility). Fluid cannot be created or destroyed.

Sub-step 2a: Calculate Divergence

Divergence measures how much fluid is being "created" at each point:

div[i,j] = −0.5 × h × ((Vx[i+1,j] − Vx[i−1,j]) + (Vy[i,j+1] − Vy[i,j−1]))

This uses central differences for first derivatives.

Sub-step 2b: Solve Poisson Equation for Pressure

We need to find pressure p such that: ∇²p = div

p[i,j] = (div[i,j] + p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1]) / 4

Again solved with 16 Gauss-Seidel iterations.

Sub-step 2c: Subtract Pressure Gradient

Remove the compressive component from velocity:

Vx[i,j] -= 0.5 × N × (p[i+1,j] − p[i−1,j])
Vy[i,j] -= 0.5 × N × (p[i,j+1] − p[i,j−1])

STEP 3: ADVECT — Semi-Lagrangian Backtracing

This step moves quantities (density, velocity) along the flow. Instead of pushing particles forward, we trace backwards to find where each value came from.

Backtrace Calculation:

x_src = i − dt × N × Vx[i,j]
y_src = j − dt × N × Vy[i,j]

This finds the "source position" — where the particle at (i,j) came from.

Bilinear Interpolation:

Since (x_src, y_src) usually falls between grid points, we interpolate:

i0 = floor(x_src);  i1 = i0 + 1;       // neighboring grid cells
j0 = floor(y_src);  j1 = j0 + 1;

s1 = x_src − i0;    s0 = 1 − s1;       // interpolation weights
t1 = y_src − j0;    t0 = 1 − t1;

// Bilinear interpolation formula:
d_new[i,j] = s0 × (t0 × d_old[i0,j0] + t1 × d_old[i0,j1])
           + s1 × (t0 × d_old[i1,j0] + t1 × d_old[i1,j1]);

d[i,j] = s0×t0×d[i0,j0] + s0×t1×d[i0,j1] + s1×t0×d[i1,j0] + s1×t1×d[i1,j1]

Display & Rendering Logic

The visualization converts the numerical grid data into pixels on the canvas.

Density Rendering

For each grid cell (i,j), we map density to pixel brightness:

for (j = 0; j < N; j++) {
    for (i = 0; i < N; i++) {
        d = density[IX(i+1, j+1)];     // Get density (skip boundary)
        color = getColor(d);            // Map to RGB based on scheme
        
        // Fill SCALE×SCALE pixels (4×4 = 16 pixels per cell)
        for (sy = 0; sy < SCALE; sy++) {
            for (sx = 0; sx < SCALE; sx++) {
                px = i × SCALE + sx;
                py = j × SCALE + sy;
                idx = (py × canvasWidth + px) × 4;  // RGBA index
                
                imageData[idx]     = color.r;
                imageData[idx + 1] = color.g;
                imageData[idx + 2] = color.b;
                imageData[idx + 3] = 255;           // Alpha
            }
        }
    }
}
ctx.putImageData(imageData, 0, 0);  // Draw all pixels at once

Color Schemes

Scheme	Mapping (density d → RGB)
Grayscale	R = G = B = d
Thermal	Black → Red → Yellow → White d<85: (3d,0,0), d<170: (255,3(d-85),0), else: (255,255,3(d-170))
Ocean	Black → Blue → Cyan → White d<85: (0,0,3d), d<170: (0,3(d-85),255), else: (3(d-170),255,255)
Plasma	Purple → Pink → Orange → Yellow t=d/255: (255×min(1,0.5+t), 255×t², 255×(1-0.7t))

Velocity Vector Rendering

When "Show Velocity" is enabled, arrows are drawn on top of the density field:

for (j = 1; j < N; j += 12) {          // Every 12th cell
    for (i = 1; i < N; i += 12) {
        vx = Vx[IX(i,j)];
        vy = Vy[IX(i,j)];
        mag = sqrt(vx² + vy²);
        
        if (mag < 0.0005) continue;       // Skip tiny vectors
        
        // Normalize direction
        nx = vx / mag;
        ny = vy / mag;
        
        // Arrow length (8-40 pixels)
        len = min(40, max(8, mag × 50));
        
        // Color by magnitude (green → yellow → red)
        intensity = min(1, mag × 5);
        r = 255 × intensity;
        g = 255 × (1 − 0.5 × intensity);
        
        // Draw shaft + arrowhead
        drawLine(i×SCALE, j×SCALE, endX, endY);
        drawArrowhead(endX, endY, angle);
    }
}

Complete Frame Pipeline

#	Stage	Description
1	Read Input	If mouse is dragging, add density and velocity at cursor location
2	Velocity Step	DIFFUSE(Vx) → DIFFUSE(Vy) → PROJECT → ADVECT(Vx) → ADVECT(Vy) → PROJECT
3	Density Step	DIFFUSE(density) → ADVECT(density)
4	Fade Density	Multiply all density by 0.995 (smoke slowly dissipates)
5	Render Density	Map density grid to canvas pixels using color scheme
6	Render Velocity	If enabled, draw velocity arrows on top
7	Next Frame	requestAnimationFrame(loop) schedules next iteration

Performance Considerations

• Float32Array: We use typed arrays instead of regular JavaScript arrays for ~10x better performance in tight loops.
• 1D Array: The 2D grid is stored as a flattened 1D array for better cache locality.
• Scaled Rendering: We compute on a 128×128 grid but render at 512×512 (4x scale) for visual quality without computational cost.
• ImageData API: Direct pixel manipulation is faster than drawing individual shapes.

Applications of Navier-Stokes

The Navier-Stokes equations are fundamental to understanding fluid behavior across many fields. Try the demo scenarios above to explore these applications:

• Weather Prediction: Atmospheric models solve NS equations on massive scales to forecast weather patterns, hurricanes, and climate.
• Aircraft Design (→ Wind Tunnel + Airfoil): Computational Fluid Dynamics (CFD) simulates airflow around wings, fuselages, and control surfaces. Our NACA 2412 airfoil demonstrates lift generation.
• Video Games (→ Fire scenario): Real-time smoke, fire, and water effects. The buoyancy-driven fire simulation shows how game engines create realistic flames.
• Movies: Visual effects for oceans, explosions, smoke — Hollywood uses similar solvers but at much higher resolution.
• Medical (→ Pipe Flow): Blood flow simulation in arteries. The parabolic velocity profile (Poiseuille flow) matches real hemodynamics in vessels.
• Engineering: HVAC systems, car aerodynamics, ship hull design, chemical mixing — anywhere fluids flow!

The NACA 2412 Airfoil

The airfoil obstacle uses the NACA 4-digit airfoil system developed by NASA's predecessor (NACA) in the 1930s. The designation "2412" encodes:

• 2 = Maximum camber is 2% of chord length
• 4 = Maximum camber occurs at 40% of chord from leading edge
• 12 = Maximum thickness is 12% of chord length

This profile is used on aircraft like the Cessna 172 — the most produced aircraft in history. At 12° angle of attack (as shown in simulation), the airfoil generates significant lift but approaches the stall angle where flow separates from the upper surface.

The NACA thickness formula:

y_t = 5t × (0.2969√x − 0.1260x − 0.3516x² + 0.2843x³ − 0.1015x⁴)

where t = thickness ratio and x = position along chord (0 to 1).

The Millennium Prize Problem

Fun fact: Proving whether smooth solutions to the 3D Navier-Stokes equations always exist is one of the seven Millennium Prize Problems in mathematics, with a $1 million prize! We can simulate fluids numerically, but mathematically proving the equations are well-behaved in all cases remains unsolved.

Tips for Experimentation

• Wind Tunnel + Airfoil: Select both and watch vortices form behind the wing. Enable velocity vectors to see the flow field.
• Wind Tunnel + Circle: Classic von Kármán vortex street — alternating vortices shed from each side.
• Fire Mode: Watch turbulent flames rise. The thermal color scheme makes temperature visible.
• Pipe Flow: Observe the parabolic velocity profile — fastest in center, zero at walls (no-slip condition).
• Create Swirls: In Free Draw mode, make quick circular motions to generate vortices.
• Compare Materials: Switch between Smoke, Water, and Honey to see how viscosity affects flow decay.
• Step Mode: Use Step ▶ button to advance frame-by-frame and analyze flow patterns in detail.
• Velocity Vectors: Enable to see the actual flow field. Colors show speed: green (slow) → red (fast).
• Adjust Wind Speed: Higher wind creates more dramatic vortices but may need smaller time step for stability.