Dynamics of simultaneously impinging drops on a dry surface: Role of inhomogeneous wettability and impact shape

Journal of Colloid and Interface Science 516 (2018) 232–247

Contents lists available at ScienceDirect

Journal of Colloid and Interface Science journal homepage: www.elsevier.com/locate/jcis

Regular Article

Dynamics of simultaneously impinging drops on a dry surface: Role of inhomogeneous wettability and impact shape K. Ashoke Raman Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 117576, Singapore

g r a p h i c a l a b s t r a c t Role of inhomogeneous wettability and impact shape on the interaction dynamics between two simultaneously impinging drops onto a dry surface.

Role of inhomogeneous wettability and impact shape on the interaction dynamics between two simultaneously impinging drops onto a dry surface. Contact angle Hysteresis

Wettability Gradient

Impact Shape

a r t i c l e

i n f o

Article history: Received 27 October 2017 Revised 16 January 2018 Accepted 17 January 2018

Keywords: Contact angle hysteresis Wettability gradient Lattice Boltzmann method Impact shape

E-mail address: [email protected] https://doi.org/10.1016/j.jcis.2018.01.063 0021-9797/Ó 2018 Elsevier Inc. All rights reserved.

a b s t r a c t Hypothesis: The quality of the printed lines in applications such as ink-jet printing and additive manufacturing is affected by the interactions between the impinging drops. Impact shape and the inhomogeneity in surface wettability govern the spreading and recoiling dynamics of the interacting drops. Hence, understanding the role of these factors on the interaction dynamics is essential to optimize these applications. Numerical experiments: Phase-field based lattice Boltzmann method solver has been employed to investigate the interaction dynamics of two simultaneously impinging drops onto a dry surface. A geometrybased contact angle scheme is used to model the moving contact line. Findings: Numerical simulations reveal that the previously identified interaction modes (Raman et al., 2017) are sensitive to the contact angle hysteresis, resulting in different impact outcomes. Two different interaction mechanisms have been discerned when drops impinge on a surface with a wettability gradient. It is shown that the deviation from the spherical geometry of the impact shape leads to different spreading behaviors and droplet morphology around the connecting region. With the increase in the cross-sectional aspect ratio, the interaction dynamics of oblate
oblate combination is similar to its spherical counterpart, albeit at a faster recoiling rate. Ó 2018 Elsevier Inc. All rights reserved.

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

1. Introduction When a liquid droplet falling onto a dry surface interacts with another wall impinging drop, the resulting spreading and recoiling dynamics of the system is considerably different from those observed from a single droplet impact scenarios. The ensuing morphological structures and hydrodynamics are much more complex. An uprising sheet is observed when two drops with high Reynolds number and Weber number impinge simultaneously on a solid surface [1]. Column-like morphological structures are observed when drops impact successively at low-velocity ratios [2]. Such flows involving these complex interactions are ubiquitous in nature and engineering applications, most notably in fuel spray injection, spray cooling, inkjet printing, rapid prototyping, touchless cleaning and many others. In spray cooling systems, the overall heat transfer rate between the droplet and the substrate is influenced by these interactions. Similarly, the geometric precision and uniformity of printed lines are dependent on the hydrodynamics of droplet interactions [3]. In general, the nature of these dropdrop interactions is subjected to several working conditions of the system. These conditions include both the physical factors such as impact velocities, droplet diameter, smoothness and rigidity of the target surface as well as chemical factors like viscosity, surface tension, density and chemical nature of the substrate. Therefore, a fundamental understanding of the effects of these physical and chemical factors on drop-drop interactions is essential for finding out optimal working conditions of these applications. Owing to its importance in several industrial application and occurrence in varied natural phenomenon, the impact of a single droplet on wet and dry surfaces have been investigated by several researchers [4,5]. Important dynamical aspects of the impact process of a single droplet on thin films and drywall are presented in the review articles provided by Yarin [6] and Rein [7]. Coalescence dynamics between two sessile drops has been studied by several researchers [8–10] to investigate the growth and expansion of the liquid bridge. The neck growth is always viscously dominated at initial times. As the neck grows in size and the Reynolds number increases, the evolution process becomes inertia driven. It is found that sessile drops from different but completely miscible do not always undergo instantaneous coalescence upon contact [11,12]. They go through a temporary phase of non-coalescence by being connected through a liquid bridge. Tong et al. [13] performed a numerical study on successive droplet impingement on a solid substrate and concluded that the contact angle plays a key role in the time of droplet interaction, spreading and recoiling characteristics of the droplet system. They observed that a dynamic contact angle model gives a better agreement with experiments. The interaction between drops is dependent on the radius of curvature, which in turn is governed by the contact angle subtended by the interface with the substrate. Unlike ideal smooth surfaces which are characterized by an equilibrium contact angle, the contact angle of real surfaces employed in industrial applications exhibit inhomogeneities in wettability. The contact angle on such surfaces varies between the advancing and receding contact angles [14]. In this range, also known as the contact angle hysteresis, the contact line appears to be immobile. The effect of contact angle hysteresis on the droplet coalescence and mixing was investigated by Nilsson and Rothstein [15]. The droplet interactions were characterized with large deformations in case of low contact angle hysteresis. However, with an increase in contact angle hysteresis, the maximum deformation was observed to be marginally above the analytical radius obtained from two drop coalescence. Lee and Son [16] found that an increase in the receding contact angle leads to shrinkage in the droplet-wall contact area and retraction. Higher advancing contact angle was found to be favorable for printing a narrow line. For successively impinging

233

droplets with different impact velocities, it is shown that an increase in the advancing contact leads to greater recoiling of the droplet system in-phase coalescence mode [2]. Increase in receding contact angle is found to promote the interaction time for droplets impacting in out-of-phase mode. Effect of advancing contact angle on the interactions multiple droplet deposition has been investigated by Zhang et al. [17]. They proposed an optimal overlap ratio to obtain uniform line formation. Apart from contact angle hysteresis, nonuniformities in the microstructure patterns and chemical composition would lead to the formation of a wettability gradient on the solid surface. Such surfaces with gradients in wettability arising due to chemical irregularities have been discerned to induce a directional motion when it interacts with a single droplet [18]. Owing to the applicability of directional motion in dropletbased devices, several experimental and numerical investigations [19–22] have been conducted and the motion of a single droplet on substrate driven by wettability gradient is well understood. However, there are limited studies in the literature investigating the role of wettability gradient on drop-drop interactions. Gong and Cheng [23] performed numerical investigations employing lattice Boltzmann model to study droplet motion and coalescence on a surface with wettability gradient. Using a finite element based phase field formulation, Ahmadlouydarab and Feng [24] explored the motion and coalescence of sessile drops driven by wettability gradient with an external flow. They observed that drop coalescence is favored by strong wettability gradient while external flow assisted separation. It has also been shown that continuation of flow field inside the droplet system promotes mixing when a single drop impacts onto a sessile droplet resisting on a surface with wettability gradient [25]. Interactions between wall impinging drops are not only influenced by the substrate characteristics and chemical properties of the droplet liquid but also governed by the physical aspects of the droplet system. For large inter-droplet spacing between drops impinging successively onto a solid surface, it is observed that the liquid in the second droplet swells up after the collision with the first one [26]. Similarly, it is found that the quality of the printed lines formed due to droplet deposition is sensitive to the center-to-center distance between the drops [27]. The lines would either remain continuous or rupture depending on the droplet spacing. Among the various thermal management technologies employed for the proper functioning of electronic devices, spray cooling is thought to be most appropriate due it’s high heat removal capability and uniformity. Owing to the inherent complexity of the sprays incoherent droplet impact, random trajectories arising from multiple interactions before impact and the large number of physical variables involved, fundamental understanding of the physical mechanisms is a challenging task to perform. Hence, in order to have a better understanding of the physical mechanisms, it is important to isolate the number of variables by employing multiple streams of mono-dispersed droplets. It has been observed that the transport phenomenon and heat transfer characteristics are not only influenced by the rate of droplet deposition, but also from the interactions between these streams of multiple periodic droplet impingements [28–31]. Horizontal droplet spacing between adjacent droplets is found to be an important factor in reducing drop-drop collision and increasing heat transfer in case of multiple droplet train impingement cooling [28,29]. Simulations performed by Batzdorf et al. [31] predicted a decrease in heat transfer characteristics due to interaction and coalescence of the droplets on the solid surface. Zhang et al. [30] performed experimental investigations to study the hydrodynamic and heat transfer characteristics of multiple droplet trains impinging on a pre-wetted solid surface. They concluded that the horizontal spacing and impingement pattern play significant roles in cooling performance. For triple droplet

234

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

train impingement, an optimal spacing was observed. However, unlike the case for double droplet arrangement, better cooling performance is observed at lower droplet spacing for hexagonal impingement pattern. It is clear from the literature that interaction modes between the droplet trains impinging on the surface and the resulting morphological structures of the droplet system have a significant effect on the transport phenomenon. In order to understand the physical mechanism involved in the interactions of double droplet train impingement, researchers have performed investigations on two simultaneously impinging drops [32– 34,31]. Two-dimensional numerical investigations performed by Wu et al. [32] indicate the generation of sparse segments when two droplets with different sizes impinge simultaneously on a dry surface. They attributed this observation to the formation of an oblique jet. When two drops with different impact velocities impinge simultaneously on a neutral wetting solid surface, two different coalescence modes, depending on the motion of the contact lines at the time of the collision, were identified by Raman et al. [33]. In addition, they showed that owing to the drag induced on the droplet system, the density of the surrounding medium would influence the interaction and coalescence outcome. Another important physical factor which would influence droplet interactions is the shape of the droplets at the time of impact. Falling drops usually oscillate between oblate to prolate shapes before impacting onto a surface. Studies on investigating the role of impact shape have been primarily limited to single droplet impact scenarios [35–37]. The present body of work aims to answer two interesting questions about the interaction dynamics when two drops impinge simultaneously on a solid surface: what is the role of wetting inhomogeneity in the form of contact angle hysteresis on the interaction modes identified by Raman et al. [33]?Subsequently, we investigate the role of chemically induced wettability gradient on droplet interaction and recoiling dynamics. Secondly, how does the geometric shape of the impacting drops affect the interaction dynamics? While the first question accounts for the chemical aspects of the droplet system, the second question probes into the physical factor accounting for the role which nonspherical drops would have on the interaction dynamics. To the best of our knowledge, the effect of chemically induced wall wettability gradient and the role of impact droplet geometry on the interaction dynamics of simultaneously impinging drops has not been investigated in prior studies. We set out in filling these gaps in knowledge through a clear analysis of the interaction mechanisms by investigating various scenarios of morphological change around the connecting region and the transitions between surface energy and droplet inertia. These questions are numerically investigated by employing a three-dimensional framework based on the lattice Boltzmann method (LBM) [38,39]. The kinetic nature of LBM facilitates in simulating multiphase flows by incorporating intermolecular forces to model macroscopic surface tension effects. In addition, compared to conventional methods, LBM does not require the solving of Poisson equation and is an explicit scheme. Owing to its local nature for computations, it is easy to parallelize. This makes it an efficient tool for three-dimensional multiphase flow simulations. The remainder of the paper is organized as follows: Section 2 describes the system of equations governing the phase-field lattice Boltzmann model. The geometry-based contact angle model which also takes into account for the contact angle hysteresis is detailed in this section. The problem setup along with the details of the computational domain and the applied boundary conditions are outlined in Section 3. We then proceed with the presentation of the results and discussions in Section 4. Concluding remarks are provided in Section 5.

2. Numerical method In contrast to the conventional numerical methods solving the Navier-Stokes equation, the lattice Boltzmann method (LBM) [38,39] is a simplified and reduced mesoscopic kinetic model based on the Boltzmann equation and has been successfully used to perform efficient computations of multiphase flow problems [40,41]. With its basis in the kinetic gas theory, it relies on the evolution of a particle distribution function and the macroscopic quantities are obtained from the hydrodynamic moments of these distribution functions. The mesoscopic nature of the LBM has been particularly useful in simulating the interfacial flows. The most popular multiphase LBM models are the color gradient model [42,43], interparticle potential model [44], free-energy model [45] and index-function model [46]. Owing to its similarities to the lattice gas automata, the color method is the earliest of the multicomponent extension of LBM. While the free energy models are derived from a free-energy functional, thereby making them thermodynamically consistent, the thermodynamics of multiphase flows emerges from the interactions between the lattice fluids. The index-function models employ two distribution function and the corresponding LBEs recover the macroscopic Cahn-Hilliard interface tracking equation and Navier-Stokes equations. A comprehensive analysis of different multiphase LB models is presented in the book by Huang et al. [47]. To understand the interaction dynamics between simultaneously wall impinging drops, we employ the phase-field based lattice Boltzmann method [48,49,46,50]. Two particle distribution functions are used in this model [48] to recover incompressible Navier-Stokes equation (g a ) and a macro interface capturing phase-field equation (f a ). The model consists of stress and potential forms of intermolecular forcing terms in the momentum equation and the phase-field model for order parameter, respectively. Stable discretization schemes are considered to discretize the forcing terms in each collision step which help in improving numerical stability for high-density ratio cases. The authors [48] have employed the treatment given in [51] on the characteristics of stress and potential forms related to the pressure tensor. The discrete Boltzmann equation (DBE) for pressure and momentum, which is obtained from the following transformation




ga ¼ f a þ

 p
q Ca ð0Þ c2s

ð1Þ

is expressed as:

@g a @g ðg
g eq ðeai
ui Þ@ i ðqc2s Þ a Þ þ ðCa ðuÞ
Ca ð0ÞÞ þ eai a ¼
a k c2s @t @xi þ

ðeai
ui Þ½j@ i ð@ k q@ k qÞ
j@ j ð@ j q@ i qÞ Ca ðuÞ c2s

ð2Þ

where k is the relaxation time due to collision and ea is the microscopic particle velocity. The corresponding (DBE) for the order parameter is given as eq

@f a @f ðf
f a Þ þ eai a ¼
a k @t @xi þ

ðeai
ui Þ½@ i qc2s
q@ i ð/
j@ 2j qÞ Ca ðuÞ: c2s

ð3Þ

The equilibrium distribution functions are given by

"

eq

f a ¼ wa q 1 þ " g eq a ¼ wa

ea :u ðea :uÞ2 ðu:uÞ þ
c2s 2c2s 2c4s

#

# p qea :u qðea :uÞ2 qðu:uÞ þ þ
c2s c2s 2c2s 2c4s

ð4Þ ð5Þ

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

pffiffiffi where cs ¼ 1= 3 and wa is the corresponding integral weights for a D3Q19 lattice velocity model:

1 wa ¼ ; a ¼ 0 3 1 ; a 2 ½1; 6 wa ¼ 18 1 wa ¼ ; a 2 ½7; 18 36 and

"

Ca ðuÞ ¼ wa 1 þ

ð6Þ # ð7Þ

These two equations are discretized along the characteristics over a time step dt. The trapezoidal rule is employed for time integration in ½t; t þ dt which is coupled with the space integration in ðx þ ea dt; t þ dtÞ. The resulting lattice Boltzmann equations are solved in three steps: (i) Pre-streaming step

ga ðx; tÞ ¼ g a ðx; tÞ


þ

1 sat sat sat /  4bðq
qsat v Þðq
ql Þðq
2 ðqv þ ql ÞÞ

D¼

ð14Þ

  dt ðeai
ui Þ½j@ i ð@ k q@ k qÞ
j@ j ð@ j q@ i qÞ  C ðuÞ a  2 c2s ðx;tÞ ð8Þ

 f ðx; tÞ ¼ f ðx; tÞ
f a
f a  a a 2s ðx;tÞ

ð15Þ

r¼

sat 3 pffiffiffiffiffiffiffiffiffi ðqsat l
qv Þ 2jb: 6

ð16Þ

The density of the fluid q, hydrodynamics pressure p and the velocity u are calculated by taking the moments of the corresponding distribution function:

X f a

ð17Þ

X dt @ qc2s ga c2s þ ui 2 @xi a 

 X dt @ @q @q @ @q @q
g a ea þ j ui ¼ 2 @xi @xk @xk @xj @xi @xj a

ð19Þ

The relaxation parameter s is related to the kinematic viscosity m ¼ sc2s dt, which can be calculated by a linear interpolation

s ¼ C sl
ð1
CÞsv

C¼

ðx;tÞ

ð9Þ (ii) Streaming

gðx þ ea dt; t þ dtÞ ¼ ga ðx; tÞ f ðx þ e dt; t þ dtÞ ¼ f ðx; tÞ a a

ð18Þ

ð20Þ

where sl and sv are the relaxation times for liquid and vapour, respectively and the parameter C is the composition defined by

 2  dt ðeai
ui Þ½@ i qc2s
q@ i ð/
j@ j qÞ  C ðuÞ  a 2  2 cs

ð10Þ ð11Þ

(iii) Post-streaming step

  1  Þ ðga
geq a  2s þ 1 ðxþea dt;tþdtÞ

2s dt ðeai
ui Þ@ i ðqc2s Þ ðCa ðuÞ c2s 2s þ 1 2   2s dt
Ca ð0ÞÞ þ 2s þ 1 2

ðq
qsat v Þ sat Þ ðqsat
q v l

ðxþe dt;tþdtÞ

a   ðeai
ui Þ½j@ i ð@ k q@ k qÞ
j@ j ð@ j q@ i qÞ  C ðuÞ  a  c2s xþea dt;tþdt

ð12Þ f a ðx þ ea dt; t þ dtÞ ¼ f a ðx þ ea dt; t þ dtÞ  1 eq  ðf a
f a Þ
2s þ 1 ðxþea dt;tþdtÞ

ðxþea ;tþdtÞ

ð13Þ where s ¼ k=dt. Detailed discussion and derivations of the discrete Boltzmann equations and the intermolecular forcing terms presented above are outlined in [48,52,46,53,51].

ð21Þ

The mixed difference scheme and the second-order central difference scheme are considered for discretizing the forcing terms in the pre-streaming and post-streaming collision steps, respectively. Further details on the discretization schemes can be found in [48]. A geometry formulation contact angle proposed by Ding and Spelt [54] is employed to model the three phase contact line. This geometric scheme can also be used to model contact angle hysteresis and its applicability to droplet impact problems has been demonstrated [55]. The procedure essentially involves updating the values of the order parameter (density in the present model) to enforce the following equation:

n  $q ¼
tan

þ

 2  2s dt ðeai
ui Þ½@ i qc2s
q@ i ð/
j@ j qÞ  þ Ca ðuÞ  c2s 2s þ 1 2

j

2b

p¼

  dt ðeai
ui Þ@ i ðqc2s Þ  ð C ðuÞ
C ð0ÞÞ a a  2 2 cs ðx;tÞ

¼ ga ðx þ ea dt; t þ dtÞ


rffiffiffiffiffiffi

a

 g a
g eq a  2s ðx;tÞ

g a ðx þ ea dt; t þ dtÞ

4 sat ðqsat l
qv Þ

where j is a constant related to the magnitude of surface tension. The surface tension force r is represented as

q¼

eq 

þ

The chemical potential / is given by

sat where b is a constant, and qsat v and ql are the saturation densities of the vapour and liquid phases, respectively. The interface thickness, denoted as D, is given by

ea :u ðea :uÞ2 ðu:uÞ þ
c2s 2c2s 2c4s

þ

235

p 2


h j$q
ðn:$qÞnj

ð22Þ

From the above form, we can achieve the desired wettability between the solid and the fluid by specifying the desired contact angle. Once the density at the boundary points is specified, the normal gradient condition in Eq. (22) is satisfied in the solver. It is important to mention at this point that a static contact angle model is employed in the current work to investigate the droplet interaction dynamics. The static contact angle is the angle subtended by the interface with the wall when the contact line is stationary. This angle satisfies Young’s equation [56]. However, when the contact line moves, the angle subtended by the interface depends on the speed of the contact line, fluid properties, and the characteristic length scale [57]. In general, modeling moving contact line is a challenging problem as it involves a large range of length scales. For the slow motion of the contact line, the interface is divided into three regions [58]. The inner region follows a wedge pattern forming the microscopic contact angle. Unlike the apparent contact angle formed at the outer region, the microscopic contact

236

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

boundaries are considered to be periodic. In the present study, the droplet diameter Di and the impact velocity of droplet 2 U 2 are set to be the characteristic length and velocity, respectively. The non dimensional time is given by T  ¼ ðtU 2 =Di Þ, where t is the simulation time in lattice units. The important dimensionless parameters governing droplet impact on solid substrate include   (i) Density ratio q ¼ qql (ii) Viscosity ratio l ¼ lll (iii) Reyg g   q U2 D nolds number Re ¼ ql Ul2 Di (iv) Weber number We ¼ l r2 i and  (v) Velocity ratio V r ¼ UU12 . The density and viscosity ratios were

angle is not influenced by the large-scale flow geometry. Models relating the dependence of dynamic angle on contact line speed are available in literature [59,60]. A comprehensive review on various contact line models is discussed in Sui et al. [61]. Employing a dynamic contact angle model would result in greater accuracy in numerical predictions when compared to corresponding experimental results for wall impinging drops. However, to capture the range of length scales near the contact region, adaptive mesh refinements techniques are required [62]. The present work does not include such modeling techniques and leaves scope for future work.

set to be 841 and 52, respectively. In order to investigate the effects of nonuniform wettability and impact shape, all the simulations are performed at Re = 300 and We = 25. The center-to-center distance between the droplets is set to be DX ¼ 1:6 in the present work. The effect of contact angle hysteresis is considered only in Section 4.1 whereas the other sections have been investigated using the equilibrium contact angle. The developed 3D solver based on the present model has been calibrated to evaluate fluid-fluid and fluid-solid interactions for both static and dynamic scenarios. For brevity, we do not discuss them in this work. Static benchmark cases such as verification of Laplace law for a stationary bubble, evaluation of equilibrium contact angle for a droplet resting on a solid surface been discussed in our previous work [55]. The experimental comparison of the

3. Problem set-up Consider two side-by-side droplets of equal size with density ql and the dynamic viscosity ll impacting simultaneously on a solid substrate as shown in Fig. 1(a). Both droplets are of equal diameters (Di). The two droplets have been initialized with an impact velocity of U 1 and U 2 for droplet 1 (left) and droplet 2 (right), respectively. The computational domain with dimensions Lx  Ly  Lz where Lx ¼ 281 and Ly ¼ 181 are the lateral dimensions of the domain and Lz ¼ 91 is the dimension along the vertical direction. The drop diameter (Di) is set to be 64 lattice units. No-slip are imposed on the top and bottom boundaries and the side

(a)

U1

Z

ΔX

Y

U2

(b)

X

Droplet 1 Droplet 2

(c)

(d) 100

Simulation

80

Radius (μm)

Experiment [63] 60

40

20

0

0

50

100

150

Time (μS) Fig. 1. (a) Schematic of the computational model for the simultaneous two droplet impact problem. (b) Distribution of the contact angle on a surface characterized by wettability gradient. The surface is divided into three zones. (i) Zone A: Region of constant contact angle h1 = 120° (ii) Zone B: Region of constant wettability gradient. (iii) Zone C: Region of constant contact angle h2 . (c) Cross sections and geometric parameters of the Oblate and Prolate shaped droplets. (d) Validation of the LB solver against the experimental data of the two-droplet interaction situation from Fig.8 in [63].

237

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

dynamic case for a droplet impacting on a solid surface with different contact angles has been outlined in our preceding article on the interactions between simultaneously impinging drops [33]. To further validated the employed solver for multiple droplet interactions, a simulation was conducted to compare with the experimental results from Yang et al. [63]. In the experiments, a water droplet of diameter 60 lm impinged with an impact velocity of 0.77 m/s onto a sessile water drop resting on a stainless steel surface with h = 90°. As observed from Fig. 1(d), a good agreement is observed between the simulations and the experimental data. An initial dip in the contact radius is also observed our simulations which results due to the coalescence of the two interacting droplets. However, the deviations in the maximum contact radius and the evolution trend are attributed to the contact angle hysteresis exhibited in real surfaces. The application of a static contact angle model instead of a dynamic model also contributes to these deviations between experimental and numerical results.

as it expedites the expansion of the connecting neck as observed at T⁄ = 3.5. The expansion of the connecting neck results due to the pumping of the liquid from both drops into the neck region. This growth in neck width, arising from unbalanced Laplace pressure between the drop and the negatively curved neck, is resisted with the increase in hA . As mentioned previously, increase in hA leads to greater accumulation of droplet liquid inside the outer rims. This results in the increased pressure difference between the convex outer rims and the flat central region of the droplet, promoting vigorous recoiling of the droplet system as observed at T⁄ = 3.5. This upward motion of the droplet liquid along its central axis resists the rate of fluid inflow into the connecting neck and hinders neck expansion. To further quantify the upward recoiling flow field inside the droplet system, we monitor the non-axial kinetic energy coefficient [36,64] defined as follows:

hR

gxy ¼ hR 4. Results and discussion In the following subsections, we would present the numerical results and elucidate the physics of the drop-drop interactions. The role of inhomogeneity in surface wettability and impact shape of the impinging droplets on the interaction dynamics are systematically investigated. 4.1. Role of contact angle hysteresis We first begin our study on the interaction dynamics by exploring the role of surface inhomogeneities in the form of contact angle hysteresis. Two different modes of interactions were identified by Raman et al. [33] for drops impacting with a given velocity ratio. These modes were classified based on the propagating direction of the interacting contact lines. While a neutral wetting surface (h = 90°) is employed in [33], the rate at which the contact line advances and recedes depends on the advancing (hA ) and receding (hR ) contact angles, respectively. When both the contact lines are advancing, the interaction mode is termed as Out-of-Phase [33]. As such, we investigate the effect of hA on this mode of interaction. The temporal evolution of the impacting drops on solid surfaces having hA = 70° (top row) and hA = 120° (bottom row) is illustrated in Fig. 2(a). The droplets are impinging with a velocity ratio of 0:5 and the receding contact angle is fixed at hR = 60°. When a drop impacts onto a solid surface, the mean curvature of the advancing contact line increases as the contact angle subtended by it is increased. This leads to greater accumulation of the droplet liquid into the spreading rims. Hence, as observed at T⁄ = 1.5, the droplet impacting for hA = 120° is characterized with a thicker spreading rim when compared to that with hA = 70°. This, in turn, delays the contact time of the propagating rims and hence the formation of the connecting neck. To further elucidate the role of the advancing contact angle, the kinetics of the propagating contact lines are shown in Fig. 2(b). It shows the time evolution of the positions of the approaching contact lines corresponding to the left (X il ) and the right (X ir ) droplet. During the initial phase of the impact process, the evolution of X il and X ir remains fairly independent of hA as the impact process is inertia driven. As time proceeds, it is noticed that an increase in hA delays the connecting time. This is inferred from both Fig. 2(b) and the droplet profiles shown at T⁄ = 1.5. For the case with hA = 70°, the rims have already undergone coalescence leading to the formation of the connecting neck while two distinct drops are seen for hA = 120°. Lower advancing contact angle facilitates stronger capillary spreading. This behavior of the moving contact line supplemented by the inertia of the drop liquid inside the advancing rims leads to earlier connecting time as well

Vdrop Vdrop

ð1=2Þqðu2x þ u2y ÞdV

i

ð1=2Þqðu2x þ u2y þ u2z ÞdV

i

ð23Þ

The observations illustrated in Fig. 2(c) clearly indicate that the oscillations in the fraction of the total kinetic energy transferred along the wall normal direction increase with the increase in hA . g = 1 implies that the flow field inside the droplet is primarily parallel to the surface. This is observed from the initial sharp rise in the non-axial kinetic energy coefficient to g = 1 for all the cases considered, as the droplet system undergoes inertial spreading, reaches maximum spread and undergoes recoiling. For hA = 70°, we notice that after an initial variation, g reaches a nearly constant value close to 1 implying that the drop system begins to recede to attain its equilibrium configuration. Most of the initial kinetic energy for this case is converted into surface energy  R  Er ¼ V Eo þ jj$qj2 dV as the droplet system undergoes greater capillary spreading, where Eo 2

2

is approximated as

sat Eo  bðq
qsat v Þ ðq
ql Þ . This is inferred from the maximum value of Er for the case with hA ¼ 70 as compared to the other two cases as illustrated from the temporal evolution of surface energy shown in Fig. 2(d). When the velocity ratio is decreased, we observe that at the time of interaction, one droplet undergoes recoiling whereas the other is in its spreading phase. In such situations, the interactions between the impinging drops occur when both the contact edges propagate in the same direction. This mode of coalescence is termed as In-phase coalescence mode [33]. In contrast to the previous case where both the contact edges collide head-on, in this case, the advancing contact edge of droplet 2 has to outreach the receding contact edge of droplet 1 for the drops to coalesce. This recoiling rate of the contact line is governed by the receding contact angle (hR ), which in turn determines the interaction outcome. We investigate the role hR when the droplets impact with a velocity ratio of V r = 0.2. The advancing contact angle is set to be hA = 100°. Fig. 3 illustrates the temporal evolution of the impact process for hR = 50° (top) and 90°. As droplet 1 reaches its maximum spread and flow reversal develops along the direction of the upward central axis of symmetry, the droplet begins to recoil as observed at T⁄ = 2.25 for both the considered cases. However, owing to contact angle hysteresis, the contact line remains pinned to the surface. This duration (Dthyst ) of contact line pinning depends on the difference between hA and hR . For the case with hR = 50°, the retracting contact edge of droplet 1 remains pinned to the surface for a longer duration as the angle subtended by it changes from 100° (advancing contact angle) to 50° (receding contact angle) before it starts moving. This increase in Dt hyst gives sufficient time for the advancing contact edge of droplet 2 to coalesce with the that of droplet 2 leading to the neck formation as shown at T⁄ = 3.5. While the

238

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

(a) T∗ = 1.5

T∗ = 3.5

(c)

(b)

T∗ = 10.0

(d) 1

2.4

2.8

2.2 0.6

Eσ

2.4

ηxy

Xil , Xir

0.8

0.4

2

1.6 0

0.5

1

1.5

2

2.5

0

2

θA = 70◦ θA = 90◦ θA = 120◦

0.2

0

2

4

T∗

6

8

10

1.8

12

θA = 70◦ θA = 90◦ θA = 120◦ 0

2

4

6

T∗

8

10

12

T∗

Fig. 2. (a) Temporal sequence of impacting droplets on a non-ideal surface with a fixed receding contact angle hR = 60° and different advancing contact angles: hA = 70° (top row) and hA = 120°. (b) Evolution of left (X il ) and right (X ir ) intermediate contact positions with respect to time denoted by hollow and solid symbols, respectively. Temporal evolution of (c) non-axial kinetic energy coefficient (gxy ). (d) Surface energy (Er ) for varying hA .

T∗ = 2.25

T∗ = 3.5

T∗ = 6.75 0.7

Z

X

0.6

Y

θR = 50◦ θR = 70◦ θR = 90◦

0.5

A∗

0.4 0.3

X -axis

θR = 50◦ θR = 70◦ θR = 90◦

0.2 0.1

140

130

(b)

120 0

0

2

4

6

8

10

12

T∗

(c) Fig. 3. (a) Temporal sequence of impacting droplets on a non-ideal surface with a fixed advancing contact angle hA = 100° and different receding contact angles: hR = 50° (top row) and hR = 90° (bottom row). (b) Interface profiles of droplet 1 near the contact region for different hR at T⁄ = 2.5. (c) Temporal evolution of the normalized contact area for different hR .

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

increase in hR leads to smaller Dt hyst and the contact edge of droplet 1 retracts earlier. This leads to two different interaction outcomes for the cases considered as observed at T⁄ = 6.75. Along with the decrease in pinning time as hR increases, the surface energy stored in the peripheral rims increases as well. Fig. 3(b) illustrates the interface profiles of droplet 1 near the contact line at T⁄ = 2.5 for three different hR . It is observed that at a given time instant during the receding phase, the interface profile for hR = 90° retracts most followed by hR = 70° and 50°. The magnitude of the stored surface energy inside the peripheral rims is proportional to the curvature subtended by the receding contact lines. Hence, from the interface profiles shown in Fig. 3(b) we can infer that with an increase in hR , surface energy increases. This surface energy is converted into the kinetic energy of the droplet liquid and a part of it is dissipated due to viscous stresses. To further elucidate the influence of this stored surface energy on the kinematics droplet retraction, we illustrate the temporal evolution of the contact area as shown in Fig. 3(c). The contact area is normalized by the initial surface area of a single droplet. Initially, during the inertial spreading phase, A⁄ is independent of hR as the contact line advances with a fixed hA = 100°. However, as time proceeds, A⁄ increases with increase in hR . We observe

239

a sharp difference in change in A⁄ between hR = 70° and 90° at a given instant when compared to that between hR = 50° and 70°. This can be attributed to the fact that for the cases with hR = 50° and 70°, the droplet system undergoes coalescence while there is no interaction between the impinging drops for hR = 90°. As illustrated in this section, we have observed that inhomogeneities in wetting characteristics in the form of contact angle hysteresis will have a significant influence on the impact outcome. Depending on the mode of interaction, this type of surface inhomogeneity could prevent droplet interaction and coalescence, resulting in the formation of nonuniform broken printed lines. Apart from contact angle hysteresis, there exists another form of heterogeneity in the surface wettability which often results in imparting a directional motion to the impinging drops. In the following section, we will investigate the role of this type of wetting inhomogeneity on the interaction dynamics of the impinging drops. 4.2. Role of wettability gradient It has been identified that on an inclined substrate with a wettability gradient, a sessile drop could move uphill against gravity

Fig. 4. Snapshots illustrating the interaction dynamics and internal flow analysis of drops impacting on a non-ideal surface characterised by a gradient in wettability: h2 = 30° (left column) and h2 = 90° (right column).

240

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

[18]. We now focus on understanding the interactions dynamics and elucidate the resulting droplet morphology when two drops impact simultaneously on such a substrate with a constant wettability gradient. Fig. 1(b) outlines the variation of contact angle on the surface along the X-direction. For the present simulations we fix h1 = 120° and vary h2 linearly. This results in a region AB = 2:5Do on the substrate which exhibits a zone of constant wettability gradient. Fig. 4 illustrates the temporal sequence of the impact process on surfaces with wettability gradient corresponding to h2 = 30° and 90° on the left and right columns, respectively. A characteristic facet of drops impacting onto surfaces exhibiting a wettability gradient is the anisotropy in the droplet morphology observed during its recoiling phase [55]. At T  ¼ 3:75, we observe that the drops have undergone coalescence and the droplet system is in its recoiling phase. Both the droplet systems are elongated in the direction of decreasing contact angle i.e., along the direction of the wettability gradient. However, the extent of this elongation and spreading depends on the magnitude of the wettability gradient. Two different interaction mechanisms for neck growth are identified from the cases illustrated in Fig. 4. For the case with h2 ¼ 30 , droplet 2 undergoes rapid asymmetric spreading after it impacts the surface. The fluid particles inside this drop nearly remain at the same position for a considerable time duration, especially those near its right end where the droplet is in contact with the region of higher wettability. This increase in the contact area of droplet 2 facilitates the expansion of the connecting neck. At the same time, impacting on a region with lower wettability, droplet 1 recoils vigorously. Assisted by the interfacial forces, droplet 1 is pulled into droplet 2, resulting in further expansion of the connecting neck as observed at T⁄ = 11.25. As time proceeds, these high

momentum fluid particles inside droplet 1 push the low momentum fluid particles inside droplet 2. To elucidate this flow behavior, Fig. 5(a) illustrates the velocity flow field inside the droplet at T⁄ = 11.5. The vectors are colored by the magnitude of flow velocity. We observe that most of the fluid velocity inside droplet 1 exceeds that inside droplet 2 by an order of magnitude. The drag induced by the motion of the droplet surface on the surrounding fluid sets up a region of high momentum fluid close to this surface. Another noticeable feature which illustrates this behavior is the flow field close to the wall region. While the velocity vectors close to the wall region are nearly zero owing to the no-slip boundary condition, a substantial difference is noticed in the velocity gradients along the wall normal direction inside the droplet bulk. This rise in velocity gradient is attributed to the increased transfer of momentum by the high momentum fluid to the droplet liquid near the wall region. Subsequently, the two drops coalesce to form an asymmetric droplet ensemble which moves along the region of decreasing contact angle (T⁄ = 14.25). In contrast to the case with h2 = 30°, droplet 2 also undergoes vigorous recoiling when the wettability gradient is lowered (h2 = 90°) as shown at T⁄ = 6.0. This rapid recoiling of droplet 2 leads to the contraction in the width of the connecting neck and formation of a slender footprint. The slenderness of the footprint indicates the lower influence of the capillary forces at the contact region. Under the influence of surface tension forces, these two constituent drops mutually pull each other leading to the growth of connecting neck height (T⁄ = 11.25). The fluid particles in the two drops move towards each other, while a slight asymmetry is observed in the shape of droplet 2 due to wettability gradient. The velocity flow field inside the drop shown in Fig. 5(b) illustrates

Z

(a)

Y

X Velocity Magnitude 0.0055 0.0045 0.0035 0.0025 0.0015 0.0005

Z

(b)

Y

X Velocity Magnitude 0.0045 0.0035 0.0025 0.0015 0.0005

Fig. 5. Velocity fields inside the droplet system at T⁄ = 11.5 for (a) h2 = 30° and (b) h2 = 90°. The vectors are colored according the velocity magnitude. The flow field is extracted along the mid y cross-sectional plane passing through the centre of drops and the simulation domain. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

241

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

this mode of neck growth. We notice that velocity vectors inside the two drops move towards each other, with the fluid momentum being higher inside droplet 1. A clear distinction in the interaction mode can be seen from the vector field shown in Fig. 5(a) where the direction of liquid flow in both constituting drops is along the same direction. These two countercurrents inside the droplet system interact and rise upwards leading to the growth in droplet height. Consequently, the drag induced by the droplet surface on the surrounding fluid leads to the formation of a vortex pattern near the top region of drop 2. This results due to the displacement of the low momentum ambient fluid by the high momentum ambient fluid near the top surfaces of droplet 2 and droplet 1, respectively. Mechanisms similar to the push and pull mechanisms have also been observed by for a droplet system interacting with a surface characterized by wettability have also been identified by Randive et al. [65]. The authors investigated the behavior of a droplet on an inclined surface with a sharp change in the surface wettability. They suggested a ‘‘push” and ‘‘pull” type of mechanisms when the droplet crosses the junction from the hydrophobic to the hydrophilic region. To further elucidate the impact kinetics of droplet interactions on surfaces exhibiting wettability gradient, we monitor the displacement of the right and the left contact edges of the droplet system in Fig. 6(a). For a given fluid-fluid system, a lower contact angle implies that the solid substrate has high surface energy [66]. This leads to greater interfacial spreading on the surface. As such, we observe that the temporal evolution of X r increases with

a decrease in h2 . The right contact edge continues to propagate downstream as time proceeds. This illustrates the influence of the high surface energy on the contact line, which in turn facilitates its rapid motion on the substrate. However, this increasing trend in the evolution of X r decreases as h2 increases, with retraction of the contact edge observed for h2 = 90°. The temporal evolution of the left contact edge remains nearly independent of h2 as it resides primarily in the region of constant wettability for the most time duration of the impact and interaction process. Therefore, the elongation of the droplet system increases as the wettability gradient increases. This observation complements with the morphology of the droplets shown in Fig. 4. To further quantify this stretching and elongation of the droplet ensemble, Fig. 6(b) shows the temporal evolution of the surface energy (Er ) as h2 varies. We notice a sharp rise in the initial values of Er for all the considered cases with the maximum values noted for lower values of h2 . This is attributed to the inertial spreading and subsequent coalescence of the two drops leading to rapid increase in Er . As mentioned previously, increase in wettability gradient leads to greater expansion of the connecting neck width and the contact area of droplet 2. This results in high maximum values of Er for h2 = 30°. As droplet 1 pulls into droplet 2, the droplet system continues to spread with the curvature of the contact line adjusting to local contact angles for the case with h2 = 30°. This pulling-in of droplet 1 facilitates the downward motion of the droplet system as observed from the temporal evolution of X cm center of mass of the system shown in Fig. 6(c). However, we observe a decreasing trend in the temporal evolution

2.6

3.6

2.4

Xr

Eσ

4

3.2

2.2

1.4

Xl

1.2 1

θ2 = 30◦ θ2 = 60◦ θ2 = 90◦

0.8

2.8

0.6 0.4 0

0

2

4

2

2.6

T∗

6

θ2 = 30◦ θ2 = 60◦ θ2 = 90◦

2

8

4

6

1.8

8

0

2

4

6

T∗

T∗

(a)

(b)

0.8

θ2 = 30◦ θ2 = 60◦ θ2 = 90◦

0.6

Xcm

η grad

2.4

2.2

0.4

θ2 = 30◦ θ2 = 60◦ θ2 = 90◦

0.2 2 0

2

8

4

6

8

0

2

4

6

T∗

T∗

(c)

(d)

8

Fig. 6. Time evolution of (a) displacement of right (X r ) and left (X l , shown inset) contact edges. (b) Surface energy (Er ), (c) location of X-center of mass (X cm ) of the droplet system, and (d) kinetic energy coefficient (gx ) along the wettability gradient direction for different h2 .

242

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

of X cm as h2 is increased. This behavior can be explained as follows: Upon impacting, when the droplets spread and coalesce on a surface exhibiting wettability gradient, the initial kinetic energy of the drops in converted into surface energy. As the drops interact and coalesce, the surface energy is converted back into kinetic energy which is primarily transferred along the axial direction, with a part of it being dissipated by viscous stresses. However, due to a surface wettability gradient, a part of this axial kinetic energy is transferred along the direction of increasing surface energy of the substrate. To quantify this transferred kinetic energy along the X-direction (direction of wettability gradient), we define the kinetic energy coefficient along the gradient direction (ggrad ):

ggrad ¼ R ½

½ Vdrop

R Vdrop

ð1=2Þqðu2x ÞdV

ð24Þ

ð1=2Þqðu2x þ u2y þ u2z ÞdV

Here we use ux as the reference velocity in the numerator as it indicates the direction of the wettability gradient. The evolution of ggrad illustrated in Fig. 6(d) clearly indicates the increase in the transfer of kinetic energy along the direction of wettability gradient as h2 decreases. As such, with the increase in the transferred momentum along the X direction, the droplet system moves more readily with the increase in wettability gradient. This lead to directional deposition of the droplet system as illustrated in Fig. 6(c). Conversely, as the magnitude of the wettability gradient is reduced, the resistance offered to this directional momentum transfer increases, resulting in vigorous recoiling of droplet 2. Therefore, variations in Er are observed as time proceeds for the cases with a lower wettability gradient. Apart from the inertia of the impinging drops, the mobility and extent of deformation of the droplet system are also governed by the viscous resistance and interfacial tension. To investigate the influence of these physical factors on the mobility of the droplet system, Fig. 7 illustrates the temporal evolution of X cm for different (a) Oh and (b) We numbers. Increase in the Oh number leads to increase in the viscous friction offered to the deforming fluid element. This, in turn, inhibits inertial spreading of the two droplets. As the droplet system moves downstream towards more hydrophilic regions, lengthening of the droplet footprint occurs. In addition, larger internal shear rates are induced inside the droplet system as the height of droplet 2 decreases when it moves downstream and droplet 1 increases due to recoiling. Therefore as Oh number increases, the internal viscous drag increases, which leads to a decrease in the temporal evolution of X cm . Fig. 7(b) shows that as

We increases, the temporal evolution of X cm decreases. This is attributed to the fact that with an increase in surface tension, the momentum transferred along the gradient direction when the recoiling droplet 1 is pushed into droplet 2 is higher. This results in greater displacement of the droplet system in the downstream direction as surface tension increases. From the preceding section, we have observed that a wetting inhomogeneity in the form of a wettability gradient leads to different mechanisms in which the droplets interact with each other. Substantial droplet elongation resulting due to the directional bias in the surface wettability would lead to the formation of uneven printed lines. For a fixed set of impact velocities and the chemical properties of the droplet fluid, the interaction dynamics would not only depend on the wetting characteristics of the target surface but also on the shape of the impinging drops. Previous studies [62,32,31] investigating the interactions between drops impinging onto dry surfaces have assumed a spherical droplet shape at the point of impact. However, in real applications, the drops undergo mid-air oscillations before impacting the target surface. This leads to deviations in the impact shape of the droplets from the assumed spherical geometry to prolate and oblate shapes. The next section addresses the role of the impact shape on the interaction dynamics of simultaneously impinging drops on dry surfaces. 4.3. Role of impact shape As the shape of the drops change from oblate to prolate while they travel in air, we consider an ellipsoidal geometry to model droplet shapes in the current study. We initialize the shape of the ellipsoidal droplets in the computational domain by the following equation:

x y z þ þ ¼1 ðao =2Þ ðbo =2Þ ðco =2Þ

where ao and co are the major and minor axes of the cross-section for the prolate ellipsoid shapes, respectively and vice versa for the oblate ellipsoid as shown in Fig. 1 (c). The aspect ratio (AR) is defined as the ratio between the major and minor axes and is varied between 1.0 and 1.8 in the current study. The principle semi axis bo is set to be equal to the minor axis and the static contact angle is fixed at 140 . Fig. 8(a) illustrates the temporal behavior of the interaction dynamics between drops with for spherical
spherical (SS, left column), oblate
oblate (OO, middle column) and prolate
prolate (PP, right column) droplet combinations. The

(a) Oh = 0.100 Oh = 0.050 Oh = 0.016

W e = 12.6 W e = 25.0 W e = 51.2

2.6

Xcm

X cm

(b)

2.8

2.4

2.2

ð25Þ

2.4

2.2 2 2 0

2

4

6

T∗

8

0

2

4

6

8

T∗

Fig. 7. Time evolution of X-center of mass (X cm ) location of the droplet system along the wettability gradient direction for different (a) Ohnesorge number (Oh) and (b) Weber number (We).

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

243

Fig. 8. Temporal sequence of droplets with an aspect ratio (AR) of 1.8 on a surface with h = 140° for Spherical-Spherical (left column), Proate-Prolate (middle column) and Oblate-Oblate (right column) impact shape combinations.

aspect ratio fort he OO and PP droplet combinations is fixed at AR = 1.8. At T⁄ = 0.75 we notice that the interaction time between the colliding rims of the two droplets is dependent on the impact shape. For the SS combination, we observe that the formation of connecting neck has been initiated. As the rims touch each other in contrast to the PP combination where the collision between rims is yet to occur. When compared to the SS combination, a similar but faster interaction process between the impinging drops is noticed for the OO combination. A distinct shape dependent feature of the interaction process is the formation of a connecting neck for the SS and OO combinations whereas a connecting ridge is formed when the impinging drops are prolate in shape as observed at T⁄ = 1.5. The thickness of the peripheral rims of the droplet system is greater for the OO combination followed by SS and PP droplet combinations. As time proceeds and the droplet system continue to recoil, we observe the characteristic concave curvature between the interacting drops for the SS and OO cases. This results due to the interplay between the uprising high momentum fluid flow inside the bulk region of the constituent drops and the low momentum flow around the connecting neck (T⁄ = 2.5–4.0). As the two lobes, driven

by the capillary pressure arising due to the negative neck curvature are drawn closer, the height of the connecting neck increases and the constituent drops coalesce (T⁄ = 6.0–9.0). However, for the PP combination we observe that the uprising motion of the high momentum fluid in the bulk region of the constituent drops is inhibited, leading to the formation of a connecting plateau as shown at T⁄ = 2.5. The plateau gradually increases its height as the droplet system undergoes coalescence. As the elongation of the constituent drops along the axial direction is restricted in this combination, the droplet system occupies a larger footprint when compared to the other two cases considered as shown at T⁄ = 9.0. For a PP droplet combination as shown in Fig. 1(c)–(i), the asymmetry in the shape of the impinging drops is such that the plane enclosing the maximum area of the cross-section of the drop is perpendicular to the target surface. Therefore, the surface energy in the droplet system at the time of impact (T⁄ = 0) is higher when compared to its spherical counterpart. Hence, to minimize its energy the droplet system would tend to squeeze-in centrally along the x-plane to a attain a spherical equilibrium configuration. When the droplets impinge onto the target surface and undergo

244

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

inertial spreading, this capillary driven squeezing fluid flow along the x-plane facilitates the motion of the droplet liquid along the x-axis. This leads to an increase in the momentum of the droplet liquid inside the interacting rims and results in the formation of a connecting ridge. Fig. 9(a) illustrates the velocity field inside the connecting ridge at T⁄ = 5.0 along the mid y plane. From the velocity vectors shown in the figure, it is evident that driven by the high momentum droplet fluid, the two rims colloid and propagates upwards leading to the formation of a connecting ridge. In contrast, the velocity field inside the connecting region for the OO drop combination shown in Fig. 9(b) is lower in magnitude indicating the formation of a connecting neck. In this case, the plane enclosing the maximum area of the cross-section of the droplets is parallel to the target surface, the capillary driven squeezing flow is centrally along the z-plane. This restricts droplet spreading and leads to greater accumulation of the droplet liquid inside the peripheral rims. This resistance to droplet spreading is observed from the velocity vectors shown inside the peripheral rims shown in Fig. 9 (d). However facilitated by capillary squeezing, enhancement of droplet spreading for PP droplet combination as observed in Fig. 9(c). The combined effect of droplet compression along the y-axis and enhancement in spreading along the x-axis results in the formation of the connecting plateau as show at T⁄ = 2.5 in Fig. 8. We next investigate the effect of aspect ratio on the interaction dynamics for the configurations considered in this work. Fig. 10(a) shows the temporal evolution of the spread factor (D ) for the OO combination with varying aspect ratios. We notice that variation in the evolution of D is nearly independent with AR. The maxi-

(a)

Z

Y

(c)

(b) X

Z

Y

mum spread factor for all the considered cases is nearly equal with a difference in the maximum spread factors of 0.32%. However, we observe a slight shift in the temporal evolution of D towards its left-hand side with an increase in aspect ratios. This leads to an important observation that when ellipsoidal drops impinge simultaneously, such that the plane enclosing the cross section with maximum area is parallel to the target surface, the interaction dynamics follows a similar trend as that of its spherical counterpart, albeit at a faster rate. The difference in this rate increases with increase in aspect ratio. To further elucidate this behavior, we monitor the temporal evolution of the surface energy for different aspect ratios as shown in Fig. 10(b). We observe that the temporal variation in Er decreases with increase in AR. This illustrates the resisting effect of the oblate geometry on the interaction dynamics. With the increase in AR, the constituent drops are more stretched and have higher surface energy at the time of impact. This initial Er leads to the inward squeezing flow discussed previously and the droplet surface is compressed more readily leading to lower Er when compared to the base case with AR = 1.0. It should be noted that even though the temporal evolution of D is nearly same, there is a substantial difference in the corresponding Er . This is attributed to the fact that the evolution of D depends on the contact footprint which is dependent on the contact angle after inertial spreading. The restoring effects of the initial increase in Er primarily influence the bulk droplet fluid leading to thicker peripheral rims and faster recoiling of the droplet system. A complete reversal in this trend is observed for the PP droplet combination. We notice an increase in the temporal evolution of both D and Er with an

Y

(d) X

Z

X

Z

Y

X

Fig. 9. Snapshots illustrating the velocity vectors along the mid Y plane at T  = 1.25 for the Prolate-Prolate (a and c) and Oblate-Oblate (b and d) droplet combinations. The aspect ratios of the impacting drops are set to be 1.8. Plots (a)–(b) and (c)–(d) show the velocity field inside the connecting region and the peripheral rims, respectively.

245

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

3.5

3.5

(a)

3

3

AR = 1.0 AR = 1.4 AR = 1.8

2.5

2.5

D∗

D∗

2

2

1.5

1.5

1

1

0.5

0.5

0 0

2

4

6

8

10

AR = 1.0 AR = 1.4 AR = 1.8

0

12

0

2

4

T∗

8

10

12

10

12

2.8

(b)

2.6

6

T∗

2.6 2.4

Eσ

Eσ

2.4 2.2

2.2

AR = 1.0 AR = 1.4 AR = 1.8

2

AR = 1.0 AR = 1.4 AR = 1.8

2

1.8 0

2

4

6

8

10

12

T∗

1.8

0

2

4

6

8

T∗

Fig. 10. Temporal evolution of (a) Spread factor (D ) (b) Surface energy (Er ) for different aspect ratios (AR) of the droplet. The left and right columns corresponds to the OblateOblate and Prolate-Prolate combinations, respectively.

increase in aspect ratio. In this case, the initial increase in Er due to the impact shape not only influences the bulk flow inside the drop but also results in compression and elongation of the contact line along the y-axis and x-axis, respectively. This difference in the interaction dynamics for the OO and PP droplet combinations is primarily attributed to the orientation of the target surface with the plane enclosing the droplet cross-section with maximum area. Moreover, the discussion on the influence of different combinations of the impact shape on the interaction dynamics has been provided in the supporting material. The results presented above elucidate the importance inhomogeneities in the wetting characteristics of the surface and the role of impact shape on the interaction dynamics of simultaneously wall impinging drops. Considerable deviations in the interaction mechanisms and impact outcome from those modeled under constant contact angle approximation are observed when nonuniform surface wettability is taken into account. Similarly, departure from the spherical droplet geometry to asymmetric impact shapes at the time of impact have shown to yield different spreading and recoiling dynamics of the droplet system. As the nature of all target surfaces used in engineering applications exhibit such nonideal surface characteristics and the impinging drops oscillate mid-air before impact, these numerical observations and discussions provide a more realistic insight on the droplet interaction dynamics for situations observed in real-time applications. 5. Conclusions In this paper, we have numerically investigated the interaction dynamics between two simultaneously wall impinging drops,

including the effects of wetting inhomogeneity and impact shape. Various dynamics and morphological characteristics were observed through numerics: interaction outcome is determined by the contact angle hysteresis, a gradient in surface wettability would result in the formation of a directional topology of the droplet system or a change in droplet shape would lead varying connecting morphological structures. A high-density ratio based three-dimensional lattice Boltzmann solver is employed for the present two-phase computations and the dynamics of the moving contact line is captured through a geometry based contact angle model. For a fixed hR , the connecting time of the interacting rims increases with increase in hA for out-ofphase collision mode. This delay in connecting time would induce sufficient irregularities in the shape of the final contact footprint which in turn would impact the quality of the printed lines. Similarly, as hR increases for the in-phase mode, there is no interaction between the impinging drops which would lead to broken printed lines. Compared to the previous 2D [32] and 3D [67,33] works employing a static contact angle model to investigate the interactions on simultaneous, the current work establishes the importance to consider the effects contact angle hysteresis for more realistic modeling. We showed that the previously identified interaction modes for simultaneously impinging drops [33] on a surface with constant contact angle are sensitive to the contact angle hysteresis. A gradient in surface wettability induces a directional asymmetry in the morphology of the droplet system. We have shown two different interaction mechanisms depending on the magnitude of wettability gradient. When the magnitude of the wettability gradient is greater, the growth in the connecting region at later stages is primarily facilitated by the motion of droplet 1.

246

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247

However, both the drops contribute to the evolution of the connecting region for surfaces characterized by lower the wettability gradients. Hence, it is recommended that for uniformity of line patterns for printing or coating applications, wettability gradient induced surface inhomogeneity should be minimized. We do not find mixing enhancement when drops impinge simultaneously onto a surface with wettability gradient in contrast to the situations of a sessile droplet impact [25]. While it is important to have a good droplet mixing in additive manufacturing applications, such resistance for mixing is favorable in graphical application to reduce pixelation and improve printing resolution. Finally, we have shown that the geometry of the impact shape has a significant effect on the interaction dynamics of simultaneously impinging drops. Depending on the orientation of the plane enclosing the maximum cross-sectional area with the target surface, different connecting structures and spreading behaviors are observed. The interaction dynamics of Oblate
Oblate droplet combination is similar to its Spherical
Spherical counterpart, albeit at a faster rate. The Prolate
Prolate combination showed an increase in the spreading and recoiling rates when compared to the base case. Therefore, for a given set of impact conditions, when compared to spherical drops interaction between prolate droplets would result in higher maximum spread factors. While the present work accounts surface homogeneity for smooth surfaces, in practice such non-uniformity is also accounted due to surface roughness. This leaves scope for future work to understand the interaction dynamics of multiple impinging drops. Owing to the high computational cost, we have restricted our study to two droplet interactions. Even though the basic interaction behaviors are captured with two drop interactions, there is a need to scale up for a large number of drops for more realistic understanding. Moreover, employing a dynamic contact angle model is essential to yield an improved quantitative agreement with experimental findings. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.jcis.2018.01.063. References [1] I.V. Roisman, B. Prunet-Foch, C. Tropea, M. Vignes-Adler, Multiple drop impact onto a dry solid substrate, J. Coll. Interf. Sci. 256 (2002) 396–410. [2] K.A. Raman, R.K. Jaiman, T.S. Lee, H.T. Low, Lattice boltzmann study on the dynamics of successive droplets impact on a solid surface, Chem. Eng. Sci. 145 (2016) 181–195. [3] D. Soltman, V. Subramanian, Inkjet-printed line morphologies and temperature control of the coffee ring effect, Langmuir 24 (2008) 2224– 2231. [4] R. Rioboo, M. Marengo, C. Tropea, Time evolution of liquid drop impact onto solid, dry surfaces, Exp. Fluids 32 (2002) 112–124. [5] R. Rioboo, M. Marengo, C. Tropea, Outcomes from a drop impact on solid surfaces, Atom. Spray 11 (2002) 155–165. [6] D.A. Weiss, A.L. Yarin, Single drop impact onto liquid films: neck distortion, jetting, tiny bubble entrainment, and crown formation, J. Fluid Mech. 385 (1999) 229–254. [7] M. Rein, Phenomena of liquid drop impact on solid and liquid surfaces, Fluid Dynam. Res. 12 (1993) 61–93. [8] J. Eggers, J.R. Lister, H.A. Stone, Coalescence of liquid drops, J. Fluid Mech. 401 (1999) 293–310. [9] L. Duchemin, J. Eggers, C. Josserand, Coalescence of liquid drops, J. Fluid Mech. 487 (2003) 167–178. [10] D.G.A.L. Aarts, H.N.W. Lekkerkerker, H. Guo, G.H. Wegdam, D. Bonn, Hydrodynamics of droplet coalescence, Phys. Rev. Lett. 95 (164503) (2005). [11] S. Karptischka, H. Reigler, Sharp transition between coalescence and noncoalescence of sessile drops, J. Fluid Mech. Rapids 743 (R1) (2014). [12] S. Karptischka, H. Reigler, Quantitative experimental study on the transition between fast and delayed coalescence of sessile droplets with different but completely miscible liquids, Langmuir 26 (2010) 11823–11829. [13] A.Y. Tong, S. Kasliwal, H. Fujimoto, On the successive impingement of droplets onto a substrate, Numer. Heat Transf., Part A 52 (2007) 531–548.

[14] L. Gao, T.J. McCarthy, Teflon is hydrophilic. comments on definitions of hydrophobic, shear versus tensile hydrophobicity and wettability characterization, Langmuir 24 (2008) 9183–9188. [15] M.A. Nilsson, J.P. Rothstein, The effect of contact angle hysteresis on droplet coalescence and mixing, J. Coll. Interf. Sci. 363 (2011) 646–654. [16] Worim Lee, Gihun Son, Numerical study of droplet impact and coalescence in a microline patterning process, Comp. Fluids 42 (2011) 26–36. [17] L. Zhang, Y. Zhu, X. Cheng, Numerical investigation of multi-droplets deposited lines morphology with a multiple-relaxation-time lattice boltzmann model, Chem. Eng. Sci. 171 (2017) 534–544. [18] M.K. Chowdhury, G.M. Whitesides, How to make water run uphill, Science 256 (6) (1992) 1539–1541. [19] J.E. Longley, E. Dooley, D.M. Givler, W.J. Napier III, M.K. Chaudhury, S. Danie, Drop motion induced by repeated stretching and relaxation on a gradient surface with hysteresis, Langmuir 28 (2012) 13912–13918. [20] A.K. Das, P.K. Das, Simulation of drop movement over an inclined surface using smoothed particle hydrodynamics, Langmuir 25 (2009) 11459–11466. [21] Z. Li, Guo-Hui Hu, Zhi-Liang Wang, Yan-Bao Ma, Zhe-Wei Zhou, Three dimensional flow structures in a moving droplet on substrate: a dissipative particle dynamics study, Phys. Fluids 25 (072103) (2013). [22] F. Varnik, P. Truman, B. Wu, P. Uhlmann, D. Raabe, M. Stamm, Wetting gradient induced separation of emulsions: a combined experimental and lattice Boltzmann computer simulation study, Phys. Fluids 20 (072104) (2008). [23] S. Gong, P. Cheng, Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows, Comp. Fluids 53 (2012) 93–104. [24] M. Ahmadlouydarab, J.J. Feng, Motion and coalescence of sessile drops driven by substrate wetting gradient and external flow, J. Fluid Mech. 746 (2014) 214–235. [25] J.R. Castrejon-Pita, K.J. Kubiak, A.A. Castrejon-Pita, M.C.T. Wilson, I.M. Hutchings, Mixing and internal dynamics of droplets impacting and coalescing on a solid surface, Phys. Rev. E 88 (023023) (2013). [26] H. Fujimoto, S. Ito, I. Takezaki, Experimental study of successive collision of two water droplets with a solid, Exp. Fluids 33 (2002) 500–502. [27] A. Dalili, S. Chandra, J. Mostaghimi, H.T.C. Fan, J.C. Simmer, Formation of liquid sheets by deposition of droplets on a surface, J. Coll. Interf. Sci. 418 (2014) 292–299. [28] T. Zhang, H.M. Tsai, J.L. Alvarado, Effects of single and double streams of droplet impingements on surface cooling, Atom. Sprays 24 (2014) 875– 893. [29] G.E. Soriano, T. Zhang, J.L. Alvarado, Study of the effects of single and multiple periodic droplet impingements on liquid film heat transfer, Int. J. Heat Mass Transf. 77 (2014) 449–463. [30] T. Zhang, J.L. Alvarado, J.P. Muthusamy, A. Kanjirakat, R. Sadr, Heat transfer characteristics of double, triple and hexagonally-arranged droplet train impingement arrays, Int. J. Heat Mass Transf. 110 (2017) 562–575. [31] S. Batzdorf, J. Breitenbach, C. Schlawitschek, I.V. Roisman, C. Tropea, P. Stephan, T. Gambaryan-Roisman, Heat transfer during simultaneous impact of two drops onto a hot solid substrate, Int. J. Heat Mass Transf. 113 (2017) 898– 907. [32] J. Wu, J.J. Huang, W.W. Yan, Lattice boltzmann investigation of droplets impact behaviours onto a solid surface, Coll. Surf. A: Physicochem. Eng. Asp. 484 (2015) 318–328. [33] K.A. Raman, R.K. Jaiman, T.S. Lee, H.T. Low, Dynamics of simultaneously impinging drops on a dry surface: Role of impact velocity and air inertia, J. Coll. Interf. Sci. 486 (2017) 265–276. [34] K.A. Raman, R.K. Jaiman, T.S. Lee, H.T. Low, On the dynamics of crown structure in simultaneous two droplets impact onto stationary and moving liquid film, Comp. Fluids 107 (2015) 285–300. [35] H. Shetabivash, F. Ommi, G. Heidarinejad, Numerical analysis of droplet impact onto liquid film, Phys. Fluids 26 (2014). [36] S. Yun, G. Lim, Ellipsoidal drop impact on a solid surface for rebound suppression, J. Fluid Mech. 752 (2014) 266–281. [37] S. Yun, J. Hong, K.H. Kang, Suppressing drop rebound by electrically driven shape distortion, Phys. Rev. E 87 (033010) (2013). [38] L.S. Luo, Theory of the lattice Boltzmann method: lattice boltzmann models for nonideal gases, Phys. Rev. E 63 (2000) 4982–4996. [39] S. Chen, G. Doolen, Lattice boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [40] Shufang Zhao, Antoine Riaud, Guangsheng Luo, Yong Jin, Yi Cheng, Simulation of liquid mixing inside micro-droplets by a lattice boltzmann method, Chem. Eng. Sci. 28 (2015) 118–128. [41] Luz Amaya-Bower, Taehun Lee, Numerical simulation of single bubble rising in vertical and inclined square channel using lattice boltzmann method, Chem. Eng. Sci. 66 (2011) 935–952. [42] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991). 4320-1127. [43] D.H. Rothman, J.M. Keller, Immiscible cellular-automaton fluids, J. Statist. Phys. 52 (1988) 1119–1127. [44] X. Shan, H. Chen, Lattice boltzmann for simulating flows with multiple phases and components, Phys. Rev. E 47 (3) (1993) 1815–1819. [45] M.R. Swift, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann for simulations for liquid-gas and binary fluid systems, Phys. Rev. E 54 (1996) 5041–5052. [46] X. He, X. Shan, R. Zhang, A lattice boltzmann scheme for incompressible multiphase flow and its application in simulation of rayleigh Taylor instability, J. Comput. Phys. 152 (1999) 642–663.

K. Ashoke Raman / Journal of Colloid and Interface Science 516 (2018) 232–247 [47] H. Huang, M.C. Sukop, X.Y. Lu, Multiphase Lattice Boltzmann Methods: Theory and Applications, Wiley-Blackwell, 2015. [48] T. Lee, C.L. Lin, A stable discretization of the lattice boltzmann equation for simulation of incompressible two-phase flows at high density ratio, J. Comput. Phys. 206 (2005) 16–47. [49] T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible two-phase flows with large density differences, J. Comput. Phys. 198 (2004) 628–644. [50] H.W. Zheng, C. Shu, Y.T. Chew, A lattice Boltzmann model for multiphase flows with large density ratio, J. Comput. Phys. 218 (2006) 353–371. [51] D. Jamet, O. Lebaigue, N. Coutris, J.M. Delhaye, The second gradient method for the direct numerical simulation of liquid vapor flows with phase change, J. Comput. Phys. 169 (2001) 624–651. [52] X. He, X. Shan, G.D. Doolen, Discrete Boltzmann equation model for nonideal gases, Phys. Rev. E 57 (1) (1998). [53] T. Lee, C.L. Lin, Pressure evolution lattice Boltzmann-equation method for two phase flow with phase change, Phys. Rev. E 67 (056703) (2003). [54] H. Ding, P.D.M. Spelt, Wetting condition in diffuse interface simulations of contact line motion, Phys. Rev. E 75 (046708) (2007). [55] K.A. Raman, R.K. Jaiman, T.S. Lee, H.T. Low, Lattice Boltzmann simulations of droplet impact onto surfaces with varying wettabilities, Int. J. Heat Mass Transf. 95 (2016) 336–354. [56] T. Young, An essay on the cohesion of fluids, Philos. Trans. R. Soc. Lond. 95 (1805) 65–87.

247

[57] E.B. Dussan V, On the spreading of liquids on solid surfaces: static and dynamic contact lines, Annu. Rev. Fluid Mech. 11 (1979) 371–400. [58] R.G. Cox, The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow, J. Fluid Mech. 168 (1986) 169–194. [59] P. Sheng, M. Zhao, Immiscible-fluid displacement: contact-line dynamics and the velocity-dependent capillary pressure, Phys. Rev. A 45 (5694) (1992). [60] Y. Sui, P.D.M. Spelt, An efficient computational model for macroscale simulations of moving contact lines, J. Comput. Phys. 242 (2013) 37–52. [61] Y. Sui, H. Ding, P.D.M. Spelt, Numerical simulations of flows with moving contact lines, Annu. Rev. Fluid Mech. 46 (2014) 97–119. [62] P.J. Graham, M.M. Farhangi, Ali Dolatabadi, Dynamics of droplet coalescence in response to increasing hydrophobicity, Phys. Fluids 24 (112105) (2012). [63] X. Yang, V.H. Chhasatia, Y. Sun, Oscillation and recoil of single and consecutively printed droplets, Langmuir 29 (2013) 2185–2192. [64] K.A. Raman, R.K. Jaiman, Y. Sui, T.S. Lee, H.T. Low, Rebound suppression of a droplet impacting on an oscillating horizontal surface, Phys. Rev. E 94 (023108) (2016). [65] P. Randive, A. Dalal, K.C. Sahu, G. Biswas, P.P. Mukherjee, Wettability effects on contact line dynamics of droplet motion in an inclined channel, Phys. Rev. E 91 (053006) (2015). [66] E.Y. Bormashenko, Wetting of Real Surfaces, De Gruyter, 2013. [67] W. Zhao, D. Loney, A.G. Fedrov, F.L. Degertekin, D.W. Rosen, Lattice boltzmann simulations of multiple-droplet interaction dynamics, Phys. Rev. E 89 (033311) (2013).