Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model

International Journal of Multiphase Flow 76 (2015) 122–143

Contents lists available at ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model Y. Ling a,b,⇑, S. Zaleski a,b, R. Scardovelli c a

Institut Jean Le Rond d’Alembert, Sorbonne Universités, UPMC Univ Paris 06, UMR 7190, F-75005 Paris, France Institut Jean Le Rond d’Alembert, CNRS, UMR 7190, F-75005 Paris, France c DIN – Lab. di Montecuccolino, Università di Bologna, Via dei Colli 16, 40136 Bologna, Italy b

a r t i c l e

i n f o

Article history: Received 13 January 2015 Received in revised form 1 June 2015 Accepted 9 July 2015 Available online 17 July 2015 Keywords: Multiphase flows Atomization Lagrangian point-particle Volume of fluid method Multiscale modeling

a b s t r a c t Modeling and simulation of atomization is challenging due to the existence of a wide range of length scales. This multiscale nature of atomization introduces a fundamental challenge to numerical simulation. A pathway to comprehensive modeling is still to be found. The present study proposes a multiscale multiphase flow model for atomization simulations, where the large-scale interfaces are resolved by the Volume-of-Fluid (VOF) method and the small droplets by the Lagrangian point-particle (LPP) model. Particular attention is focused on the momentum coupling between LPP and resolved flow and the conversion between droplets represented by VOF and LPP. A series of multiphase flow problems are considered to validate the model. The results obtained by a number of simulations are compared against direct numerical simulation (DNS) results and experimental data. In particular, the model is applied to simulate the gas-assisted atomization experiment, and the numerical results are compared to the experimental measurements for a quantitative validation. Ó 2015 Elsevier Ltd. All rights reserved.

Introduction Atomization is an important fundamental multiphase flow problem. As it is essential to many engineering applications, such as fuel injection in combustion engines, many experimental and numerical studies on the topic have been reported in the literature (Lasheras et al., 1998; Marmottant and Villermaux, 2004; Descamps et al., 2008; Fuster et al., 2009; Shinjo and Umemura, 2010; Tomar et al., 2010; Matas et al., 2011; Fuster et al., 2013). Nevertheless, the mechanisms that control spray formation and evolution are far from being understood. Numerical simulation has been shown to be a powerful tool to investigate atomization, in particular in providing physical insight into fundamental physics and aspects of the phenomenon that are difficult to be examined experimentally. However, modeling and simulation of atomization is challenging due to the existence of a wide range of length scales, as in an atomizing liquid jet they can vary from the length of the jet (several cm) to the size of the smallest droplets (submicrons). This multiscale nature of atomization introduces a fundamental challenge to numerical simulation. A pathway to comprehensive modeling is still to be found. ⇑ Corresponding author at: Institut Jean Le Rond d’Alembert, Université Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris, France. E-mail address: [email protected] (Y. Ling). http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.07.002 0301-9322/Ó 2015 Elsevier Ltd. All rights reserved.

In the past decades, sophisticated numerical schemes have been developed to capture or track the interface evolution, such as the Volume-of-Fluid (VOF) (Hirt and Nichols, 1981), Level-Set (LS) (Sussman et al., 1994), and Front-Tracking methods (Tryggvason et al., 2001). With these numerical schemes accurate direct numerical simulations (DNS) of multiphase flows with complex interfaces can be conducted (Fuster et al., 2009; Herrmann, 2010; Shinjo and Umemura, 2010; Fuster et al., 2013). In particular, the VOF method has the advantage of conserving mass and is in a good balance between accuracy and implementation complexity. In previous studies of our group, VOF has been applied to a variety of multiphase flows with complex interface topology change, like primary atomization, droplet splash, rising bubble, yielding accurate simulation results (Josserand et al., 2005; Fuster et al., 2009; Popinet, 2009; Fuster et al., 2013). Direct numerical simulation with interface capturing/tracking schemes is a very powerful tool to understand the fundamental physics of atomization. However, the range of scales that can be covered by DNS with today’s computing capability is still smaller than what would be needed for practical applications or even comparison to experiments. Adaptive Mesh Refinement (AMR) techniques have been developed and utilized in multiphase flow simulation to reduce the cost of computation (Popinet, 2009; Agbaglah et al., 2011). In AMR, a fine mesh is used only in the regions that require a higher resolution, such as near the interface

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

and in shear layers. Nevertheless, even with AMR, DNS is still often too expensive for many atomization problems. First of all, it is well known that octtree AMR adds a large overhead to computations. This is mainly due to the placement of objects in memory, resulting numerous page faults. Because 3D block-AMR requires the creation of a very large number of boxes, similar difficulties may also apply in that case. Second, the parallelization of AMR adds additional challenges. Regions in which a high density of small computations cells exist are not always distributed uniformly in space (in atomization they are close to the nozzle) which makes domain decomposition into blocks of equal size difficult. Third, the most demanding mesh requirement is usually introduced by the small droplets that are formed in atomization. Typically, to obtain reasonable results of motion and interface dynamics for a spherical droplet, at least eight to twelve grid points across the droplet diameter are required. It is known from experiments that the diameter of the atomizing liquid jet is about three orders of magnitude larger than the smallest droplets. For such a wide range of length scales, the number of grid points that is required to resolve every scale, namely a true DNS, is extremely expensive even with AMR. An alternative approach to resolve this multiscale challenge in atomization simulation is to introduce subgrid models, although the two ideas could be combined. The idea of subgrid modeling is very simple: only the physical scales larger than the cell size are resolved, while for the physical scales that are smaller suitable models are introduced to represent the relevant physics. A good example of subgrid modeling is the Large-Eddy Simulation (LES) model for turbulence (Pope, 2000). In LES, a low-pass filtering is applied and only the large eddies are resolved. As a result, a coarser mesh compared to DNS of turbulence can be used. However, to accurately capture the large-scale physics, the small unresolved turbulent eddies and their interaction with the large scales have to be represented by a subgrid model, such as the Smagorinsky model (Smagorinsky, 1963). For turbulent multiphase flows like atomization, subgrid modeling is more challenging since both the unresolved turbulence and multiphase scales, and the interaction between two unresolved scales need to be modeled properly (Lakehal and Liovic, 2011; Lakehal et al., 2012; Chesnel et al., 2011a,b; Vincent et al., 2010; Larocque et al., 2010; Toutant et al., 2009a,b; Labourasse et al., 2007; Sirignano, 2005). If a spatial filter of a constant filtering length is applied to the whole flow field, then both the flow properties and the phase-characteristic function need to be filtered (Lakehal and Liovic, 2011; Sirignano, 2005). As a result, an ideal subgrid model would need to capture a variety of multiphase dynamics, such as the formation of sheets, ligaments, and droplets, droplet dynamics, droplet impact on liquid interface, and the effect of both resolved and unresolved turbulence on these pheomena. Although progress has been made in the past decade (Lakehal et al., 2012), such comprehensive subgrid models are exceedingly complex and the current knowledge of many of the above phenomena still requires significant improvement for the development of accurate models. Nevertheless, as the droplets formed in atomization are usually very small, it is easier to only apply subgrid modeling for their motion, and the Lagrangian point-particle (LPP) approach is a perfect candidate. The LPP model has been widely utilized in particle-laden flows and sprays (Balachandar and Eaton, 2010; Apte et al., 2003). In LPP, the particles (or droplets) are viewed as point masses and the flow around individual particles is not resolved. To accurately track the particles and the backward effect to the carrier flow, physical models are needed to calculate the instantaneous force acting on the particles. Eventually, in a so-called Eulerian–Lagrangian simulation, the macro-scale fluid flow is resolved in an Eulerian framework and the point-particles

123

are evolved in a Lagrangian framework. In atomization, the typical Weber number of the droplets, Wep , is usually very small. (Precise definition of Wep will be given in Section ‘Modeling and numerical methods’.) As a result, the droplets remain close to a spherical shape, and the LPP model is suitable to capture their motion and the backward effect to the flow. From the filtering perspective, LPP can be viewed as a local spatial filter with a variable filtering length scale (Climent and Magnaudet, 2006). Since it is difficult to filter the overall flow field, filtering is only applied at the particle location with a filtering length scale which is related to the droplet size. In this work, a multiscale simulation approach that couples the LPP model and an interface-capturing scheme is proposed. In this approach, the goal is to filter the scales on the order of a small number, say nc , of grid cells, and to represent larger scales. Thus droplets of a diameter smaller than or of the order of nc Dx are represented by the LPP model, while larger droplets and fluid masses are represented by an interface tracking scheme. Thus part of the flow is represented in a DNS-like manner while both the flow inside the small droplets and the surrounding flow perturbed by these droplets are filtered. Thus the motion of the small droplets and their backward effect to the resolved flow are represented by the LPP force model. Moreover, care is taken to work with Reynolds and Weber numbers sufficiently moderate to minimize the energy present at such small scales. In recent years, several pioneering efforts have been made to couple the LPP model with interface-capturing schemes, such as LS (Herrmann, 2010) and VOF (Tomar et al., 2010). The results are quite encouraging and show the great potential of this multiscale approach for accurate large-scale simulations of atomization. However, there still exist unresolved issues in the coupling of LPP with interface-capturing schemes which require further investigation and validation. In conventional simulations with the LPP model the grid cell is much larger than the particle size, typically dp < Dx=10, where dp and Dx represent the particle diameter and cell size, respectively, as shown in Fig. 1(a). As a result, a single computational cell may contain many particles. Sometimes, the number of particles in a cell is so large that a super-particle technique is utilized, in which one computational particle (or parcel) is used to represent many physical particles. In contrast, for a fully resolved droplet (RD), typically about ten grid points per diameter are required, as in Fig. 1(d). The ratio dp =Dx has a significant jump of about two orders of magnitude between these two limits. Even if AMR is considered, it is impossible to adapt the mesh to such a spatial jump in a single time step when a resolved droplet is converted to a Lagrangian particle or vice versa. In yet another numerical setup, if the grid were stretched and the droplet moved from a high density grid to a coarse grid region, a transition region of poor resolution would again be encountered. Thus in all cases the mesh has to be adapted gradually, covering regimes where the grid resolution of the droplet lies in between the two extreme limits, as in Fig. 1(b) and (c). For these two intermediate conditions, dp is only a few times larger than Dx, and with such a resolution the flow inside and outside the droplet is not captured accurately by DNS. It is therefore reasonable to convert these poorly-resolved droplets to LPP, but because of the larger-than-grid-size droplets the required modeling is much more challenging than in conventional case. The challenge mainly comes from the coupling between the LPP droplets and the resolved fluid flow. The fluid-dynamic force acting on each individual droplet is crucial to LPP models, since it determines both the droplet motion and the backward effect to the fluid. Typically, the calculation of this force requires the undisturbed fluid velocity at the droplet location. For the simple case of a single droplet, the ‘‘undisturbed fluid velocity’’ can be considered as the fluid velocity without the local perturbation induced by the

124

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Fig. 1. Droplet representation at different grid resolutions: (a) dp  Dx=10, (b) dp  Dx, (c) dp  4 Dx, and (d) dp  10 Dx.

droplet, or the fluid velocity far away from the droplet location if a uniform flow is imposed. For the case with many droplets, the undisturbed fluid velocity ‘‘seen’’ by an LPP droplet can be considered as the fluid velocity when the droplet of interest is absent but all the other droplets are present (Balachandar, 2009). Because the droplets and the fluid are two-way coupled, not only the fluid exerts a fluid-dynamic force on the droplets, but this force is also exerted back to the fluid with an opposite sign. Conventionally, the droplet force is applied to the fluid as a point source located at the center of mass of the droplet (Climent and Magnaudet, 2006). Strictly speaking, after the droplet force is added, the undisturbed fluid velocity is not available anymore, since the fluid flow has been ‘‘perturbed’’. However, for a droplet that is much smaller than the computational cell the perturbation induced by each individual droplet to the fluid within a cell is indeed small. As a result, the fluid velocity stored in the cell is a good approximation of the undisturbed fluid velocity. On the other hand, for droplets that are larger than the cell size, if the force is also applied as a point source then the fluid velocity of the cell, where the center of mass of the droplet is located, will be influenced significantly and it cannot be considered equal to the undisturbed fluid velocity. Furthermore, if the undisturbed fluid velocity is not computed correctly, an error will be present in the droplet force and it will propagate back to the resolved fluid flow, thus leading to serious computational errors (Balachandar, 2009; Eaton, 2009). It should be reminded that, even though the local perturbation of the flow at the scale of the droplet size is ignored in the LPP model, the LPP droplets and the resolved fluid flow are still two-way coupled. Since the force exerted on each individual LPP droplet is subtracted from the resolved fluid flow, the presence of LPP droplets will influence the resolved flow. However, this influence is only at a spatial scale that is much larger than the droplet size (Balachandar, 2009). In the conversion from a resolved droplet to an LPP droplet (RD-to-LPP) the phase-characteristic function v in the region occupied by the droplet is changed from the value v ¼ 1 of the liquid phase to v ¼ 0 of the gas, as in Fig. 2. Moreover, the local perturbation induced by the droplet in the fluid velocity field should be removed carefully, as depicted in Fig. 2, so that the ‘‘undisturbed fluid velocity’’ is available to compute the fluid-dynamic force acting on the LPP droplet. The inverse conversion (LPP-to-RD) requires a similar treatment. In recent years, more sophisticated LPP models have been proposed for larger-than-grid-size solid particles, such as the Force-Coupling method (FCM) by Maxey and his colleagues (Maxey et al., 1997; Maxey and Patel, 2001; Lomholt et al., 2002; Lomholt and Maxey, 2003). Instead of computing the actual fluid-dynamic force exerted on the particle, a Gaussian distribution of external force monopoles and dipoles is applied to mimic the local perturbation. In such a case, the undisturbed fluid velocity is not required to model the interphase coupling. However, FCM is developed based on the multipole expansion of Stokes flows. Extension to droplets with finite Reynolds number is difficult.

χ=0

u

χ=1

up

χ=0

u

χ=0 up

Resolved Drop

Lag Point-Particle

Fig. 2. Schematic representation of a resolved droplet and a Lagrangian pointparticle droplet.

Even for droplets with small Reynolds number, FCM typically requires four to six grid points per particle diameter (Lomholt and Maxey, 2003). It is thus inapplicable in atomization simulations, since we need to deal with droplet resolution dp =Dx going continuously from the very small to four grid points. The goal of the present study is to propose a multiphase flow model for atomization simulations, where the large-scale interfaces are resolved by the VOF method and the small droplets by the LPP model. To focus on the development and validation of the LPP model, in particular for droplets that are larger than the grid size, we perform the study on a fixed grid. The application of the LPP model in an octree AMR framework will be considered in future works. When an adaptive grid is used, the cell size can be increased (gradually) by at least an order of magnitude, after the small resolved droplets are converted to LPP droplets. In such a case, the total computational cost can be reduced significantly. For a fixed grid, even though the reduction of computational time is not striking, there are several other important advantages to represent small poorly resolved droplets by the LPP model. First, the LPP model typically yields more accurate results on droplet dynamics than keeping the droplets resolved by VOF with poor resolution (dp < 4Dx). Second, the poorly resolved droplets are usually removed in VOF simulations (Shinjo and Umemura, 2010). An unaccounted-for removal of the small droplets formed makes the measurements of probability density function (PDF) of droplets impossible and also leads to a loss of physics. With the LPP model, these droplets can be kept as LPP droplets. Third, the poorly-resolved droplets sometimes cause numerical instability (which is the motivation to remove them), converting them to LPP improves the numerical stability of the simulation and speeds up the convergence of Poisson solver slightly. The multiphase flow models, corresponding equations and numerical methods are presented in Section ‘Modeling and numerical methods’. A series of multiphase flow problems are considered to validate the approach, and the results obtained by a number of simulations are compared against DNS results and experimental data in Section ‘Results’. In particular, the model is applied to simulate the gas-assisted atomization experiment by Descamps et al. (2008). This experiment is of particular interest because the drop trajectories and drop flying angle statistics are measured. The

125

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

numerical results are compared to the experimental measurements for a quantitative validation. Finally, we summarize the conclusions of the present study in Section ‘Conclusions’. Modeling and numerical methods

NS Solver

Our numerical approach is made up of two main components. The first one for the resolved flow contains a VOF method similar to that implemented in several open-source codes such as SURFER (Lafaurie et al., 1994), Gerris (Popinet, 2003) and Basilisk (Popinet). The other one is the Lagrangian point-particle (LPP) method, which is the main topic of this paper. The combined approach has been implemented in the more recent free code PARIS-Simulator (Arrufat et al.). The LPP model can be viewed as a local spatial filtering on top of the small droplets with a variable filtering length scale. The LPP droplets are those with a diameter smaller than the cut-off droplet size in the flow field. The cut-off droplet size is generally a function of the grid resolution, say dcut  ð4—6Þ Dx, since droplets with a diameter smaller than dcut are poorly resolved. Furthermore, droplets should also satisfy the physical requirement of an associated small Weber number, Wep , so that they remain almost spherical. The filtered flow contains both phases and is governed by the incompressible, variable-density, Navier–Stokes equations, with the interface evolution described by the VOF method. However, the LPP model and the coupling techniques here presented can be applied to other interface methods, such as Level-Set and Front-Tracking. The fluid flow around and inside the droplets is not resolved, and the LPP force model is used both to calculate the droplet motion and as a closure model to account for the effects of the droplet dynamics on the resolved or filtered flow field. The flowchart of the combined algorithm is shown in Fig. 3. For the resolved flow, the Navier–Stokes (NS) equations and the color function advection equation are first advanced in time. Then the separated liquid structures are identified and tagged by examining the volume fraction C field. The resolved droplets that are qualified for LPP are then submitted to the RD-to-LPP conversion routine. On the other side, the equation of motion for each LPP droplet is solved to update its velocity and position. The LPP droplets are then examined to check if they are still qualified to remain particles, if not they are converted back to resolved droplets. While the resolved flow properties are needed to integrate the LPP equation of motion, the force exerted by the flow on the LPP droplets should be subtracted from the flow momentum equation. Thus the NS and LPP solvers are two-way coupled. A description of the various solvers and routines is given in the following sections. Governing equations and numerical methods for the resolved flow Governing equations The one-fluid approach is employed to compute the resolved two-phase flow, where the liquid and gas phases are treated as one fluid with material properties that change abruptly across the interface. The incompressible, variable-density, Navier–Stokes equations with surface tension are

qð@ t u þ u  ruÞ ¼ rp þ r  ð2lDÞ þ rjds n  f p ; r  u ¼ 0;

Start

ð1Þ ð2Þ

where q and l are the fluid density and viscosity, u and p the velocity and pressure fields, and D the deformation tensor with components Dij ¼ ð@ i uj þ @ j ui Þ=2. The volume fraction C is introduced to distinguish the different phases, in particular C ¼ 1 in the computational cells with only the liquid phase, and its time evolution satisfies the advection equation

RF-LPP Coupling

LPP Solver

Color Function Solver (VOF) LPP-to-RD No Conversion

Tag Drops

Qualify LPP? Yes

Qualify LPP?

Yes RD-to-LPP Conversion

No No

RF: Resolved Flow RD: Resolved Drop LPP: Lagrangian Point-Particle

t>t-End? Yes End

Fig. 3. Flow chart of the combined VOF-LPP algorithm.

@ t C þ u  rC ¼ 0:

ð3Þ

The fluid density and viscosity are then defined by

q ¼ C ql þ ð1  CÞqg ; l1 ¼ C l1 þ ð1  CÞl1 or l ¼ C ll þ ð1  CÞlg ; l g

ð4Þ ð5Þ

where the variables with subscripts ‘‘g’’ and ‘‘l’’ represent properties of gas and liquid, respectively. The variables without a subscript are corresponding to the VOF representation, which denote variables of gas, liquid, and gas–liquid mixture for C equal to 0, 1, and a fractional number, respectively. The arithmetic mean is used for density, while the arithmetic or harmonic means are used for viscosity according to the physical problem (Boeck et al., 2007; Vincent et al., 2014). Furthermore, the present approach assumes that the velocity of the two phases match at the interface, which is valid for incompressible viscous fluids without phase change (Tryggvason et al., 2011; Rudman, 1998; Le Chenadec and Pitsch, 2013). The third term on the right hand side of Eq. (1) is a singular term, with a Dirac distribution function ds localized on the interface, and it represents the surface tension force. The surface tension coefficient is r, and j and n are the local curvature and unit normal of the interface. Finally, the last term of Eq. (1), f p , is the closure term that accounts for the backward effect of the LPP droplets on the resolved flow. Numerical methods Navier–Stokes solver. The Navier–Stokes equations, Eqs. (1) and (2), are solved by the projection method (Chorin, 1968). The time integration is conducted by a second-order predictor–corrector method (Tryggvason et al., 2011). The finite-volume approach on a regular, cubic staggered grid is utilized for spatial discretization. The advection term is discretized by the third-order QUICK scheme (Leonard, 1979). The viscous term is treated explicitly with the standard second-order centered difference scheme. In the projection step, an elliptic equation for the pressure is solved by a red– black successive overrelaxation (SOR) Gauss–Seidel iteration method, which is found to be competitive with a more complex

126

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

multigrid solver, probably because of the small time step constraint. Volume-of-Fluid method. The piecewise-linear geometrical Volume-of-Fluid method is used to solve the color function advection equation Eq. (3) (Scardovelli and Zaleski, 1999). The overall VOF method proceeds in two steps: interface reconstruction and then its advection. The interface reconstruction is performed using a piecewise-planar interface representation in each cell. The interface normal is computed by the Mixed-Youngs-Centered (MYC) method (Aulisa et al., 2007). After the reconstruction, direction-split geometrical fluxes are computed to advect the interface (DeBar, 1974; Rudman, 1997). The VOF method has been shown to preserve sharp interfaces and be close to second-order accuracy for practical applications (Aulisa et al., 2007; Popinet, 2009; Bornia et al., 2011). Curvature and balanced-force surface tension calculation. The height-function (HF) method is employed to calculate the local interface curvature, in a manner similar but not identical to the method implemented in Gerris (Popinet, 2009). We compute height function values whenever a 1D block of cells exists with a full cell on one side and an empty one on the other side. The block should be at most of length N H ¼ 9 cells. We then progress with a hierarchical sequence of steps: (i) when all the height values in a 3  3 stencil that includes the cell under examination are present, the curvature is computed by a finite difference calculation (Popinet, 2009); (ii) when only 7 or 8 heights are found, the resulting interface points are fitted by an elliptic paraboloid that approximates the shape of the local interface. However, the height values must be of the same kind, which means no mixing of the heights along the z direction, with widths or depths along x and y. The derivatives of the height function are then given in terms of the coefficients of the paraboloid; (iii) when less than 7 or 8 heights are found, we use a mixture of heights, widths and depths to get a set of at least 7 points. This set is trimmed to eliminate points that are too close to each other: when two points are less than half a cell size apart, only one of the two points is kept. Unlike (Popinet, 2009), the coordinate system is not rotated and the points are directly fitted to an elliptic paraboloid with the axis along the direction with the largest normal component; (iv) when the previous option fails, all heights, widths and depths are discarded and we use the centroid of the interface fragments reconstructed with the VOF method to create a set of points, which are then fitted to an elliptic paraboloid. As in several other works (Renardy and Renardy, 2002; Francois et al., 2006; Popinet, 2009), a balanced-force surface tension discretization is employed. We have checked that the method produces accurate curvatures in 2D and 3D test cases and well–behaved droplet oscillations in air–water systems.

  ~ g  up qg Du~ g C m qg Du~ g dup dup u þ g; ¼ /þ þ  dt sp qp Dt qp Dt dt

~ g the local undisturbed gas where up is the particle velocity and u velocity. The variables with a subscript p are corresponding to the LPP droplets. In the present study the LPP model is only used to represent liquid droplets, therefore physical properties denoted by the subscript p correspond to those of the liquid phase, for example qp ¼ ql . Furthermore, the liquid volume fraction is identical to zero, i.e., C ¼ 0, at the cell where an LPP droplet locates, as the liquid is represented by the LPP model. Therefore, the resolved flow properties such as u; q, and l at the LPP droplet location are all corresponding to the gas phase, i.e., ujx¼xp  ug ; qjx¼xp  qg , and

ljx¼xp  lg . The response time, density, and added-mass coefficient of the LPP droplet are denoted by sp ; qp , and C m , respectively. The modified gravity acceleration is denoted by g ¼ ð1  qg =qp Þg 0 , where g 0 is the gravity acceleration. The droplet Reynolds and Weber numbers can then be defined as

qg dp ju~  up j ; lg qg dp ju~  up j2 Wep ¼ : r

Rep ¼

ð7Þ ð8Þ

The position xp of an LPP particle is then updated by integrating

dxp ¼ up : dt

ð9Þ

The force terms on the right hand side of Eq. (6) represent the quasi-steady force, the pressure-gradient force, the added-mass force, and the gravity force, respectively, while the Basset-history force has been ignored for simplicity (Maxey and Riley, 1983; Ling et al., 2013). The relative importance in turbulent flows of the unsteady forces compared to the quasi-steady force is investigated by a scaling argument in Ling et al. (2013). The unsteady forces are necessary to accurately capture the LPP droplet motion when the LPP droplet size becomes comparable to the Kolmogorov scale and the density ratio qg =qp is significantly large. In the present study, we consider droplets with a rather small ratio qg =qp , therefore it is expected the contribution of the unsteady forces to be small. Nevertheless, the pressure-gradient and added-mass forces are retained to capture the leading order effect of the unsteady mechanisms. The quasi-steady force in Eq. (6) is written as the Stokes drag multiplied by a correction function /, which takes into account the effect of finite Rep . In the case of liquid droplets, for which ll  lg , the standard drag correlation for solid particle can be employed (Clift and Gauvin, 1970)

Lagrangian point-particle model for the unresolved small droplets

/ðRep Þ ¼ 1 þ The Lagrangian point-particle approach has been widely used in turbulent dispersed multiphase flows (Balachandar, 2009; Balachandar and Eaton, 2010). The key idea of the LPP model is to approximate the dispersed phase, made up of small solid particles, droplets or bubbles, as point masses. As a result, the flow at the scale of these particles is not resolved. To accurately track the LPP droplets in the Lagrangian framework, the force exerted on each individual droplet is calculated in terms of the undisturbed flow field properties. The closure is typically given by the force model or the so–called equation of motion (EOM) (Crowe et al., 1998; Maxey and Riley, 1983; Gatignol, 1983)

ð6Þ

0:15Re0:687 p

The response time et al., 1978)

sp ¼

qp d2p : 18 lg

þ 0:0175Rep 1 þ

42500 Re1:16 p

!1 :

ð10Þ

sp of the LPP droplet is expressed as (Clift

ð11Þ

A second-order predictor–corrector method is used for the time integration of the particle velocity and position, Eqs. (6)–(9), which is consistent with the algorithms used for the resolved flow. In order to focus on the development of the VOF-LPP model, droplet–droplet interaction (i.e., four-way coupling effect) is

127

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

neglected in the present LPP model. For the gas assisted atomization configuration of interest here, the small droplets formed in atomization are dispersed downstream of the liquid jet. As a result, their total volume fraction is quite low and the inter-droplet interaction is of secondary importance. Nevertheless, there exist atomization applications where the effects of inter-droplet interaction and secondary breakup are important (Apte et al., 2003), it is interesting to extend the present model to include these effects in future works.

The force exerted by the resolved flow on the LPP droplets is referred to as forward coupling and its computation is based on the resolved flow properties, Eq. (6). Due to Newton’s third law, the force exerted on the LPP droplets needs to be subtracted from the resolved flow and is referred to as backward coupling. This forcing term is represented by f p in the resolved flow momentum equation, Eq. (1). In the conventional LPP model, the computational cell is usually much larger than the LPP droplet size and many droplets are contained in each cell. The forcing term f p is then computed by a volume average over the LPP droplets N

p;c 1 X F fp;i ; V c i¼1

ð12Þ

where V c and N p;c are the cell volume and the number of droplets in that cell. Note that only the force due to fluid-LPP-droplet interaction needs to be subtracted from the momentum equation

F fp;i ¼ mp;i

  dup;i g : dt

fp ¼

Np X F fp;i Gðx  xp;i Þ;

ð13Þ

Although there exist issues of coupling LPP back to the resolved flow due to the number density fluctuation (Sundaram and Collins, 1996; Ling et al., 2010), the handling of two-way coupling is relatively easy. Since the fluid velocity corresponds to an average value over the cell, and the cell contains a large number of LPP droplets, then their influence is automatically filtered out (but the effects of the ambient flow and other LPP droplets outside the cell are retained). As a result, the stored fluid velocity can be used directly as the undisturbed fluid velocity in Eq. (6). However, in the present LPP model the droplets are bigger than the cell size, the two-way coupling becomes more complicated and in particular Eq. (12) is not valid any more. What’s even worse is the fact that if the force F fp;i is applied to the cell where the mass center of the LPP droplet is located, or is distributed to all the cells occupied by the LPP droplet, the resolved flow will be affected significantly. For this reason, we retain Eq. (6) and, as suggested by Maxey et al. (1997), write the resolved velocity ug as the sum of two con~g þ u ^ g . The first term u ~ g represents the undistributions, ug ¼ u ^ g is the velocity turbed fluid velocity, the second one u disturbance induced by the droplet under examination, a term which is filtered out in the conventional LPP model. When there are only a few LPP droplets in the flow field, their contribution to the momentum equation of the resolved flow is small, and a simple solution is to neglect the coupling force term f p in Eq. (1). This is the one-way coupling approximation for which ~ g ¼ ug . It has been found in previous works that the ^ g ¼ 0 and u u effect of the coupling force on changing the fluid momentum is related to the mass fraction of dispersed phase (Crowe et al., 1998; Ling et al., 2013), hence in the one-way coupling assumption ~ g increases with the mass fraction of the LPP droplets the error in u (Eaton, 2009).

ð14Þ

i¼1

where N p is now the total number of LPP droplets in the flow field. The Gaussian function Gðx  xp;i Þ is a numerical representation of the LPP droplet coupling force (Maxey et al., 1997)

GðxÞ ¼ ð2pL2 Þ

Momentum coupling between LPP and resolved flow

fp ¼

~ g , we propose to calcuIn order to obtain a better estimate of u late f p as

3=2

expðjxj2 =2L2 Þ;

ð15Þ

where L controls the size of the region where the force should be distributed. Eq. (14) is similar to that given by Maxey et al. (1997), but the interpretation is different. In the present model, L is viewed as the length scale of a spatial filter of the coupling force at the LPP droplet location. The parameter L is chosen to be larger than the actual droplet size. Depending on the volume fraction of the LPP droplets, L is picked between 5 and 10 dp . Therefore, the th

coupling force F fp;i corresponding to the i LPP droplet is smoothly distributed within a region larger than the droplet size, rather than concentrated at the droplet location. At every LPP droplet location, the contribution of the force from the present droplet to the local resolved flow is small, so that ug computed from Eq. (1) can be used ~ g in Eq. (6). as an approximation of u There are several advantages for treating the coupling force in the above way. First, the undisturbed fluid velocity is directly given by the resolved flow momentum equation to calculate the droplet force. Second, the integration of Eq. (14) over the whole domain yields the total force exerted on all the LPP droplets and the total momentum of the combined RF-LPP system is conserved. Third, since L is determined by the droplet diameter and not by the cell size, the coupling between RF and LPP droplets does not depend on the grid. This feature is quite important if the LPP model is coupled with AMR, because the calculation of the force on the LPP droplets and the backward effect on the resolved flow will not be affected when the mesh adapts itself from the well-resolved droplet regime, Dx < 0:1dp , to the conventional LPP model, Dx > 10dp . The present treatment of the coupling force will only capture the backward effect of the LPP droplets on the resolved flow on length scales larger than L. The dynamics due to fluid-droplet interaction on scales smaller than L, such as the shear layer near the droplet interface, the wake of the droplet and the interaction with the resolved flow, will be missed. Finally, since the LPP droplets do not in general coincide with the locations at which the fluid velocity is stored, interpolation is required to compute ug at the LPP droplet location, and a tri-linear interpolation scheme is used in the present study (Apte et al., 2003; Ling et al., 2010). When the spatial variation of the resolved flows is significant over the LPP droplet size, then picking the fluid velocity at the LPP droplet center is not sufficient and the Faxén correction is often used to generalize the force models (Faxén, 1922). Furthermore, if the LPP droplet moves in a shear flow, then the shear lift force will be induced (Legendre and Magnaudet, 1997). Nevertheless, the LPP model is only applied to small droplets in the present study, the effect of inhomogeneous ambient flows on the droplet force is expected to be generally small and thus the Faxén and lift forces are neglected. Two-way conversion between LPP and resolved droplet The flow chart of the combined algorithm in Fig. 3, shows the tight coupling between the two solvers for the resolved flow and the LPP droplets. In atomization simulations that consider the jet formation as well, the initial flow is completely resolved with the VOF method, since small droplets are not yet formed. Once the

128

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

liquid jet starts to atomize, the droplets away from the interface are converted to LPP droplets. Some of the LPP droplets may follow trajectories that bring them back toward the interface. These droplets may be converted back to resolved droplets. The conversion algorithm is summarized in Fig. 4. Conversion criteria The criteria to determine whether a droplet should be represented by the LPP model or remain resolved depend on its shape and position. The present LPP force model assumes the droplet to be spherical. When the droplet Weber and Reynolds are large, or when the droplet Capillary number Cap ¼ Wep =Rep becomes large the droplet may deform significantly and the LPP model will introduce an error in its motion. However, if the droplet diameter is less than four to six grid spacings, then the shape and the dynamics of the droplet are not accurately represented even by the VOF method. Therefore, droplets with moderate Weber and Reynolds numbers should also be converted to LPP droplets. In previous atomization simulations these poorly resolved droplets, with dp 6 4 Dx, were removed to improve the stability of the solver (Shinjo and Umemura, 2010). Here, instead of simply removing them, we convert them to LPP droplets, with a positive effect on the numerical stability and on the mass conservation, and with better droplet statistics during the atomization process. For convenience of the calculation, the droplet volume V p is used instead of its diameter dp to define the droplet size. The aspect ratio c is also estimated in order to avoid converting to LPP droplets thin ligaments and sheets. Therefore, only when V p is smaller than a critical value, say V crit ’ ð4DxÞ3 , and c is close to one, then the droplet is qualified for the LPP model. Furthermore, a location criterion is used to check if the droplet is away from the liquid jet interface. Since the current LPP model does not include either the formation of a droplet or a droplet impact on the interface, only droplets that are at a certain distance away from the liquid jet interface can be converted. The distance is typically chosen to be the droplet diameter dp . The overall conversion criteria can be expressed as

for icell =1,Ncell do Identify separated droplets and generate RD array Perform droplet calculation (droplet volume, etc) end for

Tagging

RD-to-LPP Conversion

for iRD =1,NRD do Check conversion criteria if Droplet iRD qualified for LPP Change non-zero color function to zero Replace perturbed flow field by undisturbed flow field Add droplet iRD into LPP array and store droplet position and velocity end if end for

LPP-to-RD Conversion

for iLPP =1,NLPP do Check conversion criteria if Droplet iLPP not qualified for LPP Reconstruct color function for the droplet Change velocity within the droplet Remove droplet iLPP from LPP array end if end for

Fig. 4. Algorithm for conversation between resolved droplets (RD) and Lagrangian point-particle (LPP) droplets.

=0

u

Reconstructed undisturbed flow field

=0 u

u p= (∫u dV)/Vp

RD-to-LPP

=1

Resolved Drop

(a)

=0

u up+u’p =1

Lag Point-Particle

=0

u =0 LPP-to-RD

up

fLPP Qualifiedg ¼ fV p < V cut g && fjcp  1j < c g && fxp 2 Rai g; ð16Þ where c is the tolerance for the aspect ratio, and Rai is the region away from the interface. RD-to-LPP conversion For the conversion from resolved droplets to LPP droplets, the separated liquid structures must be first identified by tagging the color function field. The tagging approach of Herrmann (2010) is employed here, and cells that are attached to each other, with respect to the liquid phase, will have the same tag number. During the tagging process, calculations of the droplet properties, including its volume, aspect ratio, location and velocity of the center of mass, are also performed. A schematic representation of the conversion from RD to LPP is shown in Fig. 5(a). In the conversion, the droplet volume V p and velocity up are first computed by an integration over the region X occupied by the droplet (gray area of Fig. 5(a)), namely R R V p ¼ X v dV and V p up ¼ X u v dV, the phase-characteristic function v is then changed from the value 1 to 0. Moreover, the perturbed flow field in the same region should be replaced by the undisturbed flow, so that the LPP droplet sees the proper undisturbed fluid velocity. This is done by interpolating each component of the fluid velocity along the corresponding coordinate from the surface into the interior of a cubic region centered at the particle location. The edge length of the cubic region is about

Resolved Drop

(b)

Lag Point-Particle

Fig. 5. Conversion between a resolved droplet (RD) and a Lagrangian point-particle (LPP) droplet. The fluid flow field for the LPP droplet is slightly perturbed due to the coupling force.

two times the droplet diameter (blue1 area of Fig. 5(a)). For example, for the velocity component u along the x-direction, a linear interpolation is performed between the two values uði1  1; j; kÞ and uði2 ; j; kÞ, where i1 and i2 are the indexes of the left and right boundary cells of the cubic region. The reconstructed velocity field is globally divergence-free if the velocity on the surface of the interpolation region is divergence-free. The field is also divergence-free locally if the divergence components @u=@x; @ v =@y, and @w=@z are constant within the region. Unless the particle Reynolds number is very small, a cubic region with an edge size twice the droplet diameter is sufficient to reconstruct the undisturbed flow field in a satisfactory way. The reconstruction of the undisturbed flow field is necessary even for the simple approach of removing the small poorly-resolved droplets from

1 For interpretation of color in Figs. 5, 10, 11, 14 and 16, the reader is referred to the web version of this article.

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

the simulation, since which minimizes the influence of the removal of droplets to the resolved flows. It has been observed in the previous work that if the small liquid structures are replaced by gas with the same velocity, the interface instability in the resolved flow is modified (Shinjo and Umemura, 2010). Since in the RD-to-LPP conversion, the region originally occupied by the liquid will be filled with the gas phase, it might seem that additional gas is added to the system and mass conservation is violated. As a matter of fact, the added gas is actually considered as a ‘‘ghost fluid’’ and is used only for the droplet force calculation. When an LPP droplet is converted back to a resolved droplet, the ghost fluid will be taken out from the system. Therefore, the ghost fluid should always be excluded from the calculation of the total mass of the system. Furthermore, the LPP droplet represents not only the momentum of the droplet, but also the momentum of the unresolved perturbed flow around the droplet. The latter can be viewed as the momentum of a virtual volume of fluid that moves with the LPP droplet. Therefore, although the perturbed flow around the droplet is removed in the RD-to-LPP conversion, the corresponding momentum is kept. As a result, the momentum of the overall system is conserved. LPP-to-RD conversion Several conversion criteria can be devised to examine if an LPP droplet is still qualified for the LPP model, for example the droplet volume can increase due to collision and merging of droplets or to phase change. These physical processes are not considered in the present LPP model, where the only condition for converting an LPP droplet to a resolved one is that the droplet position is very close to the liquid jet interface. A schematic picture of the conversion process is shown in Fig. 5(b), and the algorithm is described in Fig. 4. First, a spherical droplet is rebuilt around xp , by updating properly the color function in the cells that will be occupied by the resolved droplet. The ghost fluid that was added in the inverse RD-to-LPP conversion is now taken out. In the same cells containing the new resolved droplet the velocity field needs to be updated to take care of the LPP droplet velocity. In principle, the local perturbed flow inside and outside the droplet should be reconstructed, this being the inverse process of the removal of the perturbed flow in the RD-to-LPP conversion. The reconstruction of the local perturbed flow is quite a difficult task, except in the limit of zero Reynolds number for which an analytical solution exists (however, in that case the perturbed region is infinitely large). Nevertheless, if one naively ignores the effect of the local perturbed flow and only changes the velocity inside the droplet to be the LPP droplet velocity up , then it can be easily shown that the droplet velocity will have a significant transition jump right after the LPP-to-RD conversion. This is simply due to the fact that a portion of the momentum of the droplet given in the conversion is transferred to the surrounding flow and used for the local perturbed flow development. As described in the RD-to-LPP conversion, the LPP droplet actually contains the momentum of both the droplet and the virtual fluid moving with the droplet. In the naive conversion approach described above, only the portion of droplet is included. In order to compensate the perturbed flow that is not reconstructed, the momentum of the virtual fluid also needs to be added back. Therefore, the velocity in the cells occupied by the resolved droplet should be changed to u0p for the purpose of compensating the momentum in the perturbed flow field around the droplet, i.e.,

qp V p ðu0p  u~ g Þ ¼ qp V p ðup  u~ g Þ þ

Z Vp

qf ðug  u~ g ÞdV;

ð17Þ

129

The second term on the right hand side is the momentum of the perturbed flow field, which is approximated by assuming a virtual volume of fluid, aV p , moving with the droplet velocity,

Z Vp

qg ðug  u~ g ÞdV  aqg V p ðup  u~ g Þ:

ð18Þ

Then the velocity to be used in the cells of the resolved droplet in the LPP-to-RD conversion can be expressed as

!

u0p

qg  ug ¼ 1 þ a ðup  ug Þ: qp

ð19Þ

The physical meaning of a is the ratio between the volumes of the virtual fluid and the droplet. However, it is not easy to get a precise value of a. It is known that, in order to accelerate a particle, additional forces, i.e., the added-mass and history (viscous-unsteady) forces, are needed to accelerate the ambient fluid around the particle. For the added-mass effect, it is considered that the ratio between the volumes of the virtual fluid to be accelerated through the inviscid mechanism and the particle is the added-mass coefficient C M , which is equal to 0.5 for a sphere in incompressible flows. Due to finite viscous diffusion time scale, the history force usually needs to be expressed in integral form (Basset, 1888; Mei and Adrian, 1992). Nevertheless, it has been shown by Ling et al. (2013) that, if the particle and ambient fluid acceleration time scales are much larger than the viscous diffusion time scale, the history force can also be expressed as non-integral form like the added-mass force, and a viscous-unsteady coefficient C v u similar to the added-mass coefficient C M can be derived as

  0:75 þ 0:105Rep : C v u  8:51 Rep

ð20Þ

The viscous-unsteady coefficient C v u , can be considered as the ratio between the volumes of the virtual fluid to be accelerated through the viscous mechanism and the particle. It is shown that C v u decreases with the droplet Reynolds number Rep . When Rep varies from 10 to 100, C v u decreases from 1.53 to 0.96. The excess momentum added through u0p is to mimic the effects of the added-mass and history forces on accelerating the surrounding fluid around the droplet. Therefore, a can be estimated as the sum of C M and C v u and thus depends on the droplet Reynolds number as well. In an actual computation, the resolution of RD is usually 4–6 grid cells across the diameter, and the local perturbed flow may not be computed very accurately. As a result, a usually needs to be adjusted slightly. For simplicity, a constant of a ¼ 4 is used in the present study, which has ben shown to yield a smooth droplet velocity transition in the LPP-to-RD conversion.

Results To validate the combined multiphase flow model where the smallest droplets are represented by the LPP model, four test problems with increasing complexity are conducted: settling droplet, droplet in a lid-driven cavity, pulsed jet atomization, and gas-assisted atomization. The combined model considers the point-particle approximation for small droplets, that introduces errors even when the droplets are spherical. Nevertheless, these errors should be small and should not influence the features of interest to be captured by atomization simulations, such as droplet size PDF and droplet flying angle PDF. To validate the combined model, its numerical results are compared with experimental data or DNS, depending on the test problem.

130

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Settling droplet

Table 1 Parameters for the settling droplet problem.

The combined model is first tested with the simulation of a liquid droplet settling in a large rectangular cuboid filled with a stationary gas, see Fig. 6(a). The relevant physical properties of the two phases and other parameters are given in Table 1. The dimensions of the tank are Lx ¼ Lz ¼ Ly =2 ¼ 103 m. Gravity is along the negative y-direction. The center of the droplet is initially at (0:5 Lx ; 0:5 Ly ; 0:5 Lx ) and reaches the settling velocity at about time t ¼ 0:03 s. The droplet terminal velocity is about v t ¼ 0:04 m/s with a corresponding droplet Reynolds number Ret  0:4. The droplet remains about spherical during the settling due to viscous effects, and surface tension is therefore not included in this simulation. Note that the standard drag correlation Eq. (10) is for LPP in an unbounded domain. Here, the tank width and depth are only 10 times the droplet diameter and, as a result, the tank wall has a noticeable retarding effect on the settling velocity v t . The relation between the terminal settling velocity v t and the diameter-to-width ratio dp =Lx has been studied by many researchers, typically in a cylindrical tube with a circular cross section (Di Felice, 1996). To include the wall effect, the correlation of Di Felice (1996)

g¼

vt v t;unbounded

 ¼

1  dp =Lx 1  0:33 dp =Lx

ð21Þ

is used to correct the drag calculation, where g is defined by the ratio between the terminal settling velocity in bounded and unbounded domains. This correlation is here used to correct the quasi-steady force in the equation of motion (6) in the following way

  ~ g  up / qg Du ~ g C m qg Du ~ g dup dup u þ g: ¼ þ þ  dt sp g qp Dt qp Dt dt

qg (kg/m3)

1000

10

10

3

lg (Pa s) 4

10

dp (m) 10

4

g (m/s2)

0

9.8

In the second test problem, we investigate the droplet motion in a 3D lid-driven cavity. The cavity is a cube of size L and the flow is driven by a moving lid on the upper boundary with velocity U 0 ; the Reynolds number of the cavity is Rec ¼ 100. The boundary conditions along the x- and y-directions are no-slip, periodic boundary conditions are considered along the third z-direction. The relevant physical properties of the two phases and other parameters are given in Table 2. Simulations are first run without any droplet, then the induced flow remains 2D, moving towards the steady state condition shown in Fig. 7. In Fig. 7(b), the u component of the velocity (x-direction) along the vertical centerline y ¼ z ¼ L=2,

In DNS without the LPP model and with cubic cells of side h, four different resolutions, dp =h ¼ 2; 4; 8; 16, are used to resolve the droplet. With the LPP model there is only one droplet in the tank, hence the one-way coupling approximation is made and f p ~ g and Du ~ g =Dt are equal to zero. is ignored in Eq. (1). As a result, u The LPP model result is simply obtained by integrating Eq. (22) and thus it is independent of the cell size. The results for the y-velocity of the droplet with only the VOF method or the LPP model are shown in Fig. 6(b). It can be observed that the LPP result agrees well with the DNS results obtained with the finest meshes,

0

LPP DNS,dp/h= 2 DNS,dp/h= 4 DNS,dp/h= 8 DNS,dp/h=16

-0.01

-0.02

-0.03

-0.04

-0.05

0

0.005

0.01

0.015

0.02

0.025

0.03

t (s)

(a)

r (N/m)

Lid-driven cavity

ð22Þ

v (m/s)

ll (Pa s)

dp =h ¼ 8; 16, when the droplet is sufficiently resolved. When the mesh resolution is low, dp =h ¼ 2; 4, DNS yields errors of 16.2% and 7.0% in the terminal velocity, respectively. Therefore, in this case LPP model yields more accurate results than DNS of a poorly resolved droplet. Though there is no experimental data for this specific test case, there are empirical correlation for the terminal settling velocity to compare with. For a droplet settling in an unbounded domain, the terminal velocity obtained through the standard drag (Clift and Gauvin, 1970) is about v t ¼ 0:0493 m/s. This value is considerably larger than the simulation results, implying that the wall effect is significant and should be considered. As a function of the Reynolds number, in the two limits Rep ! 0 and Rep ! 1, the terminal velocity, including the wall retarding effect, can be computed as v t;0 ¼ 0:0393 m/s and v t;1 ¼ 0:0478 m/s (Di Felice, 1996). Both DNS and LPP results, with the exception of the lowest resolution dp =h ¼ 2, are bounded by these two limits. In particular, since the Reynolds number based on the droplet terminal velocity is about Rep ¼ 0:4, we expect and find that our results are closer to the lower limit. With the correlation of Di Felice (1996) the terminal velocity is about v t ¼ 0:0406 m/s, which again agrees reasonably well with our simulations. The small discrepancy may be also due to the different shape of the domain cross section in the present study and in Di Felice (1996).

2:7 ;

ql (kg/m3)

(b)

Fig. 6. A liquid droplet settling in a stationary gas tank: (a) problem setting, (b) settling velocity as a function of time.

131

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143 Table 2 Parameters for the 3D lid-driven cavity problem.

ql

(kg/m3) 3000

qg

ll

lg

L (m)

(Pa s)

(Pa s)

dp (m)

r (N/m)

(kg/m3)

U0 (m/s)

10

103

104

104

1:5  105

3:2  104

3.125

and the v-component (y-direction) along the horizontal centerline x ¼ z ¼ L=2 are compared with the classic work of Ghia et al. (1982). Good agreement is observed with the results obtained with the two grid resolutions L=h ¼ 64; 128. After the cavity flow reaches a stationary condition at about time t ¼ 1 ms, a single or multiple droplets are seeded into the cavity. We first consider the simpler case where only one droplet is seeded. As in the settling droplet problem, the feedback effect of the single LPP droplet to the resolved flow is ignored, hence f p ¼ 0. Four different computations are performed: DNS (RD), LPP, and two different implementation of the RD-LPP conversion. In the first one (RD-LPP-1), the complete conversion algorithm described in Section ‘Two-way conversion between LPP and resolved droplet’ and Fig. 4 is used; while in the second one (RD-LPP-2) the velocity remains unchanged in the conversion, namely the undisturbed and disturbed flow fields are not reconstructed in the two conversions LPP-to-RD and RD-to-LPP. The conversion criteria for this particular test case is that the droplet is represented by the LPP model inside the region ymin < yp < ymax . The initial position of the droplet is at xp ¼ ð0:61 L; 0:61 L; 0:59 LÞ, and two different initial velocities are

considered, up1 ¼ ð0:184 U 0 ; 0:0293 U 0 ; 0Þ and up2 ¼ ð0:0320 U 0 ; 0:0320 U 0 ; 0:0320 U 0 Þ. The first value is equal to the local fluid velocity, up1 ¼ uðxp Þ, the second value has a different magnitude and direction with respect to the local fluid velocity. With the initial velocity up1 the droplet will be represented by the LPP model in the slab defined by the two limits ymin ¼ 0:61 L and ymax ¼ 0:64 L, with up2 the limits are ymin ¼ 0:59 L and ymax ¼ 0:625 L. The simulations are performed with the mesh resolution L=h ¼ 128, corresponding to a ratio dp =h ¼ 4. The droplet is initially seeded into the driven cavity as an LPP droplet. For the case of DNS, the LPP droplet is converted to a resolved droplet right after it is seeded. To have a fair comparison between the LPP and DNS results, it is important the average velocity of the resolved droplet just after the LPP-to-RD conversion matches well with the LPP droplet velocity at the corresponding time. Different values of a are tested and the results are summarized Table 3. The superscript 1 and 2 denote the droplet velocity right before and 6 steps after the LPP-to-RD conversion. The droplet velocities for LPP-to-RD conversion with different a are represented by u2p;RD ðaÞ; while the LPP droplet at the corresponding time without a conversion is denoted by u2p;LPP . It can be seen that among the values tested, a ¼ 4 yields results closest to u2p;LPP . The time evolution of the droplet velocity for the two initial velocities is shown in Figs. 8 and 9. First, it should be observed that the DNS and LPP results agree reasonably well for both cases, as the difference in magnitude of the u component in Fig. 8(a) is about 10%. If this value is compared to the time variation of the v component in Fig. 8(b), then the difference between DNS and LPP is

N=128 N= 64 Ghia

(a) (b) Fig. 7. Single phase results of the lid-driven cavity problem: (a) velocity vector plot at z ¼ L=2, (b) u at y ¼ z ¼ L=2 and

v at x ¼ z ¼ L=2.

Table 3 The effect of the paremeter a on the droplet velocity in LPP-to-RD conversion. The superscript 1 and 2 denote the droplet velocity right before and 6 steps after the LPP-to-RD conversion. The droplet velocities for LPP-to-RD conversion with different a are represented by u2p;RD ðaÞ; while the LPP droplet at the corresponding time without a conversion is denoted by u2p;LPP . The lid velocity U 0 ¼ 3:125 m/s and the undisturbed fluid velocity at the LPP droplet location before conversion is (0:1839U 0 ; 0:02928U 0 ; 0).

u=U 0

v =U 0 w=U 0

u1p;LPP

u2p;LPP

u2p;RD ða ¼ 0Þ

u2p;RD ða ¼ 2Þ

u2p;RD ða ¼ 4Þ

u2p;RD ða ¼ 6Þ

0.03200 0.03200 0.03200

0.03190 0.03200 0.03198

0.02871 0.03377 0.03007

0.03013 0.03380 0.03027

0.03155 0.03383 0.03047

0.03297 0.03386 0.03067

132

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

-0.2 -0.25

0.6

DNS (RD) LPP RD-LPP-1 RD-LPP-2

0.5

-0.35

RD-to-LPP

-0.4 -0.45

RD-to-LPP

0.3 0.2

-0.5

0.1

-0.55 -0.6 0.001

DNS (RD) LPP RD-LPP-1 RD-LPP-2

0.4

v (m/s)

u (m/s)

-0.3

0.7

0

LPP-to-RD 0.00105

0.0011

0.00115

-0.1 0.001

0.0012

LPP-to-RD 0.00105

0.0011

t (s)

t (s)

(a)

(b)

0.00115

0.0012

Fig. 8. Time evolution of the velocity components when the droplet initial velocity is the same as the local fluid velocity: (a) u, (b) v.

0.1

0.3 DNS (RD) LPP RD-LPP-1 RD-LPP-2

0 -0.1

0.2

-0.3

v (m/s)

u (m/s)

-0.2

LPP-to-RD

-0.4 -0.5

DNS (RD) LPP RD-LPP-1 RD-LPP-2

0.1

0

LPP-to-RD

-0.1

-0.6 -0.7 0.001

0.0011

0.0012

0.0013

-0.2 0.001

0.0011

0.0012

t (s)

t (s)

(a)

(b)

0.0013

Fig. 9. Time evolution of the velocity components when the droplet initial velocity is different from the local fluid velocity: (a) u, (b) v.

indeed quite small. Furthermore, this error magnitude is also consistent with the droplet low resolution, dp =h ¼ 4. When we consider the conversion model with the initial velocity up1 we observe that the droplet at about time t ¼ 1 ms is on the lower boundary of the slab with a negative y-component of the velocity. The LPP particle is quickly converted to RD, then after a transient it changes the sign of this velocity component and finally it enters the slab region where it is converted back to LPP. A similar discussion holds for the initial velocity up2 . With the conversion model described in Section ‘Two-way conversion between LPP and resolved droplet’, the droplet velocity has a smooth transition after the conversions for both initial velocities up1 and up12 . In Figs. 8 and 9 this model is represented by the RD-LPP-1 lines, with the conversion instants marked by an arrow. If the perturbed flow field is not reconstructed properly in the LPP-to-RD conversion, then remarkable transition jumps are observed in the droplet velocity after the conversion (RD-LPP-2 lines). The transition jump is more significant when the droplet is not in velocity equilibrium with the local fluid. As the drag force tends to bring the droplet in velocity equilibrium with the fluid, the droplet velocity may slowly change toward the correct value, as in Fig. 8, but the transient process can be quite long. However, if the initial error in the conversion causes a significant deviation of the droplet trajectory, the droplet velocity will never go back to the correct value, as in Fig. 9(b). Therefore, it is

critical to handle properly the local flow field in the conversion process. We also consider the case with multiple droplets seeded in the cavity. In total 183 monodisperse droplets are seeded randomly after a steady state flow is reached in the cavity with a single phase. The mass fraction of the seeded droplets, defined by the ratio of their mass to the total mass of liquid droplets and gas, is about 0.75. In this case the mass fraction is large, and two-way coupling simulations are performed with the force coupling model described in Section ‘Momentum coupling between LPP and resolved flow’. The size of the parameter that controls the force distribution is L ¼ 10 dp . Fig. 10 shows the droplet distribution inside the cavity at times t ¼ 1:04; 1:08; 1:12 ms. When the droplets are resolved by the VOF method the interfaces are plotted as the intermediate contour level C ¼ 0:5, and are indicated by a red color. The droplet distribution represented by the LPP model is given by the green spheres. The DNS and LPP results overall match well on the droplet location. The differences, shown in the images on the third column, mainly occur when the droplets are subject to collisions or serious deformations, as near the top of the cavity where the shear is strong. This indicates the need for a future improvement of the model to take into account both droplet–droplet and droplet–wall interactions, such as incorporating the collision model of Schmidt and Rutland (2000).

133

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

(a)

(b)

(c)

DN S

LPP

DNS vs LPP 3

Fig. 10. Droplet distributions in the cavity at times: (a) t ¼ 1:04  10 , (b) t ¼ 1:08  103 , and (c) t ¼ 1:12  103 s.

It is also important to point out that the results obtained with the one-way and the two-way coupling models are indeed very close. Therefore, the one-way coupling with the present conversion model seems to be adequate even when the LPP droplet mass fraction is quite large. Nevertheless, further validation of the coupling model is still required in future works. Pulsed-jet atomization Next we consider the atomization of a dense cylindrical liquid jet which is injected into a stagnant gas tank, to examine the capability of the present model to simulate multiphase flows with very complex interfaces. The primary breakup of a liquid jet injected into stagnant gas has been investigated by many previous works in detail, see for example (Ménard et al., 2007; Lebas et al., 2009; Shinjo and Umemura, 2010). Discussions of the primary breakup mechanisms of the liquid jet can be found in the above reference. In the present paper, we mainly focus on testing the present model. The jet diameter is denoted by Djet and the tank dimensions are Ly ¼ Lz ¼ 4 Djet and Lx ¼ 10 Djet . The injection is along the positive x-direction. The inflow velocity u is modulated sinusoidally, u ¼ U jet ð1 þ 0:05 sinð20 p tÞÞ, to promote the growth of primary shear instabilities. The physical properties and parameters in non-dimensional form are listed in Table 4, where Rejet ¼ ql U jet Djet =ll and Wejet ¼ ql U 2jet Djet =r. On the left and right boundaries, inflow and outflow boundary conditions are considered, respectively, while free-slip conditions are applied on the other four boundaries. The grid is composed of cubic cells of side h ¼ Djet =64. Two simulations are performed for this problem: the first one is DNS with the VOF method, the second one with the combined method, with the small droplets described by the LPP model. Since the mass fraction of the LPP

Table 4 Parameters for the pulsed-jet atomization problem.

ql =qg

ll =lg

Rejet

Wejet

20

20

800

4000

droplets in the flow is small during all the simulation, and based on the knowledge of the previous test, the one-way coupling model is used. Since the amount of liquid injected into the domain increases with time, the total number of droplets and the number of those represented by LPP model also increase. The maximum number of LPP droplets is about 6000. The simulation of the pulsed-jet problem is conducted on the CALMIP machine, with a typical CPU time on 160 processors of 245 min for the combined model and 253 min for DNS. Snapshots of the interface configuration in the two simulations at time t ¼ 9 Djet =U jet are shown in Fig. 11(a) and (b), with the VOF interface represented in green and white for the contour level C ¼ 0:1, in order to make small droplets clearly visible. The orange spheres represent the liquid droplets of the LPP model, and their radius length is such to give the actual droplet volume. The two figures are in good agreement with each other, thus confirming the capability of the combined model to capture atomization. The primary breakup of the pulsed liquid jet can be clearly seen in Fig. 11. Liquid sheets are formed at the front of the liquid jet and also on the liquid core surface. At the front, a umbrella-shape sheet is formed due to the impingement of the liquid jet against the stagnant gas and the resulting lateral liquid spread. On the liquid core surface, it is seen that the perturbation near the inlet remains almost sinusoidal, but as it is advected downstream, the interface folds and thin liquid sheets are formed. The sheets incline toward the upstream direction, looking like an umbrella. As the sheets are

134

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Fig. 11. Snapshots of the liquid–gas interface of the pulsed-jet atomization at time t ¼ 9 Djet =U jet : (a) DNS, (b) the combined model, where the orange droplets are from the LPP model.

stretched, they break up into thin ligaments. The ligaments in sequence break into small droplets due to the Plateau-Rayleigh instability (Shinjo and Umemura, 2010). The size distribution (PDF) of the droplets formed in the pulsed-jet atomization at different times is shown in Fig. 12. Since more and more liquid is injected into the domain and then atomized to small droplets as time evolves, the droplet Sauter mean decreases with time, while the number of small droplets increases. It can be seen that the droplet size distributions obtained by DNS and the combined model are in good agreement. In particular, the combined model captures the peak of the PDF, while the difference of the Sauter means with the DNS results is less than 10%. When the LPP model is introduced, the number of tiny droplets is significantly higher than the DNS results without the LPP model, see Fig. 12. The discrepancy may be due to the erroneous dynamics of small droplets when they remain to be resolved by VOF. Small fragments of VOF-represented liquid tend to be reconstructed and advected erroneously, so that these small fragments stay where they are formed they are easily swept by and merged into the moving bulk liquid jet. If these small fragments are converted to LPP, then they can be advected away from the liquid jet and remain conserved. As a result, less small droplets are seen in the pure VOF simulation than the combined VOF-LPP simulation. The forward RD-to-LPP process and the backward LPP-to-RD conversion, from the ligament breakup to the droplet merging back in the liquid jet, are shown in Fig. 13. In particular, a ligament formed at the edge of the sheet is breaking due to Plateau-Rayleigh instability in Fig. 13(a). In Fig. 13(b), it can be seen that the initially formed droplet has a rather irregular shape. Under the action of surface tension, the droplet becomes closer to a spherical shape in Fig. 13(c). Then, since the droplet satisfies both the size and aspect-ratio criteria and is enough away from the jet, it is converted from RD to LPP, as shown in Fig. 13(d). The LPP droplet moves for a while in the gas, then driven by the recirculation generated on the inside of the umbrella-shape front, the droplet is carried back to the liquid jet. When the LPP droplet is getting close to the interface, it does not satisfy the criteria of LPP and then undergoes a backward conversion from LPP to RD in Fig. 13(e). Finally it is reabsorbed by the liquid jet in Fig. 13(f). The present results indicate that the two-way conversion scheme is capable to capture complex droplet dynamics in atomization. Gas-assisted atomization The combined model is finally applied to simulate the gas-assisted atomization experiment conducted by Descamps et al. (2008). This experiment is particularly interesting to the

present study, because the drop trajectories and drop flying angle statistics are measured in the experiment and thus can be used to validate the capability of the present combined model in accurately capturing the drop dynamics. Up to the authors’ knowledge, the experiment by Descamps et al. (2008) is the only one in the literature that provides detailed measurements of drop trajectories. Simulation setup The simulation setup is shown in Fig. 14. The atomizer consists of two planar jets: on top the faster gas jet with injection velocity U g , and on bottom the slower liquid jet with injection velocity U l . The inlet height for both phases is the same, and is denoted by H. As a function of this characteristic length, the width, height and thickness of the computational domain are Lx =8 ¼ Ly =4 ¼ Lz ¼ H, respectively. The domain is subdivided with cubic cells of side h ¼ H=64 for the coarse mesh and h ¼ H=128 for the fine one. At the inlet, a separator plate, indicated by the blue color in Fig. 14, has been positioned to separate the gas and liquid flows. The need for such a plate to accurately simulate atomization has been shown in previous papers (Fuster et al., 2013). The width and height of the separator plate are lx ¼ H=4 and ly ¼ H=32, respectively. Inflow and outflow boundary conditions are applied along the x-direction, slip conditions along the vertical y-direction, periodic boundary conditions along the z-direction. A snapshot of the liquid–gas interface resolved by the VOF method (in gray) and by the LPP model (in orange) is also shown in Fig. 14. The relevant physical properties of the two phases and other parameters are given in Table 5. The properties of the gas and liquid phases are chosen in such a way to guarantee that the important dimensionless parameters are as close as possible to the corresponding value of the experiment by Descamps et al. (2008). In particular, Table 6 shows that the Reynolds and Weber numbers of the liquid inflow, Rel ¼ ðql U l HÞ=ll and Wel ¼ ðql U 2l HÞ=r, and the momentum ratio between the gas and liquid inflows, M ¼ ðqg U 2g Þ=ðql U 2l Þ, are very close to the experimental value. On the other hand, while the numerical Reynolds number based on the gas inflow boundary layer thickness dg , Reg;d ¼ ðqg U g dg Þ=lg , is also very similar, the corresponding Weber number, Weg;d ¼ ðqg U 2g dg Þ=r, is about twice as big. Finally, the numerical Reynolds number of the gas jet is fairly large, Reg ¼ ðqg U g HÞ=lg  8000, but still only about 42% of the experimental value. For such high Reynolds numbers, it is extremely expensive to resolve all turbulent scales in the gas flow, and we do not expect that the current mesh resolution can capture the smallest turbulent eddies. Therefore, the term ‘‘DNS’’ here refers

135

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

10 3

10 3

DNS Model

10 1

counts 10 1

d32,DNS=2.60E-3 d32,LPP=2.88E-3

10 0 0

0.005

0.01

10 1 d32,DNS=2.55E-3 d32,LPP=2.66E-3

0.015

0.02

DNS Model

10 2

counts

10 2

counts

10 2

10 3

DNS Model

10 0 0

0.005

d32,DNS=2.33E-3 d32,LPP=2.49E-3

0.01

d

d

(a)

(b)

0.015

10 0 0

0.02

0.005

0.01

0.015

0.02

d

(c)

Fig. 12. Size distribution and Sauter mean of droplets formed in pulsed-jet atomization at times (time units Djet =U jet ): (a) t ¼ 8; d32;DNS ¼ 2:60 mm; d32;DNS ¼ 2:88 mm, (b) t ¼ 8:5; d32;DNS ¼ 2:55 mm; d32;DNS ¼ 2:66 mm, (c) t ¼ 9; d32;DNS ¼ 2:33 mm; d32;DNS ¼ 2:49 mm.

(a)

(b)

(c)

(d )

(e )

(f)

Fig. 13. The two-way RD-to-LPP and LPP-to-RD conversions for a droplet initially formed from a ligament breakup and later reabsorbed by the liquid jet, at times (time units Djet =U jet ): (a) t ¼ 2:0, (b) t ¼ 2:2, (c) t ¼ 2:4, (d) t ¼ 2:6, (e) t ¼ 2:8, (f) t ¼ 3:2.

Table 5 Parameters for the gas-assisted atomization problem.

Lz=H

Lx=8H Periodic

ql

(kg/m3)

c=0

Slip wall 1000

LPP

Inflow

qg

ll

lg

r

(Pa s)

(Pa s)

(N/m)

Ul (m/s)

Ug (m/s)

H (m)

dg (m)

(kg/m3) 10

103

104

0.2

8.3

0.1

0.01

0.0011

Separator plate Gas

H

interface

Ly=4H

H c=1

Liquid y

Slip wall z

Outflow

x Fig. 14. Problem description of gas-assisted atomization.

to simulations solving the Navier–Stokes equation without the LPP model, but it does not necessarily mean that every physical scale is fully resolved. In terms of turbulence, the numerical simulation can be considered as implicit LES instead of true ‘‘DNS’’. Nevertheless, the mesh resolution considered in this test seems to be sufficient to capture the macro-scale features of the atomization process and droplet dynamics, as shown later.

136

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Table 6 Important dimensionless parameters of the present simulation and the experiment by Descamps et al. (2008) of gas-assisted atomization.

ql =qg

ll =lg

Rel

Wel

M

Reg;d

Weg;d

Reg

Experiment

833

58.8

1000

1.45

69.1

960

5.68

1:69  104

Simulation

100

10

1000

1.45

68.9

913

11.0

8:30  103

General behavior The liquid jet atomization assisted by the fast gas stream is shown in Fig. 16 at different times. The gray and orange colors indicate the interface resolved by the VOF method and the small droplets represented by the LPP model, respectively. The vector plot on the background shows the flow velocity, where the vector magnitude is indicated by both length and color of the arrow. The flow velocity within the liquid is much smaller than in the gas, and the length of the arrows is too small to be appreciated. The simulation results are shown in Fig. 16 on a xy-plane. Due to the large velocity difference between the liquid and gas jets, a Kelvin–Helmholtz instability develops when the two streams meet at the end of the separator plate (Marmottant and Villermaux, 2004; Fuster et al., 2013). The interface wave appearing near the inlet, which is indicated by a purple arrow, is at first characterized by a small amplitude and by a shape close to a sinusoidal wave, see Fig. 16(a). As the wave propagates downstream, its amplitude grows and the wave crest rises into the fast gas stream, introducing a significant perturbation into the gas stream that strips small droplets away from the wave, see Fig. 16(b). Since the wave propagates much slower than the gas stream, it acts like an obstacle to the gas flow. As a consequence, the gas flow over the wave separates and forms a recirculation zone, indicated by a yellow arrow in Fig. 16. The recirculation zone formed at the downstream of the interface wave has also been observed experimentally and numerically by Jerome et al. (2013). The interaction between the wave and the gas stream becomes more complex as the wave further grows and propagates downstream, and the wave crest breaks up in a more pronounced way and generates a large amount of droplets, as shown in Fig. 16(b)–(d). The interface then starts to fold and liquid ligaments are formed, see Fig. 16(e). The gas stream is significantly influenced by the deformed interface and is observed to bend upward. A counterclockwise circulation is formed near the top of the gas stream, indicated by blue arrows Fig. 16(a), (b) and (e), (f). The clockwise circulation at the downstream of the wave grows in time and pushed downstream as the wave moves. These circulation regions are observed to have a significant impact on droplet formation (Jerome et al., 2013). The ligaments, stretched by the gas stream, finally break up into many droplets due to a Plateau-Rayleigh

0.14 0.12

Liquid Volume Fraction

The simulation starts at time t 0 ¼ 0 with no liquid inside the domain and it takes a long transient for the liquid jet to flow into the domain and reach an ‘‘equilibrium’’ state between the inflow and outflow of the liquid phase. The volume fraction of liquid within the whole domain as a function of time is plotted in Fig. 15. It can be seen that the amount of liquid in the domain increases gradually and finally reaches a plateau at about te ¼ 0:4 s. To save computational time only DNS is performed up to time t e . After reaching the ‘‘equilibrium’’ state, both DNS and the combined model are executed up to the final time tf ¼ 0:6 s. When running the combined model, droplets with an equivalent diameter smaller than four grid spacings are converted to Lagrangian particles if all conversion criteria are satisfied. The total number LPP droplets varies in time, the maximum value can reach up to 60,000. The simulations with a grid resolution h ¼ H=64, are conducted on 256 processors of the CINES-JADE machine: about 33 h of CPU time are required to reach time t e , and 12.5 more hours to time tf , either with DNS or the combined model.

0.1 0.08 0.06 0.04 0.02 0 0

0.1

0.2

0.3

0.4

0.5

t (s) Fig. 15. Time evolution of the volume fraction of liquid within the whole domain.

instability, as shown in Fig. 16(f). Finally, the atomized droplets are advected downstream and leave the computational domain.

Droplet dynamics and flying angle statistics Due to the strong interaction between the wavy interface and the gas stream, the formation and dynamics of the droplets are rather complex. The droplets are generated along the streamwise direction with different formation mechanisms. Near the inlet, droplets are mainly stripped off the interface by the gas stream, as highlighted by the purple arrow in Fig. 16(a)–(c). Further downstream, the droplet formation is mainly due to ligament break-up, or a combination of the two mechanisms, see Fig. 16(d)–(f). Droplets are not only generated at different streamwise positions by different mechanisms, but their dynamics changes remarkably as well, as shown in Fig. 17 by the trajectories of the droplets formed near the two streamwise positions x1 ¼ 2 H and x2 ¼ 4 H, at times t1 ¼ 0:46 s and t 2 ¼ 0:49 s, respectively. In the figure are also plotted a snapshot of the interface and the streamwise velocity component, on the plane z ¼ 0, to show the influence of the gas-interface interaction on the droplet dynamics. In particular in Fig. 17(a) at time t1 , droplets are formed at the wave crest when its amplitude is still relatively small and the interaction with the gas stream is not very important. As a result, the gas stream remains fairly horizontal and the flying angle of the droplet is generally small (a flying angle equal to 0 is along the x-direction). Except for a few large droplets that cross the whole gas stream, most of them stay in its interior as they are convected downstream. The flying angle statistics for the droplets in the sampling rectangular cuboid A of Fig. 17(a), and defined by the coordinates ranges 4:5 H 6 x  5:5 H; H 6 y 6 3 H, and 0 6 z 6 H, is shown in Fig. 18. The time interval for the sampling is 0:45 6 t 6 0:47 s, which is chosen to select the droplets that are formed close to the inlet. From Fig. 17(a), it can be observed that the droplet trajectories are quite straight after the early acceleration within the gas stream. Therefore, it is expected that the flying angles vary little while the droplets are moving within the sampling region. It can

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

(a) t=0.44s

(d) t=0.47s

(b) t=0.45s

(e) t=0.48s

(c) t=0.46s

(f) t=0.49s

137

Fig. 16. Numerical results with the combined model of the atomizing liquid jet for times in the range 0.44–0.49 s, viewed along the z-direction. The gray and orange colors indicate the interface resolved by VOF and the small droplets represented by LPP model, respectively. The vector plot on the background shows the flow velocity, with its magnitude indicated by both length and color of the arrows. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A

2D-view

3D-view

2D-view

3D-view

(a) t=0.46 s

B

(b) t=0.49 s

Fig. 17. 2D and 3D views of trajectories of the droplets formed at the top of the perturbation wave on the liquid–gas interface at times: (a) t ¼ 0:46 s, (b) t ¼ 0:49 s. The interface and streamwise fluid velocity are also shown in the background.

be clearly seen in Fig. 18(a) that most of the droplets are evenly distributed with a flying angle ranging from 30 to 30 degrees. The statistics include all droplets, hence represented by either the VOF method or the LPP model. On the other hand, in Fig. 17(b) at time t 2 the wave crest is near x2 , the interface has deformed seriously causing the gas stream to bend upward. The liquid ligaments are mainly aligned with the

stream, so when they break-up the droplets are ejected into the bended gas stream. As a result, most of the droplets are ejected with a positive angle and the corresponding flying angle distribution is shifted to the right as shown in Fig. 18(b). The sampling is done in the same volume A, but within the time interval 0:48 6 t 6 0:5 s, which mainly captures the droplets formed by the ligaments break-up, as shown in Fig. 17(b). Because of that,

138

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

0.5

0.5

t=0.45-0.47 t=4.5-4.7ss

0.4

N/sum (N)

N/sum (N)

0.4

0.3

0.2

0.1

0 -90

t=0.48-0.50 t=4.8-5.0ss

0.3

0.2

0.1

-60

-30

0

30

60

90

Flying Angle (deg)

(a)

0 -90

-60

-30

0

30

60

90

Flying Angle (deg)

(b)

Fig. 18. Numerical results of droplet flying angle distribution for the two time intervals: (a) 0:45 < t < 0:47 s, (b) 0:48 < t < 0:50 s. The sampling volume is defined by 4:5H < x < 5:5H; H < y < 3H, and 0 < z < H, and is also indicated by the letter A and the dotted yellow line in Fig. 17(a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the droplets are more aligned with the gas stream, their trajectories are less curved than for the droplets of Fig. 18(a) and their flying angle distribution in Fig. 18(b) is more peaked. The flying angle statistics for the droplets in the sampling rectangular cuboid B of Fig. 17(b), and defined by the coordinates ranges 4:5 H 6 x  5:5 H; 2:5 H 6 y 6 3:5 H and for the longer time interval 0:4 6 t 6 0:5 s, is shown in Fig. 19. The sampling volume B is similar to that of the experiment of Descamps et al. (2008) and the time period is long enough to cover the overall cycle of wave development and break up, as in Fig. 16. The simulation results agree fairly well with the experimental data, as most of the droplets are released with a positive flying angle and the peak of the distribution is at about 20 degrees. However, the experimental distribution profile is narrower than the numerical one.

Droplet Reynolds and Weber number statistics When the LPP model is applied, the relative velocity between the droplet and the surrounding fluid can be easily obtained, being ~  up , to compute the corresponding Reynolds and the term u Weber numbers. The probability density function (PDF) and mass-weighted probability density function (mPDF) of the LPP droplets as a function of the Reynolds and Weber numbers are shown in Fig. 20, for the sampling volume A of Fig. 17(a) and the two time intervals 0:45 6 t 6 0:47 s and 0:48 6 t 6 0:5 s. In the first time interval about 45% of the total number of droplets is represented by the LPP model, and about 64% in the second one. In general, the results for the two time intervals are rather similar, even if slightly more droplets are observed in the PDF of the second time interval at small Reynolds and Weber numbers. Both PDFs decrease with Rep and while most of the droplets are in a regime of very small Reynolds numbers, Rep < 1, still there is a significant amount of droplets in the flow at higher numbers, in the range 1 < Rep < 200 as shown in Fig. 20(a). If the statistics is made on mass rather than droplet number, the mass-weighted PDF is obtained. From Fig. 20(c) a large fraction of the total mass is in the droplets with 20 < Rep < 70. It should be reminded that the criterion dp < 4Dx is used for RD-to-LPP conversion. For this range of Reynolds numbers, it is expected that the droplets would not have been resolved accurately by VOF, had they remained RD. The LPP model is actually a better alternative for these poorly resolved droplets with a relatively large Reynolds number. Both PDF and mPDF decrease monotonically as a function of the Weber number, Wep , and for most droplets it is indeed quite small,

Wep < 3, in terms of both mass and particle number. The droplet deformation is then expected to be small, and the assumption of a spherical shape in the LPP model is indeed rather good. The LPP droplet Reynolds number is plotted in Fig. 21 as a function of the ratio between the droplet and cell sizes at two different time intervals. The scaling estimate of the relationship between the Reynolds number and the diameter of solid particles by Balachandar (2009) can be applied to the droplets here. When the droplet response time sp is larger than the largest ambient flow time scale, the relative velocity is mainly controlled by the ambient flow velocity and can be considered as independent of dp . In this limit it can be shown that Rep dp . On the other hand, when sp is smaller than the smallest ambient flow time scale, such as the Kolmogorov time scale, then the relative velocity can be estimated ~  up j  sp Du ~ =Dt, where by the Equilibrium Eulerian Approach as ju ~ =Dt is the acceleration of the ambient flow and is again indepenDu dent of dp (Ferry and Balachandar, 2001). Then it can be easily 3

shown that in this limit, Rep dp . It is observed from Fig. 21(a) and (b) that, the variation of the droplet Reynolds number with the droplet diameter lies in between these two limits, for both 0:45 < t < 0:47 s and 0:48 < t < 0:5 s, although more droplets with smaller Reynolds number are formed for the latter time period. Furthermore, the droplets with larger Reynolds number agree better with Rep dp ; while the droplets with smaller Reynolds number follow the Rep dp scaling. At last, Fig. 21 also shows that there are a significant amount of droplets with large Reynolds number and small diameter-to-cell-size ratio formed in atomization. It is expected the dynamics of these droplets would be incorrect if they remain to be resolved by VOF with low resolution (dp =h K 4). This confirms it is necessary to employ LPP model to represent these droplets. Comparison of droplet trajectories with experiment The trajectories of a sample of large droplets generated by the numerical simulation of the atomization process are shown in Fig. 22(a), while in the experiment the droplet trajectories are traced through high-speed cameras (Descamps et al., 2008). Due to the limited size of the computational domain, the comparison between numerical and experimental trajectories can be done only for a few droplets. The traced trajectories are characterized by a large positive flying angle, and they are likely to correspond to droplets with a big inertia. As a matter of fact, from the simulation results we can compute trajectories of large droplets, with

139

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

0.5

0.5

Exp

0.4

N/sum (N)

N/sum (N)

0.4

Model t=0.4-0.5 s

0.3

0.2

0.1

0.3

0.2

0.1

0 -30

0

30

60

0 -30

90

0

30

60

Flying Angle (deg)

Flying Angle (deg)

(a)

(b)

90

Fig. 19. Droplet flying angle distribution: (a) experimental data, (b) numerical results. The sampling region is volume B of Fig. 17(b) with 4:5H < x < 5:5H; 2:5 < y < 3:5H for the time interval 0:4 < t < 0:5 s. The sampling region is similar to the experimental one, with 4:5H < x < 5:5H and 2:5 < y < 4H (Descamps et al., 2008).

10 0

10 2

t=0.45-0.47s t=0.48-0.50s

t=0.45-0.47s t=0.48-0.50s

10 1

10 -1

10 0 PDF

PDF

10 -2 10 -1

10 -3 10 -2 10 -4

10 -5

10 -3

0

50

100

150

10 -4

200

0

0.5

1

Rep

(a) 10 -1

1.5

2

2.5

3

2.5

3

Wep

(b) 10 1

t=0.45-0.47s t=0.48-0.50s

10 -2

t=0.45-0.47s t=0.48-0.50s

mPDF

mPDF

10 0

10 -3

10 -4

10 -1

0

50

100

150

200

10 -2

0

0.5

1

1.5

Re p

We p

(c)

(d)

2

Fig. 20. Probability density function (PDF) and mass-weighted probability density function (mPDF) as a function of the Reynolds and Weber numbers, for the droplets in the sampling volume A of Fig. 17(a), with 4:5H < x < 5:5H; H < y < 3H, and 0 < z < H, in the two time intervals 0:45 < t < 0:47 s and 0:48 < t < 0:5 s.

dp ¼ 0:24—0:5 mm, which match quite well the experimental measurements. Both the numerical and experimental trajectories are rather straight. The computed trajectories for the two droplets with diameter dp equal to 0.24 and 0.31 mm, curve slightly as the droplets enters the faster gas stream and are subject to a sudden acceleration in the streamwise direction.

The corresponding Reynolds and Weber numbers along the droplets trajectory are shown in Fig. 22(b)–(c). For the largest droplet with diameter dp ¼ 0:49 mm, which belongs to the tail of the PDF of Fig. 20, the Reynolds number can be as high as 400, while the Weber number remains smaller than 5. These two numbers vary in a remarkable way as the droplets move downstream, though

140

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Rep~dp

t=0.45-0.47s

Rep~dp

t=0.48-0.50s

10 2

10 2

Re p

Re p

Rep~dp3 10 1

10 1

10 0 1.00

2.00

4.00

Rep~dp3

10 0 1.00

2.00

dp/h

dp/h

(a)

(b)

4.00

Fig. 21. Droplet Reynolds number as a function of the ratio between the droplet and cell sizes for the droplets in the sampling volume A of Fig. 17(a), with 4:5H < x < 5:5H; H < y < 3H, and 0 < z < H, in the two time intervals 0:45 < t < 0:47 s and 0:48 < t < 0:5 s.

0.04

y (m)

0.03

0.02

dp dp=0.24mm dp dp=0.28mm dp dp=0.31mm dp dp=0.49mm Exp

Gas Jet 0.01

Liquid Jet 0

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

x (m)

(a) 600

6

dp dp=0.24mm dp 5 dp=0.28mm dp dp=0.31mm dp 4 dp=0.49mm

Wep

Rep

ddp=0.24mm p p 500 ddp=0.28mm ddp=0.31mm p p 400 ddp=0.49mm 300

3

200

2

100

1

0

0

0.02

0.04

0.06

0.08

0

0

0.02

0.04

x (m)

x (m)

(b)

(c)

0.06

0.08

Fig. 22. Numerical results of a sample of large droplets as a function of the longitudinal x-direction: (a) projection of the trajectories on the xy-plane, (b) Reynolds number, and (c) Weber number. The numerical trajectories are also compared with the experimental measurements.

the trajectories are quite smooth. These fluctuations are related to the interaction with the gas turbulent flow. The turbulence of the gas flows is shown in background of Fig. 17. However, the inertia of these droplets is quite large, and the influence of turbulence on their trajectory is not profound. Droplet size distribution Finally, the size distribution for the droplets formed in the gas-assisted atomization is shown in Fig. 23. The droplet size distributions computed by the combined model and DNS are shown

in Fig. 23(a) and (b), respectively. The comparison of the probability density functions (PDF) of the droplet size for the combined model and DNS on the same size mesh is also shown in Fig. 23(d). The results of the combined model for the droplets larger than the grid size agree very well with the DNS result on the same size mesh. This indicates that the LPP coupling and conversion model do not introduce error on the droplet formation, which is mainly resolved by VOF. When the LPP model is introduced, a larger number of droplets smaller than the cell size are observed, see Fig. 23(a). The number of smaller-than-cell-size droplets is

141

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

DNS, h=H/64 4

10

10 3

10 3

counts

counts

Model, h=H/64 10 4

10 2

10 1

10 1

10 0

10 2

0

0.5

1

1.5

2

10 0

2.5

0

0.5

1

1.5

d (mm)

d (mm)

(a)

(b) 10 2

DNS, h=H/128

2

2.5

0.2

0.25

Model, h=H/64 DNS, h=H/64 DNS, h=H/128

10 5 10 1

PDF

counts

10 4

10 3

10 -1

10 2

10 1

10 0

0

0.5

1

1.5

2

2.5

10 -2

0

0.05

0.1

0.15

d (mm)

d/H

(c)

(d)

Fig. 23. Size distributions for droplets formed in atomization. (a) Counts of droplets for present model on the mesh h ¼ H=64 ¼ 0:156 mm, (b) Counts of droplets for DNS on the mesh h ¼ H=64 ¼ 0:156 mm, (c) Counts of droplets for DNS on the mesh h ¼ H=128 ¼ 0:078 mm, (d) Probability density function (PDF) (only droplets that are larger than the cell size are considered in PDF calculation.) The characteristic length scales for the original and finer meshes are k1 ¼ 0:166 mm and k2 ¼ 0:095 mm, respectively.

significantly lower for the DNS results without the LPP model as shown in Fig. 23(b). The discrepancy may be due to the erroneous dynamics of small droplets when they remain to be resolved by VOF. Indeed small fragments of VOF-represented liquid tend to be reconstructed and advected so erroneously that they do not move at all, an effect that can be easily demonstrated but is never discussed as the dynamics of such small fragments is considered largely erroneous anyway. As these small fragments stay where they are formed they are easily swept by and merged into the moving bulk liquid jet. If these small fragments are converted to LPP, then they can be advected away from the liquid jet and remain conserved. Since the size distributions for both the combined model and DNS are possibly erroneous for the size range dp < h. The droplets smaller than the cell size are thus excluded in the PDF statistics shown in Fig. 23(d). Furthermore, an exponential decay of droplet number with droplet diameter for dp 2 ½0:2; 1:1 mm is observed in both Fig. 23(a) and (b). The decay rate can be characterized by a length scale k1 , which is considered as the ‘‘average ligament size’’ in Marmottant and Villermaux (2004). A least square fit of the data of the present model for dp 2 ½0:2; 1:1 mm yields k1 ¼ 0:166 mm. The exponential decay with k1 ¼ 0:166 is shown to match very well with the DNS results on the same size mesh, see Fig. 23(b). It is noted that k1 is comparable to the cell size h ¼ 0:156 mm.

In order to see the effect of grid resolution a DNS simulation is conducted with a finer mesh, h ¼ H=128. The total number cells is about 67 million and it cost about 400 h on 2048 cores to reach similar physical time t ¼ 0:6 s. The detailed results of the finer-mesh DNS are not the focus of the present paper and discussion of which will be relegated to future work. The droplet size distribution obtained by DNS on a finer mesh is shown in Fig. 23(c). The comparison of the PDF of the droplet size with the two different meshes is shown in Fig. 23(d). With similar sampling time period (about 0.03 s), the droplet number for the finer mesh is an order of magnitude larger than that for the original mesh. It is also seen that the size distribution of small droplets dp < 0:2 mm is slightly better resolved, but not yet converged, when the finer mesh is used. Exponential decay of droplet number similar to the original mesh is observed for dp 2 ½0:1; 0:6 mm in the fine mesh results, though the decay is faster. The characteristic length scale for the exponential decay, k2 ¼ 0:095 mm, is 42% smaller than k1 ¼ 0:166 for the original mesh, indicating lack of convergence, so that yet unresolved small scales in the flow are important to droplet formation. To yield an accurate prediction of the droplet size distribution, including the expected peak at small scales and the exponential decay, the smaller scales must be resolved with sufficient grid resolution. It is likely that the computational cost for a fully-resolved DNS would be extremely high.

142

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143

Therefore, it may be unavoidable for future exploration of atomization in this parameter range to develop subgrid models accounting for the small-scale physics, such as interface instability, interface-turbulence interaction and droplet formation at the relevant scales.

acknowledge Grants x20142b6115 and x20142b7325 from GENCI. The authors acknowledge Dr. G. Tryggvason and Dr. S. Dabiri for communicating their code, which serves as the initial version of PARIS-Simulator. We also thank T.J. Arrufat, Dr. D. Fuster, L. Malan, and Dr. P. Yecko for their contributions to the development of PARIS-Simulator.

Conclusions In this paper, we introduce a multiscale simulation approach for atomization. While the interface between different phases is resolved by the Volume-of-Fluid (VOF) method, the small droplets formed in atomization are represented by the Lagrangian pointparticle (LPP) model. To accurately compute the dynamics of the LPP droplets that are larger than the grid spacing, a new model of the momentum coupling and the two-way conversion between the LPP droplets and the resolved flow is proposed. The combined model is validated by a comparison with experimental data and DNS simulations through a series of tests. A key aspect of the present momentum coupling model is to distribute the coupling force exerted back to the resolved flow over a length scale which is larger than the droplet diameter, say 5 to 10 times dp . By doing that the influence of an individual LPP droplet on the local gas flow is small and the gas flow velocity at the LPP droplet location is a close approximation of the undisturbed gas flow velocity that is required for the droplet force calculation. Attention is also paid on reconstructing the local flow field around the droplet during the two-way conversion between the resolved droplets and LPP droplets. In the RD-to-LPP conversion, the undisturbed flow field is rebuilt by interpolation from the flow properties away from the droplet while in the reverse LPP-to-RD conversion, excess momentum is added to compensate the development of the local disturbed flow field. The tests clearly show that the reconstruction of the local flow field is necessary to have a smooth transition in droplet velocity during the conversion. The combined VOF-LPP model is applied to simulate the gas-assisted atomization experiment by Descamps et al. (2008). The numerical results show complex droplet formation mechanisms. The dynamics, characterized for instance by the flying angle of droplets formed at various streamwise locations is observed to strongly depend on the observation window. The PDF and mass-weighted PDF of both Reynolds and Weber numbers of the LPP droplets (with diameter smaller than four cells) are shown. The Reynolds number of the LPP droplets changes in a wide range (up to several hundreds), while the Weber number is in general small. Resolving the droplets with high Reynolds number by VOF with less than four cells will yield catastrophic errors in the prediction of droplet dynamics, and it is a considerable improvement to represent them by the LPP model. Agreement with the experimental measurements is observed for both the the flying angle distribution and the trajectories of the droplets. The exponential decay of the PDF of droplet diameters is recovered but the rate of decay is not converged for the current grid sizes. The present model faithfully recovers the large scale features of atomization and the droplets at scales significantly larger than the grid size, though some small-scale physics is missed. Acknowledgements This project has been supported by the ANR MODEMI project (ANR-11-MONU-0011) program, and the FIRST (Fuel Injector Research for Sustainable Transport) project supported by the European Commission under the 7th Framework Programme under Grant Agreement No. 265848. The simulations of this paper are conducted on the our lab cluster, and both the IDRIS Turing machine and the CINES Jade machines for which we thankfully

References Agbaglah, G., Delaux, S., Fuster, D., Hoepffner, J., Josserand, C., Popinet, S., Ray, P., Scardovelli, R., Zaleski, S., 2011. Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method. C.R. Mec. 339, 194–207. Apte, S.V., Mahesh, K., Moin, P., Oefelein, J.C., 2003. Large-eddy simulation of swirling particle-laden flows in a coaxial-jet combustor. Int. J. Multiphase Flow 29, 1311–1331. Apte, S.V., Gorokhovski, M., Moin, P., 2003. LES of atomizing spray with stochastic modeling of secondary breakup. Int. J. Multiphase Flow 29, 1503–1522. Arrufat, T.J., Dabiri, S., Fuster, D., Ling, Y., Malan, L., Scardovelli, R., Tryggvason, G., Yecko, P., Zaleski, S. The PARIS-Simulator code. . Aulisa, E., Manservisi, S., Scardovelli, R., Zaleski, S., 2007. Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry. J. Comput. Phys. 225, 2301–2319. Balachandar, S., 2009. A scaling analysis for point particle approaches to turbulent multiphase flows. Int. J. Multiphase Flow 35, 801–810. Balachandar, S., Eaton, J.K., 2010. Turbulent dispersed multiphase flow. Annu. Rev. Fluid Mech. 42, 111–133. Basset, A.B., 1888. On the motion of a sphere in a viscous liquid. Phil. Trans. R. Soc. A 179, 43–63. Boeck, T., Li, J., López-Pagés, E., Yecko, P., Zaleski, S., 2007. Ligament formation in sheared liquid–gas layers. Theor. Comp. Fluid Dyn. 21, 59–76. Bornia, G., Cervone, A., Manservisi, S., Scardovelli, R., Zaleski, S., 2011. On the properties and limitations of the height function method in two-dimensional cartesian geometry. J. Comput. Phys. 230, 851–862. Chesnel, J., Menard, T., Reveillon, J., Demoulin, F.-X., 2011a. Subgrid analysis of liquid jet atomization. Atomization Spray 21, 41–67. Chesnel, J., Reveillon, J., Menard, T., Demoulin, F.-X., 2011b. Large eddy simulation of liquid jet atomization. Atomization Spray 21, 711–736. Chorin, A.J., 1968. Numerical solution of the Navier–Stokes equations. Math. Comput. 22, 745–762. Clift, R., Gauvin, W.H., 1970. The motion of particles in turbulent gas streams. Proc. Chemeca 1, 14–28. Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops, and Particles. Dover Publications. Climent, E., Magnaudet, J., 2006. Dynamics of a two-dimensional upflowing mixing layer seeded with bubbles: bubble dispersion and effect of two-way coupling. Phys. Fluids 18, 103304. Crowe, C.T., Sommerfield, M., Tsuji, Y., 1998. Multiphase Flows with Droplets and Particles. CRC Press. DeBar, R. (1974). Fundamentals of the KRAKEN code. Technical report UCIR-760, Lawrence Livermore National Laboratory, Livermore, California, USA. Descamps, M.N., Matas, J.-P., Cartellier, A. (2008). Gas–liquid atomisation: gas phase characteristics by PIV measurements and spatial evolution of the spray. In: Proceedings du 2nd colloque INCA, Initiative en Combustion Avancée. Di Felice, R., 1996. A relationship for the wall effect on the settling velocity of a sphere at any flow regime. Int. J. Multiphase Flow 22, 527–533. Eaton, J.K., 2009. Two-way coupled turbulence simulations of gas-particle flows using point-particle tracking. Int. J. Multiphase Flow 35, 792–800. Faxén, H., 1922. Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist. Ann. Der Phys. 373, 89–119. Ferry, J., Balachandar, S., 2001. A fast Eulerian method for disperse two-phase flow. Int. J. Multiphase Flow 27, 1199–1226. Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M., Williams, M.W., 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. J. Comput. Phys. 213, 141–173. Fuster, D., Bagué, A., Boeck, T., Le Moyne, L., Leboissetier, A., Popinet, S., Ray, P., Scardovelli, R., Zaleski, S., 2009. Simulation of primary atomization with an octree adaptive mesh refinement and VOF method. Int. J. Multiphase Flow 35, 550–565. Fuster, D., Matas, J.P., Marty, S., Popinet, S., Hoepffner, J., Cartellier, A., Zaleski, S., 2013. Instability regimes in the primary breakup instability regimes in the primary breakup instability regimes in the primary breakup region of planar coflowing sheets. J. Fluid Mech. 736, 150–176. Gatignol, R., 1983. The Faxén formulae for a rigid particle in an unsteady nonuniform Stokes flow. J. Mech. Theor. Appl. 1, 143–160. Ghia, U., Ghia, K.N., Shin, C.T., 1982. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J. Comput. Phys. 48, 387– 411. Herrmann, M., 2010. A parallel Eulerian interface tracking/lagrangian point particle multi-scale coupling procedure. J. Comput. Phys. 229, 745–759.

Y. Ling et al. / International Journal of Multiphase Flow 76 (2015) 122–143 Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201–225. Jerome, J.J.S., Marty, S., Matas, J.-P., Zaleski, S., Hoepffner, J., 2013. Vortices catapult droplets in atomization. Phys. Fluids 25, 112109. Josserand, C., Lemoyne, L., Troeger, R., Zaleski, S., 2005. Droplet impact on a dry surface: triggering the splash with a small obstacle. J. Fluid Mech. 524, 47–56. Labourasse, E., Lacanette, D., Toutant, A., Lubin, P., Vincent, S., Lebaigue, O., Caltagirone, J.-P., Sagaut, P., 2007. Towards large eddy simulation of isothermal two-phase flows: governing equations and a priori tests. Int. J. Multiphase Flow 33, 1–39. Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., Zanetti, G., 1994. Modelling merging and fragmentation in multiphase flows with SURFER. J. Comput. Phys. 113, 134–147. Lakehal, D., Labois, M., Narayanan, C., 2012. Advances in the large-eddy and interface simulation (leis) of interfacial multiphase flows in pipes. Prog. Comput. Fluid Dyn. 12, 153–163. Lakehal, D., Liovic, P., 2011. Turbulence structure and interaction with steep breaking waves. J. Fluid Mech. 674, 522–577. Larocque, J., Vincent, S., Lacanette, D., Lubin, P., Caltagirone, J.-P., 2010. Parametric study of LES subgrid terms in a turbulent phase separation flow. Int. J. Heat Fluid Flow 31, 536–544. Lasheras, J.C., Villermaux, E., Hopfinger, E.J., 1998. Break-up and atomization of a round water jet by a high-speed annular air jet. J. Fluid Mech. 357, 351–379. Le Chenadec, V., Pitsch, H., 2013. A monotonicity preserving sharp interface flow solver for high density ratio two-phase flows. J. Comput. Phys. 249, 185–203. Lebas, R., Menard, T., Beau, P.A., Berlemont, A., Demoulin, F.-X., 2009. Numerical simulation of primary break-up and atomization: Dns and modelling study. Int. J. Multiphase Flow 35, 247–260. Legendre, D., Magnaudet, J., 1997. A note on the lift force on a spherical bubble or drop in a low-Reynolds-number shear flow. Phys. Fluids 9, 3572–3574. Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng. 19, 59– 98. Ling, Y., Haselbacher, A., Balachandar, S., 2010. A numerical source of small-scale number-density fluctuations in Eulerian–Lagrangian simulations of multiphase flows. J. Comput. Phys. 229, 1828–1851. Ling, Y., Parmar, M., Balachandar, S., 2013. A scaling analysis of added-mass and history forces and their coupling in dispersed multiphase flows. Int. J. Multiphase Flow 57, 102–114. Lomholt, S., Maxey, M.R., 2003. Force-coupling method for particulate two-phase flow: Stokes flow. J. Comput. Phys. 184, 381–405. Lomholt, S., Stenum, B., Maxey, M.R., 2002. Experimental verification of the force coupling method for particulate flows. Int. J. Multiphase Flow 28, 225–246. Marmottant, P., Villermaux, E., 2004. On spray formation. J. Fluid Mech. 498, 73– 111. Matas, J.-P., Marty, S., Cartellier, A., 2011. Experimental and analytical study of the shear instability of a gas–liquid mixing layer. Phys. Fluids 23, 094112. Maxey, M.R., Patel, B.K., 2001. Localized force representations for particles sedimenting in Stokes flow. Int. J. Multiphase Flow 27, 1603–1626. Maxey, M.R., Patel, B.K., Chang, E.J., Wang, L.P., 1997. Simulations of dispersed turbulent multiphase flow. Fluid Dyn. Res. 20, 143–156. Maxey, M.R., Riley, J.J., 1983. Equation of motion for a small rigid sphere in a nonuniform flow. Phys. Fluids 26, 883–889.

143

Mei, R., Adrian, R.J., 1992. Flow past a sphere with an oscillation in the free-stream velocity and unsteady drag at finite Reynolds number. J. Fluid Mech. 237, 323– 341. Ménard, T., Tanguy, S., Berlemont, A., 2007. 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 33, 510–524. Pope, S.B., 2000. Turbulent Flows. Cambridge University Press. Popinet, S. The basilisk code. . Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. J. Comput. Phys. 190, 572–600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228, 5838–5866. Renardy, Y., Renardy, M., 2002. PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method. J. Comput. Phys. 183, 400–421. Rudman, M., 1997. Volume-tracking methods for interfacial flow calculations. Int. J. Numer. Methods Fluids 24, 671–691. Rudman, M., 1998. A volume-tracking method for incompressible multifluid flows with large density variations. Int. J. Numer. Methods Fluids 28, 357–378. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech. 31, 567–603. Schmidt, D.P., Rutland, C.J., 2000. A new droplet collision algorithm. J. Comput. Phys. 164, 62–80. Shinjo, J., Umemura, A., 2010. Simulation of liquid jet primary breakup: dynamics of ligament and droplet formation. Int. J. Multiphase Flow 36, 513–532. Sirignano, W.A., 2005. Volume averaging for the analysis of turbulent spray flows. Int. J. Multiphase Flow 31, 675–705. Smagorinsky, J., 1963. General circulation experiments with the primitive equations. Monthly Weather Rev. 91, 99–164. Sundaram, S., Collins, L.R., 1996. Numerical considerations in simulating a turbulent suspension of finite-volume particles. J. Comput. Phys. 124, 337–350. Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146–159. Tomar, G., Fuster, D., Zaleski, S., Popinet, S., 2010. Multiscale simulations of primary atomization. Comput. Fluids 39, 1864–1874. Toutant, A., Chandesris, M., Jamet, D., Lebaigue, O., 2009a. Jump conditions for filtered quantities at an under-resolved discontinuous interface. Part 1: Theoretical development. Int. J. Multiphase Flow 35, 1100–1118. Toutant, A., Chandesris, M., Jamet, D., Lebaigue, O., 2009b. Jump conditions for filtered quantities at an under-resolved discontinuous interface. Part 2: A priori tests. Int. J. Multiphase Flow 35, 1119–1129. Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S., Jan, Y.J., 2001. A front-tracking method for the computations of multiphase flow. J. Comput. Phys. 169, 708–759. Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct Numerical Simulations of Gas–Liquid Multiphase Flows. Cambridge University Press. Vincent, S., Balmigère, G., Caltagirone, J., Meillot, E., 2010. Eulerian–Lagrangian multiscale methods for solving scalar equations-application to incompressible two-phase flows. J. Comput. Phys. 229, 73–106. Vincent, S., Brändle de Motta, J.C., Sarthou, A., Estivalezes, J.-L., Simonin, O., Climent, E., 2014. A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows. J. Comput. Phys. 256, 582–614.