Study of Rayleigh Taylor Instability with the help of CFD simulation

: The purpose of this paper is to simulate a two-dimensional Rayleigh-Taylor instability problem using the classical method of Finite Element analysis of a multiphase model using ANSYS FLUENT 19.2 . The governing equations consist of a system of coupled nonlinear partial differential equations for conservation of mass, momentum and phase transport equations. The study focuses on the transient state simulation of Rayleigh Taylor waves and subsequent turbulent mixing in the two phases incorporated in the model. 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. The Cahn- Hilliard equations were used to capture the interface between two fluids systems. The objective of this article is to perform grid dependency test on Rayleigh Taylor Instability for 2 different mesh size and compare the results for the variation in Atwood Number. The results were validated with the observations from previous published literatures. 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.


Introduction
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.
Presently large number of numerical techniques were introduced to simulate RTI including boundary integral methods ( Baker et al. 1984), front tracking methods ( Liu et al. 2007, Terashima andTryggvason 2009), volume of fluid methods (Hirt and Nichols 1981, Gopala and BGMvan 2008,Raessi et al. 2010) and phase field methods (Ding et al. 2007, Chiu and Lin 2011, Lee et al. 2011, Kim 2005 Liu et al. 2015, Xu et al. 2017, Cubeddu et al. 2017). Due to numerous experimental studies conducted by (Read 1984, Marinak et al. 1998) it is understood that RTI is a spatial process. Due to the complexity of the implications of 3D simulations, 2D simulations are studied extensively in recent years. Zhou et al. 2018 implemented 2D Large Eddy Simulation(LES) for investigating the dynamics and the evolution of Rayleigh Taylor bubbles from W-shaped, sinusoidal and random perturbations. The classical numerical methods such as the Finite Difference method, Finite Volume method, Finite Element method has been extensively used to solve physical problems in science and engineering but the implementation of these classical methods has certain drawbacks. The main drawback is that these classical methods is required to build a mesh in the computational domain and on the boundaries to solve the problem.
However the present study is based on the Finite Element method to simulate the interface evolution of RTI. Here in this study the 2D geometry is itself constructed in Spaceclaim and a computational mesh is generated. The meshing file is imported in fluent and a grid dependency test was performed and the results of the interfacial evolution of RTI is captured at different flow times and compared. The governing equations for the RTI were the modified form of both Navier-Stokes equation and Cahn-Hilliard equation. The spatial discretization involves least square cell based for gradient and for time derivative, first order implicit scheme was used in order to solve the governing equations. SIMPLE scheme was employed to overcome the pressure velocity coupling. In this study 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 in the present literature.

Governing Equations
The governing equations for The RTI problem are the Navier-Stokes equations, the continuity and the Cahn-Hilliard equations which are the non-dimensional form.
The above equations can be reduced into non-dimensionalized form as follows

Preprocessing and Mesh
At first, geometry is created as square boxes 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 surface model was further processed for meshing and the results were taken for three different meshing schemes. This preprocessing phase of the workflow is extremely important because it determines the simulation potential and limits the CFD results.
• Case A.

CFD Case Setup
• Pressure based solver and Absolute velocity is used.
• Transient-state Calculation is performed.
• Compressible fluid flow is chosen.
• Gravity is enabled in -ve Y axis • 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 model-Laminar is used.
• 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 simulation.Phase used for patching was water and volume fraction of 1 is taken for water surface and 0 for air surface.
• In Pressure velocity coupling SIMPLE scheme is used. For the momentum second order upwind scheme is used and for pressure, PRESTO is selected. For the transient formulation first order implicit scheme is chosen for better results and for stability and convergence of the residuals.

Post-processing Results
The post processing was done in CFD Fluent. Case 1 is for the meshing scheme and Case 2 is for the variation in Atwood Number For Case 1A: The transient case was run for 2000 time steps with a Time step size of 0.005. A more clearer settled fluids are observed. We can observe from the water volume fraction how the air has moved upwards while the water settles down. For Case 1B: This case was ran for 4000 time steps with a time step size of 0.0005. As the mesh size is increased we get more defined results which we can see from the below plot and it took a long time for the simulation to reach appropriate results.

Plots on Water Volume Fraction
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 (at t = 0.5 to t = 9). 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 t = 0.11 vanishes as time proceed as in t = 2.253, showing the phenomena of Rayleigh-Taylor instability. At t = 0.349 and t = 2.7725, it is clear that the falling fluid is rolled up. In our result, finally the heavier fluid has settled down completely at = 7.7 or we can say it is in the verge of complete settling. Figure 3 and 4 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. The number of time steps needed in the simulation to reach 3.6 seconds flow time was 722 and 7196 for Case 1A and 1B respectively whereas the time step size is adaptive. For ≤ 0.1, we observe a broadening of the interface as the diffusion takes place. For time, ≥ 1.5, the simulation enters the non-linear regime with the increase in amplitude of perturbation. 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.

Comparison of Variation in Atwood Number with Case 1A
With the decrease in Atwood Number from Case 1A and 1B from 0.9 to 0.3,0.4 and 0.5 for 2A, 2B and 2C respectively and when compared, we can clearly observe from Figure 7 and 8 that bubbles and spikes formation are still present at the same flow time or even after completion of simulation or at the verge of settling.The formation is because of the increase in density of secondary phase with a density of ( 2 = 536.846 3 , 2 = 427.285 3 , 2 = 332.33 3 ) for Case 2A, 2B and 2C respectively. 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. 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 figures that the larger Atwood numbers, the larger amplitude ofperturbation. It takes place due to the larger density ratio of two fluids is, the stronger force to develop the interface instability. 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. The bubble height is approximately a quadratic function of time at 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.Since in our simulation for Case 1A and 1B has Atwood Number close to 1, heavier fluid takes bubble form and thus helping in validating that our simulation is correct.

Conclusions
In the present manuscript, classical method of Finite Element analysis was adopted to study the evolution of the interface during RTI. The governing equations were Navier-Stokes equations and Cahn-Hilliard equations. The initial displacement disturbance was defined by the Cosine function. Steady state solution determines the fully developed solution. This numerical simulation was validated with the observations obtained by previous researchers and this was done by the variation in Atwood Number where we obtained the same observation. In lower Atwood numbers, spikes and globules were still present at the end of the simulation which was not there in Atwood Number of 0.9 in Case 1A.
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 seen from the Case 1B with a time step size of 0.0005. 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 simuation is the ideal setup to run the calculation for Rayleigh Taylor Instability.