Week 5 - Rayleigh Taylor Instability

Aim - To conduct the Rayleigh Taylor CFD simulation.

Objectives:

• some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves
• 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
• explain why a steady-state approach might not be suitable for these types of simulation.
• Run one more simulation with water and user-defined material (density = 400 kg/m³, 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.

Introduction - Rayleigh Taylor Instability (RTI) is a phenomenon that occurs at the interface of fluids of different density 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.

Theory - 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.

CFD models based on mathematical analysis of Rayleigh Taylor waves:

A. Mass diffusion model:

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

B. Single fluid model:

a. Potential flow model

This model is valid for all instability changes from the linear stage through the early non-linear to the asymptotic but limited to Atwood number = 1. In this model, single-mode perturbation expands the flow equations to second-order around the bubble tips, hence obtaining an ordinary differential equation for the bubble velocity.

b. Buoyancy drag model

It is valid for every Atwood number but limited to the non-linear 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 2 or 3 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.

c. Besnard-Harlow-Rauenzahn(BHR)

It is a more advanced version of the single-fluid model. The BHR model is based upon 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 Reynold stresses, turbulent mass-flux velocity, density fluctuations, and the turbulent kinetic energy dissipation rate.

C. 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 Reynold averaging of the governing equations and resultant terms.

a. Two equation K-L turbulence model:

The K-L turbulence model involves equations for turbulent kinetic energy K and turbulent length scale L. The starting point for deriving the model equations are the buoyancy-drag models for the self-similar growth of RT instabilities.

b. Three-equation K-L-a turbulence model:

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

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

The four equation model was developed as an extension to the three-equation K-L model by including a fourth equation for density-specific volume covariance. The model can accurately reflect the effects of the changes in the density fluctuations.

Solution procedure :

• We wil simulate both fluid flow and also the heat transfer to the solid i.e., CHT Analysis on an Exhaust port with the steady state .Then for solving we follow the below procedure :

A)Geometry setup

B)Meshing

C) Problem setup

D)Solving

E)Post processing

Simulating the flow with the transient solver

Case-1 : Rayleigh–Taylor instability between air and water with coarse mesh of size 0.5 mm :

A) Geometry setup :

• In fluent module we have to open spaceClaim .
• At first, geometry is created as square boxes we are going to create square of dimensions 20 mm x 20 mm by using rectangle in spaceClaim and then convert it to surface by using Pull tool.After creating surface move it to a new component and name it as air .
• Then right click on the main design and create a new component and name it as water.Then again create a new square sketch of dimensions 20 mm x 20 mm by using rectangle over the previous created on on the same axis and next thing is as we need conformal mesh we need to use share in workbench and select both componens and click ok.Then you see share topology is also automatically set to share and you will see right icon modified and blue line is created at interface between them.
• The geometry was created in spaceclaim and it involves no other importing of CAD models. The surface model was further processed for meshing .This preprocessing phase of the workflow is extremely important because it determines the simulation potential and limits the CFD results.
• This completes the creation of geometry.Then we close the spaceClaim ,the progress will be saved automatically.

B) Meshing :

• For meshing we open Meshing,where we mesh our model so that we can use it to solve our model.
• In general it is always good to name the surfaces before meshing. So let us do that.
• We first select the entity by using the entity selection by using faces so that we can provide name for our entity selection.

steps : select the face you want to name by using face selection -> then press N -> name the selection.

• In the similar way we name the water side and air side.
• Then next thing to do is to mesh.We click the mesh option in the side then give element size as 0.5 mm and mesh.So this is our first case .

• Now we have to check the quality of the mesh.

steps : mesh ->quality-> mesh metric-> Element quality

• Here our minimum quality is around 15 %,so we can move further for our solution setup.

C) Problem Setup :

• To set up our problem we open Fluent and setup inital conditions , boundary conditions material etc.,
• We first start with checking the mesh quality using perform mesh check under domain.If no error is seen in console window then we can proceed further.
• Then we move move to setting up physics where we select the solver,viscous model etc.
• We use a transient state pressure based solver to solve our model.

Physics :

General :

Solver type : pressure based

time : transient state

velocity formulation : absolute

• We need to enable gravity by giving value for Y as -9.8 m/s^2 as gravity is acting in negative Y-direction.

models :

viscous model : Laminar (since the velocity is very small)

Material : Fluid -1 : Air

Fluid -2: Water

In models enable multiphase flow and in that select,

for model : select VOF method

Number of Eulerian phases : 2

Formulation : Implicit

• And in phases option as we selected two phases give names and selcect phase ,
• For primary phase :change name as air and change phase material to air
• For secondary phase :change name as water and change phase material to water.
• Then for doing patching we have to initialize first and then we do Hybrid initialization .
• In patch for phase select water and for variable you can see volume fraction now .
• Now, patch the value fraction value of 0 to air side and 1 to water side.(i.e., at t=0 the top box will be total water and bottom box will be total air.)
• Now we are almost ready to run the simulation but as to verify the patching,we create contour for volume fraction of water .
• Now we want to see the volume fraction of water contour.
• For contour -> go to solutions tab -> graphics -> contour -> in contour of : select phases and under that select volume fraction and for phase select water -> select both the boxes or surfaces and press save /display.
• Then you can see the initial values for volume fraction of water on whole domain,where you can see in the bottom it is set to 0 and for the top region it is set to 1.So we can conclude that our initialization is correct.
• Then we want to create animation so that we can see how the solution or volume fraction of water is is varying when problem is being solved.
• For that,

steps : Activities -> solution animation -> select volume fraction contour,which we have created -> preview -> set the required view and use active -> and name it as volume fraction of water animation.

• Then in solution method for 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 choosen for better results and for stability and convergence of the residuals.
• Then we are ready to solve .

D) Solving :

• Give iterations and time step such that solution is converged.
• Run the simulation until the steady state is reached.Lets run for 1300 times steps with time step size = 0.005 and the total run time =6.5 s .
• we can find the solution is converged.
• If you see error of diverging solution ,then try to under relax or try to use lower order discretization schemes so that we can ge converged solution.
• Then let us post process the results.

E)Post processing :

Residuals plot :

You can observe from the below plot the residuals has converged.

Volume fraction of water contour :

You can observe from the below contour the volume fraction of water.

Case- 2 :Rayleigh–Taylor instability between air and water with fine mesh of size 0.1 mm :

• Everything is same from the point of geometry and when it comes to mesh the only change is the mesh size,here it is quite refined with mesh size of 0.1 mm.
• Comming to set up even it is also same.
• Here we run for 1400 times steps with time step size = 0.005 and the total run time = 7s.
• So let us directly run the simulation then postprocess the results and compare in between the results of two different grid sizes .

Residuals plot :

You can observe from the below plot the residuals has converged .

Volume fraction of water contour :

You can observe from the below contour the volume fraction of water.

• You can observe the volume fraction of the same problem but the mesh is even more refined(0.07 mm).Here the numerical diffusion(due to low order scheme used for convergence)is even more less ,we get some more accurate result than above.So small the grid size more the computational effort and more the accurate solution the result is.

Case- 3 :Rayleigh–Taylor instability between user defined material and water with fine mesh of size 0.1 mm :

• Everything is same from the point of geometry and even when it comes to mesh, the mesh size is also same as refine case above i.e., 0.1 mm.
• Comming to setup the change comes here ,we change the lighter material i.e.,from air to user defined material with density = 400 kg/m³, viscosity = 0.001 kg/m-s .
• Then again every thing is same ,we patch the lighter fluid down and denser fluid upper side of the box or domain.
• Then we run for 1800 times steps with time step size = 0.005 and the total run time = 9.00s.
• So let us directly run the simulation then postprocess the results and compare the results with above case -2 with respect to Atwood number .

Residuals plot :

You can observe from the below plot the residuals has converged .

Volume fraction of water contour :

You can observe from the below contour the volume fraction of 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 probleem we are observing the instabilities when it is occuring so we are not concerened about the final answer because steady state is more into caapturing the final results but by transient we can see the behaviour of the solution and every instant such that capturing bubbles,vortec and shockwaves. so this is reason why transient is more suitable than steady state model.

Atwood Number:

It is a dimensionless number that is used to study hydrodynamic instabilities in density stratified flows.

A=ρ1−ρ2ρ1+ρ2

Where,

ρ1- is density of heavier fluid

ρ2-is density of lighter fluid

For Air-water RT instability,

ρ1 = 998.2 kgm3 [Water density]

ρ2= 1.225 kgm3 [Air density]

A = (998.2-1.225)/(998.2+1.225)

A = 0.998

For water- UD material RT instability,

ρ1 = 998.2 kgm3 [Water density]

ρ2 = 1.225 kgm3 [UD material density]

A = (998.2-400)/(998.2+400)

A = 0.428

Behaviour of RT Instability on 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 non linear growth rate.

Conclusions :
In the present challenge, classical method of FVM was adopted to study the evolution of the interface during RTI. The governing equations were Navier-Stokes equations and Cahn-Hilliard equations.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.99 in Case -2 . From the grid dependency we also observed that the interface thickness reduces with the refined mesh.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 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 simuation is the ideal setup to run the calculation for Rayleigh Taylor Instability.