Simulations of coalescence of two colliding liquid drops using lattice Boltzmann method

Abstract

In this paper, two-dimensional numerical simulations of the head-on and off-center binary collision of liquid drops are carried out using lattice Boltzmann method. The coalescence process of drops colliding is freely captured, regardless of tracking equations and drops model. Special attention is paid to investigate the effect of the Weber number, impact velocity, drop size ratio on the coalescence process. The research results demonstrate that numerical results of lattice Boltzmann method are in agreement with qualitatively experimental data in the same Weber number. It is also noted that the first oscillation period is invariably shorter than the second period, regardless of the Weber number. It is further found that variation of the first period of oscillation in off-center collision of two equal-size drops as function of Weber number are higher than the variation of the first period of oscillation in head-on collision of two equal-size drops in the same Weber number.

Introduction

The phenomenon of binary droplet coalescence is an indispensable characteristic in many industrial processes.^1,2 The understanding of this process as well as the prediction and control of the phenomenon is overwhelmingly crucial practice. Liquid–liquid extraction, emulsification, spray coating, hydrocarbon fermentation and waste treatment are many familiar applications. ³ Especially, the coalescence of oil drops is closely related to a host of industrial and environmental clean-up operations. The superiority of liquid–liquid extraction operations lies on coalescence of the discrete drops. ⁴ In addition, the coalescence of the emulsions formed during alkaline flooding operations is essential. The coalescence of fine oil mist is achieved using porous coalescences, settlers, chemicals, and electric fields to break up the emulsions in the petroleum refineries. The collision and coalescence of liquid drops in diesel engine sprays are overwhelmingly shared, because the size of the drops is affected by coalescence in the engine cylinder, which in turn may affect its performance and emission characteristics. ³

The study of binary droplet collisions has been widely utilized across different scientific areas, from understanding cloud formation in climate theory, to engineering applications such as spray coatings, turbine blade cooling, and spray combustion in diesel internal combustion engines. Modeling of multiphase flows with traditional computational fluid dynamics (CFD) is challenging as it requires the tracking of the interface between the different phases. ¹¹ This becomes increasingly arduous in regimes where two interfaces split or meet and merge together, as is the case with droplet collisions. Previous studies of lattice Boltzmann method (LBM) have emphasized potential advantages to capture interface between gas and liquid phase.^4–9 In contrast to solvers that directly compute solutions for the discretized Navier–Stokes equations, LBM is a form of cellular automaton. It relaxes the incompressibility criterion and thus does not require an additional step to compute the pressure with an iterative method such as a multi-grid solver, ⁵ or a pressure projection step. ⁷ Multiphase surface fluids can also be computed with an approach known as smoothed particle hydrodynamics (SPH),^6,7 which does not require a fixed grid and computes the fluid properties by computation kernels defined on particle neighborhoods. These properties of the algorithm make it feasible to perform interactive fluid simulations, or adapt it to other problems.^6,13

A multitude of comprehensive experimental studies of binary droplet collisions are presented in Law et al., ¹⁴ including one by Qian and Law. ¹² Various regimes of droplet collisions are confirmed by two dimensionless parameters. The first is the Weber number We, which is expressed as

$$\mathrm{We}=\frac{\rho u_{\mathrm{rel}}^{2}D}{\sigma }$$

(1)

where D is the droplet diameter, $u_{\mathrm{rel}}$ is their relative velocity, and σ is the surface tension of the liquid. The surface tension follows as the integral along a flat surface of the mismatch between the normal $P_{\mathrm{xx}}$ and transversal $P_{\mathrm{yy}}$ component of the pressure tensor. ⁷

$$\sigma ={\int }_{-\infty }^{+\infty }(P_{\mathrm{xx}}-P_{\mathrm{yy}})\text{d}x=-\frac{{\mathrm{Gc}}_s^4}{2}{\int }_{-\infty }^{+\infty }{\left|{\partial }_x\psi \right|}^2\text{d}x$$

(2)

The second is the impact parameter B

$$B=\frac{X}{D}$$

(3)

where X is the separation between the centers of the droplets, perpendicular to their direction of motion. Different regimes for the droplet collisions are found, and a general representation of these regions is given in Figure 1.The difficulty in reproducing the results of Figure 1 with the LBM is associated with reaching the relevant values for the Reynolds number. The Reynolds number (Re) is given by

$$\text{Re}=\frac{\rho u_{\mathrm{rel}}D}{\nu }$$

(4)

Figure 1.

Schematic of (a) the head-on and (b) off-center binary collision models as viewed in the reference frame in which one of the two drops is at rest.

Based on the above discussions, the research mainly focusses on the effects of the Weber number, impact velocity, drop size ratio on the coalescence process, coalescence process of drops colliding is freely captured, regardless of tracking equations and drops model. In the following, a multiphase lattice Boltzmann model will be briefly described at first. Then the problem definition and numerical verification are presented. After that, the detailed numerical studied of the head-on and off-center binary collision of van der Waals liquid drops are given. Finally, some concluding remarks are provided.

Numerical method

The starting equation of lattice Bhatnogar–Gross–Krook (LBGK) models for phase transitions read as ¹⁰

$$f_k(x+c_k\Delta t,t+\Delta t)=f_k(x,t)+(f_k^{\mathrm{eq}}(x,t))-f_k(x,t))/\tau +\Delta f_k$$

(5)

where f_k is the particle density with given velocity c_k, τ the relaxation time, and k is the index of discrete velocity. $\Delta f_k$ is a fictitious forcing term, describing interactions between two neighboring sites, this term will be specified later. c_k is defined by the following choice for the D2Q9 model. ¹⁰

The equilibrium distribution function is chosen as

$$f_k^{\mathrm{eq}}=\rho {\omega }_k[1+\frac{c_{k\alpha }{u}_{\alpha }}{c_s^2}+\frac{c_{k\alpha }{c}_{k\beta }{u}_{\alpha }{u}_{\beta }}{2c_s^4}-\frac{u^2}{2c_s^2}]$$

(6)

where ${\omega }_k$ is the associated weighting coefficient. The corresponding values are ${\omega }_0=\frac{4}{9}$ , ${\omega }_k=\frac{1}{9},k=1,2,3,4$ and ${\omega }_k=\frac{1}{36},k=5,6,7,8$ . c_s is a constant with $c_s=1/\sqrt{3}$ . ¹⁰ The density and velocity of the fluid are defined as

$$\rho =\sum _{k=0}^8{f}_k=\sum _{k=0}^8{f}_k^{\mathrm{eq}},\rho u=\sum _{k=0}^8{f}_k{c}_k=\sum _{k=0}^8{f}_k^{\mathrm{eq}}{c}_k$$

(7)

The body force $\Delta f_k$ is the difference of equilibrium distribution functions according to the mass velocity after and before the action of a force in time step at constant density ρ.

$$\Delta f_k=f_k^{\mathrm{eq}}(\rho ,u+\Delta u)-f_k^{\mathrm{eq}}(\rho ,u)$$

(8)

where $\Delta u=F·\Delta t/\rho$ , the attractive forces are introduced between two neighbor nodes. ⁷ In this paper, the van der Waals equation of state is used and its reduced form is expressed as Qian and Law. ¹²

Using the Chapman–Enskog expansion, we derive the hydrodynamic equation up to the second order.

$$\frac{\partial \rho }{\partial t}+\nabla ·(\rho \mathbf{u})=0$$

(9)

$$\frac{\partial (\rho \mathbf{u})}{\partial t}+\nabla ·(\rho \mathbf{uu})=-\nabla p+\nabla ·(2\rho \nu \mathbf{S})$$

(10)

Numerical simulations

In the first section, the effect of variation of the Weber number on the outcome of the head-on collision of two equal-size drops is investigated. Numerical simulations are carried out for several cases. In the first and second cases, two resting cylindrical droplets of identical radii of R₀ lattice units (l.u.), initially separated by 4 l.u., are placed into a periodic rectangular domain of size $512\times 512$ l.u. The first case is included two simulations with two Weber number. In the first simulation the velocity of left liquid drop is $0.05$ , the velocity of right liquid drop is −0.05. The Weber number is 2.5. In the second simulation, the velocity of left liquid drop is 0.1, the velocity of right liquid drop is −0.1. The Weber number is 10.

In the second section, the effect of variation of the Weber number on the outcome of the offer-center collision of two equal-size drops is studied. A simple off-center collision numerical simulations is obtained by shifting center of the incident drop in the vertical direction a distance equal to its radius R, resulting in a dimensionless impact parameter $X=0.5$ .The initial condition of numerical simulation is same with the first case.

The temporal evolution of the collision for the drops using LBM for We = 2.5 and 10 is presented in Figures 2 and 4. Both cases evolve in a qualitatively similar fashion in spite of the difference in their Weber values. It is seen that the former results in two opposite flows within the coalescing drop directed toward the central plane perpendicular to the x-axis. The pressure within this region soon exceeds the surface tension pressure, causing an outwardly spreading flow along the y-axis. As this motion progresses, the overall width of the coalesced drop reduces with a consequent increase of the rim pressure. The outward motion seen in the present simulations is the analogous of the outwardly radial flow envisaged in actual liquid-drop collision experiments ¹⁴ in Figures 3 and 5. The research mainly demonstrates that numerical simulations are in qualitative agreement with experimental results and the combined drop in their largest surface deformation.

Figure 2.

Sequence of times showing the evolution of the head-on collision using LBM ( $\mathrm{We}\approx 2.5$ ).

Figure 3.

Sequence of times showing the evolution of the head-on collision experiment ¹⁴ ( $\mathrm{We}\approx 2.5$ ).

Figure 4.

Sequence of times showing the evolution of the head-on collision using LBM ( $\mathrm{We}\approx 10$ ).

Figure 5.

Sequence of times evolution of the head-on collision experiment ¹⁴ ( $\mathrm{We}\approx 10$ ).

Since the magnitudes of the internal velocities change quicker along the y-axis due to the higher surface energies involved, fluid motion undergoes stronger viscous dissipation in this direction than along the x-axis. The model of Figure 4, with We = 10, formed a thinner disk compared to that of Figure 2 having lower (We = 2.5) collision values. Also, note that for the former model the disk possesses well-pronounced concave surfaces on each side around the center-to-center line of the parent drops. Due to the higher initial kinetic energy for this case, the surface around the center-to-center line continues to move inwards until a concave surface is formed which then produces adverse pressures that prevent further inward motion of the surface. This result is in excellent agreement with previous simulations by Fockea et al. ¹⁵ with the aid of different numerical techniques, who mainly found that higher collision values of Weber translate into larger surface deformation.

Development of the disk ceases as soon as the rim pressure first balances and then overcomes the stagnation pressure, which leads to the disk to contract back under surface tension. The bottom-row slides of Figures 2 and 4 manifest the reverse motion of the surface toward completion of the first oscillation period. In particular, the last slides in Figures 2 and 4 depict the shape of the coalesced drops close to the end of the first period by the time the maximum elongation happens to be along the x-axis. These forms correspond to the stretched liquid cylinder detected in drop collision experiments after the disk contraction phase. It is found that higher Weber number collision values give rise to more elongated coalesced drops by the end of the first period of oscillation. The subsequent evolution will be governed a long-term interplay between viscous dissipation and conversion into surface energy of the internal liquid movement, with a consequent damped oscillatory motion with maximal drop elongations alternating between being parallel and perpendicular to the direction of incidence. This phase will eventually end with the formation of a circular coalesced drop. Although direct quantitative comparisons cannot be done with other existing simulations, several of the observed trends seem to be similar in spite of differences in the initial parameters and numerical techniques employed, implying that the LBM is an alternative and promising numerical scheme for simulating the dynamics of colliding drops.

A sequence of slides are displayed for the time evolution of the off-center binary of collision numerical simulations in Figures 6 and 8 when Weber number is equal to 2 and 10, respectively, direct comparison can be made with the head-on evolution of Figures 5 and 6. During the initial condition, a bridge between the drops forms again at the point of contact, which then expands in radius in as much as the same way as described for the head-on the We = 10 case. While part of the sliding motion is resisted by viscous forces due to the shearing flow layer between the sliding masses, part is transformed into rotational motion by the action of the surface tension forces which tend to circularize the drop. The centrifugal forces associated with rotation of the deforming end caps sandwiching the stagnation region along with the outwardly accelerated flow within it, induced by the initial kinetic energy pertaining to the longitudinal velocity, flatten the combined drop until a stage of maximum deformation is reached. This point is achieved in the low Weber number (We = 2.5) collision (Figure 6) and for We = 10 (Figure 8), as shown in each figure by the second slide of middle row. The peanut shape has undergone transition into a rotating dumbbell mode, because of the larger momentum effects associated with the latter case. A common feature in all models is that due to the elongated shapes, the centrifugal forces first induce solid-body rotation of the coalesced drop regarding its center of mass. It can be demonstrated that numerical simulations of LBM in We = 10 qualitatively agree with experimental results.

Figure 6.

Sequence of times evolution of the off-center collision using LBM ( $\mathrm{We}\approx 2.5$ ).

Figure 7.

Sequence of times evolution of the off-center collision experiment ¹⁴ ( $\mathrm{We}\approx 2.5$ ).

Figure 8.

Sequence of times evolution of the off-center collision using LBM ( $\mathrm{We}\approx 10$ ).

Figure 9.

Sequence of times evolution of the head-on collision experiment ¹⁴ ( $\mathrm{We}\approx 10$ ).

The coalescence collision is taken into account and the effects of the Weber number, drop size ratio, impact velocity, and internal circulation on the behavior of the combined drop are presented in head-on and off-center collision of two equal-size drops. In collision of equal-size drops, an increase of the Weber number results in increasingly thinner drops during the coalescence process. When collision of non-equal-size drops is taken into account, the early stages of coalescence circumfuse considerably larger deformations for the smaller drop. A study of the motion of the drop surface along the axis of symmetry brings to light that the two ends of the consolidated drop begin to oscillate with a phase shift for collision of non-equal-size drops. At longer times, the drop approaches a spherical shape and the phase shift is removed. The period of oscillations for the combined drop is also measured and it is noted that the first oscillation period is invariably shorter than the second period, regardless of the Weber number. The effects of an initial internal circulation within the drops are also demonstrated and the differences in shape evolution and period of oscillations are highlighted. These effects are primarily contained within the first few oscillations. The longitudinal component is primarily responsible for the coalescence and subsequent deformation of the combined drop into a plate shape as in the head-on collision models. As long as the bridge connecting the drops expands in radius, the transverse component results in the bulk of the coalescing drops to slide in opposite directions in an attempt to break the bridging between them.

It is demonstrated that the period of the first oscillation based on the maximum value obtain for the parameter which is defined as the ratio of the surface location along the X direction. The variation of the first period with Weber number in head-on and off-center collision of two equal-size drops is plotted in Figure 10. The circle points represent the variation of the first period with Weber number in head-on collision of two equal-size drops, the square points indicate off-center collision of two equal-size drops. Two minimas are on the curve, one at We = 0.5 and the other at We = 10. The latter is explained by considering the slow-downs occurred in the motion of the surface point along the X direction. Therefore, this part of the surface tends to decelerate the reverse flow rather than to accelerate that which is the case with higher Reynolds numbers. The larger period of oscillation at We = 1 is due to the increase of the viscous effects as the Weber number is decreased. Plotted in Figure 10, it is also noted that variation of the first period of oscillation in off-center collision of two equal-size drops as function of We are higher than variation of the first period of oscillation in head-on collision of two equal-size drops in the same Weber number.

Figure 10.

Variation of the first period of oscillation in head-on and off-center collision of two equal-size drops as function of Weber number.

It is manifested that dissipation of rotation after coalescence in off-center collisions is closely associated with the rate of dissipation of the longitudinal motion. Although the present simulations apply to the coalescence of two infinitely long cylinders, the predicted shape evolution during and after coalescence strongly ensembles the edge-on images envisaged in real experiments involving the collision and coalescence of spherical drops. It is also observed that variation of the first period of oscillation in off-center collision of two equal-size drops as function of Weber number are higher than variation of the first period of oscillation in head-on collision of two equal-size drops in the same Weber number.

In Figure 11, the surface tension energy (SE), the kinetic energy (KE), and the total energy (TE) of the drops from Figures 4 and 8 are plotted versus time. Initially, the force accelerating the drops together makes the kinetic energy to increase. After the force has been turned off, the drops move a short distance before colliding. Since the ambient fluid owns a viscosity, kinetic energy is dissipated by friction and the drops slow down. As the drops come in contact, the kinetic energy of the small-impact-parameter drops decreases rapidly, but the large-impact parameter drops are not affected to any significant degree. Due to the deformed drops the surface tension energy of the small-impact parameter drops increases. Since the drops remain almost spherical, the surface tension energy of the other drops hardly increases at all. When the film between the drops is obtained, part of the drops’ surface is removed and the surface energy reduced. This reduction is larger for the small-impact-parameter drops, since the area removed is larger. It is noted that the filament between the drops starts to neck down, the increase in surface area stops and the kinetic energy levels off. For the small-impact parameter drops, the rupture occurs near the time of maximum deformation, and surface energy is initially transformed into kinetic energy as the drop accommodates the new shape.

Figure 11.

Energies for the drops in Figures 4 and 8: the total energy (T.E red line), the surface energy (S.E blue line), and the kinetic energy (K.E black line) of the drops are plotted versus time.

Conclusions

The two-dimensional numerical simulations of the head-on and off-center binary collision of van der Waals liquid drops are implemented using LBM. The major finding can be summarized as follows:

The coalescence process of drops colliding is freely captured, which reflects invariably the phenomena of two droplet collision regardless of tracking equations and drops model. Numerical results of LBM are demonstrated to qualitatively agree with experimental data in the same Weber number. It is also observed that an increase of the Weber number results in increasingly thinner drops during the coalescence process in collision of equal-size drops. When collision of non-equal-size drops is considered, the early stages of coalescence involve considerably larger deformations for the smaller drop. It is further noted that the first oscillation period is invariably shorter than the second period, regardless of the Weber number. Variation of the first period of oscillation in off-center collision of two equal-size drops as function of Weber number are higher than variation of the first period of oscillation in head-on collision of two equal-size drops in the same Weber number. Clearly, the possibility of finite Weber number, arbitrary surface shape, finite velocity of approach, and the inclusion of another fluid outside the drops lends a tremendous richness to the class of singularities studied here.