Week 5 - Rayleigh Taylor Instability

Aim:

To perform simulation for observing Rayleigh Taylor Instabilities (RTI) in a column of two immiscible fluids.

Simulation should be performed for 2 cases - water and air, water and a user defined material.

For water and air case, base mesh should be taken as 0.5 mm and simulation should be performed. Later, mesh needs to be refined and another simulation should be performed and results must be compared.

For case 2, refined mesh size used in case 1 should be used and simulation should be performed.

Atwood number for both the case should be calculated and how the variation in Atwood number affects the instabilities is to be found.

Introduction & Theory:

RT instability is an instability of an interface between two fluids of different densities that occurs when lighter fluid is pushing heavier fluid.

Examples are water suspended over oil, mushroom clouds from volcanic eruptions and nuclear explosions.

(Refer Pic 1)

When a heavy fluid which is on top of lighter fluid moves downward due to gravity, the potential energy of the system changes. This phenomenon is not just for fluids which are static but also applies to fluids that are accelerated.

The above picture gives a clear view of how an RTI looks. As RTI develops, the disturbances at the interface progress from linear growth phase to non linear phase eventually developing bubbles flowing upwards and spikes falling downwards. The density disparity between fluids determines the subsequent instability flow.

CFD models based on mathematical analysis of Rayleigh Taylor waves:

1. Mass diffusion model: The mass diffusion equation describes the transport of mass and energy due to model viscosity and hence drives actual mixing in the fluids. Even if any configuration is RT unstable, the diffusion process could be responsible for significant mixing if the viscosity is large.

2. Multi-Fluid models: These models use a separate set of equations for each fluid in addition to the main flow equations and provide an accurate modelling framework for demixing by correctly capturing the relative motion of the different fluid fragments. The first step for the development of turbulence models is to perform Reynold averaging of the governing equations and resultant terms.

2.1. Two equation K-L turbulence model: The K-L turbulence model involves equations for turbulent kinetic energy K and turbulent length scale L. The starting point for deriving the model equations is he buoyancy-drag models for the self-similar growth of RT instabilities.

2.2. Three-equation K-L-a turbulence model: The three-equation model was developed as an extension to the two-equation K-L model by including a third equation for the turbulent mass-flux velocity.

2.3. Four-equation K-L-a-b turbulence model: The four equation model was developed as an extension to the three-equation K-L model by including a fourth equation for density-specific volume covariance. The model can accurately reflect the effects of the changes in the density fluctuations.

3. Single fluid model:

3.1. Potential flow model: This model is valid for all instability changes from the linear stage through the early non-linear to the asymptotic but limited to Atwood number =1. In this model, single-mode perturbation expands the flow equations to second-order around the bubble tips, hence obtaining an ordinary differential equation for the bubble velocity.

3.2. Buoyancy drag model: It is valid for every Atwood number but limited to the non-linear stages of the instability. It uses 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 2 or 3 dimensions, and as a rule, do not address demixing, also known as counter gradient transport, i.e. reduction of total fluid masses within the mixing zone. To address the above problems, two-fluid(or multi-fluid)models have to be used.

3.3. Bernard-Harlow-Rauenzahn (BHR): It is a more advanced version of the single-fluid model. The BHR model is based upon the evolution equations arising from second-order correlations and gradient-diffusion approximations. Using a mass-weighted averaged decomposition, the original HR model includes full transport equations for the Reynold stresses, turbulent mass-flux velocity, density fluctuations, and the turbulent kinetic energy dissipation rate.

Geometry creation:

To perform simulation for RT instability, two different fluid columns, one with heavier fluid on top (water in this case) and lighter fluid at the bottom is required.

For this, a 2D geometry containing two rectangles one over the other is created. These two rectangles have length and breadth of 30 x 30 mm respectively.

(Refer Pic 2)

Similarly draw another rectangle with same dimensions on top the existing rectangle.

Use the pull tool to convert the sketch into volume.

Final geometry is as below

(Refer Pic 3)

Name the top surface as "Water" and the bottom one as "Air".

(Refer Pic 4)

As there are two different surfaces, conformal meshing needs to be done on this. To do that, share topology must be enabled.

In the "Workbench" module, click on "Share" and select the common edge.

(Refer Pic 5)

Also, enable "Share" in the Tree.

(Refer Pic 6)

Pre-processing:

Before meshing, assign names to the regions to which boundary conditions should be given.

Case 1 - Water and air being the two fluids with mesh size of 0.5 mm.

Water at the top and air at the bottom.

(Refer Pic 7)

Next select Mesh -> Method. Select the geometry.

Use Quad elements of linear type.

(Refer Pic 8)

Input the mesh element size which is 0.5 mm.

(Refer Pic 9)

Generated mesh is as below

(Refer Pic 11)

Next check the mesh element quality.

(Refer Pic 10)

Mesh quality is good as all elements are having element quality of more than 98%.

Total no. of elements and nodes are as below

(Refer Pic 12)

As the share topology is enabled, it can be see that the nodes are shared between the upper and lower surfaces at the interface.

(Refer Pic 13)

Case 2 - Water and air being the two fluids with a refined mesh size of 0.1 mm.

The settings used in case 1 are also applied to this case except now the mesh size is changed to 0.1 mm.

(Refer Pic 14)

Final generated mesh

(Refer Pic 15, 16)

Total no. of elements and nodes are as below

(Refer Pic 17)

Element quality is as below

(Refer Pic 18)

Mesh quality is good as all elements are having a mesh quality of more than 96%.

Case 3 - Water and user defined material being the two fluids using refined mesh of element size 0.1 mm.

In this case, instead of air, a material is created in database with a density of 400 kg/m3 and dynamic viscosity of 0.001 Pa-sec.

As the mesh element size is same as in case 2, the same mesh can be used in this case also.

Solving:

Open fluent solver and before starting to assign input values, mesh check must be done.

(Refer Pic 19)

If the mesh quality is good, it displays the below message.

(Refer Pic 20)

Cross check all the regions/surfaces in "Display" option in "Setting up Domain" module.

(Refer Pic 21)

Next in "Setting up Physics" module, select the appropriate conditions in which the simulation must be performed.

Simulation is considered transient and pressure based.

As the water moves down due to gravity, gravity option must be enabled and its value of 9.81 m/s2 must be provided as input in Y direction.

(Refer Pic 22)

Here, the positive Y axis is towards up, so the gravity value should be negative.

The flow is laminar.

(Refer Pic 23)

As the problem has two phases, multiphase model must be used to simulate. In this, "Volume of Fluid" model is selected.

Implicit approach is used for this simulation and no. of eulerian phases is given as 2.

(Refer Pic 24)

Now materials must be assigned to the top and bottom areas.

To do that, both the materials (air and water for case 1 and 2) should be shown in materials database.

Air is the default material. Water must be added to the "Fluent fluid materials" from the fluent database.

(Refer Pic 25)

For case 3, a new material must be created in the fluent material database with a density of 400 kg/m3 and dynamic viscosity of 0.001 Pa-sec.

To do this, open the material section and click on fluent database.

(Refer Pic 27)

Select any material and click copy.

(Refer Pic 28)

A dialog box pops up stating to give name for new material.

(Refer Pic 29)

Click OK and click. By doing this, all the material properties of the selected material are copied into the new material created and is displayed in material database window.

Now modify the material properties with the required values and click "Change/Create".

(Refer Pic 30)

Next in the "Phases" section, click "List/Show All" which displays all the phases in the problem.

(Refer Pic 31)

As the problem has two phases, both of them are displayed. Open each phase and assign corresponding materials to it. Air is assigned to "Air-Primary phase" for case 1 and 2 whereas user defined material is assigned in Primary phase in case 3. Water is assigned in secondary phase.

(Refer Pic 32, 33)

Next in solving module, perform standard initialization.

(Refer Pic 34)

After doing initialization, patch option is enabled. Select the patch option.

This patch option assigns the volume fraction ratio in the initial stage to each phase.

So, in initial stage, volume fraction of water in water side is 1 and at air side is 0.

After assigning, click on patch.

(Refer Pic 35, 36)

Create a volume fraction contour plot in post processing module.

(Refer Pic 37)

As observed in above picture, the top portion has a volume fraction of 1 indicating that it is filled with water as assigned in patch section.

The bottom portion's volume fraction is 0 indicating there is no water present there but only air.

Next step is to solve the problem.

Results & Conclusions:

Transient simulation is performed in this problem. Steady state simulation is not used for this because the RT instabilities depends on time. As time progresses, interaction between heavier and lighter fluid keeps changing and instabilities are observed. Steady state simulation is only used when the end result or final state is of main focus. But in this problem, due to gravity, the disturbances at the interface progress from linear phase to non linear phase with time thereby developing bubbles and spikes. These bubbles and spikes can be observed more clearly only when transient simulation is performed thereby enabling to understand the flow dynamics.

Case 1 - Air and water being two fluids with mesh size of 0.5 mm.

As the problem is transient, a fixed time step size of 0.005 seconds is taken with no. of time steps as 300.

Simulation is ran for 8000 iterations and below are the plots

(Refer Pic 38)

As observed, residuals are not increasing indicating convergence.

(Refer Pic 39)

As seen, water gets settled at the bottom (volume fraction of 1) after some time and air moves up (volume fraction of 0).

The green and yellow regions has a volume fraction of 0.5-0.7 indicating bubbles.

Link for solution animation - https://drive.google.com/file/d/1msvOn96LCK1Kivuerx9cht46rW_bRj4v/view?usp=drive_link

Case 2 - Air and water being fluids with refined mesh of 0.1 mm.

In this case, mesh is refined and simulation is performed initially with a fixed time step of 0.005 seconds.

Simulation is ran for 5000 iterations and below are the plots

(Refer Pic 40)

As seen above, the residuals are increasing indicating that the solution is diverging which gives inaccurate results.

Time step is decreased now to 0.0025 seconds and simulation is ran again for 2000 iterations.

(Refer Pic 41)

The residuals are again increasing indicating that the solution is diverging.

Now again the time step is reduced to 0.0005 seconds and simulation is ran for 1200 iterations.

This is an iterative process and until the solution achieves convergence, time step must be changed.

(Refer Pic 42)

As observed, the continuity and volume fraction equations are diverging.

As the solution is diverging for fixed time step, variable time step is used.

In the "Run calculation" option, change the time step from fixed to variable.

Minimum and maximum time step should be given as input.

In this case, minimum time step is given as 0.001 seconds and maximum time step as 0.0025 seconds.

Courant number is given as 0.2.

(Refer Pic 43)

Even though implicit approach is used to solve the problem, global courant number should be given as input. This is because the governing equation for solving volume fraction is solved using explicit approach by default. Here the global courant number should be less than 0.5.

Simulation is ran for 300 time steps with variable time steps and the results are below

(Refer Pic 44)

As seen in above image, the residuals are not increasing indicating that the solution is converged.

(Refer Pic 45)

As observed from plot, the bottom and top portions are occupied with water and air respectively.

Link for solution animation - https://drive.google.com/file/d/14bn2g5hbPWSX3RuHQE3r0xLw4LQ_Pn0e/view?usp=drive_link

As the mesh is much more refined here, the RT instabilities can be seen very clearly i.e., during initial phase at the interface , bubbles and spikes can be seen here whereas in the baseline mesh (0.5 mm), these bubbles and spikes are not observed. Also, how the mixing of both phases happen can be clearly seen. As the mesh size is less, the areas of void fractions ranging from 0.2-0.8 are clearly distinguished i.e., the bubbles of size 0.1 mm in diameter are observed and distinguished with other bubbles. If the mesh size is large, the differentiation between different bubbles is difficult to observe.

Case 3 - Water and user defined material being two fluids with refined mesh of 0.1 mm.

In this case also, variable time step with the above settings is used.

(Refer Pic 46)

As seen, residuals are not increasing indicating solution is converged.

(Refer Pic 47)

In this case, as the lighter fluid has higher density compared to air, it takes more time to reach steady state.

Link for solution animation - https://drive.google.com/file/d/1Y4_pP8MMFrU2XbeDE9cQxgIA4gY2MLO6/view?usp=drive_link

Atwood number:

Atwood number is a dimensionless number used in the study of hydrodynamic instabilities in density stratified flows.

It is given by below formula

(Refer Pic 48)

For case 1 and case 2,

density of heavier fluid (water) = 998.2 kg/m3.

density of lighter fluid (air) = 1.225 kg/m3.

Atwood number would be = (998.2 - 1.225)/(998.2 + 1.225)

= 0.997

For case 3,

density of heavier fluid (water) = 1000 kg/m3.

density of lighter fluid (user defined material) = 400 kg/m3.

Atwood number would be = (998.2-400)/(998.2+400)

= 0.427

For higher Atwood numbers, the growth of instability is very fast and non linear. Bubbles and spikes are seen but not like mushroom clouds.

whereas for lower Atwood numbers, growth of instability is relatively slow and linear. The instabilities are more developed into mushroom like clouds in the beginning and later disperse.

Higher Atwood number flows have less intermolecular forces between the fluids which makes the heavier fluids displace lighter fluid quickly thereby achieving steady state faster.

whereas the lower Atwood number flows have higher intermolecular forces between the fluids thereby taking more time to achieve steady state.