Numerical Modeling of the Breakup Mechanism of Liquid Ligaments in Two-Phase Flow Application Using Level Set Method
Mechanical Engineering Dept., Faculty of Engineering, Taif University, Al-Haweiah, Taif, Saudi Arabia
The present paper introduces a numerical modeling for one of the most important problems of two-phase flows; namely the breakup mechanism of stretched liquid ligaments. The challenge problem of such complex flow is how to deal with the breakup of liquid ligaments in natural way without any constraints or numerical instability. The present modeling is based on solving the governing equations of motion in the liquid phase by the aid of the control volume strategy on a non-sta-ggered grid system. The level set function is further applied to track and capture the moving interface by using the normal velocity components located directly at the liquid interface. Some numerical test cases are performed to indicate the ability of the proposed numerical method in predicting the breakup mechanism of liquid ligaments. The obtained results showed that the developed numerical method is capable of simulating the formation and the breakup mechanism of the ligaments formed in two-phase flows without any numerical constraints.
Breakup Mechanism, Liquid Ligaments, Level Set Method, Numerical Modeling, Two-Phase Flow Dynamics
Received: May 20, 2015
Accepted: June 6, 2015
Published online: July 15, 2015
@ 2015 The Authors. Published by American Institute of Science. This Open Access article is under the CC BY-NC license. http://creativecommons.org/licenses/by-nc/4.0/
The separation and breakup of liquid ligaments are considered as the most challenging topic of the numerical simulation of two-phase flows . That can be seen in numerous industrial and engineering applications, e.g., atomization and spray of liquid jet , bubble and droplet dynamics , and other multiphase flow systems .
Numerical simulations have recently become an important tool for predicting and understanding many engineering processes. The mean feature of such engineering processes is the moving interface between the two phases. This moving interface can be advected by the normal velocity components on the separating surface according to the existed interfacial boundary conditions. The advection process may be a linear or a nonlinear regime according to the included terms in the governing equations. In the nonlinear regimes, there is a few numbers of developed methods for solving the governing equations [5,6]
The most general methods for treating two-phase flows with moving interfaces are boundary integral formulations  and further applications in . Many incompressible flow algorithms have been used to solve the two-phase flow problems accompanied by a suitable surface tracking method. As mentioned in , both of conventional conservative and high order conservative schemes can produce excessive numerical diffusion and oscillations which destroy the smoothness of the front. That can be seen in . Another scheme is based on the projection method developed in , and it has been combined along with the level set to simulate two-phase flows, e.g. . In conclusion, a truly general method is not yet available .
A great challenge in the numerical simulations of unsteady motion of two-phase flows is how to predict the topological change of the interface because the motion is very sensitive to small numerical perturbations, such as merging and breaking of the interface. Many techniques have been used to deal with the propagation of the interface in two-phase flows. These techniques are based basically on two different approaches: the Lagrangian approach and the Eulerian approach see for more details [14-17].
Based on the Eulerian approach, another technique has been used in . This technique is based on tracking the interface using a grid that moves through a stationary Eulerian grid. Here complication occurs when one needs to add or subtract points to the moving grid. The major limitation of this method is the treatment of interfaces that interact. Moreover, these problems are become important when solving a three dimensional problem. Another example for moving grid technique is used in . In this technique the fluid is bounded by a height function varying with time along the fluid interface. This function determines the radial distance through evolution. This algorithm is limited to a case in which the interface does not overturn, see for more details .
Two types of methods are available and widely used; the Volume-Of-Fluid (VOF) methods [20,21] and the Level Set methods . These methods avoid explicitly tracking the interface, and consequently, avoiding many of the topological limitations.
The level set method, unlike the VOF method, divides the solution domain containing the two phases into a region where a scalar quantity, say, is smaller or larger than a fixed value f0. The value of the level set function f is assigned to each grid point in the calculation domain, thereby determining whether it belongs to the liquid or the gas phase. The level set method is widely applied in two-phase flows applications, e.g. [25,12] and more recently in [26-33].
2. Mathematical Formulations
The governing equations for incompressible two-phase flow are mathematically expressed by the conservation equations of the mass and momentum at each point of the flow field. These equations can be written in the primitive variables formulation in form of continuity equation and Navier-Stokes equations, respectively, as follows:
Where r, u, p, m are the density, velocity vector, pressure and viscosity of the fluid. The subscript i denotes the liquid phase (i=1) and the gas phase (i=2), respectively
The coupling between the pressure and velocity field is made through the Poisson equation for pressure:
The Poisson equation is solved by the Successive Over-Relaxation method in each phase to obtain the pressure in the considered phase.
2.1. Computational Methodology
The so-called implicit fractional step-non iterative method is applied in our present work by presuming that the velocity field reaches its final value in two stages;
Whereby, is an imperfect velocity field based on a guessed pressure field, and is the correction velocity. The 'starred' velocity will be obtained from solving of the momentum equations. After that, the solution of Poisson equation for the pressure:
where will be called the pressure correction.
The velocity correction is obtained according to:
This fractional step method described above ensures the proper velocity-pressure coupling for incompressible flow fields .
2.2. Level Set Method
In absence of the interfacial mass transfer such as evaporation or condensation, the equation of the level set method can be written as:
And from the level set function, the normal vector and the interface curvature can be calculated as:
It should be pointed out that, the curvature of the interface is calculated only on the main grid points. However, the exact curvature should be calculated at the interfacial markers located on the intersection points of the computational grid with the interface. Consequently, the interpolation scheme applied for obtaining the fluid variables at the interface is also applied for calculating the interfacial curvature, for more details one can see 
In another form of the level set method, where the normal velocity component is included, the level set equation becomes
By further manipulation, the normal velocity is replaced by the velocity field Fext known as the extension velocity, which at the zero level set, equals the given speed Vn. In other words, the level set equation can be written as:
The above equation is solved by using the second-order Runge-Kutta method and the overall domain is re-initialized each time step .
2.3. Interfacial Boundary Conditions
The boundary conditions at the interface, or jump conditions, may be written in the following general form:
Where I is the identity matrix, p is the pressure, m is the dynamic viscosity, s is the surface tension coefficient, and D is the rate of deformation tensor. The bracket means the jump of the stresses along the fluid interface G, the unit normal vector n is taken from fluid 2 to fluid 1,denotes the gradient in the local free surface coordinates, and t is an arbitrary vector perpendicular to the normal of the interface. The curvature of the interface k is computed from the following equation:
In contrast to the previous two-phase numerical methods, in which the interfacial jump conditions are embedded naturally on the momentum equations by applying CSF model , the jump conditions between the two fluids are treated here as a boundary conditions enforced explicitly at the interface G.
The above interfacial boundary equation, Eq. (12) can be written in the normal and tangential directions, respectively:
The kinematic conditions are satisfied by assuming the continuity of the normal velocity components, i.e.
In addition to that, the liquid gas interface requires equality of dynamic stresses, especially shear stresses, on both sides of the interface, i.e.
The system of the above equations and the prescribed boundary conditions should be solved simultaneously to determine the flow field in the two fluids. The normal velocity components are then used for the advection of the interface by solving the appropriate equation of the level set function.
3. Results and Discussions
The first test case performed using the above developed numerical method is the falling of a water droplet from a pipette with specified dimensions and with a constant initial velocity. The sequence of the pictures presented shows that the droplet is firstly formed from a nearly sinusoidal wave form initiated at the open side of the pipette. Under the effect of the initial velocity, the droplet is formed and a thin ligament is shown to connect the droplet with the liquid inside the pipette. At the last stages of the simulation the droplet detaches from the liquid surface and a single droplet moves under the gravity downward.
From the above results, it can be shown that the formation of the ligament and the breakup of the single droplet are carried out in natural way without any numerical instabilities or constraints.
The second performed test case shows the breakup mechanism of an elongated droplet with initial ellipsoid shape under the surface tension effects only. The deformation of the droplet is started according to the nonlinear theory of droplet deformation. At later steps of the simulation, two satellite droplets and a mother droplet can be observed. As in the previous example, the breakup mechanism is carried out without any numerical problems.
In the present paper, the numerical modeling of the formation of liquid ligament and the further breakup of the liquid droplet is carried out. The numerical method developed is based on solving the momentum equations coupled with the level set method. The obtained results showed a remarkable capability of the developed numerical method in predicting the formation and the further breakup of the ligaments and satellite droplets without any numerical constraints. This encourages us to further include a wide range of two-phase flow industrial and engineering applications.1