Multiscale simulations and experiments on water jet atomization

Accepted Manuscript

Multiscale simulations and experiments on water jet atomization Mahdi Saeedipour, Simon Schneiderbauer, Gregor Plohl, Gunter Brenn, Stefan Pirker ¨ PII: DOI: Reference:

S0301-9322(16)30745-5 10.1016/j.ijmultiphaseflow.2017.05.006 IJMF 2596

To appear in:

International Journal of Multiphase Flow

Received date: Revised date: Accepted date:

12 December 2016 28 March 2017 19 May 2017

Please cite this article as: Mahdi Saeedipour, Simon Schneiderbauer, Gregor Plohl, Gunter Brenn, ¨ Stefan Pirker, Multiscale simulations and experiments on water jet atomization, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.05.006

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We present a RANS-based multiscale model for simulation of primary

CR IP T

atomization.

• Phase-Doppler anemometry (PDA) experiments were done to characterize water sprays.

• The model is validated against experimental data using a point-by-

AN US

point comparison.

• The model predicts the size and velocity distributions of droplets reasonably.

AC

CE

PT

ED

high-Re sprays.

M

• The model requires relatively low computational effort for modeling

1

ACCEPTED MANUSCRIPT

Multiscale simulations and experiments on water jet atomization

CR IP T

Mahdi Saeedipoura,b,∗, Simon Schneiderbauera,b , Gregor Plohlc , G¨ unter Brennc , Stefan Pirkera a

AN US

Department of Particulate Flow Modelling, Johannes Kepler University, 4040 Linz, Austria. b Christian Doppler Laboratory for Multi-scale Modeling of Multiphase Processes, Johannes Kepler University, 4040 Linz, Austria. c Institute of Fluid Mechanics and Heat Transfer, Graz University of Technology, 8010 Graz, Austria.

Abstract

Numerical simulation of primary atomization at high Reynolds number is still

M

a challenging problem. In this work a multiscale approach for the numerical simulation of liquid jet primary atomization is applied, using an Eulerian-

ED

Lagrangian coupling. In this approach, an Eulerian volume of fluid (VOF) method, where the Reynolds stresses are closed by a Reynolds stress model

PT

is applied to model the global spreading of the liquid jet. The formation of the micro-scale droplets, which are usually smaller than the grid spacing in

CE

the computational domain, is modelled by an energy-based sub-grid model. Where the disruptive forces (turbulence and surface pressure) of turbulent eddies near the surface of the jet overcome the capillary forces, droplets are

AC

released with the local properties of the corresponding eddies. The dynamics of the generated droplets are modelled using Lagrangian particle tracking ∗

Corresponding author Email address: [email protected] (Mahdi Saeedipour)

Preprint submitted to International Journal of Multiphase Flow

May 19, 2017

ACCEPTED MANUSCRIPT

(LPT). A numerical coupling between the Eulerian and Lagrangian frames is then established via source terms in conservation equations. As a follow-up

CR IP T

study to our investigation in Saeedipour et al. (2016a), the present paper aims at modelling drop formation from liquid jets at high Reynolds numbers in the atomization regime and validating the simulation results against in-house

experiments. For this purpose, phase-Doppler anemometry (PDA) was used to measure the droplet size and velocity distributions in sprays produced by

AN US

water jet breakup at different Reynolds numbers in the atomization regime. The spray properties, such as droplet size spectra, local and global Sautermean drop sizes and velocity distributions obtained from the simulations are compared with experiment at various locations with very good agreement. Keywords: Jet breakup, primary atomization, multiscale modelling, VOF,

ED

1. Introduction

M

LPT, phase-Doppler anemometry

The breakup of a liquid jet is an important step in many technical and

PT

industrial processes, from spray drying to fuel injection, from ink-jet printing to surface treatment. The jet breakup results in a stream or spray of

CE

droplets with sizes depending on the breakup mechanism. Both the mean drop size and the width of the size spectrum depend on that mechanism.

AC

For a given nozzle geometry and state of the ambient medium, v. Ohnesorge (1936) identified the pair of characteristic numbers Re = ρU D0 /µ

and Oh = µ/(σρD0 )1/2 (jet velocity U , nozzle orifice diameter D0 , liquid density and dynamic viscosity ρ and µ, liquid surface tension against the ambient medium σ) as relevant for identifying the breakup mechanism. At 3

ACCEPTED MANUSCRIPT

low values of the Reynolds number Re < 48.4 Oh−0.794 , the jet breaks up due to its capillary (Rayleigh) instability, and the most frequent drop size

CR IP T

in the size spectrum is about two times the nozzle orifice diameter. At intermediate Reynolds numbers 48.4 Oh−0.794 < Re < 228.5 Oh−0.814 , the aerodynamic interaction of the jet with the ambient medium determines the breakup (Kelvin-Helmholtz instability), and the most frequent drop size is

of the order of the nozzle orifice diameter. Finally, at high Reynolds number

AN US

Re > 228.5 Oh−0.814 , the jet breaks up due to atomization, and the drops are much smaller than the orifice diameter. While the drop size resulting from capillary jet breakup may be predicted quite accurately with analytical methods (Lafrance, 1975; Yuen, 1968), this is not possible for the windinduced Kelvin-Helmholtz mechanism (Gordillo and P´erez-Saborid, 2005),

M

and even less for atomization. Numerical simulations of the two latter processes are difficult for well known reasons but allow to predict the breakup

ED

length and/or the size of the emerging drops. Presently we concentrate on the atomization regime. The primary breakup

PT

represents the initiation of the disintegration of the liquid jet. This process is then followed by secondary breakup (i.e. breakup of larger droplets into

CE

smaller ones) towards the full atomization process. Despite the importance of primary breakup for various real-world applications, the physics behind the whole process is still not fully understood. Several factors have been identi-

AC

fied to play a role in the primary breakup process, such as liquid-gas turbulence, nozzle cavitation, velocity profile relaxation and liquid-gas properties (Srinivasan et al., 2008). The phenomena involved in liquid jet atomization exhibit different length and time scales. For instance, in addition to micro-

4

ACCEPTED MANUSCRIPT

scale velocity fluctuations, which disrupt the jet surface initiating atomization, micro-scale interfacial topological changes as well as the surface tension

CR IP T

forces play a significant role in jet breakup. Especially, the latter surface tension forces considerably impact the liquid jet breakup in the micro-scale

(Herrmann, 2013). The dependency of atomization on turbulence and topol-

ogy of the interface (i.e. the liquid surface curvature) make the numerical simulation of primary breakup still a computationally challenging problem.

AN US

In other words, proper methodologies for treating the turbulence as well as resolving the small scale liquid-gas interfacial curvature, are needed for the numerical simulation of this process. This, however, entails the requirement of high grid resolution as the confinement. Over the years, a number of

efforts have been done to model liquid jet atomization with various meth-

M

ods, mostly based on solving the Navier-Stokes equations coupled with a proper interface capturing method such as volume-of-fluid (VOF) (Hirt and

ED

Nichols, 1981; Fuster et al., 2009; Ling et al., 2015), level set (LS) (Osher and Sethian, 1988; Desjardins et al., 2008; Herrmann, 2011) and combina-

PT

tions of these methods (Sussman and Puckett, 2000; M´enard et al., 2007; Lebas et al., 2009). Due to the well-known shortcomings of each method

CE

(Tryggvason et al., 2011), finding the most appropriate choice of modelling strategy, in particular for the atomization problem, is still an active field of research.

AC

Direct numerical simulation (DNS) (M´enard et al., 2006; M´enard et al.,

2007; Gorokhovski and Herrmann, 2008; Shinjo and Umemura, 2010) and large eddy simulation (LES) (Chesnel et al., 2011; Xiao, F., 2012; Desjardins et al., 2013; Pavlovic et al., 2014) of primary atomization have been investi-

5

ACCEPTED MANUSCRIPT

gated in the past and ongoing literature. The main advantage of the DNS methods is the ability to capture all liquid structures (e.g. interface defor-

CR IP T

mations, ligaments and droplets) by the computational grid with a high level of accuracy. Also in the context of LES, the large turbulent liquid structures are resolved by the grid accurately. However, these methods require

high grid resolutions and advanced numerical schemes which lead to high computational costs. To overcome the dependency with respect to compu-

AN US

tational resources, multiscale modelling strategies have been developed in recent years. They involve mesoscopic particles or droplets suspended in macroscopic flow fields (Tomar et al., 2010). In these frameworks, numerical simulation is performed by either Eulerian-Eulerian or Eulerian-Lagrangian approaches. The former resolves both scales in an Eulerian formulation, while

M

the latter includes computing the continuous phase in an Eulerian scheme, and model the discrete phase using the Lagrangian formulation. Herrmann

ED

(2005) developed a coupling method for the primary breakup of a liquid jet in which all interface dynamics and physical processes occurring on scales larger

PT

than the available grid resolution are fully resolved and all dynamics and processes occurring on sub-grid scales are modeled in a Lagrangian framework

CE

by using a level set/vortex sheet method. Ham et al. (2003) developed a hybrid Eulerian-Lagrangian technique in which they model primary breakup using the LES and a surface capturing method and account for the secondary

AC

breakup within the Lagrangian framework. Tomar et al. (2010) performed a multiscale simulation of primary jet breakup using an Eulerian-Lagrangian two-way coupling. Grosshans et al. (2014) combined the VOF method with LPT to simulate jet atomization in a hybrid framework. In this approach,

6

ACCEPTED MANUSCRIPT

an LES-VOF simulation is used until a statistical convergence point. The rest of the simulation is then done using LPT. Recently, Ling et al. (2015)

CR IP T

proposed a multiscale model for the simulation of atomizing flows, where all interfacial length scales are resolved by a DNS-VOF approach, and the small resolved droplets are then converted to the Lagrangian point-particle (LPP)

model to improve the numerical stability and computational efficiency. In this DNS-based approach, the coupling between VOF and LPP is done by

AN US

considering a local spatial filtering.

The present study follows the authors’ previous works (Saeedipour et al., 2014, 2016b,a) in multiscale modelling of liquid jet breakup for industrial applications. Unlike the fully resolved DNS or a highly resolved LES in modelling liquid atomization, which demand high grid resolutions (i.e. ∆x ∼

M

O(10−6 − 10−5 ) m) and face limitations to moderate Reynolds numbers, the

present RANS-based approach is able to model primary atomization at high

ED

Reynolds number jets (i.e. Re ∼ O(104 - 105 )) due to the use of coarser

grid spacing (∆x ∼ O(10−4 ) m). This method predicts global properties of

PT

atomizing jets, such as the Sauter mean diameter (SMD) as well as droplet size and velocity distributions and the general behavior of spray propagation

CE

with appropriate accuracy using relatively coarse grid. It is worth mentioning that this method does not intend to resolve the small-scale interfacial topological changes leading to primary breakup. Instead their effect is mod-

AC

elled using a statistical representation of the flow field. Therefore, a sub-grid breakup model is applied to bridge statistical properties obtained from the meso-scale interface capturing with a time-averaged turbulence model and deterministic properties of the droplets in the micro-scale. This bridging be-

7

ACCEPTED MANUSCRIPT

tween scales is established based on an energy-based criterion depending on the flow properties in the sub-grid scale. Accordingly, the multiscale mod-

CR IP T

elling strategy is developed by an Eulerian-Lagrangian coupling in which the general spreading of the liquid jet is simulated by an Eulerian approach

(i.e. the VOF method), while the small-scale droplets are tracked within a Lagrangian framework (i.e. LPT).

Early experimental literature on liquid jet breakup concentrated on iden-

AN US

tifying the mechanisms leading to the destabilization of the jet. One classical work is the above referenced study by v. Ohnesorge (1936). Today there exists a wide literature about spray formation by liquid jet breakup, including experimental studies on the properties of the disperse droplet phase in the sprays. Experimental data on the drop size and velocity spectra are mostly

M

acquired with phase-Doppler anemometry (PDA). Studies providing information about spatially average drop sizes may rely on laser diffraction-based

ED

measurements with instruments of the Malvern type. Experimental investigations of jet atomization mostly relate to Diesel injection. Information

PT

about the drop size spectrum in such sprays is difficult to get by optical means since the sprays are very dense. As a consequence, for the purpose of

CE

modelling jet breakup in Diesel injection, such as realized by Yi and Reitz (2004), drop sizes are often determined at the periphery of the sprays only (Wu et al., 1986). A comprehensive review of jet breakup experiments in all

AC

the breakup regimes was given by Dumouchel (2008). Our paper is structured as follows. In the next section we give a physical

description of the liquid atomization process and present the numerical model in Sec. 3. Section 4 presents the experiments on liquid jet atomization, and

8

ACCEPTED MANUSCRIPT

the results are compared with the simulations in Sec. 5. The paper ends with the conclusions in Sec. 6.

CR IP T

2. Physical description

Although the macroscopic behavior of liquid jet breakup can be described

by instability theory, the microscopic interfacial effects including the deformation of the liquid-gas interface, formation of liquid structures (ligaments)

AN US

and consequently, the breakup into small droplets must be studied at smaller scales. Therefore, liquid jet breakup is considered as an inherently multiscale problem both in time and length scales. Resolving the breakup process from large towards small scales (i.e. instability of the interface and formation of ligaments breaking up into small droplets) is computationally expensive and

M

typically considered unfeasible for industrial applications. As mentioned in the introduction, a coarse-grid simulation which, as a compromise, covers the

ED

most effective phenomenon on the larger scale, could be a choice for modelling. Therefore, a multiscale method is required to consider micro-scale

PT

processes in a coarse-grid domain. This concept is depicted schematically in Fig. 1. In this figure, the case with finer grid consists of several small-scale

CE

interfacial topological changes which cannot be captured on the coarser grid. Besides this, the small liquid structures, i.e. the droplets, are still captured

AC

within the Eulerian framework (i.e. VOF), while on the coarser grid they belong to the sub-grid scale, and the Eulerian framework is not capable of modelling their dynamics. In this context, a sub-grid model is needed in order to account for the smaller droplets, which are not resolved by the coarse Eulerian grid. 9

ACCEPTED MANUSCRIPT

An in-depth understanding of the physics on the small scale is essential for defining a sub-grid model. According to theoretical and experimental

CR IP T

works carried out by Wu and Faeth (1993) and by Faeth et al. (1995) primary breakup of turbulent liquid-gas multiphase flows has two major effects: interfacial turbulence and aerodynamic effects. Thus, both the onset of drop

formation and drop properties after turbulent primary breakup can be studied by treating interactions between these disruptive forces against surface

AN US

tension. Accordingly, the energy required for droplet formation includes the

kinetic energy of an eddy of characteristic size relative to its surroundings plus the added mechanical energy due to the pressure drop caused by acceleration of the surrounding gas over the tip of the perturbation (i.e. aerodynamic lift). If these energies can overcome the stabilizing surface energies, droplets

M

can be formed. This is schematically shown in Fig. 2. In this figure, li is the turbulent eddy size, v 0 the associated velocity fluctuation, and it is assumed

ED

that the liquid structure moves at the bulk jet velocity U¯ . At the interface, the specific turbulent kinetic energy (Ek ∼ ρl v 02 ), specific mechanical energy

PT

due to the aerodynamic lift (Ea ∼ ρg U¯ 2 ) and specific surface tension energy

(Es ∼ σ/li ) are major players in turbulent breakup. In these relations, ρl , ρg

CE

and σ denote the liquid density, gas density and surface tension coefficient of the liquid, respectively. The energy balance between these disruptive and stabilizing energies determines whether breakup occurs or not. If disruptive

AC

energies (Ek + Ea ) exceed the surface energy, then droplets can be formed, which can be expressed as the following energy-based criterion: Ek + Ea ≥ Es

(1)

In flows showing a density ratio greater than 500 the breakup process is 10

ACCEPTED MANUSCRIPT

mainly controlled by liquid turbulence and the aerodynamic effects appear negligible. This especially applies to the present study, where a liquid water

CR IP T

jet is surrounded by air unveiling a density ratio of approximately 1000. For liquid-gas density ratios less than 500, however, aerodynamic effects are significant and the time scale of primary breakup of droplets near the surface

of the jet may become large compared to the time scale needed for secondary breakup. Then these two breakup modes can merge. In other words, the

AN US

droplet immediately breaks up into smaller droplets after formation, which however, is not the case in the present study. This theory is thouroughly validated in a DNS study by Desjardins and Pitsch (2010) and stands out as a central element of turbulent primary breakup.

In order to translate this physical description to the definition of a sub-

M

grid model, an appropriate estimation of the unresolved characteristic length scale li is needed. First, it should be emphasized that in the case of large

ED

density ratios (i.e. > 500) and in the case of very small eddies (i.e. in the size of the Kolmogorov length scale), the turbulent kinetic energy may not

PT

be sufficient to overcome the surface energy (Wu and Faeth, 1993) and consequently, no droplet is formed. This, in turn, implies that inequality (1) is

CE

not satisfied for very small eddies in the flow and an appropriate length scale should, therefore, be large enough to contain sufficient kinetic energy for the droplet formation. Second, according to Faeth et al. (1995) the characteristic

AC

length must be smaller than the largest eddies in the flow. Consequently, it is proposed by Wu and Faeth (1993) that this length scale is in the inertial sub-range of the turbulence. On the basis of turbulence energy cascade theory (Pope, 2000) shown in Fig. 3, the Taylor micro-scale λ is, therefore the

11

ACCEPTED MANUSCRIPT

most appropriate measure for li (Saeedipour et al., 2016a). Then this energy-based sub-grid breakup criterion is reformulated as (2)

CR IP T

(ρl vλ2 + Csa ρg U¯ 2 )λ3 ≥ Csi σλ2

where vλ is the velocity fluctuation perpendicular to the jet surface and Csa

and Csi are two empirical constants. Since these constants are of order unity (Wu and Faeth, 1993), they were set equal to one in the present study. The

AN US

Taylor micro-scale is modeled with the root mean square of the velocity

fluctuations hv 0 irms , kinematic viscosity ν and turbulence dissipation rate  as follows (Pope, 2000):

li ' λ =

r

ν 15 hv0 irms ε

(3)

M

It is worth mentioning that the droplets formed at the onset of primary breakup are the smallest droplets that can be formed by this energy-based

ED

mechanism (Faeth et al., 1995). This is beneficial for validation purposes because it explains the increase of droplet sizes with increasing the streamwise

PT

distance from the nozzle exit (Wu et al., 1992). This, in turn, implies that the model provides the information about the minimum droplet size that can

CE

be obtained by the simulation. Thus it can be used to evaluate the limitation of the model in predicting the droplet size distribution compared to the real

AC

spray.

3. Numerical model Liquid jet atomization is described by the incompressible Navier-Stokes

equations for the two-phase immiscible liquid-gas system, together with a 12

ACCEPTED MANUSCRIPT

transport equation for capturing the interface as presented in Eqs. (4a) to (4c). (4a)

CR IP T

~ =0 ∇·U

(4b)

∂α ~ = Sα + ∇ · (αU) ∂t

(4c)

AN US

~ ∂(ρU) ~ ⊗ U) ~ = ∇ · τ¯ + ρ~g + Sσ + Smom + ∇ · (ρU ∂t

~ is the mixture velocity field shared by In this single-fluid formulation, U both phases. The density and molecular viscosity of the fluid are defined as ρ = αρl + (1 − α)ρg and µ = αµl + (1 − α)µg , where α is the volume

M

fraction of the liquid phase in the two-phase mixture. The surface tension force Sσ is treated by the Continuous Surface Force (CSF) model by Brackbill

ED

et al. (1992) using the filtered volume fraction. This method describes the surface force due to capillarity as a volumetric source term in the momentum

PT

equation acting only on the interfacial areas and reads Sσ = σκ∇α

(6)

AC

CE

where κ is the interface curvature   ∇α κ=∇· | ∇α |

(5)

Further to this, Smom and Sα are additional source terms which typi-

cally do not appear in the system of equations for incompressible two-phase flows, however, in this work they may become non-zero in case of sub-grid modelling. 13

ACCEPTED MANUSCRIPT

By using a very high grid resolution and applying correct computational schemes, especially for the surface tension term, this set of equations can

CR IP T

resolve all the liquid interfacial structures in the breakup process (e.g. ligaments, droplets etc.), but this is computationally demanding. The main objective of this study is the numerical simulation of jet breakup using a sub-grid model with less computational costs.

The sub-grid breakup model should in fact bridge time-averaged prop-

AN US

erties of the flow from the turbulence model to deterministic properties of droplets, such as size and velocity. This bridge is established by (i) evaluation of the sub-grid criterion based on the energy balance and (ii) defining a local droplet diameter using a characteristic turbulent length scale. The former is handled by inequality (2) which evaluates the energy balance between

M

disruptive and consolidating energies on the sub-grid. The latter is done by equation (3) in which the local droplet diameter can be estimated by the

ED

Taylor micro-scale defined on the basis of statistical turbulence properties (e.g. k,  and ν). Both of these steps result in an average local diame-

PT

ter for each droplet which gets tracked within the Lagrangian framework. Therefore, once these steps are fulfilled in any cell containing the liquid-gas

CE

interface, a specific number of droplets (N ) are released at that position of the computational domain with respect to the local properties of the cell. The dynamics of the generated droplets are then governed by Lagrangian

AC

equations of motion. This coupling approach is schematically shown in Fig. 4.

We propose a model to calculate N as a function of turbulence frequency

f , eddy size λ and grid spacing ∆x as follows:

14

ACCEPTED MANUSCRIPT



N˙ = f

∆x λ

2

(7)

CR IP T

where f can be estimated as the ratio of the magnitude of the sub-grid velocity fluctuation vector and the turbulent characteristic length scale λ f=

~0 | |u λ

(8)

Then the droplets are locally injected perpendicular to the jet surface, and

∇α |∇α|

  ~ (init) = U ~ + u ~0 · ~n ~n U d

(9)

~ and u ~0 are the cell is the interface normal vector and U

M

where ~n =

AN US

their initial local velocity vector upon release reads

velocity and velocity fluctuation vectors, respectively.

ED

In this RANS-based model, a relatively coarse grid is used, and the generated droplets are smaller than the grid spacing in the computational domain. Since we are not resolving the small scales, a sufficiently accurate representa-

PT

tion of turbulence in terms of velocity fluctuations and small eddy size in the sub-grid must be provided. Therefore, for closing the RANS equations for

CE

turbulent flow, we have applied the Reynolds stress model (RSM). Because RSM is the most comprehensive type of RANS-based models that can reflect

AC

the velocity fluctuations in each direction every time step by solving additional equations for Reynolds stresses. Therefore, RSM provides more accurate statistics on the flow to be used in the proposed energy-based sub-grid criterion (Saeedipour et al., 2016a). Once this Eulerian-Lagrangian coupling method is established, the governing equations are influenced by activating 15

ACCEPTED MANUSCRIPT

the source terms on the right-hand side. Based on the model mentioned above, the source term for the VOF equation is formulated as Vdrops Vcell ∆t

(10)

CR IP T

Sα = −

where Vdrops is the total volume of droplet(s) generated in the cell and ∆t is

the Eulerian time step. It should be noted that, while this source term is activated, the volume of liquid is replaced by the corresponding volume of gas

AN US

keeping the same velocity in each cell. In other words, this is a conversion in the marker function of VOF without any change in the mixture veloc-

ity. Therefore, the mixture velocity field remains divergence free. But the ~ momentum equation is affected accordingly. Although the velocity field U ~ changes due to extraction and release of does not change, the momentum ρU

M

droplets. Hence, the momentum source term has the following formulation: (init)

+ SD

(11)

ED

~ Smom = (ρl − ρg )Sα U d

where the former term is the momentum of extracted droplets in each cell

PT

containing droplet release, and the latter term SD is the source representing the interaction between the droplets and the ambient gas within the

CE

Eulerian-Lagrangian coupling. The latter is governed by Newton’s second ~ d of each droplet law. Therefore, the position vector ~xd and velocity vector U

AC

are calculated by solving the equations

md

d~xd ~d =U dt

(12)

~d X dU ~ = F dt

(13)

16

ACCEPTED MANUSCRIPT

In the equations above, md denotes the droplet mass calculated from the liquid density and droplet volume at each Lagrangian time step dt. The droplet volume is computed from data on the turbulent flow properties in

CR IP T

~ is typically the representative of all the Eulerian framework. The vector F

the forces acting on a droplet in the spray, such as drag, weight, buoyancy, forces due to pressure gradient and added mass, etc. Some of these forces are

negligibly small due to the high density ratio between droplets and air. In

AN US

the present LPT algorithm, only drag, weight and buoyancy are accounted for. Therefore, equation (13) can be detailed as

md

3 ~d dU πd2 ~ −U ~ d) | U ~ −U ~ d | + πd (ρl − ρ)~g = ρCd (U dt 8 6

(14)

M

For the drag coefficient Cd , a frequently used correlation in spray modelling is

ED

  24    Red  

2/3

24 (1 + 16 Red ) 0.1 < Red < 1000 Red      0.44 Red > 1000

(15)

PT

Cd =

Red < 0.1

CE

where Red = ρ | U − Ud | d/µ is the droplet Reynolds number (Schiller

and Naumann, 1935). Accordingly, one can deduce that the interaction of a

AC

droplet with the surrounding gaseous continuous phase results in a difference in droplet momentum between the instants it enters tin and leaves tout a

specific computational cell Vc , as formulated in the equation   1 X ~ ~ SD = md (Ud )tout − (Ud )tin Vc ∆t d 17

(16)

ACCEPTED MANUSCRIPT

For more details, readers are referred to Saeedipour et al. (2014, 2016a). All these steps act as the main elements of bridging between the statistical nature

CR IP T

of RSM and the deterministic nature of VOF-Lagrangian coupling. The flowchart of the proposed Eulerian-Lagrangian coupling method is presented

in Fig. 5. This multiscale model is developed based on OpenFOAM, and the corresponding modifications are implemented in its C++ libraries.

AN US

4. Jet Atomization Experiments

Existing literature about spray formation by liquid jet atomization is dominated by studies on Diesel injection. This process involves injector hole geometries different from the experiments underlying v. Ohnesorge’s nomogram named in the introduction. Furthermore, in Diesel spray formation

M

cavitation plays a major role. Therefore, for validating the simulation results from the present study, we carried out our own jet atomization exper-

ED

iments with an atomizer design oriented at the one by Haenlein (1931). On the experiments on liquid jet breakup by Haenlein (1931), v. Ohnesorge’s

PT

nomogram was built (v. Ohnesorge, 1936). In those experiments, plain-orifice nozzles with a constant length-to-diameter (or aspect) ratio of 10 were used.

CE

Since the jet breakup mechanism depends on this nozzle property (McCarthy and Molloy, 1974), the breakup regimes in the Ohnesorge nomogram corre-

AC

spond to jets from such nozzles only. The aim of the present experiments was to realize atomizing jets and to characterize velocity and size spectra of the spray drops formed by PDA measurements. With the proper nozzle design we can rely on the Ohnesorge nomogram in selecting the Reynolds number regime ensuring jet atomization for a given Ohnesorge number. 18

ACCEPTED MANUSCRIPT

The plain-orifice nozzle used for the experiments had a hole diameter of 1 mm and the prescribed aspect ratio of 10, as schematically depicted in

CR IP T

Fig. 6. The nozzle body was manufactured from brass and built in the body of an SDX-type Delavan pressure-swirl atomizer (without the SDX inserts).

Liquid was supplied axially to the entrance of the nozzle body by means of a high-pressure displacement pump. According to the Re-Oh nomogram

by v. Ohnesorge (1936), the atomization regime for liquid jets issuing from

AN US

such nozzles into a stagnant atmospheric environment is determined by the inequality Re > 228.5 Oh−0.814 which, for the present experiments with water jets with the Ohnesorge number of 3.7 · 10−3 , requires that the Reynolds number must be greater than 21,700. The water mass flow rate through the

nozzle was set to 127 kg/h, 170 kg/h and 212 kg/h for the experimental cases

M

1-3, leading to Reynolds numbers of 45,000, 60,125 and 75,000, respectively. These Re numbers ensure that the jets break up by atomization. The param-

ED

eters for the three cases are listed in Table 1. The calculated numbers are based on the liquid properties ρl =1000 kg/m3 , µl =1 mPa · s, and σ=0.072

PT

N/m of water.

The measurements of spray properties were carried out using phase-

CE

Doppler anemometry (PDA). This measuring technique acquires the local size and velocity of individual drops by means of a probe volume formed by intersecting pairs of laser beams. The present Dantec PDA system uses two

AC

pairs of beams from an Ar-ion laser for measuring two drop velocity components and the drop size. The probes of the PDA system used have optical apertures of 80 mm. The optical configuration of the PDA system is given by the laser-beam intersection half angle of 1.294◦ and the scattering angle

19

ACCEPTED MANUSCRIPT

of 50◦ . Together with the elevation angles of the three detectors receiving pairs of Doppler signals, the configuration resulted in the two phase factors of

CR IP T

2.528◦ /µm and 1.264◦ /µm and a maximum measurable drop size of 188µm. The dominant light scattering mechanism was refraction.

To obtain information about entire spray cross sections, the probe volume

is traversed to various locations along spray cone diameters in planes at different distances from the nozzle orifice. Assuming axial symmetry of the

AN US

spray, the traversing trajectory goes along one radius between the axis of

symmetry and the edge of the spray. The edge of the spray is defined as the position where drops cross the probe volume with a frequency of no more than 300 Hz, while the drop frequencies are much higher inside the spray. The radial distance between pairs of measurement points was 1 mm.

M

Measurements were carried out at the distances x1 = 200 mm, x2 = 300 mm and x3 = 400 mm downstream from the nozzle orifice, as shown in Fig.

ED

7. Preliminary visualization studies showed that, at these distances from the nozzle orifice, the jet breakup and drop formation was complete and

PT

PDA measurements were possible. At each point of measurement, 20,000 droplets were measured. Primary information about the spray drops consists

CE

in triplets of the two velocity components (u, v) and the size d (further to the times of arrival at and transit through the probe volume). In PDA data, drop size and velocity are correlated, so that they may be represented by

AC

scatter plots as in Fig. 8. The three sub-figures 8(a) to (c) show data at increasing distance from the spray axis in the plane x3 = 400 mm. At the largest distance from the spray axis, the positive correlation of drop size and velocity is most pronounced. This is typical for spray drops injected into a

20

ACCEPTED MANUSCRIPT

stagnant environment. At a given distance from the level of injection, the large drops have the high velocities and vice versa. The correlation is not

CR IP T

very pronounced close to the axis of the spray, since the drop velocity there is largely determined by the momentum flux of the injected liquid.

At each location in the spray cross section, local mean values of the

drop size and velocity may be derived from the PDA data. We look at

number-based mean values here, such as the local number-mean drop size

AN US

D10 , Sauter-mean diameter D32 and the mean x-component U of the drop velocity, which are defined by the equations

PT

ED

D32 =

U =

di φn,i (di )

i=1 I P

i=1 I P

i=1 I P

(17)

φn,i (di )

i=1 I P

M

D10 =

I P

d3i φn,i (di )

(18)

d2i φn,i (di ) ui φn,i (di )

i=1 I P

(19) φn,i (di )

i=1

CE

The local sizes and mean velocities rely on a weighting with the local droplet number fluxes φn,i of drops with size di , which includes the Saffman cor-

AC

rection. The Saffman correction accounts for the effect that the Gaussian intensity distribution of the laser light irradiated into the probe volume produces a larger probe volume for larger droplets seen by the detectors. This is accounted for in the drop number fluxes derived from the PDA data for the various drop sizes (Saffman, 1987). The measured drop sizes are grouped 21

ACCEPTED MANUSCRIPT

into a total of I size classes. Profiles of these mean values in the planes at the smallest and at the largest distance from the nozzle are shown in Fig. 9.

CR IP T

The profiles start at r ≈ 2 mm, since closer to the symmetry axis the sprays

are too dense and validation rates too low to acquire reliable PDA data. The radial extensions of the profiles represent the cross-sectional dimensions of the sprays, which widen with increasing distance from the nozzle. Also, the

mean drop sizes decrease with increasing mass flow rate, indicating the po-

AN US

tential of the atomization process to produce smaller droplets with increasing

jet Reynolds number. Yet, in each spray the mean drop sizes increase with the distance from the nozzle which may be due to coalescence between the drops. Evaporation of the drops, which could make the small droplets disappear earlier than the large ones and also lead to a mean drop size increase,

M

is estimated to be of minor importance in the present experiments. Coalescence, in the contrary, is very probable in sprays of the present kind which

ED

exhibit high spatial drop concentrations and high relative velocities of drops due to the large scatter of the drop velocity (see Fig. 8).

PT

It is of interest for the characterization of spray flows to associate a global mean value of a drop size to each downstream spray cross section. For this

CE

purpose, local PDA data must be converted into global data for each spray cross section by weighting the local data by the local number flux of drops in each size class and by the part of the spray cross sectional area for which

AC

the data are representative. The latter are ring-shaped parts of the cross section with the point of measurement located at the mean radial position of the ring. The equations for the global number-mean and Sauter-mean drop

22

ACCEPTED MANUSCRIPT

sizes, D10,global and D32,global , therefore read di (rj )φn,i,j (rj , di )rj

j=1 i=1 J P I P

(20) φn,i,j (rj , di )rj

j=1 i=1

D32,global =

J P I P

j=1 i=1 J P I P

j=1 i=1

CR IP T

D10,global =

J P I P

d3i (rj )φn,i,j (rj , di )rj

(21)

d2i (rj )φn,i,j (rj , di )rj

AN US

where the areas of the ring-shaped sections were calculated as 2πrj ∆r and for the present equidistant measurement points the constant ∆r, as well as the factor 2π, were cancelled from the fractions. The number flux φn,i,j of drops of each size class i at each radial position j is calculated from the PDA data

M

applying the Saffman correction, since, due to the Gaussian distribution of the laser light intensity in the probe volume, the probe volume has a different

ED

effective cross section for drops of different size. The numbers I and J are the numbers of size classes and measurement positions in the spray cross section, respectively. The resulting global mean drop size data will be shown in Sec.

PT

5.2 together with the results from the simulations.

CE

5. Results and discussion This section presents the results of simulations compared to data from the

AC

experiments. According to the proposed multiscale method, a finite volumebased solver was developed within the framework of OpenFOAM. This solver is used to solve the three-dimensional Navier-Stokes equations coupled with the VOF equation to capture the interface behavior. The PIMPLE algorithm is adopted for pressure-velocity coupling and pressure correction. 23

ACCEPTED MANUSCRIPT

The VOF transport equation is solved by a specific algorithm in OpenFOAM called MULES solver (Multidimensional Universal Limiter with Ex-

CR IP T

plicit Solution) which belongs to algebraic VOF schemes. This method guarantees the robustness of the VOF method when using an artificial compression term to improve the interface resolution (Greenshields, 2015). For turbulence modelling, the Reynolds stress model (RSM) is used to model the

velocity fluctuations in different directions more accurately than other RANS-

AN US

based models. The computational domain and the boundary conditions are

shown in Fig. 10. Assuming D0 as the nozzle diameter, the width, depth and length of the computational domain are Lx = 500D0 and Ly = Lz = 100D0 , respectively.

The computational mesh with structured grid was generated with the

M

minimum grid spacing of ∆xmin = D0 /8 to keep the grid size larger than a typical size of the droplets which is of the order of tens of microns. Since a

ED

transient solution is desired for such highly unsteady flow regimes, the selection of the time step has to be based such that the stability of the numerical

PT

simulation is ensured. Therefore the simulation time steps were determined from the CFL criterion to keep a maximum Courant number of about 0.6. A

CE

turbulent velocity boundary condition is used at the inlet to initiate the injection in the x-direction with a randomly perturbed velocity profile with the intensity of 0.05. Accordingly, the initial stress tensor for the RSM method

AC

is set at the inlet. Outflow and slip wall boundary conditions are applied at the outlet and the surroundings, respectively. All the fluid properties in the simulation were taken at standard room temperature: ρl =1000 kg/m3 , ρg =1.2 kg/m3 , µl =1 mPa · s, µg =0.017 mPa · s and σ=0.072 N/m. 24

ACCEPTED MANUSCRIPT

5.1. Comparison methodology In the characterization of sprays, in particular in view of the intensity of

CR IP T

transport phenomena, the Sauter mean diameter (SMD or D32 ) is recognized as a relevant quantity. The SMD, which is the diameter of a representative droplet whose volume-to-surface ratio is equal to the ratio of the sum of all

droplet volumes in the spray to the sum of all droplet surface areas, is defined as di 3

AN US

n P

i=1 n P

D32 =

(22)

di 2

i=1

We have performed the simulations long enough to assume quasi-steady state conditions for the simulations. The Lagrangian part of the numerical model

M

provides the droplet size distribution at every time step. Using equation (22), the global SMD can be computed at each streamwise distance from the

ED

nozzle exit. This enables us to validate the simulation results against the measurement data.

PT

Further to this, a point-by-point comparison of experiments and simulation data is possible by filtering the Lagrangian discrete data. The local

CE

SMD and local velocity for an arbitrary location of the probe can be com-

AC

puted from:

hD32 (x = (x, y, z))i =

tN P P

tj =t1 p∈P tN P P

tj =t1 p∈P

25

d3p δ(x − xd,p (tj )) d2p δ(x

− xd,p (tj ))

(23)

ACCEPTED MANUSCRIPT

tj =t1 p∈P

|Ud,p |δ(x − xd,p (tj ))

tN P P

tj =t1 p∈P

(24)

δ(x − xd,p (tj ))

CR IP T

hU (x = (x, y, z))i =

tN P P

where P denotes the set containing all droplets in each data ensemble, xd is the current position of droplet and δ is a box filter defined as

 1 | x − xd |< ∆/2

(25)

AN US

δ(x − xd ) =

  0 | x − xd |> ∆/2

where ∆ denotes the grid spacing. The time interval tj+1 − tj between two consecutive probing events should be sufficiently long to avoid double

M

detection of one specific droplet, i.e.

∆ < min Ud,p p∈P (tj+1 − tj )

(26)

ED

This approach allows for a point-by-point comparison with PDA data at given probing positions. The computed data ensembles at each probe volume

PT

location comprise the properties of around 650 individual droplets, which is a much smaller number than the 20,000 droplets measured in the experi-

CE

ment. This difference is due to the computational limitations and must be accounted for when comparing statistical spray properties between compu-

AC

tation and experiment.

5.2. Droplet size distributions At each downstream distance from the nozzle exit, the global SMD can be computed by equations (22) and (21) from simulations and experiments, 26

ACCEPTED MANUSCRIPT

respectively. Table 2 represents the SMD values from simulation and measurement in the three different downstream spray cross sections introduced

CR IP T

in Sec. 4. Figure 11 shows these results as a function of the downstream location x.

A point-by-point comparison of local SMD values was performed based on the aforementioned comparison methodology. The SMD of the spray

droplets can be obtained at different radial distances from the jet axis in

AN US

each downstream cross section. For each case, the local SMDs of droplets are plotted at the distances x1 = 200 mm and x2 = 300 mm downstream

from the nozzle orifice (Fig. 12). In view of the quite different numbers of data from simulations and measurements used for the statistics, the results reveal reasonable agreement. The SMD values decrease with increasing mass

M

flow rate (from case 1 to 3) because of the increase in jet Reynolds number. This effect is correctly represented by the simulations.

ED

As also shown by the experiments (Fig. 9), the SMD increases with the distance from the nozzle exit. This may be due to the coalescence between

PT

the droplets. In the present computational approach no droplet coalescence model is included, but the sub-grid breakup model determines the droplet

CE

size based on the turbulence near the jet surface in which the turbulent length scales increase in the streamwise direction (Sallam et al., 2002). Accordingly, the global and local SMD values increase with the distance from the nozzle

AC

as demonstrated in Fig. 12. 5.3. Drop velocity profiles in the sprays Similar to the mean droplet sizes, the profiles of mean droplet velocities at each streamwise location can be obtained using the local droplet ensembles. 27

ACCEPTED MANUSCRIPT

For each case, the radial profile of the droplet velocity x components are plotted at distances x1 = 200 mm and x2 = 300 mm as shown in Fig. 13.

CR IP T

The velocity profiles obtained by simulations reveal a trend of variation similar to that obtained from the measurement data. The droplets closer

to the spray axis have higher velocities, while they slow down at the larger radial distances from the jet axis, as also found in the experiments. But

there is a slight deviation in the inclination of the velocity profiles compared

AN US

to the experiments. In order to discuss the deviation in each case, it is worth

mentioning that, in high-speed sprays, the drag force acts as the dominant force on the droplets and is the main mechanism of droplet deceleration. The deceleration rate is inversely proportional to a power of the droplet size, starting with dUd /dt ∝ d−2 in the Stokes regime, to dUd /dt ∝ d−1

M

for the regime with constant drag coefficient, in accordance with equations (14) and (15). Therefore, the larger droplets face lower deceleration rate in

ED

their Lagrangian motion. This explains the positive correlation of drop size and velocity found in the experiments (Fig. 8). Further to the size-velocity

PT

correlations, this physical analysis is in agreement with the trend of the velocity profiles from the simulations. It is illustrated in Fig. 12 that any

CE

over- or underestimation of droplet sizes exhibits a similar corresponding behavior in the velocity profiles. In other words, when the local SMD is overestimated by the model, its local velocity is also overestimated and vice

AC

versa.

Moreover, the spatial concentration of droplets moving in the spray may

affect the droplet trajectories and velocities. As there exists no established model for the drag coefficient of the droplets moving in a cloud, we applied the

28

ACCEPTED MANUSCRIPT

drag coefficient for an individual droplet in the Lagrangian formulation (i.e. equation (15)). But the cloud effect might lead to a different drag coefficient

CR IP T

and deceleration rate of droplets in the real spray. However, accounting for this effect is beyond the scope of the current study. Additionally, the present model does not include the interaction of dispersed droplets with gas turbulence (i.e. turbulent dispersion effect) as well as the effect of droplet impact to the liquid core which may affect the droplet motion.

AN US

It should be emphasized that this model determines the spray proper-

ties during primary breakup based on steady-statistical representation of the turbulent flow, without resolving the interfacial eddies. Sometimes, this may result in an overall over-prediction of the droplet size in the modelling of relatively fine sprays (e.g. case 3). Because coarse grid may be inadequate to

M

estimate an accurate characteristic length scale in highly-frequent sub-grid events. This compromise between accuracy and efficiency, however, leads to

ED

a reasonable performance in simulating high-speed sprays with coarse grid and low computational costs. In the case of having a more accurate turbu-

PT

lence modelling, such as large eddy simulation, for instance in case 3, we would need a grid resolution of at least 10 − 20 µm to resolve the droplets.

CE

Thus, an LES would be approximately up to four orders of magnitude more expensive.

AC

5.4. Size-velocity correlation In the PDA data, droplet size and velocity are correlated due to the mea-

surement principle. Data are shown in Fig. 8. It is also possible to present this correlation for the simulation results at given probe volume locations. Figures 14(a), (c) and (e) show the size-velocity correlation measured with 29

ACCEPTED MANUSCRIPT

PDA at different radial distances from the spray axis in the plane at x2 for cases 1, 2 and 3, respectively. The corresponding size-velocity correlations of

CR IP T

filtered Lagrangian data from the simulation at the same radial positions for each case are plotted in Figs. 14(b), (d) and (f). Although the simulation

data consist of much smaller numbers of drops than those of the experiment,

both scatter plots reveal the same correlations between droplet sizes and velocities. For comparison, the two data sets are plotted in diagrams with

AN US

equal axes.

The ranges of drop size and velocity in the simulations are well within those from the experiments, but covering only parts of the experimental ranges. It can be deduced that more simulation data would follow similar trends. As discussed in Sec. 3, resolving all the interfacial turbulence length

M

scales by the coarse grid is not practically feasible. Instead, they are modelled with a characteristic length scale in the inertial sub-range of the turbulent

ED

kinetic energy spectrum (i.e. the Taylor microscale λ). Accordingly, the droplets generated in the simulation are the smallest droplets that can be

PT

formed by the energy-based sub-grid criterion. In the real spray, however, various primary breakup mechanisms are involved in the generation of the

CE

droplets, which all appear in the experimental scatter plots. The ratio between Taylor microscale and the smallest scales in the jet

(comparable to the Kolmogorov length scale) is proportional to Re1/4 (Pope,

AC

2000). It is evident that, for all the three jets, the smallest droplets in the simulation are always larger than the smallest ones in the corresponding real spray. Nevertheless, this limitation would not affect too much the statistics for computation of local D32 , since the larger droplets found in the experi-

30

ACCEPTED MANUSCRIPT

ments are captured by the simulation.

CR IP T

6. Conclusions In this paper, the primary atomization of water jets at high Reynolds

numbers was numerically investigated using a multiscale Eulerian-Lagrangian approach and measured in a series of experiments. For the simulations, un-

resolved droplet formation was modeled by a sub-grid breakup model on the

AN US

basis of a physical description of the primary atomization. A numerical hybrid framework was established to predict the general behavior of the liquid jet by employing the VOF method together with LPT to track the droplets

in small scale. In the jet atomization experiments, water sprays were produced at high jet Reynolds numbers O(104 ) using a plain-orifice nozzle of the

M

Haenlein (1931) geometry, and the drop velocity and size were measured using phase-Doppler anemometry (PDA). These experiments were carried out

ED

for validation of the simulation results and evaluation of the performance of the sub-grid drop formation model.

PT

This sub-grid model has some limitations in case of modelling small-scale interfacial phenomena, but it proposes a low-cost method for capturing the

CE

dynamics of droplets due to the breakup. A filtering operation on the Lagrangian discrete data is carried out to enable point-by-point comparisons

AC

between simulation and experimental data. The Sauter-mean drop sizes obtained from the simulations and the experiments agree reasonably well. In addition, the velocity distribution of droplets were investigated both numerically and experimentally. The mean drop velocity profiles obtained from simulations agree reasonably with the experiments, even though this mul31

ACCEPTED MANUSCRIPT

tiscale model does not resolve all the interfacial scales. Furthermore, the positive correlation of drop size and velocity is clearly demonstrated by both

CR IP T

experiment and simulation. However, there are various shortcomings to be overcome in future work, such as (i) considering a distribution of length scales in the sub-grid model at the instant of droplet formation, (ii) accounting for

turbulent dispersion of the Lagrangian droplets and (iii) considering the effect of spatial concentration of droplets on the drag coefficient during the

AN US

integration of the Lagrangian trajectories. Future studies will also focus on employing this multiscale modeling approach for different nozzle configurations and injection conditions for optimization purposes. 7. Acknowledgments

M

We acknowledge funding by Christian-Doppler Research Association, the Austrian Federal Ministry of Economy, Family and Youth, and the Austrian

ED

National Foundation for Research, Technology and Development. The authors also want to acknowledge the financial support of the K1MET center

PT

for metallurgical research in Austria, which is partly funded by the Austrian

CE

government (www.ffg.at). References

AC

Brackbill, J., Kothe, D., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of Computational Physics 100, 335–354.

Chesnel, J., Reveillon, J., Menard, T., Demoulin, F.-X., 2011. Large eddy simulation of liquid jet atomization. Atomization and Sprays 21, 711–736.

32

ACCEPTED MANUSCRIPT

Desjardins, O., McCaslin, J., Owkes, M., Brady, P., 2013. Direct numerical and large-eddy simulation of primary atomization in complex geometries.

CR IP T

Atomization and Sprays 23, 1001–1048. Desjardins, O., Moureau, V., Pitsch, H., 2008. An accurate conservative level set/ghost fluid method for simulating turbulent atomization. Journal of Computational Physics 227, 8395–8416.

AN US

Desjardins, O., Pitsch, H., 2010. Detailed numerical investigation of turbulent atomization of liquid jets. Atomization and Sprays 20, 311–336.

Dumouchel, C., 2008. On the experimental investigation on primary atomization of liquid streams. Experiments in Fluids 45, 371–422.

M

Faeth, G., Hsiang, L., Wu, P., 1995. Structure and breakup properties of sprays. International Journal of Multiphase Flow 21, 99–127.

ED

Fuster, D., Bagu´e, A., Boeck, T., Le Moyne, L., Leboissetier, A., Popinet, S., Ray, P., Scardovelli, R., Zaleski, S., 2009. Simulation of primary at-

PT

omization with an octree adaptive mesh refinement and VOF method. International Journal of Multiphase Flow 35, 550–565.

CE

Gordillo, J., P´erez-Saborid, M., 2005. Aerodynamic effects in the break-up of liquid jets: on the first wind-induced break-up regime. Journal of Fluid

AC

Mechanics 541, 1–20.

Gorokhovski, M., Herrmann, M., 2008. Modeling primary atomization. Annu. Rev. Fluid Mech. 40, 343–366.

33

ACCEPTED MANUSCRIPT

Greenshields, C. J., 2015. OpenFOAM The Open Source CFD Toolbox. Tech. rep.

CR IP T

URL http://www.openfoam.com Grosshans, H., Sz´asz, R., Fuchs, L., 2014. Development of an efficient sta-

tistical volumes of fluid Lagrangian particle tracking coupling method. International Journal for Numerical Methods in Fluids 74, 898–918.

AN US

¨ Haenlein, A., 1931. Uber den Zerfall eines Fl¨ ussigkeitsstrahles. Forschung auf dem Gebiete des Ingenieurwesens 2, 139–149.

Ham, F., Young, Y.-N., Apte, S., Herrmann, M., 2003. A hybrid EulerianLagrangian method for LES of atomizing spray. Advances in Fluid Mechan-

M

ics, WIT Press (Computational Methods in Multiphase Flow II), 313–322. Herrmann, M., 2005. A Eulerian level set/vortex sheet method for two-phase

ED

interface dynamics. Journal of Computational Physics 203, 539–571. Herrmann, M., 2011. On simulating primary atomization using the refined

PT

level set grid method. Atomization and Sprays 21, 283–301. Herrmann, M., 2013. A sub-grid surface dynamics model for sub-filter surface

CE

tension induced interface dynamics. Computers and Fluids 87, 92–101.

AC

Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39, 201–225.

Lafrance, P., 1975. Nonlinear breakup of a laminar liquid jet. Physics of Fluids 18, 428–432.

34

ACCEPTED MANUSCRIPT

Lebas, R., Menard, T., Beau, P. A., Berlemont, A., Demoulin, F. X., 2009. Numerical simulation of primary break-up and atomization: DNS and

CR IP T

modelling study. International Journal of Multiphase Flow 35, 247–260. Ling, Y., Zaleski, S., Scardovelli, R., 2015. Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model. International Journal of Multiphase Flow 76, 122–143.

AN US

McCarthy, M., Molloy, N., 1974. Review of stability of liquid jets and the influence of nozzle design. The Chemical Engineering Journal 7, 1–20.

M´enard, T., Beau, P. A., Tanguy, S., Demoulin, F. X., Berlemont, A., 2006. Primary break up modeling Part A : DNS, a tool to explore primary break up. In: 10th Triennial International Conference on Liquid Atomization and

M

Spray Systems (ICLASS 2006), Kyoto, Japan.

ED

M´enard, T., Tanguy, S., Berlemont, A., 2007. Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary

524.

PT

break-up of a liquid jet. International Journal of Multiphase Flow 33, 510–

CE

Osher, S., Sethian, J. A., 1988. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of

AC

Computational Physics 79, 12–49.

Pavlovic, Z., Scheidl, S., Edelbauer, W., Basara, B., Brenn, G., Jakirlic, S., 2014. Numerical investigation of the liquid core length in sprays with fully turbulent boundary condition. In: 26th European Conference on Liquid Atomization and Spray Systems (ILASS 2014), Bremen, Germany. 35

ACCEPTED MANUSCRIPT

Pope, S., 2000. Turbulent Flows. Cambridge University Press. Saeedipour, M., Pirker, S., Bozorgi, S., Schneiderbauer, S., 2016a. An

CR IP T

Eulerian-Lagrangian hybrid model for the coarse-grid simulation of turbulent liquid jet breakup. International Journal of Multiphase Flow 82, 17–26.

Saeedipour, M., Pirker, S., Schneiderbauer, S., 2016b. Multiscale simulation of liquid jet disintegration and primary atomization using Eulerian-

AN US

Lagrangian coupling. In: International Conference on Multiphase Flow, Firenze, Italy.

Saeedipour, M., Schneiderbauer, S., Pirker, S., Bozorgi, S., 2014. Numerical simulation of turbulent liquid jet breakup using a sub-grid criterion with

M

industrial application. In: 26th European Conference on Liquid Atomization and Spray Systems (ILASS 2014), Bremen, Germany.

ED

Saffman, M., 1987. Automatic calibration of lda measurement volume size.

PT

Applied Optics 26, 2592–2597. Sallam, K., Dai, Z., Faeth, G. M., 2002. Liquid breakup at the surface of tur-

CE

bulent round liquid jets in still gases. International Journal of Multiphase Flow 28, 427–449.

AC

Schiller, L., Naumann, A., 1935. A drag coefficient correlation. Zeitschrift des Vereins Deutscher Ingenieure 77, 318–320.

Shinjo, J., Umemura, A., 2010. Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation. International Journal of Multiphase Flow 36, 513–532. 36

ACCEPTED MANUSCRIPT

Srinivasan, V., Salazar, A. J., Saito, K., 2008. Numerical investigation on the disintegration of round turbulent liquid jets using LES/VOF techniques.

CR IP T

Atomization and Sprays 18, 571–617. Sussman, M., Puckett, E. G., 2000. A coupled level set and volume-of-fluid

method for computing 3D and axisymmetric incompressible two-phase flows. Journal of Computational Physics 162, 301–337.

AN US

Tomar, G., Fuster, D., Zaleski, S., Popinet, S., 2010. Multiscale simulations of primary atomization. Computers & Fluids 39, 1864–1874. Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct Numerical Simulations of Gas-Liquid Multiphase Flows. Cambridge University Press.

M

v. Ohnesorge, W., 1936. Die Bildung von Tropfen an D¨ usen und die Aufl¨osung fl¨ ussiger Strahlen. Zeitschrift f¨ ur Angewandte Mathematik und Mechanik

ED

16, 355–358.

Wu, K.-J., Reitz, R., Bracco, F., 1986. Measurements of drop size at the

PT

spray edge near the nozzle in atomizing liquid jets. Physics of Fluids 29, 941–951.

CE

Wu, P., Faeth, G. M., 1993. Aerodynamic effects on primary breakup of turbulent liquids. In: 31st Aerospace Sciences Meeting & Exhibition (AlAA

AC

93-0903), Reno, NV.

Wu, P., Tseng, L., Faeth, G. M., Arbor, A., Reno, J., 1992. Primary Breakup in Gas/Liquid Mixing Layers for Turbulent Liquids. In: 30th Aerospace Sciences Meeting & Exhibiti (AlAA 92-0462), Reno, NV. 37

ACCEPTED MANUSCRIPT

Xiao, F., 2012. Large eddy simulation of liquid jet primary breakup. Ph.D. thesis, Loughborough University.

Atomization and Sprays 14, 53–80.

CR IP T

Yi, Y., Reitz, R., 2004. Modeling the primary breakup of high-speed jets.

Yuen, M.-C., 1968. Non-linear capillary instability of a liquid jet. Journal of

AC

CE

PT

ED

M

AN US

Fluid Mechanics 33, 151–163.

38

ACCEPTED MANUSCRIPT

Table 1: The parameters of the three jet atomization experiments. The nozzle hole diameter is 1 mm, the jet Ohnesorge number 3.7 · 10−3 .

Re [-] 45, 000 60, 125 75, 000

CR IP T

Jet velocity [m/s] 45 60 75

AC

CE

PT

ED

M

AN US

Case 1 Case 2 Case 3

Mass flow rate m ˙ [kg/h] 127 170 212

39

AN US

CR IP T

ACCEPTED MANUSCRIPT

Table 2: The global SMD from simulations and experiments in µm at each streamwise location. Deviations in percent relate to the experiment.

ED

M

x2 Exp. Sim. ∆ [%] 87 91 5 62 64 3 38 52 37

AC

CE

PT

Case 1 Case 2 Case 3

x1 Exp. Sim. ∆ [%] 82 89 9 57 61 7 40 46 15

40

x3 Exp. Sim. ∆ [%] 89 94 6 65 68 5 42 55 31

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

ED

(a)

AC

CE

PT

Figure 1: Schematic of the multiscale approach: (a) The pure Eulerian simulation of twophase flow on fine grid. (b) The Eulerian simulation of the continuous liquid phase (in blue) on coarse grid and the Lagrangian representative particles (in gray) on the sub-grid level.

41

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

Figure 2: Schematic of energy-based mechanism for breakup near the interface (redrawn from Wu and Faeth (1993)). The liquid-gas interface is presented by solid lines.

42

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 3: Description of turbulence energy cascade theory (redrawn from Pope (2000)).

43

AN US

CR IP T

ACCEPTED MANUSCRIPT

gas

M

~d U

liquid

ED

λ

(b) t = ti + ∆t

PT

(a) t = ti

d'λ

AC

CE

Figure 4: Schematic of the sub-grid breakup model for Eulerian-Lagrangian coupling.

44

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 5: Flow chart of the algorithm for multiscale modelling of the jet. The coupling between Reynolds stress model and volume-of-fluid (RSM-VOF) and Lagrangian Particle Tracking (LPT) is established according to the sub-grid criterion.

45

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 6: A schematic cross section through the plain-orifice nozzle.

46

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

ED

P DA

x1 x2 x3

AC

CE

Figure 7: Schematic of the planes of measurement in the sprays.

47

ACCEPTED MANUSCRIPT

100

|

U

d

| [m/s]

60

40

20

0

0

50

100

150

d [ m]

CR IP T

80

200

AN US

(a)

100

80

40

20

0

0

M

|

U

d

| [m/s]

60

50

100

150

200

d [ m]

ED

(b)

100

PT

80

AC

d

U |

CE

| [m/s]

60

40

20

0

0

50

100

150

200

d [ m]

(c)

Figure 8: Size-velocity correlations of drops measured in sprays at the lowest mass flow rate m ˙ = 127kg/h at the largest distance x3 = 400 mm from the nozzle orifice, at the radial distances of (a) 2 mm, (b) 6 mm and (c) 12 mm from the spray axis.

48

ACCEPTED MANUSCRIPT

25

[ m]

15 Case 1, X

1

D

10

Case 1, X

3

10

Case 2, X

1

Case 2, X

3

Case 3, X

1

5

0

Case 3, X

3

0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 r [mm]

90 80 70

[ m]

60

32

(a)

AN US

100

D

CR IP T

20

50 40

Case 1, X

1

30

Case 1, X

3

Case 2, X

20

0

0

2

4

6

M

1

10

Case 2, X

3

Case 3, X

1

Case 3, X

3

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38

ED

r [mm]

(b)

70 65 60

Case 1, X

55

Case 1, X

50

Case 2, X

PT

1

3

1

Case 2, X

AC

U [m/s]

CE

45

3

Case 3, X

40

1

Case 3, X

3

35 30 25 20 15 10 5 0

0

2

4

6

8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 r [mm]

(c)

Figure 9: Radial profiles of (a) the number-mean drop size D10 , (b) the Sauter-mean drop size D32 and (c) the number-mean drop velocity X-component U measured at the three different liquid mass flow rates in the closest and the farthest downstream cross sections.

49

CR IP T

ACCEPTED MANUSCRIPT

inlet

slip

Lx

ED

M

slip

AN US

D0

PT

Lz

outlet Ly

AC

CE

Figure 10: The computational domain with its boundary conditions.

50

CR IP T

ACCEPTED MANUSCRIPT

110 Case 2 (exp) Case 3 (exp)

80

Case 1 (sim)

70

Case 2 (sim)

60

Case 3 (sim)

50 40

M

D32,global (µm)

90

Case 1 (exp)

AN US

100

30

10

50

100

150

200

250

x (mm)

PT

0 0

ED

20

300

350

400

450

500

AC

CE

Figure 11: Global Sauter-mean drop size D32,global of the droplets for each spray from simulation and PDA measurement. The open and filled symbols denote values from the numerical simulation and the measurement, respectively.

51

ACCEPTED MANUSCRIPT

120 100

60 40

X (exp) 1

X (sim) 1

20

X (exp) 2

X2 (sim) 0 0

2

4

70 60

D32 [µm]

6

8

10

(a)

AN US

80

r [mm]

CR IP T

D32 [µm]

80

50

X1 (exp)

40

M

X1 (sim)

30

2

4

6

X2 (sim) 8

r [mm]

10

12

14

(b)

ED

20 0

X2 (exp)

80

X1 (exp)

70

X1 (sim)

PT

X (exp)

AC

CE

D32 [µm]

60

2

X2 (sim)

50 40 30

20 0

2

4

6

8

r [mm]

10

12

14

(c)

Figure 12: Radial profiles of the local Sauter-mean drop size D32 at different downstream cross sections for (a) case 1, (b) case 2 and (c) case 3. The open and filled symbols denote values from the numerical simulation and the experiment, respectively.

52

ACCEPTED MANUSCRIPT

45

X1 (exp)

40

X1 (sim) X2 (exp)

35

U [m/s]

CR IP T

X2 (sim)

30 25 20 15 10 5 0 0

r [mm]

10

15

(a)

AN US

5

60

X1 (exp) X1 (sim)

50

X (exp) 2

X (sim)

U [m/s]

40

2

30

10

5

ED

0 0

M

20

r [mm]

10

(b)

100

X (exp) 1

90

X (sim) 1

80

PT

X2 (exp) X2 (sim)

U [m/s]

CE

70

AC

15

60 50 40 30 20 10

0 0

5

10

r [mm]

15

20

(c)

Figure 13: Radial profiles of droplet velocity x component at different downstream cross sections for (a) case 1, (b) case 2 and (c) case 3. The open and filled symbols are for values from the numerical simulation and the experiment, respectively.

53

ACCEPTED MANUSCRIPT

80

80 70 60

U |

50 40 30

20

20 10 0

0

50

100

150

200

0

d [ m]

0

CR IP T

|Ud| [m/s]

40

d

| [m/s]

60

50

100

150

200

d [µm]

(b) case 1, r = 6 mm (sim)

AN US

(a) case 1, r = 6 mm (exp) 80

80 70

60

|Ud| [m/s]

50

40

40

|

U

d

| [m/s]

60

30

20

0

0

50

100 d [ m]

M

20

150

10

200

0

0

ED

70 60

|Ud| [m/s]

U |

CE AC

0

50

200

80

20

0

150

(d) case 2, r = 10 mm (sim)

PT

40

d

| [m/s]

60

100

d [µm]

(c) case 2, r = 10 mm (exp) 80

50

50 40 30 20 10

100

150

200

0

d [ m]

0

50

100

150

200

d [µm]

(e) case 3, r = 10 mm (exp)

(f) case 3, r = 10 mm (sim)

Figure 14: Size-velocity correlations of droplets in the experiments (left) and simulations (right) for each case at the distance x2 = 300 mm from the nozzle orifice and different radial distances from the spray axis.

54