Simulating Stratified Buoyancy in Matter.js

Simulating buoyancy in a stratified liquid with multiple density gradients requires calculating the buoyant force and hydrodynamic drag for an object as it intersects distinct fluid layers. Matter.js is a rigid-body physics engine that does not have a native fluid solver, but you can achieve realistic multi-density fluid stratification by hooking into the engine's update loop, segmenting fluid layers along the Y-axis, determining the submerged area per layer, and applying corresponding upward and damping forces to the floating bodies.

Understanding Fluid Stratification

In a stratified liquid, fluids of different densities arrange themselves into stable horizontal layers, with the densest fluid at the bottom. According to Archimedes' principle, the buoyant force (\(F_b\)) exerted on a body submerged in a fluid is equal to the weight of the fluid displaced:

\[F_b = \rho \cdot V \cdot g\]

Where:

When a body spans multiple layers, the total buoyant force is the sum of the forces exerted by each layer based on the portion of the body submerged within that specific layer's vertical bounds.

Step-by-Step Implementation

1. Define Fluid Layers

Create a data structure that defines each horizontal fluid layer by its vertical boundaries, density, and viscosity (drag coefficient).

const fluidLayers = [
  { minY: 100, maxY: 250, density: 0.001, viscosity: 0.02 }, // Top layer (e.g., oil)
  { minY: 250, maxY: 400, density: 0.002, viscosity: 0.05 }, // Middle layer (e.g., water)
  { minY: 400, maxY: 600, density: 0.004, viscosity: 0.12 }  // Bottom layer (e.g., honey/syrup)
];

2. Calculate Submerged Heights

For rectangular bodies aligned with the axes, calculate the vertical overlap between the body's bounding box and each fluid layer. For arbitrary polygons, approximate the area by slicing the polygon across layer boundaries or by calculating the fraction of vertical intersection.

Here is an approximation method based on bounding box overlap:

function getSubmergedFraction(body, layer) {
  const bodyTop = body.bounds.min.y;
  const bodyBottom = body.bounds.max.y;
  const bodyHeight = bodyBottom - bodyTop;

  if (bodyHeight <= 0) return 0;

  // Find overlapping vertical region
  const overlapTop = Math.max(bodyTop, layer.minY);
  const overlapBottom = Math.min(bodyBottom, layer.maxY);

  if (overlapBottom > overlapTop) {
    return (overlapBottom - overlapTop) / bodyHeight;
  }

  return 0;
}

3. Compute and Apply Forces in beforeUpdate

Listen to the Matter.js beforeUpdate event. In this loop, evaluate each floating body against the fluid layers, compute the combined buoyant and drag forces, and apply them using Matter.Body.applyForce.

Matter.Events.on(engine, 'beforeUpdate', () => {
  const gravity = engine.gravity.y * engine.gravity.scale;

  floatingBodies.forEach(body => {
    fluidLayers.forEach(layer => {
      const submergedFraction = getSubmergedFraction(body, layer);

      if (submergedFraction > 0) {
        // Displaced mass = area * fraction * layer density
        const submergedArea = body.area * submergedFraction;
        const displacedMass = submergedArea * layer.density;

        // 1. Buoyant Force (acts upward against gravity)
        const buoyancyForceY = -displacedMass * gravity;

        // 2. Viscous Drag (opposes velocity)
        const dragForceX = -body.velocity.x * layer.viscosity * submergedFraction;
        const dragForceY = -body.velocity.y * layer.viscosity * submergedFraction;

        // Apply forces to the center of mass
        Matter.Body.applyForce(body, body.position, {
          x: dragForceX,
          y: buoyancyForceY + dragForceY
        });

        // Apply angular damping to simulate rotational resistance
        Matter.Body.setAngularVelocity(
          body,
          body.angularVelocity * (1 - layer.viscosity * submergedFraction)
        );
      }
    });
  });
});

Enhancing Realism

To increase physical accuracy:

  1. Center of Buoyancy: Instead of applying the buoyant force directly to body.position (the center of mass), compute the centroid of the submerged portion (the center of buoyancy) for each layer. Applying forces to the center of buoyancy creates natural righting moments and rotational stability.
  2. Polygon Clipping: For complex shapes, replace the bounding-box fraction with a polygon clipping algorithm (such as Sutherland-Hodgman) to clip the body's vertices to each layer's rectangular bounds, providing exact displaced polygon areas.
  3. Quadratic Drag: For high-velocity motion, apply quadratic drag (\(F_d = \frac{1}{2} \rho v^2 C_d A\)) alongside linear drag to prevent objects from instantly penetrating deep fluid layers at high speeds.