Week 5 - Rayleigh Taylor Instability

Simulation and Analysis of Rayleigh Taylor Instability

Objective-

To perform Rayleigh Taylor Instability simulation using a CFD model of two different fluids in a system.

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.

RTI occurs at very large scales; e.g., interstellar gas, pushed out of the galactic plane by magnetic fields and cosmic rays (Zweibel 1991). It also occurs at very small scales; e.g., RTI limits the driving pressure in the collapse of acoustically forced bubbles in sonoluminescence phenomena (Lin et al. 2002).

RTI has a profound effect on Earth's climate. Winds off the coasts of Greenland and Iceland drive evaporation of the ocean surface making the upper layers saltier and heavier. This generates descending plumes that drive the abyssal ocean circulation, which in turn acts as a heat pump transporting warm equatorial waters up north. Without ocean overturning from RTI, northern Europe would likely face drastic drops in temperature (Calvin 1998).

RTI plays a crucial role in all known forms of fusion, whether confinement be magnetic (Bateman 1979), inertial (Petrasso 1994) or gravitational (Burrows 2000).

CFD models based on mathematical analysis of Rayleigh Taylor waves:

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

2. Single fluid model:

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

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

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

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

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

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

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

Procedure-

Modelling-

Here, we are making two squares of size 20X20, depicting two containers in which top square is filled with water and bottom square is filled with air and they are named accordingly after meshing. After making square we have to use share topology tool. Below image shows the model -

Meshing-

Case1- For the first case we are taking element size as 0.5mm and then naming the top and bottom surface as water and air respectively. Below image shows the model after meshing-

Pre-Processing and Solver Setting:

Here, we are using solver type as pressure based and transient solver. Also, we have to enable gravity in –y direction. Then, select viscous model as laminar model. Here, we have used simple method for solving the equation.

Then for adding water as material, go to ‘materials’ and copy ‘water-liquid’ from fluent database.

Then, under multiphase model select ‘volume of fluid’ and implicit formulation for solution since in explicit formulation time step size is constraint, i.e., time step size should satisfy the Courant number (non-dimensional number that tells how fast information travels through our discretized domain)

Then under phases tab, select air and water as primary and secondary phase and rename accordingly. Below are the images which shows the tab of multiphase model-

Now, initialize using hybrid initialization. Then, we have to use patch tool to patch different values of flow variable into different cells. For this, select water under phase > volume fraction > water under zones to patch and give 1 value to it and select patch. Similarly give 0 value to air and select patch. Below image shows the patch tab-

Now to make contour, select phases under contour and select water as phase and select save/display. Below image shows the tab of contours depicting the selection-

Now, run the calculation using time step size = 0.005 and number of time steps = 1000.

Case 2-

Here, mesh size is 0.1mm and all other processes are same as the previous case. Below image shows the model after meshing.

After meshing and setting up the solver settings, we ran the calculation using time step size = 0.001 and number of time steps =1300.

Case 3-

In this case, we are using the same model and same mesh size of the elements and the only thing we are changing here is the fluid air at the bottom. We are replacing the air at the bottom with a user defined fluid material with density = 400 kg/m³ and viscosity = 0.001 kg/m-s.

For changing the fluid, we can just change the name and input density and viscosity values of air in the material tab.

Result-

The contour of water volume fraction of cases 1, 2 and 3 are shown below-

The plot of scaled residuals of cases 1, 2 and 3 are shown below-

The animation of the volume fraction of cases 1, 2 and 3 are shown below-

Comparison between cases 1 and 2- From the above animation you can easily see the details captured in the flow of fluids in case 2 which is lacking in the case 1. It is due to the coarse and fine mesh size of the elements in cases 1 and 2. The fine mesh captures the details of the flow of the fluid.

Here, we have used transient model instead of steady model for solving because we want to quantify the parameters for each and every second, i.e., time dependent parameters.

Atwood Number

It is a dimensionless number in fluid dynamics that is used to study the hydrodynamic instabilities in density-stratified flows. The Atwood number equals the ratio of the difference in mass densities to their sum as given below-

A = (ρ₁ – ρ₂) / (ρ1 + ρ₂)

where ρ₁ is density of heavier fluid

ρ₂ is density of lighter fluid.

For Air-water Rayleigh Taylor instability,

ρ₁ = 998.2 kg/m³

ρ₂ = 1.225 kg/m³

Thus, A = (998.2 – 1.225) / (998.2+1.225) = 0.997

For water- user defined material Rayleigh Taylor instability,

ρ₁ = 998.2 kg/m³

ρ₂ = 400 kg/m³

Thus, A = (998.2 – 400) / (998.2+400) = 0.427

Behavior of Rayleigh Taylor Instability on variation of Atwood number:

Rayleigh Taylor Instability 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-

For better capturing instability or better capturing details of flow of fluid in a system, we should use finer mesh size in the CFD model.

Also, the variation of Atwood number has direct impact on the perturbation amplitude and evolution of interface growth.