A Numerical Investigation on the Collision Behavior of Unequal-Sized Micro-Nano Droplets

Micro-nano droplet collisions are fundamental phenomena in the applications of nanocoating, nano spray, and microfluidics. Detailed investigations of the process of the droplet collisions under higher Weber are still lacking when compared with previous research studies under a low Weber number below 120. Collision dynamics of unequal-sized micro-nano droplets are simulated by a coupled level-set and volume of fluid (CLSVOF) method with adaptive mesh refinement (AMR). The effects of the size ratio (from 0.25 to 0.75) and different initial collision velocities on the head-on collision process of two unequal-sized droplets at We = 210 are studied. Complex droplets will form the filament structure and break up with satellite droplets under higher Weber. The filament structure is easier to disengage from the complex droplet as the size ratio increases. The surface energy converting from kinetic energy increases with the size ratio, which promotes a better spreading effect. When two droplets keep the constant relative velocity, the motion tendency of the droplets after the collision is mainly dominated by the large droplet. On one hand, compared with binary equal-sized droplet collisions, a hole-like structure can be observed more clearly since the initial velocity of a large droplet decreases in the deformation process of binary unequal-sized droplets. On the other hand, the rim spreads outward as the initial velocity of the larger droplet increases, which leads to its thickening.


Introduction
The collision dynamics of binary micro-scale liquid droplets play an important role in the development of nanocoating [1] and nano-spray technology [2][3][4]. Micro-nano droplet collisions commonly occur in microfluidics, which is widely used in bioscience, clinical medicine, and environmental science [5][6][7]. Prominent examples are the fuel-spraying process in an engine combustor and the efficient distribution of agricultural sprays and powder production processes in food and life sciences. The frequency and outcomes of droplet collisions in the combustion chamber will determine the size and spatial distribution of the droplet group, and, ultimately, affect the heat and mass transfer effect of the system [8]. The results of droplet collisions are widely used to improve product performance and industrial efficiency. For example, entrained droplets travel through the steam-water separator of the steam generator in a nuclear power station. Based on the dynamic behavior and the discipline of droplet collision in the steam, the results can be referred during the design and optimization of separation apparatuses since considering the factors affects the separation efficiency [9].
Scenarios of microscale and nanoscale droplet collision dynamics have been widely investigated by either experiments or numerical simulations. Nahla et al. [10] demonstrated a new approach for droplet coalescence in microfluidic channels based on selective surface energy alteration, which focused on the coalescence process of the micro-nano droplet. Lian et al. [11] experimentally investigated the bouncing dynamics of impact droplets on the superhydrophobic surfaces for shortening contact time, and established the relationship between Weber number and contact time to find the minimum Weber number that formed the pancake-like bouncing behavior. Tang et al. [12] developed a theoretical model based on energy balance to explain the size-ratio dependence of the Weber number for the unequal-sized micro-nano droplets collision. Experimental results showed that the transition Weber number between bouncing and permanent coalescence regimes weakly depended on the size ratio, while the transition Weber number that separated permanent coalescence and separation regimes increased significantly as the size-ratio increased. However, for binary unequal-sized droplets, they concentrated on the mode division of the collision complex, the boundary lines between collision outcomes, and the energy evolution under low Weber number. Therefore, the numerical simulation is used to analyze the collision process of micro-nano droplets such as the flow and the velocity fields at a relatively high Weber number.
For numerical investigations, Bardia et al. [13] used the volume-of-fluid method to examine the binary nanometer-sized droplets collision. It was shown that the collision process can be described as follows: a pre-impact phase that ends with the initial contact of both droplets, and a post-impact phase characterized by the merging, deformation, and coalescence of the droplets. Kim et al. [14] investigated the collision of polymer nanodroplets by molecular dynamics simulations. The results showed that initial kinetic energy was changed into potential energy of the droplets as elastic energy when they began to contact each other. Then, the kinetic energy was completely changed into the potential energy of the droplets when droplets reached maximum deformation. Nikolopoulos et al. [15] studied the effect of the Weber number on binary unequal-sized droplets' collision at the size ratios of 0.5 and 0.6. It was found that the small droplet spread on the large one's surface at a high Weber number, while the small droplet penetrated into the large one at a low Weber number. Kuan et al. [16] developed a parallel, adaptive Eulerian-Lagrangian interface-tracking method to investigate head-on collisions of water droplets at high Weber numbers. Their simulations found that the rim grew with bread-like structures, and the extent of the retraction and deformation on the rim becomes weaker as the Weber number increases. Nikolopoulos et al. [17] employed the volume of the fluid (VOF) method to investigate the effects of the size ratio on binary droplets' collision in the reflexive regime. The collision phenomenon including the formation of filaments, the maximum deformation, the penetration of one droplet into the other, and the satellite droplet formation of binary droplets were discussed at the Weber number of 56.
In the above studies, investigations have been conducted by discussing the collision behavior of droplet, the effect of the size-ratio on droplet collision, and energy change for binary unequal-sized droplets at a low Weber number, and the collision dynamics of unequal-sized micro-nano droplets under the high Weber number (larger than 200) are important in many engineering operations even though they are not detailed now. The rim of the collision complex retracts and aggregates strongly due to surface tension at the Weber number range from 210 to 270, which causes complex deformation of the droplet. In this paper, deformation processes, the flow field, and the energy evolution of binary unequal-sized droplet collisions under a higher Weber number will be discussed in detail. The mode of the droplets' breakup and the evolution of the filament structure with the trend of the rim expansion will be analyzed. It is also expected that this study will help identify the collision dynamics of binary unequal-sized droplets at We = 210. This paper focuses on the application of the CLSVOF methodology for the numerical investigation of the collision behavior of unequal-sized micro-nano Droplets. This method can capture the complex evolutionary process of the interface with an accurate normal vector and good quality conservation, which has been applied successfully in the past work [18][19][20]. At the same time, adaptive mesh refinement (AMR) will also be used to reduce the computational expenditure.
The rest of the paper is organized as follows. Section 2 gives the introduction of the CLSVOF method, numerical model, and comparison with experimental results [21] at We = 210. In Section 3, the effects of the size ratio (from 0.25 to 0.75) and different initial collision velocities on the collision process of two unequal-sized droplets are studied in terms of the droplets' deformation, the velocity fields, and streamlines as well as the energy evolution. The concluding remarks are given in Section 4.

Mathematical Formulation
Coupled Level-Set and Volume of Fluid method (CLSVOF) combines the advantages of the Volume of Fluid (VOF) method in mass conservation with the advantages of the Level-Set method in the accuracy of interface curvature. The main processes of tracking the two-phase fluid interface include initialization of the phase function, solving of the phase function governing the equation, solving of the phase function transport equation, and re-initialization of the level set function. In the present study, the CLSVOF is used to study the collision of the binary droplets, which is modeled as incompressible, isothermal, immiscible flow. The governing equations for continuity, momentum, and phase-fraction equations are shown below.
where U is the fluid velocity, P is the pressure, F σ is the volumetric surface tension force, ρ is the mixture density and τ is the stress tensor defined as τ = −µ(∇U + ∇U T ), in which µ is the mixture viscosity. The mixture density ρ and viscosity µ adopting phase fraction α as a weight coefficient can be calculated as the transition of two-phase physical properties in the interface region.
In the original VOF function (Equation (3)), the additional compression term ∇ × (U c α(1 − α) is introduced to sharpen the interface [22]. An improved VOF equation can be written as: where U c is the compression velocity to suppress diffusion of the interface.
In this case, C α is the compression coefficient, which determines the compressive magnitude. It is generally greater than or equal to 1. The term ∇α/|∇α| is introduced to represent the convection of volume fraction function normal to the interface. The term α(1 − α) can be used to ensure itself invalid outside of the interface area and the divergence terms ensure mass conservation in terms of the entire compression term.
Furthermore, the Level-Set (LS) function ψ is introduced to ensure the interface smoothness by calculating a smoother curvature. The phase-fraction α is solved by transport Equation (3) for the VOF function, and use Equation (8) as the prime guess value to initialize the LS field.
where Γ is a small non-dimensional number, which depends on the minimum mesh size ∆x around the interface, the calculation formula is as follows.
The LS function ψ 0 can be re-distanced by a reinitialized equation.
where τ is the artificial time step. S(ψ 0 ) is a sign function defined as S(ψ 0 ) = ψ 0 /|ψ 0 |. The solution of ψ converges to |∇ψ|= 1, and the interface position is defined as the contour-line of ψ = 0. The number of iterations (ψ corr ) meets the following term.
where ε is the interface thickness. The Level-Set function can be used to calculate the interface normal by n = ∇ψ/|∇ψ| accurately. The smoother curvature is derived by k = ∇ × n. Therefore, the volumetric surface tension force can be calculated below.
Where ς is the surface tension coefficient, and δ is the Dirac function used to confine the influence of the surface tension to a narrow area around the interface.

Numerical Model
The two spherical droplets are formed in a cubic computational domain, as shown in Figure 1. Pressure-outlet boundary conditions are applied to all the computational domain planes. D S and D L present the diameter of the small droplet and the large droplet, respectively. The distance (L) between two droplets is 0.5 mm. The initial impact velocities of the two droplets are V S and V L , and the relative velocity is defined as V r = V L − V S . In the present study, the computational domain size is 6.25 mm × 6.25 mm × 6.25 mm. The grid number is 50 × 50 × 50. In the present study, the adaptive mesh refinement (AMR) is used, as shown in Figure 2. The minimum cell size ∆x is defined as L D /2 Lmax , where L D is the side length of a square, and L max is the cell refinement level. The numerical resolution related to the cell refinement level of the droplet is computed as D 0 /∆x, where D 0 is the initial diameter of the droplet. The numerical simulation method and the grid independence are verified by Qian et al. [23]. We take D S as the characteristic length and the relative velocity V r of the impacting droplets as the characteristic velocity. Four dimensionless numbers are employed to describe the collision dynamics of binary droplets in this work: the Weber number (We), the Reynolds number (Re), the size ratio (∆), and the dimensionless time (T). The Weber number (We) represents the ratio between inertial forces and the surface tension coefficient. The Reynolds number (Re) is defined as the ratio of the inertial and viscous forces. The size ratio (∆) is the diameter ratio of the small droplet and the large droplet.
We =ρ l V 2 r D S /σ (12) where ρ l is the liquid density, σ is the surface tension coefficient, and µ l is the droplet fluid viscosity. In this paper, the liquid density ρ l is 1000 kg/m 3 , the surface tension σ is 0.073 N/m, and the viscosity µ is 0.001 Nm −2 × s −1 . t 0 is the initial contact time of two droplets. The parameters used in our simulation in this paper are shown in Table 1.

Comparison with Experiments
The evolution of the droplet collision deformation for the case of We = 210 is presented in Figure 3, which is compared with the experimental work of Pan et al. [21]. It can be observed that a circular sheet forms in the middle of the collision complex after the coalescence of the droplets. The circular sheet is radially extended within the collision plane and forms a circular rim at the fringe of the circular sheet due to the end-pinching effect of the surface force, which includes Taylor-Culick rim (TC rim) [24][25][26]. Compared to the simulation results, the process of the deformation matches well with the experimental work, and this numerical method can be used to investigate the binary unequal-sized droplets collision at We = 210. However, the deviation between experimental and simulated results is the extent of the retraction and deformation on the rim. The deformation of the rim was affected by the anti-symmetrical disturbance of the RT (Rayleigh-Taylor) instability and the symmetrical disturbance controlled by the RP (Rayleigh-Plateau) instability in experiments. To simplify the model, the initial disturbance that is equivalent to instability is not imposed in the numerical simulation.

The Effects of the Size Ratio (∆)
To study the effect of the size-ratio on binary droplets' collision, different size-ratios of 0.75, 0.50, and 0.25 (case A-C in Table 1) are considered under the same Weber number of 210.
When the size-ratio is 0.75 in Figure 4a, the small droplet spreads on the surface of the large one and flattens gradually at the contact surface as the rim expands outward radially. The circular rim expands outwards and ruptures with the circular sheet, evolving into the filament and satellite droplets around the complex droplet. Figure 4b shows that the small droplet attaches to the big droplet and penetrates deeply into the big one because of a pressure difference between them, which is in agreement with the work of Nikolopoulos et al. [15]. Then the complex droplet that forms a hole-like structure evolves into the band-like filament and breaks into satellite droplets at the later stage of collision. When the small droplet becomes smaller in Figure 4c, due to the low impact momentum of the small droplet, the small droplet spreads on the surface of the large one with shallower penetration. Then the rim evolves into four distinct petal-like structures, which forms the ellipse holes on the fringe of petal-like structures. Difference size-ratio collisions of unequal-sized mico-nano droplets show similar collision behavior that two droplets coalesce rapidly after the collision due to intermolecular forces and then deform on the axial and radial direction. However, the tendencies of the rim expanding during the collision processes are different for various size ratios.  Figure 5 shows the evolution of the axial and radial deformations of the droplets, d X is the maximum radial deformation, and d y is the axial deformation in which the dimensions are characteristic of the length of equivalent volume D eq , which is defined as D eq = D 3 S + D 3 L 1/3 . T = 0 means the initial contact time of two droplets, as mentioned in Equation (17). It can be observed that the larger the size ratio, the smaller the axial deformation of droplets, and the larger the droplet spreads radially outward. In order to further illustrate the effect of the size ratio on the droplets' deformation, the streamlines and the velocity magnitude from the mixture at the cross-section Z = 0.003125 are given in Figure 6. It shows that, as the droplets move toward the collision center, gas is squeezed out, a gas jet sheet is formed between two droplets, and the streamlines form two opposite pairs of rotating vortexes at T = 0.00. As shown in Figure 6a, for binary equal-sized droplets collision, two opposite pairs of rotating vortexes are produced and attached to both sides of the rim. During the collision process, the droplets merge at a significantly flattened surface, and the rim expands outward with higher radial velocity while the axial velocity is dissipated gradually, which can be obtained in Figure 5. For binary unequal-sized droplets collision, as shown in Figure 6b-d, two different opposite pairs of rotating vortexes are produced due to the size ratio. When the size ratio is 0.75, when compared with binary equal-sized droplets collision, under the action of the opposite pair of rotating vortexes generated on the contact plane, a shallower hole-like structure is visible at T = 3.63 in Figure 6b. When the size ratio is 0.5 and 0.25, at T = 0.59, due to the detachment of the head of the main vortex by forming the opposite pair of rotating vortexes on the interface. The concave interface forms at the outer edge of the contact surface and is gradually stretched to the concaved crater at T = 1.62. The same stage of the formation of the main vortex by the vorticity generation from interfacial deformation was also observed and explained by Sun et al. [27]. Under the influence of collision dynamics, the vortex rings produced at the contact surface are dissipated at T = 3.63, and the hole-like structure is more clear, which is also shown in Figure 6c,d.

The Effects of the Initial Velocity
In order to illustrate the effect of the large droplet and the small droplet on the collision dynamics, the initial velocity of the droplets is changed under the condition of We = 210 and ∆ = 0.5, and the vector fields inside and around the colliding droplets are also further given, as shown in Figures 7  and 8. Figure 7 indicates that the momentum of the small and large droplets has different effects on the motion tendency of the mixing droplet and the expansion of the rim at We = 210. When the larger droplet impacts the smaller droplet with the velocity of 5.54 m/s, the mixing droplet integrally moves upward, forming the hole-like structure. Subsequently, the larger droplet continues to move and expand outward, which leads to a thicker rim structure. When the velocity of the smaller droplet is −5.54 m/s, the smaller droplet spreads only on the surface of the larger droplet and the motion of the mixing droplet remains almost unchanged. The mixing extent of the collision droplet is more distinct and the axial deformation of the collision droplet is larger. In Figure 8, the vector fields inside and around the colliding droplets are given to explain the flow direction of the small droplet and the large droplet. It is found that, when the droplets collide with the opposite initial velocity, the mass of the small droplet tends to flow into the large one, and spreads outward from the inside of the large droplet. In Figure 8b, when the velocity of the large droplet is 5.54 m/s, the large droplet with the large momentum deforms and carries the mass of the small droplet to spread outward widely. While the small droplet takes the velocity of −5.54 m/s, as shown in Figure 8c, a little mass of the small droplet spreads on the surface of the large drop with a shallower penetration. It is shown that, during the collision process, the momentum of the large droplet has a larger effect on the formation of the hole-like structure and the vector fields inside and around the colliding droplets also indicate that, for binary unequal-sized droplet collisions, the mass of the small droplet tends to flow into the large one after the collision, so that the colliding droplets are difficult to separate.

Evolution of Dimensionless Energy
In this section, following the research studies of Yoshino et al. [28] and Sarrokaet al. [29], the energy transfer processes are also discussed.
The initial kinetic and surface energy of the droplets at T = 0 are computed as follows.
The kinetic and the surface energy at each time step during the collision process are computed as: where SA p is the polygon area and υ cell is the volume of the computational cell, f ijk is the volume fraction of the liquid phase, u ijk , υ ijk , and w ijk is the velocity component and N cell is the number of the computational cell. The kinetic energy and surface energy are non-dimensionalized by the initial kinetic energy and surface energy, respectively. The evolutions of dimensionless energy of binary droplets under different size ratios and initial impact velocity are given in Figure 9.  Figure 9a,b shows that the larger size-ratio is, the larger the surface energy is, and the smaller the kinetic energy is. Corresponding to the maximum deformation process of droplets in Figure 4, it can be seen before deformation that the surface energy and the kinetic energy remain unchanged. As binary droplets move toward each other, the kinetic energy is transformed to the surface energy. The surface energy increases gradually with the kinetic energy decreasing. The larger the size ratio is, the more kinetic energy is converted to surface energy.
In Figure 9c,d, the energy evolution of the droplets with different initial velocity shows that the kinetic energy of the droplets with the opposite initial velocity is the largest, and the kinetic energy of the large droplet with the velocity of 5.54 m/s is larger than that of the small droplet with the velocity of −5.54 m/s. Compared with the other two cases, when the velocity of the small droplet is −5.54 m/s, it can be observed that, during the collision process, more kinetic energy is converted to surface energy. For the surface energy evolution of the droplets, it can be found that, when the velocity of the large droplet is 5.54 m/s, the surface energy is the largest. On one hand, the surface of the mixing droplet contracts as much as possible due to the surface force. On the other hand, the surface energy is increased by thickening the rim. Hence, the collision complex is stable without satellite droplets forming around the mixing droplet. The surface energy of the droplets with the opposite initial velocity is larger than that of the small droplet with the velocity of −5.54 m/s in the late stage of collision. Corresponding to the droplet deformation process, like Figure 7, it can be explained that, under the relatively high Weber number, the impact momentum of the small droplet has a great effect on the mass of the outward spreading rim, and the mass of the outward spreading rim decreases as the impact momentum of the small droplet increases. The collision dynamics of the droplets are mainly dominated by the large droplet, which had also been provided by Mohammad et al. [30] at a low Weber number.

Conclusions
Binary unequal-sized micro-nano droplets collision has been numerically simulated by a coupled level-set and volume of fluid (CLSVOF) method with adaptive mesh refinement (AMR). The effects of the size ratio and different initial collision velocities on the collision process of the droplets are discussed and summarized as follows.

1.
For binary unequal-sized micro-nano droplet collisions, the maximum axial deformation of the complex droplet increases when the difference of the diameter between two droplets becomes smaller, and the complex droplet reaches its maximum value of radial deformation when ∆ = 1.

2.
The spreading effect of the collision complex becomes more clear as the size-ratio increases and conversely facilitates the mixing effect of the collision complex. As for the atomization combustion, the complex droplets contact the gaseous medium sufficiently when the collision droplets are equal-sized, which could promote the efficiency of combustion.

3.
When two droplets keep the constant relative velocity, the motion tendency of the droplets after the collision is mainly dominated by the large droplet. As for a steam separator, the velocity direction of complex droplets is maintained by increasing the velocity of large droplets due to the large momentum of large droplets. The complex droplets are hardly affected by the steam flow and much easier to remove by the wave-type vanes as a result.

4.
The mixing effect can be improved by suppressing the initial velocity of the larger droplet. Moreover, it is difficult to separate a droplet from another after the collision since the small droplets tend to penetrate into the large one.