Billow

A shear layer picking its wavelength

Abstract

Two streams of fluid sliding past one another produce, at the surface between them, the most familiar instability in fluid mechanics: a row of rolling vortices, the Kelvin-Helmholtz billows. This companion sets out what the classical theory says about which wavelength wins, why the fastest disturbance is always two-dimensional even though the eventual mixing is not, and how the playground checks a live three-dimensional simulation against a hundred and sixty years of results. The central claim is testable while you watch: the seeded disturbance should grow at a rate the Rayleigh eigenvalue problem fixes in advance, and it does, to within about ten percent, until it saturates.

---

1. The instability

1.1 The mechanism

Take an interface across which the tangential velocity jumps. Perturb it into a small ripple. Fluid must speed up passing over a crest and slow down through a trough, and by Bernoulli's principle the pressure does the opposite: low over the crests, high in the troughs. That pressure difference pushes the crests further up and the troughs further down. The ripple feeds itself.

What follows is not a bigger ripple. The same shear that drives the growth also tilts the wave over, and past a certain amplitude the crest folds forward and wraps around itself into a vortex. The result is the row of curls that gives the phenomenon its name and this playground its title.

1.2 Where it appears

• Kelvin-Helmholtz clouds, the row of breaking waves along an atmospheric shear, where a cloud layer marks the interface and makes the instability visible.
• Saturn and Jupiter, at the boundaries between counter-rotating zonal bands.
• The solar corona, where billows have been imaged on the flanks of coronal mass ejections.
• Clear air turbulence, which is frequently the breakdown of exactly this instability where an aircraft cannot see it coming.
• Anywhere fluids mix, from an estuary to a stirred cup, since a folded interface is what mixing at scale actually looks like.

---

2. What the theory says

2.1 The Rayleigh equation

A parallel flow U(y) perturbed by a normal mode

psi(x, y, z, t) = phi(y) exp( i ( alpha x + beta z - alpha c t ) )

obeys, in the inviscid limit, the Rayleigh equation

(U - c) ( phi'' - k^2 phi ) - U'' phi = 0,      k^2 = alpha^2 + beta^2

with phi decaying at both edges. The eigenvalue c is a complex phase speed, and the mode grows in time at the rate sigma = alpha * Im(c).

Rayleigh proved in 1880 that a necessary condition for instability is an inflection point in the velocity profile, a place where U'' changes sign. A shear layer has one, right in the middle. A boundary layer on a flat plate does not, which is why it is stable to this mechanism and needs viscosity to become unstable at all.

2.2 The vortex sheet, and why it is not enough

The simplest model is a discontinuity: U = -U0 below the interface, +U0 above, with no transition at all. Kelvin and Helmholtz solved it and found

sigma = k * dU / 2

Every wavelength is unstable, and shorter wavelengths grow faster without bound. That is a broken model. It has no shortest scale, so it makes no prediction about what the instability looks like, and the initial value problem is ill posed.

The playground checks this limit rather than assuming it: the eigenvalue solver, handed a very long wave on a very thin layer, returns a phase speed within one percent of the vortex-sheet answer.

2.3 Finite thickness fixes it

Give the profile a thickness. Rayleigh's piecewise-linear layer, linear across a layer of half-thickness h and constant outside it, has a closed-form eigenvalue relation:

c^2 = ( U0 / (2 k h) )^2 * [ (1 - 2 k h)^2 - exp(-4 k h) ]

The bracket is negative, and therefore the layer unstable, only below a wavenumber that solves

(1 - 2 k h)^2 = exp(-4 k h)

whose nontrivial root is k h = 0.6392. Above it, nothing grows.

The smooth version, and the one this playground uses, is the hyperbolic tangent

U(y) = U0 * tanh( y / delta )

Michalke solved its Rayleigh problem numerically in 1964 and the results are the numbers on the dispersion chart:

The neutral point is exact, and worth stating plainly: a shear layer is stable to every wavelength shorter than about 6.3 times its half-thickness. A disturbance smaller than the layer cannot tell the two streams apart, so there is nothing for it to feed on.

The consequence is the whole point. The thickness of the layer sets the size of the billows. Without a thickness there is no wavelength, no scale, and no picture.

2.4 Squire's theorem

Disturbances need not travel along the flow. An oblique wave has a spanwise wavenumber beta as well as a streamwise one. Squire showed in 1933 that such a wave satisfies the same equation as a planar wave at the resultant wavenumber, with the same c. Only the prefactor changes:

sigma_3D(alpha, beta) = alpha * c_i( sqrt(alpha^2 + beta^2) )
                      = cos(theta) * sigma_2D( sqrt(alpha^2 + beta^2) )

Obliqueness is charged twice, and the second charge is much the larger. The cosine is at most one. And the resultant wavenumber is always bigger than the streamwise one, so the oblique wave is read off further along a dispersion curve that is already falling toward zero. A mode that is merely oblique can be outright stable while the planar mode of the same streamwise wavelength grows perfectly well.

The primary billow of a parallel shear flow is therefore two-dimensional, always. The playground verifies this numerically rather than asserting it: searching over all spanwise wavenumbers for the one that maximises growth returns exactly zero.

2.5 The tension, and why three dimensions matter anyway

If the fastest mode is planar, why is a real mixing layer not planar?

Because growing fastest is not the same as being the end of the story. The planar rollers form first. They then constitute a new base flow, one which is emphatically not parallel, and that new flow has instabilities of its own. Streamwise vortices, the ribs, appear in the braids between neighbouring rollers and stretch. Moser and Rogers traced this in direct numerical simulation in 1993 and found what everyone since has found: the transition to turbulence in a plane mixing layer is a secondary instability riding on a primary one, and it is where essentially all the mixing happens.

Two-dimensional billows are beautiful and mix almost nothing. They fold the interface, but folding is reversible in principle and the sheet stays a sheet. Three-dimensional breakdown is what turns folding into mixing.

---

3. The model

3.1 Domain and base flow

The simulated domain is the unit cube, periodic in all three directions.

Periodicity in the shear-normal direction has a consequence that has to be faced rather than hidden: a profile that climbs from one velocity to another cannot return, so a single shear layer cannot be periodic. There must be an even number of them. The base flow is therefore a pair of opposite-signed interfaces:

u_x(y) = U0 [ tanh( (y - 1/4) / delta ) - tanh( (y - 3/4) / delta ) - 1 ]

which is +U0 through the middle half of the box and -U0 in the outer half. Both interfaces roll. The two of them are part of the formulation, not a visual duplicate.

3.2 Seeding

The layer is forced the way a plane mixing layer is forced in experiment: a planar primary carrying most of the amplitude, plus a weaker pair of oblique waves at the same streamwise wavenumber, plus a little broadband noise to break the exact symmetry of the seeded modes.

Seeding only the oblique pair would be a poorly designed experiment. Squire's theorem guarantees it loses the race to whatever planar disturbance the round-off supplies, so the run would be decided by arithmetic noise rather than by anything intended. Seeding both puts the competition where it can be watched, and lets the spanwise control do what it should: not change how fast the layer rolls, but change how quickly the rollers break.

3.3 Discretisation

The equations solved are incompressible Navier-Stokes,

du/dt + (u . grad) u = -grad p + nu * laplacian(u),     div u = 0

by Stam's stable-fluids splitting: advect, add forces, diffuse, project, then carry the marker through the result. Three choices deserve comment.

One-sided, adjoint divergence and gradient. The usual GPU recipe takes both by central differences. Their composition is then the Laplacian on a stencil of width two, which is blind to the checkerboard mode, and the pressure solve cannot remove the part of the divergence that lives there. Measured on the CPU reference, roughly fifteen percent of a purely divergent flow survived that projection. Forward differences for the divergence and backward differences for the gradient compose into the compact seven-point Laplacian, which is exactly the operator the solve inverts, and the residual falls to zero.

Multigrid, not Jacobi. Jacobi damps an error of wavelength L by about 1 - (pi h / L)^2 per sweep. At sixty-four cells across, the wavelength of the seeded billow needs thousands of sweeps, and the two dozen a real-time budget allows remove under a fifth of it. The velocity field is then not incompressible, and a compressible flow does not conserve what it carries: the marker measurably drains away. Coarsening cures it, because a wavelength that is long on the fine grid is short on a coarse one. Four levels of V-cycle, from sixty-four cells down to four, reduce the residual by more than a factor of four per cycle.

MacCormack transport with a monotonicity limiter. A single semi-Lagrangian step is first order and heavily diffusive; after a few hundred steps the interface has smeared into a gradient several cells thick and the billow has no edge. Advecting forward and then backward estimates that error and subtracting it recovers second order. The correction can overshoot, which for a marker between zero and one is meaningless, so it is clamped to the range of the cells the backtrace actually interpolated between.

3.4 Diagnostics

Every diagnostic is an integral over planes of constant height, so the volume is collapsed on the GPU in two reduction passes and only the collapsed profile crosses to the CPU. This is not only cheaper than reading back a coarse copy of the volume, it is more accurate: the momentum thickness

theta = integral of (1/4) ( 1 - ( U_mean(y) / U0 )^2 ) dy

is a property of the horizontally averaged profile, and averaging in the other order produces a number that counts the billows themselves. For a single tanh interface the integral collapses exactly to delta / 2, which is what makes it checkable.

The growth rate is measured as half the slope of the log of the disturbance energy, by least squares over a sliding window. The base flow carries no cross-stream or spanwise velocity at all, so that energy is entirely the disturbance, and a normal mode of growth rate sigma carries energy growing as exp(2 sigma t).

---

4. Results

4.1 The prediction holds

Run the default preset and watch the growth chart. The measurement is flat for a fraction of an e-folding, while the seeded field settles onto the eigenfunction, and then climbs at a slope that matches twice the Rayleigh eigenvalue to within about ten percent. The dashed line is not fitted to the data: its slope is the eigenvalue and nothing else.

Then it stops matching, which it must. Linear theory is a statement about small disturbances. The point where the measurement bends away from the line is the end of the linear phase and the beginning of the billow.

4.2 What each preset shows

The fourth is the sharpest of them. Nothing about the picture suggests a threshold; the layer looks exactly as unstable as the others. It just does not roll, because the wavelength asked for is on the wrong side of alpha delta = 1.

---

5. Limitations

Resolution. Sixty-four cells across for the flow and ninety-six for the marker is coarse. Above a Reynolds number of roughly one thousand the numerical dissipation of the transport scheme exceeds the physical viscosity, so the Reynolds number in the panel stops being the one the fluid experiences.

Vorticity confinement is a numerical force with no counterpart in the Navier-Stokes equations. It restores roughly what the discretisation removes, which is what keeps the ribs alive long enough to see, but it injects energy rather than conserving it. It can be set to zero, and anything that disappears when it is was partly its doing.

The marker is passive. There is no density difference between the streams and no gravity, so there is no stratification and none of the Richardson-number physics that decides when the atmosphere actually produces billows. A stratified shear layer is stable wherever the gradient Richardson number exceeds one quarter, the Miles-Howard criterion, and this model cannot reproduce that.

Marker conservation. Semi-Lagrangian transport is not conservative. The marker total is displayed and drifts by well under a percent per unit time through the linear phase, faster once the flow is turbulent at grid scale. Nothing physical depends on it, since the marker is passive and every measured quantity comes from the velocity field.

Two interfaces, not one. The classical theory describes an isolated layer. Above a half-thickness of roughly 0.12 the two interfaces are no longer independent, and the predictions should be read as indicative.

Late times are a demonstration, not a prediction. Nonlinear shear-layer evolution is famously sensitive to the initial condition and the resolution. Linear growth here is a prediction and is checked against one. Vortex pairing, breakdown and the fine structure that follows are not, and should not be treated as reproducible without a convergence study: multiple resolutions, multiple time steps, a preserved random seed, and comparison of integral quantities rather than pictures.

---

References

• Kelvin, W. Thomson (1871). Hydrokinetic solutions and observations. Philosophical Magazine 42, 362.
• Helmholtz, H. (1868). On discontinuous movements of fluids. Philosophical Magazine 36, 337.
• Rayleigh, Lord (1880). On the stability, or instability, of certain fluid motions. Proceedings of the London Mathematical Society 11, 57.
• Squire, H. B. (1933). On the stability for three-dimensional disturbances of viscous fluid flow between parallel walls. Proceedings of the Royal Society A 142, 621.
• Michalke, A. (1964). On the inviscid instability of the hyperbolic tangent velocity profile. Journal of Fluid Mechanics 19, 543. https://doi.org/10.1017/S0022112064000908
• Drazin, P. G. and Reid, W. H. (2004). Hydrodynamic Stability, 2nd edition. Cambridge University Press.
• Moser, R. D. and Rogers, M. M. (1993). The three-dimensional evolution of a plane mixing layer: pairing and transition to turbulence. Journal of Fluid Mechanics 247, 275. https://ntrs.nasa.gov/citations/19930047801
• Brown, G. L. and Roshko, A. (1974). On density effects and large structure in turbulent mixing layers. Journal of Fluid Mechanics 64, 775.
• Stam, J. (1999). Stable fluids. SIGGRAPH 99, 121.
• Selle, A., Fedkiw, R., Kim, B., Liu, Y. and Rossignac, J. (2008). An unconditionally stable MacCormack method. Journal of Scientific Computing 35, 350.
• Schroeder, P. W. et al. (2019). On reference solutions and the sensitivity of the 2D Kelvin-Helmholtz instability problem. Computers and Mathematics with Applications 77, 1010. https://doi.org/10.1016/j.camwa.2018.10.030