Week 5 - Rayleigh Taylor Instability
Objective -
- 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/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 affect the behavior of the instability.
Introduction to Rayleigh Taylor Instability -
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, Water suspended atop air is an everyday example of Rayleigh-Taylor instability, and it may be modelled by two completely plane-parallel layers of immiscible fluid, the denser fluid being on the top and both subject to the Earth's gravity.
The equilibrium here is unstable to any perturbations or disturbances of the interface: if a parcel of heavier fluid is displaced downward with an equal volume of lighter fluid displaced upwards, the potential energy of the configuration is lower than the initial state. Thus, the disturbance will grow and lead to a further release of potential energy, as the denser material moves down under the (effective) gravitational field, and the less dense material is further displaced upwards.
As the RT instability develops, the initial perturbations progress from a linear growth phase into a non-linear growth phase, eventually developing "plumes" flowing upwards (in the gravitational buoyancy sense) and "spikes" falling downwards. In the linear phase, the fluid movement can be closely approximated by linear equations, and the amplitude of perturbations is growing exponentially with time. In the non-linear phase, perturbation amplitude is too large for a linear approximation, and non-linear equations are required to describe fluid motions.
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). The difference in the fluid densities divided by their sum is defined as the Atwood number, A. For A close to 0, RT instability flows take the form of symmetric "fingers" of fluid; for A close to 1, the much lighter fluid "below" the heavier fluid takes the form of larger bubble-like plumes.
Types of Instability and some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves -
Richtmyer-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 describes the phenomenon where a falling jet or cylinder of fluid (such as water) at one-point ceases to be a jet and breaks into multiple droplets of smaller total surface area. This fluid instability is exploited in the design of a particular type of ink jet technology whereby a jet of liquid is perturbed into a steady stream of droplets.
Kelvin-Helmholtz instability - Kelvin-Helmholtz (KH) instability typically occurs between two immiscible fluids with a velocity difference across their interface. This instability can be observed, for example, as a natural phenomenon of certain cloud patterns. These patterns are characterized by the generation of small-scale motion at the interface, which later leads to the formation of an unstable interfacial vortex street. This generates spiral structures at the fluid-fluid interface. Wave formation on still water due to blowing wind is another common example of this instability.
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 empty space represents the lighter fluid.
Single Fluid Motion -
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 Motion -
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 are not modelled here as well. Additionally, the velocity slip between vapor and the liquid phase is in this model taken into account.
Turbulence model -
Turbulence models are needed to predict the average mixing behavior in flows that are on average one- or two-dimensional. The approach to the construction of tile turbulence model is guided by the experimental behavior. The equations governing turbulent flows can only be solved directly for simple cases of flow. For most real-life turbulent flows, CFD simulation use turbulent models to predict the evolution of turbulence. These turbulence models are simplified constitutive equations that predict the statistical evolution of turbulent flows.
Rayleigh—Taylor Instability Simulation -
The idea is, the dense fluid will be at the top and lighter fluid will be at the bottom and due to gravity, there will be instability at their interface, which is going to be calculated by using VOF multiphase method in ANSYS Fluent.
CASE 1 -
Geometry -
The geometry is essentially two squares, one on top of another of side 20 mm. The top square has been named 'water' and bottom square 'air'. The common edge (coloured blue) between the domains was shared using share tool in Workbench. The created geometry is depicted below

Meshing -
For case 1 we'll be using the base line mesh which has an element size of 0.5mm as shown below.


Setup -
For this simulation we'll be using a pressure based, transient solver with gravity (9.81m/s2) turned on in the negative y direction as shown below

The Method we'll be using is the Laminar Volume of Fluid (Vof) multiphase method with implicit scheme as shown below

Material Properties -


Since we're using the multiphase method, we'll have to associate a material to a phase as shown below

After Standard Initialization ,we have to patch both the domains to the respective materials, wherein we provide a volume-fraction value of 1 to water and 0 to air as shown below


In order to see if we have done the patching correctly, we'll create a phase contour plot and it should show water at top and air at bottom at our geometry.

We'll now run the simulation using the implicit scheme with time step - 0.005 and the total steps - 1500.
Results -
Residuals -

Solution Animation -
CASE 2 -
The only change we need to make in the case 2 is that we'll refine the mesh by changing the element size to 0.2mm from 0.5mm in the previous case. Except that everything else stays the same.
Refined mesh -


Now we can start the simulation but before that in the setup we'll have to do the patching again as we've changed the mesh the previous patching is no more.
Results -
Residuals -

Simulation Animation -
CASE 3 -
In this case we'll replace air with a user defined liquid whose density = 400 kg/m3 and viscosity = 0.001 kg/m-s. All the things are same as in the above cases but we'll replace air with the user defined liquid.
Now we can run the simulation.
Results -
Residuals -
Simulation Animation -
Observations -
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 the 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, this is due to the fact that we have a greater number of nodes and due to which we can see the changes happening at the smallest level and hence mesh refinement produces better result.
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 reason we are using the transient state approach and not the steady state approach is because we need to observe the small-time variation of instability which can only be done by transient state approach. The steady-state simulation is performed if we are concerned more about the final state results or the equilibrium state.
In RT-Instability analysis, we are more concerned 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 as it is occurring and are not concerned about the final answer because steady state is more into capturing the final results but by transient state, we can see the behavior of the solution at every instant such and hence see properly the formation of 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 = (ρ1 - ρ2)/(ρ1 + ρ2)`
where
ρ1 = density of heavier fluid
ρ2 = density of lighter fluid
Atwood for case 1 and 2 is calculated to be:
`A = (ρ1 - ρ2)/(ρ1 + ρ2) = (998.2 - 1.225)/(998.2 + 1.225) = 0.9975`
Atwood for case 3 is calculated to be:
`A = (ρ1 - ρ2)/(ρ1 + ρ2) = (998.2 - 400)/(998.2 + 400) = 0.4278`
For A close to 0, RT instability flows take the form of symmetric "fingers" of fluid; for A close to 1, the much lighter fluid "below" the heavier fluid takes the form of larger bubble-like plumes.
Variation in Atwood number in the above two cases affect 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 is 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.
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 light fluid is a function of acceleration time scale. A*g*t2, where is the acceleration due to gravity and 't' is the time.