Rayleigh Taylor Instability Challenge

Rayleigh Taylor instability:

• When two fluids of different densities superposed one over the other (or accelerated towards each other); the instability of the plane interface between the two fluids, when it occurs, is called the Rayleigh-Taylor instability.
• If the horizontal surface of a liquid at rest under gravity is displaced into the form of regular small corrugations and then released, standing oscillatory waves are produced. Theoretically, a liquid could exist in a state of unstable equilibrium with a flat lower horizontal surface supported by air pressure.
• The physics involved are rather intuitive: “heavy stuff falls, light stuff rises”.
• For a heavy fluid over a light fluid, the heavy stuff wants to fall and the light stuff wants to rise; however, for a perfectly flat interface there are no avenues for this vertical motion to occur. When this perfectly flat interface is displaced (perturbed) the waves generated (as seen in the stable case) provide the necessary means for the light fluid to rise (in the peak of the wave) and the heavy fluid to fall, in the valley of the wave.

CFD Models based on Rayleigh taylor instability :

1. Mass diffusion model:

The mass diffusion equation describes the transport of mass and energy due to a model viscosity and hence drives actual mixing in the fluids. Even if any configuration is Rayleigh–Taylor unstable, diffusion processes could be responsible for significant mixing if the viscosity is sufficiently large.

2. Single Fluid Model:

Potential flow model:

Used by Layzer (1955) and Hecht et al. (1994), this model is valid for all instability stages from the linear stage through the early nonlinear to the asymptotic but limited to A =1. In this model,single-mode perturbation and expand the flow equations to second order around the bubble tips, hence obtaining an ordinary differential equation for the bubble velocity.

Buoyancy drag model:

Buoyancy–drag models valid for every Atwood numbers, but limited to the nonlinear 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, buoyancy, and drag forces. These models cannot model multiple mixing interfaces; cannot be easily extended to two and three dimensions; and, as a rule, do not address demixing, also known as counter gradient transport, i.e., reduction of total fluid masses within the mixing zone. To address the above problems, two fluid (or multi-fluid) models have to be used.

Besnard-Harlow-Rauenzahn (BHR):

It is a more advanced version of a single fluid model. The BHR model is based on 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 Reynolds stresses, turbulent mass-flux velocity, density fluctuations, and the turbulence 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 Reynolds averaging of the governing equations and resultant terms.

Two-equation K-L turbulence model:

The k–L turbulence model evolves equations for the turbulent kinetic energy k and the turbulent length scale L. The K-L model was proposed by Dimonte and Tipton for describing the turbulent self-similar regime of RTI induced mixing. The starting point for deriving the model equations are the buoyancy-drag models for the self-similar growth of RT instabilities.

Three-equation K-L-a turbulence model:

The three-equation model of Morgan and Wickett was developed as an extension to the two-equation K-L model of Dimonte and Tipton by including a third equation for the turbulent mass-flux velocity.

Four-equation K-L-a-b turbulence model:

It is an extension of a K-L-a model by adding an evolution equation for the density-specific volume covariance. This model can accurately reflect the effects of the changes in the density fluctuations.

Simulation Analysis :

Geometry of simulation :

• In this simulation, we are about to analyse the Rayleigh Taylor Instability that occurs between two different fluids, namely water & air.
• We are simulating a 2D analysis over a surface on XY plane.
• Thus we have a geometry of two phases having a common edge.
• Initially each phase consists of different fluids, which will be disturbed by gravity for the simulation to take place.

Computational Meshing :

• As we are simulating in 2D for two square phases, the mesh need not have to be triangonal.
• The best suited mesh type for this simulation is multiphase meshes.
• We are glad to compare the results for different grid sizes. Thus we are going to take 3 different cases by varying the element size of the mesh.

Case (1): Case (2): Case (3):

Element size = 0.5 Element size = 0.3 Element size = 0.2

Setup & Simulation :

• As we are calculating the instability in the fluids, we are simulating the problem with transient condition.
• The Gravity force is enabled along the negative Y direction.
• Two materials are selected from the Fluent database.
• Primary & secondary phases are specified in the material phases.
• After the initialisation of the simulation, Patching is done to differentiate one phase from the other.
• Animation of the countour is captured as the itration goes on.

Initial Contour of the surface :

• The red surface denotes the water & the blue surface denotes the air.
• The two phases of fluid is seperated from each other by a common edge.
• As the seperation is disturbed (i.e) Removed , the gravitational force enabled will tend to form instability to the fluids.
• Let us compare the animatons of the instabilities caused by the gravity for different mesh sizes.

Surface Animation :

Case (1):

Case (2):

Case (3):

Comparison of the simulation :

• For better conclusions, let us compare the snaps at different itrations.
• The number of itrations completed is denoted below the snaps.
• For the same number of itrations, let us compare the contour results for all 3 cases.

Case (1):

15 50 100 200 500 700

Case (2):

15 50 100 200 500 700

Case (3):

15 50 100 200 500 700

Conclusion based on the comparison :

• As the mesh size is deceased, we can see that the instability caused is faster.
• The mixing of the fluids starts quicker when there are higher number of grids.
• We can find more number of bubble formation on the refined mesh Rather than the baseline mesh.
• As the smaller grids cover more area, even small bubble formations of instability caused is calculated.
• The settlement of the higher density fluid is quicker in refined mesh.
• But at the same time we can see some bubbles during the settlement in refined mesh, whereas in baseline mesh there are no visible bubbles during the settlement of the fluid.
• Since we can get more details from smaller grid size, we can conclude that smaller grids can provide more accurate results than the baseline mesh.
• We can also see a much smoother contour in the refined mesh.

Why steady state is not suitable? :

• As we are calculating the instabilities between two fluids, steady state cannot be used.
• This is because , if we use steady state analysis the residual values during the itration increases.
• This causes a huge error in the result.
• By taking steady state, we may neglect some important calculations hich may lead to error in the output.

Simulation Using User defined Material :

• For further clarrifications on the instability, we can try the simulation for different fluids.
• Let us have water as it is and introduce an user defined primary phase with the details given in the problem.
• The same conditions are provided for this simulation also.
• This simulation is modified from the refined mesh (i.e) from case(3).
• Now we can compare results of this simulation with the results obtained from case(3).

Animation of surface instability :

Comparison with Case (3) :

Case (3):

15 100 200 500 700 1000

User Defined Material :

15 100 200 500 700 1000

Conclusion :

• From the 15th itration we can clearly see that the instability is not started with the user defined material, where as the instability in case(3) has started to form penetration of the fluids.
• This is because, the user defined material has a larger density value than that of air.
• As the density is larger, the resistance for the penetration of water also gets large.
• Thus it takes some time to start the penetration in this case.
• We can also see that the water is being pushed further by the user defined material than that of air.
• As the water starts to penetrate, some amount of volume will be taken by the water from the primary material. Thus the primary material starts to push the water in the direction opposite to the flow of water.
• Because of higher density, we can see the push from user defined material is higher.

Atwood Number :

• Atwood number, abbreviated as A, is a dimensionless number that describes density difference between two adjacent fluids with a common interface.

where

= density of heavier fluid

= density of lighter fluid

• Atwood number is an important parameter in the study of Rayleigh–Taylor instability and Richtmyer–Meshkov instability.
• In Rayleigh–Taylor instability, the penetration distance of heavy fluid bubbles into the light fluid is a function of acceleration time scale.
• For Atwood number close to 0, RT instability flows take the form of symmetric “fingers” of fluid; for Atwood number close to 1, the much lighter fluid “below” the heavier fluid takes the form of larger bubble-like plumes.
• In our case Atwood number for water and air is

A1 = 0.9975

Atwood number for water & user defined material is

A2 = 0.4277

• As we compare the start of instability we can see the difference between the 2 cases.

User Defined Material : Air :

• As discussed we can see that with air our atwood number is closer to 1. Thus we are getting only bubbles and we cannot see any spikes poking out of the bubbles.
• In the user defined material case we can see projections in the shape of mushroom heads, as the atwood number is closer to zero.