Week 5 - Rayleigh Taylor Instability

                                                           Rayleigh Taylor Instability

 

AIM: To conduct the Rayleigh Taylor CFD simulation. 

OBJECTIVE:  

• What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, 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 (Attach them from your google drive). Observe the differences in the animation and explain the effect of mesh size on simulation. Also, explain why a steady-state approach might not be suitable for these types of simulation. [ Choose Air and Water for this simulation] 
• Run one more simulation with water and user-defined material (density = 400 kg/m3, viscosity = 0.001 kg/m-s) for refined mesh.  
• Define the Atwood Number. Find out the Atwood number for both the cases and explain how the variation in Atwood number in the above two cases affects the behavior of the instability.  

THEORY: 

Rayleigh Taylor Instability (RTI) is a phenomenon that occurs at the interface of fluids of different densities wherein the lighter fluid pushes the heavier fluid due to the density difference. It is a dynamic process wherein the system tries to reduce its combined potential energy. Initially, the system is in a state of hydrostatic equilibrium, in which the heavier fluid (ex: water) sits atop a lighter fluid (ex: air) in a constant gravitational field. As the RTI develops, the initial perturbations progress from a linear growth phase, eventually developing "plumes" flowing upwards (in the gravitational buoyancy sense) and "spikes" falling downwards. In the linear phase, equations can be linearized, and the amplitude of perturbations is growing exponentially with time. In the non-linear phase, the perturbation amplitude is too large for the non-linear terms to be neglected. 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). Examples include supernova, water suspended atop oil, mushroom clouds formed due to volcanic eruptions. 

Practical CFD models 

SPH method (Smoother particle hydrodynamics) 

Smoothed-particle hydrodynamics (SPH) is a computational method used for simulating the mechanics of continuum media, such as solid mechanics and fluid flows. Smoothed-particle hydrodynamics is being increasingly used to model fluid motion as well. This is due to several benefits over traditional grid-based techniques. First, SPH guarantees conservation of mass without extra computation since the particles themselves represent mass. Second, SPH computes pressure from weighted contributions of neighboring particles rather than by solving linear systems of equations. Finally, unlike grid-based techniques, which must track fluid boundaries, SPH creates a free surface for two-phase interacting fluids directly since the particles represent the denser fluid and space represents the lighter fluid. 
 

SINGLE-FLUID MODEL 

A typical approach used for the analysis of two-phase flows is a mixture model, i.e., the individual fluid phases are assumed to behave as a flowing mixture described in terms of the mixture properties. The applied single-fluid model is a five-equation model consisting of the mass, momentum, and energy equations for a vapor/liquid mixture, and two equations describing the formation and growth of the liquid phase. 

TWO-FLUID MODEL 

In the two-fluid model, separate sets of the governing equation for the vapor and liquid phases have been used. The interaction between the droplets and the heat exchange between the liquid phase and the solid boundary is not modeled here as well. Additionally, the velocity slip between the vapor and the liquid phase is in this model considered. 

Turbulence model 

Turbulence models are needed to predict the average mixing behavior inflows that are on average one- or two-dimensional. The approach to the construction of the tile turbulence model is guided by experimental behavior. The equations governing turbulent flows can only be solved directly for simple cases of flow. For most real-life turbulent flows, CFD simulations use turbulent models to predict the evolution of turbulence. These turbulence models are simplified constitutive equations that predict the statistical evolution of turbulent flows. 

  

SOLUTION AND RESULTS 

In this project, we consider 3 different cases in which 2 cases are based on the type of MESHING and the other is with user-defined materials. 

We consider 2 materials air and water placing one above the other and see the effect of gravity acting on them. 

GEOMETRY 

For our project, we consider a simple 2D geometry of which 

• We draw 2 squares one above the other with a length of 0.02 m. 
• Then naming the below square as air and above square as water in 2 different components. 
• And share the topology for a perfect intersection between both surfaces. 

Creating first square

 creating second square over the first square.

After creating the both squares we should do to share the topology for a perfect intersection between both surfaces.

 

BASE MESH

MESHING  

There is not much change to the fundamental mesh as it is a 2D object  

Overall Element size: 0.5mm 

 

 

SETUP 

Solver Type	Pressure-Based
Velocity Formulation	Absolute
Time	Un-Steady
Models	Viscous model: laminar  Multiphase: Volume of fluid                  : Implicit formulation                   : 2 phases
Material	Fluid type: Air                 : water liquid   Solid type: Aluminum
Cell zones	Fluid type: Air                  : water
Boundaries	Standard
Reference values	Standard
Initialization	Standard  Patch: Water-surface (1 value of volume fraction of water)           : Air-surface (0 value of volume fraction of water)

 

Set the time should be transient and the gravity should be -9.81 because y-axis shows upward direction if you sign is negative then is convert the downward direction.

viscous model should be laminar

select the multiphase option and On the voume of fluid(VOF)

creating the Phases as Air and Water

 After Initialization use the option of patch its represent the which phase as water and air. we should give the value of water is 1 and air is 0.

 

Craeting the contour of phases

creating as Animation of contour

Results 

 

Residual plot 

 

Change in the Counters phase-vise 

 

 

 

 

 Video of Baseline Mesh(0.5mm) RTI

 https://drive.google.com/file/d/1kkptjFp8HGtMEEh3NjSTvL3fPt4uhPZk/view?usp=sharing

 

CASE 1: 

MESHING  

Overall Element size: 0.3mm 

 

Residual plot 

 

Change in the Counters phase-vise 

 

 

 

 

 

Video of Mesh(0.3mm) RTI

 https://drive.google.com/file/d/1KmeTLNVbFAvLcFI91c_ZPT-YJhs36BTz/view?usp=sharing

 

 

CASE 2

MESHING

Overall Element size: 0.1mm 

 

Residual plot 

 

Change in the Counters phase-vise 

 

 

 

Video of Mesh(0.1mm) RTI 

https://drive.google.com/file/d/1XpuiLLLr1XRqdpIzPB7EDbx4Tq4HNEgJ/view?usp=sharing

 

CASE 3

MESHING

Overall Element size: 0.1mm 

 

SETUP 

Solver Type	Pressure-Based
Velocity Formulation	Absolute
Time	Un-Steady
Models	Viscous model: laminar  Multiphase: Volume of fluid                  : Implicit formulation                   : 2 phases
Material	Fluid type: User-defined-fluid (Density =400 kg/m3, viscosity = 0.001 kg/m-s)                 : water liquid   Solid type: Aluminum
Cell zones	Fluid type: Air                  : water
Boundaries	Standard
Reference values	Standard
Initialization	Standard  Patch: Water-surface (1 value of volume fraction of water)           : Air-surface (0 value of volume fraction of water)

  

Results 

Residual plot 

 

Change in the Counters phase-vise 

 

 

Video of Mesh(0.1mm) RTI 

 https://drive.google.com/file/d/1XpuiLLLr1XRqdpIzPB7EDbx4Tq4HNEgJ/view?usp=sharing

  

Observation: 

In the above cases, the simulation starts with the state of hydrostatic equilibrium and the Rayleigh Taylor instability is observed at the interference when a lower density fluid pushes a higher density fluid due to which formation of shock waves at the interface takes place. Formation of air bubbles starts taking place, which compresses the heavy fluid around it, due to which shock waves of multidimensional fashion generates and it gets more stronger as it moves upward.   

It is observed that the more we refine the mesh, the more the simulation results get smoother about the irregularities that take place at the interface of the two fluids. As higher density fluid replaces lower density fluid, the formation of air bubbles and vortex takes place, which travels towards the upward region with time. In case 3, we can see the more detail results about the formation of waves and the formation of some air bubbles that get trapped at the lower region during the initial stage and then they travel towards the upper region, generating shock waves. At the end of the simulation, it is to be observed that the two phases get separated from each other though, we can see some diffusivity between them which is a volume fraction of air and water. 

  

why a steady-state approach might not be suitable for the above types of simulation? 

The difference between the steady and transient is that you can't see the small-time variation of instability. The steady-state simulation is performed if we are concerned more about the final state results or the equilibrium state. In RT-Instability CFD models, we are more concerned to learn about the transition of the irregularities that starts developing when we pour high dense fluid upon low dense fluid under gravity effect so, by using transient solver along with refined mesh of the model, we can compute the smooth transition of irregularities that takes place at the interface of the fluids. The final state results for both the steady-state and transient state will be the same. 

In this problem we are observing the instabilities when it is occurring so we are not concerned about the final answer because steady state is more into capturing the final results but by transient, we can see the behavior of the solution and every instant such that capturing bubbles, vortex and shockwaves. so, this is reason why transient is more suitable than steady state model 

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=(rho_1-rho_2)/(rho_1+rho_2)`

`rho_1` is the density od heavier fluid`

`rho_2` is the density of lighter fluid`

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 asymmetric “finger” of fluid; for Atwood number close to 1, the much lighter fluid “below” the heavier fluid takes the form of larger bubble-like plumes. 

The calculated Atwood number is close to 1 and from the simulation results, it is found that when high dense fluid i.e. water poured upon low dense fluid i.e. air under gravity, the formation of air bubble-like plumes takes place which travels towards upward region in the form of waves and some gets trapped at the lower regions during the initial stages, which afterward try to move towards upper region and at the end, both phases gets separated with some diffusivity left at the middle portion between them. Thus, the calculated Atwood number is validated for our simulation results.   

For above cases 

case 1: 

heavier fluid(water) 

density 1=998.2 `(kg)/m^3` 

lighter fluid(air) 

density 2=1.225 `(kg)/m^3`

Atwood number(A) = (density 1-density2)/ (density 1+density2) 

Atwood number(A)=0.996 

case 2: 

heavier fluid(water) 

density 1=998.2 `(kg)/m^3`

lighter fluid (user-defined material) 

density 2=400 `(kg)/m^3`

Atwood number(A) = (density 1-density2)/ (density 1+density2) 

Atwood number(A)=0.427 

variation in Atwood number in the above two cases affects the behavior of the instability.  

For Atwood number close to 0, RT instability flows take the form of symmetric fingers of fluid. 

When Atwood number close to 1, the much lighter fluid below the heavier fluid takes the form of larger bubble-like plumes. 

Therefore, from the value of Atwood number it is seen from the simulations that air forms larger bubbles when water is poured down the air pushes through and forms large bubbles. Therefore, the results can be validated by obtaining the value of Atwood number. 

The behavior of RT Instability on the variation of Atwood number: 

RTI evolves in three distinct stages: Linear stability, mushroom head(falling) or bubble(rising) formation, long term evolution due to bubble merging and mixing. For low Atwood number, bubble and mushroom head shape is relatively symmetrical, the disturbances or perturbations are linear for a larger time. In contrast, for large Atwood number, bubble and mushroom head shape is less symmetrical, also it occurs very rapidly at an accelerating rate which leads to a nonlinear growth rate. 

CONCLUSION: 

• From the above cases, we observe that finer mesh captures the instability better than the coarse mesh, but it takes more time for the heavier fluid to settle down. As computational time is long but produces results with better accuracy. 
• Since the simulation involves capturing instabilities that are changing in an unpredictable manner as a function of time, the transient approach is preferred over the steady-state approach.  
• The Atwood no validates the CFD simulation as it is observed. In the case of water and air mixing, bubble formation takes place, whereas in the other case no bubble takes place.