🌊 Tutorial: WebGL Fluid Simulation with Metaballs

📚 Table of Contents

1. Introduction
2. What is a Metaball?
3. Code Structure
4. WebGL Shaders
5. SPH Physics Simulation
6. The Floating Boat
7. Smoke Particle System
8. Optimizations and Tips

1. Introduction

This tutorial explains how to create a real-time fluid simulation using WebGL and the metaball technique. Our project includes:

100 fluid particles simulated with SPH pressure forces
Metaball rendering in a fragment shader
A boat that floats on the fluid surface
Smoke particles rising from the chimney
Mouse interaction to push the fluid away

Prerequisites: Basic knowledge of JavaScript, HTML5 Canvas, and GLSL (shader language) concepts.

2. What is a Metaball?

Metaballs are a rendering technique where multiple shapes "merge" together organically, like water droplets joining together.

The Mathematical Principle

For each pixel on the screen, we calculate a sum of influences from all particles:

sum = Σ (r² / d²)

where r = effective radius, d = distance to particle center

If this sum exceeds a threshold (usually 1.0), the pixel is considered to be "inside" the fluid.

Shader Code

// For each particle, add its influence
float sum = 0.0;
for(int i = 0; i < NUM_PARTICLES; i++) {
  vec2 p = u_particles[i];
  float d2 = dot(pos - p, pos - p);  // squared distance
  sum += (u_radius * u_radius) / (d2 + 0.0001);
}

// If sum > 1.0, we're inside the fluid
float t = step(1.0, sum);

Tip: We use d² instead of d to avoid an expensive square root calculation!

3. Code Structure

3.1 Simulation Parameters

const NUM_PARTICLES = 100;    // Number of particles
const RADIUS = 0.04;          // Metaball size
const BUBBLE_RADIUS = 0.4;    // Container bubble radius
const GRAVITY = 0.0008;       // Downward gravity
const PRESSURE_RADIUS = 0.08; // Interaction radius between particles
const STIFFNESS = 0.004;      // Pressure stiffness
const STIFFNESS_NEAR = 0.01;  // Near pressure stiffness

3.2 Particle Initialization

// Initial distribution in a circle
for(let i = 0; i < NUM_PARTICLES; i++){
  let angle = Math.random() * Math.PI * 2;
  let r = Math.random() * 0.3;
  particles.push({
    x: r * Math.cos(angle),
    y: r * Math.sin(angle),
    vx: 0,  // X velocity
    vy: 0   // Y velocity
  });
}

3.3 Animation Loop

function animate() {
  // 1. Update particle physics
  // 2. Calculate SPH pressure forces
  // 3. Apply boundary constraints
  // 4. Calculate boat position
  // 5. Update smoke particles
  // 6. Send data to shaders
  // 7. Draw
  
  requestAnimationFrame(animate);
}

4. WebGL Shaders

4.1 Vertex Shader

The vertex shader is very simple - it draws a fullscreen quad (two triangles covering the entire screen):

attribute vec2 a_position;
varying vec2 v_uv;

void main() {
  // Convert position [-1,1] to UV [0,1]
  v_uv = a_position * 0.5 + 0.5;
  gl_Position = vec4(a_position, 0.0, 1.0);
}

4.2 Fragment Shader

The fragment shader does all the rendering pixel by pixel:

precision highp float;
varying vec2 v_uv;
uniform vec2 u_resolution;
uniform vec2 u_particles[100];
uniform float u_radius;

void main() {
  // Convert UV to centered coordinates [-1, 1]
  vec2 pos = v_uv * 2.0 - 1.0;
  
  // Correct aspect ratio
  float aspect = u_resolution.x / u_resolution.y;
  pos.x *= aspect;
  
  // Calculate metaball sum
  float sum = 0.0;
  for(int i = 0; i < 100; i++) {
    vec2 p = u_particles[i];
    p.x *= aspect;
    float d2 = dot(pos - p, pos - p);
    sum += (u_radius * u_radius) / (d2 + 0.0001);
  }
  
  // Render: blue if inside fluid, background otherwise
  float t = step(1.0, sum);
  vec3 fluidColor = vec3(0.0, 0.8, 1.0);
  vec3 bgColor = vec3(0.53, 0.81, 0.92);
  vec3 col = mix(bgColor, fluidColor, t);
  
  gl_FragColor = vec4(col, 1.0);
}

WebGL 1.0 Warning: Loops with u_particles[i] (dynamic indexing) don't work on all devices. Our code "unrolls" the loop manually for compatibility.

5. SPH Physics Simulation

SPH (Smoothed Particle Hydrodynamics) is a method for simulating fluids with particles. Here's how it works:

5.1 Applied Forces

// Gravity
p.vy -= GRAVITY;

// Random nudge (natural wave motion)
if(Math.random() < 0.05) {
  p.vx += (Math.random() - 0.5) * 0.02;
  p.vy += (Math.random() - 0.5) * 0.015;
}

// Mouse repulsion
let dx = p.x - mouseX;
let dy = p.y - mouseY;
let dist2 = dx * dx + dy * dy;
if(dist2 < 0.05) {
  p.vx += dx * 0.02;
  p.vy += dy * 0.02;
}

5.2 SPH Pressure Forces

Particles repel each other to prevent compression:

// Step 1: Calculate local density
let density = 0, nearDensity = 0;
for(let j = 0; j < particles.length; j++) {
  if(i === j) continue;
  let d = distance(pi, pj);
  if(d < PRESSURE_RADIUS) {
    let q = 1 - d / PRESSURE_RADIUS;  // [0, 1]
    density += q * q;
    nearDensity += q * q * q;
  }
}

// Step 2: Apply repulsion force
let pressure = STIFFNESS * density;
let pressureNear = STIFFNESS_NEAR * nearDensity;

// Push particles apart
let force = pressure * q + pressureNear * q * q;
pj.x += dx * force * 0.5;
pj.y += dy * force * 0.5;
pi.x -= dx * force * 0.5;
pi.y -= dy * force * 0.5;

Why two pressures?

pressure: Standard force that prevents compression
pressureNear: Additional force for very close particles (prevents "clumping")

5.3 Circular Boundary Constraint

// Check if particle exits the circle
let px = p.x * aspect;
let py = p.y;
let dist = Math.sqrt(px * px + py * py);
let radiusScaled = BUBBLE_RADIUS * aspect;

if(dist > radiusScaled) {
  // Bring particle back to the edge
  let nx = px / dist;
  let ny = py / dist;
  p.x = (nx * radiusScaled) / aspect;
  p.y = ny * radiusScaled;
  
  // Bounce (invert velocity)
  p.vx *= -0.3;
  p.vy *= -0.3;
}

6. The Floating Boat

6.1 Finding the Fluid Surface

We use a binary search to find where the fluid surface (sum = 1.0) is located:

function getSurfaceY(atX) {
  let low = -0.5;   // Inside the fluid
  let high = 0.5;   // Above the fluid
  
  // Binary search: 20 iterations
  for(let i = 0; i < 20; i++) {
    let mid = (low + high) / 2;
    let sum = getMetaballSum(atX, mid);
    
    if(sum > 1.0) {
      low = mid;   // Still inside fluid
    } else {
      high = mid;  // Outside fluid
    }
  }
  return (low + high) / 2;
}

6.2 Calculating the Tilt

// Surface height at left and right of the boat
let leftY = getSurfaceY(boatX - 0.06);
let rightY = getSurfaceY(boatX + 0.06);

// Angle = arctan(difference / distance)
let boatAngle = Math.atan2(rightY - leftY, 0.12);

6.3 Drawing the Boat (Shader)

// Transform coordinates to boat local space
vec2 relPos = pos - boatPos;
relPos = rotate(relPos, -u_boatAngle);

// Hull (trapezoid)
float widthAtY = boatWidth * (1.0 - (hullTop - relPos.y) / boatHeight * 0.5);
if(relPos.y > hullBottom && relPos.y < hullTop && abs(relPos.x) < widthAtY) {
  col = vec3(1.0, 1.0, 1.0);  // White
}

// Cabin (rectangle)
if(relPos.y > hullTop && relPos.y < hullTop + cabinHeight && abs(relPos.x) < cabinWidth) {
  col = vec3(1.0, 1.0, 1.0);
}

// Porthole (circle)
float hublotDist = length(relPos - hublotPos);
if(hublotDist < hublotRadius) {
  col = vec3(0.3, 0.6, 0.8);  // Light blue
}

7. Smoke Particle System

7.1 Parameters

const NUM_SMOKE = 30;           // Maximum number of particles
const SMOKE_SPAWN_RATE = 0.15;  // Spawn probability per frame
const SMOKE_BUOYANCY = 0.00008; // Upward lift
const SMOKE_FADE = 0.008;       // Alpha decrease per frame
const SMOKE_GROW = 0.0003;      // Size increase per frame

7.2 Creating a New Smoke Particle

// Chimney position in world coordinates
let cosA = Math.cos(boatAngle);
let sinA = Math.sin(boatAngle);
let chimneyWorldX = boatX + (chimneyLocalX * cosA - chimneyLocalY * sinA);
let chimneyWorldY = surfaceY + (chimneyLocalX * sinA + chimneyLocalY * cosA);

// Spawn with some random variation
smokeParticles.push({
  x: chimneyWorldX + (Math.random() - 0.5) * 0.005,
  y: chimneyWorldY + Math.random() * 0.005,
  vx: -0.002,     // Initial velocity leftward
  vy: 0.004,      // Initial velocity upward
  alpha: 0.8,     // Initial opacity
  size: 0.012     // Initial size
});

7.3 Updating Particles

for(let s of smokeParticles) {
  s.vy += SMOKE_BUOYANCY;   // Buoyancy
  s.vx -= 0.00005;          // Wind leftward
  s.x += s.vx;
  s.y += s.vy;
  s.alpha -= SMOKE_FADE;    // Gradual fade
  s.size += SMOKE_GROW;     // Expansion
}

7.4 Smoke Rendering (Shader)

// Soft circle with smoothstep
float d = length(pos - smokePos);
float puff = smoothstep(smokeSize, smokeSize * 0.3, d);

// Blend with existing color
vec3 smokeColor = vec3(1.0, 1.0, 1.0);  // White
col = mix(col, smokeColor, puff * smoke.z * 0.7);

8. Optimizations and Tips

8.1 Avoid Expensive Calculations

Use d² instead of sqrt(d) when possible
Pre-calculate aspect ratio once
Add 0.0001 to denominators to avoid division by zero

8.2 WebGL Compatibility

// ❌ Doesn't work everywhere (dynamic indexing)
for(int i = 0; i < NUM; i++) {
  sum += u_particles[i];
}

// ✅ Solution: unroll the loop in JavaScript
let shaderCode = '';
for(let i = 0; i < NUM_PARTICLES; i++) {
  shaderCode += `sum += u_particles[${i}]; `;
}

8.3 Anti-aliasing

// Hard edge
float mask = step(radius, dist);

// Soft edge (anti-aliased)
float aa = 0.002;  // Gradient width
float mask = smoothstep(radius - aa, radius + aa, dist);

8.4 Performance Improvements

Going Further:

Use a spatial grid to limit pressure calculations to nearby neighbors
Move the simulation to the GPU using data textures
Reduce particle count on mobile devices

🎉 Conclusion

You now understand the basics of a fluid simulation with metaballs!

Key Concepts:

Metaballs create organic shapes that merge together
SPH physics simulates realistic fluid behavior
The fragment shader renders in parallel on the GPU
Binary search efficiently finds the fluid surface

Ideas for Improvement:

Add waves on the surface
Create reflections and refractions
Add multiple boats
Implement water flow (fountain)

← Back to Demo