Rayleigh–Taylor instability

Aim:-

To simulate a two-dimensional Rayleigh-Taylor instability problem using the classical method of Finite Element analysis of a multiphase model. Then perform grid dependency test on Rayleigh Taylor Instability for 2 different mesh size and compare the results for the variation in Atwood Number with the inclusion of user defined material in secondary phase.

Introduction:-

The Rayleigh Taylor instability is an instability of an interface between two fluids of different densities which occurs when the lighter fluid is pushing the heavier fluid in a gravitational field. The problem was governed by the Navier-Stokes and Cahn-Hilliard equations in a primitive variable formulation. In Multiphase flow, the dynamic variables like the velocity, viscosity, pressure and density are generally used to explain the movement of fluids under the gravitational field. Multiphase flows takes place in nuclear geothermal power plants, food processing industries and many others. One of the basic example is Rayleigh Taylor Instability. The Rayleigh Taylor Instability is a hydrodynamic instability arises when a high density fluid is placed over a low density fluid in a gravitational field. The RTI phenomenon was initially introduced by Rayleigh 1883 and later was observed for all accelerated fluids by Taylor 1950. When the density interface is disturbed, the hydrostatic pressure is generated, so the heavier fluid moves downwards and mixes with the lighter one, with the amplitude of perturbation initially growing exponentially but after sometime saturates and forms characteristic mushroom shape. When a heavier fluid is displaced downward with an equal volume of lighter fluid displaced upwards, it results in decrease of potential energy than it’s initial state. This RTI has been applied to a wide range of problems including the mushroom clouds like volcanic eruptions, nuclear explosions and plasma fusion reactors instability.

Theory:-

One of the classic examples of hydrodynamic instability is the mixing of two fluids that takes place if a heavy fluid initially lies above a lighter one in a gravitational field. The first to investigate this problem was Lord Rayleigh, who late in the last century considered the linear stability problem for inviscid fluids with various stratification profiles, including sharply stratified fluids. 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.

Water suspended atop oil is an everyday example of Rayleigh–Taylor instability, and it may be modelled by two completely plane-parallel layers of immiscible fluid, the more dense on top of the less dense one and both subject to the Earth's gravity. The equilibrium here is unstable to any perturbations or disturbances of the interface: if a parcel of heavier fluid is displaced downward with an equal volume of lighter fluid displaced upwards, the potential energy of the configuration is lower than the initial state. Thus the disturbance will grow and lead to a further release of potential energy, as the more dense material moves down under the (effective) gravitational field, and the less dense material is further displaced upwards. This was the set-up as studied by Lord Rayleigh. The important insight by G. I. Taylor was his realisation that this situation is equivalent to the situation when the fluids are accelerated, with the less dense fluid accelerating into the more dense fluid. This occurs deep underwater on the surface of an expanding bubble and in a nuclear explosion.

As the RT instability develops, the initial perturbations progress from a linear growth phase into a non-linear growth phase, eventually developing "plumes" flowing upwards (in the gravitational buoyancy sense) and "spikes" falling downwards. In the linear phase, the fluid movement can be closely approximated by linear equations, and the amplitude of perturbations is growing exponentially with time. In the non-linear phase, perturbation amplitude is too large for a linear approximation, and non-linear equations are required to describe fluid motions. In general, the density disparity between the fluids determines the structure of the subsequent non-linear RT instability flows (assuming other variables such as surface tension and viscosity are negligible here). The difference in the fluid densities divided by their sum is defined as the Atwood number, A. For A close to 0, RT instability flows take the form of symmetric "fingers" of fluid; for A close to 1, the much lighter fluid "below" the heavier fluid takes the form of larger bubble-like plumes.

This process is evident not only in many terrestrial examples, from salt domes to weather inversions, but also in astrophysics and electro-hydrodynamics. For example, RT instability structure is evident in the Crab Nebula, in which the expanding pulsar wind nebula powered by the Crab pulsar is sweeping up ejected material from the supernova explosion 1000 years ago. The RT instability has also recently been discovered in the Sun's outer atmosphere, or solar corona, when a relatively dense solar prominence overlies a less dense plasma bubble. This latter case resembles magnetically modulated RT instabilities.

Atwood number:-

The Atwood number is a dimensionless number in fluid dynamics used in the study of hydrodynamic instabilities in density stratified flows. It is a dimensionless density ratio related to the density of heavier fluid and the density of lighter fluid. (The difference in the fluid densities divided by their sum). Atwood number is an important parameter in the study of Rayleigh–Taylor instability. 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

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

A	Atwood number (dimensionless)
ρ1	Density of heavier fluid (kg/m3)
ρ2	Density of lighter fluid (kg/m3)

Pre-processing and meshing:-

The geometry is created as two square boxes 20 mm× 20 mm and converted to surface using Pull tool and share topology is turned on. The geometry was created in Spaceclaim and it involves no other importing of CAD models. The geometry is shared using share tool to compute the interference region between two square as shown in Figure 1. After completing geometry section it is further processed to meshing section.

Figure 1, Geometry

In the following section, we analysed the R-T instability phenomena for two immiscible incompressible fluids having different densities./ The change is that the density of the secondary phase is varied with the inclusion of user defined material.

We create a simple hexa grid mesh with two different meshing schemes as shown in Figure 2. This pre-processing phase of the workflow is extremely important because it determines the simulation potential and limits the CFD results.

Figure 2, Mesh

Mesh dependency test:-

For every numerical simulation we have to check for a Mesh dependency of solution. Making a mesh too much fine can lead to consume too much resource and round off error in simulation. But also making mesh course can make us compremise with simulation accuracy. So this test to define the limit which up to we refine our mesh until there be no change in result we ae can observe. For our present case we choose a two mesh case as mentioned below.

• Case A. Mesh details: Baseline mesh- 0.5mm Element size, Elements – 3200, Nodes – 3321
• Case B. Mesh details: Mesh - 0.2mm Element size, Elements – 20000, Nodes – 20301

As the mesh size is increased we get more defined results which we can see from the video link and it took a long time for the simulation to reach appropriate results. The shape of the globules deforming as flow time progresses actually varies with the meshing. The volume fraction in simulations on different meshes at a particular flow time shows the interface thickness is smaller on finer mesh.

Solver setup:-

Based on the mesh convergence suitable mesh size is selected for all the cases performed . now as we reach at defining numerical section in Fluent. Here we choose a pressure based, transient state is selected with gravity enable in -Y direction with constant 9.81 m/s. Multiphase Model is turned on along with Volume of Fluid and number of Eulerian phases is chosen as 2 and Implicit scheme is selected. Viscous laminar model is selected.

Material is selected as air and water then for Atwood number variation we define one user defined material their properties are defined in below Table 1.

Material	Density	Viscosity
Air	1.225 𝑘𝑔/𝑚3	1.7894e-05 𝑘𝑔/𝑚-s
Water liquid	998.2 𝑘𝑔/𝑚3	0.001003 𝑘𝑔/𝑚-s
User defined	400 𝑘𝑔/𝑚3	0.001 𝑘𝑔/𝑚-s

Table 1

Case 1-In the following section, we analysed the R-T instability phenomena for two immiscible incompressible fluids having different densities (Primary phase, 𝜌1 = 1.225 𝑘𝑔/𝑚3 as air and Secondary phase, 𝜌2 = 998.2𝑘𝑔/𝑚3 as water liquid) with different viscosities(𝜇1 = 1.7894e-05 𝑘𝑔/𝑚-s ,𝜇2 = 0.001003 𝑘𝑔/𝑚-s).

Case 2- In the following section, we analysed the R-T instability phenomena for two immiscible incompressible fluids having different densities (Primary phase, 𝜌1 = 400 𝑘𝑔/𝑚3 as air and Secondary phase, 𝜌2 = 998.2𝑘𝑔/𝑚3 as user defined material) with different viscosities(𝜇1 = 0.001 𝑘𝑔/𝑚-s ,𝜇2 = 0.001003 𝑘𝑔/𝑚-s).

After setting primary and secondary phases we select Pressure velocity coupling scheme as Coupled because couples scheme shows good result in transient case simulation. For pressure, PRESTO is selected. For the momentum second order upwind scheme is used. For the transient formulation first order implicit scheme is

chosen for better results and for stability and convergence of the residuals. then solution is initialized with hybrid initialization and in water volume fraction we choose 0 as we initially want air to fill in entire volume. After that we need to patch the volume fraction of the two phases is performed before the start of simulation. Phase used for patching was water and volume fraction of 1 is taken for water surface and 0 for air surface. After that we start calculation for time step size of 0.005 second and number of time step be 900 so total we run the simulation for 4.5 seconds. Max iteration time step be 20. For both cases 1 and 2 the video link for volume fraction is given below.

For Case 1- https://drive.google.com/file/d/1NYYvves2zlZXgU0SZmJtQBwrm9f0iuIW/view?usp=sharing

For Case 1 Refined mesh-https://drive.google.com/file/d/1UOOV-ohZdzJKfiAC08IGvSCNaX4Xd1rS/view?usp=sharing

For Case 2-https://drive.google.com/file/d/1WwGyQkRIH8r53qBKh_8bj9kKBkYf4sR5/view?usp=sharing

Result:-

At first the fluids behave linearly and we can observe small linearity sinking going on between them, slowly the chunky bubbles from the higher density fluid on top settles down. This can be identified from the highest mesh case resulting in finest behaviour of the instability for this condition. At t = 0, the small enough, initial perturbation in flow can be described by a smooth function. We observe the disturbance in flow from smooth to turbulent growth has been observed. We also observed the transition of a single globule into lots of deformed globules with with lots of turbulence along with some breakup and penetration of heavier fluid into the lighter one. These breakups which seems at start vanishes as time proceed as to last stage, showing the phenomena of Rayleigh-Taylor instability. From Figure 3 the residual plot shown a jump its because when nest timestep has been started nd then residual becomes fof same amplitude ,the difference between their change is resduced as solution goes on,in transient solution will not converge so from this stste we can conclude that solution has converged to satisfactory level. At the middle stage of video, it is clear that the falling fluid is rolled up. In our result, finally the heavier fluid has settled down completely at last or we can say it is in the verge of complete settling. IN above link it shows the 2D-Rayleigh Taylor Instability simulation having the contours of two phase mixture model. It shows the penetration of heavier fluid into lighter fluid with the progress of flow time.At later time we observe the formation of low density bubbles and its rise into the high density phase whereas the high density globules sinks at the edges of the domain of the simulation into the light fluid. With the time evolution, the interface gets smeared due to numerical diffusion within the interfaces. The shape of the globules deforming as flow time progresses actually varies with the meshing. The volume fraction in simulations on different meshes at a particular flow time shows the interface thickness is smaller on finer mesh.

Figure 3, Residual Plot

Atwood number calculation:-

The Atwood Number is a dimensionless number in Fluid Dynamics. It is used in the study of instabilities of the dynamics of the fluid flow .It gives the relation between density of higher fluid and the density of lighter fluid. For Atwood no. close to 0, RT instability flows takes the form of symmetric "fingers" of fluid. It was seen that with larger Atwood Number, larger amplitude of perturbations are observed. It is due to the larger density ratio of the two fluids, the stronger is the force to develop the interface instability.

For our Case 1 Atwood number be equal to 0.99, and for our Case 2 Atwood number be equal to 0.42. 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 a little influence. As we see for low Atwood number the interference start breaking slowly and then spikes of small size start penetrating less denser fluid ,the size of plumes are smaller. Solution that large time to stabilize. 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 seen in the Video 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. The spikes of denser fluid start penetrating the lighter fluid spikes are of small size but the plume are of larger size and then it start rotating expanding and due to density the dense fluid rotates downward and plume vanish as times pass on. The flow structure of bubbles and globules are quite different both for the variation in meshing and also with the variation in Atwood Number. 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 RTI in the steady growth stage is shorter. We can observe that the compressibility on the bubble velocity is strong at larger Atwood Number.

Conclusion:-

In above simulation we adopted Volume of fluid method to study the evolution of the interface during Rayleigh Taylor Instability. This numerical simulation is analysed by the variation in Atwood Number. For our case 1 we get Atwood number equal to 0.99 and for Case 2 it is 0.42 In lower Atwood numbers, spikes and globules were still present at the end of the simulation which was not there in Atwood Number 0.9 and it take much time to get settled . From the grid dependency we also observed that the interface thickness reduces with the refined mesh and when run with smaller time steps. For Transient simulation, compared to steady state, has higher demand in terms of computational time. Transient solver better captures the flow and can be used to run at smaller time steps which is with a time step size of 0.005. Smaller time steps are required to observe the instabilities and can be only possible through transient simulation. 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.