Rayleigh Taylor Instability Challenge
Aim - To conduct the Rayleigh Taylor CFD simulation.
Introduction - Rayleigh Taylor Instability (RTI) is a phenomenon that occurs at the interface of fluids of different density wherein the lighter fluid pushes the heavier fluid due to the density difference. It is a dynamic process wherein the system tries to reduce its combined potential energy.
Theory - Initially the system is in a state of hydrostatic equilibrium, in which the heavier fluid(ex: water) sits atop a lighter fluid(ex: air) in a constant gravitational field. As the RTI develops, the initial perturbations progress from a linear growth phase, eventually developing "plumes" flowing upwards(in the gravitational buoyancy sense) and "spikes" falling downwards. In the linear phase, equations can be linearized and the amplitude of perturbations is growing exponentially with time. In the non-linear phase, the perturbation amplitude is too large for the non-linear terms to be neglected. In general, the density disparity between the fluids determines the structure of the subsequent non-linear RT instability flows (assuming other variables such as surface tension and viscosity are negligible here).
Examples include supernova, water suspended atop oil, mushroom clouds formed due to volcanic eruptions.
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 Rayleigh-Taylor unstable, the diffusion process could be responsible for significant mixing if the viscosity is sufficiently large.
2. Single fluid model:
2.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.
2.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.
2.3. Besnard-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 BHR model includes full transport equations for the Reynold stresses, turbulent mass-flux velocity, density fluctuations, and the turbulent kinetic energy dissipation rate.
3. Multi-Fluid model:
These models 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.
The first step for the development of turbulence models is to perform Reynold averaging of the governing equations and resultant terms.
3.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 are the buoyancy-drag models for the self-similar growth of RT instabilities.
3.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.
3.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.
Modeling Approach:
Pre-Processing:
Step 1: Preparation of geometry
The geometry is essentially two squares , one on top of another of side 20 mm. The top square has been named 'water' and bottom square 'air' . The common edge between the domains was shared using share tool in Workbench. The created geometry is depicted below:
Step 2: Meshing
2.I. Baseline mesh
Mesh properties:
Element size = 0.5 mm
Statistics:
The mesh created is depicted below:
Mesh Quality:
1. Element Quality:
2. Aspect Ratio:
3. Skewness:
4. Orthogonal Quality:
2.II. Refined mesh
Mesh properties:
Element size = 0.2 mm
Statistics:
The mesh created is depicted below:
Mesh Quality:
1. Element Quality:
2. Aspect Ratio:
3. Skewness:
4. Orthogonal Quality:
Step 3: Solving using Fluent
Solver - Pressure based , transient solver with Absolute velocity formulation with Gravity of value 9.81 `m/s^2` enabled
Viscous model - Laminar
Material properties of air, water and user-defined material:
Multi-phase model properties:
Solution Methods:
After Standard Initialization ,we have to patch both the domains to the respective materials, wherein we provide a volume-fraction value of 1 to water and 0 to air and user-defined material.
After solving for 1000 Time steps with 5e-3 dt (time step size), the following solution/animation was obtained:
Case 1: Baseline mesh(Air and water):
Case 2:Refined mesh(Air and water):
Case 3:Refined mesh(Water and user-defined material):
Observation:
Finer mesh captures the instability better than the coarse mesh , but it takes more time for the heavier fluid to settle down. So it is computationally expensive but produces results with better accuracy.
Since the simulation involves capturing instabilities which are changing in an unpredictable manner as a function of time , transient approach is preferred over steady state approach.
Atwood Number
It is a dimensionless number that is used to study hydrodynamic instabilities in density stratified flows.
`A=(ρ_1-ρ_2)/(ρ_1+ρ_2)`
Where,
`ρ_1` is density of heavier fluid
`ρ_2` is density of lighter fluid
For Air-water RT instability,
`ρ_1` = 998.2 `(kg)/m^3` [Water density]
`ρ_2` = 1.225 `(kg)/m^3` [Air density]
A = (998.2-1.225)/(998.2+1.225)
A = 0.998
For water- UD material RT instability,
`ρ_1` = 998.2 `(kg)/m^3` [Water density]
`ρ_2` = 1.225 `(kg)/m^3` [UD material density]
A = (998.2-400)/(998.2+400)
A = 0.428
Behaviour of RT Instability on variation of Atwood number:
RTI evolves in three distinct stages: Linear stability , mushroom head(falling) or bubble(rising) formation , long term evolution due to bubble merging and mixing. For low Atwood number, bubble and mushroom head shape is relatively symmetrical, the disturbances or perturbations are linear for a larger time. In contrast, for large Atwood number, bubble and mushroom head shape is less symmetrical, also it occurs very rapidly at an accelerating rate which leads to a non linear growth rate.
Conclusion:
Atwood number forms a basis to compare and study research papers based on RT Instabilty. Thus we can use these papers to validate our results for beetter accuracy.