Abstract
In the geophysical context, there are a wide variety of mechanisms which may lead to the formation of unstable density stratification, leading in turn to the development of the Rayleigh-Taylor instability and, more generally, interfacial gravity-driven instabilities, which involves moving boundaries and interfaces. The purpose of this work is to study the level set method and to apply the process to study the Rayleigh-Taylor instability experimentally and numerically. With the help of a simple, inexpensive experimental arrangement, the R-T instability has been visualized with moderate accuracy for real fluids. The same physical phenomenon has been investigated numerically to track the interface of two fluids of different densities to observe the gravitational instability with the application of level set method coupled with volume of fraction replacing the Heaviside function. Good agreement between theory and experimental results was found and growth of instability for both of the methods has been plotted.
1. Introduction
The Rayleigh-Taylor instability is instability of an interface of two fluids of different densities which occurs when the interface between the two fluids is subjected to a normal pressure gradient with direction such that the pressure is higher in the light fluid than in the dense fluid. This is the case with an interstellar cloud and shock system. A similar situation occurs when gravity is acting on two fluids of different density—with the denser fluid above a fluid of lesser density—such as water balancing on light oil. Considering two completely plane-parallel layers of immiscible fluid, the heavier on top of the light one and both subject to the Earth’s gravity, the equilibrium here is unstable to certain perturbations or disturbances. An unstable disturbance will grow and direct to a release of potential energy, as the heavier material moves down under the gravitational field and the lighter material is displaced upwards. Such instability can be observed in many situations including technological applications as laser implosion of deuterium-tritium fusion targets, electromagnetic implosion of a metal liner and natural phenomena as overturn of the outer portion of the collapsed core of a massive star, and the formation of high luminosity twin-exhaust jets in rotating gas clouds in an external gravitational potential.
Various numerical and experimental works have been done by many researchers concentrating on the growth of single wavelength perturbations as well as considering different wavelength modes. Sharp [1] presented some of the critical issues concerning Rayleigh-Taylor instability. The importance to carry out the three-dimensional study of Taylor instability and the role of statistically distributed heterogeneities on the growth of instability have been analyzed in his work. Read [2] experimentally investigated the turbulent mixing by Rayleigh-Taylor instability and the results showed that if the instability arises from small random perturbations, the width of the mixed region grows in proportion to . The same investigation has been done numerically by Youngs [3] to simulate the growth of perturbations at an interface between two fluids of different density. If the mixing process evolves from small perturbations then the growth of instability is controlled by the non-linear interaction between bubbles of different sizes. Dalziel [4] investigated the Rayleigh-Taylor instability experimentally using a simple apparatus of novel design where the initial nonlinear perturbations to the flow have been introduced by the removal of the barrier separating the two fluid layers and a good agreement between the results of this work and a previous one has been achieved. Velocity measurements have been done by particle tracking using the method of particle image velocimetry. Voropayev et al. [5] experimentally analyzed the evolution of gravitational instability of an overturned, initially stable stratified fluid. In the present analysis, the instability is initiated by overturning the experimental setup such that the heavy fluid lies over the lighter one. The present study is mainly concerned about the propagating interface between the two fluids and its formation and growth rate. A propagating interface is a closed surface in some space that is moving under a function of local, global, and independent properties. A variety of numerical algorithms are available to track propagating interfaces, and in the present numerical simulation level set method coupled with volume of fraction has been used. Level set method is a computational technique for tracking moving interfaces which rely on an implicit representation of the interface whose equation of motion is numerically approximated using schemes built from those for hyperbolic-conservation laws. The consequential techniques are able to handle problems in which the speed of the evolving interface may sensitively depend on local properties such as curvature and normal direction, as well as complex physics of the front and internal jump and boundary conditions determined by the interface location.
The volume of fluid (VOF) technique has been presented by Hirt and Nichols [6] as a simple and efficient means for numerically handling free boundaries in a calculation mesh of Eulerian or arbitrary Lagrangian-Eulerian cells. It works extremely well for a wide range of complicated problems and this process is very much conservative in nature, but the appropriate tracking of the interface is not possible by this method.
Sethian [7] presented a case of the evolution of a front propagating along its normal vector field with curvature-dependent speed. Numerical methods based on finite difference schemes for marker particles along the front are shown to be unstable in regions where the curvature builds rapidly. And then the front tracking based on volume of fluid techniques has been used together with the entropy condition.
Various numerical methods were developed to study the propagating interfaces. Osher and Sethian [8] devised new numerical algorithms, called PSC algorithms, for fronts propagating with curvature-dependent speed. Merriman et al. [9] extended the Hamilton Jacobi formulation of Osher and Sethian and proposed a level set method for the motion of multiple junctions where the diffusion equation was shown to generate curvature-dependent motion. Zhu and Sethian [10] considered hydrodynamic problems with cold flame propagation by merging a second-order projection method for viscous Navier stokes equations with modern techniques for computing the motion of interfaces propagating with curvature-dependent speed. A new method was presented by Unverdi and Tryggvasan [11] to simulate unsteady multifluid flows in which a sharp interface or a front separates incompressible fluids of different density and viscosity. Chopp and Sethian [12] studied hyper surfaces moving under flow that depends on the mean curvature. The approach was based on a numerical technique that embeds the evolving hypersurface as the zero Level Set of a family of evolving surfaces. Sussman et al. [13] combined a level set method with a variable density projection method for capturing the interface between two fluids to allow for computation of two-phase flow where the interface can merge or break considering a high Reynolds number. Chang et al. [14] presented a level set formulation for incompressible, immiscible multi fluid flow separated by a free surface and the interface was identified as the zero Level Set of a smooth function.
Theory and algorithms of level set method were reviewed by Sethian [15] for the evaluation of the complex interfaces. Topological changes, corner and cusp development, and accurate determination of geometric properties such as curvature and normal direction were obtained by the method. Few years later, Sethian [16] summarized the development and interconnection between narrow band level set method and fast marching method, which provides efficient techniques for tracking moving fronts. At another paper, Sethian [17] reviewed past works on fast marching method and level set method for tracking propagating interfaces in two or three space dimensions.
Kaliakatos and Tsangaris [18] studied the motion of deformable drops in pipes and channels using a level set approach in order to capture the interface of two fluids. The shape of the drop, the velocity field, and the additional pressure loss due to the presence of the drop, the relative size of the drop to the size of the pipe or channel cross-section, the ratio of the drop viscosity to the viscosity of the suspending fluid, and the relative magnitude of viscous forces to the surface tension forces were computed. Son and Hur [19] combined a level set method with the volume of fluid method to calculate an interfacial curvature accurately as well as to achieve mass conservation. They developed a complete and efficient interface reconstruction algorithm which was based on the explicit relationship between the interface configuration and the fluid volume function.
Sethian and Smereka [20] provided an overview of level set methods, introduced by Osher and Sethian [8], for computing the solution to fluid-interface problems. They discussed the essential ideas behind the computational techniques that rely on an implicit formulation of the interface and the coupling of these techniques to finite-difference methods for incompressible and compressible flow. Majumder and Chakraborty [21] developed a novel physically based mass conservation model in the skeleton of a level set method, as a substitute to the Heaviside function based formulation. The transient evolution of a rounded bubble in a developing shear flow and rising bubbles in a static fluid, the Cox angle, and the deformation parameter characterizing the bubble evolution were critically examined. Carlès et al. [22] used a magnetic field gradient to draw down a low density paramagnetic fluid below a more dense fluid in a Hele-Shaw cell. An extended level set method for classical shape and topology optimization was proposed based on the popular radial basis functions by Wang et al. [23]. The implicit level set function was approximated by using the RBF implicit modeling with multiquadric splines. Sun and Tao [24] presented a coupled volume of fluid and level set (VOSET) method for computing incompressible two-phase flows. VOF method was used to conserve the mass and level set method was used to get the accuracy of curvature and smoothness of discontinuous physical quantities near interfaces. Sussman and Puckett [25] presented a coupled level set/volume-of-fluid (CLSVOF) method for computing 3D and axisymmetric incompressible two-phase flows and Sussman [26] presented a coupled level set and volume of fluid method for computing growth and collapse of vapor bubbles. A level set method was combined with the volume of fluid method by Son [27] for computing incompressible two-phase flows in three dimensions where the interface configurations were much more diverse and complicated. A passive scalar transport model has been studied by Wang et al. [28] to study the 3D Rayleigh-Taylor instability. The characteristic behavior and the principle of the interfacial motion from both sinusoidal and random perturbations have been achieved. Youngs [29] numerically simulated three-dimensional turbulent mixing of miscible fluids of RT instability which concluded that significant dissipation of turbulent fluctuations and kinetic energy occurs via the cascade to high wave numbers. The chaotic stage of Rayleigh-Taylor instability is characterized by the evolution of bubbles of the light fluid and spikes of the heavy fluid. Gardner et al. [30] proposed a statistical model to analyze the growth of bubbles in a Rayleigh-Taylor unstable interface. The model using numerical solutions based on the front tracking method has been compared to the solutions of the full Euler equations for compressible two-phase flow. Later, Glimm et al. [31] numerically studied the dynamics of the bubbles in chaotic environment and their interactions with each other as well as the acceleration of the bubble envelope.
The Rayleigh-Taylor instability is a gravity driven instability of a contact surface and this growth of this instability is sensitive to numerical or physical mass diffusion. Li et al. [32] addressed this problem using a second-order TVD finite difference scheme with artificial compression. They numerically simulated the 3D Rayleigh-Taylor instability using this scheme. A new model was proposed by Chen et al. [33, 34] for the momentum coupling between the two phases. The Rayleigh-Taylor instability of an interface separating fluids of distinct density is driven by acceleration across the interface. Two-phase turbulent mixing data were analyzed, which have been obtained from direct numerical simulation of the two-fluid Euler equations by the front tracking method. Direct numerical simulation of three-dimensional Rayleigh-Taylor instability (RTI) between two incompressible, miscible fluids has been presented by Cook and Dimotakis [35]. Mixing was found to be even more sensitive to initial conditions than growth rates. The flow structure and energy budget for Rayleigh-Taylor instability using the results of a high resolution direct numerical simulation have been examined by Cook and Zhou [36]. Later Cook et al. [37] described large eddy simulation for computing RT instability. A relation has been obtained between the rate of growth of the mixing layer and the net mass flux through the plane associated with the initial location of the interface.
In this present work, the level set methodology has been applied to visualize theoretically the RT instability using a triangular distribution of initial disturbance. The fraction of volume in the interface control volumes has been successfully incorporated for identifying the interface very accurately. The topological changes with time have been captured accurately and this has been matched effectively with the experimental results. The instability growth rate which is predicted by the theory is confirmed by the experimentation with the initial incipience of linear distribution of disturbances as already stated. This is a positive contribution along with the theoretical topological visualization of the RT effects.
On the other hand, the merging and consequent breaking up of the interfaces has been captured while the RT instability growth takes place. These results are important as they provides the probable trapping, merging, and consequent breaking of the oil and natural gas pools trapped between the formation of salt domes and overlying sedimentary rocks. These effects of the geothermal RT instabilities and deformation of the rocks above the salt domes are important as they provide the possibility of exploration of oil and gas pools, thus coagulated and subsequent fragmented in huge mass under the earth for million of years. These results are encouraging and can bridge our knowledge of RT to apply to the oil and gas industry.
2. Geometrical Description
The geometry of the problem is shown in Figure 1. A fluid layer with a thickness and density overlies a second layer of thickness and density . The upper boundary and lower boundary are assumed to be rigid surfaces. Here, is greater than . The undisturbed interface between fluid layers is taken to be at . Due to gravitational instability, the interface between the fluids distorts and motions occurs in the fluid layers. The displacement of the disturbed fluid layers is denoted by . When the heavy fluid lies above the light fluid, the configuration becomes unstable. The time to grow the instability depends on the viscosity of the fluid and the density difference of the fluids. When the viscosity of the fluid is high and the density difference is smaller, the instability takes longer to grow.
3. Experimental Setup
The experimental setup consists of a closed rectangular box made of Perspex of 20.4 cm × 10.2 cm × 15 cm dimension. There are two openings at the top surface with valve arrangement for the purpose of filling the box with the required liquids. The two side handles are provided for convenience turning of the setup to upside down or vice versa in quick time. The setup is placed on a preleveled surface and lower half of the box is filled with glucose solution and upper half is filled with colored refined soya bean oil, with the help of funnels. The viscosities of both the liquids were measured in the laboratory at room temperature by “Falling Sphere method” and density of the fluids was measured by simply measuring their mass and volume (see Figure 2).
The viscosity and specific gravity of the liquids have been measured as follows: Viscosity of glucose syrup = 350 Pa-S, Viscosity of oil = 0.0791 Pa-S, Specific Gravity of Glucose syrup = 1.4, Specific Gravity of oil = 0.92.
4. Experimental Technique
In the experiment, first the setup rests at position 3 where the light fluid lies over the heavy one. In this position, it is totally balanced and stable. Then the setup is turned upside down quickly so that heavy liquid lies in the upper half and thus instability is initiated. The instability can also be initiated by keeping the setup at position 2 where the heavy and light liquids stand vertically side by side in an unbalanced and unstable condition. Naturally all these configurations want to return to position 3 to minimize the potential energy and to gain a stable and balanced position. The whole process is captured to track the moving interface and to study the growth rate of the instability with time (see Figure 3).
5. Formulation of Two-Phase Flow with Surface Tension
The term two-phase flow refers to the motion of two different interacting fluids or with fluids that are in different phases. In the present analysis, only two immiscible incompressible fluids have been considered and a low enough Reynolds number is assumed so that the flow can be considered as laminar flow. Level set method may be applied to track the interface efficiently in case of incompressible, immiscible fluids in which steep gradient in viscosity and density existed across the interface. In these problems, the role of surface tension is crucial and formed an important part of the algorithm.
6. Numerical Modeling
For mathematical analysis, we assume a system of two-fluid phases constituting a two-dimensional domain. The individual fluid phases are assumed to be incompressible but deformable in shape on account of shear stresses prevailing between various fluid layers as well as fluid-solid interfaces. We assume the flow field to be two dimensional and laminar.
Navier-Stokes equation is given as Assume a sharp fluid interface between two fluids with different densities, and also the flow is incompressible, and thus
The surface tension term acts normal to the fluid interface and is proportional to the curvature, due to balance of force argument between the pressure on each side of the interface. This leads to the relation, Thus, surface tension acts as an additionally forcing term in the direction normal to the fluid interface.
Now replacing normal by and when distance is approximated by , we have The curvature can be expressed by and its derivatives as follows, As in [14], regularized delta function can be defined as This recasts the surface tension in the level set framework. If is always reinitialized to the distance function, the Dirac delta function itself can be smoothed.
Thus the equation of motion become The governing equations can be written as follows.
6.1. Continuity
Consider
6.2. Momentum
Consider
A scalar variable, level set function is used to identify the interface between two fluids and also acts as a distance function. The equation transporting the interface can be written as where is the level set function prescribing position of the interface at any specified time instant. If the value of the at the interface is taken as zero, it effectively becomes a distance function satisfying
But at all instant of times must remain a distance function, to ensure that another scalar variable needs to be introduced and solved. This variable must be constrained to constitute a distance function having the same interface value as . This can be achieved by obtaining a pseudo-steady-state solution for the following transient transport equation of : where with being a pseudo-time step.
Equation (12) is subjected to the following initial condition
The reinitialization process is iteration of (12) with a pseudo-time step, and within a few iterations it comes to a steady state solution. Then the reinitialization procedure ends leading to reassignment of the level set value.
It is evident that pseudo-steady-state value of is the value of at the time instant . Success of the mass correction is affected by (12) which depends on the accuracy of the interpolation of physical properties such as density, and viscosity across the interface. This can be achieved by calculating a property within a control volume as where is called Heaviside function.
The equation for the one-dimensional volume fraction is given by
and for two-dimensional volume fraction the concept has been taken from [21].
At the solid boundary, the Neumann boundary condition for the level set function has been utilized.
7. Solution Procedure
The governing differential equations, coupled with appropriate boundary conditions, are solved using a pressure based finite volume method, as per the SIMPLER algorithm [38]. Convection-diffusion terms in the conservation equations are discretized using the power law scheme [38].
The location of the interface at time has been specified and then the normal distance for all nodes from the interface is calculated. The properties at all nodes have been specified using (15). The continuity and momentum conservation equations at time instant are solved. Then using the velocities obtained in a previous step, using (10), has been solved. Next, using the values of from preceding step as initial values, the pseudo-steady-state (12) has been solved. Setting , the procedure is going on until the desired convergence is achieved.
The temporal term of the momentum equation has been discretized as follows. Equation (9) in two-dimensional form is discretized to get algebraic linear simultaneous equations as follows: where represents general variables and and are the total (convection plus diffusion) fluxes defined by where and denote the velocity components in the and directions, is the source term, and represents the diffusion coefficient. The integration of (17) over the control volume (Figure 4(a)) gives
(a)
(b)
The source term is linearized in the usual manner anticipating negative slope while the unsteady terms and are assumed to prevail over the whole control volume. In a similar fashion, the continuity equation is also linearized.
Similarly the pressure gradient term is discretized considering the staggered control volume as: where .
This is for the equation as shown in Figure 4(b). The corresponding other equations are discretized in a similar fashion. Finally, guessing the velocity field, the pressure equation is solved, and consequently by correcting the velocity field the variables are solved. This method has an essence physically possible solution by removing unrealistic checker board results.
8. Results and Discussions
8.1. Gravitational Instability due to Density Difference with Initially Horizontal Layers of Fluids
If the box is rotated in the plane quickly so that the heavy liquid occupies the upper portion, then instability will initiate at once in the presence of sufficient unavoidable perturbations and the interface starts moving. The position of the interface at different times, especially in the initial stage of growing instability, has been analyzed in the present investigation (see Figure 5).
8.2. Theoretical Growth of Instability
From theoretical analysis of the problem, the growth rate is given by
Here it has been assumed that the flow is laminar owing to the fact that the viscosities are of very high order in a two-dimensional, incompressible flow.
The solution of this equation is where is the initial displacement of that point of the interface from the undisturbed interface and is the growth time.
Now, it can be seen from the growth equation that growth rate varies linearly with displacement of that point at a particular time.
Now, for comparison purpose, a point at a distance of 6.1 cm, that is, approximately distance from the left vertical wall, is considered, where is the wavelength of the applied perturbation sine curve which in our case is the length of the box = 20.3 cm.
The displacement of the considered point is measured at different times from the undisturbed interface by proper measurement in the series of snapshots presented in Figure 11 and the graph between growth rate and time has been plotted (see Figure 6).
It can be seen that the best fitted curve is unbounded exponential in nature, which agrees very much with the theory demanding exponential growth of the instability. The curve is of the form , or we can write where growth time is 9.09 seconds. We also see that the initial perturbation of the observed point is 0.468 cm.
Theoretically, the value of can be calculated as where is equivalent coefficient of dynamic viscosity and it can be expressed as
where = density of lighter liquid = 0.92 × 103 kg/m3, = density of heavier liquid = 1.4 × 103 kg/m3, = viscosity of lighter liquid = 0.0791 Pa-s, and = viscosity of Heavier liquid = 350 Pa-s. So becomes 211.23 Pa-s.
Here, = height of the upper or lower rigid boundary from the undisturbed, liquid interface = 7.5 cm = 0.075 m, = acceleration due to gravity = 9.8 m/s2, and = wave length of the perturbation sine curve = 20.3 cm. So, becomes 7.846 seconds.
Now, from the experimental study, the initial perturbation of the specified point is 0.468 cm.
So the theoretical growth equation becomes, —theoretical growth equation, —experimental growth equation.
It can be observed from the above two expressions that the characteristics of the development of growth of instability are quite similar, with a slight difference in the growth time. The growth time is slightly higher in case of experimental observation than the numerical investigation.
Figure 7 shows the growth of instability with time and Figure 8 shows the comparison between the theoretical and experimental results. Figure 9 depicts the theoretical and experimental comparisons of growth rate of instability. From the figures, it is observed that, at the early stage of growth of instability, the experimental and theoretical results matche considerably while the growth rate differs with the increase in time. This may be due to the fact that theoretically the flow has been assumed to be two dimensional, but in case of experimentation the three-dimensional characteristics come into consideration and due to this effect of three dimensionality the theoretical results differ with the experimental result and it increases with increase of time.
The same problem is numerically analyzed considering a rectangular two-dimensional domain. Two arrays of and grid points in axial and radial directions, respectively, have been used. It has been observed that the grid independent study has shown 0.001% change and the results are almost unaffected considering both the grid meshes. The grid array of has been used for all subsequent results reported here with uniform mesh size and time step s. A disturbance of the vertical component of velocity having triangular distribution has been introduced as shown in Figure 10.
The variation of the propagating interface with time has been shown. Both the experimental and numerical results are presented here (see Figure 11). In the numerical results, red color represents the lighter fluid and blue color represents the heavy fluid, whereas in the experimental results white solution is the heavy fluid and the red colored fluid is the light fluid.
It can be observed from the above figures that the experimental results are in good agreement with the numerical results. The interface between the two fluids shows similar pattern during the study for both experimental and numerical analysis. However, three-dimensional features are observed to affect the results as seen in the experimental study. Figure 12 shows the comparison of growth rate of instability. The numerical result matches the linear theory at the initial period, whereas with the increase of time it varies with the linear theory.
In Figure 13, two-bubble merging and consequent breaking up have been evaluated with time. The matrix is the heavier fluid. This is the result of numerical experimentation.
9. Conclusions
The nature of the development of instability was experimentally found as a function of sine curve as predicted by theoretical model. A numerical methodology was devised and validated with experimental results so that the methodology can handle any gravitational interfacial instability. It was found that, in the early stages of the growth of instability, the growth rate is proportional to the instantaneous growth in a particular position, that is, growth rate varies linearly with growth at that moment at a particular point on the interface. But, at the later stage of development of instability, substantial deviation from the linear theory was observed. The pictorial views of the interface between the two fluids have been studied both theoretically and experimentally and they have matched satisfactorily.
Nomenclature
| : | Atwood number (−) |
| : | Height of one fluid layer (cm) |
| : | Speed function (m/s) |
| : | Acceleration due to gravity (m/s2) |
| : | Heaviside function (−) |
| : | Pressure (N/m2) |
| : | Time (s) |
| : | Velocity (m/s) |
| : | Growth at a particular time (cm) |
| : | Initial growth (cm) |
| : | Small time step (s). |
| : | Level set function (m) |
| : | Dummy variable for level set function (m) |
| : | Dirac delta function (−) |
| : | Fluid properties such as density and viscosity (−) |
| : | Surface tension coefficient (N/m) |
| : | Coefficient of viscosity (Pa-s) |
| : | Curvature (m−1) |
| : | Density of the liquid (kg/m3) |
| : | Growth time (s) |
| : | Wavelength of the perturbation (cm). |
Acknowledgment
This work is supported by the Council of Scientific and Industrial Research (CSIR), Government of India.