PHYSICAL REVIEW E 89, 033311 (2014)

Lattice Boltzmann simulations of multiple-droplet interaction dynamics Wenchao Zhou,1 Drew Loney,1 Andrei G. Fedorov,1,2 F. Levent Degertekin,1,2 and David W. Rosen1 1

The George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0405, USA 2 Parker H. Petit Institute for Bioengineering and Bioscience, Georgia Institute of Technology, Atlanta, Georgia 30332-0405, USA (Received 29 October 2013; published 21 March 2014) A lattice Boltzmann (LB) formulation, which is consistent with the phase-field model for two-phase incompressible fluid, is proposed to model the interface dynamics of droplet impingement. The interparticle force is derived by comparing the macroscopic transport equations recovered from LB equations with the governing equations of the continuous phase-field model. The inconsistency between the existing LB implementations and the phase-field model in calculating the relaxation time at the phase interface is identified and an approximation is proposed to ensure the consistency with the phase-field model. It is also shown that the commonly used equilibrium velocity boundary for the binary fluid LB scheme does not conserve momentum at the wall boundary and a modified scheme is developed to ensure the momentum conservation at the boundary. In addition, a geometric formulation of the wetting boundary condition is proposed to replace the popular surface energy formulation and results show that the geometric approach enforces the prescribed contact angle better than the surface energy formulation in both static and dynamic wetting. The proposed LB formulation is applied to simulating droplet impingement dynamics in three dimensions and results are compared to those obtained with the continuous phase-field model, the LB simulations reported in the literature, and experimental data from the literature. The results show that the proposed LB simulation approach yields not only a significant speed improvement over the phase-field model in simulating droplet impingement dynamics on a submillimeter length scale, but also better accuracy than both the phase-field model and the previously reported LB techniques when compared to experimental data. Upon validation, the proposed LB modeling methodology is applied to the study of multiple-droplet impingement and interactions in three dimensions, which demonstrates its powerful capability of simulating extremely complex interface phenomena. DOI: 10.1103/PhysRevE.89.033311

PACS number(s): 47.11.−j, 68.03.−g, 68.08.−p, 47.55.D−

I. INTRODUCTION

Droplet impact is a century-old problem [1] that has fascinated scientists for its natural beauty and many industrial applications. While single-droplet impact has been extensively studied [2,3], interactions of multiple droplets upon impingement on the substrate have rarely been reported due to the lack of proper research tools. The computational cost for simulating droplet impingement in three-dimensional settings prohibits the study of multiple-droplet impingement. The lattice Boltzmann method (LBM) is a promising alternative to conventional computational fluid dynamics (CFD) simulations for its computational efficiency and has experienced rapid development during the past two decades [4,5]. Originating from lattice-gas automata for simulating gas dynamics [6], the LBM is a particle-based method that treats a group of molecules as a fluid particle rather than continuum and solves the dynamics of particle populations that comes from a microscopic description of the fluid behavior. The LBM has been widely used to model multiphase flow and interface phenomena. A number of approaches for modeling interparticle force have been proposed, of which three models have achieved significant success, including the Shan-Chen model [7,8], which mimics the intermolecular interactions with an empirical forcing potential function to correct the velocity field after each time step; the free-energy-based model [9,10], which incorporates the thermodynamic effects of complex fluids into a modified equilibrium distribution function resulting from a modified momentum flux tensor using the concept of free-energy functional; and the He-ShanDoolen model [11,12], which provides a solid theoretical 1539-3755/2014/89(3)/033311(13)

foundation for a transition from the continuous Boltzmann equation to the lattice Boltzmann equation and a different perspective of viewing the LBM as a special finite-difference approximation of the Boltzmann equation. The well-known numerical instability problem arising from the spurious flow around the interface in multiphase flows, especially for modeling interactions between phases with high-density and/or -viscosity ratio, has also been tackled by different research groups [13–15]. Briant et al. developed free-energy-based wetting boundary conditions to simulate contact line dynamics [16,17]. Yan and Zu studied spontaneous water droplet spreading on both homogeneous and heterogeneous partially wetting surfaces [15]. Although the equlibrium profiles of the spreading droplets matched well with predictions, the transient profiles were not compared with either experiments or theory. Lee and Liu [18] reported transient LBM simulations of the droplet impingement and used a definition of the spreading factor different from that of the experiments in order to match the experimental data. In spite of all the progress, several challenges remain to be addressed in order to achieve satisfactory computational efficiency and accuracy. Although the LBM originated from a microscopic description of the fluid behavior, it could be linked to the macroscopic description through the Chapman-Enskog multiscale analysis [19] recovering the continuous Navier-Stokes (NS) equations. From this perspective, the LBM can be viewed as a special approach of solving the macroscopic Navier-Stokes equations. The phase-field model is a well-established approach for simulating two-phase flows [20] and has been implemented by many researchers in custom-designed computer code as well as in commercial software, such as COMSOL [21]. There is a long

033311-1

©2014 American Physical Society

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014)

history of successful application of the phase-field model in simulating droplet impingement dynamics, including its validation in comparison to experiments [3,22,23]. Therefore, this paper provides an alternative perspective to lattice Boltzmann (LB) formulation and implementation by maintaining consistency between the well-established phase-field NS-based model and the still-developing LB simulation methodology to ensure accurate treatment of the interparticle interactions and boundary conditions in the LBM. First, we derive an expression for the interparticle forcing term for the LBM from the phase-field model. Second, we demonstrate the inconsistency between the existing LB schemes and the phase-field model in calculating the relaxation time at the phase interface and propose a different approximation for the relaxation time as a function of interfacial phase composition function to ensure the consistency with the phase-field interface treatment. Third, we show that the commonly used bounceback scheme [15,16,24] and equilibrium bounceback scheme [18] for the velocity boundary condition with the binary fluid LB scheme do not conserve momentum at the wall boundary and we propose a modified scheme to ensure momentum conservation. Fourth, contrary to the popular surface energy formulation of the wetting boundary condition in the existing LB schemes [16–18], we propose a geometric formulation of the wetting boundary condition to be consistent with the phase-field model. Results show that the geometric formulation enforces the prescribed contact angle in simulations better than the surface energy formulation for both static and dynamic wetting cases. The proposed LB formulation was applied to simulating droplet impingement dynamics in three dimensions and the simulation results were compared with those predicted by the phasefield model, previously reported LBM simulations [18], and experimental data from the literature [25]. Upon incorporation of all above-stated improvements, we show that the proposed LB simulation approach not only yields a significant computational speed improvement over the phase-field model in simulating droplet impingement dynamics on a submillimeter length scale, but also has better accuracy than both the phase-field model and other LB techniques. Finally, the case of multiple interacting droplets, including interfacial dynamics upon droplet impact, spreading, coalescence, and breakup, is successfully simulated and used to demonstrate the capabilities enabled by the proposed LBM simulation algorithm. The rest of the paper is organized as follows. In Sec. II the proposed LBM formulation is presented. Section III presents an approach for calculation of the relaxation time at the interface and the implementation of the boundary conditions. In Sec. IV the LBM simulation results are compared with predictions of the phase-field model and previously reported LBM simulations and validated against the experimental data. Section V demonstrates the capability of the proposed LB numerical model to simulate the multiple droplet impingement and interactions in three dimensions. A summary is given in Sec. VI.

II. LATTICE BOLTZMANN FORMULATION

the behavior of the single-particle distribution function [12], fi − fi(0) ei − u (0) ∂fi + ei · ∇fi = − +F· f , ∂t λ ρcs2 i

(1)

where fi  fi (x, ei , t) is the single-particle distribution function, x is the spatial coordinates, ei is the local particle velocity in the ith direction of a discretized velocity space, λ is the time scale for local particle distribution relaxing back to its equilibrium state, F is the external body force exerted on a particle, ρ is the fluid density, cs is the lattice speed of sound, which is a scaling factor that depends on the specific lattice structure [26], u is the macroscopic velocity of fluid, and fi (0) is the equilibrium distribution function. The equilibrium distribution is approximated by expanding the Maxwellian distribution to the second order of macroscopic velocity u:   ei · u (ei · u)2 u·u (0) , (2) − fi = ti ρ 1 + 2 + cs 2cs4 2cs2 where ti is the weighting factor that depends on the specific lattice model. In order to directly track the pressure field rather than using the equation of state, He et al. [11] introduced a distribution function using the transformation gi = fi cs2 + ψ(ρ)i (0), where i (u) = fi(0) /ρ, ψ(ρ) =p − ρcs2 , and p is the thermodynamic pressure for nonideal fluids, which is only a function of density ρ under isothermal conditions. Equation (1) then becomes gi − gi(0) ∂gi + ei · ∇gi = − + (ei − u) ∂t λ ·[i (u)F + i (0)∇ψ(ρ)],

(3)

where gi(0) = fi(0) cs2 + ψ(ρ)i (0). Next we apply the Chapman-Enskog multiscale analysis to Eq. (3) by expanding gi = gi(0) + gi(1) +  2 gi(2) + O( 3 ),

(4)

where  is identified as the Knudsen number, which is introduced in the Chapman-Enskog theory to keep track of the order of the terms in the series. Higher-order terms represent the influence from a longer time scale and larger spatial scale. Then we can recover the macroscopic equations ∂p + u · ∇p + ρcs2 ∇ · u = 0, ∂t

(5)

  ∂ 1 (1) (ρu) + 2 ∇ ·
(0) = F + ∇ψ(ρ), g +
g ∂t cs

(6)

where
(0) g =
(1) g =





gi(0) ei ei = cs2 pI + ρcs2 uu,

(7)

gi(1) ei ei = −λcs4 [∇ρu + (∇ρu)T ],

(8)

with I the identity matrix and  gi(1) = −ti λ Qi : ∇ρu + ei · ∇ψ(ρ) − ei ∇ : ρuu

The LBM solves macroscopic motion of a fluid by following the evolution of a lattice Boltzmann equation that governs 033311-2

 1 + 2 [(ei · ∇)(Qi : ρuu)] , 2cs

(9)

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

PHYSICAL REVIEW E 89, 033311 (2014)

where Qi = ei ei − cs2 I. It is interesting to point out that
(1) g is recovered from the O( 1 ) order component of gi and is dependent on velocity gradients, which is consistent with the microscopic view of gi(1) . Because the material derivative of the pressure is negligible in simulating incompressible fluids, the first two terms in Eq. (5) can be dropped and Eq. (5) becomes ∇ · u = 0, which recovers the mass conservation equation for incompressible fluids and enforces the incompressibility. Equation (6) can be made equivalent to the momentum equation of the phase-field model by letting ν = λcs2 ,

(10)

F = −∇ψ(ρ) + μ∇C + ρg,

(11)

where ν is the kinematic viscosity, C is the phase composition of the fluid, g is the gravity, and μ is the chemical potential defined as the derivative of the Gibbs free energy with respect to C at constant temperature and pressure. The free energy takes the Ginzburg-Landau form for two-phase flow [27]: E(C,∇C) = β(C − Cl )2 (C − Ch )2 + κ |∇C|2 /2, where the first term on the right-hand side is the bulk free energy E0 that acts to separate the phases, the second term is the gradient energy that favors mixing the phases together, β is a constant relating to bulk free energy, Cl and Ch are two constants representing two different phases (here we choose Cl = 0 for surrounding air and Ch = 1 for liquid), and κ is a parameter related to surface tension σ . Minimizing the free energy of the system leads to a constant chemical potential from the calculus of variation: μ = 4β(C − Cl )(C − Cm )(C − Ch ) + κ∇ 2 C, where Cm = (Cl + Ch )/2. For a plane interface under equilibrium, the composition profile across interface can thus be obtained:   Ch + C l C h − Cl 2z C(z) = + tanh , (12) 2 2 ξ where z is the coordinate normal to the interface √ and ξ is the interface thickness, which is given by ξ = 4 κ/2β/(Ch − Cl ). The surface tension can be calculated by integrating the √ free energy across the interface σ = (Ch − Cl )3 2κβ/6. Because there is large variation in C but little difference in μ across the interface under equilibrium, we replace the surface tension term μC in F [Eq. (11)] with −Cμ using the identity μC = –Cμ [20] to reduce numerical errors around the interface. Note that this forcing term has been used by Lee and Liu [18], but here we show that it can be directly derived from the phase-field model. A second distribution function hi = Cfi /ρ [11,18] is required to track the evolution of C. The partial differential equation governing the evolution of composition C in the phase-field model is the Cahn-Hilliard diffusion convection equation ∂C + u · ∇C = ∇ · M∇μ, ∂t

(13)

where M is mobility. We can derive the evolution equation for hi with Eqs. (1) and (13) to make sure it recovers to the

Cahn-Hilliard equation hi − h(0) ∂hi i + ei · ∇hi = − + (ei − u) ∂t λ   C C ∇ρ i (u) F + ∇C − · ρcs2 ρ + i (u)∇ · M∇μ,

(14)

(0) where h(0) i = Cfi /ρ. Two points should be noted. First, the diffusion process driven by the chemical potential difference is assumed to occur on a longer time scale than the particle collision. Second, to recover the Cahn-Hilliard equation, only the zeroth- and first-order velocity moments of the distribution function hi are used to link the microscopic description of the particle motion to the macroscopic fluid motion, while to recover the mass and momentum conservation equations [Eqs. (5) and (6)] the zeroth-, first-, and second-order velocity moments of the distribution function gi are used. This has important implications for boundary conditions. That is, the boundary conditions for the distribution functions should not change the relationships between the velocity moments of the distribution functions and the macroscopic quantities up to first order for hi and second order for gi . To solve the evolution equations for gi [Eq. (3)] and hi [Eq. (14)], the trapezoidal rule can be applied to integrate over each time step δt along the characteristic direction given by ei . Modified distribution functions g¯ i and h¯ i can be introduced to maintain an explicit scheme in a similar way as done in [11,18]. The density is taken as a linear function of the composition C across the interface: ρ = (ρh − ρl )(C − Cl )/(Ch − Cl ) + ρl . This section has shown the derivation of the LB evolution equations of the distribution functions and the forcing term F from the phase-field model. Although a similar forcing term and LB evolution equations have been used in the literature before (e.g., [18]), we have shown a different approach of deriving the LB evolution equations from the phase-field model by comparing the recovered macroscopic equations with the phase-field equations. The physical meaning of different components of the distribution functions for the proposed LB model is identified. It is found that the secondand higher-order components of distribution functions gi and hi do not affect the macroscopic fluid dynamics because there are no corresponding terms in the recovered macroscopic equations, that is, the macroscopic NS equations only contain the information of the evolution dynamics of zeroth and first order of the distribution functions. That is the zeroth- and first-order distribution functions gi(0) and h(0) i are responsible for the local conservation of mass and of momentum and the composition C, respectively, and gi(1) represents the momentum flux from the neighboring region due to the velocity gradients. The explicit expressions for gi(0) , h(0) i , and gi(1) are given in Eqs. (3), (14), and (9), respectively, and it has no impact on the macroscopic fluid is found that h(1) i dynamics, as there is no corresponding term in the recovered macroscopic equations. Importantly, this understanding of physical meaning of different components of the distribution functions helps formulate the boundary conditions in a selfconsistent manner.

033311-3

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 1. (Color online) (a) Relaxation time and (b) dynamic viscosity as a function of C for low-viscosity liquid (τ l > τ h ) calculated from different equations: (i) Eq. (15), (ii) Eq. (16), (iii) Eq. (17), and (iv) Eq. (18). III. RELAXATION TIME AND BOUNDARY CONDITIONS

In order to apply the LBM to simulate droplet impingement dynamics, we need to address three additional issues that were insufficiently addressed in the literature, including the treatment of the relaxation time τ across the liquid-air interface, the wetting boundary condition, and the velocity boundary for the solid wall. A. Relaxation time

A treatment popular in the LBM literature for the relaxation time τ at the interface is to assume the same linear form as is used for density: C − Cl (τh − τl ) + τl , τ= C h − Cl

of the fluid C across the diffuse interface:   1 1 C − Cl 1 1 + . = − τ Ch − Cl τh τl τl

However, there are inconsistencies between these treatments of the relaxation time in the existing LB schemes and interface modeling in the phase-field model as well as other well-established traditional CFD algorithms (e.g., level-set method [28]). In the phase-field model, the dynamic viscosity η (η = ρτ cs2 δt) variation across the interface is treated to be proportional to C, which leads to τ=

(15)

where τh and τl are the relaxation times for the liquid and the surrounding air, respectively. One deviation from this typical approach is that due to Lee and Liu [18], who proposed an inverse linear form, arguing that the collision frequency should be linearly proportional to variation in the phase composition

(16)

C−Cl (ρh τh − ρl τl ) + ρl τl Ch −Cl . C−Cl (ρh − ρl ) + ρl Ch −Cl

(17)

The profiles of relaxation time and dynamic viscosity calculated from these three different approaches are plotted in Figs. 1 and 2, respectively, for two different situations (low liquid viscosity τ l > τ h and high liquid viscosity τ l < τ h ). As one can see, there are clear differences between the relaxation profiles calculated using these three different approaches.

FIG. 2. (Color online) (a) Relaxation time and (b) dynamic viscosity as a function of C for high-viscosity liquid (τ l < τ h ) calculated from different equations: (i) Eq. (15), (ii) Eq. (16), (iii) Eq. (17), and (iv) Eq. (18). 033311-4

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

PHYSICAL REVIEW E 89, 033311 (2014)

Lee and Liu [18] reported a much slower droplet spreading, not matching experimental data, for low liquid viscosity with the adoption of the linear form [Eq. (15)] than with the inverse linear form [Eq. (16)]. As Fig. 1 shows, the inverse linear form [Eq. (16)] yields relaxation time behavior qualitatively similar to that produced using the linear dynamic viscosity form [Eq. (17)] when τ l > τ h. On the other hand, the inverse linear form agrees poorly with calculations using the linear dynamic viscosity form when τ l < τ h , as shown in Fig. 2. The simplest linear form [Eq. (15)] yields predictions that deviate substantially in both trend and magnitude from those by the other two models for both viscosity limits (Figs. 1 and 2). In addition, the radical change of τ for the inverse linear form when C undershoots (C < 0) as shown Fig. 1 or when C overshoots (C >1) as shown in Fig. 2 during the transient simulations causes an instability issue when the viscosity ratio is high. It is compelling to use Eq. (17) because of its fundamental consistency with the phase-field model. However, as shown in Figs. 1 and 2, the relaxation time exhibits drastic change due to high-density variation when C is around 0. Not only are nonphysical negative τ values observed at low fluid viscosities, but also severe numerical instabilities in solving the LBM equations occur when C varies slightly near 0 during the transient simulation [18]. The numerical instability has been confirmed by our numerical tests for a wide range of conditions. Therefore, as an approach justified on the basis of improving robustness for numerical computations, yet without a significant deviation from the fundamental form inspired by the phase-field model, we propose a quadratic form to alleviate the stability issue associated with oscillations and the negative values produced when Eq. (17) is used:  C−Cl 2 (ρh τh − ρl τl ) + ρl τl C −C τ = h l 2 . (18) C−Cl (ρ − ρ ) + ρ h l l Ch −Cl The construction of the function form ensures that the calculated relaxation profile always remains close to that calculated with the phase-field approach [i.e., Eq. (17)] because both [(C − Cl )/(Ch − Cl )]2 and (C − Cl )/(Ch − Cl ) vary from 0 to 1 with only a slight difference in the rate of change with varying C. The predictions of the relaxation time using Eq. (18) for both situations (τ l > τ h or τ l < τ h ) are also shown in Figs. 1 and 2, respectively, to demonstrate that the quadratic form generally follows the results of Eq. (17) in the off-zero range of composition C values, yet avoids the oscillation and negative values in the immediate vicinity of C = 0, allowing us to ensure robustness of numerical implementation. B. Partial wetting boundary conditions

In most previous studies, a surface energy formulation for computation of the contact angle is employed by taking into account the wall free energy through linear, quadratic, or cubic polynomial approximations [16,18]. A geometric formulation for the wetting condition is mathematically equivalent to the surface energy formulation, but computationally more efficient and more accurate in numerical discretization for the phase-field model [29] because the surface energy formulation

FIG. 3. Illustration of the geometric relationships of a contact angle.

usually causes a discrepancy between the prescribed contact angle and the computed contact angle. Inspired by the adoption of the geometric formulation in the phase-field model, we also adopt the geometric formulation for the proposed LB formulation in order to be consistent with the phase-field implementation of the contact angle treatment. The geometric relationships are illustrated in Fig. 3: ns =

∇C , |∇C|

(19)

where ns is the outer normal to the interface, and π

ns · n n · ∇C tan −θ = = , 2 |ns − (ns · n)n| |∇C − (n · ∇C)n| (20) where θ is the equilibrium contact angle and n is the normal to the wall. A second boundary condition is required to ensure there is no mass flux across the nonpermeable wall: n · ∇μ|wall = 0.

(21)

For a given direction of the wall (i.e., n), Eqs. (20) and (21) specify the Neumann boundary conditions for C and μ, respectively. In order to test its impact on the numerical accuracy, we compared the proposed geometric formulation and the existing surface energy formulation [18] for both static and dynamic wetting. In static liquid wetting of the surface, an equilibrium droplet is generated sitting on a solid wall with different prescribed contact angles ranging from 15° to 150°. The relevant simulation parameters are ρh = 1, ρl = 0.0237, τl = 0.205, τh = 0.075, σ = 0.001 56, and M = 0.02/β. It is assumed that the droplet reaches equilibrium when the kinetic energy per unit volume converges to zero. The equilibrium contact angle is then measured at the contact line using the geometry of the droplet profile (at the interface phase function C = 0.5) obtained from simulations. We tested the difference between the proposed geometric formulation and the commonly used surface energy formulation in enforcing the prescribed contact angle in LB simulations. As shown in Fig. 4, the proposed geometric formulation enforces the prescribed contact angle more accurately than the surface energy formulation in most cases. They both show zero numerical error when the prescribed contact angle is 90°. For impingement dynamics, a droplet is generated with an initial velocity impinging on a solid surface with a prescribed equilibrium contact angle of 107°. The Weber number and the Reynolds number are set to 12.8 and 238, respectively. The contact angle is measured at the contact line using the droplet profile at every time step during the impinging process. As we can observe from Fig. 5, the measured contact angle from

033311-5

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 6. Illustration of particle populations at the wall boundary. The dashed vectors stand for incoming unknown particle populations after the streaming step. FIG. 4. (Color online) Difference between the proposed geometric formulation and the existing surface energy formulation in enforcing the prescribed contact angle in static liquid wetting (θ m is the measured contact angle and θ the prescribed contact angle).

simulations with the surface energy formulation gradually relaxes to the prescribed contact angle during the droplet impinging process while the proposed geometric formulation can enforce the prescribed contact angle consistently during the entire impinging process. C. Velocity boundary conditions

It is straightforward to apply a Dirichlet velocity boundary condition for incompressible Navier-Stokes equations because Navier-Stokes equations explicitly solve for the velocity field. On the other hand, the LBM solves the evolution of particle populations that have more degrees of freedom than what can be imposed by the constraints of the velocity boundary condition. Therefore, it can be complicated to translate the velocity boundary conditions in the continuous formulation into the equivalent boundary conditions for the particle populations. Many previous studies chose a popular bounceback scheme for the nonslip wall boundary conditions

due to its simplicity [15,16], but it has a fundamental deficiency of not conserving momentum at the boundary. Inspired by the boundary analysis for single-phase flow in [30], we follow a similar analysis for the proposed two-distributionfunction LB formulation for two-phase incompressible flow to circumvent a problem with momentum conservation of the existing equilibrium bounceback scheme [18]. The boundary conditions are implemented by finding the appropriate particle population distribution at each computational node near the wall after the streaming step, as illustrated by the dashed vectors in Fig. 6, such that the constraints imposed by the specified velocity at the boundary are enforced. It should be noted that the proposed LBM formulation solves for evolution of the modified distribution functions g¯ i and h¯ i , which are introduced for convenience of calculations and do not have clear physical meanings. However, all the constraints from the boundary conditions only apply to the original distribution functions gi and hi . Therefore, to implement the boundary conditions in the proposed approach to LBM simulations, we need to first obtain gi and hi from g¯ i and h¯ i , apply all constraints to find out the unknown gi and hi , and then calculate the unknown g¯ i and h¯ i . Note that the information from the previous time step is used to obtain gi and hi from g¯ i and h¯ i in order to avoid implicitness, which still gives second-order accuracy in time [31]. We will explain how to determine the unknown gi and hi in the following. After the streaming step, one knows the prescribed velocity u at the boundary, the outgoing particle populations g+ and h+ , and the particle populations parallel to the wall g0 and h0 . The constraints between the particle populations and macroscopic quantities that are required to recover macroscopic equations are   gi ei , (22) p= gi , ρucs2 =
g = 

FIG. 5. (Color online) Difference between the proposed geometric formulation and the existing surface energy formulation in enforcing the prescribed contact angle during the droplet impinging process [θ m (deg) is the measured contact angle and θ (deg) the prescribed contact angle].



(1) gi ei ei =
(0) g +
g ,

hi = C,



hi ei = Cu.

(23) (24)

The unknown incoming particle populations g− and h− are the quantities that need to be computed to satisfy the boundary conditions. First, we can find the macroscopic quantity composition C with known particle populations and Eq. (24):

033311-6

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

PHYSICAL REVIEW E 89, 033311 (2014) TABLE I. List of droplet impingement conditions. Parameter

Case 1 Case 2 Case 3 Case 4 Case 5

impact speed U (m/s) 4.36 4.36 4.36 12.2 12.2 48.8 48.8 50.5 50.5 droplet diameter D0 (μm) 48.8 Weber number We 12.8 12.8 12.8 103 103 Ohnesorge number Oh 0.0151 0.0151 0.0151 0.0148 0.0148 Reynolds number Re 238 238 238 689 689 equilibrium contact angle 31° 90° 107° 31° 107°

FIG. 7. (Color online) Difference between the proposed velocity boundary condition and the existing equilibrium bounceback boundary condition.

C=



Cun =

hi = h+ + h0 + h− ,



hi ei = h+ − h− ,

(25) (26)

where un is the velocity normal to the wall. Therefore, we can obtain 1 C= (2h+ + h0 ). (27) 1 + un Then the density ρ on the boundary can be obtained from C and pressure p can be found from Eq. (22) using the same procedure that we used to find C. With the macroscopic quantities C, ρ, p, and u known on the boundary, one needs to find a set of particle populations gi and hi that satisfies Eqs. (22)–(24). Therefore, we can reconstruct gi and hi as gi = gi(0) (p,ρ,u) + gi(1) , hi = h(0) i (C,ρ,u),

(28)

where gi(1) can be obtained from Eq. (9). This reconstruction of gi and hi makes sure that the constraints (22)–(24) are satisfied and thus the macroscopic equations recovered at the boundary are identical as those recovered in the bulk [i.e., Eq. (6)]. However, the existing equilibrium bounceback boundary [18] only takes into consideration the local momentum conservation, that is, the zeroth order of the distribution functions gi(0) and h(0) i . The first-order distribution function gi(1) is missing, which represents the momentum flux due to velocity gradient and viscosity effects. Therefore, the recovered macroscopic momentum equation will be missing the term
(1) g in Eq. (6) and thus the momentum will not be conserved on the boundary. We have tested the effects of the term gi(1) on droplet spreading, as shown in Fig. 7. The Weber number and the Reynolds number are set to 12.8 and 238, respectively, for the simulation with a prescribed contact angle of 107°. The experimental data are obtained from [25]. One can see clearly that the droplet spreads and retracts more slowly and not as far with the equilibrium bounceback boundary than with the proposed velocity boundary due to the loss of momentum at the boundary. The results illustrate the effects for neglecting the gi(1) component in the boundary formulation and further

confirm the importance of the identification of the physical meanings of different components of the distribution functions and the explicit expression of gi(1) . IV. RESULTS AND VALIDATION

We have implemented a three-dimensional (3D) numerical solver based on the proposed LBM scheme with a D3Q19 lattice model to simulate droplet impingement on a solid surface. The simulations have shown a significant speed improvement for computing droplet impingement dynamics over traditional CFD algorithms based on the phase-field model, as implemented, for example, by the commercial software COMSOL. We have tested five different cases for a single-droplet impingement to compare the LBM simulations with the phase-field simulation results and experimental data from Ref. [25]. The impingement conditions are listed in Table I and are taken from experimental studies, which are used for validation of simulations. The fluid properties of the liquid droplet and the surrounding gas are taken to be those of water and air, respectively, at 1 atm and 25 °C. For all LBM simulations, the computational domain was set to 100 × 100 × 70 and the droplet radius was equal to 25 LB grid units. Although the LBM has the advantage for massive parallelization due to the locality of the computations, all reported simulations were run on a single thread on a laptop PC with a memory requirement slightly over 1 GB for each simulation and were completed within 20 h for each case. For a similar mesh density (i.e., spatial resolution of computations), it is estimated it would take over 1 month for COMSOL to run the same 3D simulation on a 16-core cluster with over 100 GB memory requirement based on our experience. Therefore, all phase-field simulations presented here were performed with a 2D axisymmetrical model and implemented using the COMSOL solver [3,22]. The comparison of the dimensionless spreading factor D*, droplet height H *, and shape coefficient, as predicted by our LBM and phase-field simulations as well as the LBM simulations by Lee and Liu [18], against the experimental data are presented in Figs. 8(a)–12(a) for combinations of parameters summarized in Table I. The spreading factor is defined as the ratio of the diameter of the wetted area to the droplet diameter, the dimensionless droplet height is defined as the ratio of the height of the droplet above the substrate to the droplet diameter, and the shape coefficient is a metric that quantifies the droplet shape achieved relative to a desired shape [22,32]. Note that in Lee and Liu’s results [18], the spreading factor is defined as the ratio of the spreading diameter to the droplet diameter where

033311-7

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 8. (Color online) (a) Validation of the spreading factor D* and dimensionless droplet height H * for case 1 and (b) comparison of the shape coefficient change between LBM and COMSOL simulations for case 1.

the spreading diameter is the blob diameter when the droplet spreads and the diameter of the wetted area when the droplet retracts, which is different from the definition of the spreading factor used in experiments. The graphs in Figs. 8(b)–12(b) show the evolution of droplet shapes using the shape coefficient metric for the LBM simulations based on the algorithm proposed in this work and the benchmark phase-field simulations. We can see that overall the results agree very well in all five cases. There is excellent agreement in the time evolution of the spreading factor between our LBM simulations and the phase-field simulations. The phase-field simulations and the previously reported LBM simulations fail to capture the details of the droplet height change, which is particularly evident in cases 2, 3, and 5 as shown in Figs. 9,10, and 12, respectively, for the times between 10 and 20 μs, which leads to a discrepancy in the shape coefficient, as shown in Figs. 8(b)–12(b). However, overall, the evolution of the shape coefficient, which characterizes the overall droplet shape, matches well with our LBM simulations and the phase-field model, establishing the consistency of our LBM formulation with the phase-field model of continuum NS equations. We can also see that for cases 1 and 2, as shown in Figs. 8 and 9, respectively, both the LBM and phase-field simulations fail to capture the oscillations of the droplet height. One possible reason is that the contact line pinning due to possible

surface contamination could cause the discrepancy of the final equilibrium spreading factor during experiments and affect the overall interface dynamics that is responsible for the droplet height oscillation. The droplet impingement dynamics is dictated by three competing forces due to inertia, surface tension, and viscous stresses, respectively, Fp = ρU 2 D02 , Fσ = σ D0 ,

Fvis = ηU D0 .

(29)

The Weber number can be interpreted as the ratio of Fp to Fσ and the Reynolds number can be considered as the ratio of Fp to Fvis . For example, the inertial, surface tension, and viscous forces are equal to 379.5, 3.68, and 0.55 μN, respectively, for case 5 before the droplet impacts the substrate. Therefore, the droplet spreading is completely dominated by the inertial force in the initial stage of spreading, which results in an increase of the spread factor and a decrease of the droplet height. Because of the small viscous force compared to the surface tension force, very little kinetic energy is dissipated and most of the kinetic energy is converted into interface energy in the form of a significantly deformed droplet interface and extended wetted area of the substrate. It should be noted that a hydrophilic surface converts surface energy into kinetic energy during the wetting process, while a hydrophobic surface requires kinetic energy input to become wetted. The time

FIG. 9. (Color online) (a) Validation of the spreading factor D* and dimensionless droplet height H * for case 2 and (b) comparison of the shape coefficient change between LBM and COMSOL simulations for case 2. 033311-8

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 10. (Color online) (a) Validation of the spreading factor D* and dimensionless droplet height H * for case 3 and (b) comparison of the shape coefficient change between LBM and COMSOL simulations for case 3.

scale of the initial wetting stage can be estimated by D0 /U , which is 4.14 μs for case 5. As shown in Fig. 12, the droplet reaches the maximum spread radius before 10 μs. After the droplet reaches the maximum spread radius, the inertial force significantly decreases due to the conversion of kinetic energy into interface energy and the surface tension force takes control and starts to pull the droplet back to its equilibrium shape dictated by the equilibrium contact angle. Because the equilibrium spread factor is much smaller than the maximum spread factor, the spreading radius starts to retract and converts the interface energy back to kinetic energy. The decrease of the spreading radius also starts to affect the droplet height as the droplet tries to minimize the interface energy by adjusting its shape during the retracting process. The kinetic energy gained during the retracting process may result in an inertial force greater than the surface tension, which can lead to the droplet shape transcending its equilibrium state and further deformation. This would lead to oscillation of the droplet shape with √ the time scale estimated as ρD03 /σ , which is 42 μs for case 5. As shown in Fig. 12, the half cycle of the oscillation is 30 μs, which is consistent with the order of magnitude estimate presented above. There is little dissipation of imparted mechanical energy during the simulation time because of the relatively small viscous force. The time scale of the viscous

dissipation can be estimated by ρD02 /η, which is 2.8 ms for case 5. Therefore, the oscillation of the droplet will continue for many cycles before the energy is fully dissipated and the droplet is stabilized in its equilibrium state. An extensive study of the droplet shape evolution dynamics during impingement can be found in [22]. V. MULTIPLE-DROPLET IMPINGEMENT

The computational efficiency of the implemented LBM solver, based on an improved algorithm described and validated in previous sections, in simulating an individual droplet impingement in three dimensions opens up the possibility for studying multiple-droplet interactions. We exemplify this capability with three different simulation cases, including (i) an impingement of two interacting droplets with a droplet spacing of 80 μm, as shown in Fig. 13; (ii) a-line-of-droplets impingement with a droplet spacing of 65 μm, as shown in Fig. 14; and (iii) a-square-array-of-droplets impingement with droplet spacing of 65, 80, 90, and 100 μm, as shown in Figs. 15–18, respectively. The droplet spacing is defined as the distance between the centroids of the droplets. The droplet diameter is set to 50 μm, the impact velocity is 10 m/s, and the contact angle is 90° for all cases. The Weber number is 100 and

FIG. 11. (Color online) (a) Validation of the spreading factor D* and dimensionless droplet height H * for case 4 and (b) comparison of the shape coefficient change between LBM and COMSOL simulations for case 4. 033311-9

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 12. (Color online) (a) Validation of the spreading factor D* and dimensionless droplet height H * for case 5 and (b) comparison of the shape coefficient change between LBM and COMSOL simulations for case 5.

Ohnesorge number is 0.04, so the simulations are focusing on inertia-dominated hydrodynamics, with surface tension being a more important effect on droplet shape evolution than viscous forces. In the two-droplet interaction case, as shown in Fig. 13, upon impingement on the substrate, the droplets initially behave independently of each other; the dynamics start to change compared to single-droplet impingement when the spreading droplets contact each other. As studied in [22], there are three competing forces that drive the shape evolution of a droplet upon impingement on the substrate, including the inertial force from kinetic energy, surface tension from interface energy, and viscous force from viscosity effects. Upon contact and merging of two or more droplets, a total interfacial area of the droplet ensemble is reduced and the resulting release of excess surface energy is converted mainly into kinetic energy of spreading with minimal viscous energy dissipation (in the case of high We and low Oh numbers). On a longer time scale after the coalescence of the droplets, the kinetic energy of the evolving droplet ensemble is dissipated by viscous forces and the shape of the droplet relaxes back

to the equilibrium shape (i.e., a spherical cap sitting on the substrate) of a single large droplet driven by surface tension. For the case of the line-droplet interactions as shown in Fig. 14, each droplet is neighbored by two droplets along a line with symmetric boundary conditions on the two ends, which effectively simulates the interaction dynamics of an infinitely long linear array of droplets. One notable observation is that the final equilibrium shape is very different from the two-droplet interaction case (Fig. 13) due to the constraint on the lateral spreading imposed by the presence of additional neighboring droplets in an infinite line array. However, similar to the case of two-droplet interactions, conversion of surface to kinetic energy upon initial coalescence of droplets is the most critical event that determines the outcome of the interaction dynamics, which is mainly determined by the droplet spacing under the given impingement conditions. Since the droplet spacing in this case (65 μm) is smaller than that in the twodroplet impingement case (80 μm), the droplets contact each other at an earlier time and thus have greater kinetic energy at the merger of the interfaces. Therefore, droplet shapes deform greatly with the emergence of fluid ridges and fingers during coalescence, with significant differences between the shapes

FIG. 13. Demonstration of two-droplet impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 80 μm). See Ref. [33] for an animation of the two-droplet interaction dynamics.

FIG. 14. Demonstration of a-line-of-droplets impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 65 μm). See Ref. [33] for an animation of the a-line-of-droplets interaction dynamics.

033311-10

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

PHYSICAL REVIEW E 89, 033311 (2014)

FIG. 16. Demonstration of an-array-of-droplets impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 80 μm). See Ref. [33] for an animation of an-arrayof-droplets interaction dynamics with a droplet spacing of 80 μm.

FIG. 15. Demonstration of an-array-of-droplets impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 65 μm). See Ref. [33] for an animation of an-arrayof-droplets interaction dynamics with a droplet spacing of 65 μm.

of the droplets in Figs. 13 and 14 as the droplet ensembles evolve towards their equilibrium states. In the array-droplet interaction, each droplet is surrounded by eight droplets, which are arranged into a square array with symmetric boundary conditions on all the edges, as shown in Figs. 15–18. Therefore, an infinitely large two-dimensional array of droplets is effectively simulated. Due to the additional neighboring droplets, the interaction dynamics become much more complex than in the previous cases and the final equilibrium shape is different for arrangements with different droplet spacing (65, 80, 90, and 100 μm), as illustrated in Figs. 15–18. Unlike in the previous two cases, in which only two droplets coalesce at one location, four droplets meet at the same location in this case. Therefore, during droplet coalescence, even more interface energy is converted into kinetic energy that drives dramatic deformation of coalescing droplet shapes, including splash and bounceback phenomena, as evident in comparing Figs. 14 and 15. Indeed, as shown in Fig. 15, the kinetic energy gained from the coalescence of the droplets is so large that it drives the droplet shape to deform to such an extent that the droplets break up and generate a stream of satellite droplets. If the droplet spacing is increased to 80 μm as shown in Fig. 16, periodic depressions in the liquid film are formed at the locations where the droplets meet

due to trapping of air bubbles during droplet coalescence. If the droplet distance is further increased to 90 μm as shown in Fig. 17, the droplets contact each other when the kinetic energy is already fairly small and coalescence of the droplets is largely driven by surface tension. As a result, a uniform film is formed as the final equilibrium shape. If the droplet distance is further increased to 100 μm, the droplets meet each other at a time when the droplets almost start to retract and the retracting kinetic energy is large enough to overcome an energetic advantage of reduced surface energy associated with coalescence of the droplets, causing the droplets to separate at equilibrium, as shown in Fig. 18. These results further illustrate the importance of droplet spacing and impingement velocity in multiple-droplet interaction dynamics. Of course, the dynamics will change in the case of fluids with different properties, such as viscosity and surface tension, as different

FIG. 17. Demonstration of an-array-of-droplets impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 90 μm). See Ref. [33] for an animation of an-array-of-droplets interaction dynamics with a droplet spacing of 90 μm.

033311-11

ZHOU, LONEY, FEDOROV, DEGERTEKIN, AND ROSEN

PHYSICAL REVIEW E 89, 033311 (2014) VI. CONCLUSION

FIG. 18. Demonstration of an-array-of-droplets impingement in the unit of dimensionless time (We = 100, Oh = 0.04, and the droplet spacing is equal to 100 μm). See Ref. [33] for an animation of an-array-of-droplets interaction dynamics with a droplet spacing of 100 μm.

energetic interactions may be more dominant during droplet impingement and coalescence events. This would require an in-depth investigation and is a subject left for future work. As can be seen from the evolution of droplet shapes presented in Figs. 13–18, which capture the events of droplet coalescence and breakup, the interface topology and dynamics are highly complex, especially with an increase in the number of interacting droplets. Importantly, the proposed numerical LBM solver effectively captures the computational complexity of the physical problem and preserves some inherent symmetries of each simulation case, even though the simulations were performed in a general 3D format. These results demonstrated the LBM capability for handling highly complicated interface dynamics with the potential to become a powerful simulation tool to deepen our understanding of interfacial phenomena and, when implemented in massively parallel fashion, become a digital design tool for dynamic fluid interfaces of interest to many practical applications.

[1] M. Worthington, Proc. R. Soc. London 25, 261 (1876). [2] M. Rein, Fluid Dyn. Res. 12, 61 (1993). [3] W. Zhou et al., Virtual Phys. Prototyping 7, 49 (2012). [4] K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010). [5] S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech. 30, 329 (1998). [6] G. R. McNamara and G. Zanetti, Phys. Rev. Lett. 61, 2332 (1988). [7] X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993). [8] X. Shan and H. Chen, Phys. Rev. E 49, 2941 (1994). [9] M. R. Swift et al., Phys. Rev. E 54, 5041 (1996).

An approach to formulate and implement the lattice Boltzmann method for the fluid mechanics problem involving interfacial flows was proposed based on reconciling the deficiencies of existing methodologies and establishing fundamental consistency with the continuous phase-field Navier-Stokes formalism. An interparticle forcing term was derived from the phase-field model. A robust expression for calculating the relaxation time at the phase interface was developed to remedy the inconsistency between the existing LB schemes and the phase-field model. In addition, a geometric wetting boundary condition was proposed and demonstrated to be numerically more robust and accurate than the commonly used surface energy formulation. Moreover, a careful examination of the equilibrium bounceback velocity boundary condition revealed that it does not conserve momentum at the wall for twophase flow simulations. This critical deficiency resulting in underprediction of the interface evolution rate in conventional LBM implementation was reconciled using a different velocity boundary condition that ensured momentum conservation on the boundary. A numerical solver was implemented based on the improved LBM formulation for simulating droplet impingement dynamics. Simulation results were compared against the predictions using the continuous phase-field CFD simulations, the LBM simulations reported in the literature, and experimental data. It was found that the proposed LBM formulation not only yields significant improvement in computational efficiency, but also produces more accurate predictions than the phase-field model and the previously reported LBM models. The capability of the proposed LBM solver in handling highly complex interface dynamics in three dimensions provides a foundation for its utility in many practical applications involving droplet deposition on a solid substrate. Several examples of simulating multiple-droplet interaction dynamics were given to demonstrate the powerful capability of the developed solver. ACKNOWLEDGMENTS

We gratefully acknowledge the US National Science Foundation for support through Grant No. DMI-0900322. We would also like to thank Dr. Hongming Dong for providing the original experimental data and Chad Hume for discussion.

[10] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett. 75, 830 (1995). [11] X. He, S. Chen, and R. Zhang, J. Comput. Phys. 152, 642 (1999). [12] X. He, X. Shan, and G. D. Doolen, Phys. Rev. E 57, R13 (1998). [13] T. Inamuro et al., J. Comput. Phys. 198, 628 (2004). [14] T. Lee and C.-L. Lin, J. Comput. Phys. 206, 16 (2005). [15] Y. Y. Yan and Y. Q. Zu, J. Comput. Phys. 227, 763 (2007). [16] A. J. Briant, A. J. Wagner, and J. M. Yeomans, Phys. Rev. E 69, 031602 (2004). [17] A. J. Briant and J. M. Yeomans, Phys. Rev. E 69, 031603 (2004). [18] T. Lee and L. Liu, J. Comput. Phys. 229, 8045 (2010). [19] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory

033311-12

LATTICE BOLTZMANN SIMULATIONS OF MULTIPLE- . . .

[20] [21] [22] [23] [24]

[25] [26] [27]

PHYSICAL REVIEW E 89, 033311 (2014)

of Viscosity, Thermal Conduction and Diffusion in Gases (Cambridge University Press, Cambridge, 1991). D. Jacqmin, J. Comput. Phys. 155, 96 (1999). C. Zhou, P. Yue, and J. J. Feng, Int. J. Multiphase Flow 34, 102 (2008). W. Zhou et al., AIChE J. 59, 3071 (2013). K. Khatavkar et al., J. Fluid Mech. 581, 97 (2007). H. Kusumaatmaja and J. Yeomans, in Simulating Complex Systems by Cellular Automata, edited by J. Kroc, P. M. A. Sloot, and A. G. Hoekstra (Springer, Berlin, 2010), pp. 241–274. H. Dong et al., AIChE J. 53, 2606 (2007). J. Latt, Ph.D. thesis, University of Geneva, 2007. P. Yue et al., J. Comput. Phys. 219, 47 (2006).

[28] E. Olsson and G. Kreiss, J. Comput. Phys. 210, 225 (2005). [29] H. Ding and P. D. M. Spelt, Phys. Rev. E 75, 046708 (2007). [30] J. Latt et al., Phys. Rev. E 77, 056703 (2008). [31] T. Lee, C.-L. Lin, and L.-D. Chen, J. Comput. Phys. 215, 133 (2006). [32] W. Zhou et al., in Proceedings of the ASME 2012 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference Chicago (ASME, New York, NY, 2012). [33] See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevE.89.033311 for the animations.

033311-13