Rayleigh Taylor Instabiliy Challenge

AIM:

1. 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]
2. Run one more simulation with water and user-defined material(density = 400 kg/m³, viscosity = 0.001 kg/m-s) for refined mesh.

INTRODUCTION:

The purpose of this problem is to simulate a two-dimensional Rayleigh-Taylor instability problem using the classical method of Finite Element analysis of a multiphase model using ANSYS FLUENT. The governing equations consist of a system of coupled nonlinear partial differential equations for conservation of mass, momentum and phase transport equations. The study focuses on the transient state simulation of Rayleigh Taylor waves and subsequent turbulent mixing in the two phases incorporated in the model. 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 in a gravitational field. The problem was governed by the Navier-Stokes and Cahn-Hilliard equations in a primitive variable formulation. The Cahn- Hilliard equations were used to capture the interface between two fluids systems.

Q: 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.

A:

Kevin-Helmholtz Instability:

The Kelvin–Helmholtz instability typically occurs when there is velocity shear in a single continuous fluid, or additionally where there is a velocity difference across the interface between two fluids. A common example is seen with wind blowing over water, the instability constant is able to manifest itself through waves on a water surface. The Kelvin-Helmholtz instability is not only restricted to a water surface as clouds, but is evident through other natural phenomena as the ocean, Saturn's bands, Jupiter's Red Spot, and the sun's corona.

In the case of a short wavelength, if surface tension is ignored, two fluids in parallel motion with different velocities and densities yield an interface that is unstable for all speeds. However, surface tension is able to stabilize the short wavelength instability and predict stability until a velocity threshold is reached. The linear stability theory, with surface tension included, broadly predicts the onset of wave formation as well as the transition to turbulence in the important case of wind over water.

The following illustration shows the flow pattern for two-phase medium:

Fig. 1 Kevin-Helmholtz Instability

In the case of a short wavelength, if surface tension is ignored, two fluids in parallel motion with different velocities and densities yield an interface that is unstable for all speeds. However, surface tension is able to stabilize the short wavelength instability and predict stability until a velocity threshold is reached. The linear stability theory, with surface tension included, broadly predicts the onset of wave formation as well as the transition to turbulence in the important case of wind over water.

It was recently discovered that the fluid equations governing the linear dynamics of the system admit a parity time symmetry, and the Kelvin-Helmholtz instability occurs when and only when the parity-time symmetry breaks spontaneously.

For a continuously varying distribution of density and velocity (with the lighter layers uppermost, so that the fluid is RT-stable), the dynamics of the Kelvin-Helmholtz instability is described by the Taylor Goldstein Equation

`(U-c)[psi -k^{2}psi ]+ (N^{2}/{U-c})-U]psi =0`

and its onset is given by the Richardson Number `Ri`. Typically the layer is unstable for `Ri<0.25`. These effects are common in cloud layers. The study of this instability is applicable in plasma physics, for example in inertial confinement fusion and the plasme-beriyllium interface. At times a situation in which there in a state of static stability, evident by heavier fluids found below than the lower fluid, the Rayleigh-Taylor instability can be ignored as the Kelvin-Helmholtz instability is sufficient given the conditions.

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. It is related to the RTI and is part of a greater branch of fluid dynamics concerned with fluid thread breakup. 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.

The driving force of the Plateau–Rayleigh instability is that liquids, by virtue of their surface tensions, tend to minimize their surface area. A considerable amount of work has been done recently on the final pinching profile by attacking it with self-similar solutions.

The following figure shows the diagram for Plateau-Rayleigh instability:

Fig. 2 Plateau-Rayleigh Instability

The flow begins with tiny perturbations in the stream. Some grow and some decay over time. The growth of some may be faster and that of others might be slower. The growth rate of the perturbations can be mapped as a function of wave number and radius of the cylindrical system. Assuming initially that all component amplitudes are roughly equal, the wave numbers are checked. As time progresses the component with maximum growth rate will be the one that pinches the stream into drops. Young’s Laplace equation uses the concept of Cylinder radius and Radius of Curvature (ROC) to assess the pressure applied to the stream. As the cylinder radius decreases (at trough), pressure due to surface tension increases, and as the Radius of Curvature is negative at the trough, the pressure decreases. This might lead to a small cancellation of the pressure effects due to both phenomenons. If ROC dominates the Cylinder Radius, the component results in a decay and vice-versa. The highest growth rate component satisfies the following relation to pinch the droplets:

`kR_0~=0.697` where `k` is the wave number and `R_0` is the initial radius.

Richtmyer–Meshkov instability

The Richtmyer–Meshkov instability (RMI) occurs when two fluids of different density are impulsively accelerated. Normally this is by the passage of a shock wave. The development of the instability begins with small amplitude perturbations which initially grow linearly with time. This is followed by a nonlinear regime with bubbles appearing in the case of a light fluid penetrating a heavy fluid, and with spikes appearing in the case of a heavy fluid penetrating a light fluid. A chaotic regime eventually is reached and the two fluids mix. This instability can be considered the impulsive-acceleration limit of the RTI.

The following figure shows a depiction of the instability:

Fig. 3 Richtmyer Meshkov Instability

This instability involves Shock interaction with an interface where vorticity is deposited at interface by means of baroclinic torque. The mechanism of the vorticity involves misalignment of the pressure gradient across the shock and the local density gradient at the contact during shock passage. The rotational flow associated with this localized vorticity causes the interface to distort strongly in areas of maximal misalignment, forming characteristic structures such as bubbles and spikes.The respective fluids `a` and `b` possess respective densities (`rho_(a_0),rho_(b_0)`) and specific weights (`gamma_(a),gamma_(b)`) at the initial condition. The initial density ratio :

`R_0=rho_(a_0)/rho_(b_0)`

The shock is reflected if `R_0>R_0^text(tt)`. The shock undergoes total transmission if `R_0=R_0^text(tt)`, and Reflected Rarefaction if `R_0`<`R_0^text(tt)`. where,

`R_0^text(tt)=(gamma_b(gamma_b+1)M_i^2)/(gamma_a-gamma_b+gamma_b(gamma_a+1)M_i^2)`

Initial Shock Mach= `M_i=D_i/(sqrt(gamma_bp_0/(rho_(b_0))))`

where `p_0` is the pressure in fluid b before shock, `D_i` is shock initial velocity, and the shock is induced in fluid b.

Q: 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.

A:

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_T=(rho_(a_0)-rho_(b_0))/(rho_(a_0)+rho_(b_0))=(R_0-1)/(R_0+1)`

with the acronyms as discussed before.

In Rayleigh–Taylor instability, the penetration distance of heavy fluid bubbles into the light fluid is a function of acceleration time scale, `Atext(g)t^2` where g is the gravitational acceleration and t is the time.

For Case 1:

`A_T=(rho_(water)-rho_(air))/(rho_(water)+rho_(air))=(998.2-1.225)/(998.2+1.225)=0.9975`

For Case 2:

`A_T=(rho_(water)-rho_(user-material))/(rho_(water)+rho_(user-material))=(998.2-400)/(998.2+400)=0.4278`

It is found that the initial density stratification caused by compressibility plays a stabilizing role, while the expansion-compression effect of flow plays a destabilizing role. For the case of small Atwood number, the density difference between the two sides of the interface is small, and the density distribution of the upper and lower layers is nearly symmetrical. The initial density stratification plays a dominant role, and the expansion-compression effect has little influence. With the increase in the Atwood number, the stabilization effect of initial density stratification decreases, and the instability caused by the expansion-compression effect becomes more significant. The flow structures of bubbles and spikes are quite different at medium Atwood number. The effect of compressibility on the bubble velocity is strong at large `A_T`.

For Atwood number close to 0, the flow takes the form of symmetric finger type protrusions and when closer to 1, the heavier fluid takes large bubble like plumes.

Q: Explain why a steady-state approach might not be suitable for these types of simulation?

A:

Steady State simulations are run for equillibrium flows where the final state of the flow is required. Transient state analysis is run to understand the variation in the flow at every time step. The instability perturbations and fluctuations can be captured in Transient State.

SOLVING AND MODELLING APPROACH:

The following flow chart is followed for the Solving Approach:

Fig. 4 Solving Approach Flow Chart

ANSYS involves selection of Fluent solver, input and verification of Geometry, Meshing, Setup, Solutions and Results. The Geometry is input in ANSYS Space Claim. The mesh sizing and face naming is done while Meshing. The Setup section involves Mesh Checking, Turbulence Model selection, Solver type, velocity formulation and time (steady or transient state), Output report files, Material selection, Initial and Boundary conditions. The Results section is used for the visualiztion of the output in a GUI.

PRE-PROCESSING AND SOLVER SETTINGS:

GEOMETRY:

Fig. 5 Geometry for Rayleigh Taylor Instability experiment

The above geometry was drwan in spaceclaim. The sides of each square are 20mm and the Pull tool was used to create a small extrusion for 2D simulation. The two drawings are moved into a new component and the Share Prep tool is used to Share the information along the connecting surfaces.

MESHING:

The baseline mesh of element size 0.5mm is shown below:

Fig. 6 Baseline mesh for geometry

The refined mesh is created for an element size of 0.25mm. The Basic Meshing is used. To improve the details Inflation layers can be added at the boundaries.

SETUP:

Since the problelm statement is considered for a low Mach number, a Pressure Based Solver is considered for the simulation. The Laminar model was used for the respective test case calculations. A Transient State simulation is run.The Multiphase option is selected and the Volume of Fluids Solving method is employed. The respective Materials were copied from Fluent Database and the phases were input with the materials. Standard Initialization is used for t=0 and 500 iterations were run for each case. Each phase was patched with a volume fraction of water (1 for water and 0 for air/user material).

RESULTS:

Case 1: Baseline Mesh

Case 2: Refined Mesh

Case 3: Refined Mesh- User Defined Material

CONCLUSION:

For the first test case a bubble like plum was observed. When the mesh was refined, some changes in the flow created a smoother plot over time.

For the user input material, where the Atwood number is reduced, the flow progressed with symmetric finger like protrusions and took larger time to stabilize.

The Transient State analysis was used to capture the fluctuations in the instability at every time step.

In the firsy two cases, for higher Atwood number the bubble like protrusions were uniform and significant in the refined case. Most of the denser fluid traversed along the wall because of the bubble plumes.

A more unstable mixing was observed for the third case at lower Atwood number. The multiple protrusions led to a better mixing of fluids in this case.

Due to the lower density difference between the user defined material and water, the instability was observed at a later time step and the mixing and stabilizing of the fluids took longer time.