Rayleigh Taylor Instability Challenge

Introduction:

The Rayleigh–Taylor instability, or RT instability (after Lord Rayleigh and G. I. Taylor), is an instability of an interface between two fluids of different densities which occurs when the lighter fluid is pushing the heavier fluid. It is a dynamic process whereby the two fluids seek to reduce their combined potential energy. The ensuing turbulence and mixing has far-reaching consequences in many natural and man-made flows.

Examples:

• Behavior of water suspended above oil in the gravity of Earth,
• Mushroom clouds
• Volcanic eruptions
• Atmospheric nuclear explosions,
• Supernova explosions in which expanding core gas is accelerated into denser shell gas.
• Instabilities in plasma fusion reactors
• Inertial confinement fusion.

________________________________________________________________________________________________________________________

Theory:

Simplest Explanation of the Occurrence of Rayleigh-Taylor Instability

• Imagine the ceiling of a room plastered uniformly with water to a depth of 1 meter. The layer of water will fall. However, It is not through lack of support from the air that the water will fall, the pressure of the atmosphere is equivalent to that of a column of water 10 meter thick, quite sufficient to hold the water against the ceiling, But in one respect the atmosphere fails as a supporting medium and fails to constrain the air-water interface to flatness, No matter how carefully the water layer was, prepared to begin with, it will deviate from planarity by some small amount. That portion of fluid which lies higher than the average, experience more pressure than is needed for their support. It begin to rise, pushing aside water as they do so, A neighboring portion of the fluid, where the surface hangs a little lower than average, will require mow than average pressure for it, support, It 4 begins to fall. The air cannot supply the variations in pressure from place to place necessary to prevent the interface irregularities from growing; the initial irregularities therefore increase in magnitude, exponentially with time at the beginning. The water which is moving downward concentrates in spikes. The air moves upward through the the floor, the same layer of water stable. Irregularities die out. Water in round topped columns. The water falling on the floor would have been perfectly thus we can infer a simple criterio.1 for the onset of Taylor instability at the interface between two fluids of different densities: If the heavy fluid pushes the light fluid, the interface is stable. If the light fluid pushes the heavy fluid, the interface is unstable. A criterion of comparable simplicity governs the onset of Kelvin-Helmholtz instability: The interface between two fluids is unstable if there is a jump in the tangential component of the velocity across the interface. These two criteria are among the most basic principles’ in the subject of interface instability.

Phenomenology

There is a complex phenomenology associated with Tailor unstable interface. This includes the formation of spikes, curtains and bubbles. The development of Helmholtz instability on the side of the spikes, competition among the bubbles leading to their amalgamation, formulation of droplets, entrainment and turbulent mixing and possible chaotic limit with a fractionized interface.

Organized description of the growth of instability in to number of stages can be as follows:

• Stage: 1

> If the initial perturbations of the interface or the velocity are extremely small, the early stages in growth of instability can be analyzed using the linearized form of dynamical equations for the fluid. The result is that small amplitude perturbations of wavelength grow exponentially with time. When the amplitude of the initial perturbations grows to a size of the order 0.1 to 0.4, substantial deviations from the linear theory are observed.

• Stage: 2

> During the second stage while amplitude of perturbation grows nonlinearly to a size of the order, the development is strongly influenced by three dimensional effects and the value of the density ratio, or Atwood number A = (_{H – L}) / (_{H + L}), where _{H =} Density of heavy fluid and _L = Density of light fluid.

> IF A 1, light fluid moves in to heavy fluid in the form of round topped bubbles with circular cross sections. Note that two-dimensional plane bubbles are unstable to perturbations along the axis perpendicular to the plane of the bubble and a trough having a plane bubble as a cross section will break up in to three-dimensional bubbles. The heavy fluid will form spikes and walls or curtains between the bubbles, so that a horizontal section would show a honeycomb pattern. If A 0, one will instead find a pattern more like two sets interpenetrating bubbles.

> Note that this will differ from simple pattern of bubbles and spikes that we think of two-dimensional. I will nevertheless often refer to bubbles and spikes for the sake of simplicity.

• Stage: 3

> The next stage is characterized by development of structure of spikes and interactions among the bubbles. These phenomena can originate from several sources. There is a non-linear interaction among initial perturbations of different frequency. Also Helmholtz instability along the sides of the spike can cause it to mushroom, increasing the effect of drag forces on the spike. This effect is more pronounced at low density ratios. There is some experimental evidence for bubble amalgamation, a process in which large bubbles grow larger and move faster.

> The process of heterogeneities in various physical quantities can modify the shape and speed of the bubbles and spikes to a degree which depends on the strength and length of the scale of heterogeneity.

• Stage: 4

> In the final stage, we encounter the breakup of spikes by various mechanisms. The penetrations of bubbles through a slab of fluid of finite thickness and other complicated behavior that leads to a regime of turbulent or chaotic mixing of the two fluids.

Vorticity explanation:

Here is the image to describe Rayleigh-Taylor instability by means of velocity.

• It is Visualization of unstable Rayleigh-Taylor instability configuration, where baroclinic torque at the interface creates vorticities and induces a velocity field that increases the baroclinic torque.

Here omega is vorticity, p is pressure, rho is density, u is velocity and g is gravity.

Thick circular arrows represents velocity field created by vortex.

________________________________________________________________________________________________________________________

CFD models:

There are some practical CFD models have been used based on mathematical analysis of Rayleigh-Taylor waves.

Models used for Rayleigh-Taylor Instability analysis:

• Early Time Growth
• Potential Flow
• Buoyancy-Drag Balance

1) Early-time growth

It was noted by Taylor (1950) that the two superimposed liquids with and interface normal to the acceleration could exist in a state of unstable equilibrium. However, if small corrugation were to exist on the interface, a growing standing wave pattern can be observed.

2) Potential flow

Development of unstable interfaces beyond the linear phase was first studied by Davies & Taylor (1950) both papers begin by considering a two-dimensional incompressible potential in a vertical cylinder with an initially flat-free under surface. This is Rayleigh Taylor instability with Atwood No A = 1. The velocity potential must satisfy Bernoulli’s equation in cylinder-polar coordinates.

3) Buoyancy-drag Balance

The complexity of directly evaluating a potential flow to calculate growth behavior can be avoided from the first principle, the following equation describes how buoyancy competes to accelerate a body, which is, in this case, is the bubble of lower layer fluid penetrating the upper layer.

Following are the types of Instability:

• Kelvin–Helmholtz instability
• Rayleigh–Taylor instability
• Plateau-Rayleigh instability (similar to the Rayleigh–Taylor instability)
• Richtmyer-Meshkov instability (similar to the Rayleigh–Taylor instability)
• Shock Wave Instability
• Benjamin-Feir Instability (also known as modulational instability)
• Bénard instability

________________________________________________________________________________________________________________________

AIM:

Objective of this project is to survey Rayleigh-Taylor instability; describing phenomenology occurs at Taylor unstable interface with two-dimensional CFD model and reviewing attempts with different mesh sizes to understand these phenomena quantitatively.

________________________________________________________________________________________________________________________

Modeling Approach:

• For modeling Rayleigh-Taylor instability, we have to take two different phases, one as water and another as Air.
• Here, we will create a model in a way that heavy fluid will fall on the light fluid due to gravity and accelerations between two fluids.
• IT can be modeled by two completely plane parallel layer of miscible fluid, the more dense on the top of the less dense and both subjects to earth gravity.
• The equilibrium here is unstable to any perturbation are disturbances of the interface.
• If a parcel of heavier fluid is displaced downward with an equal volume of lighter fluid displaced upward.
• The potential energy of the configuration is lower than the initial stage.
• Thus the disturbance will grow to further release potential energy as the more dense material moves downward under the effective gravitational field.
• And the less dense material is further displaced upward.

Geometry:

• Create two rectangle attached to one another in X-Y plane.
• Create Surfaces from these two rectangles as different component in SpaceClaim.
• Give the name to upper surface as Water and lower surface as Air.
• Using “SharePrep” or “Share” option, share the boundaries between two fluids for conformal Mesh.
• Make sure that “Share Topology” is enabled.

See the detailed geometry as given below:

________________________________________________________________________________________________________________________

Case Setup Details:

• Solver Type: Pressure Based
• Time Domain: Transient / Unsteady
• Viscous Model: Laminar
• 2D Space: Planar
• Gravity: - 9.81 m/s² in Y-Direction
• Multiphase Model: Volume of Fluid
• of Eulerian Phases: 2
• Formulation for Volume Fraction: Implicit

• > Here we will discuss three different cases with different Mesh sizes.
• > Baseline Mesh case M1 and than refined mesh cases will be as case M2 & M3.

________________________________________________________________________________________________________________________

Meshing Details:

Linear Mesh size according to size of elements for different cases is as given below:

> Case - M1: 0.5 mm

> Case - M2: 0.3 mm

> Case - M3: 0.2 mm

Images of Meshing for above three cases are as given below:

________________________________________________________________________________________________________________________

Pre-Processing & Solver Setting:

• Set the solver to pressure based with transient and absolute velocity formulation.
• We have to give the Gravitational acceleration to the fluid as 9.81 m/s² In the –Y direction.
• Create materials as Air and Water (liquid) for phases.
• Now after creating material, Multiphase method as “Volume of Fluid”.
• In list of the Phases, we have to assign primary material as “Air” and Secondary material as “Water” for changing of the Phases for analysis.
• Then after, initialize with Standard Initialization.
• After initialization, we need to Patch the Phases to Zones as given below:

• For that process, Patch the Water with Variable for Water Volume Fraction called and keep the value =1 for water and Patch it.
• After Patching the Water, Patch the Air with =0.
• Here =0 will be given for Primary fluid and =1 for the Secondary Fluid.
• Then create contours of Phases as Volume Fraction and Water for all methods.
• Also create Solution Animation for viewing the animation of contours.
• Then Run the simulation with appropriate Time Step value and Iterations up to convergence.
• IF we modify any settings then after initialization every time Patching Process mentioned above is compulsory for proper running simulation without any errors.

NOTE:

• If we use Explicit Method, then we have to ensure the Time Step should be selected in a manner that CFL (Courant-Friedrichs-Levy) number is

Because, SRS (Scale-Resolving Simulations) model operates at Grid limit.

• In Implicit method also we have to ensure that value of Ushould not exceed the value of. Otherwise there may be errors regarding floating point and other simulation errors, as the Time Step will be higher than the elemental length.

________________________________________________________________________________________________________________________

Results:

Results for different cases as M1, M2 and M3 are as given below:

CASE: M1

Residuals:

• It converges at near about 1600 iterations.

Contours (Volume Fraction: Water):

CFD Post:

CASE: M2

Residuals:

• It converges at near about 6000 iterations.

Contours (Volume Fraction: Water):

CFD Post:

CASE: M3

Residuals:

• It converges at near about 14000 iterations.

Contours (Volume Fraction: Water):

CFD Post:

___________________________________________________________________________________________________________

User Defined Fluid:

• Now we will see another analysis with water and User defined fluid.

Material Specifications:

• Here heavy fluid will be water and UDF will be light fluid instead of Air, because Water density is much higher than UDF.
• We give the name of this case as Case: M4 (UDF-Water)
• Following above pre-process, we have used UDF as less dense fluid instead Air and run the simulation.

Results:

Results for case M4 are as given below:

CASE: M4

Residuals:

• It converges at near about 21500 iterations.

• Contours (Volume Fraction: Water):

CFD Post:

___________________________________________________________________________________________________________

Comparison of RT. Instability with Time progress:

• Comparison of RT. Instability with time progress for case M3 and M4 as given below.

NOTE: Mesh size in case of M3 and M4 is same for better comparison.

CASE: M3 (Water – Air)

CASE: M4 (Water – UDF):

• We can see that in case of UDF (M4), uniform generation of spikes can be observed.
• Also mixture of both fluids accelerates in equal manner.
• Bubbles like plums in case of M3 are more than M4.
• In case of M4, we can see the starting interface creates spikes like mushroom clouds or finger forms.
• In case of M3 less dense material is AIR and in case of M4, less dense material id UDF, but density of UDF is much higher than AIR, so it is taking more time for getting stability.

___________________________________________________________________________________________________________

Atwood Number:

Here we will see the Atwood number and its effect on the behavior of RT. Instability.

Atwood Number:

• It is a dimensionless number which describes density difference between two adjacent fluids with common interface.

Atwood number for case M1, M3, M3 common as (1) and for M4 (Water-Air) as (2), respectively:

Variation of Atwood Number and Behavior of Instability:

• IF Atwood number is positive, interface is unstable.
• If Atwood number is negative, interface is stable.
• If Atwood number is close to 0, RT. Instability flow takes form of fingers of fluid.
• If Atwood number id close to 1, lighter fluid below heavy fluid takes the form of large bubbles like plumes.

• Here Atwood number is near to 1, in case of M1, M2 and M3, that is why we can see the bubble like plumes in the refined case.

• In case of M4, Atwood Number is near to 0.5, so we can see the very good mixing effect of both the fluid and finger forms as well as spikes.

• In case of M4 spikes and bubbles are small compared to other case but the stability creates very smooth flow of both the fluids, as density ratio is low.

________________________________________________________________________________________________________________________

Conclusion:

• Here we have taken Transient mode for simulation, which is because steady mode only gives the final results and for observation of instability of the fluids, we have to track it during progress of time. So for behavior of Rayleigh-Taylor Instability with respect to time, transient mode is most needed.
• Here also instability of the fluid is highly dependent on small time steps that are why we have taken implicit method for solution.
• Also implicit method is more complex to program and require more computational effort in each time step size, so numerical stability is increased.
• We can see the difference in mesh and results, that with more refined mesh we can see the formation of bubbles and spikes with increase in instability at the interface.
• Also with increase in cell count or with refinement of mesh, much smaller time step size is needed to solve.
• As mesh is more refined, results will be accurate and clear contours we can get.
• High Atwood number or High Density Ratio becomes easily stable compared to low Atwood number or low Density Ratio.
• Instability of the flow in high Atwood number has larger bubbles and spikes compared to low Atwood number.

________________________________________________________________________________________________________________________

NOTE: Links for the video animations for M1, M2, M3 and M4 case, are given as below:

M1:https://drive.google.com/open?id=14JXQOARoCBJWUhZbJq1-mvIZKB6cin_4

M2:https://drive.google.com/open?id=1ucJM6p_QYfju-GB2UYQDUSG-Yj9X6wft

M3:https://drive.google.com/open?id=10f2YHi6HrScxt-w3HMVSNleTcVnvsAyJ

M4:https://drive.google.com/open?id=1XRYJ9aVAUgqHEbKTeO9xmTOaoYQlLUi6

Also link for the Report in Ducument and Pdf version given as below:

Doc:https://drive.google.com/open?id=14jq7Svd6BZ02Cb8bZD4THtBX2UxWP6kv

Pdf:https://drive.google.com/open?id=1eHgncMHInjfXuMMi4HPN3qA2rawxc5f1

________________________________________________________________________________________________________________________

Thanks & Regards,

Dharmesh Joshi