Numerical modeling of turbulence-induced interfacial instability in two-phase flow with moving interface

Applied Mathematical Modelling 36 (2012) 3593–3611

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Numerical modeling of turbulence-induced interfacial instability in two-phase flow with moving interface Ashraf Balabel ⇑ Mechanical Power Engineering Dept., Faculty of Engineering, Menoufia University, Shebin Elkom, Menoufia, Egypt

a r t i c l e

i n f o

Article history: Received 2 May 2011 Received in revised form 20 October 2011 Accepted 2 November 2011 Available online 9 November 2011 Keywords: Interfacial instability Level set method Nonlinear turbulence modeling Numerical methods Two-phase flow

a b s t r a c t The present paper introduces a new interfacial marker-level set method (IMLS) which is coupled with the Reynolds averaged Navier–Stokes (RANS) equations to predict the turbulence-induced interfacial instability of two-phase flow with moving interface. The governing RANS equations for time-dependent, axisymmetric and incompressible two-phase flow are described in both phases and solved separately using the control volume approach on structured cell-centered collocated grids. The transition from one phase to another is performed through a consistent balance of kinematic and dynamic conditions on the interface separating the two phases. The topological changes of the interface are predicted by applying the level set approach. By fitting a number of interfacial markers on the intersection points of the computational grids with the interface, the interfacial stresses and consequently, the interfacial driving forces are easily estimated. Moreover, the normal interface velocity, calculated at the interfacial markers positions, can be extended to the higher dimensional level set function and used for the interface advection process. The performance of linear and non-linear two-equation k–e turbulence models is investigated in the context of the considered two-phase flow impinging problem, where a turbulent gas jet impinging on a free liquid surface. The numerical results obtained are evaluated through the comparison with the available experimental and analytical data. The nonlinear turbulence model showed superiority in predicting the interface deformation resulting from turbulent normal stresses. However, both linear and nonlinear turbulence models showed a similar behavior in predicting the interface deformation due to turbulent tangential stresses. In general, the developed IMLS numerical method showed a remarkable capability in predicting the dynamics of the considered two-phase immiscible flow problems and therefore it can be applied to quite a number of interface stability problems. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Turbulent two-phase flow can be seen in various industrial and engineering applications such as liquid atomization [1], internal engine fuel spraying [2] and droplet dispersion in turbulent gas–liquid flow [3], inkjet printing [4] and flow coating [5]. One of the important issues in studying turbulent two-phase flow is the stability and the topological evolution of the moving interface between the two phases. Great attention should also be directed to the nonlinear behavior in the vicinity of the singular point where the interface separates. According to the high complexity of the physics encountered in turbulent two-phase flow and the intrinsic features of the interface, no accurate and satisfying analytical models exist up to this date.

⇑ Tel.: +20 48 2311655/0106768952; fax: +20 48 2235695. E-mail address: [email protected] 0307-904X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.11.006

3594

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

Moreover, it is often difficult to use the experimental measurements to visualize and quantify accurately the interface motion in transitional and turbulent flows. Recently, the huge evolution of the numerical methods and computer speed has allowed the accurate prediction of many complex turbulent two-phase flows of engineering and industrial relevance. However, the numerical study of such complex flows has been somewhat limited because it poses several great challenges [6,7]. The numerical simulations of the complex turbulent two-phase flow with moving interfaces need both refined turbulence models and accurate numerical methods with higher order numerical discretization schemes, usually coupled to a high accurate and robust tracking method in order to follow the complex topological changes of the interface with accurate prediction of the interface normal and curvature. The main objective of most previous numerical methods concerning the complex turbulent two-phase flows was to perform reliable and accurate simulations in order to investigate ligament formation, stretching and fragmentation into small scale droplets and to further study the influence of the Reynolds and the Weber numbers on the spray-formation dynamic process [7,8]. However, such computations in the presence of large density and viscosity ratio and the singular nature of the interfacial forces provide unique challenges, so that it is not surprising that numerical simulations of two-phase flow in the nonlinear regimes are few and require some ad-hoc prescriptions for continuation through singularity [9,10]. There are many approaches to turbulence modeling of two-phase flows. Previous studies in turbulent two-phase flows are carried out in the context of direct numerical simulations (DNS). Some examples of numerical studies of the evolution and structure of turbulent shear layers at free surface can be found in [11–13]. Clearly, the advantages of DNS lie in its independence from any modeling or assumptions. However, for high Reynolds number flow, which is the prevailing circumstances in industrial and engineering applications, DNS could not resolve all physical processes encountered and then DNS is beyond possibility. Moreover, the DNS of the atomization process is challenged by involvement of high complexity of physics, surface instabilities, stretching and breakup of ligaments and droplets producing a multi-scale problem that requires high resolution to undertake. To overcome these computational limitations of DNS, the governing equations of motion are subjected to either Reynolds averaging or spatial filtering, resulting in unclosed terms that require closure relations. The governing equation are then closed, however, the choices of such closure models is not straightforward and has some limitations. Recently, the increase in computational resources has made large eddy simulation (LES) an attractive approach for modeling of turbulent flows. In such context, the large scale eddies are simulated, while the small scales called the sub-grid scales (SGS) are modeled [14]. In industrial level, LES is applied in various cases like, internal flows with complex geometries, flows with large separated area, atmospheric boundary layers, simulation in a combustion chamber, and complex flows with moderate Reynolds number. However, LES is extremely computationally demanding for high Reynolds number flows. Moreover, LES is probably realistic away from a solid boundary, however, it has all the limitations of a simple eddy viscosity models near to a boundary. Nevertheless of the previous comments, LES can be directly applied in free-surface flows, provided that the free surface does not exhibit any significant deformation. However, the presence of the complex topological changes in turbulent two-phase flows has relatively a dynamic interaction with both the resolved and modeled turbulence scales in the flow. This can be referred to the interfacial processes that are not described explicitly in the governing equations. Consequently, the modeling of turbulent flow by LES has consequences in the numerical simulation of turbulent two-phase flows. It should be pointed out that, the modeling of turbulent interfacial flow using RANS or LES in the Eulerian formulation is limited. Over the past decade, efforts have focused mainly on the construction of two-equation eddy-viscosity models and Reynolds-stress-transport closure [15]. The two-equation k–e eddy viscosity model still represents a good compromise between accuracy and computational efficiency. Therefore, the two-equation eddy viscosity models have been the subject of much research in the last years as they still the most widely used in industrial and engineering applications, even though they fail to predict correctly a number of complex flows. It was calibrated and validated for many kinds of high-Reynolds number turbulent flows [16]. The conventional eddy-viscosity models are based on a set of linear Boussinesq stress-strain relations. This approach is attractive from a computational point of view, especially in terms of numerical robustness, and it has a large popularity with Computational Fluid Dynamics (CFD) practitioners. However, this approach is known to be affected by major weaknesses which can be considered as the source of substantial errors where the gradient of the normal stresses contribute significantly to the momentum balance, or what is called complex strain. Although some of ad-hoc corrections have been made, usually to the length-scale equation, and/or the formulation of alternative equations for different length-scale parameters, none of these corrections can address the fundamental limitations arising from unrealistic constitutive relations, see for more details [17,18]. The alternative route thus pursued extensively over the past decade has been second-momentum closure which is mathematically complex and numerically challenging and often computationally expensive. Accordingly, it has some limitations in the context of industrial CFD. This has thus motivated efforts to construct models which combine the simplicity of the eddy-viscosity formulations and the superior fundamental strength of second-moment closure. These efforts have given rise to the group of nonlinear eddy-viscosity models. These models have to handle turbulent flow accurately while keeping a high level of accuracy as required by DNS/LES approach. Nonlinear eddy-viscosity models can be traced back to the work of [19]. Such approach is generally found to marginally improve the predicted turbulence intensities. However, relative to the linear models, convergence is mostly difficult to achieve. Most nonlinear eddy-viscosity models utilize constitutive equations which are functions of two turbulence scales (usually k and e) as well as strain and vorticity invariants, e.g. Craft–Launder–Suga (CLS) model [20]. The nonlinear models

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3595

have been previously applied to a wide range of flows in order to assess the ability of these models to predict anisotropy and streamline curvature effects, e.g. plane and fully-developed curved channel flow and complex two-dimensional flows [20,21], and more recently for different jet flow applications [18]. However, the nonlinear models were not widely applied for turbulent two-phase flows. In such situations, the flow exhibits the combined effects of three-dimensionality and streamlines curvature, interfacial normal stresses and surface instabilities. In the computations of two-phase flows, the preservation of the interface shape and the mass conservation properties of individual fluids are of great importance. Recently, a large number of numerical methods have been developed for computing the interface movement and predicting its time-dependent topological changes. The literature review of such methods is quite extensive and comprehensive, therefore only a brief account of some important tracking methods is given. In general, the available numerical method for computing two-phase flows with moving interface can be categorized as interface tracking and interface capturing methods according to the scheme used for identifying the interface. The interface tracking methods involve essentially explicit computational element or Lagrangian particles moving through an Eulerian grid. Such Lagrangian-based methods can predict the interface location precisely. Also, these methods have the ability to resolve features of the interface that are smaller than the grid spacing of the regular Eulerian mesh [22,23]. However, these methods could not deal with the complex topological changes of the interface such as merging and breaking-up which require a specified algorithm in order to add or remove marker particles as they get too far apart or too close together. Moreover, the geometrical quantities of the interface such as normal and curvature are difficult to calculate. These problems are amplified when solving a three-dimensional problem. The interface capturing methods, in contrast to the interface tracking method, implicitly locate the interface through a definition of a separate phase function that discretized on the fixed Eulerian grid. Among several kinds of interface capturing methods, the volume-of-fluid (VOF) [24] and the level set (LS) method [25] are the most popular interface capturing methods. Although the VOF method has been widely applied for predicting different complex two-phase flows, it suffers from several numerical problems such as interface reconstruction algorithms and the difficult calculation of the interface curvature [6]. These numerical problems can, in particular, limit the accuracy and the stability of the numerical method adopted for calculation of two-phase flows, especially when the surface tension is included. A comprehensive review for the different VOF methods and their numerical constraints especially in turbulent flow can be found in [26]. In contrast to the VOF methods, the level set methods offer highly robust and accurate numerical technique for capturing the complex topological changes of moving interfaces under complex motions. The basic idea of LS method is the use of a continuous, scalar and implicit function defined over the whole computational domain with its zero value is located on the interface. Unlike the VOF method, which divided the spatial domain into cells that contain material function, the LS method divided the domain into grid points that contain the value of the scalar function; therefore, there is an entire family of contours. The interface is then describes as a signed distance function at any time, and consequently, the geometric properties of the highly complicated interfaces are calculated directly from LS function. Moreover, the complex topological changes of interfaces such as merging and breaking-up are handled automatically in a quite natural way without any additional procedure. In addition, the extension of the LS method to three-dimensional problems is easy and straightforward. Referring to the previous discussion, the LS methods have seen much use in different CFD-applications of diverse areas, e.g. two-phase flows, turbulent atomization, grid generation and premixed turbulent combustion [27]. However, the LS methods suffer from numerical diffusion which may cause a smoothing out of sharp edges of interface. The LS function is usually evolved by a simple Eulerian scheme, and consequently, the final implementation of LS methods does not provide full volume conservation, so highly accurate transport schemes are required. From the numerical point of view, both of the VOF and LS methods has its own advantages and disadvantages. However, the great attention that turned recently on the development of new capturing-LS-based methods has given the superiority to the LS methods in the computations of two-phase flows with moving interfaces. In more recent comparison between LS and VOF methods [26], it is concluded that when turbulence is taken into consideration, the VOF technique could capture mean effects. However, difficulties are encountered when computing problems with instabilities or when capturing interfaces with strong deformation, e.g. [28]. Therefore, in order to alleviate some of the geometrical problems of VOF method and to improve the mass conservation in the LS method, a large number of coupled LS and VOF methods have been developed, see for more details [29,30]. The resulting approach is completely Eulerian and does not include any characteristics of the front tracking methods. Recently, another class of hybrid methods has been developed and applied for different two-phase flows applications, e.g. the particle level set method [31], the Eulerian level set-vortex sheet method [32], the hybrid Lagrangian–Eulerian particlelevel set method [6], the coupled level set-ghost method [7] and the coupled level set-boundary integral method [33]. All these methods have been quite successful, however, each method has its own advantages and disadvantages. Consequently, it is difficult to affirm which one is generally superior. The success of each method is dependent on the available computational resources and the description of the problem considered. According to the special features of such hybrid methods, the computational cost is much greater and complex than the classical level set method. Moreover, as a result of the different types of the coupled techniques in hybrid methods, there is usually a time step restrictions to obtain stable schemes. Although the ghost method introduces an attractive way for handling the density and pressure jump of two-phase interface, the discretization of the viscous term is challenging and complex [7]. In general, the ghost method results in additional computational cost that becomes very demanding in large three-dimensional simulations.

3596

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

In the present work, we are concerned with further development of standard the level set method proposed in [34], which belongs to the group of Eulerian interface-capturing methods. An exhaustive review of such method can be found in [35]. Although the standard level set formulation provides a good solution to the problem of interface advection in two-phase flow, however, it suffers from various numerical problems. The inherent problem resulted from introducing the transition region and the smoothing of the interface discontinuity in a continuum formulation remains challenging [30], and results in smearing of flow properties and variables, forcing them to be continuous across the interface in spite of the appropriate jump conditions. Even if smooth initial data are considered, the interface can lose smoothness and develop singularities in finite time. As a consequence, the numerical simulation using the standard LS approach cannot be continued past the point of singularity without some sort of artificial regularization being applied. Since the surface tension effects play a significant role in the most important two-phase flow problems such as droplet deformation, bubble motion and liquid ligament breakup, a great attention has been given for the numerical modeling of the surface tension force in two-phase flow simulations. In the context of the standard level set method, the continuum surface force (CSF) model [36] has been widely used to model surface tension. In the CSF model, surface tension effect is treated as a body force or as a source term in the momentum equation. A striking feature of this model, as reported by [22], is the so called spurious currents resulted from the inconsistent modeling of the surface tension force. In some cases, these numerical artifacts may lead to catastrophic instabilities in interface calculations and inaccurate description of steep gradients occurring at the interface. In such cases, the effect of the unbalanced forces acting on the interface can reduce the accuracy of the calculation. Therefore, it is not surprisingly that till now new surface tension models are continuously developed [30]. In order to try alleviating the above numerical problems, in the present paper, new modeling concepts are considered. Firstly, instead of using the continuum formulation of the single phase model [36], the momentum equations are solved in both phases separately. The normal and the tangential stresses as well as the kinematic conditions are satisfied on the interface separating the two phases. The splitting of the two phases presents, from our numerical point of view, several advantages over the standard level set approach. The Navier–Stokes equations are solved in both fluids with constant properties. Consequently, the computation of the pressure can be done in a standard way without pressure and velocity oscillations at the interface that are common in two-phase level set methods with large density ratios. Moreover, the continuity equation is always enforced in both fluids, thus allowing the use of standard collocated methods [35]. The interfacial stresses are modeled by calculating their effects on a number of interfacial markers located on the interface intersections points with the computational grids. Therefore, the presented numerical method is referred to interfacial marker level set method (IMLS). The rest of the paper is organized as follows. In Section 2, the RANS equations for time-dependent, incompressible turbulent two-phase flows are presented, along with the implemented linear and nonlinear turbulence models. Following, the associated interfacial stresses modeling and the interface capturing LS method are explored. Then, in Section 3, the numerical IMLS method applied to solve the complete set of equations is described in brief. In Section 4, the developed IMLS is tested numerically and its accuracy is assessed by performing one of the most challenging industrial turbulent two-phase flow problem; namely, the impinging of gaseous jet on a liquid interface. Finally, in Section 5, conclusions of the present work are drawn. 2. Mathematical formulations Over the last years a large number of computational methods have been developed for solving turbulent two-phase flows with moving interfaces. However, no pretence has been made that any of the previous numerical methods can be applied to all complex turbulent two-phase flows: such as ‘universal’ method may not exist. In the IMLS method, described here, the two immiscible fluids (such as liquid/gas) are treated separately; therefore, two sets of governing equations are introduced. Both fluids are assumed to be governed by the incompressible Reynolds form of the continuity and Navier–Stokes equations with constant properties. In order to extend the governing equations to the existing model of two immiscible fluids, two fluid domains X1 ðtÞ and X2 ðtÞ are considered. Both domains are evolving with time t and satisfying:

X ¼ X1 ðtÞ [ X2 ðtÞ;

X1 ðtÞ \ X2 ðtÞ ¼ ;:

ð1Þ

Each fluid has its own material properties qa and la (a = 1, 2), where the subscript a = l or g indicates the liquid and the gas phase presented at a given point in space. The discontinuity of the mass density and viscosity arise along the interface C, which denoted by:

CðtÞ : @ X1 ðtÞ \ @ X2 ðtÞ:

ð2Þ

2.1. Governing equations Taking the interaction of the two fluids at the interface into account, the instantaneous local equations of motion (known as RANS equations), in conventionally averaged variables, within each fluid after neglecting body forces can be written in tensor notation for two dimensional, unsteady, incompressible and axisymmetric flows, as follows:

3597

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

@ ðr qU i Þ ¼ 0; r@xi   @ @ @P @ @U i in Xa ða ¼ 1; 2Þ; ðqU i Þ þ U j ðr qU i Þ ¼  þ r l þ qRij @t r@xj @xi r@xj @xj

ð3Þ

where U is velocity vector, q is density, l is viscosity, P is pressure and Rij is the apparent turbulent stresses or Reynolds stresses (Rij ¼ ui uj ) where an overbar is used to denote Reynolds averaging. In the above set of equations, we consider that the upper case denotes ensemble-mean quantities and lower case indicates fluctuating or turbulence quantities. Different turbulence models have been previously developed in order to model the Reynolds stresses and to close the set of the momentum equations. In the following, the definition of linear and nonlinear k–e turbulence models is introduced. 2.2. Linear and nonlinear turbulence models The turbulent characteristics of two-phase flows can be obtained by solving the transport equation for kinetic energy k and its dissipation rate e. In such context, the linear standard two-equation k–e turbulence model (STD k–e) [16] applied in X can be written as follows:

  @ @ @ @k þ qðPk  eÞ; ðqkÞ þ U j ðr qkÞ ¼ r ðl þ lt =ck Þ @t r@xj r@xj @xj

ð4Þ

  @ @ @ @e e þ q ðC e1 Pk  C e2 eÞ; ðr qeÞ ¼ r ðl þ lt =ce Þ ðqeÞ þ U j @t r@xj r@xj @xj k

ð5Þ

where the turbulent viscosity lt = Clk2/e and ck ; ce ; C e1 ; C e2 ; C l are the model constants given as 1, 1.3, 1.44, 1.92, 0.09, respectively. The turbulent kinetic energy production rate Pk is given as follows:

Pk ¼ Rij

@U i : @xj

ð6Þ

The modeling of the production term in the turbulent flow simulation is one of the most important tasks as it affected by the approximation of the Reynolds stresses. Models following the Boussinesq assumption (1877) in approximating the Reynolds stresses are referred to linear turbulence models. In such models, a linear stress-strain relation can be defined as follows:

2 Rij ¼ mt Sij  kdij ; 3 where

ð7Þ

mt is the eddy viscosity, dij is the Kronker delta and Sij is the mean stress tensor defined as: Sij ¼

  @U i @U j : þ @xj @xi

ð8Þ

It is more conveniently to re-express the turbulent stress in terms of the dimensionless anisotropy tensor as follows:

aij ¼

ui uj 2 mt  dij ¼  Sij : 3 k k

ð9Þ

It is known that the linear eddy-viscosity models could not represent adequately the turbulence physics, particularly in respect of the different rates of production of the different Reynolds stresses and the resulting anisotropy, see for example [18]. The nonlinear models mimic the response of turbulence to complex strain by using a nonlinear constitutive relation between turbulent stress and mean rate of strain, i.e. roughly of form:

aij ¼ 

mt k

Sij þ Q 1 ðSij ; fij Þ þ Q 2 ðSij ; fij Þ;

ð10Þ

where the Q1 and Q2 represent quadratic and cubic functions of the mean strain Sij and the vorticity tensor fij that defined as:

fij ¼

  @U i @U j :  @xj @xi

ð11Þ

The quadratic eddy viscosity model, i.e. (Q2 = 0), has showed little range of applicability [37]. However, the cubic eddy viscosity model, (Q2 – 0), exhibits advantages over the quadratic model in allowing successful inclusion of normal-stress anisotropy and streamline curvature effects on turbulence. This model has been proposed in [20,37] and further implemented in [38], however, for single phase jet impinging problem. In such model, a nonlinear relation for the anisotropy tensor is used:

aij ¼ 

mt k

Sij þ









mt k 1 1 C S S  S S d þ C 2 ðSik fkj þ Sjk fki Þ þ C 3 fik fjk  flk flk dij 3 e 1 ik kj 3 kl kl ij

þ C 6 ðSkl Skl Sij Þ þ C 7 ðfkl fkl Sij Þ:

 þ

mt k2 ½C ðS f þ S f ÞS e2 4 ki lj kj li kl ð12Þ

3598

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

The model coefficients C1, C2, C3, C4, C6, C7 are given as 0.1, 0.1, 0.26, 10C 2l ; 5C 2l ; 5C 2l , respectively. The coefficient C l is now given by the expression:

 C l ¼ min 0:09;

1:2 1 þ 3:5g þ fRS

 ð13Þ

with

qffiffiffiffiffi S2I ;

ð14Þ

g ¼ ðk=eÞmax

nqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffio 0:5Sij Sij ; 0:5fij fij ;

ð15Þ

SI ¼ Sij Sjk Ski =ð0:5Snl Snl Þ1:5 :

ð16Þ

fRS ¼ 0:235ðmaxð0; g  3:333ÞÞ2

In relation to the flow of interest here, particular weakness of the linear eddy viscosity models related to turbulence production and insensitivity to streamline curvature can be avoided. In the present study, the performance of linear and nonlinear k–e turbulence models variants is investigated when these models are applied to predict a turbulent gas jet impinging on a free surface as one of the most challenging problem of two-phase flow for the first time according to our knowledge. 2.3. Near-wall modeling The above nonlinear model, as reported by [20], should contain additional low-Reynolds-number and near-wall terms in order to account for the wall damping effects. In such context, a very large number of nodes are required in order to resolve the flow right down to the wall. In addition, special viscosity-dependent modifications to the turbulence model are required. Instead of that, the wall effects on turbulence are modeled here via the wall-function and near-wall-modeling approaches [39] when applied either the linear or nonlinear turbulence models. Consequently, the low-Reynolds-number terms in the applied nonlinear model are neglected here; see for more details [38]. 2.4. Level set method In order to describe the complex topological changes of the interface separating the two phases with a robust and efficient method, the level set approach is applied. Since the introduction of the level set method in two-phase calculations [34], a large amount of bibliography on this subject has been published and several types of problems have been tackled with this method; for instance one can see the cited reviews of [35,40,41]. The basic idea of the LS method is to embed a moving interface C as the zero level set of a smooth phase function /, defined over the whole computational domain. Consequently, the two phases are distinguished by the level set function /(x, t), which can be defined as negative in the dispersed phase and positive in the bulk fluid. Fig. 1 shows a moving interface between two immiscible fluids with the definition of the level set function and the projection of 2D velocity components into a normal velocity Vn and tangential velocity Vt according to the normal vectors nx and ny in the direction normal and tangential to the interface, respectively. Mathematically, one can split the velocity components U in the direction of the normal and the tangent vector to the interface in order to get the normal and the tangential velocity components, i.e.



Vn Vt



 ¼

þni nj

   Ui  : ni Uj nj

Fig. 1. Definition of level set function and velocity components of a moving interface.

ð17Þ

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3599

In the numerical formulation of the level set, a smooth function / is typically initialized as a signed distance function from the interface, i.e. its value at any point is the distance from the nearest point on the interface and its sign is positive on one side and negative on the other. Let us set / as positive in liquid and negative in gas. The location of the interface is then given by the zero level set of the function /. In absence of the interfacial mass transfer such as evaporation or condensation, the equation of the level set method can be written in a form that is equivalent to the Hamilton–Jacobi equation:

@/ þ V n jr/j ¼ 0; @t

ð18Þ

where Vn is the normal velocity at the interface, defined as Vn = U  n. The unit normal on the interface n, drawn from the liquid into gas, and the curvature of the interface j can be defined in terms of / as:

n¼

r/ jr/j

j ¼ r  n:

and

ð19Þ

In general, the level set function is advected by using the normal velocity component defined for all the level sets in the computational domain. However, strictly speaking, the normal velocity has a physical meaning only at the interface location and it is independent of the others away from the interface. Therefore, the thinking about building of extension velocity field in the context of level set methods has received a great attention. One of the most recently developed numerical techniques for constructing of extension velocity field is proposed in [42], where the fast marching method (FMM) proposed by Sethian [43] is applied. In this method, the normal velocity Vn is replaced by some velocity field Fext known as the extension velocity, which at the zero level set, equals to the given speed Vn. In other words;

@/ þ F ext jr/j ¼ 0: @t

ð20Þ

The extension velocity Fext is calculated over the all computational domain X by solving the following equation:

rF ext  r/ ¼ 0;

where F ext ¼ V n jC :

ð21Þ

An important step in the solution algorithm of the level set function is to maintain the level set function as a distance function within the two fluids at all times, especially near the interface region, i.e. the Eikonal equation, jr/j ¼ 1, should be satisfied in the computational domain [34]. By the transient advection of the level set function, it ceases to be the signed distance function, although this property is tacitly employed by applying the FMM [42]. However, that is not accurately achieved in case of complex topological changes of the interface. Consequently, in the present work, the re-initialization algorithm described in [34] is applied for a specified small number of iterations in the all computational domain. 2.5. Interfacial boundary conditions Generally, the numerical simulations of two-phase flow could become quite a challenging problem and face significant difficulties. In particular, accurate tracking of complex interfaces, where large deformations and inter-penetration of phases occur, is required. In addition, the basis of improving such numerical simulations lies in satisfying the interface- normal and tangential stress boundary conditions more accurately. The accurate boundary conditions at the interface C separating two immiscible fluids are the balance of the dynamic stresses and the kinematic conditions. Assuming the stress tensor for turbulent flow of an incompressible fluid is given by:

Iij ¼ pI þ lD þ qRij ;

ð22Þ

where I is the identity matrix, D is the deformation tensor. This equation is applied at each point in the flow field; however, this equation loses its generality at the interface between the two phases as a result of the discontinuous fluid properties. The surface tension effect is known to balance the jump of the normal stress along the fluid interface [36]. In case of two immiscible fluids, the interfacial boundary conditions or jump conditions can be written as:

½Iij   n ¼ ðreff jn þ rt reff Þ;

ð23Þ

where the unit normal vector n is taken from liquid phase to gas phase and t is an arbitrary vector perpendicular to the normal to the interface. The bracket denotes the jump of the stresses along the fluid interface C, reff is the effective surface tension coefficient; equals to the molecular surface tension coefficient r and the turbulent surface tension coefficient rt and j is the average interface curvature similar to the laminar calculation. The modeling of the turbulent surface tension coefficient is discussed in more detail in [44], where it is concluded that the turbulent pressure fluctuation near the interface can result in a large surface tension force at the interface due to the large fluctuation of the local interface curvature. Such interfacial force can be modeled, in the context of RANS, either by considering a turbulence-induced surface tension effect with average interface curvature or a turbulence-induced curvature with the molecular surface tension coefficient. In general, it is suggested ffi pffiffiffiffiffiffiffiffiffi that the surface tension force in turbulent flow is increased by a factor of mt =m. In the present work, and analogy to [44], we assume that turbulent surface tension coefficient can be modeled as: 0:5 reff ¼ rð1 þ C r C 0:5 ke0:5 Þ; l m

where C r is the model constant to be determined due to calibration.

ð24Þ

3600

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

An alternative model for the jump condition described in Eq. (23), analogy to the above discussion, can be introduced by considering the turbulence-induced curvature with the molecular surface tension coefficient unchanged, i.e.

½Iij   n ¼ ðrjeff n þ rt rÞ;

ð25Þ

where jeff is the effective interface curvature; equals to the average local interface curvature j and the turbulent curvature term jt . The modeling of the turbulent curvature can be described as follows:





0:5 jeff ¼ j 1 þ C j C 0:5 ke0:5 ; l m

ð26Þ

where C j is the model constant to be determined due to calibration. By comparing Eqs. (24) and (26), it is observed that both equations are basically similar except for the different model coefficients. In the present study, we have used the turbulent curvature formulation of Eq. (26), with C j ¼ 1 [44]. The second term on the right hand side of Eq. (23) is the stress due to gradients on surface tension or Marangoni effect [45], usually important when a temperature gradient is applied parallel to the interface, e.g. thermo-capillary convection. Taking the projections of the jump conditions in the directions normal and tangential to the interface, considering no thermo-capillary convection effect, one obtains the following two equations in the normal and tangential directions, respectively:

½p þ leff ðrU þ rT UÞn ¼ rjeff n;

ð27Þ

½leff ðrU  nÞ  t þ leff ðrU  tÞ  n ¼ 0;

ð28Þ

where leff is the effective viscosity (leff = l + lt). It is clear from the above equation that the normal stresses are balanced by the surface tension force; however, the tangential stresses are continuous through the interface, or in other words, the tangential direction, an equality of the shear stress on both sides of the interface should be satisfied, i.e.

sl ¼ sg ;

ð29Þ

where s represents the total shear stress in the liquid phase and the gas phase, respectively. In addition to the equality of the dynamically interfacial stresses described above, the kinematic conditions should also be considered. When there is no mass transfer through the interface, the kinematic conditions is satisfied at the interface by assuming the continuity of the normal velocity component, i.e.

V n jl ¼ V n jg :

ð30Þ

2.6. Modeling of interfacial forces Based on the Continuum Surface Force Model [36], the standard level set formulation for incompressible and viscous twophase flows is based on expressing the surface tension effect in term of a singular source function defined in the momentum equations [46,47]. Consequently, the interfacial jump condition at the interface is integrated into the Navier–Stokes equations, resulting in a body force concentrated on the interface. A striking feature of the CSF model is the so called ‘‘spurious’’ currents. These numerical artifacts result from inconsistent modeling of the surface tension force and the associated pressure jump. In some cases, the effect of the unbalanced forces acting on the interface can reduce the accuracy of the calculation. Therefore, new surface tension models are still developed for the numerical simulation of two-phase flows [48]. Moreover, the CSF model considered that the effect of surface tension is to balance the jump of the normal stress along the fluid interface between two inviscid fluids having a constant surface tension coefficient. Consequently, the interfacial jump condition is reduced to Laplace’s formula for the surface pressure [49] due to the dropping of the viscous stresses. The effect of neglecting the interfacial shear stress can be seen through increasing of the normal stress magnitude at the phase boundary [50,51]. Moreover, the formation of the capillary ripples on a thin liquid film is dependent on receiving their energy essentially from the pressure and shear stress exerted at the interface by the flowing concurrently high velocity gas stream [52]. In contrast to CSF model, the jump conditions at the interface between the two immiscible fluids are treated here as boundary conditions for pressure enforced explicitly at the interface. This numerical technique has been implemented in our previous work [53,54], and good results have been obtained for a number of challenging problems. The idea of our modeling is straightforward. By introducing a number of so called ‘‘interfacial markers’’ on the intersection points of computational grids with the interface, the local average curvature of the interface is calculated using an interpolation technique from the regular grid point to the ‘‘interfacial markers’’ in order to represent the interface values. Referring to Fig. 2, the average interface curvature can be expressed by / and its derivatives at the regular grid points according to the formula

jð/Þ ¼

ð2/x /y /xy  /2x /yy  /2y /xx Þ ð/2x þ /2y Þ1:5

;

ð31Þ

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3601

Fig. 2. Computational domain and grid points.

where the first and the second derivatives of the level set function is calculated according to the central difference approximation. The local interface curvature at the interfacial marker point b, for example, can be calculated according to:

jb ð/Þ ¼ ð1  f ÞjP ð/Þ þ f jB ð/Þ; f ¼

/P : ð/P  /B Þ

ð32Þ

Consequently, the above equation can be applied for the calculation of the interphase boundary value for any arbitrary variable. The calculation of the local curvature at the interfacial markers enables us to calculate the surface tension force, i.e. the surface pressure, and then it can be used to drive the liquid phase through the pressure gradient seen in the momentum equations. Moreover, at the grid points near the interface, the gradient of variables can be calculated by means of the appropriate values of the properties at the interfacial marker and by using the actual distance to the interface (ha, hb, hc, hd). The present surface tension model ensures that both the pressure calculated within the liquid phase and the surface tension pressure is consistent and dynamically similar, as their effect is determined in the same way. Accordingly, the pressure drop across the interface cancels exactly the surface tension potential at the interface. Such modeling could, to some extent, reduce the spurious current problem that might be associated with the inconsistent modeling of the surface tension force in twophase flow simulation, as the interfacial curvature and, consequently, the interfacial driving forces are calculated accurately. For more generality of the present model, it is considered that the interfacial pressure inside the liquid phase pl is determined by evaluating the interfacial pressure in the gas phase pg and other different ‘‘pseudo pressure’’ terms, i.e.

pl ¼ pg þ pr þ pl þ pgr ;

ð33Þ

where pr ; pl and pgr represent the interfacial ‘‘pseudo pressure’’ due to the surface tension effect, the viscous normal stresses and the buoyancy effect. This formulation has been successfully applied in [53–55]. In such context, the present model enables us to simply include any external interfacial driving force in the pressure boundary condition, e.g. magnetic or electric field, instead of incorporating it as a body force in the momentum equations. The pressure values calculated from the above equation is then used as Dirichlet boundary conditions for solving the Poisson equation for the pressure. Consequently, this model can, in principle, be applied to any immiscible fluids problem as long as the stress on the second phase can be specified. The satisfying of the previous interfacial boundary conditions is an important task in the numerical simulation of twophase flows as the pressure and velocity field inside the liquid phase are caused by the normal and tangential stresses of the external gas field. Therefore, the exact pressure level inside the liquid phase, which considers as the driving force, should be accurately defined. 3. IMLS numerical method The numerical methods for two-phase flows can be distinguished by their ways of determining the pressure field or their techniques to capture the interface separating the two immiscible fluids. In general, the calculation of the accurate pressure field has the most important effect on the topological changes of the interface rather than other fluid variables; see for more details [56]. Moreover, in the numerical simulation of turbulent two-phase flow, the inertial effects are predominant [57]. According to the difficulties of measuring flow field variables such as pressure and velocity fields near free surface flow, numerical simulation is a powerful way to gain an understanding of the interfacial dynamics. So far, much previous researches have been performed using different numerical techniques [58–60].

3602

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

Fig. 3. Illustration of the irregular mesh caused by an interface separating two phases.

The standard level set method for incompressible two-phase flows, originated nearly twenty years ago [34], has become very popular. In this method, the instantaneous local momentum equations within each fluid are solved with appropriate numerical scheme. However, the direct solution of the momentum equations in such case is naturally difficult, since the fluid properties change discontinuously across the interface, and the concentrated surface tension also becomes infinite in an infinitesimal volume. To overcome these difficulties, the interface is smoothed across a finite thickness region, usually a few grid points thick. Moreover, the fluid properties are smoothed according to the Heaviside function that depends on a parameter describing the thickness of the interface or in other words a ‘‘transition region’’, see for more details [47]. Although the problem of the transition region is solved by using a sharp interface immersed boundary formulation [61] and a level-set/ghost-fluid method [62], the resulting approach introduces an additional computational cost that can be very demanding in 3D numerical simulation. In contrast to the standard level set formulation, in the present paper, the two phases are solved separately by using the control volume approach on structured cell-centered collocated grids, shown in Fig. 3. A brief, descriptions of the developed numerical method is described in the following. The governing equations are discretized and solved on the basis of the control volume approach proposed by Patankar [63]. As a result of the complexities introduced by the staggered grid system in two-phase flow calculations and the need to track the liquid surface, a non-staggered grid system is applied here, which requires the use of a single cell network for all variables and a collocated specification of variables at the center of each cell. As there is no pressure transport equation necessitates the consideration of the continuity equation, Poisson equation for pressure is solved by means of the Successive Over-Relaxation method. 3.1. Computational procedure In our algorithm, the implicit fractional step-non iterative method is applied to obtain the velocity and pressure field in each phase according to its properties by presuming that the velocity field reaches its final value in two stages; that means

Unþ1 ¼ U þ Uc ;

ð34Þ



where U is an imperfect velocity field based on a guessed pressure field, and Uc is the corresponding velocity correction. Firstly, the ‘starred’ velocity will result from the solution of the momentum equations. The second stage is the solution of Poisson equation for the pressure:

r2 pc ¼

qa Dt

r  U ;

ð35Þ

where pc is called the pressure correction and Dt is appropriate time step. The prescribed time step is considered to be constant and small enough to ensure a stable and oscillation free solution. Once the Poisson equation is solved, one gets the appropriate pressure correction, and consequently, the velocity correction is obtained according to:

Uc ¼ 

Dt

qa

rpc :

ð36Þ

This fractional step method described above ensures the proper velocity-pressure coupling for incompressible flow field. It should be pointed out that, more recently; this method has been successfully applied for different single phase flow configurations [18].

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3603

3.2. Solution of Poisson equation in irregular domain In two-phase flow problems, the interface shape may be of complicated geometry. In such context, the accurate solution of the surface pressure occurring at transient fluid interfaces of arbitrary and time dependent topology is a key element of the proposed method, which involves a special treatment of the discretization at the interface [64]. Assuming that a regular mesh is used for the calculation, the curved shape of the interface causes unequal spacing between the interface and some internal grid points, as shown in Fig. 2. The Dirichlet pressure boundary conditions are assigned at the interfacial marker positions by applying Eq. (33) for the calculating the interfacial liquid pressure according to the case considered. The source term of Eq. (35) is calculated at the main grid points using the imperfect velocity field calculated from the momentum equations and approximated at the irregular grid. In order to solve Poisson equation, it should be first generalized for an arbitrary computational grid as shown in Fig. 2, in which the four neighboring points are at different distances from the central point of the control volume according to the interface position. Referring to Fig. 2, the Poisson equation can be approximated at point p by expanding the pressure value at the neighboring points using Taylor’s series as follows, e.g.

!   2 @p ha @ 2 p 3 pa ¼ pij  ha þ  Oðha Þ; @x ij 2 @x2

ð37Þ

!   2 @p hb @ 2 p 3 pb ¼ pij þ hb þ þ Oðhb Þ: @x ij 2 @x2

ð38Þ

ij

ij

By dropping the cubic and higher-order terms and by applying some mathematical manipulations and rearrangements, the desired approximation of Poisson equation can be expressed as:

pij ¼



 pa pb pc pc þ þ þ þ Sp ; ha ðha þ hb Þ hb ðha þ hb Þ hc ðhc þ hd Þ hd ðhc þ hd Þ

ð39Þ

where Sp is the pressure source term calculated according to the right-hand side of Eq. (35). It can easily be shown that the above formula is equivalent to the regular grid formula if the distances ha = hb = Dx, and hc = hd = Dy. Our developed strategy applied for obtaining the dynamics of the interface separating the two immiscible fluids is as follows; the complete set of the governing equations are solved for the gaseous field to obtain the velocity, pressure, viscous stresses and any existing force around the interface. By using the velocity and pressure values on the gas phase as boundary conditions defined on the interface, the solution of the momentum equations and Poisson equation in the liquid phase can be carried out. After the whole computational domain is calculated, the turbulent equations are solved on both phases simultaneously. The normal velocity at the interface is then used to advect the interface using the level set approach and to obtain its topological changes. Consequently, the whole algorithm is repeated until it would reach the specified time or the statistically steady state condition. It should be pointed out that, the basic idea of the present approach has been applied previously by the present author for calculating the dynamics of the interface separating two phases in simple as well as complex two-phase flow applications [53,54,65]. However, the present approach with the implementation of nonlinear turbulence models has not been applied previously for predicting the interaction dynamics of turbulent flow with a moving free surface. 4. Numerical simulation In this section, the validation of our developed IMLS numerical method is carried out through performing one of the most challenging problems of two-phase with moving interface; namely, the deformation of a liquid interface upon the impinging of turbulent gaseous jet. In such case, the interfacial boundary conditions are affected by the surface tension force, the viscous stresses, the gravity force and the turbulence quantities. A new dimensionless number, called Impinging number (IM), which describe the impinging process and the related operating parameters is introduced. The presented numerical results for the considered case are compared with the available analytical results and experimental measurements. 4.1. Liquid surface deformation due to impinging process The impingement of an axisymmetric gaseous jet through lances upon a liquid surface or a molten bath is relevant in industrial and metallurgical processes, e.g. dust removal from the producer gas obtained from biomass gasification and electrical arc furnace steel production. As a result of the kinetic energy of the jet stream, the surface of the liquid will be depressed and a cavity is created on the surface of the liquid. The experimental visualization of such cavity formation is a difficult task, so that the structure of the jet impingement on a liquid surface is not completely understood [66]. Moreover, the analytical and theoretical treatments are based on great simplifications and assumptions [67]. Consequently, a clear understanding of such problem is still required through numerical simulation [68]. The important characteristics in such two-phase impinging process are the interface shape, cavity depth and width, lip formation and the recirculation pattern in the liquid. An extended review of the most important articles concerned such

3604

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

Fig. 4. (a) The features of the jetting system. (b) The force analysis exerted upon the liquid surface.

impinging process up to 2010 can be found in [68]. The problem configuration is shown in Fig. 4a, while Fig. 4b shows the analysis of forces acting upon the liquid surface during the impinging process and as the cavity shape is developed. The parameters described in the Fig. 4a for the jetting system are given as follows; the nozzle diameter Dn, vessel diameter Dc, the height of the undisturbed liquid H, distance between nozzle tip and liquid surface h, cavity depth hc, cavity diameter dc, lip height hL, and the momentum of the jet is given by M ¼ 0:25pD2n qg V 2j , where qg and Vj is the gas jet density and velocity, respectively. The Reynolds number, consequently, is given by Re = qgVjDn/lg. On the other hand, the schematic of the force analysis can be seen in Fig. 4b, where Fm, Fs, Fp, Fr, and Fg are the associated forces due to jet momentum, tangential shear stress, pressure force due to recirculatory flow, surface tension force and gravity force. In the present algorithm and in contrast to the CSF model [36], all the interfacial forces are modeled as pressure boundary conditions at the interfacial marker points on the interface. As a measure of the algorithm accuracy, a statistically steady state cavity shape that formed due to jet impingement should be obtained after a specified period of time corresponding to the operating parameters. This can only be achieved if the interfacial forces are dynamically balanced, i.e.

X

~ F i ¼ 0:

ð40Þ

The problem of how to numerically formulate these interfacial boundary conditions has, however, not been completely solved. When the associated modeling of such interfacial forces is not accurate, consequently, it is expected that no statically steady state cavity will be obtained. Therefore, the obtained steady state cavity shape over long time intervals is considered as a challenge problem of such simulations and a challenge validation for our developed numerical method. 4.2. Developing of a new dimensionless number (IM-number) The characteristics of the cavity formed in the impinging process are dependent on the properties of the impinging jet as well as the liquid properties. When the jet momentum is small, a steady state cavity is just formed. On the other hand, a lip can appear at larger jet momentum. Further large jet momentum causes splash and ripples on the liquid surface. The present computation does not consider splash and ripples. In order to combine the effect of the operating parameters as well as the fluids properties in such jetting system, the dimensional analysis theory is applied to obtain a new dimensionless number; we call here the Impinging number (IM), which can describe the impinging process in terms of the gaseous jet and the liquid bulk properties. The IM number is expressed as:

IM ¼

qg V j Dn  lg

pffiffiffiffiffiffiffiffiffiffiffiffi ql rH

ll

:

ð41Þ

The above equation can be written in form of the known dimensionless number, Reynolds number (Re) and Ohnesorge number (Oh), as follows:

IM ¼

Re : Oh

ð42Þ

The Reynolds number is considered as a measure of the jet momentum while the dimensionless Ohnesorge number measures the ratio of the viscous force on an element of fluid to the surface tension and the inertia force. The combination of both numbers could describe the strength of the impinging process adequately.

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3605

In the present numerical experiment, the Ohnesorge number is taken constant and equals 0.94, while the jet velocity is changed to obtain a wide range of Re number and consequently a wide range of IM number. 4.3. Computational grid and boundary conditions The numerical simulation for the considered problem is carried out on 91  93 grid with Dx = 1.14e3, Dy = 2.1e3 m and a constant time step of Dt = 5e4 s. Different types of boundary conditions were used to describe the flow field within the computational domain used to solve the jetting system shown in Fig. 4a. The gas jet inlet boundary conditions imposed at the nozzle exit were described by uniform profiles for the inlet jet velocity, kinetic energy of turbulence and its dissipation rate are described; as follows:

U ¼ U inlet ;

V ¼ 0;

k ¼ kinlet ¼ 1:5ðIU inlet Þ2 ; I ¼ u0 =U inlet ¼ 0:16ðReÞ0:125 ;

e ¼ einlet ¼ ðC l Þ

0:75

ð43Þ

1:5

ðkinlet Þ =l;

l ¼ Dn ; where the turbulence length scale l, at the inlet is assumed to be the nozzle diameter. The turbulence intensity I is considered to be a function of the Reynolds number based on the nozzle diameter and inlet velocity. At the wall boundaries, the no-slip boundary conditions were imposed. Moreover, the wall effects on turbulence were evaluated via the wall-function and nearwall-modeling approaches [16]. Since the flow is axisymmetric, the axis boundary was used along the nozzle centerline, at which the radial velocity component and the gradients of the other dependent variables were equal to zero. For other boundaries, the slip boundary conditions were applied. The gravity is assumed to be lumped with the pressure jump across the interface. A piezometric pressure has to be defined at the interface; pg = qgz, where z is the vertical coordinate at the interface. For a more detailed discussion one can see Ref. [55]. 4.4. Results and discussion Fig. 5 shows the effect of the IM number on the obtained statistically steady cavity profiles for different IM numbers, ranging from small to medium and large numbers for jet spacing of (h = 5 cm). The axial velocity contours of the impinging jet is superimposed over the bulk liquid and plotted in the same figure. The formation of dimples is clearly evident only at relatively high IM number. At low IM number, the dimples formed die out due to dissipation of their energy by viscous friction. The circulation of the liquid and the velocity vectors in both phases are shown in Fig. 6a, while the level set contours are illustrated in Fig. 6b at IM = 9771 and a jet distance of h = 5 cm. The tangential flow and the development of the recirculation regions are clearly visible with a center of the circulating flow in the liquid phase close to the centerline. It can also be seen that, in the opposite direction to the impinging jet, upward flow is built up due to the circulation initiated in the liquid phase. The existence of such upward flows, to some extent, can be considered as a stabilizing force that could prevent further increase of the cavity depth. The level set contours obtained in Fig. 6b reveals the capability of the present developed level set technique in capturing the smooth topological changes of the interface, especially the lip formation at the ends of the formed cavity. An important feature of the present case, which also considered as a challenging problem, is the obtaining of the steady cavity profile over a large period of time. The statistically steady cavity profile obtained reveals the accurate modeling implemented for the interfacial forces presented at the interface in the case considered. Fig. 7a shows the evolution of the cavity depth for a wide range of IM numbers during the jetting process until it reaches to an equilibrium position with observable small oscillation of the cavity depth about this equilibrium position. It can be seen that, for all IM numbers, the liquid interface is stable within the first few seconds until the gas jet starts to impinge the liquid surface and then the cavity is formed. The pressure force of the impinging jet which is provided by the momentum flux of the gaseous jet causes the liquid surface

Fig. 5. The axial jet velocity contours superimposed over the statistically steady state cavity profile at different IM number.

3606

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

Fig. 6. (a) Velocity vectors in the jetting system at IM = 9771 and jet spacing h = 5 cm. (b) Level set contours.

to deform whereby a cavity is formed with a specified depth. Consequently, lateral flow is developed and moved over the liquid surface generating a shearing force along the surface that drives the liquid towards the walls of the cylinder and producing symmetrically recirculation cells on both sides of the cavity, as shown previously in Fig. 6a. In Fig. 7b, the time evolutions of the topological changes of the liquid interface, starting from a flat interface to the final cavity shape, are illustrated. The equilibrium position is obtained when all the driving forces are dynamically balanced. In almost all cases, the cavity depth reaches a maximum value after a specified time as the liquid moves downward and then it rises again in a direction opposite to the gravity in order to satisfy the continuity by refilling the cavity created by the air jet impingement, then the cavity depth remains almost constant with an oscillation around an average value. The most important parameters of the formed cavity are the interface shape, the width and the depth of the cavity and the height of the peripheral lip. Therefore, the numerical prediction of such parameters should be evaluated by comparing the obtained results with those predicted either analytically [67] or measured experimentally [69]. Consequently, a series of

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3607

Fig. 7. The time evolution of the cavity depth and shape during the jetting process: (a) for a wide range of IM numbers; (b) at different time intervals.

numerical experiments has been performed over a wide range of the IM number by using both the STD k–e turbulence model and the nonlinear turbulence model.

3608

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

Fig. 8. Comparison of the numerical results with the experimental and theoretical data for cavity depth (a) and cavity radius (b).

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3609

The comparison of the cavity depth and the cavity radius with the available analytical and experimental data is illustrated in Fig. 8a and b. The jet spacing is considered to be constant (h = 5 cm) for all numerical experiments. The comparison of the numerically predicted cavity depth, shown in Fig. 8a, reveals that the numerical results obtained by using both turbulence models indicate fairly general agreement with the experimental measurements rather than the analytical prediction. This can be attributed to the fact that the theory did not take into account the effects of surface tension forces and the existence of recirculation in the liquid phase. By comparing both turbulence models adopted, it can be seen that the numerical results obtained from the nonlinear turbulence model are in closer agreement with the experimental date than the results obtained from the STD k–e model. The overestimated cavity depth, obtained from the STD k–e model, can be referred to the excessive generation of the turbulent production at the impinging point with the liquid surface. This is a general problem of linear turbulence models which produces incorrect prediction of anisotropic turbulence in many industrial important flows such as impinging jets [18]. Consequently, the accurate flow field prediction could not be well predicted at the stagnation point where gas jet impinging the liquid surface. The comparison of the cavity radius, shown in Fig. 8b, illustrates that the present numerical results from both linear and nonlinear turbulence models overpredict both of the experimental measurements and the analytical results. However, the experimental trend was well numerically predicted. It should be pointed out that a part of the disagreement may be attributed to both of the associated error in the measurement of the cavity width due to the error in measurements of the cavity radius [69]. On the other hand, the lip formation phenomena could not be predicted analytically, and therefore, there is uncertainty in the theoretical definition of the cavity width. The present numerical prediction of the cavity radius indicates that no significant difference has been obtained by using either STD k–e or nonlinear turbulence models. This implies that the tangential stresses predicted from both turbulence models are nearly similar. In general, the increase in the cavity depth and radius with increasing the IM number is clearly evident. 4.5. Prediction of centerline velocity The adopted turbulence models are now used for predicting the centerline velocity, VCL, starting from the nozzle exit. The numerical results are compared with the analytical equation developed in [70] and based on the theoretical analysis provided by [71]:

V CL B ¼ ; z=Dn þ B VJ

ð44Þ

Fig. 9. Comparison of the predicted centerline velocity with the theoretical analysis.

3610

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

where B is a function of the gas jet Reynolds number and is given a value of 7.6 in [70]; for the specified case where Dn = 6 mm, DC = 290 mm, h = 154 mm, H = 111 mm and VJ = 56.2 m/s. It should be pointed out that, the above equation, Eq. (44), has been also used for validating the numerical results obtained from the commercial CFD code Fluent [72]. Fig. 9 shows the comparison between the centerline jet velocity, predicted by using the linear and the nonlinear turbulence models adopted, and the theoretical results of Eq. (44). It can be seen that the predicted numerical results from both turbulence models agree well with the analytical solution nearly up to the point at which the jet enters the cavity region (at z/Dn  50). The agreement between the linear and the nonlinear turbulence models adopted in predicting the mean centerline velocity reveals that the nonlinear constitutive relation has no effect over the prediction of the mean variables. However, the nonlinear formulation has a significant effect on the cavity depth which influenced by the correct prediction of the anisotropic normal stresses at the impinging region (cf. Fig. 8a). Consequently, there is a need for refined experimental measurements of the cavity radius in such applications in order to validate the numerical prediction from both turbulence models adopted. 5. Conclusions A new interfacial marker-level set method for predicting the dynamics of two-phase flows with moving interface has been developed. The governing unsteady RANS-equations are coupled with the level set method and solved in each phase separately on structured cell-centered collocated grids using the control volume approach on the physical domain of the problem considered. By fitting a number of interfacial markers on the intersection points of computational grids with the interface, the interfacial stresses and consequently the interfacial driving forces are easily estimated. Moreover, the normal interface velocity, calculated at the interfacial markers position, are extended to the higher dimensional level set function and applied for the advection algorithm of the level set method via the fast marching method. The developed numerical method has been applied to simulate one of the most challenging problems in two-phase flow applications; namely, the impinging of turbulent gaseous jet onto a liquid surface. Both linear and nonlinear turbulence models have been applied to predict the turbulence effects of the gaseous jet and the interaction with the free surface. The numerical results obtained show advantages of the nonlinear turbulence model over the linear STD k–e model in predicting the parameters associated with the anisotropic normal stresses. However, a similar prediction of the mean variables has been obtained from both turbulence models adopted. In general, the developed IMLS method demonstrates a remarkable capability to predict the dynamical characteristics of the considered two-phase immiscible flow problems and therefore it can be applied to quite a number of two-phase applications related to different industrial and engineering applications. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

A.H. Lefebvre, Atomization and Sprays, Hemisphere Publishing Corporation, 1989. H. Won, N. Peters, Investigation of cluster-nozzle concepts for direct injection diesel engines, Atomization Sprays 19 (10) (2009) 983–996. M. Garbero, M. Vanni, U. Fritsching, Gas/surface heat transfer in spray deposition processes, Int. J. Heat Fluid Flow 27 (2006) 105–122. P.A. Torpey, Prevention of air ingestion in a thermal ink-jet device, in: Proceedings of the 4th International Congress on Advances in Non-impact Print Technologies, Springfield, VA, 1988. D. Dressler, L. Li, S. Green, M. Davy, D. Eadie, Newtonian and non-newtonian spray interaction with a high-speed moving surface, Atomization Sprays 19 (2009) 19–39. Li. Zhaorui, Farhad A. Jaberi, Tom I.-P. Shih, A hybrid Lagrangian–Eulerian particle-level set method for numerical simulations of two-fluid turbulent flows, Int. J. Numer. Methods Fluids 56 (2008) 2271–2300. O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost fluid method for simulation turbulent atomization, J. Comput. Phys. 227 (2008) 8395–8416. M. Herrmann, A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids, J. Comput. Phys. 227 (2008) 2674– 2706. J. Eggers, Nonlinear dynamics and breakup of free-surface flows, Rev. Mod. Phys. 69 (3) (1997) 1–65. R.M. Schulkes, The evolution and bifurcation of a pendant drop, J. Fluid Mech. 278 (1994) 83–100. G. Tryggvason, A. Esmaeeli, J. Lu, S. Biswas, Direct numerical simulations of gas/liquid multiphase flows, Fluid Dyn. Res. 38 (2006) 660–681. R. Lebas, T. Menard, P.A. Beau, A. Berlemont, F.X. Demoulin, Numerical simulation of primary break-up and atomization: DNS and modeling study, Int. J. Multiphase flow 35 (3) (2009) 247–260. G. Luret, T. Menard, A. Berlemont, F.X. Demoulin, A DNS study ranging from dense to dilute turbulent two-phase flows, in Seventh International Conference on Multiphase Flow, ICMF, FL, USA, 2010. P. Liovic, D. Lakehal, Interface–turbulence interaction and bubble dynamics, in: Seventh International Conference on CFD in the Minerals and Process Industries, CSIRO, Melbourne, Australia, 2009, pp. 9–11. A. Balabel, A.M. Hegab, M. Nasr, S.M. El-Behery, Assessment of turbulence modeling for gas flow in two-dimensional convergent–divergent rocket nozzle, Appl. Math. Model. 35 (2011) 3408–3422. B.E. Launder, D.B. Spalding, Lectures in Mathematical Models of Turbulence, Academic Press, London and New York, 1972. p. 169. W. El-Askary, A. Balabel, Prediction of reattachment turbulent shear flow in asymmetric divergent channel using linear and non-linear turbulence models, Eng. Res. J. (ERJ), Faculty Eng. Menoufiya Univ. 30 (4) (2007) 535–550. A. Balabel, W. El-Askary, On the performance of linear and nonlinear k–e turbulence models in various jet flow applications, Eur. J. Mech. B/Fluids 30 (2011) 325–340. S.B. Pope, A more general effective-viscosity hypothesis, J. Fluid Mech. 72 (1975) 331–340. T.J. Craft, B.E. Launder, K. Suga, Development and application of a cubic eddy-viscosity model of turbulence, Int. J. Heat Fluid Flow 17 (1996) 108–115. D.D. Apsley, M.A. Leschziner, A new low-Reynolds-number nonlinear two-equation turbulence model for complex flows, Int. J. Heat Fluid Flow 19 (1997) 209–222.

A. Balabel / Applied Mathematical Modelling 36 (2012) 3593–3611

3611

[22] S. Popinet, S. Zaleski, A front tracking algorithm for the accurate representation of surface tension, Int. J. Numer. Methods Fluids 30 (1999) 775–793. [23] S. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [24] B.D. Nichols, C.W. Hirt, Methods for calculating multi-dimensional, transient free surface flows past bodies, in: Proceedings of the First International Conference on Numerical Ship Hydrodynamics, Gaithersburg, 1975, pp. 20–23. [25] S. Osher, J. Sethian, A Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [26] P. Trontin, S. Vincentm, J.L. Estivalezes, J.P. Caltagirone, Detailed comparisons of front-capturing methods for turbulent two-phase flow simulations, Int. J. Numer. Methods Fluids. 56 (2008) 1543–1549. [27] N. Peters, Turbulent Combustion, Cambridge University Press, Cambridge, UK, 2000. [28] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling, merging and fragmentation in multiphase flows with SURFER, J. Comput. Phys. 113 (1994) 134–147. [29] Th. Schneider, R. Klein, Overcoming mass losses in level set-based interface tracking schemes, in: Proceedings of the Finite Volume for Complex Application II, 1999, pp. 41–50. [30] Y.T. Albert, W. Zhaoyuan, A numerical method for capillary-dominant free surface flows, J. Comput. Phys. 221 (2007) 506–523. [31] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys. 183 (2002) 83–116. [32] M. Herrmann, A Eulerian level set/vortex sheet method for two-phase interface dynamics, J. Comput. Phys. 203 (2005) 539–571. [33] M. Garzon, L.J. Gray, J.A. Sethian, Numerical simulation of non-viscous liquid pinch-off using a coupled level set-boundary integral method, J. Comput. Phys. 228 (2009) 6079–6106. [34] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flows, J. Comput. Phys. 114 (1994) 146– 159. [35] J.A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annu. Rev. Fluid Mech. 35 (2003) 341–372. [36] J.U. Brackbill, D.B. Kothe, A. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [37] T.J. Craft, H. Iacovides, J.H. Yoon, Progress in the use of non-linear two equation models in the computation of convective heat-transfer in impinging and separated flows, Turbul. Combust. 63 (1999) 59–80. [38] T.J. Craft, H. Iacovides, N.A. Mostafa, Modelling of three-dimensional jet array impingement and heat transfer on a concave surface, Int. J. Heat Fluid Flow 29 (2008) 687–702. [39] B.E. Launder, D.B. Spalding, The numerical computation of turbulent flows, Comput. Methods Appl. Mech. Eng. 3 (2) (1974) 269–289. [40] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759. [41] S. Zahedi, K. Gustavsson, G. Kreiss, A conservative level set method for contact line dynamics, J. Comput. Phys. 228 (2009) 6361–6375. [42] D. Adalsteinsson, J.A. Sethian, The fast construction of extension velocities in level set methods, J. Comput. Phys. 148 (1999) 2–22. [43] J.A. Sethian, A marching level set method for monotonically advancing fronts, Proc. Natl. Acad. Sci. 93 (4) (1996) 1591–1595. [44] E. Shirani, A. Jafari, N. Ashgriz, Turbulence models for flow with free surfaces and interfaces, AIAA J. 44 (7) (2006) 1454–1462. [45] C. Marangoni, Ueber die Ausbreitung der Tropfen einer Flüssigkeit auf der Oberfläche einer Anderen, Ann. Phys. Chem. 143 (1871) 337–354. [46] Y.C. Chang, T.Y. Hou, B. Merriman, S. Osher, A level set formulation of Eulerian interface capturing method for incompressible fluid flows, J. Comput. Phys. 124 (1996) 449–464. [47] T. Watanabe, Numerical simulation of oscillations and rotations of a free liquid droplet using the level set method, Comput. Fluids 37 (2008) 91–98. [48] S. Zhongguo, Xi Guang, Chen Xi, Numerical simulation of binary collisions using a modified surface tension model with particle method, Nucl. Eng. Des. 239 (2009) 619–627. [49] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, Pergamon, London, 1959. [50] T.B. Benjamin, Shearing flow over a wavy boundary, J. Fluid Mech. 6 (1959) 161–205. [51] A.M. Sterling, C.A. Sleicher, The instability of capillary jets, J. Fluid Mech. 68 (3) (1975) 477–495. [52] J.C. Asali, T.J. Hanratty, Ripples generated on a liquid film at high gas velocities, J. Multiphase Flow 19 (2) (1993) 229–243. [53] A. Balabel, B. Binninger, M. Herrmann, N. Peters, Calculation of droplet deformation by surface tension effects using the level set method, J. Combust. Sci. Technol. 174 (11–12) (2002) 257–278. [54] A. Balabel, A. Hegab, Numerical modeling of the interfacial instability in two-phase flows considering buoyancy and capillary phenomena, in: Proceedings of the 12th ASAT Conference, 29–31 May, CFD-07, 2007, pp. 1–16. [55] P.M. Carrica, R.V. Wilson, F. Stern, An unsteady single-phase level set method for viscous free surface flows, Int. J. Numer. Methods Fluids 53 (2007) 229–256. [56] M. Seifollahi, E. Shirani, N. Ashgrizb, An improved method for calculation of interface pressure force in PLIC-VOF methods, Eur. J. Mech. B/Fluids 27 (2008) 1–23. [57] S. Vincent, J. Larocque, D. Lacanette, A. Toutant, P. Lubin, P. Sagaut, Numerical simulation of phase separation and a priori two-phase LES filtering, Comput. Fluids 37 (2008) 898–906. [58] K.A. Pericleous, M. Cross, G. Moran, P. Chow, K.S. Chan, Free surface Navier–Stokes flows with simultaneous heat transfer and solidification/melting, Adv. Comput. Math. 6 (1996) 295–308. [59] C. Bailey, D. Wheeler, M. Cross, An integrated modeling approach to solder joint formation, IEEE Trans. Comp. Packag. Technol. 2 (1999) 497–502. [60] E. Samiei, M. Shamsa, R. Ebrahimi, A novel numerical scheme for the investigation of surface tension effects on growth and collapse stages of cavitation bubbles, Eur. J. Mech. B/Fluids 30 (2011) 41–50. [61] J. Yang, F. Stern, Sharp interface immersed-boundary/level-set method for wave body interaction, J. Comput. Phys. 228 (17) (2009) 6590–6616. [62] R.P. Fedkiw, T. Aslam, B. Merriman, S.J. Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the Ghost fluid method), J. Comput. Phys. 154 (1999) 393–427. [63] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. [64] F. Gibou, R. Fedkiw, L.T. Cheng, M. Kang, A second order accurate symmetric discretization of the Poisson equation on irregular domains, J. Comput. Phys. 176 (2002) 205–227. [65] A. Balabel, Numerical prediction of droplet dynamics in turbulent flow, using the level set method, Int. J. Comput. Fluid Dyn. 25 (5) (2011) 239–253. [66] D. Nakazono, K. Abe, M. Nishida, K. Kurita, Supersonic O2-jet impingement on liquid iron with surface chemistry, ISIJ Int. 44 (1) (2004) 91–99. [67] F.R. Cheslak, J.A. Nicholls, M. Sichel, Cavity formed on liquid surfaces by impinging gaseous jets, J. Fluid Mech. 36 (1969) 55–64. [68] A. He, A. Belmonte, Deformation of a liquid surface due to an impinging gas jet: a conformal mapping approach, Phys. Fluids 22 (2010) 1–7. 042103. [69] S. Eletribi, D.K. Mukherjee, V. Prasad, Experiments on liquid surface deformation upon impingement by a gas jet, in: Proceedings of the ASME Fluids Engineering Division, vol. 244, 1997, pp. 235–247. [70] D.E. Wakelin, The Interaction Between Gas Jet and the Surface of Liquids Including Molten Metals, PhD, University of London, UK, 1996. [71] R.B. Banks, D.V. Chandrasekhara, Experimental investigation of the penetration of a high-velocity gas jet through a liquid surface, J. Fluid Mech. 15 (1963) 13–34. [72] A.V. Nguyen, G.M. Evans, Computational fluid dynamics modelling of gas jets impinging onto liquid pools, Appl. Math. Model. 30 (2006) 1472–1484.