Week 5 - Rayleigh Taylor Instability

Aim: The aim is to understand the Rayleigh Taylor instability phenomena and how it occurs.

Objective :

• Some practical CFD models have been based on the mathematical analysis of Rayleigh Taylor waves.
• Explain how these mathematical models have been adapted for CFD calculations.
• Perform the Rayleigh Taylor instability simulation for 2 different mesh sizes with the base mesh being 0.5 mm. Compare the results by showing the animations
• Explain why a steady-state approach might not be suitable for these types of simulations.
• Run one more simulation for refined mesh with water and user-defined material(density = 400 kg/m3, viscosity = 0.001 kg/m-s).
• Define the Atwood Number. Find out the Atwood number for both cases and explain how the variation in the Atwood number in the above two cases affects the behavior of the instability.

RAYLEIGH TAYLOR INSTABILITY

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. 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, and Supernova explosions in which expanding core gas is accelerated into denser shell gas. There are various models of Rayleigh Taylor Instability in which fluid interfaces behave differently.

Kelvin-Helmholtz or Rayleigh Taylor Hybrid Model.

• Kelvin-Helmholtz Instability (KHI) is a hydrodynamic instability in which immiscible, incompressible, and in-viscid fluids are in relative and irrotational motion.
• In KHI, the velocity and density profiles are uniform in each fluid layer, but they are discontinuous at the (plane) interface becomes an unstable vortex sheet that rolls up into a spiral.
• KHI is observed in a wide range of nature. Oceankelvin - Helmholtz waves, Kelvin-Helmholtz Instability in clouds, Clumping in supernova remnants.

Potential Flow models of Rayleigh-Taylor and Richtmyer-Meshkov for bubble fronts.

• In the potential Flow model of Rayleigh-Taylor and Richtmyer-Meshkov bubbles of an interface between an incompressible fluid and a constant supporting pressure (Atwood number A= 1) are presented.
• In this model, which extends the work of Layzer, the ordinary differential equations for the bubble height and curvatures are obtained by considering the potential flow equations near the bubble tips.
• The model is applied to two-dimensional single-mode evolution as well as two-bubble competition, for both the Rayleigh-Taylor(RT) and the Richtmyer-Meshkov (RM) instabilities, the latter treated in an impulse approximation.
• The model is applied to the evolution of three-dimensional modes with an initial rectangular geometry. The model yields the aspect ratio dependence of the early nonlinear stages, in agreement with third-order perturbation theory.

Granular Rayleigh-Taylor instability

• A granular instability is driven by gravity in which, the instability arises as grains fall in a closed Hele-Shaw cell where a layer of dense granular material is positioned above a layer of air.
• The initially flat front defined by the grains subsequently develops into a pattern of falling granular fingers separated by rising bubbles of air.
• A transient coarsening of the front is observed right from the start by a finger at the center of the rising bubbles which we can observe in the below picture.

Lattice Boltzmann Model:

• The Lattice Boltzmann method was originated from Ludwig Boltzmann’s kinematic theory of gases.
• The fundamental idea of LBM is that gases/fluids can be imagined as consisting of a larger number of small particles moving with random motions. The exchange of momentum and energy is achieved through particle streaming and billiard-like particle collision.
• The Lattice Boltzmann method (LBM) has emerged as a promising tool for modeling the Navier-Strokes equations and simulating complex fluid flows.
• LBM is based on microscopic models and mesoscopic kinetic equations. In some perspectives, it can be viewed as a finite difference method for solving the Boltzmann transport equation.
• Moreover, the Navier-Strokes equations can be recovered by LBM with a proper choice of the collision operator.

K L Turbulence Model:-

• The k-L turbulence model, where k is the turbulent kinetic energy and L represents the turbulent eddy scale length, is a two-equation turbulence model that has been proposed to simulate turbulence induced by Rayleigh-Taylor (RT) and Richtmyer Meshkov (RM) instabilities, which play an important role in the implosions of inertial confinement fusion (ICF) capsule targets.
• The model describes the dominant eddies in the mixing zones with evolutionary equations for their characteristic dimension (L) and energy per unit mass (K=(v^2)/2).
• The equations are based on buoyancy-drag models for RT and RM flows.

WORKBENCH SETUP

After getting the theoretical overview of the Rayleigh Taylor Instability, let us discuss the numerical setup of two-dimensional simulation on Ansys Fluent.

There were Four types of simulations performed in this project. The first two are for mesh dependency and to understand the differences in simulation concerning mesh grids.

The third case is performed to study the effect of Atwood number in the Rayleigh Taylor method. In this case, the density of the secondary phase is altered with the incorporation of user-defined material resulting in variation of Atwood Number.  

The results of the grid dependency study and the variation of Atwood Numbers are carefully observed and reported after simulation in the present literature.

PREPROCESSING

Geometry setup

At the starting point, geometry, as two square 2D boxes are created in the Space claim and converted into surfaces using the Pull tool. Also, to share mesh in-between two faces share topology is executed on the geometry.

Meshing

The surface model was further processed for meshing and the final results were taken for three different simulations. Mesh for the First simulation was created using 0.5mm mesh grid size.

The nodes are 3321 and elements 3200

Here 0.2 is taken into consideration for a finer mesh.

Fluent Case setup

• A pressure-based solver and Absolute velocity are used.
• Transient-state Calculation is performed.
• Gravity is enabled in -ve Y-axis

Cell zone selected as Fluid-air;water-surface:water-liquid. 

Patching of the volume fraction of the two phases is performed before the start of the simulation. The phase used for patching was water and a volume fraction of 1 is taken for water surface and 0 for air surface.

 POST-PROCESSING

The evolution of the Rayleigh Taylor Instability follows four main stages. In the first stage, the perturbation amplitudes are small when compared to their wavelength, the equations of motion can be linearized, resulting in exponential instability growth. In the early portion of this stage, a sinusoidal initial perturbation retains its sinusoidal shape.

 

 After the end of the first stage, when non-linear effects begin to appear, one observes the beginnings of the formation of the ubiquitous mushroom-shaped spikes (Fluid structures of heavy fluid growing into light fluid) and bubbles (fluid structures of light fluid growing into heavy fluid). The growth of the mushroom structures continues in the second stage and can be modeled using buoyancy drag models, resulting in a growth rate that is approximately constant in time. At this point, nonlinear terms in the equations of motion can no longer be ignored.

 The spikes and bubbles then begin to interact with one another in the third stage. Bubble merging takes place where the nonlinear interaction of mode coupling acts to combine smaller spikes and bubbles to produce larger ones. Also, bubble competition takes place, where spikes and bubbles of smaller wavelengths that have become saturated are enveloped by larger ones that have not yet saturated.

This eventually develops into a region of turbulent mixing, which is the fourth and final stage in evolution. It is generally assumed that the mixing region that finally develops is self-similar and turbulent, provided that the Reynolds number           

Case 1 - Animation - Water and Air (Coarse Mesh)

Case 2 - Animation - Water and Air (Fine Mesh)

Observation:-

Observation:-

At the very beginning, the small enough, initial perturbation in a flow can be described by a smooth function. Then the disturbance in a flow from smooth to turbulent growth has been observed in further time steps.

We also observed the transition of a single globule into lots of deformed globules with lots of turbulence along with some breakup and penetration of heavier fluid into the lighter one. These breakups which seem after some time steps, show the phenomena of Rayleigh-Taylor instability.

Before settlement, one can clearly observe that the falling fluid is rolled up. In our result, finally, the heavier fluid has settled down completely at or we can say it is on the verge of completely settling.

 Differences in Animation of coarse and fine mesh

First before making an animation to simulate the results the fine mesh required a very small time step (0.0005) than the coarse mesh(0.002) which required more time to simulate the results for the fine mesh case.

The interface thickness between two fluids is reduced with the fine mesh which provides a little sharper view in animation.

Fine mesh capture more detail in animation than a coarse one. Coarse mesh just shows the flow pattern at the interfaces of two fluids but in fine mesh animation, we can see all the small bubble formations with sharp flow patterns in a very detailed manner.

 

Steady State simulation

In the steady-state simulation, it solves the problem iteratively to get final results. Therefore, it never shows the changes in flow patterns with respect to time. Study of Raleigh Taylor Instability is the observation of instability of an interface between two fluids with a different density which occurs when the lighter fluid is pushing the heavier fluid. For this, we need to capture the frames with the changes in the fluid as time passes which we can only observe in the Transient state.

 

Following animation is an observation of Raleigh Taylor Instability between two fluid simulated using Steady-state solver.

 

 

Case 3 - Animation - Water and User Define material (Fine Mesh)

 

Observation:-

For the Rayleigh Taylor Instability simulation with the user define the material, we can observe this phenomenon with a much slower time frame which means we can run this simulation with a greater number of initial time steps than the above two cases.

we observe a broadening of the interface as the diffusion takes place. For a time, the simulation enters the non-linear regime with the increase in the amplitude of perturbation. At a later time, we observe the formation of low-density bubbles and it rises into the high-density phase whereas the high-density globules sink at the edges of the domain of the simulation into the lighter fluid. Finally, like in previous cases, the heavier fluid settled down at the bottom.

 

Case 4 - Animation - Steady-State - Water and Air (Fine Mesh)

 

 

We can see, steady-state models try to achieve a condition of an equilibrium state which is also independent of time. Hence it is not suitable for is types of problems.

Atwood Number:-

The 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=ρ1-ρ2ρ1+ρ2</math>">A=ρ1−ρ2ρ1+ρ2�=�1-�2�1+�2

Where

rho1 = density of the heavier fluid

rho2 = density of the lighter fluid

 

For Case1 and Case2

ρ1</math>">ρ1�1 = density of heavier fluid (Water) = 998.2 kg/m3

ρ2</math>">ρ2�2 = density of lighter fluid (air) = 1.225 kg/m3

 

                                                   A=998.2-1.225998.2+1.225</math>">A=998.2−1.225998.2+1.225�=998.2-1.225998.2+1.225

                                                   A=0.99754</math>">A=0.99754�=0.99754

It was seen that with a larger Atwood Number, a larger amplitude of perturbations is observed. It is due to the larger density ratio of the two fluids, the stronger is the force to develop the interface instability.

For Case3

ρ1</math>">ρ1�1 = density of heavier fluid (Water) = 998.2 kg/m3

ρ2</math>">ρ2�2 = density of lighter fluid (user Define) = 400 kg/m3

 

                                                   A=998.2-400998.2+400</math>">A=998.2−400998.2+400�=998.2-400998.2+400

                                                   A=0.42783</math>">A=0.42783�=0.42783

Observation:-

If we compare the high Atwood number with the decreased Atwood number in case 3, we can clearly observe that bubbles and spikes formation are still present but not for the same time flow. For the case of smaller Atwood Number, the density difference between the two sides of the interface is small and the density distribution of the upper and lower layers is nearly symmetrical. The initial density stratification plays a dominant role and the expansion compression effect has little influence.

With the increased Atwood Number close to 1 the stabilization effect of initial density stratification decreases and the instability caused by the expansion-compression effect becomes more significant.

It can be observed in the animation that the larger Atwood numbers, the larger amplitude of perturbation. It takes place due to the larger density ratio of two fluids is, the stronger force to develop the interface instability.

With the increase of Atwood number, the time of spikes or globules of the lighter phase fluid arriving at the bottom of the computational domain is accelerated, then the duration of Rayleigh Taylor instability in the steady growth stage is shorter. We can observe that the compressibility on the bubble velocity is strong at a larger Atwood Number.

The bubble height is approximately a quadratic function of time at the potential flow growth stage. The average bubble acceleration is nearly proportional to the square of the Mach Number of Atwood Number close to 1. For Atwood No. close to 1, the much lighter fluid below the heavier fluid takes the form of larger bubble-like plumes.

 

 

 

Conclusion

 

In the report, the classical method of Finite Element analysis was adopted to study the evolution of the interface to observe Rayleigh Taylor Instability.

From the grid dependency, we also observed that the interface thickness reduces with the refined mesh and when run with smaller time steps.

Transient simulation, compared to steady-state, has higher demand in terms of computational time. The transient solver better captures the flow and can be used to run at smaller time steps.

The steady-state solution only simulates the problem, using the iterative technique, to get the solution at its equilibrium state, and therefore, it is not preferable to use a steady-state solver to simulate the solution of Rayleigh Taylor Instability.

Since it is a time marching solution, we expect more accuracy from transient simulation for which Transient simulation is the ideal setup to run the calculation for Rayleigh Taylor Instability.

Smaller time steps are required to observe/capture the instabilities and can be only possible through transient simulation.

• For the case of smaller Atwood Number, the density difference between the two sides of the interface is small and the density distribution of the upper and lower layers is nearly symmetrical. The initial density stratification plays a dominant role and the expansion compression effect has little influence and vice versa.

References:-

[1] The Rayleigh- Taylor instability of viscous fluid layers - https://aip.scitation.org/doi/10.1063/1.869283

[2] Potential Flow models of Rayleigh-Taylor - https://aip.scitation.org/doi/10.1063/1.868391

[3] kelvin-helmholtz instability- https://www.sciencedirect.com/topics/earth-and-planetary-sciences/kelvin-helmholtz-instability

[4] Granular Rayleigh-Taylor instability - https://www.researchgate.net/publication/337621555_Granular_Rayleigh-Taylor_instability