Numerical simulation of droplet detachment from solid walls under gravity force using lattice Boltzmann method

Journal of Molecular Liquids 212 (2015) 544–556

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Numerical simulation of droplet detachment from solid walls under gravity force using lattice Boltzmann method S.E. Mousavi Tilehboni a, Ehsan Fattahi b, Hamid Hassanzadeh Afrouzi a,⁎, Mousa Farhadi a a b

Faculty of Mechanical Engineering, Babol Noshirvani University of Technology, P.O. Box 484, Babol, Islamic republic of Iran Department of Computer Science 10 (System Simulation), Erlangen-Nurnberg University, Department of numerical mathematics, Technical University of Munich, Germany

a r t i c l e

i n f o

Article history: Received 17 August 2015 Received in revised form 26 September 2015 Accepted 3 October 2015 Available online xxxx Keywords: Lattice Boltzmann method Droplet Contact angle Eotvos number

a b s t r a c t In this study, interparticle potential model (Shan and Chen) of lattice Boltzmann method is employed to simulate the dripping and detachment of immiscible droplet under gravity force. In order to validate the simulations, free deformation of initialized square drop, coalescence of two adjacent drops and the Laplace law of two dimensional drops have been investigated. Dynamic behaviors of droplet, i.e. droplet detachment from ceiling of channel and droplet detachment from left solid wall of channel with different conditions of wall wetting properties are õinvestigated. Simulation results are presented for different Eotvos numbers (3≤Eo ≤48) at the small Ohnesorge number (Oh = 0.07), as well. The results show that the droplet will not detach from wettable wall in low Eotvos numbers (Eo ≤ 6), but in hydrophobic walls, the droplet detaches from the wall in total range of examined Eotvos number. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Droplet motion and its dynamic behavior is crucial in engineering applications and scientific research. Design of variety of machinery, such as liquid rocket motors and blood pumping machines, requires detailed knowledge of dynamics behavior of droplet in a two-phase flow in macroscopic level [1] and molecular level [2]. The formation, deformation and breakup mechanism of liquid droplets as well as surface tension and wetting related phenomena, can be found in many industrial supplies and are important in many applications such as micro emulsion [3,4], paint sprays [5], ink-jet printers [6], fuel injectors [7], microcapsules [8], DNA arraying [9,10], the manufacture of particles [11], MEMS and microfluidics [12], dispensing processes for imprint lithography [13], and spin etching [14]; an encyclopedic review can be found in Ref. [15]. When an immiscible droplet is placed in contact with a solid wall, there is a contact line between the fluids and the solid surface. The contact angle between the fluids and the wall depends on the properties of all three materials (two fluids and solid phases) coming together. This contact angle can be calculated through Young's equation [16, 17]. There are many works in the field of a two phase flow that investigated the dynamic behavior of droplets with different geometries and conditions. Gunn and Kinzer [18] accurately examined the dynamic behavior of drops. In their study, the terminal velocities for distilled water ⁎ Corresponding author. E-mail addresses: [email protected] (S.E. Mousavi Tilehboni), [email protected] (E. Fattahi), [email protected] (H. Hassanzadeh Afrouzi), [email protected] (M. Farhadi).

http://dx.doi.org/10.1016/j.molliq.2015.10.007 0167-7322/© 2015 Elsevier B.V. All rights reserved.

drops that are falling through static air are experimentally determined. Arecchi et al. [19] published an experimental investigation of the fragmentation procedure of a heavy droplet falling in a lighter miscible fluid. Instability of the droplet is determined based on Schmidt and fragmentation numbers. The results show that variation of Schmidt and fragmentation numbers change droplet stability. They also showed that by increasing of fragmentation number, droplet becomes completely unstable and breakup in to small fragments and finally dissolves in light fluid. Krishna et al. [20] experimentally investigated the effects of the size of transverse cross-sections of vertical channel, on the rising velocity and deformation of bubbles in water. They found that the effect of side walls of the channel on the bubble rising can be ignored while the ratio of bubble diameter to channel width is less than 0.125. Beside the experimental works (some of them are mentioned above), in recent years numerical methods are considered in many studies of multi-phase flows.Han and Tryggvason [21] investigated the unsteady motion of falling liquid droplets for different mixtures of two fluids, using a finite difference/front tracking numerical method. Results are classified for two density ratios, a density ratio equal to 1.15 and a density ratio of ten. The breakup of droplets were controlled by the Eotvos number (Eo), the Ohnesorge number (Oh), and the viscosity and density ratios between the two different fluids. Backward-facing and forward-facing bags were shown in the simulations for the range of Eo b 144; also deformation and breakup regimes are summarized in six Eo–Oh shape maps. Annaland et al. [22] reported the volume of fluid (VOF) method to simulate dynamic behavior of multi bubbles three-dimensionally (3D). It can be seen that due to buoyancy effects, bubble rises and moves to the top of the channel. Interaction effects between two rising bubbles, collision and coalescence together for

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

different initial conditions are simulated as well. Ni et al. [23] presented the direct numerical simulation of falling droplet in a different quiescent liquid. The effects of the walls on the droplet movement have been studied by settling a circular droplet through a static fluid at different lateral positions between two parallel walls. By changing the Reynolds and Weber numbers, the effect of inertial force on the oscillatory motion and deformation of the droplets were investigated and different deformation mechanisms were presented. Jalaal and Mehravaran [24] simulated the fragmentation of falling droplets in a quiescent fluid, using direct numerical simulations. An adaptive VOF technique was employed. Three modes of deformation, disintegration, and dispersion of fragments of a falling droplet were performed by changing Eotvos numbers. In addition, the mechanism of bag breakup was investigated in detail. In the last two decades, the lattice Boltzmann method (LBM) has emerged as a powerful competitor for conventional computational fluid dynamics (CFD) methods. LBM has been applied to a wide range of problems, furthermore it is used to develop new algorithms for the complex fluid and flow simulations, because it is easy to implement, and to parallelize [25–28]. LBM is a powerful method for simulation of complex fluid flow phenomena such as: Nano fluids [29,30], dispersion and deposition of micro-particles [31] and Nano-particles [32], melting of materials [33,34], and porous media [35,36]. Because the lattice Boltzmann equation (LBE) model is based on microscopic models and mesoscopic kinetic equations, it offers many advantages for the study of multi-component or multi-phase flow problems. The four major lattice Boltzmann-based approaches for analyzing two-phase flows are: color-fluid model developed by Gunstensen et al. [37]; interparticle potential model proposed by Shan and Chen [38,39]; free-energy model proposed by, Swift et al. [40] and mean field theory model presented by He et al. [41]. Fakhari and Rahimian [42,43],have presented lattice Boltzmann simulations of falling droplets for a variety of Ohnesorge, Archimedes, and Eotvos numbers. All simulations were performed in 2D or axisymmetric domains. Mousavi Tilehboni et al. [44] simulated the deformation and breakup of a falling droplet under the gravity force, using interparticle potential model of LBM. Different modes of droplet breakup for small Ohnesorge numbers (Oh b 0.1) at several Eotvos numbers were studied. Xing et al. [45] developed a lattice Boltzmann methodbased single-phase free surface model to study the interfacial dynamics of coalescence, droplet formation and detachment phenomena related to surface tension and wetting effects. The limitations of this model are inability of analysis of multicomponent-multiphase flows and incapacitation on the investigation of the actual physics of multiphase flows. Kang et al. [46] studied the displacement of a 2D immiscible droplet under gravity force in a vertical channel, using the lattice Boltzmann method. The dynamic behavior of the droplet was shown and the effects of the contact angle, Bond number, droplet size, density and viscosity ratios were investigated. In their study, the contact angles considered ranged from only 60° to 120° and no explicit relationship between contact angle and adhesion parameter has been proposed. Huang et al. [47] proposed a method for approximating the adhesion parameters in the Shan and Chen multicomponent-multiphase lattice Boltzmann model that leads to the desired fluid–solid contact angle. Instead of Young's equation, this method can be used to involve the interaction parameters (cohesion parameter and a density factor for the fluid–fluid interfacial tension, and the adhesion parameters for the corresponding fluid– solid interfacial tensions). Thus, straightforward and explicit relationship between Young's equation, cohesion parameters and adhesion parameters in the Shan–Chen model has been proposed. The interparticle potential model is based on the concept of nearestneighbor interaction. The model uses two (for two phase flows) different particle distribution functions to model two different components. To model interactions between these components, Shan and Chen modified the collision operator by using an equilibrium velocity, which

545

includes an interactive force. This force guaranties phase separation and introduces surface tension effects. In order to apply the interparticle potential model to non-ideal fluids, Yuan and Schaefer [48] and Qin [49] proposed modified interparticle potentials to obtain a suitable equation of state (EOS). In this paper, the interparticle potential multicomponent-multiphase (Shan and Chen) model has been employed to simulate the dynamic behavior of droplet in contact with the solid wall. The main reasons for using the Shan and Chen model are the following advantages. First, this model is not required to applied complex equation for tracking interface. Second, due to the separated force term on the momentum equation, the external forces such as gravity force and adhesion forces (between solid and fluid) can easily be applied in this model. Issues such as different static contact angle between droplets and solid walls, dripping and detachment of droplets (with different wettabilities) from the solid surfaces, deformation and breakup of droplets (with different wettabilities) after detachment from solid walls are investigated. In the previous works in the Shan and Chen multicomponent– multiphase lattice Boltzmann model determining contact angles from adhesion parameters is strictly empirical and limited to a relatively small range of contact angles. Huang et al. [47] presented an equation for evaluation of contact angle that acquires systematic, consistent contact angle measurements. But the dynamic behavior of droplet and gravity effect is not investigated and results provided for static contact angles. In this study the contact angles of a static droplet have been evaluated for a large range (range of contact angle measured in this paper is 19.6 b θ b 157.6). Besides that, dynamic behaviors of droplet in contact with side walls due to the effect of gravity force have been investigated. 2. Methodology In the multicomponent-multiphase interparticle potential model (Shan and Chen model) [38]discrete Boltzmann equations apply for each component. The standard D2Q9 model with 9 velocity directions is used in this study. Thus, the distribution function of each component is governed by: σ

σ

f i ðx þ ei δt ; t þ δt Þ−f i ðx; t Þ ¼ −


 1
σ σ ðeqÞ f ðx; t Þ− f i nσ ; uσeq τσ i

ð1Þ

σ

where f i ðx; tÞ is the number density distribution of component σ with lattice velocity ei at position x and time t. τσ is the relaxation parameter of component σ and f i function evaluated as: " fi

σ ðeqÞ

¼ wi nσ 1 þ

ei :uσeq c2s

þ

σ ðeqÞ

  ei :uσeq 2 2c4s

−

is the equilibrium distribution

uσeq 2 2c2s

# :

ð2Þ

Discrete lattice velocities in the above equation are: ½e1 ; e2 ; e3 ; e4 ; e5 ; e6 ; e7 ; e8  0 1 0 −1 0 ¼ 0 0 1 0 −1

1 1

−1 1

−1 −1

1 −1

 ð3Þ

and weighting factors are: 8 4 > > > ; > > <9 1 wi ¼ ; > 9 > > > 1 > : ; 36

i¼0 i ¼ 1; 2; 3; 4 i ¼ 5; 6; 7; 8

ð4Þ

546

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

Relaxation parameter (τσ) is related to the kinematic viscosity by: τσ ¼ c2s ðτσ −0:5δtÞ

ð5Þ

The total intermolecular force between particles (two fluids) that acting on component σ at the location x is given by [50] X σ X F σf ðxÞ ¼ −∇V ðx; x0 Þ ¼ −ψσ ðxÞ Gσ σ ψ ðxþei δtÞei : σ ;σ

where υσ is the kinematic viscosity of component σ. δt = 1 and c2s ¼ RT ¼ 1=3 are chosen, where cs is speed of sound in lattice scale. The number density and velocity of component σ can be obtained by nσ ðx; t Þ ¼

X

After collision step, the momentum of σcomponent for each lattice cell should be modified as [50] ρσ uσeq ¼ ρσ u0 þτ σ F σ

σ

f i ðx; t Þ i X σ nσ uσ ðx; t Þ ¼ ei f i ðx; t Þ:

ð6Þ

i

σ

ρ ¼m n

X

σ

ð7Þ

where mσ is the molecular mass of component σ. In order to import nonlocal interaction between particles, Shan and Chen [38] Introduced an interaction potential between components σ and σ V ðx; x0 Þ ¼ Gσσ ψσ ðxÞψσ ðx0 Þ

ð8Þ

where x′ = x + eiδt is the location of the neighbor sites, ψσ and ψσ denotes the intermolecular potential among components σ and σ, respectively. Gσσ is the interaction strength between two fluid components. For two immiscible fluids, when σ and σ are different (indicating two different phases), Gσ σ is positive. As a result, interaction forces between the particles of two phases are repulsive. When σ and σ are the same (particles are for one fluid), Eq. (8) presents the interaction between the particles of one phase, which implies that the nature of the forces are attractive and therefore the value of Gσ σ must be negative. Surface tension between two fluids will control by changing the magnitude of Gσ σ . For nearest neighbor interactions, Gσ σ is given by: Gσ σ ðx; x0 Þ ¼



0 Gσ σ

jx−x0 jNjei δtj jx−x0 j ¼ jei δtj

ð9Þ

The interaction strength parameters G11 and G22 describe the interaction between particles of each component with the neighbors. If we assume component σ as non-ideal gas, then equation of state will be [50]: P ¼ ρσ RT þ Gσσ

RT 2 ½ψðρσ Þ 2

ð10Þ

where P represents the pressure in the equation of state. The original σ

σ

form of the intermolecular potential functions (ψ andψ ) proposed by Shan and Chen [38] is widely used and defined as follows:    −ρσ : ψðρ Þ ¼ ρ0 1− exp ρ0 σ

ð14Þ

where uσeq is the macroscopic equilibrium velocity of component σ, u′ is combinational velocity common to the various components defined as

The mass density is σ

ð13Þ

i

ρσ uσ =τ σ u ¼ X σ σ ρ =τ 0

σ

ð15Þ

σ

and Fσis the total interaction force on fluid component σ including fluid–fluid interaction (Fσf ) fluid–solid interaction (Fσs ) and gravity force (Fσg ), i.e. F σ ¼ F σf þ F σs þ F σg :

ð16Þ

Interaction force between solid wall and fluid (Fσs ) acting on component σ at location x is calculated as [53] X Gads;σ wi Sðx þ ei δt Þei F s σ ðxÞ ¼ −ψσ ðxÞ

where F s σ ðxÞ is the fluid–solid force acting on component σ that is adjacent to the solid wall. Gads, σ is the adsorption coefficient (interaction parameter between solid and fluid) for component σ, and w is the weighting factor as shown in Eq. (4). For the D2Q9 lattice structure, S is a switch, which takes on the value 1 when the lattice node at x + ciδt is a solid node, otherwise it is zero. Interaction of the solid and fluid nodes determines the contact angle of the fluid and solid wall, the adsorption coefficient, Gσads determines the magnitude and type of the interaction force. Based on the magnitude of the adsorption coefficient, three different types of fluids can be specified, namely: non wetting (hydrophobic), neutral, and wetting (hydrophilic) fluids. Usually it is assumed that for Gads , c − Gads , d N 0, (‘c’ denotes the continuous phase and ‘d’ denotes the droplet phase) component of droplet is assumed wetting phase and for Gads , c − Gads , d b 0component of droplet is assumed non-wetting phase. Total mass density and macroscopic momentum of flow are obtained by [50] ρðx; t Þ ¼

X

mσ

σ

X

σ

f i ðx; t Þ

i

ρðx; t ÞU ðx; t Þ ¼

X

mσ

X

σ

ð11Þ

ð17Þ

i

σ

ei f i ðx; t Þ þ

i

1X σ F ðx; t Þ: 2 σ

ð18Þ

Based on equation of state and using the Chapman–Enskog expansion for multicomponent-multiphase model, macroscopic total pressure of flow will be obtained as [50]

In this study ρ0 is set equal to 1. If component σ is assumed to be ideal gas then Gσσ = 0 andψ(ρσ) = ρσ (based on equation of state) [51]. Thus, for this case the interaction strength between two components and the immiscibility of the mixture are described byGσ σ (σ≠σ ). By choosing density ratio of two components close to one, based on the equation of state for both fluids following equation would be satisfied [38,52]

2.1. Implementation of gravity force

ψðρσ Þ ¼ ρσ :

Falling (rising) of droplets (bubbles) occurs due to the gravity force and density difference between two fluids (components or phases).

ð12Þ

P ðx; t Þ ¼ RT

X σ

mσ nσ þ

RT X Gσσ ψσ ψσ : 2

ð19Þ

σσ

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

547

3. Numerical validation 3.1. Relaxation of liquid drops

Fig. 1. Representation of the buoyancy force versus the gravity force. ρdis the density of droplet fluid and ρc is the density of surrounding fluid (continuous phase).

Fig. 1 shows the schematic of gravity and buoyancy effects on droplet suspended in the other fluid. In order to consider the buoyant effect associated with density difference between two fluids, an effective buoyant force, (Fe,σ) was introduced:

The most important and complex part of the two phase flow analysis is concerning the dynamic behavior of interface between two fluids. Thus, as a first step to demonstrate the accuracy of the simulation, investigation of interface behavior that shows the surface tension effect is properly considered. Hence, two tests have been performed. First, an initial two dimensional square drop has been investigated to freely deform to a circular drop as shown in Fig. 2. This phenomenon occurs due to surface tension effects and it is a benchmark validation criterion [16,42,54]. In the second test the coalescence of two circular static adjacent drops are studied (Fig. 3). It has been indicated that when two initial drops put near each other with negligible velocities (assumed to initial static drops), the van der Waals force pulls two drops together and a liquid bridge is formed. The bridge expands rapidly due to surface tension effects. Eventually, the two drops are merged and form a single larger drop with a minimum surface area (circular shape for 2D domain). Details of numerical simulations for two mentioned tests are as follows: In the first test, a square drop with a length of 60 lattice units is placed in the center of the square computational domain with 120 × 120 lattice units. In the second test, two circular drops with Radius 20 lattice units are placed in the center of square computational domain with 120 × 120 lattice units and the central distance between two drops is equal to 42 lattice units. In both tests, periodic boundary conditions are applied in all directions. Initial densities inside the drops are ρd =ρ0 (ρd is density of droplet) and ρc = 0 (ρc is density of surrounding fluid); and outside the drops ρc = ρ0 andρd = 0. The kinematic viscosity of each fluid is 0.1667 in lattice unit. Therefore it can be concluded that in the absence of external forces, a two phase system under surface tension effects, tends to minimize its interface area and a droplet deform to circular shape for a 2D domain (and spherical shape for 3D domain) as shown in Figs. 2 and 3. 3.2. Laplace law

F e;σ ¼ ρσ  g applied :

ð20Þ

In fact, due to the density difference between two fluids, the buoyancy and gravity resultant is:
 F net;g ¼ g eff Δρeff ¼ g eff ρd;eff −ρc;eff :

ð21Þ

Considering Eqs. (20) and (21) the following equation is obtained:

g applied ¼


 g eff ρd;eff −ρc;eff ρσ

ð22Þ

where ρd is the density of droplet and ρc is the density of surrounding fluid, and ρσ = d is the density of the droplet as practically used in the simulation. gapplied indicates the gravitational acceleration that is used in the simulation and applied to density of droplet. Thus, Δρeff and geff are the effective physical density difference and gravitational acceleration quantities. This approach of placing the buoyant effects into a body force greatly expands the range of effective density difference between droplet and surrounding fluid that can be readily simulated with the interparticle potential (Shan and Chen) model.

The surface tension between two fluids (for 2D drops or bubbles) is evaluated based on Young–Laplace's equation given by the following explanation [38,39] ΔP ¼ P in −P out ¼

γ R

ð23Þ

where R is the radius of the drop at the equilibrium state, γ is the surface tension coefficient and ΔP is the pressure difference between the inside (Pin) and the outside (Pout) of the drop. In order to verify the Laplace law, the Laplace drop test checks that the effect of surface tension induces a correct pressure difference at the interface between the two fluids. Initially static drops with different radii are placed inside the horizontal computational domain with 120 × 120 lattice units. Periodic boundary conditions are used for all directions in the computational domain and kinematic viscosity for both fluids is 0.1667 in lattice scale. The densities must be converted to pressures via the Equation of States (Eq. (19)) and compute theΔP. The results show (Fig. 4) a linear relation between pressure difference of the drops and their inverse of radius that satisfies the Laplace law. The slope of the linear line is 0.19 that represents the coefficient of surface tension (in the lattice scale) between two fluids for this study. 3.3. Validation with front tracking method Front tracking method is selected as a CFD method for further validation of the present numerical procedure in multicomponent multiphase

548

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

Fig. 2. Display of free deformation of a static drop in a horizontal plate at different lattice time steps.

problem. Resulting comparisons of falling droplet in a vertical channel are shown in Fig. 5. As shown, results of present study are in good agreement with the results of Han and Tryggvason's [21] which used CFD with front tracking technique. For present LBM results that are shown in Fig. 5, non-dimensional numbers are based on the equivalent spherical droplet.

Huang et al. [47] with respect to the result of Sukop and Thorne [52], presented the following equation for calculation of contact angle in the Shan–Chen model (equivalent to the Young's equation) cosθ1 ¼

Gads;2 −Gads;1 : ρ1 ðxÞ−ρ2 ðxÞ G12 2

ð26Þ

4. Calculation of static contact angle between droplet and solid wall To calculate the contact angle between the droplet and solid surface Young's equation [16] is used that can be written in the following form: cosθ1 ¼

γS2 −γ S1 γ S2 −γ S1 : ¼ γ 12 γ

ð25Þ

This equation determines the contact angle θ1 which is measured in fluid 1 (see Fig. 6). Young's equation for computing the contact angle contains interfacial tension values between the two fluids (γ12) and between each fluid and the solid wall (γS1andγS2). Using interaction parameters between two fluids and between each fluid and solid wall, the Young's equation is equivalent to the interparticle potential model.

According to the Eq. (26), when Gads,2 − Gads,1 N 0 contact angle (θ1) is less than 90° and for Gads,2 − Gads,1 b 0 contact angle is more than 90°. In this study, for different interaction parameters (8 selected cases) between solid and fluids, the contact angle of droplets is considered (Fig. 7). Simulation domain is 200 × 100 lattice cells. Initial condition is a rectangle droplet with 41 × 20 lattice cells in contact with the walls. Relaxation time (τ) is 1 for both fluids. Fig. 7 represent the steady shape of droplet after 24,000 lattice time steps. In order to calculate the contact angle using Eq. (26), each lattice site in the computational domain is occupied by both components while one of the two components is dominant. The minor components can be thought of as dissolved within the dominant component. Therefore, to calculate the contact angle of Fig. 7 (θ1) it is assumed that inside the

Fig. 3. Collision and coalescence of two identical circular drops and free deformation to a circular shape at different lattice time steps.

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

549

Fig. 6. Contact angle.

If the base and height of a droplet are L and H, respectively, the radius of the droplet (R) and contact angle (Ө) are obtained as follows [47]

R¼

Fig. 4. Verification of Laplace's law. Pressure difference between the inside and the outside of a series of drops as a function of inverse drop radius with linear fit of the points.

droplet ρ1(x) is the main density and ρ2(x) is the dissolved density. In the simulations of this studyθ = θ1 = θd represents the contact angle θ1 that is shown in Fig. 6. Also Gads, 1 and Gads, 2 represents the Gads, d and Gads,c respectively. Contact angle between simulated droplet and the solid wall, can be computed by measuring the base and height of droplet on a solid surface (see Fig. 8).


 4H 2 þ L2 8H

;

θ¼

L : 2ðR−HÞ

ð27Þ

Calculated contact angles from young's equation [16] and measured contact angles from present simulations (Fig. 7) are shown in Fig. 9 and Table 1 for different interaction parameters between droplets and solid wall. It can be seen that the lattice Boltzmann simulation results are in very good agreement with Young's equation. For contact angles more than 90° the simulation result approximately fitted on the Young's equation curve. But for the contact angles less than 45° deviation from this curve is increased. For example, in case 1 of Fig. 7, measured contact angle from simulation is 157.6° and calculated contact angle from Young's equation is 156.4°. In addition, in case 7 of Fig. 7 measured contact angle from simulation is 41.4° and calculated contact angle from Young's equation is 46.6° that shows that the deviation is more than that of the previous case. As shown in Fig. 9, relationship between contact angles and adhesion parameters is not linear. 5. Results and discussion 5.1. Dripping and detachment of droplet from solid wall In presence of solid wall, two more boundary conditions must be determined; the wetting properties of the solid wall and the zero velocity at the wall. In order to realize the no slip condition at the solid wall the bounce back scheme is applied. The wetting properties of the solid wall are defined by static contact angle between droplet of fluid and this wall. For contact angle (θd) less than 90°, wall is defined as wettable or hydrophilic (the fluid tends to wet the solid wall and called wetting droplet). For (θd) more than 90°, these properties are defined as unwettable wall or hydrophobic solid surface (the fluid has less affinity for the solid wall and called nonwetting droplet). In the next sections, dripping and detachment of droplets from solid walls are simulated in different wetting properties (different contact angles). An immiscible droplet (denote fluid “d”) attached to the wall is put into the channel filled with surrounding fluid (denote fluid “c”). Initial attached droplets shapes are semicircular with the equivalent diameter (De) equal to 30 Lattice units (initial area of droplet is equal to a circle with diameter De). The length and width of the channel are 360 and 720 Lattice units, respectively. Half way bounce back boundary condition is implemented for all directions. After applying the initial condition, two fluids are affected by gravitational force along the vertical direction of channel. For all following cases the kinematic viscosity is 0.1667 for both fluids. For wetting casesGads, c = 0.4, and Gads, d = − 0.4 are chosen. For nonwetting cases values of Gads , c and Gads , d are − 0.2 and 0.2, respectively. Simulation results are characterized by five dimensionless parameters that are defined as follows: Eotvos number gives the ratio between the gravitational force and the surface tension force:

Fig. 5. Validation of droplet dynamic during falling from vertical channel at Oh = 0.046, ρρd o

μ

¼ 1:15, μ d ¼ 1 and Eotvos number of 12. (A) Present simulation. (B) Front tracking method. o

Eo ¼

gðρd −ρc ÞD2e : γ

550

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

Fig. 7. Simulations of different contact angles between droplets and solid walls. The density of droplet is shown in dark gray. Images represent the steady shape of droplet after 24,000 lattice time steps. For all cases Gads,c = −Gads,d. For cases 1 to 8 Gads,c are −0.4, −0.3, −0.2, −0.1, 0.1, 0.2, 0.3 and 0.4, respectively.

Ohnesorge number is the ratio of the viscous force to the surface tension force; two different Oh ¼

μd 1

ðρd De γ Þ2

based on the droplet fluid:

Time is nondimensionalized with respect to the droplet diameter and the acceleration of gravity t t  ¼
1 : 2 De g

The density and the viscosity ratios, can be selected as the other two parameters ρd μ ; d ρo μo where g is the acceleration gravity, ρd is the density of droplet fluid, ρc is the density of surrounding fluid, De is the equivalent diameter of circular droplet that is set to 30 lattice unit, γ is the surface tension between two fluids, μd and μc are the dynamic viscosity of droplet and surrounding fluid, respectively, and t is the lattice time step of the simulation. In this study Oh = 0.07 is chosen (Oh is very small). Eotvos number is the main controlling parameter to study the dynamic behavior of droplet in small Oh numbers [21]. For all the simulations of this paper, density and viscosity ratio between two fluids are equal to 1.

Fig. 8. Base and height of droplet.

5.2. Droplet detachment from ceiling with different wetting conditions Droplet detachment from hydrophilic (wettable) solid ceiling (θd b 45∘) affected by gravity force, and formation suspended droplet have been shown in Fig. 10 for different Eotvos numbers. Initial droplet attached to the wall is shown at t⁎ = 0 for all cases (a–e) and different shapes of the droplet (for each case) are shown in one line at several selected times. It can be seen (Fig. 10), at the Eotvos numbers equal to 3 and 6, that gravity force cannot overcome the adhesion between solid wall and droplet fluid. Thus, the fluid does not detach from solid ceiling. When Eo = 12, neck region is created on the droplet. By decreasing the pressure on the neck region inside the droplet, diameter of this region decreases and finally droplet detaches from wall (ceiling). Due to the high adhesion force between the wall and droplet, a small volume of the liquid does not detach from ceiling. The remaining liquid retracts back to the ceiling and then spreads on the wettable ceiling since the surface tension is dominant compared with the effect of gravity force on such a small volume of fluid. As shown in Fig. 10(d and e) at higher Eotvos numbers (Eo = 24, 48), droplet will detach from the neck region and fall more easily. At Eo = 48 effects of inertial forces are dominant. Therefore, detachment of droplet from very thin neck region occurs faster and deformation of falling droplet is more than that for low Eotvos numbers. Fig. 11 shows the simulated results of the droplet formation and

Fig. 9. Contact angles of droplets interacting with solid walls versus the magnitude of the adhesion parameters Gads,c. Physical setup of simulation is the same as expressed in the above and results are same as Fig. 7.

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

551

Table 1 Contact angle of droplet for eight selected cases (Fig. 7) with different adhesion parameters (Gads,c = −Gads,d). Different cases of Fig. 7

Gads,c

Contact angle of droplet (θd) calculated from Young's equation(degree)

Contact angle measured from present simulation as shown in Fig. 7 (degree)

Difference between simulation and Young's equation(degree)

1 2 3 4 5 6 7 8

−0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4

156.4 133.4 117.3 103.2 76.8 62.7 46.6 23.6

157.6 134.8 118.1 103.4 75.7 60.3 41.4 19.6

1.2 1.4 0.8 0.2 1.1 2.4 5.2 4

Fig. 10. Simulated snapshots of droplet detachment from wettable ceiling for different Eotvos numbers and nine selected dimensionless times at the fixed Ohnesorge number equal to 0.07. μd ρd ∘ ρ ¼ μ ¼ 1 and contact angle is less than 45° (θd b 45 ) with Gads,c =0.4 (Gads,d = −Gads,c). c

c

552

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

Fig. 11. Simulated snapshots of droplet detachment from unwettable ceiling for different Eotvos numbers and nine selected dimensionless times at the fixed Ohnesorge number equal to μ 0.07. ρρd ¼ μd ¼ 1 and contact angle is more than 90° (90∘ b θd) with Gads,c = −0.2 (Gads,d = −Gads,c). c

c

detachment from a Hydrophobic (unwettable) ceiling wall (θd N 900). In the nonwetting cases, effect of adhesion force between droplet and solid wall is less than wettable wall case. In this case (compared to the wettable ceiling walls) droplet can be detached at the lower gravitational forces that are shown by small Eotvos numbers (Fig. 11, a and b). This means that by increasing the contact angle (θd) effect of gravity force will increase compared to adhesive force. As shown in Fig. 11, for different Eotvos numbers from 3 to 48 (Fig. 11, a–e), effect of gravity force is dominant and no necking step occurs. Similar to the wettable ceiling wall, the fluid tends to minimize its surface area as it detaches from the ceiling. The fluid falling down consists of one whole droplet, and no liquid remains attached to the ceiling since the droplet fluid is nonwetting to the ceiling wall. By increasing Eotvos number, formation of droplet from ceiling and deformation of falling droplet change. Also,

due to the effect of inertial force, it can be seen that the rate of droplet deformation increases by increasing Eotvos number (Fig. 12). When Eotvos number is equal to 3 (Fig. 11a), the droplet deforms into an oblate ellipsoid and moves with a steady state shape. When the Eo is equal to 6 (Fig. 11b), after droplet detachment, first the droplet shape is similar to the previous case, but afterwards the back of the droplet becomes more convex and ultimately the droplet deforms into a thin disklike shape that moves nearly with a steady state shape. When the Eotvos number increases further to 12 (Fig. 11c), the suspended droplet deforms more and eventually forms a backward-facing bag. By increasing Eotvos number to 24, in the falling droplet no sign of bag formation has been shown. At Eo = 48, droplet detaches from wall similar to the previous cases of (Fig. 11,a–d), but for this case the droplet becomes slender at the poles. Finally, fragments of droplet are sheared from the edges

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

553

Fig. 12. Simulated snapshots of droplet detachment from wettable left wall of cannel for different Eotvos numbers and nine selected dimensionless times at the fixed Ohnesorge number equal to 0.07. ρρd ¼ μμd ¼ 1 and contact angle is less than 45° (θd b 45∘) with Gads,c =0.4 (Gads,d = −Gads,c). c

c

and the shear breakup mechanism occurs. Similar phenomena can be seen in Reference [21] for falling droplet. 5.3. Droplet detachment from left wall of channel with different wetting conditions Figs. 12 and 13, show the initial droplets attached to the left solid wall of the channel and move under gravity force. Results are shown for wettable solid wall with contact angle less than 45° in Fig. 12. In this case, the adhesion forces between droplet and solid wall are Significant. Under the influence of adhesion forces, the droplet tends to stick on the wall. Gravity forces pull the droplet down. Due to gravity and adhesion effects between droplet and solid walls, shear forces are created. As shown in Fig. 12, at Eo = 3 (Fig. 12a) the effect of adhesion forces is

dominant. This means that gravity force and therefore shearing effect between solid and fluid cannot overcome the wetting effect. It can be seen that the sticking droplet is spreading and no detachment occurs. When the Eotvos number increases to 6(Fig. 12b) the effect of gravity following the shearing forces increases, however still wetting effects are dominant and results are similar to the previous case. As shown in mentioned cases (Fig. 11a, b) for low Eotvos numbers (Eo b 6) The simulation goes on until steady state is reached. The steady state here is defined as the state that the droplet moves along the wall with a constant velocity and its shape does not significantly change with time. Thus, eventually this shear force overcomes the adhesion forces and after the spreading of droplet along the wall, a portion of droplet detached and is entirely attached onto the bulk fluid. When Eotvos number is 24 (Fig. 12d) initial stages of droplet behavior are similar to the

554

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

Fig. 13. Simulated snapshots of droplet detachment from unwettable left wall of channel for different Eotvos numbers and nine selected dimensionless times at the fixed Ohnesorge number equal to 0.07. ρρd ¼ μμd ¼ 1 and contact angle is more than 90° (90∘ b θd) with Gads,c = −0.2 (Gads,d = −Gads,c). c

c

previous cases, but the detachment of droplet occurs earlier because the shear forces are increased. Simulation results for Eo = 48 shows that in addition to the process considered at the lower Eotvos numbers, after detachment of droplet from solid wall a fine fragment of droplet is sheared from the right edge and the shear breakup mechanism occurs. The rate of deformation after detachment of droplet increases as well. This behavior of the droplet is expected, because at higher Eotvos numbers surface tension force is smaller and gravitational forces are dominant, and results in high deformation of the droplet and eventually breakup. It is clear from Fig. 12, in the wetting condition because of high adhesion forces, droplet tends to stick to the wall. As a result, a portion of droplet doesn't detach and spreads on the wall. Also it can be seen that the detached portion of the droplet at higher Eotvos numbers is larger than that at lower Eotvos numbers. According to Fig. 13, for nonwetting cases (contact angle is more than 90°) because of high inertial and gravitational effects the interface between droplet and solid wall is reduced and the contact line (between the droplet and the

wall) eventually goes to zero. That means that the shear forces between the droplet and the wall are very large and eventually cause the entire droplet to detach from the wall. For these cases droplet does not spread along the wall because the adhesion forces between wall and fluid is negligible compared to the shear forces, therefore at very low Eotvos numbers (Eo = 3 or 6) an entire droplet is detached (Fig. 13a, b). When Eotvos number increased to 12 (Fig. 13c), droplet detachments occurs more rapidly. By increasing the Eotvos number to 24 (Fig. 13d) and then to 48(Fig. 13e) the speed of droplet detachment increases in comparison with lower Eotvos numbers. In addition, because of high inertial effects the rate of droplet deformation increases before detachment from wall. Therefore, because of high-speed droplet detachment and increasing inertia a small portion of droplet is not possible to detach and attach to the wall, due to its low gravity force. As it has shown in Fig. 13, due to wall repulsion, the detached droplet deforms and moves to the centerline of channel. When the droplet approaches the centerline, the pressure difference and the rotation force on the droplet

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

will push the droplet back to the wall. When the droplet is close to the wall, the wall repulsion and the rotation force will drive the droplet away to the channel centerline. Thus, the detached droplet has an oscillating motion. By increasing Eotvos number, the inertia of droplet increases and therefore the surface tension effect decreases. For these reasons at the higher Eotvos numbers the oscillating amplitude is stronger than lower Eotvos numbers. When Eotvos number is very low (Eo = 3) the inertia of droplet and oscillating amplitude is very small due to the large surface tension (Fig. 13a). Simulation results show that the effects of contact angle magnitude and Eotvos number are against each other; as the contact angle increases then the value of Eotvos number for droplet detachment from the wall is decreased compared to the lower contact angle. For the perception concept of this sentence see Figs. 11–13, and compared them.

6. Conclusion In this paper, interparticle potential model (Shan and Chen) of lattice Boltzmann method is employed to simulate the dripping and detachment of immiscible droplet under gravity force. Different static contact angles between droplet and solid wall are simulated and the explicit relationship between Young's equation (including three surface tensions between fluids and solid) and the cohesion and adhesion parameters of Shan and Chen model are presented. Present lattice Boltzmann simulation results are shown for range of contact angle between 19.6 and 157.6 and are in very good agreement with Young's equation. Two phenomena related to the contact angles have been studied: 1) droplet detachment from ceiling of channel and 2) droplet detachment from left solid wall of channel. In the first phenomenon, if the contact angle less than 45° (wettable solid wall), for low Eotvos numbers (3 and 6) adhesion force between droplet and solid wall is dominant (compared to the gravity force) and droplet will not be detached. But by increasing the Eotvos number (12 ≤ Eo), adhesion effect decreases in comparison with the gravity force and droplet detaches from ceiling. Since the wall is very hydrophilic, a portion of droplet does not detach and spreads on the wall. For the contact angle of more than 90°, wettability is very low and the entire droplet will be detached from ceiling in all Eotvos numbers (3 ≤ Eo ≤ 48). In the second phenomena, physical conditions are same as previous cases with the difference that due to shearing force between fluid and solid, contact line increases and the droplet tends to wet the solid left wall.

References [1] R. Clift, J.R. Grace, M.E. Weber, Bubbles, Drops, and Particles, Courier Corporation, 2005. [2] S.M. Langroudi, M. Ghassemi, A. Shahabi, H.R. Nejad, A molecular dynamics study of effective parameters on nano-droplet surface tension, J. Mol. Liq. 161 (2011) 85–90. [3] V. Lisy, B. Brutovsky, A. Zatovsky, A. Zvelindovsky, On the theory of light and neutron scattering from droplet microemulsions, J. Mol. Liq. 93 (2001) 113–118. [4] N.A. Vodolazkaya, Y.A. Kleshchevnikova, N.O. Mchedlov-Petrossyan, Differentiating impact of the AOT-stabilized droplets of water-in-octane microemulsions as examined using halogenated fluoresceins as molecular probes, J. Mol. Liq. 187 (2013) 381–388. [5] K.R. Ellwood, J.L. Tardiff, S.M. Alaie, A simplified analysis method for correlating rotary atomizer performance on droplet size and coating appearance, J. Coat. Technol. Res. 11 (2014) 303–309. [6] S. Hosseini, R. Orsi Koga, N. Ashgriz, S. Chandra, Inkjet printer drop impact on coated and uncoated papers, APS Division of Fluid Dynamics Meeting Abstracts 2013, p. 3005. [7] P.-C. Chen, W.-C. Wang, W.L. Roberts, T. Fang, Spray and atomization of diesel fuel and its alternatives from a single-hole injector using a common rail fuel injection system, Fuel 103 (2013) 850–861. [8] J. Zhang, R.J. Coulston, S.T. Jones, J. Geng, O.A. Scherman, C. Abell, One-step fabrication of supramolecular microcapsules from microfluidic droplets, Science 335 (2012) 690–694. [9] T. Yasui, Y. Inoue, T. Naito, Y. Okamoto, N. Kaji, M. Tokeshi, Y. Baba, Inkjet injection of DNA droplets for microchannel array electrophoresis, Anal. Chem. 84 (2012) 9282–9286.

555

[10] M.A. Burns, B.N. Johnson, S.N. Brahmasandra, K. Handique, J.R. Webster, M. Krishnan, T.S. Sammarco, P.M. Man, D. Jones, D. Heldsinger, An integrated nanoliter DNA analysis device, Science 282 (1998) 484–487. [11] H.M. Shewan, J.R. Stokes, Review of techniques to manufacture micro-hydrogel particles for the food industry and their applications, J. Food Eng. 119 (2013) 781–792. [12] J. Berthier, Micro-drops and digital microfluidics, William Andrew, 2012. [13] I. Reinhold, M.S. Shafran, W. Longsine, M.C. Traub, Y. Srinivasan, V.N. Truskett, W. Zapka, High-speed, low-volume inkjet and its role in Jet and Flash™ Imprint Lithography, NIP & Digital Fabrication Conference, Society for Imaging Science and Technology 2014, pp. 408–412. [14] A. Graf, D. Sonnenberg, V. Paulava, A. Schliwa, C. Heyn, W. Hansen, Excitonic states in GaAs quantum dots fabricated by local droplet etching, Phys. Rev. B 89 (2014) 115314. [15] V. Cristini, Y.-C. Tan, Theory and numerical simulation of droplet dynamics in complex flows—a review, Lab Chip 4 (2004) 257–264. [16] T. Young, An essay on the cohesion of fluids, Philosophical Philos. Trans. R. Soc. Lond. (1805) 65–87. [17] R. Finn, The contact angle in capillarity, Phys. Fluids 18 (2006) 047102. [18] R. Gunn, G.D. Kinzer, The terminal velocity of fall for water droplets in stagnant air, J. Meteorol. 6 (1949) 243–248. [19] F. Arecchi, P. Buah-Bassuah, F. Francini, S. Residori, Fragmentation of a drop as it falls in a lighter miscible fluid, Phys. Rev. E 54 (1996) 424. [20] R. Krishna, M. Urseanu, J. Van Baten, J. Ellenberger, Wall effects on the rise of single gas bubbles in liquids, Int. Commun. Heat Mass Transfer 26 (1999) 781–790. [21] J. Han, G. Tryggvason, Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force, Phys. Fluids 11 (1999) 3650–3667. [22] M. van Sint Annaland, N. Deen, J. Kuipers, Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method, Chem. Eng. Sci. 60 (2005) 2999–3011. [23] M.-J. Ni, S. Komori, N.B. Morley, Direct simulation of falling droplet in a closed channel, Int. J. Heat Mass Transf. 49 (2006) 366–376. [24] M. Jalaal, K. Mehravaran, Fragmentation of falling liquid droplets in bag breakup mode, Int. J. Multiphase Flow 47 (2012) 115–132. [25] S. Ghafouri, M. Alizadeh, S. Seyyedi, H.H. Afrouzi, D. Ganji, Deposition and dispersion of aerosols over triangular cylinders in a two-dimensional channel; effect of cylinder location and arrangement, J. Mol. Liq. 206 (2015) 228–238. [26] G. Kefayati, Natural convection of ferrofluid in a linearly heated cavity utilizing LBM, J. Mol. Liq. 191 (2014) 1–9. [27] G.R. Kefayati, Simulation of non-Newtonian molten polymer on natural convection in a sinusoidal heated cavity using FDLBM, J. Mol. Liq. 195 (2014) 165–174. [28] G.R. Kefayati, Mesoscopic simulation of magnetic field effect on double-diffusive mixed convection of shear-thinning fluids in a two sided lid-driven cavity, J. Mol. Liq. 198 (2014) 413–429. [29] H. Nemati, M. Farhadi, K. Sedighi, E. Fattahi, A. Darzi, Lattice Boltzmann simulation of nanofluid in lid-driven cavity, Int. Commun. Heat Mass Transfer 37 (2010) 1528–1534. [30] A.A. Mehrizi, M. Farhadi, H.H. Afroozi, K. Sedighi, A.R. Darz, Mixed convection heat transfer in a ventilated cavity with hot obstacle: effect of nanofluid and outlet port location, Int. Commun. Heat Mass Transfer 39 (2012) 1000–1008. [31] H.H. Afrouzi, M. Farhadi, A.A. Mehrizi, Numerical simulation of microparticles transport in a concentric annulus by lattice Boltzmann method, Adv. Powder Technol. 24 (2013) 575–584. [32] H. Pormirzaagha, H. Hassanzadeh Afrouzi, A. Abouei Mehrizi, Nano-particles transport in a concentric annulus: a lattice–Boltzmann approach, J. Theor. Appl. Mech. (2015) (Article in Press). [33] C. Yan, L. Hao, A. Hussein, D. Raymont, Evaluations of cellular lattice structures manufactured using selective laser melting, Int. J. Mach. Tools Manuf. 62 (2012) 32–38. [34] M. Jourabian, M. Farhadi, A.A.R. Darzi, Outward melting of ice enhanced by Cu nanoparticles inside cylindrical horizontal annulus: lattice Boltzmann approach, Appl. Math. Model. 37 (2013) 8813–8825. [35] A.A. Mehrizi, M. Farhadi, K. Sedighi, M.A. Delavar, Effect of fin position and porosity on heat transfer improvement in a plate porous media heat exchanger, J. Taiwan Inst. Chem. Eng. 44 (2013) 420–431. [36] A.A. Mehrizi, M. Farhadi, K. Sedighi, A.L. Aghili, Lattice Boltzmann simulation of heat transfer enhancement in a cold plate using porous medium, J. Heat Transf. 135 (2013) 111006. [37] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320. [38] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815. [39] X. Shan, H. Chen, Simulation of nonideal gases and liquid–gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (1994) 2941. [40] M.R. Swift, W. Osborn, J. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (1995) 830. [41] X. He, X. Shan, G.D. Doolen, Discrete Boltzmann equation model for nonideal gases, Phys. Rev. E 57 (1998) R13. [42] A. Fakhari, M.H. Rahimian, Simulation of falling droplet by the lattice Boltzmann method, Commun. Nonlinear Sci. Numer. Simul. 14 (2009) 3046–3055. [43] A. Fakhari, M.H. Rahimian, Investigation of deformation and breakup of a falling droplet using a multiple-relaxation-time lattice Boltzmann method, Comput. Fluids 40 (2011) 156–171. [44] S. Tilehboni, K. Sedighi, M. Farhadi, E. Fattahi, Lattice Boltzmann simulation of deformation and breakup of a droplet under gravity force using interparticle potential model, Int. J. Eng. Trans. A 26 (2013) 781.

556

S.E. Mousavi Tilehboni et al. / Journal of Molecular Liquids 212 (2015) 544–556

[45] X.Q. Xing, D.L. Butler, S.H. Ng, Z. Wang, S. Danyluk, C. Yang, Simulation of droplet formation and coalescence using lattice Boltzmann-based single-phase model, J. Colloid Interface Sci. 311 (2007) 609–618. [46] Q. Kang, D. Zhang, S. Chen, Displacement of a two-dimensional immiscible droplet in a channel, Phys. Fluids 14 (2002) 3203–3214. [47] H. Huang, D.T. Thorne Jr., M.G. Schaap, M.C. Sukop, Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models, Phys. Rev. E 76 (2007) 066701. [48] P. Yuan, L. Schaefer, Equations of state in a lattice Boltzmann model, Phys. Fluids 18 (2006) 042101. [49] R. Qin, Mesoscopic interparticle potentials in the lattice Boltzmann equation for multiphase fluids, Phys. Rev. E 73 (2006) 066703.

[50] X. Shan, G. Doolen, Multicomponent lattice-Boltzmann model with interparticle interaction, J. Stat. Phys. 81 (1995) 379–393. [51] Z. Yu, A Novel Lattice Boltzmann Method for Direct Numerical Simulation of Multiphase Flows, The Ohio State University, 2009. [52] M. Sukop, D.T. Thorne Jr., Lattice Boltzmann Modeling Lattice Boltzmann Modeling, 2006. [53] N.S. Martys, H. Chen, Simulation of multicomponent fluids in complex threedimensional geometries by the lattice Boltzmann method, Phys. Rev. E 53 (1996) 743. [54] H. Bararnia, S. Seyyedi, D. Ganji, B. Khorshidi, Numerical investigation of the coalescence and breakup of falling multi-droplets, Colloids Surf. A Physicochem. Eng. Asp. 424 (2013) 40–51.