Analysis of Rayleigh taylor instability using VOF multiphase method

Objective: To perform numerical analysis of Rayleigh taylor instability(air-water) using Volume of fluid(VOF) multiphase method in ANSYS fluent.

Rayleigh taylor instability:

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.

It may be modeled by two completely plane-parallel layers of immiscible fluid, the more dense on top of the less dense one 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 more dense material moves down under the (effective) gravitational field, and the less dense material is further displaced upwards.

Examples of Rayleigh taylor instability:

1. behavior of water suspended above oil in the gravity of Earth

2. mushroom clouds like those from volcanic eruptions and atmospheric nuclear explosions

3. supernova explosions in which expanding core gas is accelerated into denser shell gas,

4. instabilities in plasma fusion reactors and inertial confinement fusion.

CFD models based on mathematical analysis rayleigh taylor waves:

1. Single fluid model:

1.1 Potential flow model:

It provides a well-known framework for the study of surface modes and has been the basis of many classical stability results. Important aspects of this approach are that the most severe interfacial perturbations can be isolated by this model, giving usually an upper bound on instability growth. Furthermore, it pertains to nonlinear evolution in a relatively simple manner.

1.2 Buoyancy drag models:

buoyancy drag models use 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 two and three dimensions; and, as a rule, do not address demixing, also known as countergradient transport, i.e., reduction of total fluid masses within the mixing zone. To address the above problems, two-fluid (or multi-fluid) models have been proposed.

1.3 Besnard-Harlow-Rauenzahn (BHR) model:

A more advanced version of the single-fluid models. The BHR model is based on the evolution equations arising from second-order correlations and gradient-diffusion approxima-tions. Using a mass-weighted averaged decomposition, the original BHR model includes full transport equations for the Reynolds stresses, turbulent mass-flux velocity, density fluctuations, and the turbulence kinetic energy dissipation rate.

2. Multi-fluid model:

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

2.1 Two equation K-L turbulence model:

The K-L model was proposed by Dimonte and Tipton for describing the turbulent self-similar regime of RT and RM induced mixing. The starting point for deriving the model equations are the buoyancy-drag models.In this model K represents turbulence kinectic energy and L represents turbulent length scale.

2.2 Three equation K-L-a turbulence model:

The three-equation model of Morgan and Wickett was developed as an extension to the two-equation K-L model of Dimonte and Tipton by including a third equation for the turbulent mass-flux velocity.

2.3 Four equation K-L-a-b turbulence model:

This model extendes the K-L-a model by adding an evolution equation for the density-specific volume covariance which governs the primary production mechanism of the turbulent mass-flux and, therefore,needs to be modeled accurately to reflect the effects of the changes in the density fluctuations.

2.4 BHR-2 model:

This model uses an algebraic closure for the Reynolds stresses and the gradient diffusion approximation for the turbulent fluxes. The model introduces several additional terms in the governing equations of the turbulence length scale(L), turbulent mass-flux velocity (ai), and density-specific volume covariance (b).

Computational fluid domain:

For our present case, we have considered two dimensional fluid domain of dimensions 20x 40 mm where initially in the lower domain will be filled with air and upper domain will be filled with water separated by a interface line(can be seen in blue color).Share topology option was used to define the interface.

Mesh:

Case 1:

Element size: 0.5 mm

Number of elements: 3200

Case 2:

Element size: 0.3 mm

Number of elements: 8978

Case 3:

Element size: 0.1 mm

Number of elements: 80000

Multi-phase method: Volume of fluid(VOF) method

Volume of fluid(VOF) method: The VOF model is a surface-tracking technique applied to a fixed Eulerian mesh. It is designed for two or more immiscible fluids where the position of the interface between the fluids is of interest. In the VOF model, a single set of momentum equations is shared by the fluids, and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain. Applications of the VOF model include stratified flows , free-surface flows, filling, sloshing , the motion of large bubbles in a liquid, the motion of liquid after a dam break, the prediction of jet breakup (surface tension), and the steady or transient tracking of any liquid-gas interface.

Fluid Properties:

Primary Fluid: AIR

Secondary fluid: water

Numerical schemes:

Solver: Pressure-based solver

Pressure-based solver: The pressure-based solver employs an algorithm which belongs to a general class of methods called the projection method . In the projection method, wherein the constraint of mass conservation (continuity) of the velocity field is achieved by solving a pressure (or pressure correction) equation. The pressure equation is derived from the continuity and the momentum equations in such a way that the velocity field, corrected by the pressure, satisfies the continuity. Since the governing equations are nonlinear and coupled to one another, the solution process involves iterations wherein the entire set of governing equations is solved repeatedly until the solution converges.

SIMPLE: The SIMPLE algorithm uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field.

PRESTO: The PRESTO! (PREssure STaggering Option) scheme uses the discrete continuity balance for a \"staggered\'\' control volume about the face to compute the \"staggered\'\' (i.e., face) pressure. This procedure is similar in spirit to the staggered-grid schemes used with structured meshes . Note that for triangular, tetrahedral, hybrid, and polyhedral meshes, comparable accuracy is obtained using a similar algorithm. The PRESTO! scheme is available for all meshes.

PATCH:

After Standard intialization, patch option is used to allocated the fluid domain with their repective fluids(i.e water is allocated in upper domain and air in lower domain).Therefore in patch option volume fraction of water in lower doamin was set to 0 and volume fraction of water in upper domain is set to 1 as shown in the figure below.

After patch the fluid domains with their respective fluids, the same can be rechecked by ploting the contour of volume fraction of water as shown in figure below.From the figure below, we can observe that the upper domain is red(meaning volume fraction of water is 1 in upper domain) and lower domain is blue(meaning volume fraction of water is 0 in lower domain).

Residual plots:

Case 1:

Case 2:

Case 3:

For case 1,case 2 and case 3, the residuals varaition display repeated patterns,Therefore the solutions can be considered as converged.The absolute convergence criteria used is 1e-06.Since transient simulation, the calculation was run until steady solution is obtained(water settling in lower domain).

RESULTS:

Volume fraction of water contour animation:

Case 1:

Time step: 0.005

Case 2:

time step: 0.002

Case 3:

Time step: 0.001

Observations:

1. For all cases, rayleigh taylor instability can be observed at the interface of two fluids due to their difference in densities.As the mesh gets finer(case 3),the instability can be accuractely observed but takes larger time for heavier fluid to settle at the bottom.

2. Initially the lighter fluid tends to move up in a bulk manner(bubbles) as heavier fluid pass through them.

3.Since the difference in densities between the fluids is large(i.e. atwood number is large), appropriate turbulence modellling is required in order to accrately observed the behaviour.

Atwood number(A):

Atwood number, abbreviated as A, is a dimensionless number that describes density difference between two adjacent fluids with a common interface. Atwood number is an important parameter in the study 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=frac{rho_1-rho_2}{rho_1+rho_2}`

where, `rho_1`= density of heavier fluid(kg/m3)

`rho_2`= density of lighter fluid(kg/m3)

For our case,

`A= frac{998.2-1.225}{998.2+1.225}=0.9975`

Effects of Atwood number on rayleigh taylor instability:

Compressibility has two effects on Rayleigh-Taylor instability. The initial density stratification caused by compressibility plays a stabilizing role, while the expansion-compression effect of flow plays a destabilizing role. The overall impact of compressibility on the evolutions of RT instability depends on competition between two effects.

1. For the case of small Atwood number, the stabilizing effect of initial density stratification plays a dominant role. The density difference between the two sides of the initial interface is very small, and the density stratification of the upper and lower layers is significant at At = 0.1. Therefore, the development of RT instability is suppressed with the increase in the Mach number.

2. For medium Atwood numbers, the spike velocity is generally suppressed, but the bubble velocity has different possibilities depending on the competition between the two effects.

3. With the increase in the the Atwood number, the stability effect of density stratification becomes weaker. The effect of compressibility on RTI is dominated by the expansion-compression effect at large Atwood number.

Why steady state approach might not be suitable on these simulation?

staedy state approach is not suitable because we are not considered about the end results of the simulation(settling down of fluids) rather we are considered in observing the instabilities caused with respect to time which is why transient approach is implemented.

Google drive link: https://drive.google.com/drive/folders/1iOkcY4GW-gUkCBmsve0jjpmWzvsMjYbR?usp=sharing