Experimental Physics Final Project (Group)

Modeling and Analysis of Kelvin–Helmholtz (KHI) and Rayleigh–Taylor (RTI) Instabilities

AdvisorNTHU Physics

PeriodSpring 2025

A numerical study of two classical hydrodynamic instabilities common in geophysical and astrophysical settings: the Kelvin–Helmholtz instability (KHI), which arises at fluid interfaces with velocity shear, and the Rayleigh–Taylor instability (RTI), which occurs when a denser fluid is supported above a lighter fluid in a gravitational field. Both are responsible for a wide range of phenomena from ocean waves and cloud formations to supernova shock fronts and inertial confinement fusion. This project implements a 2D compressible Euler solver from scratch, using the HLLC Riemann solver for shock-capturing and a dynamic CFL-based time-step scheme. Simulations are run across a parameter space of density ratios and velocity shear amplitudes to study how these factors control instability growth, vortex roll-up, and energy redistribution. The solver code (RTRT2.py) is also accelerated with Numba JIT compilation for high performance on large grids.

▪Implemented a full 2D compressible Euler solver: spatial HLLC fluxes (x and y directions) with Numba @jit(nopython=True) acceleration on Nx×Ny grids up to 512×512.
▪Boundary conditions: periodic in x (for KHI); reflective extrapolation on top/bottom walls (for RTI) to model a closed gravitational box.
▪Initial conditions for RTI: hyperbolic-tangent smoothed density interface at y=0.5L with cosine perturbation amplitude 0.05L and a gravitational source term g=2 applied each timestep via operator splitting.
▪Initial conditions for KHI: two fluid layers with horizontal velocity shear (v_upper − v_lower) and a thin tanh-smoothed interface to avoid numerical artifacts.
▪Parameter studies: density ratios ρ_upper/ρ_lower ∈ {0.3, 2.7, 6} for KHI; {3, 5, 10} for RTI. Velocity shear amplitudes |Δv| ∈ {0.2, 0.5, 0.8} for KHI.
▪Visualization: density field and kinetic energy density (½ρv²) snapshots exported at regular intervals, assembled into .mp4 simulation videos.

Methodology

1. Compressible Euler Equations

The simulation solves the 2D compressible Euler equations in conservation form:

\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = \mathbf{S}

where the conserved state vector and fluxes are:

\mathbf{U} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ E \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho u \\ \rho u^2 + P \\ \rho u v \\ (E+P)u \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \rho v \\ \rho u v \\ \rho v^2 + P \\ (E+P)v \end{pmatrix}

with the ideal gas equation of state: $P = (\gamma - 1)\left(E - \frac{1}{2}\rho(u^2+v^2)\right)$ , $\gamma = 5/3$ .

2. HLLC Riemann Solver

Intercell fluxes are computed via the HLLC (Harten–Lax–van Leer–Contact) approximate Riemann solver, which captures both shocks and contact discontinuities. For each cell interface, three wave speeds are estimated:

S_L = \min(u_L - c_L,\; u_R - c_R), \quad S_R = \max(u_L + c_L,\; u_R + c_R)

S^* = \frac{P_R - P_L + \rho_L u_L (S_L - u_L) - \rho_R u_R (S_R - u_R)}{\rho_L(S_L-u_L) - \rho_R(S_R-u_R)}

The HLLC flux selects between the left, star-left, star-right, and right states depending on the sign of $S_L$ , $S^*$ , and $S_R$ . This correctly resolves contact waves that the simpler HLL solver smears out, preserving sharp density interfaces critical for resolving instability morphology.

3. Gravitational Source Term (RTI)

For RTI simulations, a uniform gravitational acceleration $g$ in the $-y$ direction is applied via operator splitting after each flux update:

\rho v \leftarrow \rho v - g \rho \Delta t

E \leftarrow E - g \rho v \Delta t

This maintains conservative form and correctly models the buoyancy-driven instability driving force.

4. CFL-Based Adaptive Time-Stepping

To ensure numerical stability, the time step is dynamically bounded by the CFL condition:

\Delta t = C_{\text{CFL}} \cdot \min_{i,j} \left( \frac{\Delta x}{|u_{ij}| + c_{ij}},\; \frac{\Delta y}{|v_{ij}| + c_{ij}} \right)

where $c = \sqrt{\gamma P/\rho}$ is the local sound speed and $C_{\text{CFL}} = 0.1$ for stability.

5. Initial Conditions with Smoothed Interface

To avoid spurious numerical noise at sharp density jumps, the interface density profile is initialized with a hyperbolic tangent transition of width $\delta$ :

\rho(y) = \rho_L + (\rho_U - \rho_L) \cdot \frac{1 + \tanh\left(\frac{y - y_0(x)}{\delta}\right)}{2}

The interface itself carries a small cosine perturbation:

y_0(x) = 0.5L + A \cos\left(\frac{2\pi x}{\lambda}\right)

with amplitude $A = 0.05L$ and wavelength $\lambda = L$ , seeding the instability growth.

Results & Analysis

RTI: Density Finger and Mushroom Cap Evolution

The RTI density field evolves in three distinct phases. Initially, the heavy fluid (red) sits above the light fluid (blue) with only a gentle cosine perturbation at the interface. As time progresses, the gravitational buoyancy force amplifies the perturbation into a downward-plunging spike ("finger") of heavy fluid. By late times, the head of the finger broadens into the characteristic mushroom-cap shape, while the shoulders develop secondary KHI roll-ups driven by the velocity shear at the finger edges.

Effect of Density Ratio on RTI Growth

Comparing runs at density ratios ρ_upper/ρ_lower ∈ {3, 5, 10}, higher density contrast produces a narrower, faster-falling spike with sharper KHI curls at the edges. At ρ=10:1, the spike tip reaches the bottom boundary earlier and the mushroom cap is significantly thinner, consistent with the Atwood number scaling:

A_t = \frac{\rho_U - \rho_L}{\rho_U + \rho_L}, \quad \gamma_{\text{RTI}} = \sqrt{A_t g k}

where $k$ is the wavenumber of the perturbation.

Kinetic Energy Density Evolution

The kinetic energy density maps (½ρv²) capture how energy is spatially redistributed as the instability evolves. At early times, kinetic energy is concentrated as a thin bright band along the interface. As the spike descends, it forms a concentrated jet of high kinetic energy. By late time, energy spreads along the full length of the spike and into the KHI vortices at the shoulders.

KHI: Velocity Shear and Roll-Up

For KHI, smaller density contrast and larger velocity shear both accelerate roll-up formation. The vortex sheet at the interface wraps into a coherent spiral structure. At higher shear velocities (Δv=0.8), roll-up is complete by earlier times, and secondary instabilities are visible along the vortex arms.

Figures

Figure 1: RTI initial state (ρ=3:1) — cosine-perturbed interface between heavy (red) and light (blue) fluid, just before instability growth

Figure 2: RTI mid-stage (ρ=3:1, t=130) — heavy fluid finger plunging downward, mushroom cap beginning to form

Figure 3: RTI late stage (ρ=3:1, t=250) — fully developed mushroom cap with secondary KHI roll-ups at the shoulders

Figure 4: RTI at ρ=5:1 (t=169) — higher density contrast produces a narrower spike with tighter KHI spiral vortices at the edges

Figure 5: Kinetic energy density map (t=100) — energy concentrated along the descending spike interface as the RTI finger accelerates

Figure 6: Kinetic energy density map (t=250) — energy redistributed along the full spike length and into the KHI vortex arms

Simulation Videos

Video 1: Rayleigh–Taylor instability — density field evolution showing finger formation, mushroom cap development, and secondary KHI roll-ups

Video 2: Kelvin–Helmholtz instability (ρ=2.7:1) — velocity shear at the interface drives vortex roll-up and spiral entrainment of the two fluid layers

Video 3: KHI parameter study (shear velocity Δv=0.5) — intermediate shear showing early-stage roll-up dynamics and interface deformation

Back to ResearchExperimental Physics Final Project (Group)

Modeling and Analysis of Kelvin–Helmholtz (KHI) and Rayleigh–Taylor (RTI) Instabilities

Advisor

NTHU Physics

Research Period

Spring 2025

Overview

Implemented a full 2D compressible Euler solver: spatial HLLC fluxes (x and y directions) with Numba @jit(nopython=True) acceleration on Nx×Ny grids up to 512×512.
Boundary conditions: periodic in x (for KHI); reflective extrapolation on top/bottom walls (for RTI) to model a closed gravitational box.
Initial conditions for RTI: hyperbolic-tangent smoothed density interface at y=0.5L with cosine perturbation amplitude 0.05L and a gravitational source term g=2 applied each timestep via operator splitting.
Initial conditions for KHI: two fluid layers with horizontal velocity shear (v_upper − v_lower) and a thin tanh-smoothed interface to avoid numerical artifacts.
Parameter studies: density ratios ρ_upper/ρ_lower ∈ {0.3, 2.7, 6} for KHI; {3, 5, 10} for RTI. Velocity shear amplitudes |Δv| ∈ {0.2, 0.5, 0.8} for KHI.
Visualization: density field and kinetic energy density (½ρv²) snapshots exported at regular intervals, assembled into .mp4 simulation videos.

Methodology

1. Compressible Euler Equations

The simulation solves the 2D compressible Euler equations in conservation form:

\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = \mathbf{S}

where the conserved state vector and fluxes are:

\mathbf{U} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ E \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho u \\ \rho u^2 + P \\ \rho u v \\ (E+P)u \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \rho v \\ \rho u v \\ \rho v^2 + P \\ (E+P)v \end{pmatrix}

with the ideal gas equation of state: $P = (\gamma - 1)\left(E - \frac{1}{2}\rho(u^2+v^2)\right)$ , $\gamma = 5/3$ .

2. HLLC Riemann Solver

S_L = \min(u_L - c_L,\; u_R - c_R), \quad S_R = \max(u_L + c_L,\; u_R + c_R)

S^* = \frac{P_R - P_L + \rho_L u_L (S_L - u_L) - \rho_R u_R (S_R - u_R)}{\rho_L(S_L-u_L) - \rho_R(S_R-u_R)}

3. Gravitational Source Term (RTI)

For RTI simulations, a uniform gravitational acceleration $g$ in the $-y$ direction is applied via operator splitting after each flux update:

\rho v \leftarrow \rho v - g \rho \Delta t

E \leftarrow E - g \rho v \Delta t

This maintains conservative form and correctly models the buoyancy-driven instability driving force.

4. CFL-Based Adaptive Time-Stepping

To ensure numerical stability, the time step is dynamically bounded by the CFL condition:

\Delta t = C_{\text{CFL}} \cdot \min_{i,j} \left( \frac{\Delta x}{|u_{ij}| + c_{ij}},\; \frac{\Delta y}{|v_{ij}| + c_{ij}} \right)

where $c = \sqrt{\gamma P/\rho}$ is the local sound speed and $C_{\text{CFL}} = 0.1$ for stability.

5. Initial Conditions with Smoothed Interface

To avoid spurious numerical noise at sharp density jumps, the interface density profile is initialized with a hyperbolic tangent transition of width $\delta$ :

\rho(y) = \rho_L + (\rho_U - \rho_L) \cdot \frac{1 + \tanh\left(\frac{y - y_0(x)}{\delta}\right)}{2}

The interface itself carries a small cosine perturbation:

y_0(x) = 0.5L + A \cos\left(\frac{2\pi x}{\lambda}\right)

with amplitude $A = 0.05L$ and wavelength $\lambda = L$ , seeding the instability growth.

Results & Analysis

RTI: Density Finger and Mushroom Cap Evolution

Effect of Density Ratio on RTI Growth

A_t = \frac{\rho_U - \rho_L}{\rho_U + \rho_L}, \quad \gamma_{\text{RTI}} = \sqrt{A_t g k}

where $k$ is the wavenumber of the perturbation.

Kinetic Energy Density Evolution

KHI: Velocity Shear and Roll-Up

Visualizations & Trajectories

Figure 1: RTI initial state (ρ=3:1) — cosine-perturbed interface between heavy (red) and light (blue) fluid, just before instability growth

Figure 2: RTI mid-stage (ρ=3:1, t=130) — heavy fluid finger plunging downward, mushroom cap beginning to form

Figure 3: RTI late stage (ρ=3:1, t=250) — fully developed mushroom cap with secondary KHI roll-ups at the shoulders

Figure 4: RTI at ρ=5:1 (t=169) — higher density contrast produces a narrower spike with tighter KHI spiral vortices at the edges

Figure 5: Kinetic energy density map (t=100) — energy concentrated along the descending spike interface as the RTI finger accelerates

Figure 6: Kinetic energy density map (t=250) — energy redistributed along the full spike length and into the KHI vortex arms

Simulation Videos

Video 1: Rayleigh–Taylor instability — density field evolution showing finger formation, mushroom cap development, and secondary KHI roll-ups

Video 2: Kelvin–Helmholtz instability (ρ=2.7:1) — velocity shear at the interface drives vortex roll-up and spiral entrainment of the two fluid layers

Video 3: KHI parameter study (shear velocity Δv=0.5) — intermediate shear showing early-stage roll-up dynamics and interface deformation

Modeling and Analysis of Kelvin–Helmholtz (KHI) and Rayleigh–Taylor (RTI) Instabilities

Methodology

1. Compressible Euler Equations

2. HLLC Riemann Solver

3. Gravitational Source Term (RTI)

4. CFL-Based Adaptive Time-Stepping

5. Initial Conditions with Smoothed Interface

Results & Analysis

RTI: Density Finger and Mushroom Cap Evolution

Effect of Density Ratio on RTI Growth

Kinetic Energy Density Evolution

KHI: Velocity Shear and Roll-Up

Figures

Simulation Videos

Modeling and Analysis of Kelvin–Helmholtz (KHI) and Rayleigh–Taylor (RTI) Instabilities

Overview

Methodology

1. Compressible Euler Equations

2. HLLC Riemann Solver

3. Gravitational Source Term (RTI)

4. CFL-Based Adaptive Time-Stepping

5. Initial Conditions with Smoothed Interface

Results & Analysis

RTI: Density Finger and Mushroom Cap Evolution

Effect of Density Ratio on RTI Growth

Kinetic Energy Density Evolution

KHI: Velocity Shear and Roll-Up

Visualizations & Trajectories

Simulation Videos

[ peer-evaluation log ]