NUMERICAL ANALYSIS OF RAYLEIGH-TAYLOR INSTABILITY

Performing 2D Multiphase Simulation of Rayleigh Taylor Instability for Air and Water as Its Two Phases.

AIM:

• Observing the differences in the effect of mesh size on simulation.
• Performing simulation with water and user-defined material (density = 400 kg/m3, viscosity = 0.001 kg/m-s) for refined mesh.

OBJECTIVE:

• Determining Atwood number for both the cases and explaining how the variation in Atwood number in the above two cases affects the behavior of the instability.

CFD Models Based on Mathematical Analysis of Rayleigh Taylor Waves:

Rayleigh–Taylor instability Model:

The Rayleigh–Taylor instability, or RT instability (after Lord Rayleigh and G. I. Taylor), is an instability of an interface between two fluids of different densities which occurs when the lighter fluid is pushing the heavier fluid. Examples include the behavior of water suspended above oil in the gravity of Earth, mushroom clouds like those from volcanic eruptions and atmospheric nuclear explosions, supernova explosions in which expanding core gas is accelerated into denser shell gas, instabilities in plasma fusion reactors and inertial confinement fusion.

RT instability fingers evident in the Crab Nebula:

Atwood number is an important parameter in the study of Rayleigh–Taylor instability and Richtmyer–Meshkov instability. Atwood number (A) is a dimensionless number in fluid dynamics used in the study of hydrodynamic instabilities in density stratified flows. It is a dimensionless density ratio defined as:

A = (ρ₁ – ρ₂) / (ρ₁ + ρ₂)

Where,

ρ₁ = density of heavier fluid

ρ₂ = density of lighter fluid

Penetration distance of heavy fluid bubbles into the light fluid is a function of acceleration time scale, given by: A×g×t²

Where, g is the gravitational acceleration and t is the time.

Richtmyer–Meshkov instability Model:

The Richtmyer–Meshkov instability (RMI) occurs when two fluids of different density are impulsively accelerated. Normally this is by the passage of a shock wave. The development of the instability begins with small amplitude perturbations which initially grow linearly with time. This is followed by a nonlinear regime with bubbles appearing in the case of a light fluid penetrating a heavy fluid, and with spikes appearing in the case of a heavy fluid penetrating a light fluid. A chaotic regime eventually is reached and the two fluids mix. This instability can be considered the impulsive-acceleration limit of the Rayleigh–Taylor instability,

Kelvin–Helmholtz instability Model:

The Kelvin–Helmholtz instability (after Lord Kelvin and Hermann von Helmholtz) typically occurs when there is velocity shear in a single continuous fluid, or additionally where there is a velocity difference across the interface between two fluids. A common example is seen with wind blowing over water, the instability constant is able to manifest itself through waves on a water surface. The Kelvin-Helmholtz instability is not only restricted to a water surface as clouds, but is evident through other natural phenomena as the ocean, Saturn's bands, Jupiter's Red Spot, and the sun's corona.

The dynamics of the Kelvin-Helmholtz instability is described by the Taylor–Goldstein equation and its onset is given by the Richardson number Ri.

A KH instability on the planet Saturn, formed at the interaction of two bands of the planet's atmosphere:

Kelvin–Helmholtz - Rayleigh–Taylor Hybrid Model:

An improved spray atomization model is presented for use in both diesel and gasoline spray computations. The KH-RT hybrid atomization model consists of two distinct steps: primary and secondary breakup. The Kelvin-Helmholtz (KH) instability model was used to predict the primary breakup of the intact liquid core of a diesel jet. The secondary breakup of the individual drops was modeled with the Kelvin-Helmholtz model in conjunction with the Rayleigh-Taylor (RT) accelerative instability model.

Plateau–Rayleigh instability Model:

Plateau–Rayleigh instability (or Rayleigh instability) of a liquid jet. This instability, sometimes called the hosepipe (or firehose) instability, occurs due to surface tension, which acts to break a cylindrical jet into a stream of droplets having the same total volume but higher surface area. This 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 concerned with fluid thread breakup.

The driving force of the Plateau–Rayleigh instability is that liquids, by virtue of their surface tensions, tend to minimize their surface area.

NUMERICAL ANALYSIS:

FLOW DOMAIN:

The bottom 2D surface is for lighter fluid and the top 2D surface s for the denser fluid. The fluid domain is set to share at the interface between the two 2D surfaces.

TRANSIENT SIMULATION 1:

Materials used: Air (for the bottom surface), Water (for the top surface) (properties below)

The simulation is conducted with only gravity (g = 9.81 m/s²) as the driving force.

Hence, Laminar Model. Since, for a flow domain length of 40 mm the Reynolds number will not exceed the threshold of turbulence.

Solution Method: SIMPLE

M < 0.3 (Here, we do not need to take the compressibility factor into account; Incompressible flow)

Hence, Pressure based solver is used.

Multiphase Model: Homogenous Model - Volume of Fluid

Number of Eulerian Phases – 2

Volume Fraction Parameters – Formulation – Implicit

Multiphase Models `->` Phases

Primary Phase – Air

Secondary Phase – Water (Liquid)

Time Step Size – 0.005 s

Number of Time Steps - 1500

BASELINE MESH:

Patch control is used to assign materials to their respective domains after standard initialization.

SIMULATION ON BASELINE MESH:

REFINED MESH:

Patch control is used to assign materials to their respective domains after standard initialization.

Solution Method `->` Spatial Discretization `->` Momentum = First Order Upwind

Solution Controls `->` Under-Relaxation Factors: Pressure = 0.2

Momentum = 0.6

SIMULATION ON REFINED MESH:

INFERENCE:

CASES	ELEMENT SIZE	ELEMENT COUNT
BASELINE MESH	0.5 mm	3200
REFINED MESH	0.1 mm	80000

When we increased the mesh size more nodes were made available for interpolation. Hence, the simulation became more accurate to the physics involved with the Rayleigh-Taylor Instability.

Here, Steady-state is not preferred as we need to capture the collapse of the interface between the distinct phases (i.e, Top Phase- Denser Water liquid and Bottom Phase- Lighter Air) as a function of time using volume of fluid method to study about the Rayleigh-Taylor Instability.

TRANSIENT SIMULATION 2:

Materials used: User Defined (for the bottom surface), Water (for the top surface) (properties below)

The simulation is conducted with only gravity (g = 9.81 m/s²) as the driving force.

Hence, Laminar Model. Since, for a flow domain length of 40 mm the Reynolds number will not exceed the threshold of turbulence.

Solution Method: SIMPLE

M < 0.3 (Here, we do not need to take the compressibility factor into account; Incompressible flow)

Hence, Pressure based solver is used.

Multiphase Model: Homogenous Model - Volume of Fluid

Number of Eulerian Phases – 2

Volume Fraction Parameters – Formulation – Implicit

Multiphase Models `->` Phases

Primary Phase – User Defined

Secondary Phase – Water (Liquid)

Time Step Size – 0.005 s

Number of Time Steps - 1500

MESH:

Patch control is used to assign materials to their respective domains after standard initialization.

Solution Method `->` Spatial Discretization `->` Momentum = First Order Upwind

Solution Controls `->` Under-Relaxation Factors: Pressure = 0.2

Momentum = 0.6

SIMULATION WITH DENSER LIQUID WATER AT TOP AND LIGHTER USER DEFNED MATERIAL AT BOTTOM:

INFERENCE:

Atwood Number, A = (ρ₁ – ρ₂) / (ρ₁ + ρ₂)

Where,

ρ₁ = density of heavier fluid

ρ₂ = density of lighter fluid

For Transient Simulation 1, for air and water:

ρ₁ = 998.2 Kg/m³

ρ₂ = 1.225 Kg/m³

A = `(998.2 -1.225) / (998.2 + 1.225)` = 0.997548904

For Transient Simulation 2, for user defined material and water:

ρ₁ = 998.2 Kg/m³

ρ₂ = 400 Kg/m³

A = `(998.2 - 400) / (998.2 + 400)` = 0.4278357889

CONCLUSION:

The Atwood number for the first transient simulation is higher than the second transient simulation. After the application of gravity into the simulation, it took longer for the second simulation to start the instability. It also took a longer time for the second simulation to attain equilibrium. This may be due the lower difference in densities of the two phases in the second simulation compared to the first simulation. Hence, this disparity in densities showcased in the Atwood number indicates how long an instability lasts before reaching equilibrium.