Rayleigh Taylor Instability of water with air and a user-defined-material.

OBJECTIVES:

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

2) 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.

3) Run one more simulation with water and user-defined material(density = 400 kg/m3, viscosity = 0.001 kg/m-s) for refined mesh.

4)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 densities where 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.

From Newton's law of viscosity, we know that even the fluid cannot resist the slightest shear force. Density is one of the main properties of the fluid which the shear force can incorporate with it. Depending upon the density we have lighter fluid like water etc and denser fluid like mercury, oil, honey, etc.

Lord Rayleigh observed that when two liquids i.e lighter and heavier fluids are mixed together where the denser fluid is top on the lighter fluid, due to the gravity and the density variation the interface of the liquid is not stable. (i.e) water suspended on the oil.

So, the instability of the interface between the different densities of fluid when the lighter fluid is accelerated towards the heavier fluid is knowns as RAYLEIGH TAYLOR INSTABILITY.

Types of instabilities

Ritchmyer-Meshkov instability: The Richtmyer-Meshkov instability arises when a shock wave interacts with an interface separating two different fluids. It combines compressible phenomena, such as shock interaction and refraction, with hydrodynamic instability, including nonlinear growth and subsequent transition to turbulence, across a wide range of Mach numbers.

Plateau Rayleigh instability: The Plateau Rayleigh instability, often just called the Rayleigh instability, explains why and how a falling stream of fluid breaks up into smaller packets with the same volume but less surface area. This fluid instability is exploited in the design of a particular type of inkjet technology whereby a jet of liquid is perturbed into a steady stream of droplets.

Kelvin instability: can occur when there is velocity shear in a single continuous fluid, or where there is a velocity difference across the interface between two fluids. An example is a wind blowing over water: The instability manifests in waves on the water surface. More generally, clouds, the ocean, Saturn's bands, and the sun's corona show this instability.

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 supernovas, water suspended atop oil, mushroom clouds formed due to volcanic eruptions.

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

MODELING APPROACH:

• Choose Fluent from the workbench setup, edit SpaceClaim to create geometry.
• Choose the XY plane to sketch 2Dgeometry, create a square of 20mm, and select pull tool to make it a surface of that.
• Move the geometry to a new component and rename it as Air.
• Now create new geometry under the main design as water and create a square of 20mm right on top of the air square which creates an interface between the 2 boxes.
• Select the entire geometry for share prep under the workbench and select share so that the interface can face information and conformal meshing can be achieved between the 2 phases.

As the simulation is 2D and geometry is very simple we can adapt the conformal mesh with hexahedral element type. To make sure conformal mesh is present for simulation it is important to have geometry topology as shared. In which the common edge between two surfaces is shared with each other. This can be done by using the share topology option in space-claim which will create a common interface between the two geometry.

Fig 1: Geometry

MESHING

The simulation is run for three cases with mesh sizes of 0.5mm, 0.3mm and 0.2mm. Hexahedral mesh is applied by default.

Fig 2: (a)mesh size=0.5 (b)mesh size=0.3 (c)mesh size=0.2

SOLVER SETTINGS

Type-pressure-based

Velocity formulation- Absolute

Time- Transient

2D space-planer

Gravity=`-9.81m/s^2` in y-direction

The setting of physics Viscous model-Laminar Multi phase model-Volume of fluid Courant number-0.25s

Solution Initialization Standard time step size-0.0005 No of time step-1500

Patch: Water-surface (1 value of volume fraction of water)

: Air-surface(0 value of volume fraction of water)

Residuals

Case (a) mesh size 0.5(air+water) Case (b) mesh size 0.3(air+water)

Case (c) mesh size 0.2(air+water) Case (d) mesh size 0.2(use defined material+water)

observation:

Case (a)

Case (b)

Case (d) mesh size 0.2(use defined material+water)

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. It is observed that by refining the mesh, the more the simulation results get smoother about the irregularities that take place at the interface of the two fluids and capture more instability. 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 more detailed 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 travel towards the upper region. 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.

Steady vs Transient

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 the 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. So this is the reason why transient is more suitable than the 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 = (1000 - 1.25) / (1000 + 1.25) = 0.9975`

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 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 the 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 the 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

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.

• In the case of water and air mixing, bubble formation takes place, whereas in the other case no bubble takes place.

• Therefore, the results can be validated by obtaining the value of the Atwood number. Atwood number is an important parameter in the body of Rayleigh Taylor instability and Richtmyer-Meshkov instability. In Rayleigh-Taylor instability, the penetration distance of heavy fluid bubbles into the lighter fluid is a function of the acceleration time scale. `A*g*t^2`, where 'g' is the acceleration due to gravity and 't' is the time.