Week 5 - Rayleigh Taylor Instability

Aim - To conduct the Rayleigh Taylor Instability CFD simulation.

Introduction - Rayleigh Taylor Instability (RTI) is a phenomenon that occurs at the interface of fluids of different
density wherein the lighter fluid pushes the heavier fluid due to the density difference. It is a dynamic process
where the system tries to reduce its combined potential energy.

Theory - Initially the system is in a state of hydrostatic equilibrium, in which the heavier fluid (ex: water) sits atop a
lighter fluid (ex: air) in a constant gravitational field. As the RTI develops, the initial perturbations progress from a
linear growth phase, eventually developing "plumes" flowing upwards (in the gravitational buoyancy sense) and
"spikes" falling downwards. In the linear phase, equations can be linearized and the amplitude of perturbations is
growing exponentially with time. In the non-linear phase, the perturbation amplitude is too large for the non-linear
terms to be neglected. In general, the density disparity between the fluids determines the structure of the subsequent
nonlinear RT instability flows (assuming other variables such as surface tension and viscosity are negligible here).
Examples include supernova, water suspended atop oil, mushroom clouds formed due to volcanic eruptions.

CFD models based on mathematical analysis of Rayleigh Taylor waves:
1. Mass diffusion model:
The mass diffusion equation describes the transport of mass and energy due to model viscosity and hence drives
actual mixing in the fluids. Even if any configuration is Rayleigh-Taylor unstable, the diffusion process could be
responsible for significant mixing if the viscosity is sufficiently large.

2. Single fluid model:

2.1. Potential flow model
This model is valid for all instability changes from the linear stage through the early non-linear to the asymptotic but
limited to Atwood number = 1. In this model, single-mode perturbation expands the flow équations to-second-order
around the bubble tips, hence obtaining an ordinary differential equation for the bubble velocity.

2.2. Buoyancy drag model
It is valid for every Atwood number but limited to the non-linear stages of the instability. It uses ordinary differential equations to evolve the width of the mixing layer. The bubble, or spike, amplitudes are described by balancing the inertia, buoyaney and drag forces. These models cannot model multiple mixing interfaces; cannot be easily
extended to 2 or 3 dimensions, and as a rule, do not address demixing, also known as counter gradient transport, ie.
reduction of total fluid masses within the mixing zone. To address the above problems, two-fluid (or multi-
fluid models have to be used.

2.3. Besnard-Harlow-Rauenzahn(BHR)
It is a more advanced version of the single-fluid model. The BHR model is based upon the evolution equations arising
from second-order correlations and gradient- diffusion approximations. Using a mass-weighted averaged
decomposition, the original BHR model includes full transport equations for the Reynold stresses, turbulent mass-flux
velocity, density fluctuations, and the turbulent kinetic energy dissipation rate.
3. Multi-Fluid model:

These models use a separate set of equations for each fluid in addition to the main flow equations and provide an
accurate modeling framework for demixing by correctly capturing the relative motion of the different fluid fragments.
The first step for the development of turbulence models is to perform Reynold averaging of the governing equations
and resultant terms.

3.1. Two equation K-L turbulence model:
The K- turbulence model involves equations for turbulent kinetic energy K and turbulent length scale L. The starting
point for deriving the model equations are the buoyancy-drag models for the self-similar growth of RT instabilities.

3.2. Three-equation K-L-a turbulence model:
The three equation model was developed as an extension to the two-equation -L model by including a third
equation for the turbulent mass-flux velocity.

3.3. Four-equation K-L-a-b turbulence model:
The four equation model was developed as an extension to the three-equation K-L model by including a fourth
equation for density-specific volume covariance. The model can accurately reflect the effects of the changes in the
density fluctuations.

Modeling approach :-

CAD clean up/Pre-process-SpaceClaim - Prepare the surface geometry in spaceclaim. And turn on shareprep.

Meshing :-Apply the elements size about 0.5 mm.

Setup :- In this section, we set up the transient simulation. Apply gravity in negative Y axis. laminar model is applied for simulation. Apply multiphase implicie method. Apply phase name and types accordingly. Prepare phase contour and prepare the animation .Give the patch after standard initialization.Certain No of time steps are set for simulation.

Solution :- we can see the animations of instability of interface between water and air.

CAD clean up/Pre-process-SpaceClaim :- To prepare geometry get two equal square (20 x 20 mm)surface. And turn on shareprep. We get the colored line between these two surfaces.

Meshing :-

Elemet size = 0.5 mm

No of elements = 3200

Mesh quality:-

Mesh metric :-

Here we can see that all the elements have high quality as the range between 0.99 - 1 .

Setup:-

Set the simulation for tansient flow. Consider gravity in negative Y direction as 9.81 m/s^2. Apply laminar model. Take volume to fluid homogeneous model. Apply implicit method as well.

Apply water and air phase accordigly.

Apply patch after standard initialization. Take value 1 for water and for air take 0.

Create phase contour. Prepare animation of phase contour.

To run the simulation , time step size = 0.005 s and No of time steps = 1500

Result:-

Here is the link of the animation how instability occurs in rayleigh taylor instability problem. We clearly see that the bubble formation is unsymmetric but it takes less time to simulate.

https://drive.google.com/file/d/1x0Tw8fU-74lEDtSVFCiES1F7TkeGiPZw/view?usp=sharing

Refined case :-

i)

Meshing :-

Elemet size = 0.35 mm

No of elements = 6498

Mesh quality :-

mesh metric

Here we can see that all the elements have high quality as the range between 0.98 - 1 .

Setup:- It remains same as previous case. Take water and air material.

Result:-

Here is the link of the animation how instability occurs in rayleigh taylor instability problem. We clearly see that the bubble formation is unsymmetric but it takes less time to simulate.

https://drive.google.com/file/d/1QzSfUfDk6xIlR7XiBjo_9lV4nlJagAeh/view?usp=sharing

ii) Water and user define material

Meshing :-

Elemet size = 0.35 mm

No of elements = 6498

Setup:- It remains same as previous case. Take WATER and USER DEFINE material (density = 400 kg/m³, viscosity = 0.001 kg/m-s).

Result:-

Here is the link of the animation how instability occurs in rayleigh taylor instability problem. We clearly see that the bubble formation is symmetric but it takes larger time to simulate.

https://drive.google.com/file/d/1BXM0oWfX5CIvlqf-hy68hcsGBeYRiNPM/view?usp=sharing

Observation:

Finner mesh captures the instability better than the coarse mesh, but it takes more time for the heavier fluid to settle down. So it is computationally expensive but produces results with better accuracy.

Since the simulation involves capturing instabilities which are changing in an unpredictable manner as a function of time, transient approach is preferred over steady state approach.

Atwood Number:-

It is a dimensionless number that is used to study hydrodynamic instabilities in density stratified flows.

An = ρ1-ρ2/(ρ1+ρ2)

where, ρ1 is density of heavier fluid
ρ2 is density of lighter fluid

For Air-water RT instability,

ρ1 = 998.2 kg/m^3 [Water density]

ρ2 = 1.225 kg/m^3 [Air density]

An = (998.2-1.225)/(998.2+1.225)

An = 0.998

For water- User define material RT instability,

ρ1 = 998.2 kg/m^3 [Water density]

ρ2 = 400 kg/m^3 [User define material density]

An = (998.2-400)/(998.2+400)

An = 0.428

Behaviour of RT Instability on variation of Atwood number:

RT evolves in three distinct stages: Linear stability, mushroom head (falling) or bubble (rising) formation, long term
evolution due to bubble merging and mixing. For low Atwood number, bubble and mushroom head shape is
relatively symmetrical, the disturbances or perturbations are linear for a larger time. In contrast, for large Atwood
number, bubble and mushroom head shape is less symmetrical, also it occurs very rapidly at an accelerating rate
which leads to a non linear growth rate.

Conclusion:

As we can see that the Atwood number for water - UD material is low so bubble and mushroom head shape is relatively symmetrical, the disturbances or perturbations are linear for a larger time. It takes too much time for simulation. In the other hand Atwood number for water - air is high , bubble and mushroom head shape is less symmetrical, also it occurs very rapidly at an accelerating rate. But it tskes less time for simulation as compared to low Atwood number.