PHYSICAL REVIEW E 90, 033305 (2014)

Consistent simulation of droplet evaporation based on the phase-field multiphase lattice Boltzmann method Hesameddin Safari,1,* Mohammad Hassan Rahimian,1,† and Manfred Krafczyk2 1

2

Department of Mechanical Engineering, College of Engineering, University of Tehran, Tehran, Iran Institute for Computational Modeling in Civil Engineering, Technische Universit¨at Braunschweig, Braunschweig, Germany (Received 10 March 2014; published 10 September 2014) In the present article, we extend and generalize our previous article [H. Safari, M. H. Rahimian, and M. Krafczyk, Phys. Rev. E 88, 013304 (2013)] to include the gradient of the vapor concentration at the liquid-vapor interface as the driving force for vaporization allowing the evaporation from the phase interface to work for arbitrary temperatures. The lattice Boltzmann phase-field multiphase modeling approach with a suitable source term, accounting for the effect of the phase change on the velocity field, is used to solve the two-phase flow field. The modified convective Cahn-Hilliard equation is employed to reconstruct the dynamics of the interface topology. The coupling between the vapor concentration and temperature field at the interface is modeled by the well-known Clausius-Clapeyron correlation. Numerous validation tests including one-dimensional and two-dimensional cases are carried out to demonstrate the consistency of the presented model. Results show that the model is able to predict the flow features around and inside an evaporating droplet quantitatively in quiescent as well as convective environments. DOI: 10.1103/PhysRevE.90.033305

PACS number(s): 47.11.−j, 64.70.fm, 44.35.+c, 47.55.−t

I. INTRODUCTION

Droplet evaporation is a well-known phenomenon occurring in various natural processes and engineering applications. Liquid sprays and jets are amongst the most significant industrial applications driven by the dynamics of evaporating (or nonevaporating) droplets. A spray consists of droplets and ligaments interacting with its surrounding gas and exchanging momentum, heat, and mass with ambient gas. The investigation of a single evaporating droplet would be the first step to gain in-depth insight into the different physical phenomena taking place inside a spray. Thus, the problem of modeling droplet evaporation has always been an active research area since the pioneering work of Spalding [1] in the mid-1950s who proposed the well-known d 2 law. Nowadays, a great body of research can be found on the subject of evaporating droplets and sprays simulation, and different textbooks (e.g., [2,3]) as well as review articles (e.g., [4]) have been published. According to [2], the models of evaporating droplets can be divided into three major categories: (1) Models based on some simplifying assumptions for droplet temperature distribution, conductivity, and flow recirculation inside the droplet. This class contains different models such as the constant droplet temperature model, the infinite liquid conductivity model, the conduction limit model, and the effective conductivity model. (2) Vortex model for droplet heating which simulates the physical processes more directly. (3) The full solution of the coupled continuity, NavierStokes, and energy equations. The main difference between the models in the third class is how to simulate or resolve the interface between the liquid and vapor phase. The interface may be considered with a defined shape (say spherical). In this strategy, either the

* †

[email protected] [email protected]

1539-3755/2014/90(3)/033305(13)

droplet cannot be deformed (applicable at low Weber numbers) or the deformation is limited to a predefined shape (e.g., ellipsoid shape). Another approach is to directly simulate the interface, so it can be deformed freely in accordance with the flow field around the interface. Two different strategies can be employed for direct simulation of the phase interface: interface-tracking and interface-capturing methods. Tracking methods directly track the location of the interface and are based on a Lagrangian description for interface motion. On the other hand, capturing methods capture the interface between a number of the computational grids. The volume-of-fluid (VOF) method, the level set method, and phase-field modeling are amongst the most common capturing methods introduced so far. Below we present a brief overview on the progress achieved in direct numerical simulation of multiphase flows including phase change for incompressible fluids (such as droplet evaporation, boiling flows, etc.) in the last two decades. Haywood et al. [5] employed a moving mesh method for simulation of the deformation and evaporation of droplets at medium Reynolds numbers. They developed a twodimensional axisymmetric model and investigated the effect of transient droplet deformation on the evaporation rates and liquid phase heating. However, their model was limited to small deformations and resulted in inconsistent prolate shapes of the droplet as compared with experimental results. Juric and Tryggvason [6] presented a front-tracking method including phase changes for film boiling in which the mass, momentum, and energy equations are solved on a fixed grid; however, the interface is tracked by a moving front represented by linked lines or surface elements. In this method, the interface is smoothed with a finite thickness and the fluid properties are smeared out over a number of grid points to avoid numerical difficulties and enhance the stability. Thus, their method can be regarded as a numerically diffuse interface method. Despite the ability of interface-tracking methods in predicting the exact location of the phase interface, the application of these methods in situations with substantial topological changes and strong

033305-1

©2014 American Physical Society

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

interface deformations such as breakup and coalescence is difficult. Interface-capturing methods have been demonstrated to be applicable in modeling strongly deformed interfaces in a simple manner. Son and Dhir [7] used the level set method for computation of nucleate boiling on a horizontal surface. They also represented a diffuse description for the fluid properties across the interface using the so-called smeared-out Heaviside function. Welch and Wilson [8] adopted the VOF method to account for the phase-change process and employed it for modeling film boiling. Both classical level set and VOF methods have their own drawbacks: the level set method suffers from the lack of mass conservation while in the VOF method the accurate calculation of the interface curvature is difficult. In order to overcome these deficiencies, sometimes these two methods are combined. Tormar et al. [9] and Guo et al. [10] used this approach for simulation of film boiling. Schlottke and Weigand [11] presented a model for a three-dimensional evaporating and deforming droplet based on the VOF method. They showed that the correct calculation of both phases’ velocities at the interface is crucial for cases with high evaporative fluxes. They proposed an iterative procedure to accurately compute these velocities. Their model was validated against the experimental results for various Sherwood numbers. Besides the above models which can be considered as numerically diffuse interface methods, sharp interface-capturing methods have been extended for phase-change problems. Pioneering works in this field have been published by Tanguy et al. [12] and Gibou et al. [13] in which the level set method and the ghost fluid method are coupled to successfully capture the interface motion as well as to accurately impose the correct jump conditions at the sharp interface. Tanguy et al. [12] simulated the evaporation of a two-dimensional axisymmetric droplet and obtained a consistent Stefan flow around the evaporating droplet. However, their results have not been compared quantitatively with any experimental data. Gibou et al. [13] used various one-dimensional test cases to validate their solutions. They also simulated the two-dimensional film boiling. The VOF method, level set method, and front tracking of Tryggvason et al. are often called numerical-diffuse-interface methods. Another approach for the direct numerical simulation of multiphase flows is based on the so-called physical-diffuseinterface methods. In this approach the liquid-vapor interface is postulated as a region with strong but smooth transition of physical properties between the bulk values. The second gradient theory [14] and the phase-field method [15,16] are well-known physical-diffuse-interface approaches employed in computational fluid dynamics calculations. Jamet et al. [14] presented a method based on the second gradient theory for modeling two-phase flows with phase change. In this method the transition at the interfacial zone is obtained by relating the internal energy of the fluid to the density gradient. One interesting feature of this method is its thermodynamically consistent behavior, so that the phase-change process as well as the displacement and topology change of the interface can be determined by solving a single set of partial differential equations without any special treatment. The major drawback of this method is that the interface thickness and the surface tension are functions of the mass flux and the temperature

gradient across the interface. Thus, limited heat and mass transfer rates across the interface can be simulated. In the phase-field method [16], the convective Cahn-Hilliard equation is used to track the order parameter which determines the distribution of each phase. For multiphase flows with a density contrast between phases, the volume fraction can be selected as the order parameter (similar to the volume fraction variable in the volume-of-fluid method). All of the aforementioned models are incorporated into Navier-Stokes equation solvers. Recently, the lattice Boltzmann method (LBM) has attracted growing attention as an alternative approach for computational fluid dynamics [17–19]. One of the most significant advantages of LBM over classical Navier-Stokes solvers is its versatility in computation of the isothermal multicomponent and multiphase flows. Since LBM is based on the mesoscopic kinetic equations, interfacial physics of multiphase flows can be modeled from a more fundamental basis. Different lattice Boltzmann (LB) multiphase models have been proposed to date [20]. Most of them can be classified as physical-diffuse-interface methods. In spite of two decades of active research on the LB simulation for multiphase flows, a comprehensive model capable of simulating phase change quantitatively in a consistent manner (especially droplet evaporation) is still lacking. Some limited works can be found in literature concerning phase-change modeling using multiphase LB models [21–26]. The majority of them such as Refs. [22,24,25] dealt with boiling flows. More details on the drawbacks of available LB models in the scope of phase change have been discussed in our previous publication [26]. Moreover, some interesting works on the numerical simulation of the multiphase microflows, droplet breakup, and coalescence (applicable to spray modeling) using LBM have been published in recent years [27–29]. The present work aims to generalize our previous model [26] in a way that evaporation can take place at every temperature. In our previous work [26], the liquid temperature was assumed to be the saturation temperature while the gas phase was superheated. Thus, the liquid temperature remained unchanged during the vaporization and the temperature gradient at the interface considered as the driving force for vaporization process. Now, in the current article, vaporization is assumed to occur by the diffusion of vapor molecules from the liquid-vapor interface to the ambient gas. If the vapor concentration of the surrounding gas is less than that of the interface, there exists a driving force for vaporization. So the evaporation from the phase interface may happen at every temperature. From thermodynamics, one can easily find that the vapor concentration (or vapor mass fraction) at the interface is highly dependent on the interface temperature. Hence, in addition to the flow field, the energy equation and vapor species transport equation has to be solved. We show that for high vaporization rates, employing passive scalar lattice Boltzmann solver for the temperature and vapor mass fraction transport equations leads to overprediction of the calculated evaporating flux. Therefore, the transport equations of the energy and vapor species are solved using a finite-difference scheme. The coupling between the temperature and vapor mass fraction at the interface is modeled via the famous Clausius-Clapeyron equation. Furthermore, the computation of the two-phase flow field is carried out by the multiple-relaxation-time (MRT)

033305-2

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . .

PHYSICAL REVIEW E 90, 033305 (2014)

collision operator in order to stabilize the simulations for the higher Reynolds numbers.

the macroscopic velocity (volume averaged velocity), and cs2 eq is a constant representing the lattice speed of sound. fα , the equilibrium distribution function, is given by [30]   ea · u (ea · u)2 (u · u) , (4) − fαeq = wα ρ 1 + 2 + cs 2cs4 2cs2

II. GOVERNING EQUATIONS AND MATHEMATICAL MODEL A. Computation of the flow field using phase-field LBM

Our model is based on the phase-field lattice Boltzmann approach of Lee [30,31]. In this method the interface between the two immiscible and incompressible phases is governed by the convective Cahn-Hilliard equation. The liquid phase composition denoted by C = ρρ˜ll is considered as the order parameter, where ρ˜l and ρl are local and bulk densities of the liquid phase. The original Cahn-Hilliard equation is derived without the consideration of phase-change phenomena within the interface. Including the contribution of evaporative mass flux across the interface, the modified Cahn-Hilliard equation is given as [26] ˙  m ∂C + ∇ · (uC) = ∇ · (M∇μ) − , ∂t ρl

(1)

˙  represents the volumetric mass flux due to phase where m ˙  depends on the gradient of change. In our previous work, m the temperature; however, in this article the volumetric mass flux depends on the gradient of the vapor mass fraction. M is the mobility factor and μ denotes the chemical potential calculated as in Refs. [30,32]. The divergence of the velocity field within the interface is consequently obtained as [26]   1 1  ˙ ∇·u=m . (2) − ρg ρl

in which wα is the weight factor determined from lattice structure. The intermolecular force (F) can be written in the potential form to eliminate the parasitic currents: F = ∇ρcs2 − ∇p − C∇μ,

where p is the hydrodynamic pressure to enforce incompressibility. The DBE of Eq. (3) is transformed into the DBE for the pressure evolution and momentum equations by introducing a new distribution function (gα ) as   gα = fα cs2 + p − ρcs2 α (0), (6) eq

with α (u) = fρα . Another distribution function, hα = (C/ρ)fα , is also defined for interface capturing and recovering the modified convective Cahn-Hilliard equation. As shown in [26], the DBEs for gα and hα can be derived as   ∂gα + eα · ∇ga = − gα − gαeq + (eα − u) ∂t  · ∇ρcs2 [α − α (0)] − C∇μα   1 1 ˙  α (0), − (7) + ρcs2 m ρg ρl  ∂hα 1 + eα · ∇ha = − hα − heq α + (eα − u) ∂t λ   C · ∇C − 2 (∇p + C∇μ) α ρcs   ˙  m + M∇ 2 μ − (8) α , ρl

In order to compute the flow field, the discrete Boltzmann equation (DBE) of a nonideal fluid with a generalized collision operator  is written as Dfα ∂fα = + eα · ∇fα Dt ∂t   1 = − fα − fαeq + 2 (eα − u) · Ffαeq , ρcs

(3)

where fα is the particle distribution function, eα is the αdirection microscopic particle velocity, ρ is the density, u is

(5)

where λ is the relaxation time. It should be noted that in [26] we used a single-relaxation-time (SRT) operator for both DBEs. In the current study, the DBE for the pressure evolution and velocity is obtained by employing a generalized collision operator instead of a SRT operator. Applying the trapezoidal integration along characteristics over the time step δt and employing the modified particle distribution function as in Ref. [30], the lattice Boltzmann equations can be written as



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



   



1 1 δt 1 δt 1 ˙  ˙  α (0)

α (0)

+ ρcs2 m − + ρcs2 m − , 2 ρg ρl 2 ρg ρl (x,t) (x+eα δt,t+δt) 

1 ¯

hα − h¯ eq α (x,t) τ + 0.5  

C M M M + δt(eα − u) · ∇ C − 2 (∇ p + C∇ μ) α

ρcs (x,t)



    



˙ ˙ δt m δt m 2 2

α

α

M∇ μ − M∇ μ − + + , 2 ρl 2 ρ l (x,t) (x+eα δt,t+δt)

(9)

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

033305-3

(10)

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

where S = δt, I is the identity matrix, τ = λ/δt is nondimensional relaxation time, and ∇ M denotes the mixed difference approximation. A detailed explanation about gradient operators calculation is found in [31]. We utilize the MRT collision operator presented by Fakhari and Lee [33] in order to stabilize the computations at lower viscosities. For brevity we avoid mentioning the details of the collision operator (for details see Ref. [33]). The macroscopic properties can be computed using the zero- and first-order moments of the modified distribution functions as  C= h¯ α , (11)

fraction gradient at the interface can be seen as the driving force for evaporation. Applying the vapor species conservation at the interface gives the local vaporizing mass flow rate per unit surface as [11,12] ˙  = m

p=

 α

g¯ α +

δt u · ∇ρcs2 . 2

˙  = m (12)

(13)

B. Vapor mass fraction and temperature field

In order to obtain the vapor mass fraction field, the vapor species transport equation should be solved. The continuity equation of the vapor species inside an incompressible gas flow in the absence of chemical reactions is expressed as ∂ (ρg Yv ) + ∇ · (ρg Yv u) = ∇ · (ρg Dvg ∇Yv ), (14) ∂t where Yv is the vapor mass fraction, ρg represents the bulk density of the gas phase, and Dvg denotes the binary mass diffusivity of the vapor inside the gas phase. Since in our study both liquid and gas phases are assumed to be incompressible, Eq. (14) can be rewritten as ∂Yv + u · ∇Yv = Dvg ∇ 2 Yv . (15) ∂t Likewise, by neglecting the viscous dissipation, the energy equation for both phases in the incompressible limit with constant properties is given by ˙  hfg m ∂T + u · ∇T = α∇ 2 T − , ∂t ρcp

with

    ea · u ea · u sαeq = wα Yv 1 + 2 and qαeq = wα T 1 + 2 . cs cs (21)

Applying the multiscale Chapman-Enskog expansion of Eq. (20) shows that the following convection-diffusion equations would be recovered with second-order accuracy in time and space [36], Dvg ∂ ∂Yv + ∇ · (uYv ) = Dvg ∇ 2 Yv + 2 ∇ · (Yv u), ∂t cs ∂t

(16)

˙  hfg m ∂T α ∂ + ∇ · (uT ) = α∇ 2 T − + 2 ∇ · (T u), ∂t ρcp cs ∂t (22)

(17)

in which the subscripts l and g denote the bulk values of the liquid and the gas phase, respectively. The last term on the right-hand side of Eq. (16) accounts for a part of energy which is absorbed by the interface to accomplish the vaporization with hfg being the latent heat of evaporation. Previously [26], we assumed that the liquid temperature remains constant at the boiling temperature during vaporization and the existence of temperature gradient across the interface leads to phase change. Now we remove this assumption by considering the vapor mass fraction field. In this concept the evaporation is considered as a continuous diffusion of vapor from the liquidvapor interface to the surrounding gas. Thus, the vapor mass

(19)

 ∂sα 1 + eα · ∇sα = − sα − sαeq , ∂t λv  ˙  hfg m ∂qα 1  + eα · ∇qα = − qα − qαeq − wα , (20) ∂t λT ρcp

α = Cαl + (1 − C)αg , (1 − C)ρg Cρl cp,l + cp,g , ρ ρ

ρg Dvg ∇Yv · ∇C. 1 − Yv

The advection-diffusion equations [Eqs. (15) and (16)] can be solved either in a lattice Boltzmann framework or by applying finite-difference schemes. Using a passive scalar approach, another set of distribution functions is employed to recover the advection-diffusion equation. One can introduce two sets of distribution functions, namely, qα and sα , to recover the energy and vapor mass fraction transport equations, respectively. In general, the DBEs for recovering Eqs. (15) and (16) can be expressed respectively as [35–37]

where α, the thermal diffusivity, and cp , the specific heat capacity at constant pressure, are calculated by

cp =

(18)

where nˆ is the unit vector normal to the interface. As pointed ˙  with the local interface content out in [26], by multiplying m (|∇C|) [34] and having nˆ = ∇C/|∇C|, the volumetric mass flow rate of evaporation can be computed by

α

δt 1  eα g¯ α − C∇μ, ρu = 2 cs α 2

ρg Dvg ˆ ∇Yv · n, 1 − Yv

where the binary mass diffusivity and thermal diffusivity are related to the relaxation times by Dvg = λv cs2 and α = λT cs2 . The last terms on the right-hand sides are considered as the unwanted terms recovered by the multiscale expansion. At steady-state conditions, the unwanted terms would eventually go to zero. For transient incompressible flow, the unwanted 2 terms are of order O( uc2 ), which can be neglected in low s Mach number limit ( cus  1) [37]. In the case of neglecting the unwanted terms, comparing Eq. (22) with Eqs. (15) and (16) reveals that these two sets of equations will be identical if the flow field is incompressible (i.e., ∇ · u = 0). Since in the present model, bulk phases are considered to be incompressible, utilizing a LB scheme for solving the temperature and vapor mass fraction seems to be appropriate.

033305-4

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . .

However, we have a diffuse interface between two phases across which the phase change takes place. As shown in previous sections, we impose a virtual compressibility by assigning Eq. (2) and incorporate it into the flow field solver. It has been demonstrated in [26] that this approach generates a consistent Stefan flow around the phase interface in a simple manner without the necessity of imposing separate mass, momentum, and energy jump conditions at the interface. On the other hand, as the divergence of the velocity has a nonzero value (∇ · u = 0) within an evaporating front, Eq. (22) and Eqs. (15) and (16) would not be identical to each other inside a diffuse interface area which undergoes vaporization. So using the standard LB scheme for solving the advection-diffusion equation will introduce some errors in evaluating the temperature and vapor mass fraction within the evaporating interface. The solutions of Eq. (22) deviate more from those of Eqs. (15) and (16) as the evaporation rate increases. In order to minimize the effects of virtual compressibility on the calculation of the temperature and vapor mass fraction fields, a finite-difference scheme for Eqs. (15) and (16) is developed rather than a lattice Boltzmann scheme. More comparisons about this issue will be presented in the following section. In order to have second-order accuracy in both time and space, Eqs. (15) and (16) are discretized employing the Crank-Nicholson scheme as 1 Yvn+1 − Yvn 1 = Dvg ∇ 2 Yvn+1 + Dvg ∇ 2 Y,n − un · ∇Yvn , δt 2 2 (23) 1 T n+1 − T n 1 = α∇ 2 T n+1 + α∇ 2 T n δt 2 2 ˙ n hfg m −un · ∇T n − . ρcp

(24)

Equation (23) is only solved for the gas phase, while Eq. (24) is considered for the whole computational domain. So the vapor mass fraction on the liquid-vapor interface should be known in order to be applied as a boundary condition of Eq. (23). The Clausius-Clapeyron correlation is utilized to estimate the amount of vapor mass fraction at the interface [38]:    hfg WVM 1 1 pv = pamb exp − , (25) − R Ti TB where pv is the saturated vapor pressure on the interface, pamb is the ambient pressure, R is the universal perfect gas constant, WVM is the vapor molecular weight, Ti is the interface temperature, and TB is the liquid boiling temperature corresponding to the ambient pressure. It is known that the vapor mole fraction is related to the vapor partial pressure as [38] pv . (26) Xv = pamb In a gas mixture, the vapor mass fraction can be simply related to the vapor mole fraction by Xv WVM . Yv = M i Xi Wi

PHYSICAL REVIEW E 90, 033305 (2014)

FIG. 1. Schematic of the one-dimensional evaporation in an open end container. III. RESULTS AND DISCUSSION

In vaporization context, we encounter a physical phenomenon called Stefan flow. Stefan flow is a kind of convective flow which is generated internally and can be present even if there is not any externally imposed flow. This flow is induced as a result of the gasification of a liquid. A reliable and consistent model for the evaporation process should correctly account for the mass transfer as well as heat transfer across the interface. So the validity of the model in the above-mentioned aspects should be examined. In order to validate the computed evaporative mass flux, a one-dimensional simulation is conducted. We assume an open container partially filled with liquid, say water, is located in a surrounding gas, say air, as illustrated in Fig. 1. If the vapor concentration in the surrounding gas is less than the vapor concentration at the liquid-gas interface, the vapor will diffuse from the interface to the open end of the container and evaporation occurs. The liquid temperature is assumed to be constant, so the vapor concentration at the interface will remain unchanged during vaporization. The height of L and the environment vapor concentration are kept constant, too. Furthermore, if the surrounding gas is insoluble in the liquid phase, there will be no net transport of the surrounding gas in the container. From the continuity equation at the steady-state condition [d(ρu)/dx = 0], we have ˙  (x) = ρu = m ˙ v + m ˙ gas = constant. m

(28)

˙ gas = 0, the conservation of the vapor mass in the Since m gaseous phase can be written as ˙ v − ρDvg ˙ v = Yv m m

dYv . dx

(29)

Integrating Eq. (29) (with the assumption of ρDvg = constant) and applying two boundary conditions as Yv (x = 0) = Yv,s ,

Yv (x = L) = Yv,L ,

(30)

one can easily find the evaporative mass flux in the form of ˙ v = m

ρDvg ln(1 + BM ), L

(31)

in which the mass transfer number is defined as

(27)

033305-5

BM =

Yv,s − Yv,L . 1 − Yv,s

(32)

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

2.4

4 LB Solver for Transport Equations

2.2

Finite Difference for Transport Equations

3.5

Nondimensional Evaporating Flux

Nondimensional Evaporating Flux

2 1.8 1.6 1.4 1.2 Analytical Solution

1

Simulation (Δx / L = 0.005)

0.8

Simulation (Δx / L = 0.010)

0.6

Simulation (Δ x / L = 0.020)

0.4

3 2.5 2 1.5 1

Simulation (Δx / L = 0.050)

0.5

0.2 0 0

2

4 6 Mass Transfer Number

8

0

10

0

1

2 3 Mass Transfer Number

4

5

FIG. 2. (Color online) Comparison between analytical solution and simulation results of nondimensional evaporative mass flux for different grid sizes.

FIG. 3. (Color online) Comparison between different strategies in solving advection-diffusion transport equations during evaporation.

A D1Q3 lattice Boltzmann model is developed for simulating the above constant temperature evaporation problem. The liquid and gas phases are selected as water and air, respectively. The vapor mass fraction at the interface is set between 0.05 and 0.90. Also the value of Yv,L = 0 is kept constant during the simulation. Figure 2 compares the analytical solution with simulation results of the nondimensional evaporative mass ˙ v /(ρDvg /L)] for different mass transfer numbers and flux [m various grid sizes ( x ). As shown, the simulation results L match perfectly with the analytical solution as the grid is refined, which confirms that the model predicts the generated Stefan flow consistently and satisfactorily accounts for the mass transfer across the phase interface. As stated in the preceding section, in the case of phasechange occurrence across the interface, employing the LB solver for the heat and species advection-diffusion transport equations may cause some errors in evaluating the evaporative flux. Since the evaporative mass flux depends on the gradient of the vapor mass fraction across the interface, accurate estimation of the vapor mass fraction field is crucial for consistent calculation of the vaporization rate. As discussed before, if the LB solver is used for solving the advectiondiffusion transport equation of vapor mass fraction, two terms would affect the accuracy of the results. (1) Yv ∇ · u, which leads to the contribution of the virtual compressibility in evaluation of the vapor mass fraction across the interface and D (2) cvg2 ∂t∂ ∇ · (Yv u), which is considered as the unwanted term s recovered from the Chapman-Enskog expansion. Figure 3 demonstrates the differences between the results of the LB solver and the finite-difference (FD) solver for the test case of Fig. 1. As shown, the LB solver overpredicts the calculated evaporating flux for mass transfer numbers higher than 0.5. In this test case, due to the assumption of the constant length between the interface and the container opening as well as the constant vapor concentration at the

interface, evaporation occurs at steady-state rate. For steadystate simulations, the unwanted term would be zero. Thus, the discrepancies between the results of the LB solver and FD solver in Fig. 3 are owing to the nonzero value of the divergence of the velocity within the interface. In order to quantify the source of these differences, the average value of Yv ∇ · u through the interfacial region is nondimensionalized 2 using the mass diffusion time scale ( DLvg ). Figure 4 shows the dimensionless value of Yv ∇ · u against the mass transfer number. It can be seen that for mass transfer numbers below 0.25, the value of Yv ∇ · u is of order less than O(10−1 ) and the difference between the LB and FD solver may be 12 11

Nondimensional Value of Y v∇⋅u

10 9 8 7 6 5 4 3 2 1 0 0

1

2

3

4 5 6 7 Mass Transfer Number

8

9

10

FIG. 4. (Color online) The average nondimensional value of Yv ∇ · u inside the phase interface for various mass transfer numbers.

033305-6

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . .

(a) 0.4 Dvg ∂ ∇⋅(Yv u) 0.3 c2s ∂t 0.2

Yv∇⋅u

0.1 0 -0.1

0

0.0025 0.005 0.0075 Nondimensional Time

0.01

0.0005 (b) 0.00025 Dvg ∂ ∇⋅(Yv u) c2s ∂t 0 Yv∇⋅u

-0.00025

-0.0005 0.25

0.4

0.55 0.7 0.85 Nondimensional Time Dvg ∂ 2 ∂t

1

∇·(Yv u)

FIG. 5. (Color online) The average ratio of cs Yv ∇·u across the interface at (a) the initial stages of the simulation and (b) the developed phase of evaporation.

negligible. However, increasing the mass transfer number will substantially increase the divergence of the velocity field within the interfacial region. Therefore, for problems with high mass transfer numbers, the finite-difference scheme should be employed in order to achieve consistent simulations. For transient cases, in addition to Yv ∇ · u, the unwanted D term [ cvg2 ∂t∂ ∇ · (Yv u)] may further affect the deviation between s the LB and the FD solver. To perform a deeper analysis on this matter, the test case of Fig. 1 is adapted so that the interface can recede as a result of evaporation. So the evaporation rate would not be constant during the vaporization process. In order to examine the order of magnitude of the unwanted term relative to the value of Yv ∇ · u, the average ratio of Dvg ∂ cs2 ∂t

∇·(Yv u)

Yv ∇·u

the relative magnitude of the unwanted term for simulation times in which the evaporation process reaches the quasisteady state or developed phase. In this phase, the interface has been accelerated and continues receding with a nearly constant (but not constant) rate. Since the length between the interface and the open end of the container increases with time, the interface velocity has a decreasing trend. However, the rate of decreasing the interface velocity is so small that the temporal variations in the flow field can be ignored as shown in Fig. 5(b). The overshoots observed in the Fig. 5(b) are due to the displacement of the numerical interface location [defined by the cutoff value for the composition (C) as mentioned in [26]] between two neighboring computational cells. These overshoots are of order O(10−4 ) and, hence, can be neglected. Thus, it can be concluded that except for the initial time steps, the unwanted term does not change the discrepancies observed between the LB and FD solvers. In order to evaluate the capability of the model in predicting the combined heat and mass transfer across the interface, twodimensional static water droplet evaporation in the quiescent air is considered. If the relative humidity of the surrounding air is less than 100%, liquid water will evaporate from the interface and the temperature of the interface and its neighboring air-vapor mixture will decrease consequently. The liquid temperature decreases until the wet-bulb temperature is reached. At this steady-state condition, the coupled heat and mass transfer continues with constant rates across the interface and further cooling of the liquid phase will not occur. A two-dimensional D2Q9 lattice Boltzmann model is developed for the simulation of the nonisothermal water drop evaporation in the air. The geometry of the computational domain is shown in Fig. 6. A quadrant of a two-dimensional droplet is initialized at the corner of the domain. Since no external flow is imposed and only Stefan flow will be generated internally from the liquid-gas interface, the initial symmetry with respect to the axis of coordinates is not violated. Thus,

Symmetric boundary condition

0.5

PHYSICAL REVIEW E 90, 033305 (2014)

within the interface is plotted in Fig. 5 against the

nondimensional time (t ∗ = Lvg2 ) for the mass transfer number of BM = 9.0. As observed in Fig. 5(a), at initial time steps, the magnitude of the unwanted term is comparable to Yv ∇ · u, but it is rapidly damped. The high levels of the unwanted term at the initial stages of the simulation are mainly due to the fact that the interface starts accelerating from the stationary state, and the velocity magnitude as well as the divergence of the velocity field within the interface is changed rapidly. For t ∗ > 0.005, the magnitude of the unwanted term seems to be negligible compared to Yv ∇ · u. Figure 5(b) demonstrates

∂p h = 0 , T =T g , b , Yv =Y v , b ∂n

Gas phase

Yv =Y v,s

(determined by Clausius-Clapeyron equation)

D t

Liquid phase y Symmetric boundary condition x

FIG. 6. (Color online) Schematic of the computational domain for nonisothermal static droplet evaporation.

033305-7

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

TABLE I. Physical properties of the liquid and gas phase.

Gas (air) Liquid (water)

ρ (kg m−3 )

μ (kg m−1 s−1 )

λ (W m−1 K−1 )

cp (kJ kg−1 K−1 )

Molecular weight (kg kmol−1 )

hfg (kJ kg−1 )

1.165 997

1.85 × 10−5 8.55 × 10−4

0.026 0.613

1.005 4.18

28.85 18

2400

for the reduction of the computational cost, a quadrant of the whole physical domain can be solved instead. The length of the domain in both the x and the y direction is L = 5R0 , where R0 is the initial droplet radius. In our simulations, we assume the constant physical properties and the related values for both phases at 25 ◦ C, which are listed in Table I. The temperature and humidity of the surrounding air is set to the desired value by applying the constant temperature of Tg,b and the constant vapor mass fraction of Yv,b as boundary conditions on the outflow boundaries. The vapor mass fraction on the outflow boundaries is chosen in a way that a desired relative humidity is obtained on the boundaries. For a given dry-bulb temperature and relative humidity, the mv value of the humidity ratio (ω = m ) can be determined from a the psychrometric charts, where mv and ma are the mass of water vapor and dry air, respectively. One can easily relate the ω humidity ratio to the vapor mass fraction as Yv = 1+ω . For ∂ph the flow field, the boundary condition of ∂n = 0 is applied on the outflow boundaries and second-order extrapolation is employed for evaluation of the unknown distribution functions (h¯ α ) of the phase composition (C). More details about the implementation of the boundary conditions are found in our previous article [26]. As the first test case, a droplet with the initial radius of R0 = 50 μm and the temperature of 25 ◦ C is placed inside the stagnant air with the temperature of 25 ◦ C and the relative humidity of 40%. Since the thermal and mass diffusivities of the gas phase are larger than the liquid thermal diffusivity, droplet evaporation is controlled by heat transfer inside the liquid droplet. So we use the liquid heat diffusion R2 time scale ( α0l ) to nondimensionalize the simulation time (t) by t ∗ = αRl2t . The following figures represent the simulation 0 results at t ∗ = 1.7 in which the droplet reaches its equilibrium temperature (wet-bulb temperature) and the heat and mass are transferred with a steady rate across the interface. The velocity vectors around the interface and the droplet interface shape are shown in Fig. 7. It can be seen that the droplet shape remains circular during the vaporization which is consistent with theory. Moreover, the radial direction of the velocity vectors normal to the interface shows the ability of the model to internally generate a consistent Stefan flow around the interface of a vaporizing drop as expected. Figure 8 represents the temperature field inside the domain and the vapor mass fraction field around the interface. It can be observed that no artificial heating or cooling of the drop interface occurs and the temperature and vapor mass fraction are coupled successfully at the interface using Clausius-Clapeyron relation [Eq. (25)]. In order to obtain better insight into the evolution of the temperature field inside the droplet and the surrounding air, the temperature profiles for different nondimensional

Dvg (s−1 )

σ (N m−1 )

2.6 × 10−5 0.071

times are depicted in Fig. 9. As shown, the temperature in the vicinity of the interface decreases rapidly at the initial stages of the evaporation, while more time is needed for the inside temperature of the drop to be reduced (because of the = 0.0065). low thermal diffusivity of the liquid phase, ααwater air As the evaporation proceeds, heat flows outwards from the interface causing the temperature of the interface and inside the droplet to decrease continuously until the equilibrium temperature (wet-bulb temperature) is reached. At this state, the temperature inside the drop is uniform and lower than the ambient temperature. The amount of heat transferred from the surrounding to the interface is equal to the amount of heat required for vaporization and the drop temperature remains constant. To examine the predicted wet-bulb temperature, different simulations are performed on the above computational domain. First, for the definite relative humidity of 40%, the far-field air temperature (considered as the dry-bulb temperature) and the initial droplet temperature are varied between 288.15 and 318.15 K. Each simulation proceeds until the droplet temperature becomes uniform and stationary. Figure 10 compares the computed wet-bulb temperatures with the values obtained from the psychrometric charts. Excellent agreement can be observed in the chart. Now we set the dry-bulb temperature to 298.15 K and vary the relative humidity in the range of 10%–90%. By increasing the relative humidity in a constant temperature field, the vapor mass fraction on the boundaries will increase and the driving force for evaporation (which is the gradient of the vapor mass fraction at the interface)

FIG. 7. (Color online) Velocity vectors around the droplet interface (t ∗ = 1.7).

033305-8

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . .

PHYSICAL REVIEW E 90, 033305 (2014) 300

T T (K)

299 298 297

Temperature ( K )

298 297.5 297 296.5 296 295.5 295 294.5 294 293.5 293 292.5 292 291.5 291 290.5 290

296 295 X X X 294

X

X X X X

293

X

292

XX

X

X

X

X

X

X

X

X

X

X

XX

X

X XX

XX

XX

X XX

XX

0.0112 0.011 0.0108 0.0106 0.0104 0.0102 0.01 0.0098 0.0096 0.0094 0.0092 0.009 0.0088 0.0086 0.0084 0.0082 0.008

X

XX

XX

X

t* = 0.02 t* = 0.17 X

t* = 0.23

291

t* = 0.47

290

t* = 0.93

288

Y Xv V

XX

XX

t* = 0

t* = 1.7

289

(a)

XX

X XX

0

0.2

0.4 0.6 Nondimensional Distance

0.8

1

FIG. 9. (Color online) The temperature profile at different times for the drop and surrounding air.

examined. The evaporation of a two-dimensional droplet in a laminar forced convective air environment is considered here. The droplet is at rest initially. As the air flows over the drop, it will start accelerating and according to the Reynolds and Weber numbers (defined based on the relative velocity and gas properties), the droplet can be either deformed or remains circular. The schematic of the problem is illustrated in Fig. 12. The water and air is selected for the liquid and gas phase, respectively. We set the initial droplet temperature and the inflow air temperature to 300 K. Furthermore, the inflow air is assumed to be dry (Yv,∞ = 0) and the surrounding air contains no vapor content at the beginning of the simulation. Physical

(b) FIG. 8. (Color online) (a) Temperature field and (b) Mass fraction field (t ∗ = 1.7).

305 Psychrometric Chart

will decrease consequently. Hence, vaporizing mass flux will decrease and less heat will flow out of the droplet resulting in a higher wet-bulb temperature. The comparison between the numerical results and the psychrometric chart values for varying relative humidity is presented in Fig. 11. The agreement between the results is satisfactory. The maximum relative error of our numerical results in Figs. 10 and 11 is about 5%. Thus, from the above test cases it can be concluded that the coupling between the heat and mass transfer across the interface is fully consistent. The accuracy of the volumetric mass flux evaluated by Eq. (19) and the capability of the Clausius-Clapeyron equation to relate the temperature to the vapor mass fraction at the interface in our framework is also verified. The previous test cases do not involve external flows. However, in most practical situations a relative velocity exists between the droplet and the background gas. So the consistency of our framework in a convective flow should be

Wet Bulb Temperature ( K )

Simulation

300

295

290

285

280 285

290

295 300 305 310 Dry Bulb Temperature ( K )

315

320

FIG. 10. (Color online) Comparison between the computed wetbulb temperature and the psychrometric chart value at the constant relative humidity of 40%.

033305-9

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

300

Psychrometric Chart

297

Wet Bulb Temperature ( K )

Simulation

294

291

288

285

282

0

20

40 60 Relative Humidity ( % )

80

100

FIG. 11. (Color online) Comparison between the computed wetbulb temperature and the psychrometric chart value at the constant dry-bulb temperature of 298.15 K.

properties of water and air are those given in Table I. The computations are performed on three different grids (coarse: 100 × 300, medium: 200 × 600, and fine: 400 × 1200). In order to quantify our simulations in a convective environment, the Sherwood number on the drop interface is defined as hm D Sh = , (33) Dvg where D is the droplet diameter and hm represents the mass transfer coefficient which can be evaluated as ˙  m , (34) hm = ρ(Yv − Yv,∞ ) where ρ is the bulk density of the gaseous phase. Different empirical correlations are available in the literature [39] which relate the Sherwood number to the Reynolds number, Schmidt number, and mass transfer number (BM ) for evaporating droplets. Convective evaporation of a droplet can be considered in two conditions: low vaporization rates (usually occurs at temperatures sufficiently below the boiling temperature) and

high vaporization rates (occurs at near boiling and/or boiling temperatures associated with high temperature gradient at the drop interface). At low vaporization rates, the effect of the Stefan flow on the thickness of the boundary layer around the droplet can be neglected, so the Sherwood number can only be considered as a function of the Reynolds and Schmidt numbers. However, at high vaporization rates, the Stefan flow affects the flow field and causes the boundary layer thickness to increase which results in the decrease of the Sherwood number. In this situation, the effect of the mass transfer number should be taken into account. Since in the present work, a two-dimensional droplet is considered, one should look for a correlation applicable to the mass transfer over a two-dimensional droplet (liquid cylinder). Heat transfer study around a solid cylinder is a well-documented case and numerous correlations have been presented to date for the average Nusselt number over a solid cylinder in a convective cross flow. Due to the analogy between the heat and species transport equation, experimental heat and mass transfer correlations obtained for a particular geometry can be interchangeable. This can be simply achieved by replacing the Nusselt (Nu) with the Sherwood (Sh) number and the Prandtl (Pr) with the Schmidt (Sc) number. Although the physics of the flow between the convective evaporation of a two-dimensional drop and the convective heat transfer over a solid cylinder are not exactly identical to each other, at low vaporization rates, where the Stefan flow is negligible, the average Sherwood number over an evaporating droplet is comparable to the average Nusselt number over a solid cylinder. The Churchill-Bernstein correlation is a comprehensive equation which can be employed to estimate the average Nusselt number for a cylinder in the cross flow [40]. This equation is recommended for Re · Pr  0.2 and given by

5/8 4/5  0.62Re1/2 Pr1/3 Re Nu = 0.3 + . 1+ [1 + (0.4/ Pr)2/3 ]1/4 282 000 (35) As mentioned above, by replacing Nu with Sh and Pr with Sc, one can obtain a rough estimate for the average Sherwood number over a two-dimensional evaporating droplet in the convective cross flow at low-vaporization-rate regimes:

5/8 4/5  0.62Re1/2 Sc1/3 Re Sh = 0.3 + . 1+ 282 000 [1 + (0.4/Sc)2/3 ]1/4 (36)

Inflow Boundary Free Slip Boundary u x =u ∞

Outflow Boundary

T =T ∞

Yv = Y v ,∞

Droplet Symmetric Boundary

y x

FIG. 12. (Color online) Schematic view of the computational domain.

Figure 13 shows the comparison between the numerical results and empirical correlation of Eq. (36). As expected, by increasing the Reynolds number, the boundary layer thickness decreases and the Sherwood number increases consequently. Numerical predictions are presented for different grid sizes. As shown, by refining the computational grids, the simulation results become closer to the values obtained by the ChurchillBernstein correlation. However, as seen, there exists a definite deviation between the numerical and empirical results over the considered Reynolds numbers and as the Reynolds number increases, this deviation becomes larger. Our simulation

033305-10

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . . 5 Churchill - Bernstein Correlation Coarse Grid Medium Grid

Sherwood Number

4

Fine Grid

3

2

1

3

4

5

Re1/2Sc1/3

6

7

8

FIG. 13. (Color online) Comparison of the predicted Sherwood number and the empirical correlation.

always predicts a lower value for the average Sherwood number than the empirical correlation. This matter can be explained by the fact that the employed correlation has been obtained for a solid cylinder where the gas-solid interface respects the no-slip condition and cannot move under the shear. On the other hand, in our simulation framework the liquid surface moves under the shear stress and this movement decreases the relative velocity between the liquid-vapor interface and the

PHYSICAL REVIEW E 90, 033305 (2014)

surrounding gas which results in the reduction of the Sherwood number. Thus, it can be concluded that our approach also gives consistent quantitative results under convective environments. In order to demonstrate the potential of the present model in more detail, temporal variations of the motion, deformation, and evaporation of the water droplet with the initial diameter of 50 μm in the air stream is investigated in more detail here. The Reynolds and Weber numbers are defined based on the gas-phase properties and the relative velocity between the droplet and inflow air stream. The results are presented for Reg = 78.5 and Weg = 0.5. The characteristics time scale (τ = R/u∞ ) is used to describe the nondimensional time (t ∗ = t/τ ), where R is the initial droplet radius, u∞ is the inflow air velocity, and t is the simulation time. Figure 14 illustrates the relative velocity streamlines around the droplet for different dimensionless times as well as the droplet shape represented by the contour of C = 0.5. At initial times, the droplet is accelerated gently from rest due to the drag force exerted on its surface by the flowing air. As the drop starts accelerating, aerodynamic forces cause it to deform. This deformation is highly dependent on the ratio of the inertial forces to the surface tension force (the Weber number). It has been mentioned that the onset of deformation occurs at We ≈ 0.3 [41]. Thus, as expected, the simulated droplet at We = 0.5 experiences a slight deformation. As a result of the sudden impact of the air flow to the stationary droplet, it takes an oblate shape at the beginning of the simulation (t ∗ = 5). However, since the Weber number is too low, the surface tension force prevents substantial deformation. As can be seen, interactions between the deforming aerodynamic forces and restorative surface tension force lead the droplet to oscillate around its circular shape. The oscillation decays as the drop speeds up and

FIG. 14. (Color online) Relative velocity streamlines and droplet shape for different dimensionless times. 033305-11

SAFARI, RAHIMIAN, AND KRAFCZYK

PHYSICAL REVIEW E 90, 033305 (2014)

FIG. 15. (Color online) Temperature field for surrounding air and droplet at various dimensionless times.

the relative velocity decreases which results in the reduction of the Weber and Reynolds numbers. Figure 14 obviously demonstrates the generation and evolution of internal flow inside the droplet. The shear at the liquid surface causes a motion and circulation of the liquid inside the drop. Figure 15 depicts the temperature fields for the ambient air and droplet at different nondimensional times. As shown, at initial times (t ∗ = 5) that the internal circulation is not formed, the front section of the droplet shows a lower temperature than other parts. This is due to the fact that the minimum thickness of the vapor mass fraction boundary layer around the droplet is located at the forward stagnation point. Thus, the gradient of vapor mass fraction reaches its maximum at this point and the evaporative mass flux is maximized accordingly. Hence, more energy flows out of the liquid-vapor interface leading to a higher reduction in temperature. As time passes, the shear at the liquid-vapor interface makes the liquid move along its surface. Thus, the liquid temperature field is advected with the liquid surface velocity. As a result of this motion, the minimum temperature point pushes to the aft section of

TABLE II. Computational performance of the code in terms of MLUPS (million lattice-node updates per second). Domain size 50 × 50 100 × 100 200 × 200 400 × 400 800 × 800

the droplet. It can be simply seen that eventually the cooled liquid is drawn into the droplet interior due to the internal circulation. These results clearly show how internal circulation can enhance the heat and mass transfer rates over a droplet and decrease the characteristics length and time scales of this phenomenon. In the current article, some fundamental test cases have been employed to prove the ability of the proposed model to account for droplet evaporation phenomena consistently. Our focus was on physical modeling rather than the optimization of the computations and high performance computing. However, it would be advantageous to examine the computational performance of the code. A serial double-precision FORTRAN code has been developed and the simulations have been performed on a personal computer with an Intel Core-i5 CPU 3 GHz with 4 GB RAM. The performance of the developed code is reported in Table II for the different domain sizes for the test case of Fig. 6. The computational performance of Lee’s multiphase model is relatively inferior compared with other lattice Boltzmann multiphase models due to special treatments for the directional and spatial derivatives. Thus, in order to extend the current model to three-dimensional practical engineering applications, making use of modern high performance computing techniques such as graphics processing unit computing would be inevitable.

MLUPS 1.69 1.66 1.56 1.23 1.23

IV. CONCLUSION

Evaporation of a single droplet in quiescent ambient gas and convective gas flow has been modeled and simulated by our extended LBM phase-field multiphase model. In this article, we extended our previous work of Ref. [26] in two

033305-12

CONSISTENT SIMULATION OF DROPLET EVAPORATION . . .

PHYSICAL REVIEW E 90, 033305 (2014)

main aspects. First, evaporation is driven by the vapor species concentration gradient rather than the temperature gradient so that vaporization can occur at any temperature. (Previously we assumed that the liquid is at its boiling temperature [26].) Secondly, an MRT collision operator has been implemented for evaluating the pressure-momentum distribution function. Thus, the stability has been improved and we were able to set all material properties in accordance with their physical values for the water-air system in moderate Reynolds numbers. Various validation tests including one- and twodimensional cases were conducted to validate the consistency of the method. For evaporation of a static water droplet in quiescent air, the calculated wet-bulb temperature was compared to the psychrometric chart results for different air temperatures and relative humidity. Excellent agreement between the simulation and psychrometric chart results has

been observed. In convective test cases, the Sherwood number was used to validate the simulation results. The presented results prove that the added source terms can satisfactorily account for the Stefan flow caused by the vaporization. Moreover, the accurate coupling between the temperature and vapor concentration at the phase interface via the ClausiusClapeyron equation was also verified. Here we have addressed the quantitative simulation of droplet evaporation in the LBM framework with realistic material properties. In this work, we confirmed the consistency and ability of our method in the modeling of vaporization in two-phase fluid flow by considering basic test cases with real physical properties. So it is expected that the presented method can be applied to practical engineering problems. The extension of our approach to three dimensions is technically straightforward and will be the subject of future work.

[1] D. B. Spalding, Symp. (Int.) Combust., [Proc.] 4, 847 (1953). [2] W. A. Sirignano, Fluid Dynamics and Transport of Droplets and Sprays, 2nd ed. (Cambridge University Press, Cambridge, 2010). [3] N. Ashgriz, Handbook of Atomization and Sprays, Theory and Applications (Springer, New York, 2011). [4] S. S. Sazhin, Prog. Energy Combust. Sci. 32, 162 (2006). [5] R. J. Haywood, M. Renksizbulut, and G. D. Raithby, Int. J. Heat Mass Transfer 37, 1401 (1994). [6] D. Juric and G. Tryggvason, Int. J. Multiphase Flow 24, 387 (1998). [7] G. Son, V. K. Dhir, and N. Ramanujapu, J. Heat Transfer 121, 623 (1999). [8] S. W. J. Welch and J. Wilson, J. Comput. Phys. 160, 662 (2000). [9] G. Tormar, G. Biswas, A. Sharma, and A. Agrawal, Phys. Fluids 17, 112103 (2005). [10] D. Z. Guo, D. L. Sun, Z. Y. Li, and W. Q. Tao, Numer. Heat Transfer, Part A 59, 857 (2011). [11] J. Schlottke and B. Weigand, J. Comput. Phys. 227, 5215 (2008). [12] S. Tanguy, T. Menard, and A. Berlemont, J. Comput. Phys. 221, 837 (2007). [13] F. Gibou, L. Chen, D. Nguyen, and S. Banerjee, J. Comput. Phys. 222, 536 (2007). [14] D. Jamet, O. Lebaigue, N. Coutris, and J. M. Delhaye, J. Comput. Phys. 169, 624 (2001). [15] D. Jacqmin, J. Comput. Phys. 155, 96 (1999). [16] H. Ding, P. D. M. Spelt, and C. Shu, J. Comput. Phys. 226, 2078 (2007). [17] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 145 (1992). [18] X. He and L. S. Luo, Phys. Rev. E 56, 6811 (1997). [19] R. R. Nourgaliev, T. N. Dinh, T. G. Theofanous, and D. Joseph, Int. J. Multiphase Flow 29, 117 (2003). [20] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010).

[21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]

[41]

033305-13

B. J. Palmer and D. R. Rector, Phys. Rev. E 61, 5295 (2000). R. Zhang and H. Chen, Phys. Rev. E 67, 066711 (2003). T. Lee and C.-L. Lin, Phys. Rev. E 67, 056703 (2003). Z. Dong, W. Li, and Y. Song, Int. J. Heat Mass Transfer 53, 4908 (2010). A. M´arkus and G. H´azi, Phys. Rev. E 83, 046705 (2011). H. Safari, M. H. Rahimian, and M. Krafczyk, Phys. Rev. E 88, 013304 (2013). D. Chiappini, G. Bella, S. Succi, and S. Ubertini, Int. J. Mod. Phys. C 20, 1803 (2009). D. Chiappini, G. Bella, S. Succi, F. Toschi, and S. Ubertini, Commun. Comput. Phys. 7, 423 (2010). G. Falcucci, S. Ubertini, D. Chiappini, and S. Succi, IMA J. Appl. Math. 76, 712 (2011). T. Lee, Comput. Math. Appl. 58, 987 (2009). T. Lee and L. Liu, J. Comput. Phys. 229, 8045 (2010). T. Lee and P. F. Fischer, Phys. Rev. E 74, 046709 (2006). A. Fakhari and T. Lee, Phys. Rev. E 87, 023304 (2013). S. Hardt and F. Wondra, J. Comput. Phys. 227, 5871 (2008). Z. Guo, B. Shi, and C. Zheng, Int. J. Numer. Methods Fluids 39, 325 (2002). B. Chopard, J. L. Falcone, and J. Latt, Eur. Phys. J.: Spec. Top. 171, 245 (2009). H. Huang, X. Lu, and M. C. Sukop, J. Phys. A: Math. Theor. 44, 055001 (2011). S. R. Turns, An Introduction to Combustion, Concepts and Application (McGraw Hill, New York, 2000). J. Schwarz and J. Smolik, Int. J. Heat Mass Transfer 37, 2139 (1994). T. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. Dewitt, Fundamentals of Heat and Mass Transfer, 7th ed. (Wiley, New York, 2011). L. Hsiang and G. Faeth, Int. J. Multiphase Flow 21, 545 (1995).