Rayleigh Taylor instability simulation

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 against any external force(gravity here).

The lighter fluid exerts a force on the denser fluid in the presence of gravity.

As the perturbation grows, the potential energy of the total system reduces as the heavier fluid displaces the lighter fluid at the bottom, because any system practically strives to attain a stable position by losing its energy.

Some practical observations of Rayleigh Taylor instability:

Examples include the behaviour 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.

Models used for Rayleigh-Taylor instability

1. Early-time growth
2. Potential flow
3. Buoyancy-drag balance
4. Energy budget

1. Early-time growth

It was noted in Taylor (1950) that two superposed fluids with an interface normal to the acceleration field could exist in a state of unstable equilibrium. However if small corrugations were to exist on the interface, a growing standing wave pattern could be expected. By considering an unbounded two-dimensional domain, incompressible potential flow, and a sinusoidal interface, the velocity potentials φ in upper and lower layers are given by

1. Potential flow

Development of unstable interfaces beyond the linear phase was first studied by Davies & Taylor (1950), and generalised somewhat in Layzer (1955). Both papers begin by considering a two-dimensional incompressible potential flow in a vertical cylinder with a initially flat free under-surface. This is Rayleigh-Taylor unstable with an Atwood number A → 1. The velocity potential must satisfy Bernoulli’s equation in cylindrical polar ordinates: where α (t) is arbitrary

1. Buoyancy-drag balance

The complexity of directly evaluating a potential flow to calculate growth behavior can be avoided. From first principles, the following equation describes how buoyancy and drag compete to accelerate a body, which in this case is a bubble of lower layer fluid penetrating the upper layer:

4. Energy-budget method

It is an application of the law of conservation of energy. The energy available for evaporation is determined by considering the incoming energy, outgoing energy and energy stored in the water body over a known time interval.

5.Lattice Boltzmann Scheme

The lattice Boltzmann scheme is based on distribution functions. The dynamic functions at the interface such as phase segregation and surface tension. These functions are derived from the continuous Boltzmann functions with many approximations.

It is only suitable for incompressible flows. An index function is used to track the interference between the phases. However, this model is only suitable for 2D simulation.

6. K L Turbulence model: This method is used to simulate the self-similar growth of the RT instability and Richtmyer-Meshkov instabilities. This CFD model captures the dominant eddies with the parameters K and L.

Where,

K = V²/2 energy per unit mass

L = Characteristic length.

7. Kelvin-Helmholtz and Rayleigh Taylor Hybrid model: the KH RT hybrid model is a CFD model that efficiently simulates the atomization of fuel spray. This model consists of two distinct steps, The KH instability model determines the primary break -up of the jet. The RT instability determines the effect of acceleration of the different phases.

This model is captures the penetration of fluid in different phases where temperature is taken into consideration. This model predicts the spray shape. Additionally, Rosin-Rammler distribution is used to determine the size of children drop after the disintegration of the parent drop.

Mathematical model:

A 2D model of size 20mm wide and 40mm high is used to visualize the phenomenon. Standard initialisation is done and patch is done so as to distinguish the two phases. Since two fluids are involved multiphase is used with a transient solver. In the multiphase model, a standard implicit volume cut-off fraction value of 1 e-6 is set as default.

Mesh Independency Study

The simulation is carried out for three different grid size. The base mesh is of 0.5 mm, it is further refined to 0.3 mm and 0.2 mm respectively as shown in table 1.

Table 1. Simulation for different grid size

Size	Nodes	Elements
0.5 mm	3321	3200
0.3mm	9180	8978
0.2mm	20301	20000

Case 1:

In this case a Rayleigh instability is studied for element size of 0.5 mm. Multiphase model is chosen and this case is made to run in Transient step simulation.

Explicit solver is used with variable time stepping method,with the minimum time step being equal to 0.05 and the courant’s number=0.25.

No. of time steps for which the calculation is run= 1000

Fig 1. The volume fraction contour after 1000 time steps

v Fig 2. The scaled residuals

Case 2:

The simulation is made to run with a finer mesh in this case. In this case the mesh size is 0.3 mm, and the other solving parameters were unaltered.

The contours obtained for Case 2 are as follows.

Fig 3. The volume fraction contour for case 2.

Fig 4. The residuals plot after 1000 time steps

Case 3:

In this 3^rd case, a finer mesh is applied to the model. The size of the mesh chosen was 0.2mm.

Fig 6. The volume fraction contour for Case 3 simulation.

Fig 7. The residual plot for case 3.

Significance of Atwood’s number:

The Atwood’s number is a dimensionless number used in fluid dynamics to study hydrodynamic instabilities in density stratified flows.

Its denoted by the letter ‘A”.

A = (ρ₁ - ρ₂) / ( ρ_{1 ­}+ ρ₂)

Where,

ρ₁ = the density of the heavier fluid

ρ₂ = density of the lighter fluid.

Atwood number serves as an important parameter in the study of Rayleigh Taylor instability and Richtmyer-Meshkov instability.

Particularly, in Rayleigh-Taylor instability, the penetration distance is seen as a function of acceleration time scale Agt² where g is the gravitational constant and t is the time.

Calculation of Atwood number:

In our case, the heavier fluid is water and the lighter fluid is air.

ρ₁ = 997 Kg/m³

ρ₂ = 1.225 Kg/m³

A = (997 – 1.225) / (997 + 1.225) = 0.99754

Implication of Atwood number:

For instabilities at RT configuration for large density differences, the molecular mixing develops in the shape of large spikes on the lighter fluid zone and bubbles on the heavy fluid zone as seen from the contour plot. For smaller values of Atwood no. in the range of less than 0.5, the nuances in the interface won’t be appreciably noticeable.

In this case since the density difference between the air and water is very large, the formation of spikes and bubbles.

Why a steady state solver is not suitable?

A steady state solver is used only when the final state (equilibrium) is taken into account. The time step variation of the parameters involved can’t be determined accurately using a steady state solver. On the other hand, a transient state simulation allows us to use a custom time step, either linear or variable with predefined volume cut-off ratio.

By using a transient solver, the time step can be varied so as to obtain the change of state from the initial state to the steady state. While, a steady state solver only focuses on the final result that is obtained at the end of the simulation.

In a transient state solver both constant time step as well as variable time step function can be used.

YOUTUBE Links for the animations of the 3 cases

https://youtu.be/Cz4r3ehwZ10 mesh size= 0.5mm

https://youtu.be/_KjIq6vaoL4 mesh size= 0.3mm

https://youtu.be/5sP0Rg_YcJQ , mesh size= 0.2mm

uploader name= Karan Gupta([email protected])