Rayleigh Taylor Instability Challenge

1. Aim

To simulate Rayleigh-Taylor instability using VOF method in ANSYS Fluent

2. Introduction

Rayleigh-Taylor instability is a type of instability occurring at the interface of two fluids with difference in densities in the presence of gravitational field. The fluid with higher density is located above the less dense fluid. Due to the gravitational pull, the higher density fluid is pushed into the lower density fluid resulting in peculiar patterns of spikes and bubbles transitioning into unstable mixture. RT instability can be observed in several real-life applications such as the behaviour of water suspended above oil, mushroom clouds volcanic eruptions and atmospheric nuclear explosions, supernova explosions.

Fig 1- RT instability in nuclear explosions and supernova (Ref- Google images)

3. CFD models based on mathematical analysis of Rayleigh Taylor waves

To understand the numerical modelling methods for simulating RT instability, it is essential to define the free surface flows. A free surface flow is the surface of a fluid flowing that is subjected to both zero perpendicular normal stress and parallel shear stress. In presence of surface tension, the instability formed is called Rayleigh instability, an example of which could be water jet from hose pipe breaking into stream of droplets.

The different models to simulate free surface flows are:

3.1. Marker and Cell method

It is one of the first developed numerical models based on Eulerian grid of control volume, the grid being fixed in space throughout the simulation. Using a set marker of particles to flag the cell zones present in the domain under study, it classifies cells into full, empty and surface cells. Full cells are the ones with marker particles, empty cells are ones without marker particles and surface cells are the adjacent ones to the empty cells. Note that the cells defined do not have mass or volume.

3.2 Langragian grid methods

In this method, the properties of the fluid are calculated using a common grid which moves along with the fluid. Since both the grid and fluid flow together, the effects of surface tension and shear stress are negated and we get the characteristics of free stream flow. However, the presence of uneven surfaces in the domain would result in accurate result as the common grid cannot capture flow effectively.

3.3 Level set method (LSM) method

Level set method measures the characteristics using level set functions. The idea is to intersect the moving fluid flow with a planar surface and plot it contour using a PDE. Level sets functions are defined such that positive contour are obtained for one fluid and negative contours for the other fluid.

3.4 Volume of flow (VOF) method

Volume of flow method is an Euler method where the grid can be stationary or non-stationary. It is based on introducing a fractional function C. This function is determined using integration of fluid characteristics in the domain. The VOF method is not a stand-alone method and requires solution of Navier-Stokes continuity equation for the domain under study.

The function C is defined as a continuous function with 0 ≤ C ≤ 1. The physical significance of C being cell is empty for C=0 and filled completely for C=1.

Due to average values calculated across cell zones, the results obtained show diffusion at the interfaces. This is resolved using reconstruction techniques such as:

• Line techniques
• Donor Acceptor Techniques
• Higher order Differencing Techniques

4. Methodology

The Rayleigh-Taylor instability is simulated for a case of two fluids- water and air in 2D. The model is generated using Spaceclaim and meshed for three different mesh element sizes- 0.5 mm, 0.25 mm and 0.1 mm.

The setup is patched and run for multi-phase laminar model to simulate the fluid flow for a transient state. The volume fraction of water contour is studied over the course of time to understand the stages of development

5. Pre-processing & Solver settings

5.1 Geometry

Two components with surface pull of dimensions 20x20 mm are generated as shown in figure. The common edge is shared between the two components.

Fig 2- Geometry

5.2 Meshing

The baseline mesh is generated for element size 0.5 mm and is refined for two cases of element sizes- 0.25 mm and 0.1 mm.

Fig 3- Mesh of coarse, moderate and fine type

Table 1- Mesh parameters

Mesh	Coarse	Moderate	Fine
Element size	0.5 mm	0.25 mm	0.1 mm
Number of Nodes	2415	9453	58315
Number of cells	2312	9248	57804

5.3 Setup

Time: Transient state

Type: Pressure based

Multi-phase: Volume of fluid (VOF) method- Implicit type

Number of Eulerian phases: 2 (air and water)

Viscous: Laminar

Initialization: Standard

Patch: Volume fraction of water is zero in lower component and 1 in upper component

Time step size is varied from 0.005, 0.001 and 0.0005 seconds for the coarse, moderate and fine case, respectively, to avoid floating point error. The case is run for about 2.5 seconds.

5.3.1 Why steady state analysis cannot capture RT instability

As steady state analysis is time independent, the solution state simulated captures stability (equilibrium) of the two phases without capturing the fluid phenomena occurring at the interface of the two phases. Hence, steady state analysis cannot be used to observe the Rayleigh-Taylor instability

6. Results

6.1 Animation

It can noted from the animation for the three cases that mesh refinement helps in effective capturing of spikes and bubbles formed during the mixture of two phases- water and air.

6.1.1 RT Instability coarse mesh (size 0.5 mm)

6.1.2 RT Instability moderate mesh (size 0.25 mm)

6.1.3 RT Instability fine mesh (size 0.1 mm)

6.2 Stages of development

Fig 4- Stages of development

The phenomenon of RT instability involves four stages of development as presented in the figure. In the first stage, the perturbation amplitudes are less compared to their amplitudes. This generates a sinusoidal shape at the interface. At the next stage, we observe formation of bubbles and spikes in the domain based on the densities of the two phases. For this case, the lower density fluid pushes the higher density fluid to form bubbles at the interface. This result in spikes of higher fluid as it enters the lower component. These mushroom-like structures grow in size over time. In the third stage, bubbles compete and coalesce to form a larger mushroom-like structure before reaching the final stage of instability. At the final stage, the fluid flow involves inter-mixing which may become turbulent based on the Reynolds number. After several time steps, the domain attains the state of equilibrium.

6.3 Atwood number and validation of results

As the nature of the fluid flow at the interfaces is dependent on the densities of the two phases, it is quantified by defining the Atwood number (A). It is calculated as the difference in the fluid densities divided by their sum.

a) For A ~ 0, RT instability flows take the form of symmetric "fingers" of fluid;

b) For A ~ 1, the lower density fluid below the heavier fluid takes the form of larger bubble.

For the given case, ρ₁ = 997 kg/m³ and ρ₂ = 1.225 kg/m^3. Hence, Atwood number is calculated to be

A= 0.998 ~ 1

Since Atwood number is close to 1, we observe large bubbles of lower density fluid- air and spikes of higher density fluid-water.

7. Conclusion

The study successfully studies the Rayleigh-Taylor Instability in two-phase domain by using the volume of fluid method in ANSYS Fluent. The effect of mesh size refinement in improved flow structures is established using animation and the stages of development are explained. The observed phenomenon is validated by calculating the Atwood number for the given case.