Kelvin-Helmholtz instability

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

Kelvin-Helmholtz instability

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 ( $\gamma = 1.4$ ). The domain $[0,1]^2$ has a central layer ( $0.25 < y < 0.75$ ) with density $\rho = 2$ moving right, and outer layers with $\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.

References

Source code scenario.hpp

Powered by the samurai-euler solver engine - the code below is this case's scenario (initial & boundary conditions); the flux, time integration and mesh adaptation come from the shared engine.

scenario.hpp 72 lines

// 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)