Richtmyer-Meshkov

The Richtmyer-Meshkov instability develops when a shock wave interacts with a disturbed boundary between materials of different densities. An example would be a shock traveling through a tube containing two different gases, with a small ripple in the boundary separating them. The shock compresses the gases and sets them in motion. Rather than simply carrying the ripple along unchanged, the interaction can cause it to grow, fold, and eventually roll up into more complicated structures.

If you have tried the Rayleigh-Taylor app, some of these shapes will probably look familiar. The main difference is how the motion gets started. In the RT example, acceleration continues to drive the instability as the heavy material falls and the light material rises. In RM, the shock provides a brief, strong acceleration. The motion created during that interaction continues after the shock has passed. You can think of it as giving the interface a sudden push instead of continuously pulling on it.

The wavy boundary is important because the pressure change across the shock and the density change across the interface do not point in the same direction everywhere. This misalignment generates rotation in the fluid, which helps the interface develop its characteristic rolling shapes. The technical name for this is baroclinic vorticity generation. As the materials stretch and fold around each other, the flow can become increasingly complicated and may eventually develop turbulence. Whether the shock travels from heavy to light material or light to heavy material also affects the direction of the generated rotation. In the heavy to light case, the ripple can first flatten and then grow in the opposite direction, with its original peaks becoming troughs and its troughs becoming peaks. This reversal of the ripple is called phase inversion. (Review of the Richtmyer-Meshkov instability)

This app is meant to make that sequence easier to see. You can watch the shock approach, cross the interface, and leave a growing disturbance behind. Changing the initial conditions lets you explore which parts of that development depend on the shock, which depend on the materials, and which depend on how well the simulation resolves the flow.

I would start by running the default case and watching what happens before changing anything. You can use "Advance to impact" to calculate forward to the estimated time when the shock reaches the average interface position. Since the interface is wavy, some parts encounter the shock earlier than others. The time labeled t* is measured relative to that estimated arrival, so negative values mean before impact and positive values mean after it.

During the early development of a small disturbance, the canonical Richtmyer-Meshkov instability model predicts that its amplitude grows approximately linearly with time. This is different from the exponential growth associated with the ideal RT growing mode. A common way to express the initial growth rate is:

η̇ ≈ k A ΔU η

Here, η is the ripple amplitude just after compression, and η̇ is its initial rate of growth. The wavenumber is k = 2π / λ, where λ is the wavelength. A, the Atwood number, describes the density contrast after the shock, with its sign determined by the material ordering, and ΔU is the change in interface velocity. These are not all the same as the initial values selected in the controls. The shock changes the conditions as it passes through. The Atwood slider specifies the initial contrast as a nonnegative value, while Shock direction selects the material ordering. This equation is useful for understanding the trends, but the app calculates the fluid motion rather than using this expression to animate the interface. (Richtmyer’s original analysis)

The main settings let you explore different parts of this interaction:

• Shock strength (Mach) sets the shock speed relative to the sound speed in the material ahead of it. A higher Mach number produces stronger compression and generally gives the interface a larger velocity change. (Normal shock equations)
• Atwood number sets the initial density difference relative to the sum of the densities. Zero means the materials have equal densities, which removes the density contrast driving this instability. The shock still travels through and compresses the fluid. The default value of 0.5 makes the heavy material three times as dense as the light material.
• Shock direction chooses whether the shock travels from light material into heavy material or the other way around. The incident shock always travels from left to right in this app. Reversing the material order changes the interaction and can cause the ripple to initially flatten and reverse before growing.
• Waves sets how many ripples fit across the domain. More waves means shorter wavelengths. Because amplitude is specified relative to wavelength, adding waves also reduces their absolute height.
• Amplitude / wavelength sets how large each ripple is compared with its wavelength. Larger values start with a more distorted interface and can reach the folding stage sooner.
• Resolution controls how much detail the grid can represent. A coarse grid may smooth away smaller structures. Try repeating the same case at different resolutions to see which features remain consistent.

Material fraction is the clearest display for watching the interface deform. Density and Pressure help show what the shock does to the surrounding fluid. "Show shocks" highlights pressure gradients, so it can also reveal pressure features that are not shocks. Intermediate material colors represent the relative material masses in a cell.

The general motivation and setup for these examples are described on the Instability page. This app uses a different physics solver based on the compressible Euler equations, allowing it to represent shock waves and compression. Both materials use an ideal gas equation of state with a fixed ratio of specific heats of 1.4. (Euler equations and ideal gas flow)

How the calculations work

The domain is divided into a fixed grid of cells. Each cell stores the mass of each material, momentum, and total energy, all per unit volume. Pressure is calculated from the internal energy remaining after subtracting the kinetic energy. The solver advances these quantities by calculating their exchanges across cell boundaries. Neighboring cells use the same exchange at their shared boundary, so what leaves one cell enters the next. This is the basis of the conservative finite volume method.

The materials begin on opposite sides of a wavy interface with a smooth transition defined by a hyperbolic tangent function. An already formed planar shock starts to the left of the interface. The conditions behind the shock are calculated from the selected Mach number using the normal shock jump relations, which connect density, pressure, and velocity through conservation of mass, momentum, and energy. (Normal shock equations)

Calculating the exchanges between cells

The solver estimates conditions at cell boundaries using a linear reconstruction of the values stored in neighboring cells. As in the other instability apps, a monotonized central slope limiter reduces artificial oscillations near sharp changes. This gives second order spatial reconstruction in smooth regions, with more limited reconstruction near discontinuities. (An illustrated explanation of slope limiting)

The estimated states on either side of a boundary form a local Riemann problem. The HLLC solver approximates the resulting wave structure using two outer signal speeds and a middle contact wave. Including the contact wave helps preserve changes in density and material composition. If its intermediate states fail the validity checks, the solver uses the more dissipative HLLE approximation at that boundary. Both methods introduce numerical smoothing, even though physical viscosity is not included. (HLLC solver and wave structure)

Moving forward in time

The solution advances using the same two stage SSP RK2 method described on the Instability page. The first stage produces a trial update. The second advances that trial state again and averages the result with the original state. Here, the material masses, momentum, and energy all advance through these stages together. (SSP RK2 time integration, also described as Heun’s method)

The timestep depends on fluid velocity, sound speed, and cell size. The CFL number defaults to 0.35 and can be adjusted from 0.05 to 0.5 in Advanced options. Changes apply on reset. Updates are checked for finite values, positive density and pressure, and acceptable material fractions. Failed steps are retried using first order reconstruction with HLLE fluxes, followed by smaller timesteps if needed. If recovery fails, the last accepted state remains available for inspection. (CFL condition)

Boundaries and reference frame

The top and bottom boundaries are periodic. At the left and right boundaries, the solver copies the adjacent interior state. This permits flow through the boundaries but does not eliminate all wave reflections. The wider domain provides additional space for checking their influence. (Boundary conditions)

The moving frame uses an estimated contact velocity from a separate planar shock interaction calculation. Initial velocities and total energies are transformed consistently into that frame, keeping the interface approximately in place after impact. This is a fixed reference speed rather than continuous tracking of the distorted interface.

How this runs in the browser

WebGL 2 hardware acceleration is implemented similarly to the other instability apps, allowing grid calculations to run in parallel on the GPU. The NVIDIA fluid simulation explanation linked on the Instability page provides background on this approach. (GPU fluid simulation)