RAYLEIGH TAYLOR INSTABILITY

Rayleigh Taylor Instability

objective

1. To find out some of the practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves and explain how these mathematical model have been adapted for CFD calculation.
2. To Perform the Rayleigh Taylor instability simulation for 2 different mesh sizes with the base mesh being 0.5 mm. Observe the differences in the animation and explain the effect of mesh size on simulation.
3. To explain why a steady-state approach might not be suitable for these types of simulation.
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.

Rayleigh Taylor Instability:

The Rayleigh-Taylor instability, or RT instability, is instability of an interface between two fluids of different densities which occurs when the lighter fluids is pushing the heavier fluid. Example - water suspended above oil in the gravity of Earth, mushroom clouds like those from volcanic erruptions and atmospheric nuclear explosions, supernova explosions in which expanding core gas is accelerated into denser shell gas, instability in plasma fusion reactors and inertial confinement fusion.

Practical CFD model based on Rayleigh Taylor waves:

1. 1D RAYLEIGH-TAYLOR MIXING:

The turbulece models have been aplied to the simulation of simple 1D RT mixing cases with 3:1 density ratio(`rho1=3g/(cm^3)` and `rho2=1g/(cm^3)`). The computational domain extend [-8,20] cm with the heavy fluid placed on the left side of the domain and the initial interface at x=0. Unless otherwise stated , the computational grid consists of 100 cells and the adiabatic exponent is `gamma=5/3` for both fluids. The following relation is satisfied in all cases `Aog=1`; thus, for the `3:1` density ratio the gravitational acceleration is `g=2(cm)/s^2`.

2. TILTED-RIG RAYLEIGH-TAYLOR MIXING:

`he tilted-rig test case originates from a series of experiments performed at the Atomic weapons establishment in the United Kindom in the late 1980 to study the mixing between two variable density fluids included by the Rayleigh-Taylor instability. In the experiment, a tank containing two fluids of different densities, a heavy fluid placed in a tank below a lighter fluid(a stable configuration), is accelerated downwards between two parallel guide rods by firing rocket motors. The downward acceleration caused to the tank effectively changes the direction of "gravity"(external body force) so that the system becomes RT unstable, causing the two fluids to mix. The acceleration from the rocket motors was not conctant but averaged approximately 35 times normal gravity. It eventually attains a roughly constant value, but the time during which the acceleration is varied is significant and still needs to be considered. One approach is to directly incorporate the measured values into the simulation[i.e. a variable gravity g(t)]. This works well for incompressible Navier-Stokes solver; however, it can create problem for compressible codes. An alternative approach is to make use of a constant acceleration for which a nondimensional time, `tau`, is used for comprerison with experiments result:

`tau=int{(A*g)/(L)}^(1/2)dt+delta`

3. SPRAY BREAKUP MODELING:

i. CHU breakup model:

The Rayleigh-Taylor instabilities from the basis of the model formulated by the Chu. The theoretical correlation was constructed in order to predict droplet sizes based on an exponential function which was uesd in fuel-coolant interaction models. The model implemented in this represents a modified version of Chu's original model in that it includes a Taylor-seires expansion of the exponential term to first order, and subsequently an explicit discretization in time in which drop diameter correlation is written as:

`(dD)/dt=-Co*C1*We^(C2)(Urel*epsilon^0.5)tau^((C1-1)),`

where Co is a function of the density ratio `epsilon`

`Co=0.1708-0.149epsilon^0.5, epsilon=(rhog)/(rhol)`

and `tau**` is a non-dimensional acculated time within the droplet integration time

step when using sub-cycling method for solving Eq.(1):

`tau**=(epsilon^0.5Urel*t)/(Do)`

The default values for the cnstants are C1= 0.772, C2=0.246.

ii. Wave Model:

The Wave model(KH model) was proposed by Reitz. The model is bassed on a first-order linear analysis of a Kelvin-Helmholtz instability growing on the surface of a cylindrical liquid jet with initial diameter `ro` that is penetrating a stationary incompressible gas with a relative velocity `Urel` .

The instability at the surface is a small axisymmetric displacement `eta`.

`eta=eta(i)e^(ikz+wt),`

whre `eta(i)` is an infiniteesimal initial perturbation, k is the wavenumber, `omega` is the growth rate. Wavelength `lambda` is related to the wavenumber, as follows:

`lambda=(2pi)/k`.

iii. KHRT model:

It is known that liquid droplets in diesel sprays have very high velocities and decelerate rapidly due to drag forces. In this case, Rayleigh-Taylor instability may also play an important role in droplet breakup mechanism in addition to the kelvin-Helmholtz instability as shown in fig.

In the KHRT model, the aerodynamics force on the drop flattens it into the space of a liquid sheet, and the decelerating sheet breaks into large-scale fragments by means of RT instability. KH waves with a much shorter wavelength originate at the edges of the fragments, and these wave break up into micrometer-sized drops.

Now to visualize the rayleigh taylor instability phenomenon we have uased ANSYS 18 -

To simulate the physics we have used the Transient, Multiphase modelling (VOF) technique with two different fluid -water and air. The geometrical domain is shown below-

fig1: geometry over which simulation is done. water is in upper half and air is in lower

Case1

The simulation is done for two different grid size with baseline mesh be 0.5mm of length.

SIZE	ELEMENT	NODE
0.5mm	3200	3321
0.13mm	47432	47895

above table shows the details of the meshes of two cases.

Mesh 1

fig2: baseline mesh with 0.5mm grid length

After transient simulation with time step of 0.025 sec , the animation of volume fraction contour is below-

The time steps should be very small so that it can capture the smallest scale length possible at the interface of the instability.

Mesh 2

fig3: Mesh 2 is further refined to refine the results more. grid size- 0.13mm

After the transient simulation we get the animation of the flow below-

after both simulation we can see that Mesh 2 is showing more details and capture the interfacial instabilty more accurately. This is because in Mesh 2 both time scale and Length scale was smaller then Mesh 1, that is why it capture more details. This is also called Mesh InDependency test, where we refine the mesh in search of more accurate results.

Case2 (User-defined material)-

In this simulation we changed the material density and investigate the behaviour. Heavier density (`rho=400(kg)/m^3 and mu=0.001(kg)/(m-s)`) being in the upper counter and air being in the bottom.

Atwood Number

A dimensionless number in fluid mechanics that tells the behaviour of two fluid interacting with different density.

`A=(rho1-rho2)/(rho1+rho2)=(997-1.225)/(997+1.225)=0.9955` .......... Case1

`A=(400-1.225)/(400+1.225)=0.99389`...............................................Case2`

we can see that 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 development of an R-T mixing region as a three-step process-

1. Initially an exponential growth of infinitesimal perturbations that correspond to linear stanility analysis.
2. At amplitudes of about one-half of tge wavelength,the instability saturates and longer wavelengths (that have yet to reach saturation) take over. Emmons et al. (1960) coined the term 'bubble competition' to describe this regime.
3. Eventually, through nonlinear mode interaction and successive wavelength saturation, a self-similar R-T mixing layer is formed with a scale similarity as `g*t^2`.