Rayleigh Taylor Instability Challenge
OBJECTIVE
- Explain different CFD models used for solving Rayleigh Taylor Instability problems
- Perform a Rayleigh Taylor instability analysis using the Volume of Fluid (VOF) multiphase model for air and water.
- Repeat the process for water and user-defined material with density = 400 kg/m3, viscosity = 0.001 kg/m-s, and compare the behavior
- Define the Atwood number and derive it for each case. Explain the effect of Atwood number on instability.
THEORY
The Rayleigh–Taylor instability or RT instability is the instability of an interface between two fluids of different densities which occurs when the lighter fluid is pushing the heavier fluid. As the RT instability develops, the initial perturbations progress from a linear growth phase into a non-linear growth phase, eventually developing "plumes" flowing upwards (in the gravitational buoyancy sense) and "spikes" falling downwards. In the linear phase, the fluid movement can be closely approximated by linearly, and the amplitude of perturbations is growing exponentially with time. In the non-linear phase, perturbation amplitude is too large for a linear approximation, and non linearity is required to describe fluid motions.
Water suspended atop oil is an everyday example of Rayleigh–Taylor instability, and it may be modeled by two completely plane-parallel layers of immiscible fluid, the denser fluid on top of the less dense one and both subject to the Earth's gravity. The equilibrium is unstable to any disturbances of the interface. If a parcel of heavier fluid is displaced downward with an equal volume of lighter fluid displaced upwards, the potential energy of the configuration is lower than the initial state. Thus the disturbance will grow and lead to a further release of potential energy, as the denser material moves down under the effect of the gravitational field, and the less dense material is displaced upwards. The image shows a typical RT instability case.
Different examples of RT instability in real life are a Mushroom cloud in volcanic eruptions, atmospheric nuclear explosion, supernova explosion in outer space, instabilities in plasma fusion and inertial confinement fusion.
Different types of Instabilities
1) Richtmyer–Meshkov instability
The Richtmyer–Meshkov instability (RMI) occurs when two fluids of different density are impulsively accelerated. Normally this is by the passage of a shove wave. The development of the instability begins with small amplitude perturbations which initially grow linearly with time. This instability can be considered the impulsive-acceleration limit of the Rayleigh-Taylor Instability
2) Kelvin Helmholtz indtability
The Kelvin–Helmholtz instability can occur when there is velocity shear in a single continuous fluid or where there is a velocity difference across the interface between two fluids. An example is wind blowing over water.
3) Plateau Rayleigh instability
The Plateau–Rayleigh instability, often just called the Rayleigh instability, explains why and how a falling stream of fluid breaks up into smaller packets with the same volume but less surface area. It is related to the Rayleigh-Taylor Instability and is part of a greater branch of fluid dynamics concerned with fluid thread breakup. This fluid instability is exploited in the design of a particular type of ink jet technology whereby a jet of liquid is perturbed into a steady stream of droplets. The driving force of the Plateau–Rayleigh instability is that liquids, by virtue of their surface tension, tend to minimize their surface area.
4) Staffman Taylor instability
The Saffman–Taylor instability, also known as viscous fingering, is the formation of patterns in a morphologically unstable interface between two fluids in a porous medium. This situation is most often encountered during drainage processes through media such as soils. It occurs when a less viscous fluid is injected, displacing a more viscous fluid.
CFD model based on mathematical analysis of RT instability
1) Volume of Fluid (VOF)
In CFD, the volume of fluid (VOF) method is a free surface modeling technique for tracking and locating the free surface or fluid - fluid surface. It belongs to the class of Eulerian methods which are characterized by a mesh that is either stationary or is moving in a certain prescribed manner to accommodate the evolving shape of the interface
2) K-L turbulence model
The K-L model was proposed by Dimonte and Tipton for describing the turbulent self-similar regime of RT and RM induced mixing. The model describes the dominant eddies in the mixing zone with evolutionary equations for their characteristic dimension (L) and energy per unit mass (K).
3) Kelvin-Helmholtz / Rayleigh Taylor Hybrid Model
This model was used to predict the primary break up of the intact liquid core of the jet and the Raleigh Taylor model determines the effects of acceleration of the different phases. This model predicts the shape of the spray considering the temperature at different phaces.
PROCEDURE
Following procedure is followed in setting up a simulation
Geometry - 2 components created using surface pull option each of 20mm x 20mm sharing a common edge are created

Mesh - Basic quad mesh is created. Baseline mesh size = 0.5mm. Refined = 0.25 mm and 0.1 mm

| Mesh - Configuration | No of Nodes | No of Elements |
| Mesh - baseline - 0.5 mm | 3321 | 3200 |
| Mesh - 0.25 mm | 13041 | 12800 |
| Mesh - 0.1 mm | 80601 | 80000 |
Setting up Physics - Transient, Pressure based, Multi-phase -> VOF -> 2 Eulerian fluid -> Implicit, Create 2 phase -> Air (Primary phase) & Water (secondary phase)
Initialization - Standard initialization, Patch -> water fraction -> water mixture = 1 -> air mixture = 1
Output - Contour -> Phase -> water -> select all -> save/display (as shown below), Create solution animation -> create every 2 timesteps.
Run Calculation - different timesteps size like 0.005 and 0.0025 were tried to create a successful simulation and avoid any floating point.
RESULTS
1) Case I Baseline mesh, mesh size = 0.5mm, timestep size = 0.005, Air + Water
Residuals Plot

Simulation
2) Case II Refined mesh 1, mesh size = 0.25 mm, timestep size = 0.005, Air + Water
Residuals Plot

Simulation
3) Case III Refined mesh 2, mesh size = 0.1 mm, timestep size = 0.0025, Air + Water
Residuals Plot

Simulation
4) Case IV Refined mesh 2, mesh size = 0.1 mm, timestep size = 0.0025, user defined material(density = 400 kg/m3, viscosity = 0.001 kg/m-s) + Water
Residuals Plot

Simulation
OBSERVATIONS
We can observe following things from each simulation
Stage 1 - smaller perturbations are generated resembling a sinusoidal waves at the interface
Stage 2 - bubbles and spikes are formed based on the domain density difference. At interface lower density bubbles pushes into the higher density fluid to form bubbles and spikes are formed when high density fluid pushes low density fluid
Stage 3 - spikes and bubbles starts interacting with one another. As a result large mushroom like structures are formed before complete instability
Stage 4 - mixing of fluids happen and the flow becomes turbulent due to high reynolds number.
Equilibrium - After sometimes the flow settles down and attains the state of equilibrium with high density fluid at bottom and low density fluid at top and variable fluid mixture in the middle and no clear interface between them.





Why Steady State simulation is not suitable ?
- We know that steady state simulation is a time independent approach.
- Based on the simulation we know there is high turbulence involve and the change in the flow occurs between very small time variation. This is not possible to capture in steady state solution
- Also steady state solver ignore the time differential terms involved in the simulation and thus even if it may be able to capture the state of equilibrium, it skips to capture the non linearities in the simulation precisely.
- Below video represents a case when steady state solver is used to simulate RT instability
ATWOOD Number
It is a dimensionless number in fluid dynamics used to study hydrodynamic instabilities in density stratified flows
A = `(ρ1-ρ2)/(ρ1+ρ2)`
Case A - Air + Water, ρ1 = 998 and ρ2 = 1.225=> A = 0.99754
Case B - User defined material + Water, ρ1 = 998 and ρ = 400 => A = 0.42775
Effect of Atwood Number
- When A is near 1, there is clear formation of bubble of low density fluid and spikes of high density fluid.
- When A is away from 1, it is difficult to differentiate between the bubbles and spike due to less density difference.
- Also the simulation attains equilibrium quickly when A ~ 1 where as the equilibrium time increases considerably when A is distant from 1 as the fluid mixes too well.