samurai gallery

Kelvin-Helmholtz instability

Two shear layers roll up into the iconic billows; the mesh tracks the growing vortices.

ρ 0.887 to 2.18
3 mesh levels 46,279 cells 61 frames

When two fluid layers slide past each other, the shear interface is unstable: any small ripple grows and rolls up into a train of vortices - the Kelvin-Helmholtz billows seen in clouds, ocean currents and countless astrophysical flows. It is a favorite showcase for adaptive solvers because the action is confined to thin, evolving shear layers.

This case is a new scenario contributed to samurai-euler; the code shown is its definition (initial shear profile and the seeded perturbation).

Equations

Compressible Euler for an ideal gas (γ=1.4\gamma = 1.4). The domain [0,1]2[0,1]^2 has a central layer (0.25<y<0.750.25 < y < 0.75) with density ρ=2\rho = 2 moving right, and outer layers with ρ=1\rho = 1 moving left, at uniform pressure. A small vertical velocity localized at each interface seeds the instability.

Numerical method

  • Flux: HLLC approximate Riemann solver.
  • Time: explicit Euler.
  • Adaptation: multiresolution with a positivity-preserving prediction.

The image shows the density ρ\rho.

What to look for

Two rows of billows form along the interfaces and roll the light and heavy fluids into interlocking spirals. The adaptive mesh refines precisely along the braided shear layers, leaving the uniform interiors coarse.

Source

scenario.hpp · 72 lines

This case runs on the samurai-euler solver. The code below is its scenario: initial and boundary conditions. Fluxes, time integration and mesh adaptation come from the shared solver. 

// Copyright 2025 the samurai team
// SPDX-License-Identifier:  BSD-3-Clause

#pragma once

#include <cmath>
#include <numbers>

#include <samurai/bc.hpp>

#include "../variables.hpp"
#include "registry.hpp"

// Kelvin-Helmholtz instability: two horizontal layers in shear. A small
// vertical-velocity perturbation localized at the two interfaces grows into
// the characteristic rolled-up billows.
namespace test_case::kelvin_helmholtz
{
    constexpr double pi = std::numbers::pi;

    // Layer states
    double rho_in  = 2.0;  // inner layer (0.25 < y < 0.75)
    double rho_out = 1.0;  // outer layers
    double v_shear = 0.5;  // horizontal shear velocity (+/-)
    double p0      = 2.5;  // uniform pressure

    // Perturbation
    double amp   = 0.1;    // amplitude of the seeded vertical velocity
    double sigma = 0.05;   // interface thickness of the seed
    int    mode  = 2;      // number of billows (wavenumber = 2*mode)

    auto init_fn = [](auto& u, auto& cell)
    {
        auto c        = cell.center();
        const double x = c[0];
        const double y = c[1];

        const bool inner = (y > 0.25 && y < 0.75);
        const double rho = inner ? rho_in : rho_out;
        const double vx  = inner ? v_shear : -v_shear;

        const double seed = std::exp(-(y - 0.25) * (y - 0.25) / (2 * sigma * sigma))
                          + std::exp(-(y - 0.75) * (y - 0.75) / (2 * sigma * sigma));
        const double vy = amp * std::sin(2 * mode * pi * x) * seed;

        PrimState<2> state{
            rho,
            p0,
            xt::xtensor_fixed<double, xt::xshape<2>>{vx, vy}
        };
        u[cell] = prim2cons<2>(state);
    };

    void bc_fn(auto& u, double /*t*/)
    {
        samurai::make_bc<samurai::Neumann<1>>(u, 0., 0., 0., 0.);
    }

    template <std::size_t dim>
    auto box_fn()
    {
        xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
        xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1.};
        return samurai::Box<double, dim>(min_corner, max_corner);
    }
}

REGISTER_TEST_CASE(kelvin_helmholtz,
                   test_case::kelvin_helmholtz::box_fn,
                   test_case::kelvin_helmholtz::init_fn,
                   test_case::kelvin_helmholtz::bc_fn)