Large Eddy Simulation of liquid sheet breakup using a two-phase lattice Boltzmann method

Computers and Fluids 160 (2018) 93–107

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Large Eddy Simulation of liquid sheet breakup using a two-phase lattice Boltzmann method Hadi Amirshaghaghi a, Mohammad Hassan Rahimian a,∗, Hesameddin Safari b, Manfred Krafczyk b a b

School of Mechanical Engineering, College of Engineering, University of Tehran, Iran Institute for Computational Modeling in Civil Engineering, Technische Universität Braunschweig, Braunschweig, Germany

a r t i c l e

i n f o

Article history: Received 27 January 2017 Revised 29 July 2017 Accepted 23 October 2017 Available online 25 October 2017 Keywords: Liquid sheet Breakup Two phase Lattice Boltzmann method Turbulence, high density ratios

a b s t r a c t A two-phase flow Large Eddy Simulation (LES) model is developed within the framework of phase-field lattice Boltzmann methods (LBM). The proposed model is based on a standard Smagorinsky model for subgrid closure and uses the Multiple-Relaxation-Time (MRT) collision operator. Relevant test cases are considered to demonstrate the capability and accuracy of the new model in dealing with two-phase flows at wide range of density contrasts. These include the Rayleigh–Taylor instability at two distinct density ratios, the breakup of a turbulent liquid column, and assisted breakup of a liquid layer. The validated model is then utilized to examine a liquid sheet emerging from a nozzle into a stagnant gas at different Reynolds numbers ranging from 10 0 0 to 60 0 0. The influence of jet injection velocity on the flow structure, variation of the jet centerline velocity, and jet area is also investigated. This work reports for the first time the extension of the Smagorinsky model for LES to the binary fluid model of Fakhari & Lee [1] and its application in simulation of turbulent liquid sheet breakup at the density ratios as high as 250. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Disintegration of liquid jets is a fundamental topic of fluid mechanics and has numerous engineering applications. Among them the breakup of liquid fuels in combustion systems of gas turbines is of great importance, since combustion efficiency is directly influenced by the atomization of the fuel jet [2]. Therefore, to improve the design of these devices better understanding of the breakup process through experimental investigations and numerical simulations have to be provided. Liquid jet breakup is a result of an intricate process that includes the development of instabilities along the gas-liquid interface, formation of liquid ligaments and their subsequent stretching and break-up into dispersed droplets of various sizes. The breakup is potentially affected by a variety of factors: cavitation inside the nozzle, gas compressibility, aerodynamic forces, surface tension, and turbulence [3,4]. Due to this complexity, one frequently encounters technical barriers to perform accurate experimental measurements. Thus, development of numerical methods is a powerful complementary approach to obtain a deeper insight into the physics of liquid jet disintegration. ∗

Corresponding author. E-mail address: [email protected] (M.H. Rahimian).

https://doi.org/10.1016/j.compfluid.2017.10.023 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

Over the past decade, numerous studies have been devoted to establishing accurate numerical schemes to simulate the liquid jets breakup. Most of these studies adopt the Direct Numerical Simulation (DNS) approach for turbulence modeling along with front tracking [5] or front capturing methods [6,7] to describe the interface dynamics. Klein [8] studied a spatially developing liquid film exhausting into gaseous atmosphere in two and three dimensions, using direct numerical simulation (DNS) and a Volume of Fluid (VOF) [9]. Results qualitatively validated against available experimental data. The author addressed several important issues in his study. It was argued that both 2D and 3D simulations result in an identical optimum wavelength (most amplified wavelength on the interface) which is a key parameter to predict jet breakup. In addition, it was stated that except the growth rate of surface waves and shear stresses, all other remaining quantities (such as turbulence statistics) at least qualitatively follow the same trends for both 2D and 3D simulations. Furthermore, it was mentioned that the VOF method has its own disadvantages. These include the uncertainties in the interface shape and thus in calculation of the surface tension forces and relatively high computational costs especially in 3D. To alleviate these issues, different interface tracking/capturing approaches have been proposed including the study of Menard et al. [10] on simulation of the primary breakup using the coupled VOF/Level Set (LS)/ghost fluid (GF), Herrmann’s

94

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

method of interface tracking using an auxiliary high resolution grid [11], the conservative LS method of Olsson et al. [12], and the DNS of liquid jet atomization by Lebas et al. [13]. It is worth mentioning that utilization of DNS provides an indepth insight into the fundamental physics and mechanisms of liquid jet atomization. However, the extreme computational costs of DNS approach, prevents its application for most engineering problems. Evidently, the high computational costs of DNS studies originate from the multitude of spatial and temporal scales present in two-phase turbulent systems. Note that, for turbulent multi-phase flows the computational domain should be sufficiently large to resolve the most energetic length scales which are on the order of the nozzle diameter. Simultaneously, the resolution should be high enough to capture the smallest dispersed liquid droplets. Therefore, the use of DNS in practical applications due to the excessive computational demands is not feasible. This drawback is somewhat relieved in Large Eddy Simulation (LES) approaches, where the smallest length scales are ignored via low pass filtering of the Navier–Stokes equations [14]. While LES models have been widely used in single phase flows, its application in simulation of turbulent flows in presence of interfaces is relatively new. Moreover, a majority of the studies was concerned with large eddy simulation of separated flows such as breaking waves [15,16] and gas-liquid bubbly flows [17,18], while the subject of LES of liquid jet disintegration is rather unexplored. The relevant studies found in the literature can be summarized as the pioneering work of De Villier et al. [19] on LES of diesel spray, Chesnel et al.’s work [20] on subgrid scale (SGS) analysis of liquid jet atomization, Desjardins et al. [21] study on LES/LS/GF simulation of the turbulent liquid jet breakup, and Xiao et al. [22] investigation on air-blast atomization of round water jets using the LS/VOF/LES methodology. The phase interface in the latter study was directly resolved without modeling the SGS terms, while the regular LES filter was applied in the single phase regions. Besides, the SGS stress tensor was modeled via the standard Smagorinsky approach [23]. From these references we conclude that the simulation of liquid jet disintegration using Navier–Stokes-based solvers is still demanding. The main problems include the laborious numerical procedures, prohibitive computational costs, sophisticated treatment of discontinuities at the gas-liquid interface and numerical stability issues for large density ratios. Over the last few years, the LBM has been considered as a powerful alternative to classical methods to solve the Navier–Stokes equations and has attracted attention in various disciplines of computational fluid dynamics [24]. Since the LBM is based on mesoscopic kinetic equations, it can be considered as a promising approach to handle systems in which an interface is present [25]. A variety of multiphase LB models have been proposed to date [26– 30]. Obviously, the preliminary versions suffered from numerous drawbacks including numerical instability for large density ratios and the presence of parasitic currents at the interface. These issues, however, are alleviated to a great extent in the model proposed by Lee and Lin [31]. Using a two-distribution discrete Boltzmann equation as well as an improved representation of the surface tension force, they enabled the LBE formulation to treat problems with a density ratio of up to 10 0 0. Furthermore they introduced the potential form of forcing terms along with a compact isotropic discretization to control the interface thickness and to avoid an overly diffusive interface as well as a reduction of spurious currents. With the increasing popularity of LB methods, turbulent flow simulations have been successfully carried out for single phase flows [32–41]. It has been demonstrated that LES within the LBM framework not only can accurately describe the turbulence hydrodynamics, but also has distinct computational advantages over its traditional continuum counterparts [42]. In spite of these interest-

ing efforts, a comprehensive LB model capable of simulating turbulent multiphase flows is still lacking. The only relevant studies are those of Fakhari and Lee [1] and Banari et al. [43,44]. In the first study the binary-fluid LB model of Lee [45] was extended to high Reynolds number flows using the multiple-relaxation-time (MRT) collision operator [46]. It was shown that the MRT technique has a superior performance compared to the single-relaxation-time model in terms of numerical stability at high Reynolds numbers. However, since no subgrid scale model is employed in Fakhari and Lee’s formulation [1], fully resolved simulations capturing all contributing scales (from small scales up to the largest structures) are required. Therefore, to reduce the prohibitive computational costs, it is advisable to only resolve the large structures and account for the missing small scales by employing an appropriate subgrid closure model, in a similar manner to LES of single phase flows. A rare example is the work of Banari et al. [43,44], where the effect of sub-scale terms has been taken into account by modifying the molecular viscosity in calculation of the viscous volume forces. This approach which is based on Inamuro et al.’s work [47], while promising has its own drawbacks; it is computationally less efficient since it requires solving a Poisson equation for pressure correction, parasitic currents are present at the interface, and it impairs the inherent simplicity of the LBM, as argued by Lee and Lin [31]. In the present paper we address the development of a suitable SGS closure approach for binary immiscible fluids within the framework of phase field lattice Boltzmann methods. The proposed method is able to simulate the moderate and even high Re flows, leads to stable solutions at high density and viscosity ratios and does not require solving an additional pressure correction equation. To this end, we incorporate the Smagorinsky subgrid closure [23] into the binary-fluid model of Lee [45]. The philosophy behind the developed formulation is similar to that proposed by Xiao et al. [22]; the interface is simulated without any SGS modeling and spatial filtering is applied in the single phase regions of both gas and liquid. Using the filtered LB equations of both the flow field and the order parameter, we simulate the primary breakup of a turbulent liquid sheet emerging from a nozzle into a stagnant gas at transitional Re numbers ranging from 10 0 0 to 60 0 0 and liquidto-gas density ratios up to 250. Note that, to the best of the authors’ knowledge, no study deals with the simulation of the turbulent liquid jets using LBM and the published researches are only related to the Lattice Boltzmann simulation of single phase jets [24,33,48]. Therefore, the current paper is the first reported application of LBM to simulation of turbulent liquid sheets breakup. The rest of the paper is organized as follows: In the next section the theoretical framework is described and details of the spatially filtered LB equations as well as the employed subgrid closure are given. In Section 3, the implemented model is assessed and validated through relevant test cases. In Section 4 the simulation results for a liquid sheet introduced into the quiescent gas are given. A summary of the obtained results are presented in Section 5. 2. Mathematical models and governing equations As mentioned earlier, the main idea behind the LES is to resolve the large scale eddies and mimic the presence of small scale structures via the modeling of subgrid scales [49]. In this method to simulate the resolved scales, a space filtering is applied to the different quantities. The filtering operation is defined as:

w¯ (x ) =



w(x )G(x, x )dx ,

(1)

where w represents the physical quantity under consideration and G is the filter function. Note that, various choices for the filter function are available, however as recommended in [32], we use a box

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

filter defined as [50]:







G x, x =

|xi − x i | < i /2 , otherwise

1 / i 0

(2)

with i being the filter width in direction i. Once the filter function is selected, the first step is to apply the filtering operation to the velocity, pressure, composition and also the particle distribution functions. The filtering decomposes the instantaneous velocity field u(x, t) into two components: the spatially filtered velocity u¯ (x, t), and the fluctuating component u (x, t):

u(x, t ) = u¯ (x, t ) + u (x, t )

(3)

Furthermore, in the current approach (similar to Xiao et al. [22]) spatial filtering is not applied to the kinematic viscosity and density and these quantities are obtained from the properties of one or both phases depending on the local value of the filtered composition C¯:

ρ









= C¯ ρl + 1 − C¯

ρg

(4)

ν = C¯νl + 1 − C¯ νg

(5)

Note that, the temporal and spatial evolution of the phase interface and flow field are determined by applying the filtering operation to the transport equations of the composition and momentum. The modified governing equations are presented in the following sections. 2.1. Filtered Cahn–Hilliard and phase field Navier–Stokes equations The numerical framework utilized in the current study implements a phase field approach as a class of diffuse interface models where the interface has a thin but nonzero thickness. Across this interface the composition or the order parameter (C) continuously varies from the value of 1 in the liquid phase to 0 in the gas phase. In the proposed model of Lee, the transport of composition is modeled as [45]:

∂C + ∇ · (uC ) = ∇ · (M∇μ ), ∂t

(6)

where u is the macroscopic velocity, M is the mobility factor, μ denotes the chemical potential and M∇μ describes the volume diffusive flow rate. In order to determine the mobility factor (M), we fix the M × β at 0.04, where the β is explicitly given by [51]:

β=

12σ , D

(7)

with σ being the surface tension and D the interface thickness. Note that other choices for M × β are also possible. However, based on our numerical tests the value of 0.04 was found to be optimal with respect to accuracy and stability. Now, we apply a space filtering operation box filter [50] to Eq. (6):

∂C + ∇ · (uC ) = ∇ · (M∇μ ), ∂t

95

total subgrid scale terms of the LES transport equation, is very small. Thus, one can simply omit this term without considerable loss of accuracy. This is equivalent to assume that the composition is primarily transported by the resolved velocity through the domain and the role of SGS velocity is insignificant. Finally, the filtered Cahn–Hilliard equation can be written as:

  ∂ C¯ + ∇ · u¯ C¯ = ∇ · (M∇ μ ¯ ), ∂t

(11)

where μ ¯ is the filtered chemical potential. It is important to note that, as a consequence of the above mentioned assumptions, the resolved Cahn–Hilliard equation cannot track the interfacial distortions due to under-resolved dynamic features. In other words, the filtered composition forms the dominant part of interface geometry and the contribution of distortions due to the small scale eddies disappears. Thus, in the above equation the filtered chemical μ¯ is given by

μ¯ = μ¯ 0 − κ∇ 2C¯,

(12)

where κ is the gradient parameter and μ ¯ 0 is the filtered classical part of the chemical potential which is given by μ ¯0 = 2β C¯ (C¯ − 1 )(2C¯ − 1 ). In addition to the Cahn-Hilliard equation that describes the evolution of the phase interface; extra equations are needed to compute the flow field. Considering both liquid and gas as incompressible, the phase field modified momentum equation that includes the surface forces can be written as [52]:



ρ

   ∂u + u.∇ u = −∇ p + ∇ · μ ∇ u + ∇ uT + μ∇ C, ∂t

(13)

where p is the pressure. Applying the filtering operation to momentum equation, using the eddy viscosity type models for subgrid closure, and neglecting the SGS terms arising from the filtering of the surface forces, the spatially filtered Navier–Stokes becomes:



   ∂ u¯ ρ + u¯ .∇ u¯ = −∇ p¯ + ∇ · μ ∇ u¯ + ∇ u¯ T + μ ¯ ∇ C¯ − ρ∇ .τSGS , ∂t (14) where τSGS = uu − u u = −2νt S¯ is the residual stress tensor with S¯ = ∇ u¯ + ∇ u¯ T being the large scale strain rate tensor and ν t being the subgrid kinematic eddy viscosity. 2.2. Filtered form of discrete boltzmann equations Considering the intermolecular interactions and non-ideal gas effects, the Discrete Boltzmann equation (DBE) for the transport of density and momentum with a generalized collision operator  is written as [53]

D fα = Dt

  ∂ fα 1 + eα · ∇ fα = − fα − fαeq + (eα − u ) · F fαeq ∂t ρ cs2

(8)

(15)

recalling the fact that the incompressibility hypothesis holds if operator (. )is used, i.e.∇ · u¯ = 0, and rearranging Eq. (8), one can write:

where fα is the particle distribution function, eα is the microscopic particle velocity, ρ is the density, u is the volume averaged veloceq ity. In addition fα , the equilibrium distribution function, is represented by

  ∂ C¯ + ∇ · u¯ C¯ + τC = ∇ · (M∇ μ ¯ ), ∂t

(9)

where τ C is the interfacial subgrid term and is represented by:

τC = ∇ · (uC ) − u¯ .∇ C¯.

(10)

Note that, as elaborated by Chesnel et al [20], except for the atomization area the contribution of the interfacial terms to the



fαeq = wα ρ 1 +



ea · u (ea · u )2 (u · u ) + − , 2 cs 2cs4 2cs2

(16)

with wα being the weight factor [54]. In Eq. (15), the intermolecular forcing (F) is described by

F = ∇ ρ cs2 − ∇ p + μ∇ C,

(17)

96

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

wherep is the hydrodynamic pressure. Note that in above equation the −C ∇μ term is replaced with μ∇ C expression. As explained by Lee and Connington [55], either versions are valid. However, based on the authors experience the new version is much more stable for large density ratios. In above relations, the collision operator is considered as [1]

 = (M−1

SM − 2I )/δt,

(18)

where I is the identity matrix. In addition, M and Sˆ are the orthogonal transformation matrix and the diagonal relaxation matrix, respectively and are given by [1]

⎡

1 ⎢−4 ⎢0 ⎢ ⎢0 ⎢ M=⎢ 0 ⎢0 ⎢ ⎢0 ⎣ 0 0

⎡

0 ⎢0 ⎢0 ⎢ ⎢0 ⎢

S = ⎢0 ⎢0 ⎢ ⎢0 ⎣ 0 0

1 −1 −2 1 −2 0 0 1 0 0 1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0

1 −1 −2 0 0 1 −2 −1 0 0 0 0 0 0 0 0 0 0

1 −1 −2 −1 2 0 0 1 0 0 0 0 0 1.7 0 0 0 0

1 −1 −2 0 0 −1 2 −1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1.7 0 0

1 2 1 1 1 1 1 0 1

1 2 1 −1 −1 1 1 0 −1 0 0 0 0 0 0 0 s8 1

1 2 1 −1 −1 −1 −1 0 1

⎤

0 0⎥ 0⎥ ⎥ 0⎥ ⎥ 0 ⎥, 0⎥ ⎥ 0⎥ ⎦ 0 s8

⎤

1 2⎥ 1⎥ ⎥ 1⎥ ⎥ 1 ⎥, (19) −1⎥ ⎥ −1⎥ ⎦ 0 −1

(20)

cs2 . ν + 0.5cs2

   νt = Cs 2 S = Cs 2 2SS,

1

(21)

(27)

with Cs and  being the Smagorinsky constant and filter width, respectively. Note that, tuning the Smagorinsky constant in an appropriate way is crucial. Since the present approach is based on the use of the spatially filtered DBE in the single phase regions, we employ the recommended value for LBE-LES of single phase jets, i.e., Cs = 0.1 [56]. Furthermore, the filter width is the grid spacing which is normally set to unity in LBM simulations. As pointed out earlier, another critical problem in solving Eq. (22) arises from the fact that there is no simple method to close the second term on the right hand side, the intermolecular forcing term. A feasible way could be to neglect the impact of subgrid scale on the intermolecular forcing term. First of all it is important to notice that the intermolecular forcing term comprises the difference between the DBE for single phase flows and that of binary fluids. Thus, the forcing term is expected to be zero throughout the flow field, except close to the interface. Meanwhile, the filtered interface evolution equation, Eq. (11), already lacks the contribution of SGS velocity and basically is unable to capture the small scale fluctuations of the interface. As a consequence, omitting the subgrid terms in the forcing term which is only present close to the interface, will not introduce additional errors to our model. We thus propose to approximate the filtered form of the intermolecular forcing term by

ρ cs2

in which s8 is related to the kinematic viscosity by

s8 =

where ν t is the turbulent viscosity. In the current study to model the unresolved scales we implement the standard Smagorinsky approach [23], in which the turbulent viscosity is calculated via,

(eα − u ) · F fαeq ≈

1

ρ cs2

(eα − u¯ ) · F¯ fαeq (ρ , u¯ ) with,

(28)

¯ ∇ C¯ . F¯ = ∇ ρ cs2 − ∇ p¯ + μ

(29)

Applying the filtering operation to Eq. (15), the coarse-grained, f¯α , satisfies:

In summary, the filtered Discrete Boltzmann equation (DBE) for transport of density and momentum for binary fluids can be expressed as

  D fα 1 = − fα − fαeq + (eα − u ) · F fαeq . Dt ρ cs2

  D fα 1 = −∗ fα − fαeq (ρ , u¯ ) + (eα − u¯ ) · F¯ fαeq (ρ , u¯ ). Dt ρ cs2

(22)

There are two major issues in developing a solution for Eq. (22): the filtered non-linear collision term needs to be closed; and the non-linear filtered forcing term should be approximated. Note that, in classical approaches it is a common practice to take into account the effect of unresolved scales by adding the SGS viscosity ν SGS to the viscos stress. The similar idea has been implemented within the single phase LES-LB framework by modifying the relaxation time or matrix [32,33,37]. Here, we utilize an identical method and modify the relaxation matrix to mimic the presence under-resolved scales, so that:

     fα − fαeq = ∗ fα − fαeq (ρ , u¯ ) ,

(23)

where fα (ρ , u¯ ) has the same formulation as Eq. (16), except that the velocity vector is replaced by the filtered velocity. Note that, in eq general, the fα (ρ , u¯ ) is not equal to the filtered equilibrium distribution function. In the above equation, Ʌ∗ is the modified collision operator and is determined by eq

−1 ∗

 = (M S M − 2I )/δt, ∗

(24)

with Sˆ ∗ being the modified diagonal relaxation matrix, in which the matrix elements are identical to Eq. (20), except that the molecular viscosity ν is replaced by the efficient viscosity ν eff :

s8 =

νe f f

cs2 , + 0.5cs2

νe f f = ν + νt ,

(25) (26)

(30)

It can be proven that the Chapman–Enskog expansion of Eq. (30) recovers the spatially filtered nearly incompressible Navier–Stokes equation in the presence of an interface, i.e., Eq. (14). As discussed in [57], in order to transform the DBE for the evolution of the density and momentum to the DBE for the evolution of pressure and momentum, a new distribution function (g¯ α ) is introduced;





¯ α (0 ) with g¯ α = f¯α cs2 + p¯ − ρ cs2

(31)

 e .u¯ (e .u¯ )2 u¯ .u¯ ¯ α (u ) = fαeq (ρ , u¯ )/ρ = ωα 1 + α2 + α 4 − 2 . cs

2 cs

2 cs

(32)

Now by taking the total derivative of g¯ α and assuming an incompressible flow, dropping the terms with higher order than Ma2 and finally integrating using the trapezoidal method, the following LBE is obtained:





 g¯¯α (x + eα δt, t + δt ) − g¯¯α (x, t ) = −(S + 2I ) g¯¯α − g¯¯eq α (x,t )  M 2   ¯α +δt (eα − u¯ ) · ∇ ρ cs ¯ α − ¯ α (0 ) + μ∇ ¯ MC¯ , (x,t )

(33)

in which S = ∗ δt, I is the identity matrix, ∇ M denotes the mixed difference approximation, and g¯¯ α and g¯¯ α are the modified distribution function and modified equilibrium distribution function, respectively, defined as:

 δt S g¯ α − g¯eq g¯¯ α = g¯ α + (eα − u¯ ) α − 2 2

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

·



¯eq

   ∇ C ρ cs2 ¯ α − ¯ α (0 ) + μ∇ CC¯ ¯ α . δt

eq

g¯ α = g¯ α −

2



(eα − u¯ ) · ∇ ρ C

cs2



(34)

  ¯ α − ¯ α (0 ) + μ∇ CC¯ ¯ α with, (35)





eq ¯ α (0 ), ¯ )cs2 + p¯ − ρ cs2 g¯eq α = f α (ρ , u

(36)

∇C

where stands for the central difference approximation. Details about the gradient operators can be found in [31,45]. To determine the evolution of the interface, an additional particle distribution function needs to be defined. The second distribution function should recover the filtered Cahn–Hilliard equation, Eq. (11),









h¯ α = C¯/ρ f¯α ,

(37)

eq ¯ ). h¯ eq α = C¯ /ρ f α (ρ , u

Taking the Eq. (11) yields:

total

(38) derivative

of

the

h¯ α

and

utilizing

 1 ∂ h¯ α ¯) + eα .∇ h¯ α = − h¯ α − h¯ eq α + ( eα − u ∂t λ     C¯  ¯ α + M∇ 2 μ ¯ · ∇ C¯ − − μ∇ ¯ ¯ α , p C ∇ 2 ρ cs (39) with λ being the single relaxation time. Applying the trapezoidal integration, employing a modified form of particle distribution, defining τ = λ/δt, and doing some simplifications the LBE for composition evolution can be written as:

   τ (x,t )      ¯ C + δt (eα − u¯ ) · ∇ MC¯ − ¯ MC¯ α  ∇ M p¯ − μ∇ ρ cs2

h¯¯ α (x + eα δt, t + δt ) − h¯¯ α (x, t ) = −

+

δt  2



1 h¯¯ α − h¯¯ eq α + 0.5

    δt  M∇ 2 μ ¯ α  + M∇ 2 μ ¯ α  (x,t ) ( x+e 2

(x,t )

α δt ,t +δt )

,

(40)

eq where h¯¯ α , h¯¯ α are defined as:





 ¯ ¯ eq  ¯h¯ = h¯ + hα − hα − δt e − u¯ · ∇ C C¯ − C¯ ∇ C p¯ − μ∇ C ¯ ¯ C α . ( α ) α α 2τ 2 ρ cs2 (41) ¯ eq δt (eα − u¯ ) · h¯¯ eq α = hα − 2

  C¯  ∇ CC¯ − 2 ∇ C p¯ − μ∇ ¯ C C¯ α . ρ cs

(42)

In solving the above LBE, the single relaxation time, λ, is fixed at 0.5δ t. 2.3. Estimation of macroscopic properties of flow

¯ h¯ α ,

(43)

α

ρ u¯ = p¯ =

1  ¯ δt eα g¯α + μ∇ ¯ C¯, 2 cs2 α

 α

g¯¯α +

δt 2

u¯ · ∇ρ cs2 .

3. Model validation As mentioned earlier various researches are dedicated to understanding the mechanisms leading to liquid jet breakup and atomization. Lasheras and Hopfinger [58] elaborated that depending on the jet Re number either the surface tension or Kelvin-Helmholtz (KH) instabilities have the prominent role in amplification of interface perturbations and formation of longitudinal waves. Further, it has been argued [21] that the Rayleigh–Taylor type instabilities are responsible for the formation of bulges at the top of the waves crest and subsequent destabilization of liquid jets. Consequently, to demonstrate the capability of the implemented phase field model in capturing of interfacial instabilities, the oscillation of a liquid drop induced by the initial perturbation, RTI and KHI benchmarks seem to be appropriate. Note that, to the authors’ knowledge, there is no reported simulation of Rayleigh-Taylor instability at high density ratios in presence of surface tension using LB method. Thus, the details of this test case are given in the following section. In contrast, the oscillating drop benchmark as well as the KH instability have been previously well studied by the authors [30]. Hence, for brevity we avoid repeating the details and results of these verifications. Additionally, to collect more evidences on validity of the numerical framework in simulation of liquid sheet breakup we compare the predicted mean velocity profiles for a turbulent infinitely long liquid column with the previously published data. Furthermore, using the developed model, assisted breakup of a liquid layer is examined and the obtained potential core length is validated against the empirical data. 3.1. Rayleigh–Taylor instability The RTI is considered to evaluate the accuracy of the method in calculating the interface dynamics for a finite density ratio. Two distinct cases are studied: In the first case, the density difference between two fluids is relatively small and the location of interface tip is checked with the results reported in [59]. In the second case, the density ratio between the heavier and lighter phase is app. ∼ 10 0 0, and the growth rate of interface disturbances is verified against the analytical solution presented by Schott et al. [60]. Note that in both cases the surface tension is taken into account. Thus, these test cases also confirm the ability of the model in dealing with large density jumps across the phase interface in the presence of surface tension. In the first benchmark, the computational domain is a W × 4 W rectangular box, with W being the domain width and is covered by 256 × 1024 grid points. The domain is occupied with two immiscible fluids, where initially the denser fluid (ρ l = 1) is located on top and the lighter phase whose density is ρ g = 0.143 occupies the bottom of the domain. The interface of two fluids is initially perturbed using a sinusoidal wave with the amplitude of 5% over the system width. Both upper and lower boundaries are considered as freeslip, while the side walls are periodic. For the entiresimulations gravity is pointing downward and is selected so that W g = 0.04.



The macroscopic quantities of the flow including the composition, density, velocity, and pressure is calculated as follows [45]:

C¯ =

97

(44)

(45)

Additionally, similar to [59] the Reynolds number (Re = W gW/ν ) is 1225.19, the viscosity ratio is set to unity, and the Atwood number, A = (ρL − ρg )/(ρL + ρg ), is 0.76. In presenting the results the



time is normalized by the T ∗ = W/g. Fig. 1 shows the temporal evolution of the interface shape. In this type of instability, the denser phase gradually penetrates into the light fluid due to gravity and an umbrella-like interface is formed. As seen, this behavior is well captured by the implemented model and the obtained results compare quantitatively well with those reported previously [38,59]. Fig. 2 presents the comparison between the location of the bubble front and spike tip

98

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

Fig. 1. Evolution of the interface shape with time. The system size is 256 × 1024, Re = 2048 and A = 0.5.

Fig. 2. Comparison with Terashima [59] and study on time histories of nondimensionalized position of bubble front and spike tip. The system size is 256 × 1024, Re = 1225.19 and A = 0.76.

(normalized by the system width: W) with those of Terashima and Tryggvason [59]. As can be seen, we obtain a quantitative agreement. It is worth noting that Terashima and Tryggvason’s study [59] is conducted at a density ratio of 7. Thus, to collect more evidence on reliability of the current numerical framework for high density ratios, we compare our simulation to the analytical data reported in [60]. The configuration of this case is somehow similar to the former one, however differs from the first case in several aspects. In the study of Schott et al. [60] Re = 10 0 0, ρ L /ρ g = 10 0 0 and hence A = 0.995, and the dynamic viscosity ratio (μL /μg ) is 10 0 0. Further, the boundary condition on the top and down is changed to no-slip. The interface initial perturbation is defined by −a0 cos(2π x/W ), where a0 denotes the disturbance amplitude which is set to be 1% of the characteristic length (e.g. system width).

Fig. 3. Comparison of the non-dimensional growth rate obtained from simulation with those reported by Schott et al. [60] for different dimension-less wave numbers. The system size is 1024 × 4048, Re = 10 0 0, ρ L /ρ g = 10 0 0, μL /μg = 10 0 0, and A = 0.995.

Once the numerical setup is determined, we run the simulations to examine the dependence of dimensionless growth factor on non-dimensional wave number and compare the results with those predicted by linear theory [60]. In presenting the results, the maximum interface displacement is normalized with the domain width (W), the wavenumber is specified as k = 2π /W and scaled as k∗ = k(ν 2 /g)1/3 , and the dimensionless growth factor is defined by:

α∗ =

ln(an /a0 ) 1   T 3 g2 /ν

(46)

where, ν = νl = νg is the kinematic viscosity, T is the simulation time and an is the amplitude of interface displacement. Fig. 3 displays the growth of initial disturbances at different wavenumbers.

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

99

3.3. Assisted breakup of a liquid sheet As the final test case we consider the simulation of the assisted breakup which is characterized by the liquid injection into a highspeed co-flowing gas. We aim to match the predicted potential cone length with those obtained from the experiments [62]. The flow configuration is shown in Fig. 7, similar to that reported by Rayana et al. [62]. The liquid sheet enters the domain through an inlet with the height of Hl with a bulk velocity of Ul , while the gas layer has a thickness of Hg and velocity of Ug . The liquid and high speed gas streams are separated with a plate of thickness θ . The inflow velocity profile for both streams is determined by:

y ux (y ) = U0 tanh( ),

(47)

δ

where δ is the vorticity of the gas stream; δ = 6Hg /Reg 0.5 , with Reg = Ug Hg /νg being the gas Reynolds number, y is the vertical co is specified for each stream such that ordinate, and U 0

Ul =

Fig. 4. Schematic representation of the temporally evolving liquid sheet.

The graph is compared with linear theory presented in Schott et al. [60] and obviously the agreement is very good. 3.2. Infinite liquid sheet breakup As another test case, the disintegration of an infinitely long turbulent liquid sheet is investigated. Here we consider a thin liquid sheet of thickness H inside a gaseous medium, as shown in Fig. 4. The surface tension is present at the interfaces and both liquid and gas phases considered as incompressible. The computational setup consists of a rectangular box of size [0,4H] × [0,6H]. The results are verified against those presented by Desjardins [61] using a combined LS/GF approach. The symmetric boundary condition is imposed to the top and bottom boundaries with periodic boundaries at the left and right sides of the domain. The sheet has an initial bulk velocity of U0 , which is obtained from a precursor periodic channel flow simulation. Simulations are carried out on a domain of size 800 × 1200 for the density ratios (ρ l /ρ g ) of 40, viscosity ratios (μl /μg ) of 40, Re = 30 0 0.0 (Re = ρl U0 H/μl ), and Weber numbers (We = ρl U20 H/σ ) of 500 and 2000. Note that, the time in dimensionless form is defined by t = T U0 /H. Incorporating these data to our simulation, time evolution of the interface shape for the Weber number of 20 0 0 is depicted in Fig. 5. As seen by marching in time, the tangential velocity difference triggers the modulations across the material interface. Amplification of these disturbances, leads to shrinking and widening in different parts of the liquid sheet and its final disintegration into ligaments. Fig. 6 presents a comparison between the mean profiles of streamwise velocity at t = 30 obtained from the simulations with those reported by Desjardins [61] at two different Weber numbers of 500 and 2000. In this figure the y-coordinate is normalized with jet half width (δ U ), which is defined as the distance from the jet centerline where the axial velocity falls below the half jet centerline velocity (Ucl ). The mean velocity profile is obtained by averaging the velocity profile at three different cross sections: x = H, 2H, and 3H. As seen in this figure, the predicted trends are similar, although the results do not exactly coincide.

1 Hl

Hl

ux (y )dy,

Ug =

0

1 Hg

Hl + θ +Hg

ux (y )dy.

(48)

Hl +θ

The simulations are carried out on a domain of size 573 × 267 (x × y) for the density ratios (ρ l /ρ g ) of 100 and viscosity ratios (μl /μg ) of 20. The separating plate thickness, θ , is set to 4. A summary of simulation parameters is presented in Table 1. Note that, the time in dimensionless form is defined by t = TUg /Hg . Snapshots of the simulation results at two different momentum ratios (M = ρgUg2 /ρl Ul2 ) of M = 4 and M = 16, are displayed in Fig. 8. The shear between high-speed gas and low-speed liquid streams and KH type instabilities lead to formation of longitudinal waves at the interface. Further amplification of the interface perturbations is followed by the development of liquid ligaments that gradually rise into the gas flow. These ligaments are stretched and break into dispersed liquid droplets resulting in the dilution of the liquid phase. This phenomena occurs within a certain distance from the nozzle outlet which is called the breakup length or liquid potential cone length [63]. It was argued that the normalized potential cone length is only a function of momentum ratio M according to the relation [64]:

Lc 12 = √ . Hl M

(49)

In order to attain the potential cone length we use a similar approach to [21], and sample the values of

H (x, t ) = x



Ly 0

C (x, y, t )dy,

(50)

over space and time, where x is the grid spacing, C is the order parameter, and Ly is the domain height. Note that, the order parameter C can be also regarded as the liquid volume fraction and therefore Integration of this quantity in the vertical direction can provide a good measure of interface height and its evolution over time. A comparison between the predicted values at two different momentum ratios of 4 and 16 and those obtained from Eq. (49) is presented in Fig. (9). As seen the predicted trends are consistent with experimental data. 4. Results and discussion In the previous section, several classical two-phase flow test cases as well as the breakup of an infinite liquid sheet were studied. The results confirm the capability of the modified phase field model in accurate capturing of the interfacial instabilities. Therefore, one can deduce that the numerical framework is also suitable for simulations of liquid jet breakup. The flow configuration

100

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

Fig. 5. Time Evolution of the liquid sheet for ρl /ρg = 40, Re = 30 0 0, We = 20 0 0. Table 1 A summary of simulation parameters for LES of assisted breakup of a liquid sheet. M (ρgUg2 /ρl Ul2 )

Rel (Ul Hl /ν l )

Reg (Ug Hg /ν g )

Wel (ρl Ul2 Hl /σ )

Weg (ρgUg2 Hg /σ )

4 16

1250 1250

50 0 0 10,0 0 0

4.53 4.53

18.1 72.46

is shown in Fig. 10, which consists of a two-dimensional planar liquid jet exhausting through a nozzle into the stagnant gas. The domain size is defined based on nozzle height (D), assumed to be H = 10D and L = 14D in normal and streamwise directions, respectively. On the lateral direction, open boundaries are specified. Such a choice allows the flow to stream freely from and into the domain, minimizing the effect of side walls on the jet during its spreading phase. In streamwise direction, the right boundary is treated as convective boundary condition in a similar procedure elaborated by Lou et al. [65]. On the left boundary, assuming that the x and y components of velocity as well as the composition are known (e.g. ux = uinlet , uy = 0, ρ = 1.0, and C = 1.0), the unknown distribution functions are replaced according to the method presented in [66] which has been derived for single phase flows. As in the inlet boundary only the liquid phase enters into the domain, that approach can be simply adapted to our framework. On the remaining boundaries the no slip wall is applied. Before starting the simulations, the sufficient spatial resolution should be determined. Here, to find the proper grid size, we examine the sensitivity of

the calculated penetration length to the system size. As listed in Table 2 we start from the coarsest grid (Grid I), in which the nozzle height is resolved with 100 lattices, and then successively increase the system size by a factor of 1.5 in each direction to obtain higher resolutions (i.e. Grid II, III, and IV). To perform the grid independency test, we choose four dimensionless groups that describe the flow physics:

ρLU0 Ds μL μL Liquid Jet Ohnesorge number : Oh =  ρL D s σ μL Dynamic viscosity ratio : κ= μg ρL Density ratio : φ= ρg Liquid Jet Reynolds number :

Re =

(51) (52) (53) (54)

The simulation parameters are ρl /ρg = 250, μl /μg = 27, Re = 20 0 0 and Oh = 0.0258. Results are presented in Fig. 11. As

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

101

Fig. 8. Snapshots of interface shape at two different momentum ratios: a) M = 4, b) M = 16.

seen beyond the Grid III, further increasing the system size has no substantial effect on the calculated distance. Hence, in the rest of this section, we perform the simulations on a fixed grid size of 112 × 225 and 3150 × 2250 in the nozzle and chamber, respectively. Note that in performing grid sensitivity analysis, it was alternatively possible to use other parameters such as the shape of the interface. However, based on the simulation results it is observed that the interface shape is slightly grid dependent. Similar finding was previously reported in the study of He et al. [57], in which the shape of spikes in Rayleigh–Taylor instability was indicated to be somewhat grid dependent. Now that the sufficient spatial resolution is provided, the breakup of liquid jets can be examined in details. The dimensionless quantities are presented in Table 3. As seen, for all cases, the density ratio is ρl /ρg = 100, the viscosity ratio is μl /μg = 10, and Ohnesorge number (oh) is 0.0258. Further, for all present simulations the gas-liquid interface is assumed to be spread over 4 lattice cells (i.e. D = 4). Fig. 6. Comparison of mean streamwise velocity profile with those reported by Desjardins [61] at t = 30 for Re = 30 0 0 and We numbers of 500 and 2000.

Fig. 7. Schematic representation of the examined geometry.

102

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

Fig. 9. Comparison of normalized H(x, t) with Eq. (49), dash-dotted line, at two different momentum ratios: a) M = 4, b) M = 16.

Fig. 12 shows a liquid jet emanating into quiescent gas at the density ratio of 100 and Re = 6000 (Case 2 of Table 3). The evolution of the interface is depicted at different time steps. Here, the evolution time is normalized as t ∗ = tU0 /L, where L is the domain length and t is the time steps elapsed after simulation starts. As seen, at the early stages of simulation, due to impact of the moving liquid to the still gas an umbrella-like head shape is developed. Shinjo and Umermura [67] argued that from the jet head, rotating vortices are shed into the upstream and a strong recirculation zone is formed behind the jet head. The presence of this region with non-zero axial velocity and dominant disturbance leads to the escalation of the interface modulations, which primarily appear as a result of the shear forces between the rapid moving liquid and its surrounding gas. As the time passes the oscillations of the in-

terface are amplified, that results in formation of ligaments which are then dilated and finally break up into droplets. 4.1. Influence of Re number on jet disintegration Various test cases are simulated to determine the effect of jet Re number on liquid sheet breakup. For all cases, dimensionless numbers are identical except for the Re number which is changed by altering the inflow velocity. The jet inflow velocities are 0.01, 0.03, 0.04, 0.05 (in lattice units) and their corresponding Re numbers are 10 0 0, 30 0 0, 40 0 0, and 50 0 0, respectively. Simulation parameters are summarized in Table 3, cases 2–5. Fig. 13 depicts the instantaneous snapshots of liquid interface at different Re numbers. Further, the velocity field obtained by averaging the com-

Fig. 10. Schematic representation of the liquid sheet simulation.

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

103

Fig. 13, these disturbances are then amplified in streamwise direction that cause higher spreading rates further downstream. Obviously, these dynamical features become stronger by increasing the Re number, which promote the interface perturbations and consequently hasten the liquid disintegration. Fig. 14 depicts the evolution of jet centerline velocity along the streamwise direction at different Re numbers. As seen, for all cases up to around X/D = 8, the velocity almost maintains its initial value, which resembles the self-similar state of the single phase jets at the far field. Note that, for the free turbulent single phase jets, under certain conditions the velocity reaches a self-preserving shape and the evolution of the centerline velocity can be expressed by [8]:

 U 2 0

Ucl

Fig. 11. Penetration length versus time at four different spatial resolutions for ρl /ρg = 250, Re = 20 0 0, Oh = 0.0258, and μl /μg = 27. Table 2 Grid Independency. Grid

System Size (i × j)

I

Nozzle Chamber

50 × 100 1400 × 1000

II

Nozzle Chamber

75 × 150 2100 × 1500

III

Nozzle Chamber

112 × 225 3150 × 2250

IV

Nozzle Chamber

168 × 337 4718 × 3370

Table 3 Presentation of examined test cases (Note: quantities are in lattice units).

νl

ID Case Case Case Case Case

1 2 3 4 5

0.0019 0.0023 0.0023 0.0023 0.0023

νg 0.0188 0.0225 0.0225 0.0225 0.0225

σ −5

9.37 × 10 8.44 × 10−5 9.37 × 10−5 9.37 × 10−5 9.37 × 10−5

U0

ρ l /ρ g

Re

Oh

0.05 0.01 0.03 0.04 0.05

100 100 100 100 100

60 0 0 10 0 0 30 0 0 40 0 0 50 0 0

0.0258 0.0258 0.0258 0.0258 0.0258

ponents of velocity vector over a time period of 5 (i.e. t∗ = 5) is also shown in the same figure. As seen, the jet breakup appears different; with the cases with higher Re number the liquid core being much more perturbed, and the liquid column being fully shrouded by the droplets. This is a clear indication of the aerodynamic forces that play a significant role in the atomization of liquid jets [68]. In fact, by increasing the inlet velocity we necessarily increase the relative velocity between liquid and gas phases. Higher relative velocities lead to stronger Kelvin–Helmholtz instabilities that amplify the waves on the surface of the jet which in turn enhances the disintegration of the liquid column. Additionally, the liquid jet breakup mechanism can be seen as the transformation of kinetic energy at inflow to surface energy [69]. In other words, for the jet with higher inlet velocity due to a growing amount of kinetic energy more intense interface distortions will be achieved. Alternatively, the interface perturbation can be related to the unrolling of vortices emanating from shear layer. The vortices trigger the first interface modulations at upstream. As shown in

=A

X D



−B ,

(55)

where Ucl is the jet centerline velocity, and A and B are constants. However, except for the case of Re = 10 0 0, beyond the X/D = 8 the variation of the second power of the (U0 /Ucl ) with the normalized distance deviates from Eq. (55). This can be related to the fact that downstream of the potential core region, the dissipation becomes faster, the breakup is accelerated and hence the liquid column is easily deflected in transverse direction, which leads to diminished jet centerline velocity. The decrease in centerline velocity is more evident at higher Re numbers, which is a clear sign of promoted liquid jet breakup. In contrast, for Re = 10 0 0, the potential core length falls outside the computational domain and thus the jet retains its nozzle exit velocity within the depicted range. A closer scrutiny of this figure also reveals that the decay of the centerline velocity from the nozzle exit value occurs earlier at higher Re numbers, indicating that the potential core length varies inversely as the Re number. We close the discussion in this section by providing a quantitative analysis of the effect of Re number on liquid jet disintegration. Here, we examine the temporal evolution of liquid sheet surface by measuring the amount of liquid surface within the computational domain. The computed surface is then normalized by the unperturbed liquid column surface (similar to Fig. 13, Re = 10 0 0) which essentially is A0 = LD. The results are presented in Fig. 15. As seen in the figure, the trends confirm that the liquid area is a function of Re number. Since the liquid area can be regarded as a measure of spreading rate, one may conclude that stronger liquid breakup can be achieved at higher Re numbers. 5. Conclusion A phase field lattice Boltzmann model, capable of handling large density contrasts has been extended to LES of two-phase turbulent flows. To this end, the standard Smagorinsky subgrid-scale model has been introduced into the MRT-LB method for immiscible fluids [1]. Once the capability of the developed LB scheme in capturing of two important instabilities was checked, the model has been employed to examine the breakup of liquid sheets exhausting through a nozzle into the still gas. Further, the influence of the inlet velocity on disintegration pattern, jet centerline velocity and jet spreading area were investigated. The results indicated that formation of stronger vortical features within the shear layer at higher jet inflow velocities causes promoted interface distortions. It was shown that these amplified perturbations that grow further downstream lead to higher jet spreading area, stronger disintegration of the liquid sheet and thus better mixing rates. Besides, regarding the results it was evident that the jet potential core length is reduced at higher Re numbers. To conclude, although the results are encouraging, there are still some remaining issues to be alleviated in future. Most notable of these, are the limited range of parameters, high computational

104

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

Fig. 12. Temporal evolution of interface location for case 1 with Re = 60 0 0. From top the non dimensionalized time, t∗ is 0.25, 0.5, 0.75, 1.0, and 1.25, respectively.

costs, lack of mass conservation, numerical diffusion at the material interface, and the lack of a model for the turbulence/interface interaction and the inability to track the interfacial distortions due to under-resolved dynamic features. In our experience, the minimum gas viscosity achievable in computations is limited by the occurrence of instabilities especially at high density ratios. This issue restricts the implementation of the model in engineering problems where the Re number, density ratio and viscosity ratio are typically large. Furthermore, in our simulations we have observed the issue of droplet dissolution. This problem arises from both the lack of mass conservation in the Lee’s model as shown in [70] and the diffuse nature of the method. It is worth mentioning that we have ex-

tended our numerical code to the mass-conservative version of the Lee’s model [71] very recently. Our preliminary simulations show that not only the issue of drop dissolution is relieved to some extent, but also the new model has a superior performance in terms of efficiency and computational speed-up. It is also important to mention that, not only the present approach but also all other developed formulations within the framework of Navier–Stokes equations, confront with profound difficulties when consider the LES of two phase flows in the presence of the interface. This problem arises from the lack of a generalized model to take into account the influence of subgrid flow properties on the interfacial deformations, regardless of the flow configuration. One possible way is

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

105

Fig. 13. Spatial evolution of liquid sheet at t∗ = 3 (Cases 2–5) with different injection velocities and Re numbers (left) - Temporal averaged flow field after t∗ = 5 (right).

Fig. 14. Effect of Re number on the evolution of the jet centerline velocity (Cases 2–5).

Fig. 15. Effect of Re number on temporal evolution of the dimensionless liquid sheet surface within the computational domain (Cases 2–5).

106

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107

to modify the volumetric surface tension force to include the subgrid terms due to two phase flow. Finally, we acknowledge that both the turbulence phenomena and the liquid jets breakup are 3D in nature. Thus, to obtain the results that accurately represent the flow features it is inevitable to perform the simulations in 3D. Extending the simulations to 3D, use of adaptive mesh refinement techniques such as [27] that fully resolve the interface and keep the computational demands low, and implementation of parallel computing algorithms will be the subject of future studies. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.compfluid.2017.10.023. References [1] Fakhari A, Lee T. Multiple-relaxation-time lattice Boltzmann method for immiscible fluids at high Reynolds numbers. Phys Rev E 2013;87(2):023304. [2] Arienti M, Soteriou MC. Time-resolved proper orthogonal decomposition of liquid jet dynamics. Phys Fluids (1994-present) 2009;21(11):112104. [3] Hsiao FB, Huang JM. On the evolution of instabilities in the near field of a plane jet. Phys Fluids A(1989-1993) 1990;2(3):400–12. [4] Ashgriz N. Handbook of atomization and sprays. Springer; 2010. [5] Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys 1992;100(1):25–37. [6] Gueyffier D, Li J, Nadim A, Scardovelli R, Zaleski S. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows. J Comput Phys 1999;152(2):423–56. [7] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114(1):146–59. [8] Klein M. Direct numerical simulation of a spatially developing water sheet at moderate Reynolds number. Int J Heat Fluid Flow 2005;26(5):722–31. [9] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39(1):201–25. [10] Ménard T, Tanguy S, Berlemont A. Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet. Int J Multiphase Flow 2007;33(5):510–24. [11] Herrmann M. A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids. J Comput Phys 2008;227(4):2674–706. [12] Olsson E, Kreiss G, Zahedi S. A conservative level set method for two phase flow II. J Comput Phys 2007;225(1):785–807. [13] Lebas R, Menard T, Beau P, Berlemont A, Demoulin F-X. Numerical simulation of primary break-up and atomization: DNS and modelling study. Int J Multiphase Flow 2009;35(3):247–60. [14] Pope SB. Turbulent flows. IOP Publishing; 2001. [15] Hieu PD, Katsutoshi T, Ca VT. Numerical simulation of breaking waves using a two-phase flow model. Appl Math Modell 20 04;28(11):983–10 05. [16] Lubin P, Vincent S, Abadie S, Caltagirone J-P. Three-dimensional large eddy simulation of air entrainment under plunging breaking waves. Coastal Eng 2006;53(8):631–55. [17] Deen NG, Solberg T, Hjertager BrH. Large eddy simulation of the gas–liquid flow in a square cross-sectioned bubble column. Chem Eng Sci 2001;56(21):6341–9. [18] Liovic P, Lakehal D. Interface–turbulence interactions in large-scale bubbling processes. Int J Heat Fluid Flow 2007;28(1):127–44. [19] De Villiers E, Gosman A, and Weller H, 2004, Large eddy simulation of primary diesel spray atomization, No. 0148-7191, SAE Technical Paper. [20] Chesnel J, Menard T, Reveillon J, Demoulin F-X. Subgrid analysis of liquid jet atomization. Atomization Sprays 2011;21(1). [21] Desjardins O, McCaslin J, Owkes M, Brady P. Direct numerical and large-eddy simulation of primary atomization in complex geometries. Atomization Sprays 2013;23(11). [22] Xiao F, Dianat M, McGuirk JJ. LES of turbulent liquid jet primary breakup in turbulent coaxial air flow. Int J Multiphase Flow 2014;60:103–18. [23] Smagorinsky J. General circulation experiments with the primitive equations: I. the basic experiment∗ . Mon Weather Rev 1963;91(3):99–164. [24] Geller S, Uphoff S, Krafczyk M. Turbulent jet computations based on MRT and Cascaded Lattice Boltzmann models. Comput Math Appl 2013;65(12):1956–66. [25] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Ann Rev Fluid Mech 1998;30(1):329–64. [26] Aidun CK, Clausen JR. Lattice-Boltzmann method for complex flows. Ann Rev Fluid Mech 2010;42:439–72. [27] Tölke J, Freudiger S, Krafczyk M. An adaptive scheme for LBE multiphase flow simulations on hierarchical grids. Comput Fluids 2006;35(8-9):820–30. [28] Safari H, Rahimian MH, Krafczyk M. Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow. Phys Rev E 2013;88(1):013304. [29] Safari H, Rahimian MH, Krafczyk M. Consistent simulation of droplet evaporation based on the phase-field multiphase lattice Boltzmann method. Phys Rev E 2014;90(3):033305.

[30] Amirshaghaghi H, Rahimian M, Safari H. Application of a two phase lattice Boltzmann model in simulation of free surface jet impingement heat transfer. International communications in heat and mass transfer; 2016. [31] Lee T, Lin C-L. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J Comput Phys 2005;206(1):16–47. [32] Hou S, Sterling J, Chen S, and Doolen G, 1994, A lattice Boltzmann subgrid model for high Reynolds number flows, arXiv preprint comp-gas/9401004. [33] Yu H, Girimaji SS. Near-field turbulent simulations of rectangular jets using lattice Boltzmann method. Phys Fluids (1994-present) 2005;17(12):125106. [34] Girimaji SS. Boltzmann kinetic equation for filtered fluid turbulence. Phys Rev Lett 2007;99(3):034501. [35] Jafari S, Rahnama M. Shear-improved Smagorinsky modeling of turbulent channel flow using generalized Lattice Boltzmann equation. Int J Numer Methods Fluids 2011;67(6):700–12. [36] Malaspinas O, Sagaut P. Consistent subgrid scale modelling for lattice Boltzmann methods. J Fluid Mech 2012;700:514–42. [37] Krafczyk M, Tölke J, Luo L-S. Large-eddy simulations with a multiple-relaxation-time LBE model. Int J Modern Phys B 2003;17(01n02):33–9. [38] Chen S, Tölke J, Krafczyk M. Simple lattice Boltzmann subgrid-scale model for convectional flows with high Rayleigh numbers within an enclosed circular annular cavity. Phys Rev E 2009;80(2):026702. [39] Stiebler M, Krafczyk M, Freudiger S, Geier M. Lattice Boltzmann large eddy simulation of subcritical flows around a sphere on non-uniform grids. Comput Math Appl 2011;61(12):3475–84. [40] Geier M, Schönherr M, Pasquali A, Krafczyk M. The cumulant lattice Boltzmann equation in three dimensions: Theory and validation. Comput Math Appl 2015;70(4):507–47. [41] Pasquali A, Schonherr M, Geier M, Krafczyk M, Simulation of external aerodynamics of the DrivAer model with the LBM on GPGPUs, 2015. [42] Sagaut P. Toward advanced subgrid models for Lattice-Boltzmann-based Large-eddy simulation: theoretical formulations. Comput Math Appl 2010;59(7):2194–9. [43] Banari A, Janßen C, Grilli ST, Krafczyk M. Efficient GPGPU implementation of a lattice Boltzmann model for multiphase flows with high density ratios. Comput Fluids 2014;93:1–17. [44] Banari A, Janßen CF, Grilli ST. An efficient lattice Boltzmann multiphase model for 3D flows with large density ratios at high Reynolds numbers. Comput Math Appl 2014;68(12):1819–43. [45] Lee T, Liu L. Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces. J comput phys 2010;229(20):8045–63. [46] dHumieres D, Ginzburg I, Krafczyk M, Lallemand P, Luo L-S, Bushnell DM, Multiple-relaxation-time Lattice Boltzmann models in 3D, 2002. [47] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput phys 2004;198(2):628–44. [48] Yang Y-T, Chang S-C, Chiou C-S. Lattice Boltzmann method and large-eddy simulation for turbulent impinging jet cooling. Int J Heat Mass Transfer 2013;61:543–53. [49] Sagaut P. Large eddy simulation for incompressible flows: an introduction. Springer Science & Business Media; 2006. [50] Ferziger JH, 1981, Higher-level simulations of turbulent flows. [51] Lee T. Effects of incompressibility on the elimination of parasitic currents in the lattice Boltzmann equation method for binary fluids. Comput Math Appl 2009;58(5):987–94. [52] Gurtin ME, Polignone D, Viñals J. Two-phase binary fluids and immiscible fluids described by an order parameter. Math Models Methods Appl. Sci 1996;6(06):815–31. [53] He X, Shan X, Doolen GD. Discrete Boltzmann equation model for nonideal gases. Phys Rev E 1998;57(1):R13. [54] He X, Luo L-S. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Phys Rev E 1997;56(6):6811. [55] Connington K, Lee T. Lattice Boltzmann simulations of forced wetting transitions of drops on superhydrophobic surfaces. J Comput Phys 2013;250:601–15. [56] Yu H, Luo L-S, Girimaji SS. LES of turbulent square jet flow using an MRT lattice Boltzmann model. Comput Fluids 2006;35(8):957–65. [57] He X, Chen S, Zhang R. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability. J Comput Phys 1999;152(2):642–63. [58] Lasheras J, Hopfinger E. Liquid jet instability and atomization in a coaxial gas stream. Ann Rev Fluid Mech 20 0 0;32(1):275–308. [59] Terashima H, Tryggvason G. A front-tracking/ghost-fluid method for fluid interfaces in compressible flows. J Comput Phys 2009;228(11):4012–37. [60] Schott B, Rasthofer U, Gravemeier V, Wall W. A face-oriented stabilized Nitsche-type extended variational multiscale method for incompressible two-phase flow. Int J Numer Methods Eng 2015;104(7):721–48. [61] Desjardins O. Numerical methods for liquid atomization and application in detailed simulations of a diesel jet. Stanford University; 2008. [62] Rayana FB, Cartellier A, Hopfinger E, Assisted atomization of a liquid layer: investigation of the parameters affecting the mean drop size prediction, Proc. proceedings of the international conference on liquid atomization and spray systems (ICLASS), Kyoto, Japan. [63] Villermaux E. Mixing and spray formation in coaxial jets. J Propul Power 1998;14(5). [64] Marty S, 2015, Contribution à l’étude de l’atomisation assistée d’un liquide: instabilité de cisaillement et génération du spray, Université Grenoble Alpes.

H. Amirshaghaghi et al. / Computers and Fluids 160 (2018) 93–107 [65] Lou Q, Guo Z, Shi B. Evaluation of outflow boundary conditions for two-phase lattice Boltzmann equation. Phys Rev E 2013;87(6):063301. [66] Mohamad AA. Lattice boltzmann method: fundamentals and engineering applications with computer codes. Springer Science & Business Media; 2011. [67] Shinjo J, Umemura A. Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation. Int J Multiphase Flow 2010;36(7):513–32. [68] Wu P-K, Faeth G. Aerodynamic effects on primary breakup of turbulent liquids. Atomization Sprays 1993;3(3).

107

[69] Sander W, Weigand B. Direct numerical simulation and analysis of instability enhancing parameters in liquid sheets at moderate reynolds numbers. Phys Fluids (1994-present) 2008;20(5):053301. [70] Guo Z, Zheng C, Shi B. Force imbalance in lattice Boltzmann equation for two-phase flows. Phys Rev E 2011;83(3):036707. [71] Geier M, Fakhari A, Lee T. Conservative phase-field lattice Boltzmann model for interface tracking equation. Phys Rev E 2015;91(6):063309.