A numerical analysis of air entrapment during droplet impact on an immiscible liquid film

A Numerical Analysis of Air Entrapment During Droplet Impact on an Immiscible Liquid Film

Journal Pre-proof

A Numerical Analysis of Air Entrapment During Droplet Impact on an Immiscible Liquid Film Firoozeh Yeganehdoust , Reza Attarzadeh , Ida Karimfazli , Ali Dolatabadi PII: DOI: Reference:

S0301-9322(19)30536-1 https://doi.org/10.1016/j.ijmultiphaseflow.2019.103175 IJMF 103175

To appear in:

International Journal of Multiphase Flow

Received date: Revised date: Accepted date:

25 July 2019 14 November 2019 20 November 2019

Please cite this article as: Firoozeh Yeganehdoust , Reza Attarzadeh , Ida Karimfazli , Ali Dolatabadi, A Numerical Analysis of Air Entrapment During Droplet Impact on an Immiscible Liquid Film, International Journal of Multiphase Flow (2019), doi: https://doi.org/10.1016/j.ijmultiphaseflow.2019.103175

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

1

2

3

4

5

6

7

Highlights • Comprehensive numerical study of droplet impact and spontaneous bouncing from an immiscible liquid film surface using VOF method • Characterizing the liquid film and the air cushion hydrodynamics for various Weber numbers and film thicknesses • Increasing the probability of droplet detachment due to the combined effect of air entrapment and the liquid film interface deformation

1

9

A Numerical Analysis of Air Entrapment During Droplet Impact on an Immiscible Liquid Film

10

Firoozeh Yeganehdoust, Reza Attarzadeh, Ida Karimfazli, Ali Dolatabadi∗

8

11 12

13

Department of Mechanical, Industrial, and Aerospace Engineering, Concordia University, Montreal, Quebec, Canada

Abstract The air entrapment during droplet impingement is responsible for spontaneous droplet bouncing on an arbitrary solid surface at low Weber numbers. However, for the impact on liquid film surfaces, the outcome would significantly change, making it more favorable for the fabrication of non-wetting lubricant impregnated surfaces (LIS/SLIPS). In this paper, we describe a problem associated with the impact on a liquid surface using a three-phase flow model that captures the details of the gas layer thickness and dynamics of fluid motions. The numerical model was based on the finite volume solution coupled with the volume of fluid method to track the phases. The model was validated with the analytical solution. Consequently, the numerical tool was utilized to investigate the thickness of the entrapped air during the impact process while the behavior of droplet and the immiscible liquid film was quantitatively measured. The morphology of the interfacial gas layer was analyzed for key parameters including impact velocity and film thickness. It was observed that the presence of liquid film can reduce the probability of rupturing the gas layer. The results for the profile of liquid film during droplet impact illustrated that the effect of film thickness can considerably influence the bouncing behavior.

14

Keywords: Droplet impact, Liquid film, LIS, Spontanous bouncing

∗ Corresponding

Author: [email protected]

Preprint submitted to Journal of LATEX Templates

November 23, 2019

15

1. Introduction

16

The phenomenon of droplets impacting surfaces has been studied due to

17

various interesting physical phenomena and its importance in applications such

18

as thermal spraying, coatings and inkjet printing [1, 2, 3, 4]. Depending on the

19

application, the impact of droplet can be classified according to the impacting

20

surface which may be a dry, thin liquid film or a deep liquid pool surface. In all

21

cases, droplet impact consists of several discernible regimes including splashing,

22

bouncing, and coalescing and understanding the physics behind the phenomena

23

are of critical importance [5, 6, 7].

24

The outcome of droplet impact dynamics is mostly the inteplaying between

25

surface tension, inertia and viscous force. Accordingly, these regimes have been

26

characterized with the definition of dimensionless numbers, such as the Reynolds

27

number (Re), the Weber number (W e), and the Ohsenorge number (Oh). In

28

the case of a wet surface, there are several key parameters such as liquid film

29

thickness, liquid properties, and impact velocity that modulate the droplet dy-

30

namics behavior. The interacting dynamics of both liquids (droplet and liquid

31

film) increases the complexity of the droplet impact process. Numerous exper-

32

imental, theoretical, and numerical studies [8, 9, 10] have been conducted on

33

droplet impacts on dry and wet surfaces (thin, thick and deep pool liquid film)

34

to further define these different phenomena under various conditions..

35

An important associated phenomenon that can occur during the bouncing of

36

a droplet from a solid surface is the formation of a thin gas layer underneath the

37

center of the droplet upon the impact. This thin gas layer forms at low values

38

of the Weber number and there have been few experimental observation using

39

optical techniques [11, 12, 13]. De Ruiter et al. [12] reported that a droplet

40

can bounce off a hydrophilic surface due to the creation of the thin air cushion,

41

which they observed using reflection interference contrast microscopy (RICM).

42

Some studies have illustrated that a high-pressure air pocket forms between the

43

droplet and the substrate. This air pocket deforms the bottom surface of the

44

droplet into a dimple-like shape and may endure throughout the impact until

3

45

the droplet bounces. Further experimental studies also measured the thickness

46

of the gas layer for a droplet impacting on a solid surface [11, 12, 14, 15].

47

Numerical simulation have also confirmed the presence of the air layer and have

48

provided further details on its dynamics [16, 17, 18].

49

Although most studies have been dedicated to analyzing the air cushion

50

during the impact process on a solid surface, the mechanism for droplet impact

51

on a liquid surface is still far from being fully understood.

52

Droplet impact on liquid surfaces is important in several industrial applica-

53

tions and recently in lubricant impregnated surfaces (LIS) due to their promising

54

applications in omniphobic surfaces [19, 20, 21]. Complicated phenomena may

55

occur during impact on a liquid surface such as crown formation and splashing

56

that are influenced including different parameters of Weber and Reynolds num-

57

bers, dimensionless film thickness and viscosity ratios. Regarding the splashing

58

phenomena, several studies introduced threshold splashing correlations for dif-

59

ferent conditions [22, 23, 24]. Chen et al. [25] performed experimental study of

60

impact onto thin liquid film of miscible and immiscible liquids. They observed

61

that the miscibility can affect the impact dynamics such as the crown formation

62

and splashing, while the effect of film thickness is more complex. Amongst com-

63

plicated phenomena during impact on liquid film, the air entrapment has been

64

received less attention. For instance, characterizing the behavior of liquid film

65

in LIS is vital for designing the slippery surface. The impact on finite liquid

66

film thickness results in merging or bouncing, Tang et al. [26] demonstrated

67

that for impact inertia less than the critical value, increasing the film thickness

68

leads to a non-monotonic transition from merging to bouncing to merging and

69

finally to bouncing again. This illustrated the role of the film thickness in the

70

impact dynamics. Tran et al. [27] observed that for drops impacting a liquid

71

bath of the same material, the liquidair interface of the bath deforms. They

72

deduced that the ruptures of the air film are related to the liquid viscosity and

73

the impact velocity. Recently, few experimental studies have been performed

74

to observe the interfacial gas layer during the impact of droplet on immiscible

75

liquid film [28, 27, 29]. Tang et al. [28] illustrated the temporal dynamics of the 4

76

gas layer for droplet bouncing on the liquid surface using high-speed imaging

77

and color interferometry. They measured the thickness of the gas layer under

78

several conditions and demonstrated the importance of the liquid film on droplet

79

dynamics. Che and Matar [29] observed that at W e = 3.84, the droplet im-

80

pacting on a liquid surface bounces several times due to the gas layer separating

81

the droplet and the immiscible liquid. Experimental studied demonstrated that

82

due to the microscopic nature of the entrapped gas layer, characterizing the

83

dynamics of air entrapment, measuring its thickness and interface morphology

84

is quite challenging [28, 29]. Hence, additional comprehensive investigations are

85

required.

86

For this process to be understood more fully, the goal of the present study

87

is using numerical analysis to investigate the profile of liquid film and droplet

88

during impingement, measuring the air thickness and capturing the phenomena

89

during the droplet impact dynamics. Although there are a significant number of

90

studies in the literature on two-phase flows [30, 31, 32, 33], much less attention

91

has been devoted to three-phase flows [34, 35, 36]. This is mostly due to the

92

complexity of modeling these systems because of their phase interactions and the

93

potential for triple-junctions to be present. To model the interaction between

94

different phases in three-phase flow, Smith et al. [37] also proposed a method to

95

decompose the surface tension coefficient. Another method is using the volume

96

of fluid (VOF). In the current study, a volume of fluid (VOF)-based numerical

97

tool for modeling the three-phase flow system was developed and validated using

98

the analytical solution. Temporal characterization of the gas dynamics and fluid

99

interface for different film thicknesses and impact velocities were investigated.

100

Meanwhile, we characterized the role of the air layer on the behaviors of both the

101

droplet and the liquid film using the two dimensionless parameters: normalized

102

film thickness h∗ =

103

and σ are droplet density, impact velocity, droplet initial diameter, and surface

104

tension, respectively, and h is the thickness of the liquid film.

h D0 ,

and Weber number W e =

5

2 ρw Uw D0 , σ

where ρw , Uw , D0 ,

105

2. Numerical methodology

106

Among different surface tracking methods that have been reported in the lit-

107

erature for computations of flows with moving interfaces [38, 39, 40, 41], the vol-

108

ume of fluid (VOF) method has proven to capture interacial changes accurately

109

due to its ability to conserve mass and handle sharp changes at the interface. In

110

this work, the VOF method combined with the continuum surface force (CSF)

111

model [42] was implemented in the open-source platform OpenFOAM [43]. Be-

112

cause of the symmetrical nature of this problem, a 2 dimensional axyisymmetric

113

coordinate system (r, z) was used to model the flow dynamics, where r and z

114

correspond to the radial and axial coordinates, respectively.

115

2.1. Governing equations

116

117

For an incompressible, newtonian, laminar flow, the equations governing the physics are given as follows, → − ∇· U =0

(1)

  → − → − → −T  → −→ − ∂ρ U ∇ U + ∇ U + ∇ · ρ U U = −∇P + ∇ · [µ( ) + ρ · g + Fs ∂t 2

(2)

118

where ρ, µ, U and P are the density, dynamic viscosity, velocity field, and

119

pressure of the flow, respectively. The surface tension force is denoted by Fs ,

120

which is computed at the interface region based on the continuum surface force

121

(CSF) [42]. This force is computed as Fs = σκˆ nδ , where σ is the surface tension

122

coefficient and κ is the interface curvature defined as −∇ · n ˆ , where n ˆ is a unit

123

surface vector.

124

In order to capture the dynamics of each phase, the volume of fluid (VOF)

125

method was used for tracking the interface, which is expressed by the volume

126

fraction αi , in which a value of i was assigned to each phase. The value of the

127

volume fraction indicator αi ranges between 0 and 1, where αi = 0 corresponds

128

to cells filled with gas, αi = 1 corresponds to cells filled with liquid, and 6

129

0 < αi < 1 represents the liquid/gas interface. A schematic of the volume

130

fraction distribution for a two-phase flow is shown in Figure 1a. The evolution

131

of the interface was computed using the following transport equation, → →  −  − ∂ραi + ∇ · U αi + ∇ · U r αi (1 − αi ) = 0 ∂t

(3)

132

where Ur is an artificial compression relative velocity at the free surface that

133

is applied normally into the interface to counteract the numerical diffusion and

134

maintain a sharp interface defined as, → − ∇αi U r = (min (Cαi |U | , max (|U |))) |∇αi |

135

(4)

where Cαi is the compression coefficient, which determines the degree of

136

compression.

137

2.2. Numerical method in three-phase flows

138

Since a three-phase flow system of a water droplet, an immiscible liquid film,

139

and the surrounding air was the focus of the current work, a multiphase flow

140

solver was required for producing a numerical solution.

141

142

143

For a three-phase immiscible flow system (water–oil–air), the volume fraction distribution is shown in Figure 1b. The volume fraction of each phase is n P distributed according to the transport equation, provided that αi = 1, where i=1

144

n is the the number of phases. In this study, n = 3 as we have three different

145

immiscible flows. In regions where all three phases exist, the volume fraction

146

is distributed as shown in Figure 1b. In the case of three-phase flows, all vari-

147

ables of fluid properties such as density and viscosity are defined by their phase

148

fraction and their values of ρi and µi , which define the density and viscosity

149

of each phase, respectively. Accordingly, the physical properties of fluids are

150

calculated according to the average weight of the volume fraction of each fluid.

151

152

The physical properties of fluids are obtained as a volumetric mixture value, NP =3 NP =3 ρ= αi ρi and µ = αi µi where N is the number of phases. i=1

i=1

7

153

In regions where all three phases exist (triple points), the complexity of three-

154

phase modeling is much higher than for two-phase modeling. Simultaneous

155

modeling of three-phase interaction is challenging and there are few studies

156

carried out to model the interactions [37, 36]. For instance, Smith et al. [37]

157

proposed a method to decompose the surface tension coefficient based on the

158

binary surface tension between phases.

(a)

(b)

Figure 1: Schematic of volume fraction distribution in VOF method for (a) two-phase: air αa )/water (αw ) (b) three immiscible phase: air (αa )/water (αw )/oil (αo ). Values in each cell corresponds to the volume fraction contained in the cell

159

In the current work, the surface normal and surface curvature are determined

162

using the Dirac delta function [35]. In the case of three-phase flows, this method 3 P 3 P calculated the surface tension force with Fs = σij κij δij where σij and κij

163

and j, respectively. In fact, all the domain follows the two-phase flow algorithm

164

except the surface tension at the triple region (where all three phases meet)

165

which is the summation of all binary interfacial tension force. The Dirac delta

166

function as implemented was,

160

161

i=1 j6=1

are the physical interfacial tension and the curvature between the two phases of i

Fs = − 167

3 X 3 X i=1 j6=1

σij (∇ · (

αj ∇αi − αi ∇αj ))(αj ∇αi − αi ∇αj ) |αj ∇αi − αi ∇αj |

(5)

The governing equations were then discretized by using the finite volume

8

168

method and following the schemes detailed in Table 1. Finally, the simulation

169

of the 2D axisymmetric geometry was performed in parallel using the PIM-

170

PLE algorithm, which is a combination of the PISO (Pressure Implicit with

171

Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked

172

Equations) algorithms. Table 1: Discretization schemes

Term ∂ ∂t

173

Discretization scheme

(ρU )

Euler

∇ · (ρU U )

Gauss Limited Linear V1 (TVD)

∇ · (αU )

Gauss VanLeer

∇ · (Ur α (1 − α))

Gauss interfaceCompression

2.3. Adaptive time step To optimize the computational time, an adaptive time-step control feature was applied. The time-step control feature also ensured the stability of the solution procedure as it was adjusted based on the Courant Friedrichs Lewy number (Co =

ui ∆t ∆x )

at the beginning of the time iteration loop. Using the

values of the velocity of the phase fractions and ∆t from previous time step, the maximum local Courant number (Co0 ) was calculated and a new time-step was iteratively initiated.     Comax Comax ∆t = min ∆t0 ; 1 + λ1 ∆t0 ; λ2 ∆t0 ; ∆tmax Co0 Co0

(6)

174

where Comax and ∆tmax were prescribed values for the Courant number and

175

time step, respectively. Moreover, in the beginning of the simulation at a very

176

small initial time-step (∆tini ), an intermediate time-step value was calculated

177

by using the formula ∆t∗ini = min



Comax ∆tini ; ∆tmax Co0 9



(7)

178

Then, this intermediate value (∆t∗ini ) was used as ∆t0 in the new time-step

179

calculation, providing a Co0 as close to the maximum Courant number. The

180

size of time steps varied during the calculations because of this Courant value

181

maximization, which allowed the time step to be smoothly adjusted.

182

2.4. Validation of the three-phase model

183

In this section, following the method of Boyer et al. [36], the model was

184

validated for an equilibrium droplet lens problem, which is a well-suited case for

185

a three-phase flow model verification. The droplet lens problem is the situation

186

of a spherical water droplet spreading at the middle of a straight surface between

187

two immiscible fluids (oil and air). A schematic of the initial conditions are

188

shown in Figure 2a. The balance of interfacial forces leads the droplet to reach

189

a final equilibrium condition (Figure 2b).

(a)

(b)

Figure 2: Schematic of droplet spreading between two immiscible phases (a) Initial condition (b) Equilibrium condition

190

In this model, we used a 1 mm droplet in the initial spherical position,

191

and neglected the effects of gravity. Therefore, the interfacial surface tension

192

force between phases was entirely responsible for the final configuration of the

193

droplet. At the equilibrium state, the final diameter of the droplet D (shown in

194

Figure 2b) and the contact angle of each fluid were obtained according to the

195

geometry and mathematical balancing of the interfacial tensions. The details of

196

calculating equilibrium diameter and contact angle are presented in Appendix.

10

197

To assess the accuracy of this model compared to the analytical solution,

198

two cases for different interfacial surface tension were simulated [33]. All details

199

of quantitative and qualitative comparison of validation results are illustrated in

200

Appendix. The results of the final droplet diameter and contact angles of phases

201

were compared to that obtained in the analytical model and it was observed that

202

they were in good agreements.

203

2.5. Computational set-up

204

A 2D axisymmetric geometry was performed for all cases. The droplet was

205

located in the center of the computational domain, which was sized about 5D0 ×

206

10D0 to make sure that the hydrodynamics of the droplet were not affected by

207

the domain boundaries (Figure 3). For high liquid film thickness (h∗ > 1),

208

since the liquid film deformation is more significant, the computational domain

209

is expanded to 10D0 × 10D0 . As the thickness of the air layer is on the order of

210

a few micrometers, the size of the cells in the computational domain is around

211

1 micrometer.

Figure 3: Schematic of computational domain of the numerical model

11

212

A no-slip boundary condition applied to the bottom surface of the geometry

213

with the gradient effect of the pressure is considered to be zero. Since the top

214

surface is subjected to the atmosphere, the boundary conditions of atmospheric

215

pressure with no gradient effect of velocity has been considered.

216

Regarding the fluids properties, a water droplet with a diameter of D0 =

217

2mm and a perfluorinated liquid (Dupont Krytox 103) were used as the droplet

218

and liquid film materials for this study, since this system is used in various

219

applications. The reason behind using this liquid film is due to its wide range

220

applications specially in the lubricant impregnated surfaces (LIS). This system

221

also has the simplifying advantage of low adhesion of the water droplet to the

222

surface. The physical properties of the droplet and the lubricant are presented

223

in Table 3. Table 2: Physical properties of fluids

Fluids

Density(kg /m3 )

Viscosity(P a. s)

Surface tension(N /m )

Dupont Krytox103

1900

0.15

0.017

Water

1000

0.001

0.07

224

For droplet impact on liquid surface, according to what have been observed

225

in experimental study, the presence of the liquid film influences the behavior of

226

the droplet dynamics and air layer since the liquid film is able to deform during

227

the impact process [28, 29]. In this study, we carried out a numerical study of

228

water droplet impact on a liquid film for film different thicknesses and impact

229

velocities, to identify this air layer thickness (δ) quantitatively. Additionally,

230

to better understand how the droplet and liquid film behave during the im-

231

pingement, we used two green and red points on the droplets and liquid film,

232

as shown in Figure 4, that better represent the similarities and differences for a

233

visual comparison.

12

Figure 4: Schematic of droplet impact on a liquid film of lubricant

234

In the next sections, we discuss comprehensively the effects that impact

235

velocity and film thickness had on the mechanism of the impact events.

236

3. Results and discussions

237

3.1. I. Role of the droplet impact velocity on the air layer

238

To identify the behavior of the droplet and the liquid film, the impact of

239

a water droplet on a liquid film with a thickness of h∗ = 0.0025 was studied

240

at three different Weber numbers (W e = 1.5, 5, and 10). The liquid film can

241

be treated as the thin film surface (h∗ < 1) [44, 28] which can be suitable

242

candidate for LIS applications [45, 46]. For a droplet at W e = 1.5, numerical

243

and experimental studies have shown the presence of an air layer between the

244

droplet and a solid surface [12, 17]. Analogously, we found that an air layers was

245

entrapped during the impact of a droplet on a liquid surface, which results in the

246

droplet bouncing without touching the liquid film. At low Weber number, the

247

base surface (liquid film) was slightly affected by the droplet inertia, however

248

the droplet follows the same hydrodynamics as was observed on solid surface

249

(i.e. spreading, rebounding, bouncing). Figure 5 illustrates this situation and

250

compares the results from the numerical modelling (droplet impact on thin

251

liquid film) performed in this study to experimental results (droplet impact on

252

solid surface) conducted by De ruiter [12].

13

Figure 5: Comparison of droplet impact at W e = 1.5 (a) On a solid surface, experimental results [12] and (b) On a thin liquid film surface (h∗ = 0.0025), numerical results

253

It can be seen that although both the droplet and liquid film interface was

254

affected, droplet hydrodynamics was similar to the case of droplet impact at

255

low Weber number on a solid surface in the experimental study.

256

As discussed, the liquid film was deformed during the impingement process.

257

To investigate the deformation in the liquid film and the droplet, the profile of

258

the central points of both the water droplet and the liquid film were labelled

259

to clearly identify their behavior during the impingement. Figure 6 illustrated

260

that the liquid film (lubricant) was deformed as the droplet started moving

261

downward, and the air cushion was generated beneath the droplet. The velocity

262

vector field in fluids at few millisecond of impingement demonstrated that how

263

the thin liquid film would react to droplet impact and the squeezed air can

264

deform the liquid film. It can be seen that the central point of liquid film profile

265

at low Weber number of 1.5 was almost constant after the first deformation of

14

266

the liquid film occurred.

Figure 6: Profile of water droplet and liquid film interface at the center of impact over time (W e = 1.5 and h∗ = 0.0025)

267

To investigate the effect of Weber number on droplet and the liquid film

268

hydrodynamics, impacts at two Weber numbers of W e = 5 and 10 were also

269

modeled. Since the droplet detached from the solid surface by W e ' 7 [12],

270

it was expected that the droplet would display spontaneous rebounding from

271

the liquid surface below that threshold. However, the presence of a liquid film

272

affected the impingement hydrodynamics even at a very small liquid film thick-

273

ness. By increasing the Weber number to W e = 5, droplet spread over a layer of

274

entrapped air beneath the droplet, as shown in Figure 7. The droplet reached

275

its maximum diameter (t ' 3 ms ) and retracted due to the surface tension

276

until it bounced from the liquid surface (t ' 10 ms ).

Figure 7: Evolution of water droplet impact on a thin liquid film (W e = 5 and h∗ = 0.0025)

15

277

The air cushion can produce a dimple-like shape in the water droplet as

278

reported for the impact of a droplet on a solid surface [11, 27]. The interface

279

profile of both the droplet and liquid film changed significantly during the few

280

milliseconds of droplet impact, as shown in Figure 8. In fact, during the spread-

281

ing stage, the pressure build-up would change both the interface profile of water

282

droplet and liquid film, even for a thin liquid film. This would indicate that the

283

deformation occur in both water droplet and liquid film, especially during the

284

few milliseconds of impingement until it spread to maximum diameter. After

285

that, during the rebounding stage, more phenomena happen which would be

286

illustrated in quantitative details.

Figure 8: Profile of the water droplet and liquid film interface (W e = 5 and h∗ = 0.0025)

287

The profile of central points of the droplet and the liquid film was measured

288

over time for the impact at W e = 5 (Figure 9). A significant deformation in the

289

profile of both the droplet and the liquid film was observed. The deformation

290

during the spreading stage until t ' 4 ms remained constant.

291

It is found that an air bubble is formed within the droplet (Figure 7) during

292

the retraction process (t ' 4 ms to t ' 5 ms). Due to the kinetic energy of

293

droplet, droplet deformed into the pyramid shape during the spreading stage due

294

to the capillary waves. Accordingly, the air cavity was formed at the center of the

295

droplet which results in air bubble entrapment within the droplet that has been

16

296

also observed in experimental studies, such as impact on a superhydrophobic

297

surface, soft and viscoelastic surfaces [47, 48, 49]. This air bubble was not formed

298

at W e = 1.5, since the droplet deformation was not large enough to generate

299

the bubble. Then, the bubble was ejected from the droplet (t ' 7.5 ms), and,

300

consequently, the height of droplet profile increased slightly (at t ' 8 ms) as the

301

liquid film height decreased. Finally, droplet bounced completely at t ' 10 ms.

Figure 9: Profile of the water droplet and the liquid film interface at the center of impact over time (W e = 5 and h∗ = 0.0025)

302

Further increase of the Weber number to W e = 10 was also investigated.

303

The temporal evolution of the shape of the water droplet and thin liquid film

304

is shown in Figure 10. Despite the failure of the droplet to bounce from a solid

305

surface at Weber numbers greater than around 7, we found that the bounc-

306

ing behavior of a water droplet from a liquid surface to be different. Figure

307

10 illustrates the qualitative comparison of the numerically calculated results

308

of droplet impact on a thin liquid film with the time lapse images from the

309

experimental droplet impacts performed by Hao et al. [45]. There was good

310

agreement in the hydrodynamics of the droplet between the numerical model

311

and the experimental results.

17

Figure 10: Comparison of water droplet impact on a thin liquid film (a) Experimental results [45] and (b) Side-view of numerical results (W e = 10 and h∗ = 0.0025)

312

The spreading factor, defined as the ratio of maximum droplet contact length

313

(Dmax ) to the initial droplet diameter (D0 ), was quantified and validated against

314

experimental results [45]. Figure 11 illustrates the result that the normalized

315

droplet spreading length in the numerical model agreed well with the experi-

18

316

mental results.

Figure 11: Comparison of numerical study of spreading factor (Dmax /D0 ) of droplet impact at W e = 10 on a thin liquid film of h∗ = 0.0025 with experimental results [45]

317

The spreading factor alone cannot describe all of the details behind the

318

droplet impact behavior. For this reason, we characterized the deformation of

319

the central point of the droplet and the liquid film over time. For W e = 10

320

we could subdivide the droplet bouncing mechanism into four distinct sections

321

to better clarify the process, as displayed in Figure 12. In section I, there

322

was a slight increase in the droplet height and a decrease in the liquid film

323

height due to the dimple formation, which was also seen in the lower Weber

324

numbers (W e = 1.5 and 5). After this region, the droplet started spreading until

325

t ' 4 ms and no significant deformation in the liquids interface was seen during

326

the spreading process. Once the water droplet began to rebound, the liquid

327

film moved downwards (section II). In section III, the droplet height increased

328

during the rebounding process. In fact, at the start of the droplet rebounding, a

329

small bubble was entrapped within the water droplet. It was then ejected from

330

the water droplet, which led to an increment in droplet height and a decrease

331

in the liquid film height. Finally, the effect of the bubble on the fluid height

332

dissipated and the droplet bounced from the liquid surface in section IV. 19

Figure 12: Profile of the water droplet and liquid film interface at the center of impact over time (W e = 10 and h∗ = 0.0025)

The interface profile of droplet and liquid film is plotted for W e = 10 and

333

∗

334

h = 0.0025 in Figure 13. Similarly, both of the droplet and liquid film interface

335

were changing upon the impact.

336

We also calculated the thickness of the entrapped air layer (δ as shown in

337

Figure 3) between the droplet and the liquid film for models run at different

338

Weber numbers, and plotted them in Figure 14a. The air layer was maintained

339

during almost all of the impingement process at the lowest Weber number tested,

340

W e = 1.5, while considerable deformation of this layer was observed at the

341

higher Weber numbers. Additionally, the effect of the inertia force in the W e =

342

5 and W e = 10 case resulted in a squeezed air layer, especially during the

343

rebounding stage, except the rise in some cases because of the ejecting bubble

344

inside the droplet.

345

A scaling law for the dimensionless height of dimple formation for impact

346

on a solid surface has been suggested by Bouwhuis et al. [50], Hd /R ∼ St−2/3 ,

347

where St is the Stokes number which is defined by St = ρl RU/µg . For We= 5

348

and 10, it is observed that the dimple height on thin liquid film is slightly lower

349

(∼ 3 − 5%) due to the damping caused by the thin liquid film comparing to the 20

350

predicted height for impact on a solid surface [50].

Figure 13: Profile of the water droplet and the liquid film interface (W e = 10 and h∗ = 0.0025)

Figure 14: Evolution of air layer thickness for different Weber numbers

351

It can be observed in all cases that although the deformation of the very thin

352

liquid film was largely restricted by the solid surface, the presence of the liquid

353

film and its movement during impingement substantially changed the conven-

354

tional mechanism of impact on a solid surface and increased the probability of

355

the droplet bouncing.

21

356

3.2. II. Role of the liquid film thickness on the air layer

357

In this part, we looked at the effect of the liquid film thickness at the same

358

Weber number (W e = 5) for a range of three film thicknesses. It should be

359

noted that different classifications of the liquid surface have been used for dif-

360

ferent thicknesses, as they behave differently [8]. For instance, h∗ < 0.1 [51] is

361

recommended as a thin liquid film, while several investigations suggest that thin

362

films are associated with h∗ < 1. In this study, h∗ = 1 was used as the thresh-

363

old to distinguish the liquid thin film regime from the liquid pool regime. The

364

results for liquid thin film of (h∗ = 0.0025) for W e = 5 were demonstrated in

365

previous section. The effect of the liquid film thickness on a thin film (h∗ = 0.1)

366

and liquid pool film (h∗ = 2.5) was also studied. Although the case of h∗ = 0.1

367

belongs to the group of liquid thin film (h∗ < 1), this treatment could have a

368

significant impact on the hydrodynamic changes of the liquid and the gas. For

369

the case of h∗ = 0.1, the time lapses of the impact process are displayed in Fig-

370

ure 15. The process of the droplet rebounding from the liquid surface occurred

371

due to the entrapped air beneath the droplet, which was a similar mechanism

372

to what has been explained in the previous section.

Figure 15: Evolution of water droplet impact on a thin liquid film surface (W e = 5 and h∗ = 0.1)

373

The spreading factor (Dmax /D0 ) of the h∗ = 0.1 case was similar to the

374

case of a thin liquid film of h∗ = 0.0025 (as shown in Figure 20). The profile

375

of the droplet and the liquid film at the center point of fluids were calculated

376

and are shown in Figure 16. Both the droplet and the liquid film experienced 22

377

deformation through the impingement process. The plot of interface profile in

378

Figure 17 is also provided for the W e = 5 and h∗ = 0.1). It was observed that

379

although the thickness of liquid film is in the thin film regime, the profile of the

380

liquid film interface was tending to deform more significantly during the initial

381

impingement process comparing to the h∗ = 0.0025.

Figure 16: Profile of water droplet and liquid film interface at the center of impact over time (W e = 5 and h∗ = 0.1)

Figure 17: Profile of the water droplet and the liquid film interface (W e = 5 and h∗ = 0.1)

23

382

For the case of a deep liquid pool of liquid film where h∗ = 2.5 at W e = 5,

383

the evolution of the droplet impact is shown in Figure 18. A more significant

384

deformation in the deep liquid pool during the impact process has been observed.

385

The air layer was significantly squeezed during the impingement while the large

386

deformation in the liquid film occurred. In fact, at this Weber number, the large

387

pool thickness was less constrained by the solid surface than was the thin film,

388

and, consequently, the liquid surface adapted in shape to the droplet impact.

Figure 18: Evolution of water droplet impact on a deep liquid film surface (W e = 5 and h∗ = 2.5)

24

389

In Figure 20, the spreading factor of the droplet indicated that the droplet

390

experienced less spreading in the case of the liquid pool compared to impinge-

391

ment at lower values of h∗ . This is because the liquid film was less restricted by

392

the solid surface, and, consequently, a large deformation of the liquid resisted

393

the droplet spreading outwards compared to the case of the thin liquid film.

394

The delay in the rebounding of the droplet from a deep liquid surface was also

395

attributed to a higher penetration of both the droplet and the liquid film, which

396

had been also observed experimentally by Tang et al. [28].

397

In the case of deep pool of liquid film, the profile of fluids interface are

398

presented in Figure 19. Since the inertia force is not high and the liquid film is

399

not constraint by the solid surface, it can be seen that the motion of liquid film

400

is alongside the water droplet while the air film is captivated within the droplet

401

and liquid film. Thus, the dimple-shape is less significant in this case.

Figure 19: Profile of the water droplet and the liquid film interface (W e = 5 and h∗ = 2.5)

25

Figure 20: Spreading factor over time at W e = 5 and for different liquid film thicknesses ( h∗ = 0.0025, h∗ = 0.1 and h∗ = 2.5)

402

In order to quantify and analyze the dynamics of the deep pool of an immis-

403

cible liquid film on droplet impingement, the variation of the droplet and liquid

404

film profile at the center points were evaluated (Figure 21). Both of the fluids

405

deformed during the impact process. The droplet deformed substantially and

406

rebounded after reaching its maximum spreading, and a small bubble formed

407

within the inside of the water droplet during the rebounding stage. A slight in-

408

crease in the height of the droplet and a decrease in the height of the liquid film

409

were demonstrated this incident according to eject of bubble inside the water

410

droplet. Finally, when the droplet detached from the liquid surface, the height

411

of both fluids increased.

26

Figure 21: Profile of the water droplet and the liquid film at the center of impact over time (W e = 5 and h∗ = 2.5)

412

To better understand the effects of the entrapped air layer and its behavior,

413

the air thickness (δ) was measured for the mentioned three liquid film thick-

414

nesses of h∗ = 0.0025, h∗ = 0.1 and h∗ = 2.5 at W e = 5, which is shown

415

in Figure 22a. It is interesting to note that there is a decrease in the air

416

layer thickness during the process of droplet approaching the surface (initial

417

impingement), and the trend of this thickness change between runs is that

418

δh=50µm < δh=200µm < δh=5 mm . It is also observed that the dimple height

419

in higher liquid film thicknesses would be larger comparing to the introduced

420

numerical correlation for the case of impact on solid surface [50]. This is likely

421

because the thin liquid film was constrained by the solid surface, and the air

422

layer was squeezed to a greater extent compared to the deep liquid pool case

423

during the spreading stage. Afterward this approach, the thickness of the in-

424

terfacial gas layer in the deep liquid pool decreased during the retracting stage.

425

This is due to the easily deformed liquid surface that was adaptive to the impact

426

during the droplet impingement.

27

Figure 22: Air layer thickness over time for different liquid film thicknesses (h∗ = 0.0025, h∗ = 0.1 and h∗ = 2.5)

427

Additionally, Tang et al.[28] observed a similar trend in their recent experi-

428

mental study using high-speed imaging and color interferometry. They showed

429

that the air layer thickness decreased almost linearly with time for both thin and

430

thick liquid layers. We found that the rebounding stage was dramatically dif-

431

ferent from the approaching stage. It is observed a substantial rapid increment

432

in the air layer thickness for a deep pool of liquid film, which started slowing

433

down after rebounding until the droplet detached from the liquid surface.

434

4. Conclusion

435

A three-phase flow solver has been developed to investigate the hydrody-

436

namics change of a liquid film (lubricant) and the entrapped air layer during

437

droplet impact and bouncing. The numerical model was validated against the

438

analytical model based on the balance of the interfacial forces. In order to un-

439

derstand and characterize the nature of the air layer underneath the droplet,

440

droplet impact on an immiscible liquid film at different film thicknesses and

441

Weber numbers have been studied.

442

It was observed that at low Weber number (W e = 1.5 and h∗ = 0.0025), 28

443

despite the liquid film deformation, the hydrodynamics change of the droplet on

444

a thin liquid film was similar to that on a solid surface which is caused by the

445

presence of the air layer. For higher Weber number (W e = 10), droplet sticks

446

to the solid surface whereas on the liquid film, droplet detaches due to the local

447

deformation of the liquid film interface caused by the entrapment of the air. For

448

the liquid film in the deep pool regime, liquid deformation was more significant

449

due to the diminished constraints of the solid surface on this thicker liquid film.

450

Furthermore, in some cases, an air bubble was formed within the droplet during

451

the retracting process, which also facilitated the droplet rebounding process.

452

The numerical results can provide valuable guidelines for exploring the behavior

453

of three-phase flow, especially in LIS.

454

5. Acknowledgment

455

The authors gratefully acknowledge the financial support from Natural Sci-

456

ences and Engineering Research Council of Canada (NSERC) and EU partner-

457

ship project Phobic2Ice. Computing resources were provided by Calcul Qubec

458

and Compute Canada

459

6. References

460

References

461

[1] A. McDonald, M. Lamontagne, C. Moreau, S. Chandra, Impact of plasma-

462

sprayed metal particles on hot and cold glass surfaces, Thin Solid Films

463

514 (1-2) (2006) 212–222.

464

[2] S. D. Aziz, S. Chandra, Impact, recoil and splashing of molten metal

465

droplets, International journal of heat and mass transfer 43 (16) (2000)

466

2841–2857.

467

[3] B. Derby, Inkjet printing of functional and structural materials: fluid prop-

468

erty requirements, feature stability, and resolution, Annual Review of Ma-

469

terials Research 40 (2010) 395–414. 29

470

471

472

473

474

475

476

477

478

479

480

481

[4] A. Davanlou, R. Kumar, Thermally induced collision of droplets in an immiscible outer fluid, Scientific reports 5 (2015) 9531. [5] F. Blanchette, T. P. Bigioni, Partial coalescence of drops at liquid interfaces, Nature Physics 2 (4) (2006) 254. [6] A.-B. Wang, C.-C. Chen, Splashing impact of a single drop onto very thin liquid films, Physics of fluids 12 (9) (2000) 2155–2158. [7] T. Gilet, J. W. Bush, Droplets bouncing on a wet, inclined surface, Physics of Fluids 24 (12) (2012) 122103. [8] C. Tropea, M. Marengo, The impact of drops on walls and films, Multiphase Science and Technology 11 (1). [9] S. Howison, J. Ockendon, J. Oliver, R. Purvis, F. Smith, Droplet impact on a thin fluid layer, Journal of Fluid Mechanics 542 (2005) 1–23.

482

[10] S. Chandra, C. Avedisian, On the collision of a droplet with a solid surface,

483

Proceedings of the Royal Society of London. Series A: Mathematical and

484

Physical Sciences 432 (1884) (1991) 13–41.

485

[11] J. de Ruiter, D. van den Ende, F. Mugele, Air cushioning in droplet impact.

486

ii. experimental characterization of the air film evolution, Physics of fluids

487

27 (1) (2015) 012105.

488

[12] J. De Ruiter, R. Lagraauw, D. Van Den Ende, F. Mugele, Wettability-

489

independent bouncing on flat surfaces mediated by thin air films, Nature

490

physics 11 (1) (2015) 48.

491

[13] J. de Ruiter, R. Lagraauw, F. Mugele, D. van den Ende, Bouncing on thin

492

air: how squeeze forces in the air film during non-wetting droplet bouncing

493

lead to momentum transfer and dissipation, Journal of fluid mechanics 776

494

(2015) 531–567.

495

[14] J. M. Kolinski, L. Mahadevan, S. Rubinstein, Drops can bounce from per-

496

fectly hydrophilic surfaces, EPL (Europhysics Letters) 108 (2) (2014) 24001. 30

497

498

499

500

[15] C. Josserand, S. T. Thoroddsen, Drop impact on a solid surface, Annual review of fluid mechanics 48 (2016) 365–391. [16] V. Mehdi-Nejad, J. Mostaghimi, S. Chandra, Air bubble entrapment under an impacting droplet, Physics of fluids 15 (1) (2003) 173–183.

501

[17] H. Shetabivash, A. Dolatabadi, Numerical investigation of air mediated

502

droplet bouncing on flat surfaces, AIP Advances 7 (9) (2017) 095003.

503

[18] F. Yeganehdoust, I. Karimfazli, A. Dolatabadi, Numerical analysis of

504

droplet impact on a smooth slippery surface, in: ASME 2018 5th Joint

505

US-European Fluids Engineering Division Summer Meeting, American So-

506

ciety of Mechanical Engineers, 2018, pp. V001T15A005–V001T15A005.

507

[19] T.-S. Wong, S. H. Kang, S. K. Tang, E. J. Smythe, B. D. Hatton,

508

A. Grinthal, J. Aizenberg, Bioinspired self-repairing slippery surfaces with

509

pressure-stable omniphobicity, Nature 477 (7365) (2011) 443.

510

511

[20] C. Lee, H. Kim, Y. Nam, Drop impact dynamics on oil-infused nanostructured surfaces, Langmuir 30 (28) (2014) 8400–8407.

512

[21] M. Cao, D. Guo, C. Yu, K. Li, M. Liu, L. Jiang, Water-repellent properties

513

of superhydrophobic and lubricant-infused slippery surfaces: A brief study

514

on the functions and applications, ACS applied materials & interfaces 8 (6)

515

(2015) 3615–3623.

516

[22] F. Marcotte, G.-J. Michon, T. S´eon, C. Josserand, Ejecta, corolla, and

517

splashes from drop impacts on viscous fluids, Physical review letters 122 (1)

518

(2019) 014501.

519

[23] H. M. Kittel, I. V. Roisman, C. Tropea, Splash of a drop impacting onto

520

a solid substrate wetted by a thin film of another liquid, Physical Review

521

Fluids 3 (7) (2018) 073601.

522

[24] H. M. Kittel, E. Alam, I. V. Roisman, C. Tropea, T. Gambaryan-Roisman,

523

Splashing of a newtonian drop impacted onto a solid substrate coated by a 31

524

thin soft layer, Colloids and Surfaces A: Physicochemical and Engineering

525

Aspects 553 (2018) 89–96.

526

527

[25] N. Chen, H. Chen, A. Amirfazli, Drop impact onto a thin film: Miscibility effect, Physics of Fluids 29 (9) (2017) 092106.

528

[26] X. Tang, A. Saha, C. K. Law, C. Sun, Nonmonotonic response of drop

529

impacting on liquid film: mechanism and scaling, Soft Matter 12 (20) (2016)

530

4521–4529.

531

532

[27] T. Tran, H. de Maleprade, C. Sun, D. Lohse, Air entrainment during impact of droplets on liquid surfaces, Journal of fluid mechanics 726.

533

[28] X. Tang, A. Saha, C. K. Law, C. Sun, Bouncing drop on liquid film: Dy-

534

namics of interfacial gas layer, Physics of Fluids 31 (1) (2019) 013304.

535

[29] Z. Che, O. K. Matar, Impact of droplets on immiscible liquid films, Soft

536

matter 14 (9) (2018) 1540–1551.

537

[30] R. Attarzadeh, A. Dolatabadi, Coalescence-induced jumping of micro-

538

droplets on heterogeneous superhydrophobic surfaces, Physics of Fluids

539

29 (1) (2017) 012104.

540

[31] J. Klostermann, K. Schaake, R. Schwarze, Numerical simulation of a single

541

rising bubble by vof with surface compression, International Journal for

542

Numerical Methods in Fluids 71 (8) (2013) 960–982.

543

[32] G. Son, Efficient implementation of a coupled level-set and volume-of-fluid

544

method for three-dimensional incompressible two-phase flows, Numerical

545

Heat Transfer: Part B: Fundamentals 43 (6) (2003) 549–565.

546

[33] M. Tembely, R. Attarzadeh, A. Dolatabadi, On the numerical modeling of

547

supercooled micro-droplet impact and freezing on superhydrophobic sur-

548

faces, International Journal of Heat and Mass Transfer 127 (2018) 193–202.

32

549

[34] N. Tofighi, M. Yildiz, Numerical simulation of single droplet dynamics in

550

three-phase flows using isph, Computers & Mathematics with Applications

551

66 (4) (2013) 525–536.

552

[35] N. Weber, P. Beckstein, W. Herreman, G. M. Horstmann, C. Nore, F. Ste-

553

fani, T. Weier, Sloshing instability and electrolyte layer rupture in liquid

554

metal batteries, Physics of Fluids 29 (5) (2017) 054101.

555

[36] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar, M. Quintard, Cahn–

556

hilliard/navier–stokes model for the simulation of three-phase flows, Trans-

557

port in Porous Media 82 (3) (2010) 463–483.

558

[37] K. A. Smith, F. J. Solis, D. Chopp, A projection method for motion of

559

triple junctions by level sets, Interfaces and free boundaries 4 (3) (2002)

560

263–276.

561

[38] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi,

562

W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the

563

computations of multiphase flow, Journal of computational physics 169 (2)

564

(2001) 708–759.

565

[39] C. W. Hirt, B. D. Nichols, Volume of fluid (vof) method for the dynamics

566

of free boundaries, Journal of computational physics 39 (1) (1981) 201–225.

567

[40] F. Yeganehdoust, M. Yaghoubi, H. Emdad, M. Ordoubadi, Numerical study

568

of multiphase droplet dynamics and contact angles by smoothed particle

569

hydrodynamics, Applied Mathematical Modelling 40 (19-20) (2016) 8493–

570

8512.

571

[41] M. Ordoubadi, M. Yaghoubi, F. Yeganehdoust, Surface tension simulation

572

of free surface flows using smoothed particle hydrodynamics, Scientia Iran-

573

ica 24 (4) (2017) 2019–2033.

574

575

[42] J. U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of computational physics 100 (2) (1992) 335–354.

33

576

[43] Openfoam, http://openfoam.com (2011).

577

[44] G. Liang, I. Mudawar, Review of mass and momentum interactions dur-

578

ing drop impact on a liquid film, International Journal of Heat and Mass

579

Transfer 101 (2016) 577–599.

580

[45] C. Hao, J. Li, Y. Liu, X. Zhou, Y. Liu, R. Liu, L. Che, W. Zhou, D. Sun,

581

L. Li, et al., Superhydrophobic-like tunable droplet bouncing on slippery

582

liquid interfaces, Nature communications 6 (2015) 7986.

583

[46] J. D. Smith, R. Dhiman, S. Anand, E. Reza-Garduno, R. E. Cohen, G. H.

584

McKinley, K. K. Varanasi, Droplet mobility on lubricant-impregnated sur-

585

faces, Soft Matter 9 (6) (2013) 1772–1780.

586

[47] L. Chen, J. Wu, Z. Li, S. Yao, Evolution of entrapped air under bouncing

587

droplets on viscoelastic surfaces, Colloids and Surfaces A: Physicochemical

588

and Engineering Aspects 384 (1-3) (2011) 726–732.

589

590

[48] L. Chen, E. Bonaccurso, P. Deng, H. Zhang, Droplet impact on soft viscoelastic surfaces, Physical Review E 94 (6) (2016) 063117.

591

[49] S. Lin, B. Zhao, S. Zou, J. Guo, Z. Wei, L. Chen, Impact of viscous droplets

592

on different wettable surfaces: Impact phenomena, the maximum spreading

593

factor, spreading time and post-impact oscillation, Journal of colloid and

594

interface science 516 (2018) 86–97.

595

[50] W. Bouwhuis, R. C. van der Veen, T. Tran, D. L. Keij, K. G. Winkels, I. R.

596

Peters, D. van der Meer, C. Sun, J. H. Snoeijer, D. Lohse, Maximal air

597

bubble entrainment at liquid-drop impact, Physical review letters 109 (26)

598

(2012) 264501.

599

[51] R. L. Vander Wal, G. M. Berger, S. D. Mozes, Droplets splashing upon films

600

of the same fluid of various depths, Experiments in fluids 40 (1) (2006) 33–

601

52.

34

602

Appendix

603

According to the Neumanns relation and the summation of contact angles

604

at a phase triple junction, the contact angles of each phase can be written as sin θw sin θa sin θo = = σoa σwo σwa

(8)

θw + θo + θa = 2π

(9)

where θw , θo and θa are the contact angles of the water, oil, and air phases, respectively. By balancing the volume of the droplet at the initial and final states (which contains two spherical caps), the final diameter of the droplet is given by D = D0

p 3

4/ a

605

where D0 is the initial diameter of the droplet and a is a coefficient that is a

606

function of the contact angles of two other phases, defined as 

 2 2 1 1 1 a= + − sin(π − θo ) tan(π − θo ) sin(π − θo ) tan(π − θo )   2 1 1 2 1 + − + sin(π − θa ) tan(π − θa ) sin(π − θa ) tan(π − θa )

(10)

607

Two cases were considered: Case 1 (σwa ; σoa ; σwo )=(1; 1; 1) and Case

608

2 (σwa ; σoa ; σwo )=(0.8; 1; 1.4). In both cases, the interfacial tensions were

609

scaled by the interfacial tension between oil and air, σoa . Parameters that were

610

compared included the final diameter of the droplet and the contact angle of

611

each phase at the triple line. In both cases 1 and 2, the interfacial tensions were

612

scaled by the interfacial tension between oil and air, σoa . The axisymmetric

613

numerical simulation was performed and the resulting states were compared

614

with the analytical solutions. Parameters that were compared included the final

615

diameter of the droplet and the contact angle of each phase at the triple line.

616

For a droplet with a diameter of 1 mm, the final diameters (D) for cases 1 and

617

2 can be analytically calculated to be 1.276 mm and 1.074 mm, respectively. In 35

618

our numerical model, the equilibrium diameters in cases of 1 and 2 reached 1.2

619

mm and 1.08 mm, respectively, giving errors of about 6% and 1%, respectively,

620

that can be treated as negligible. In order to assess the grid independence of the

621

results for a three-phase flow, numerical simulations were performed at different

622

levels of cell size and the error of final states of steady equilibrium diameter

623

of droplet comparing to the analytical results have been measured as shown in

624

Figure 23.

4.5 4 3.5 3 2.5

2

1.5

20

40

60

80

100

120

140

Figure 23: Error of numerical results for different level of refinement

625

For the numerical calculation performed at the cell size of 40 µm, the error

626

of the final diameter solution did not compatible with other cell sizes. For finer

627

cell sizes, following the regime of droplet oscillations, the final diameters reached

628

a size approximately equal to the analytical solution. This showed that cell sizes

629

smaller than 40 µm possessed grid independence, and were acceptable cell sizes

630

for use in the model.

631

The contact angle of phases at the stable conditions are shown for two cases

632

in Figure 24. The contact angles were measured and compared with the ana-

633

lytical force balance at the triple line, as summarized in Table 3.

36

Table 3: The numerical and analytical comparison of contact angle for a droplet lens, case 1: (σwa ; σoa ; σwo )=(1; 1; 1) and case 2: (σwa ; σoa ; σwo )=(0.8; 1; 1.4)

Contact angle Case

θw 1

θo 2

◦

1 ◦

Numerical

120

128

analytical

120◦

136◦

120

θa 2

◦

120◦

145

1 ◦

146◦

2 ◦

87◦

120◦

78◦

120

634

To assess the grid independence of the final steady results of contact angles,

635

three-phase flow numerical simulations were performed at different levels of

636

cell size.The errors of final states have been measured (as shown in Figure25)

637

comparing to the analyrical values.

638

It can be seen that the cell size of 40 µm brings significant error to the

639

model, especially in the second case. Increasing the number of cells resulted

640

in reducing the error of the contact angle for different phases comparing to the

641

analytical results.

37

Conflict of Interest and Authorship Conformation Form Please check the following as appropriate: 642

o

All authors have participated in (a) conception and design, or analysis and interpretation of the data; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version.

o

This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue.

o

The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript

o

The following authors have affiliations with organizations with direct or indirect financial interest in the subject matter discussed in the manuscript:

Author’s name Firoozeh Yeganehdoust Reza Attarzadeh Ida Karimfazli Ali Dolatabadi

Affiliation Concordia University Concordia University Concordia University Concordia University

Author’s statement 643 The contributions of each authors are summarized as follows,

-Firoozeh Yeganehdoust: Methodology, programming, investigation, visualization, validation, writing (first draft) and writing (review and editing)

-Reza Attarzadeh: Methodology and formal analysis

-Ida Karimfazli: Conceptualization, supervision, project administration and writing (review and editing)

-Ali Dolatabadi: Conceptualization, supervision and project administration and writing (review and editing)

(a)

(b)

(c)

Figure 24: Droplet lens at (a) t=0 ms, initial conditions, (b) t=40 ms, similar interfacial tensions (σwa ; σoa ; σwo )=(1; 1; 1), (c) t=40 ms, different interfacial tensions (σwa ; σoa ; σwo )=(0.8; 1; 1.4)

40

3

w

2.5

o a

2

1.5

1

0.5 20

40

60

80

100

120

140

(a)

w o a

10 1

10 0 20

40

60

80

100

120

140

(b) Figure 25: Error of predicting contact angle comparing to analytical results for (a) case 1: Same interfacial tensions, (b) case 2: Different interfacial tension

41