Analysis of Bubble Coalescence Dynamics and Postrupture

Subscriber access provided by READING UNIV

Article

Analysis of bubble coalescence dynamics and postrupture oscillation of capillary-held bubbles in water Yesenia Saavedra Moreno, Ghislain Bournival, and Seher Ata Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b03197 • Publication Date (Web): 21 Nov 2017 Downloaded from http://pubs.acs.org on November 26, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1

1

Analysis of bubble coalescence dynamics and post-

2

rupture oscillation of capillary-held bubbles in water

3

Yesenia Saavedra Moreno, Ghislain Bournival, and Seher Ata*

4

School of Mining Engineering, University of New South Wales, Sydney, New South Wales,

5

2052, Australia

6

Keywords: bubble coalescence; neck expansion; post-rupture oscillation; VOF method

7

Corresponding author

8

*

E-mail: [email protected]

9

10

ABSTRACT

11

This paper presents a three-dimensional computational study of coalescence dynamics of two

12

capillary-held air bubbles in water using the volume of fluid (VOF) method. The interface

13

motion of the newly formed bubble indicated that smaller initial separation distances resulted in

14

a slightly faster expansion of the neck. The velocity vectors showed that the inward motion of

15

the air at the contact point of the two bubbles favors the generation of small eddies at the top and

16

bottom part of the neck at the early stage of coalescence and bulges at the edges of the bubble at

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 39

2

17

a later stage. The release of free-surface energy drove the bubble contraction and expansion and

18

consequently the oscillation motion, which created differences in pressure across the bubble

19

interface. The computationally simulated bubbles followed the oscillatory motion observed

20

experimentally, but with lower damping constant and lower angular frequency.

21

1. INTRODUCTION

22

The coalescence of air bubbles is an important phenomenon in a number of industrial

23

processes that use bubbling in a liquid flow such as in froth flotation, wastewater treatment and

24

paper recycling. Bubble coalescence is initiated by the thinning of the liquid separating the two

25

bubbles followed by a rupture of the liquid film 1. Once coalescence occurs, a neck is formed

26

joining the two bubbles. As the neck opens, the surface area and the curvature of the two initial

27

bubbles decrease. The newly formed bubble expands horizontally and contracts vertically and

28

vice versa until a stable spherical shape is achieved. The release of free-surface energy due to

29

this coalescence process is imparted to the surrounding liquid as kinetic energy 2.

30

The effects of inorganic electrolytes and surfactants in inhibiting bubble coalescence have

31

been widely studied experimentally and theoretically. Many dissolved salts have a critical

32

concentration beyond which the likelihood of bubble coalescence is significantly reduced

33

Lessard and Zieminski 3 introduced a transition concentration for each electrolyte studied where

34

the probability of coalescence was 50%, using 100% as the coalescence of a bubble pair in pure

35

water. They found that electrolytes with valence combinations of 3-1 and 2-2 were the most

36

effective electrolytes with transition concentrations of 0.032 M and 0.035 M, respectively. These

37

findings were corroborated by Craig et al. 4, who observed that highly charged salts were more

38

effective as bubble coalescence inhibitors at lower transition concentrations. More recently,


3-6

.

Page 3 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

39

Yaminsky et al. 7 reported that the speed at which the bubbles were brought together played an

40

important role in bubble coalescence, finding that instantaneous coalescence of bubbles occurred

41

at greater approach speed (150 µm/s). However, they found that high concentration of NaCl

42

electrolyte increased film stability by immobilizing the interface and therefore inhibiting

43

coalescence. For non-ionic surfactants, it is believed that the adsorption at the liquid-air interface

44

is mainly driven by the hydrophobic forces

45

surface properties of the interface, causing the Gibbs-Marangoni phenomena, which is

46

commonly used to explain the effect of surface tension on bubble coalescence. However, surface

47

forces have also been recognized as an important contributor in preventing film rupture in some

48

non-ionic surfactant systems

49

increasing surfactant concentration and coalescence time until a maximum coalescence time was

50

reached. However, all surfactants investigated had a maximum coalescence time until no change

51

with increasing surfactant concentration was observed.

10

8,9

. The amount of surfactant adsorbed affects the

. Bournival et al.

11

observed a direct relationship between

52

High-speed cameras and laser technologies are major optical technique developments which

53

have been used as tools to improve understanding of bubble coalescence dynamics. Thoroddsen

54

et al.

55

approximately 2.4 mm in average diameter using a high speed imaging technique with a capture

56

rate of 1×106 frames/second (fps). The frame sequences taken allowed the authors to analyze and

57

quantify the motion of the neck and bubble surface shape. More recently, Bournival et al. 11 used

58

a high speed camera to study the coalescence dynamics of bubbles in non-ionic surfactant

59

solutions. The authors generated two bubbles attached to the tips of a pair of fine capillaries to

60

examine the relationship between coalescence time, defined as the time interval from the

61

moment at which the bubbles are in contact to the rupture of the liquid film separating the two

12

investigated the motion of the neck during the merging of two air bubbles of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 39

4

62

bubbles, and surfactant concentration. Coalescence and oscillation of the entire bubble were

63

captured and analyzed by tracking the newly formed bubble’s change in projected surface area

64

over time. A direct relationship between increasing surfactant concentration and coalescence

65

time was observed and the post-rupture dampening of the oscillation of the newly formed bubble

66

increased with surfactant concentration. The study suggested that surface elasticity partly

67

dampened the oscillation of the newly formed bubble while other factors like surface viscosity

68

may contribute to the dissipation of energy. The experimental work provided an analysis of the

69

surface area of the oscillated bubble in two dimensions, quantifying the expansion and

70

contraction of the coalesced bubble as a projected area. However, there are still questions on how

71

fluid dynamics properties, such as flow velocity and dynamic pressure, control the bubble

72

coalescence and post-rupture oscillation. The objective of the paper is to explain bubble

73

coalescence dynamics by characterizing the flow velocity field vectors and local dynamic

74

pressure inside the bubble and in the surrounding liquid across the bubble interface using

75

computer simulation.

76

Several numerical investigations have used computational fluid dynamics (CFD) as a tool for 13

77

studying multiphase flows

. Most of these studies use multiphase methods such as volume of

78

fluid (VOF) method 14-22, front tracking method 23,24 and level set method 25-27. The VOF method,

79

first introduced by Hirt and Nichols

80

volume of each fluid phase in a specific computational cell. Using this method, Tomiyama et al.

81

14

82

rising in a laminar flow. The computational results were in reasonable agreement with the

83

experimental data suggesting that VOF may predict the qualitative and quantitative behavior of

84

bubble motion under different conditions. Nevertheless, further development of the numerical

28

, tracks the liquid-air interface by defining a fractional

investigated the effect of fluid properties and surface tension on the motion of a single bubble


Page 5 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

85

model to consider a larger computational domain and a three-dimensional model is required for a

86

precise quantitative prediction of the bubble rising motion.

87

The study of coalescence dynamics of a bubble pair rising freely in a liquid phase using the 15,16,18

18

88

VOF method has been well documented

89

dimensional model of two co-axial free bubbles rising in a liquid phase and used the VOF

90

method to track the motion, shape, and velocity of the two bubbles during the interaction and

91

coalescence processes. Although the two-dimensional numerical study did not fully represent the

92

experimental conditions, the findings corroborated the experimental data available in the

93

literature such as Chen et al. 15. The study of Hasan and Zakaria 18 showed that the VOF method

94

may be considered an accurate tool for modeling bubble coalescence time and the subsequent

95

oscillation behavior. Some authors have carried out qualitative analysis of the interaction of a

96

pair of bubbles rising side by side

97

investigate the effect of approach and rise velocity on bubble coalescence, considering liquid

98

film drainage before the film rupture.

17,19,29,30

. Hasan and Zakaria

developed a two-

. Such computational studies were intended only to

99

These studies show that no research has numerically investigated the oscillation following the

100

merging of two captive bubbles. The present study provides insight into bubble coalescence

101

dynamics by mapping the flow velocity vectors inside the bubble and in the surrounding liquid,

102

and local dynamic pressure across the bubble interface. This objective was achieved by modeling

103

two air bubbles fixed to adjacent capillaries in three dimensions using the VOF method. The

104

developed model was validated by comparing the overall oscillatory motion of the bubbles with

105

experimental data produced in a well-controlled system. The analysis provided further

106

understanding of post-rupture bubble oscillation and the fluid dynamics properties regulating the

107

coalescence process.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 39

6

108

2. METHODOLOGY

109

2.1 Computational Solution Method

110

The VOF method developed by Hirt and Nichols

28

was used to simulate the coalescence of

111

two air bubbles fixed to adjacent capillaries in water. The computational model was implemented

112

in the commercial code ANSYS-Fluent v 17.2.

113

2.1.1 Governing equations

114

In the VOF method, both phases were assumed to be insoluble and all computational cells in

115

the domain were either occupied by the water or air phase or a combination of both 28. The water

116

phase was defined as the primary phase and the air phase as the secondary phase. A VOF

117

function F was defined to track the fractional volume of a particular fluid in any specific cell of

118

the computational domain at time t. If the cell was completely empty of air, F = 0; if the cell was

119

completely full of air phase, F = 1; and if the cell contained the interface between the water and

120

air phase, 0 < F < 1. Assuming that the mass of the fluid was preserved, the VOF method tracked

121

the water-air interface by solving a continuity equation for the secondary phase (air) given by: ∇·V=0

122 123

1

where V is the velocity vector of the air phase in the entire domain. A momentum equation was solved for both water and air phases by the following equation: ∂ρV V ∇ρV V∙V V - ∇p ρg g ∇∙μ ∇V V∇V VT F Fσv ∂t

2

124

where ρ and t represent the density and time; p, g, µ, Fσv, are pressure, gravity force, viscosity

125

and surface tension force per unit volume, respectively; and T is the matrix transpose operator.

126

The fluid properties, ρ and µ, were determined by the properties of each phase as:


Page 7 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7

127 128

ρ= Fair ρair +(1-Fair )ρwater

3

µ= Fair µair +(1-Fair )µwater

4

The surface tension force per unit volume (Fσv) was calculated by adopting the continuous surface force model developed by Brackbill et al. 31 and assuming a constant surface tension: Fσv = σkn

129 130

5

In eq 5, σ is the surface tension, k is the surface curvature and n is the normal vector to the interface. The surface curvature of the interface, k, was calculated by: k=

1 n ∇ |n|-∇n |n| |n|

6

131

The normal vector, n, was evaluated by the following equation, which took into account the

132

three-phase contact perimeter formed between the bubble interface (i.e. water and air) and the

133

capillary tube, known as the wall adhesion effect: n= nwall cos θeq +twall sin θeq

134 135

7

where nwall and twall are the unit vectors normal and tangent to the capillary tube respectively, and θeq is the equilibrium contact angle between the water-air interface and the capillary tube.

136

2.1.2 Geometry and mesh design

137

A rectangular structure of dimensions 18×30×8 mm with two cylindrical capillaries of 20 mm

138

long was created in ANSYS-Design modeler. The geometry used was sufficiently large with the

139

ratio of the average diameter of the two initial bubbles (Db,ave) to the domain diameter being less

140

than about 0.25 to eliminate any wall effects on the bubble coalescence dynamics 32. A schematic

141

representation of the geometry used in ANSYS-Fluent is presented in Figure 1. The two vertical



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 39

8

142

capillaries were chosen to have an outer diameter of 1.07 mm and inner diameter of 0.69 mm to

143

replicate the experimental conditions of Bournival et al.

144

offered a well-controlled experimental system, where the neck expansion and the bulging of the

145

extremities of bubbles were not as constrained by the capillary tubes as other configurations.

146

This capillary arrangement also allowed the (almost free) surface motion of the newly formed

147

bubble, minimizing the effect of the capillary tubes on the horizontal expansion of the coalesced

148

bubbles. Practically, the vertical capillaries effectively allowed the coating of the bubbles with

149

hydrophobic particles. Thus, that configuration offered practical advantages in a stirred system

150

and for the packing of particles on the surface of the bubbles 33-36.

11

, which for experimental purposes

151

The domain was discretized in ANSYS-Meshing using a multizone mesh method to generate

152

hexahedral meshes. A local mesh refinement was implemented by increasing the number of cells

153

per average diameter of the bubbles to improve the resolution of the interface curvature around

154

the bubbles and reduce the computational error of the interface tracking. The mesh refinement

155

was created in ANSYS-Fluent by marking a hexahedron with dimensions of 2.5Db,ave × 2.5Db,ave

156

× 2.5Db,ave, as depicted in Figure 2a. A case with 23 cells per 2 mm was compared against two

157

cases with mesh refinement of 46 and 92 cells per 2 mm. The mesh independence solution was

158

evaluated by extracting the mean velocity profiles at the cross section plane. A small difference

159

of 4.1% was observed in the mean velocity for the finest mesh (92 cells per 2 mm) with respect

160

to the case with no mesh refinement (23 cells per 2 mm), indicating a mesh independent solution.

161

Consequently, for computational efficiency and numerical accuracy, a mesh density of 92 cells

162

per 2 mm was selected for all subsequent analysis in this study. A high level of grid quality

163

throughout the domain was achieved where the aspect ratio of the cells was below 5.2 and the

164

skewness was 0.51 37.


Page 9 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9

Isometric view

Lateral view

Top view, capillary inner diameter (ID) and outer diameter (OD)

165

Figure 1. Schematic representation of rectangular structure with two cylindrical capillary tubes

166

for the 3D computational simulations.

167

Time-step sizes below 1.5×10-6 s were selected to satisfy a courant number (CN= |V|∆t/∆X)

168

below 0.6. This dimensionless number quantifies the velocity magnitude at a specific cell and the

169

ratio of the time-step size to the element size. To compare the computational simulations with the

170

experimental results, a simulation time of 0.025 s was chosen.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 39

10

a

b

171

Figure 2. (a) Cross section view of the region prior to mesh refinement and after the first and

172

second mesh refinement. (b) Boundary conditions used in the computational domain.

173

2.2 Computational Parameters in ANSYS-Fluent

174

Two air bubbles, each with 2 mm diameter, were patched into the domain and placed at the

175

tips of the capillary tubes. The other ends of the capillary tubes were closed to the atmosphere.


Page 11 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11

176

The immersion depth of the capillary tubes in water was 20 mm and the air phase was patched

177

into the capillary tubes. A contact angle (θeq) of 160° between the capillary tube and the water-

178

air interface was defined. It was measured inside the air phase and calculated from the

179

experimental images by using an image processing technique (ImageJ). The rectangular structure

180

was open to the atmosphere by a pressure outlet set at the top of the rectangular structure as

181

shown in Figure 2b. Both water and air phases were assumed to be Newtonian fluids and

182

incompressible flow, which, by definition, assumed the speed of the flow was significantly lower

183

than the speed of sound

184

that these assumptions are valid for simulating the motion of coalesced bubbles

185

flow was considered laminar and unsteady. The physical properties for the water and air phases

186

specified in the computational simulations are shown in Table 1.

187

Table 1. Physical properties of water and air phases used in the computational simulations.

38

. Previous computational studies on bubble coalescence have shown

Physical property

Water phase

Fluid density (ρ), kg m-3

998.2

Dynamic viscosity (µ), mPa s

1.002 39

Surface tension (σ), mN m-1

72.8 11

Temperature (T), °C

20 11

39

18,27,29,30

. The

Air phase 1.204 39 1.805×10-2 39

20 11

188 189

2.2.1 Solution methods

190

A pressure-base solver and a time dependent solution were chosen as a solver approach for the

191

model. The pressure implicit with splitting operator (PISO) was applied as the pressure velocity



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 39

12

192

coupling in eq 2. A geo-reconstruction scheme was selected for tracking the water-air interface

193

as it is the most precise reconstruction scheme available in ANSYS-Fluent

194

based on a piecewise linear approximation of the water-air interface at a specific cell by a plane

195

in three dimensions 41. All the parameters used in ANSYS-Fluent are shown in Table 2.

196

Table 2. Parameters used in ANSYS-Fluent.

Parameter

Setting

Solver

Pressure-base

40

. This scheme is

Transient Solution method

Pressure velocity coupling: PISO

Discretization

Momentum: second order upwind Volume fraction: geo-reconstruction

197 198

2.2.2 VOF computational cases

199

The liquid film thickness between two bubbles at the rupture could not be measured in the

200

experimental case. Thus, two computational cases with initial separation distances (h) between

201

the two bubbles of 0.05 mm and 0.02 mm were generated. These separation distances were in the

202

range of 0.1 mm, which is the film thickness of the initial liquid drainage estimated by

203

Kirkpatrick and Lockett 42, and 10-5 mm which is the critical film thickness in water 43. It was not

204

possible to simulate smaller separation distances, as the cell size was 0.02 mm.

205

2.3 Experimental data

206

The experimental work of Bournival et al. 11 was used for validating the computational results

207

of the present study. These authors studied the coalescence dynamics of capillary bubbles in non-

208

ionic surfactant solutions by evaluating the oscillation of the resultant bubble. They generated


Page 13 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13

209

two bubbles attached to the tips of a pair of fine, vertical capillaries. The two air bubbles were

210

brought together in a controlled environment to allow coalescence using an electronic linear

211

actuator (T-LA28A, Zaber Technologies Inc.). A resolution time of 3.33 ms, which was the time

212

taken for the two bubbles to coalesce, was set for all experiments to evaluate cleanliness of the

213

system. The bubbles were left ageing for 90 s in Milli-Q water (resistivity of 18.2 MΩ cm)

214

before being brought in contact. The system was deemed cleaned if the bubbles coalesced within

215

one frame interval (i.e. between two frames). The details of the cleaning procedure may be found

216

in Bournival et al.

217

was maintained during all experiments. Five independent trials were run for replicating the same

218

experimental conditions. The authors used a high speed camera (Phantom 5, Vision Research

219

Inc., USA) to study the coalescence dynamics of the bubbles with a capture rate of 6024 frames

220

per second (fps). The analysis of the oscillation was performed as described in section 3.1.2.

221

3. RESULTS AND DISCUSSION

11

. All the experiments were conducted at a temperature of 20±2 °C, which

222

3.1 Analysis of bubble coalescence dynamics

223

Two capillary-held bubbles were simulated using the VOF method in ANSYS-Fluent to 11

224

evaluate the bubble coalescence dynamics, similar to Bournival et al.

225

computational study were firstly used to evaluate the effect of the initial separation distance

226

between the two bubbles on the neck expansion by calculating and comparing the Weber (We),

227

Reynold (Re) and Ohnesorge (Oh) numbers for both computational cases and the experimental

228

trial. The post-rupture oscillation behavior of coalesced bubbles was assessed by measuring the

229

projected surface area of the coalesced bubbles and comparing it to the experimental case.


. The results of the


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 39

14

230

3.1.1 Neck expansion during coalescence

231

The neck expansion during coalescence can be characterized by the Weber (We), Reynold (Re)

232

and Ohnesorge (Oh) numbers. The We number quantifies the relative magnitude of the inertial

233

force with the surface tension, as follows: Inertial force ρwater Vr 2 λ We Surface tension force σ

8

234

The Weber number takes into account the velocity of the neck radius (Vr), the density of the

235

water phase (ρwater), the horizontal width of the two coalesced bubbles (λ), which is determined

236

by the neck radius (r), and the surface tension (σ).

237 238

The Reynold number (Re) describes the ratio of inertial force to viscous force and it is given by the following equation: Re

Inertial force ρwater V+ λ Viscous force μ,-./+

9

239

where µwater is the dynamic viscosity of the water phase. The importance of the surface tension,

240

viscous and inertial forces on the motion of the neck expansion can be measured by Ohnesorge

241

number (Oh), defined as: Oh

μ,-./+ √We Re 3ρwater σλ

10

242

Frames of the VOF computational cases and experimental trial were produced and then

243

analyzed using ImageJ. The neck expansion of the two coalesced bubbles is illustrated in Figure

244

3a. The neck grows rapidly in both vertical directions until a maximum neck radius is achieved,

245

at 2.0 ms. Figure 3b illustrates the contour of the bubbles prior to coalescence compared with the

246

bubble shape at t = 0.8 ms. The initial contact point between the two bubbles was established as

247

the origin of the Cartesian coordinates system. The neck radius (r) was measured from the initial


Page 15 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


15

248

contact in a downward vertical direction to eliminate the neck growth restriction due to the

249

capillaries, as shown in Figure 3b. The initial shape of the two bubbles at t = 0 ms was used to

250

measure the horizontal width of the coalesced bubbles (λ). a

b

251

Figure 3. (a) the neck expansion contours as a function of time and (b) comparison of the shape

252

of the two bubbles (black continuous line) prior to coalescence with the bubble shape at t = 0.8

253

ms (blue dashed line), for the experimental trial conducted by Bournival et al. 11.

254

The motion of the neck radius (r) for two liquid drops in an inviscid system has been described

255

as a power-law dependence on time, r ∝ tA, where A is the power exponent

256

bubble coalescence have also reported that the growth of the neck radius follows a power-law

257

dependence on time during the early stage of coalescence 2,12,47.

44-46

. Studies on

258

In this study, the neck radius (r) was normalized by the average radius of the two initial

259

bubbles (Rb,ave=(Rb1+Rb2)/2). The time (t) was normalized by the capillary-inertial time

260

(τ=3ρwater Rb,ave3 / σ) which compared the surface tension force per unit area (σ/Rb,ave) with the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 39

16

261

inertial force per unit area 47. A prefactor, C, was calculated through the best fit of the data to the

262

power-law scaling: t A = C Rb,ave τ r

11

263

The growth of the neck for the experimental and simulated data was fitted to eq 11 and is

264

plotted in Figure 4. All three cases followed the power-law scaling with exponents in a range

265

between 0.40 and 0.55 for r/Rb,ave < 0.85. The observed initial neck radius in all three cases

266

deviated slightly, leading to different power-law exponents. However, these exponents were

267

relatively close to 0.5, which is the theoretical power exponent suggested by Eggers et al.

268

bubbles in water. Duchemin et al. 46 pointed out that the growth of the neck radius resulting from

269

the coalescence of two identical drops in an inviscid flow had a power-law exponent of 0.5. The

270

experimental and computational results were consistent with the experiments of Thoroddsen et

271

al. 12 who also showed an overall agreement with a power-law exponent of 0.5 for the expansion

272

of the neck radius of two bubbles in ethyl alcohol. The prefactors C in the present study were

273

slightly smaller, i.e., 1.10 – 1.25, than the prefactor observed by Thoroddsen et al. 12 who found a

274

value of 1.39 ±0.03. Therefore, the results presented in Figure 4 not only represented well the

275

motion of the neck radius as a power-law dependence but were also in good agreement with the

276

other theoretical and experimental findings.


44

for

Page 17 of 39

17

1.0 0.8

r/Rb,ave

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.6 0.4 0.2 0.0 0.0

0.1

277

0.2

t/τ

0.3

0.4

0.5

278

Figure 4. Neck radius (r) normalized by average radius of the two initial bubbles (Rb,ave) as a

279

function of t/τ for (solid red triangle) experimental trial with A=0.55 and C=1.10, (solid green

280

square) VOF computational case with h = 0.02 mm, A=0.40 and C=1.23, and (solid blue circle)

281

VOF computational case with h = 0.05 mm, A= 0.46 and C=1.25. The solid lines show the best

282

fit to eq 11. The initial radii of the bubbles for the experimental trial were Rb1 = 0.971 mm and

283

Rb2 = 0.982 mm. Both computational cases had two bubbles of equal radii of 1 mm.

284

Figure 5 illustrates the horizontal width of the two coalesced bubbles as a function of the

285

normalized neck radius for all three cases. Although the initial horizontal width observed in the

286

computational case with an initial separation distance of 0.05 mm was smaller than that for the

287

computational case with a separation distance of 0.02 mm, the two cases had nearly the same

288

horizontal width-neck radius curves. Based on this result, the computational case with h = 0.05

289

mm was selected for further analysis.



18

1.2 1.0 0.8

λ (mm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 39

0.6 0.4 0.2 0.0 0.0

290

0.2

0.4

0.6

0.8

1.0

r/Rb,ave

291

Figure 5. The horizontal width of the two coalesced bubbles (λ) as function of r/Rb,ave for (solid

292

red triangle) experimental trial, (solid green square) VOF computational case with h = 0.02 mm

293

and (solid blue circle) VOF computational case with h = 0.05 mm.

294

Snapshots of the experimental trial and contours of air volume fraction for the VOF case with

295

h = 0.05 mm are presented in Figure 6. Three representative r/ Rb,ave values were selected for

296

comparison. As the neck radius grew, the horizontal width between the coalesced bubbles

297

increased and therefore more surrounding liquid was pushed outwards. This motion created a

298

surface wave across the newly formed bubble surface. The VOF case captured the bubble surface

299

deformation very well as illustrated in Figure 6. However, the experiments had wider horizontal

300

widths than the computational case as shown graphically in Figure 5. The deviation could be

301

expected due to the piecewise linear scheme used for calculating the interface, which may have

302

compromised the reconstruction of the sharp edges during the neck expansion 48.

303 304


Page 19 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


19

a

b r/Rb,ave = 0.39

r/Rb,ave = 0.65

r/Rb,ave = 0.76

305

Figure 6. Comparison of the neck radius for (a) experimental trial conducted by Bournival et al.

306

11

307

the bubbles before coalescence while the green arrows indicate the horizontal width between the

308

two coalesced bubbles. The capillary tubes are 1.07 mm in diameter and act as a scale.

and (b) computational case with h = 0.05 mm. The dashed red lines represent the contour of

309

The neck growth velocity during the early stages of bubble coalescence was calculated from

310

the derivative of the neck expansion equation, Vr = dr/dt, and used to calculate the Weber and

311

Reynold numbers (eqs 8 and 9). The time at which the rupture of thin liquid film occurred was

312

chosen as time zero. For both computational cases and the experimental trial, the Weber number

313

was compared with the Reynold number (Figure 7). It is important to note that the experiments

314

had a time resolution of 0.2 ms from the rupture of the thin liquid film. For the experimental

315

data, the motion of the neck in the early stages of expansion was not fully counter-balanced by



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 39

20

316

the inertial force, although the value was close to unity as stressed by the Reynold number

317

(Re≤0.24). This finding is corroborated by the results of Stover et al. 2 who found that the motion

318

of the resultant bubbles is the result of the reduction in total surface area, which causes a

319

decrease in the surface energy of the system. As such the motion of the interface is driven by the

320

surface tension, which determines the amount of surface energy to be released. The excess

321

energy is passed to the surrounding liquid where it is opposed by the viscosity and the inertia of

322

the liquid.

323

On the other hand, the Weber number for the computational cases differed slightly from the

324

experimental results, showing the inertial force as superior in magnitude to the driving force (i.e.

325

surface tension) followed by surface tension for longer times although all results are close to

326

unity. Since the motion is driven by the surface tension, the small discrepancy could be attributed

327

to the poor capability of VOF to simulate the initial complex phenomenon of the liquid film

328

rupture. Previous numerical studies have pointed out the limitation of VOF to accurately

329

calculate the local curvature near the interface, which is then used for calculating the surface

330

tension force per unit volume 49,50. Moreover, water, in the absence of any surfactant, could have

331

small variations in surface tension

332

tension in the computational cases could also lead to small deviation in the results specifically at

333

the early stage of the neck expansion.

51

, which means that the assumption of a constant surface

334

Figure 8 shows that the Ohnesorge number was between 15 and 3. It indicates that the neck

335

motion driven by surface tension was initially restrained by a greater magnitude of the fluid

336

viscosity than the inertial force, followed by an exponential decay of viscous force to inertial

337

force ratio. Therefore, it appears that the simulation of the coalescence dynamics by VOF is

338

relatively accurate since the Weber numbers are relatively close to unity and the Ohnesorge


Page 21 of 39

21

numbers are consistent. However, the early stage of coalescence (first 1 ms) represented by the

340

initial expansion of the neck may lack in computational accuracy.

We

339

1.6

1.6

1.4

1.4

1.2

1.2

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0.0 0.0

0.5

1.0

1.5

2.0

Re

2.5

t (ms)

341 342

Figure 7. Comparison of (closed markers) Weber number and (open markers) Reynold number

343

as function of time for (red triangle) experimental trial, (green square) VOF computational case

344

with h = 0.02 mm and (blue circle) VOF computational case with h = 0.05 mm.

16.0 14.0 12.0 10.0 Oh

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8.0 Viscous force

6.0 4.0 2.0 0.0 0.0 345

0.5

1.0

1.5 t (ms)


2.0

2.5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 39

22

346

Figure 8. Ohnesorge number as function of time for (solid red triangle) experimental trial, (solid

347

green square) VOF computational case with h = 0.02 mm and (solid blue circle) VOF

348

computational case with h = 0.05 mm.

349 350

3.1.2 Oscillation behavior of coalesced bubbles

351

After obtaining the VOF results, the oscillating behavior was quantified following the 52

352

approach which has its basis in the equation proposed by Schulze

353

approximated the surface oscillation of a coalescing bubble using the following equation: A' A0 ' × e-δt sinω0 tφ B

. Bournival et al.

11

12

354

where A' and A0' denote the normalized relative projected area and initial amplitude,

355

respectively, δ is the damping constant in ms-1, φ is phase shift, t is time in ms, ω0 is the angular

356

frequency in ms-1 and B is an integration constant.

357

It is important to state that the normalized relative projected area A' in eq 12 (which was

358

calculated by the instantaneous projected area minus the average projected area towards an

359

infinite time divided by the average projected area towards an infinite time) is characterized by

360

five parameters, with three highly significant parameters. The initial amplitude describes the

361

magnitude of the deviation of the initial projected area of the bubble from the final equilibrium

362

state. The damping constant, which depends on the fluid and interfacial properties, shows if the

363

oscillation is slowly (low damping constant) or rapidly dampening (high damping constant). The

364

angular frequency gives an indication of how fast the oscillation completes one period and is

365

important for regulating the violent nature of the process at extreme points (peaks and troughs)

366

of the oscillation.


Page 23 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


23

367

The phase shift describes only the starting position of the bubble’s oscillation. It is not a very

368

critical parameter in this study and therefore not meaningful in comparing the simulations and

369

the experimental data. The integration constant is generally very close to zero and therefore not

370

relevant in the present analysis.

371

The calculation assumed equilibrium surface tension at the start of the simulation, which did

372

not change after coalescence. As such, no surface tension gradient caused by a change in surface

373

area was generated, limiting any Gibbs elasticity effect although relaxation may in fact occur 51.

374

High-frame-rate animations of the VOF numerical cases were produced and then analyzed

375

using ImageJ to track the changes in the projected surface area over time. The five parameters in

376

eq 12 were fitted using the solver command in Microsoft Excel and minimizing the sum of the

377

square error. An oscillation period of 1.5 was selected for fitting the simulation data to eq 12. A

378

previous study determined that 1.5 – 2 periods are sufficient to detect small deviations for the

379

damping constant and angular frequency 53.

380

The projected area curves for both computational case and experimental trials are illustrated in

381

Figure 9. In the VOF case, the average projected surface area towards an infinite time was

382

calculated by assuming that the total volume of the bubbles was preserved during the post-

383

rupture oscillation. From these two curves, it can be seen that the computational case followed

384

the oscillatory motion of the bubble after coalescence. The VOF results showed a higher initial

385

amplitude peak of the bubble oscillation compared to the experimental trial and indicated that the

386

coalesced bubbles in the VOF case were initially further from the final equilibrium state at which

387

a stable shape was reached. In the experiments, the damping surface oscillation was rapid

388

suggesting that the coalesced bubbles expanded and contracted quickly. The VOF computational

389

case showed a lower resistance of the bubble surface on the neck expansion after the film



24

390

rupture, leading to a lower damping constant (error of 13.3 %) and slower angular frequency

391

(error of 18.6 %) as shown in Table 3. 0.4

Normalized relative projected area

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 39

0.3 0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 0

5

10

15

20

25

t (ms)

392 393

Figure 9. Comparison of modeled data (eq 12) for (solid red triangle) experimental data and

394

(solid blue circle) VOF computational case with h = 0.05 mm.

395

Table 3. Calculated parameters from eq 12 for experimental and VOF case with an initial

396

separation distance of 0.05 mm.

Parameter

Experimentsa

VOF case

Initial amplitude, A0'

0.259

0.346

Damping constant, δ (ms-1)

0.060

0.052

Angular frequency, ω0 (ms-1)

0.586

0.477

Phase shift, φ

0.938

0.798


Page 25 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25

Integration constant, B 397

a

1.5×10-3

5.8×10-3

Obtained from averaging 5 trials in Bournival et al. 11

398 399

3.2 Flow dynamic analysis

400

Flow dynamic analysis provides an insight into the coalescence dynamic experienced by two

401

bubbles during the first cycle of post-rupture oscillation. Analysis of velocity vectors and

402

dynamic pressure gradients are presented below.

403

3.2.1 Velocity field analysis

404

Simulation snapshots of velocity vectors predicted at a vertical plane along the central axis of

405

the capillary tubes are illustrated in Figure 10a. The motion of the interface was captured in the

406

change of the projected surface area over time. Nine points in time during the initial stage of

407

bubble coalescence were chosen. As can be seen in Figure 10a, the velocity vectors indicate the

408

direction and the velocity magnitude of the airflow inside the coalesced bubbles and the

409

surrounding liquid. After the rupture of the liquid film separating the two bubbles, the air phase

410

flowed towards the bubbles’ initial contact point and then it moved in the vertical opposite

411

direction, as shown during the first 1.55 ms

412

producing the neck expansion. Small eddies with velocity magnitudes between 0.25 m/s and 0.75

413

m/s were observed at the top and bottom part of the neck as the curvature of the bubble surface

414

decreased. The neck expansion generated a surface wave that propagated horizontally in the

415

opposite direction, creating bulges at the edges of the bubble, as illustrated at 3.16 ms. The air

416

phase in the bulging area moved at much higher velocity, in the range of 1.75 m/s to 2.5 m/s,

417

than the rest of the bubble surface. The bubble contracted vertically as the surface wave recoiled

418

and the air phase flowed perpendicular towards the top of the neck creating eddies with small

54

. This airflow pushed the surrounding liquid,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 39

26

419

velocity magnitude at the bottom of the neck and leaving a stagnant void in the center of the

420

bubbles. It can be seen from 4.76 ms to 5.57 ms that the air phase near the bubble interface–

421

capillary tube flowed towards the top of the neck as the bubble continued to vertically contract. a

b

t = 0.14 ms

t = 0.74 ms

t = 1.55 ms

t = 2.35 ms


Page 27 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27

t = 2.96 ms

t = 3.16 ms

t = 3.96 ms



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 39

28

t = 4.76 ms

t = 5.57 ms

422

Figure 10. (a) Velocity vectors colored by velocity magnitude and (b) contour of dynamic

423

pressure gradients for VOF computational case with h = 0.05 mm. A logarithmic scale was

424

selected for the color map of the dynamic pressure contours. The capillary tubes are 1.07 mm

425

outer diameter and 0.69 mm inner diameter as indicated by the black contour lines. The

426

capillaries act as a scale.

427

3.2.2 Dynamic pressure gradient analysis

428

In addition to the velocity field analysis, the dynamic pressure gradients of fluid were also

429

assessed in ANSYS-Fluent to quantify the pressure associated with the difference of velocities at

430

the bubble interface. Dynamic pressure (pd=0.5×ρ|V|2), which represents the difference between

431

the total and static pressure in the system, is proportional to the density and square of the velocity


Page 29 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


29

432

vector magnitude of the fluid. Snapshots of dynamic pressure contours along a cross section

433

plane during the first 5.57 ms are shown in Figure 10b. The dynamic pressure scale was set

434

between 50 and 800 Pa for better comparison, and the white color in the snapshots represents

435

values in the range of 0 to 50 Pa. It is evident from the snapshot contours that the neck expansion

436

created dynamic pressure gradients at the interface where the motion was constrained by the

437

surrounding liquid. There were contours with higher dynamic pressure gradients in the bubble

438

interface, where the air phase flowed with high velocity, creating eddies as shown in Figure 10a.

439

Consequently, vertical bubble contraction was rapidly driven by the high-pressure gradients

440

around the bulged area that pushed the air phase inward, as shown at 3.16 ms. The high dynamic

441

pressure could be useful in understanding the detachment of particles when two bubbles coated

442

with hydrophobic particles coalesce. Bournival et al. 35 studied the coalescence dynamics of two

443

bubbles coated with silica particles in water.

444

Some snapshots of the experimental study were compared against the VOF computational case

445

to obtain an insight into the velocity fields and dynamic pressure gradients. The comparison was

446

based on the oscillation time after coalescence. The selected time of the surface contours for the

447

two uncoated bubbles presented in Figure 11 differed slightly from the coated bubbles. Apart

448

from the numerical error observed in the VOF results, one reason is the particle layer around the

449

loaded region, which may restrict the bubble surface deformation

450

detachment of a large number of particles from the bubble surface as the bulging area contracted.

451

At 3.36 ms, the bulged area moved at much higher velocity magnitude, in the range of 1.75 m/s

452

to 2.5 m/s, than the rest of the bubble surface. Therefore, a dynamic pressure greater than 800 Pa

453

is observed at the bulging area. The vertical contraction also created a difference in dynamic

454

pressure at the interface, which may explain the results by Bournival et al.


33

. Figure 11 illustrates the

36

, who pointed out


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 39

30

455

that the difference in pressure force between the air and the surrounding liquid exerts a force on

456

the particles, leading to particle detachment. Therefore, it could be argued that the pressure

457

gradients generated during the post-rupture oscillation of the coalesced bubble affect the stability

458

of particles attached at the interface. However, other factors such as particle size, density and

459

hydrophobicity may also help quantify the detachment of particles. t = 2.67 ms

t = 2.96 ms

t = 3.0 ms

t = 3.36 ms

Dynamic pressure > 800 Pa

t = 3.17 ms

t = 3.56 ms

t = 4.0 ms

t = 4.16 ms


Page 31 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


31

460

Figure 11. Velocity vectors and contours of dynamic pressure gradients for coalescence

461

dynamics of two uncoated bubbles, using the VOF computational case with h = 0.05 mm,

462

compared with an experimental study of bubble coalescence dynamics of two bubbles coated

463

with hydrophobized silica particles carried out by Bournival et al. 35. The capillary tubes may be

464

used to scale the images where capillary tubes of 1.04 mm in diameter were used in the

465

experiments and 0.69 mm inner diameter as indicated by the black contour lines on the

466

simulation.

467

4. CONCLUSIONS

468

The flow field arising from the merging of two equal-size bubbles was numerically

469

investigated for two separation distances using a finite volume methodology and was validated

470

against experimental data. The use of the VOF method allowed the detailed investigation of the

471

flow field and pressure gradients in the interface regions of the coalescing bubble. The results

472

indicated that there was a small deviation of the neck expansion motion showing a faster neck

473

growth at the smaller separation distance.

474

The dynamics of the coalesced bubbles were studied by evaluating the velocity fields and the

475

dynamic pressure gradients. The velocity vector showed that as the neck grew the surrounding



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 39

32

476

liquid was pushed outwards. The chaotic motion created a surface wave across the newly formed

477

bubble surface, which propagated in the opposite direction towards the edges of the bubble. The

478

neck growth created high dynamic pressure gradients at the top and bottom of the bubble

479

interface, where the air phase flowed with higher velocity than the surrounding liquid, which

480

created eddies. Overall, it was found that the computational case followed the oscillatory motion

481

very well, as observed in experiments, but with lower damping constant and lower angular

482

frequency. Further studies of the velocity fields and pressure gradients during the merging of two

483

bubbles with different types of surfactants are needed. These studies should be complemented by

484

investigation of the changes in surface energy during post-rupture oscillation to further elucidate

485

the oscillation mechanism. The distribution of surfactants along the interface during bubble

486

coalescence should also be studied numerically. The variation in the distribution of surfactant

487

may be resolved by using a non-constant surface tension force, which would take into account

488

surface elasticity. A similar approach may be needed to reconcile the initial growth of the neck,

489

which computer simulations considered to be inertia driven as opposed to surface tension driven.

490

Overall, the VOF method successfully tracked the interface motion of the water-air phases in the

491

current system.

492

SUPPORTING INFORMATION

493

Supplementary video of the VOF simulation with an initial separation distance of 0.05 mm can

494

be found in the online version.

495

ACKNOWLEDGEMENT

496

The authors wish to thank the following people: Mr Collin Turner of McGill University,

497

Canada for his contribution to the settings of the preliminary VOF computational case; Dr


Page 33 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


33

498

Yuqing Feng of CSIRO Mineral Resources (Australia) for helpful discussion on the model

499

settings and validation; and Dr Francisco Trujillo of the School of Chemical Engineering,

500

University of New South Wales (Australia) for facilitating access to a computer for running

501

preliminary models and vital discussion on the model settings. Mrs Yesenia Saavedra Moreno

502

would also like to thank the Australian government and University of New South Wales

503

(Australia) for financial support through the Australian Government Research Training Program

504

(RTP) scholarship.

505

ABBREVIATIONS 3D

three-dimensional

A

power exponent

A0'

initial amplitude

A'

normalized projected area

B

integration constant

C

prefactor

CFD

computational fluid dynamics

CN

courant number

Db,ave

average diameter of the two initial bubbles, mm

F

function of volume fraction

Fair

volume fraction of air phase

Fσv

surface tension force per unit volume, N m-3

g

gravity, m s-2

h

initial separation distance between the two bubbles, mm



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 39

34

n

normal vector to the interface

nwall

vector normal to the capillary tube

Oh

Ohnesorge number

p

pressure, Pa

pd

dynamic pressure, Pa

PISO

pressure implicit with splitting operator

r

neck radius, mm

Rb,ave

average radius of the two bubbles, mm

Rb1

radius of bubble 1, mm

Rb2

radius of bubble 2, mm

Re

Reynold number

t

time, ms

twall

vector tangent to the capillary tube

T

matrix transpose operator

T

temperature, °C

V

velocity vector

|V|

velocity vector magnitude, m s-1

Vr

velocity of the neck radius, m s-1

VOF

volume of fluid

We

Weber number

Greek letters ∆t

time-step size, s

∆X

element size, m


Page 35 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


35

δ

damping constant, ms-1

θeq

equilibrium contact angle, °

k

surface curvature of the interface

λ

horizontal width of the two coalesced bubbles, mm

µ

dynamic viscosity, mPa s

µair

dynamic viscosity of air phase, mPa s

µwater

dynamic viscosity of water phase, mPa s

ρ

density, kg m-3

ρair

density of air phase, kg m-3

ρwater

density of water phase, kg m-3

σ

surface tension, mN m-1

τ

capillary-inertial time, ms

φ

phase shift

ω0

angular frequency, ms-1

∇

gradient operator

506 507 508 509 510 511 512 513 514 515 516

REFERENCES (1) Marrucci, G. A theory of coalescence. Chemical Engineering Science 1969, 24, (6), 975-985. (2) Stover, R. L.; Tobias, C. W.; Denn, M. M. Bubble coalescence dynamics. AIChE journal 1997, 43, (10), 2385-2392. (3) Lessard, R. R.; Zieminski, S. A. Bubble coalescence and gas transfer in aqueous electrolytic solutions. Industrial & Engineering Chemistry Fundamentals 1971, 10, (2), 260269. (4) Craig, V. S.; Ninham, B. W.; Pashley, R. M. The effect of electrolytes on bubble coalescence in water. The Journal of Physical Chemistry 1993, 97, (39), 10192-10197.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 39

36

517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560

(5) Pashley, R.; Craig, V. Effects of electrolytes on bubble coalescence. Langmuir 1997, 13, (17), 4772-4774. (6) Christenson, H.; Bowen, R.; Carlton, J.; Denne, J.; Lu, Y. Electrolytes that show a transition to bubble coalescence inhibition at high concentrations. The Journal of Physical Chemistry C 2008, 112, (3), 794-796. (7) Yaminsky, V. V.; Ohnishi, S.; Vogler, E. A.; Horn, R. G. Stability of aqueous films between bubbles. Part 1. The effect of speed on bubble coalescence in purified water and simple electrolyte solutions. Langmuir 2010, 26, (11), 8061-8074. (8) Prud'home, R. K.; Warr, G. G. Foams in mineral flotation and separation process. In Foams: Theory, Measurements, and Applications, Prud'home, R. K.; Khan, S. A., Eds. Marcel Dekker New York, 1995; Vol. 57, pp 511-553. (9) Penfold, J.; Staples, E.; Tucker, I.; Thompson, L.; Thomas, R. Adsorption of nonionic mixtures at the air–water interface: effects of temperature and electrolyte. Journal of Colloid and Interface Science 2002, 247, (2), 404-411. (10) Wang, L.; Yoon, R.-H. Role of hydrophobic force in the thinning of foam films containing a nonionic surfactant. Colloids and Surfaces A: Physicochemical and Engineering Aspects 2006, 282, 84-91. (11) Bournival, G.; Ata, S.; Karakashev, S. I.; Jameson, G. J. An investigation of bubble coalescence and post-rupture oscillation in non-ionic surfactant solutions using high-speed cinematography. Journal of Colloid and Interface Science 2014, 414, 50-58. (12) Thoroddsen, S.; Etoh, T.; Takehara, K.; Ootsuka, N. On the coalescence speed of bubbles. Physics of Fluids 2005, 17, (7), 071703-4. (13) Yeoh, G. H.; Tu, J. Introduction. In Computational Techniques for Multiphase Flows, Elsevier: Oxford, 2010; pp 1-20. (14) Tomiyama, A.; Zun, I.; Sou, A.; Sakaguchi, T. Numerical analysis of bubble motion with the VOF method. Nuclear Engineering and Design 1993, 141, (1), 69-82. (15) Chen, L.; Li, Y.; Manasseh, R. The coalescence of bubbles-a numerical study. In Third International Conference on Multiphase Flow: ICMF, Lyon, 1998; Vol. 98, pp 1-8. (16) van Sint Annaland, M.; Deen, N. G.; Kuipers, J. A. M. Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method. Chemical Engineering Science 2005, 60, (11), 2999-3011. (17) Bothe, D. VOF-simulation of fluid particle dynamics. In Proc. 11th Workshop on Twophase Flow Predictions, Merseburg, 2005; pp 1-13. (18) Hasan, N.; Zakaria, Z. B. Computational approach for a pair of bubble coalescence process. International Journal of Heat and Fluid Flow 2011, 32, (3), 755-761. (19) Deyranlou, A.; Passandideh-Fard, M. The study of bubble coalescence in coaxial and side-by-side motions. In 20th Annual International Conference on Mechanical EngineeringISME 2012, Shiraz, 2012; pp 1-4. (20) Keshavarzi, G.; Pawell, R. S.; Barber, T. J.; Yeoh, G. H. Transient analysis of a single rising bubble used for numerical validation for multiphase flow. Chemical Engineering Science 2014, 112, 25-34. (21) Deising, D.; Marschall, H.; Bothe, D. A unified single-field model framework for Volume-Of-Fluid simulations of interfacial species transfer applied to bubbly flows. Chemical Engineering Science 2016, 139, 173-195.


Page 37 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


37

561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605

(22) Feng, J.; Li, X.; Bao, Y.; Cai, Z.; Gao, Z. Coalescence and conjunction of two in-line bubbles at low Reynolds numbers. Chemical Engineering Science 2016, 141, 261-270. (23) Muradoglu, M.; Tryggvason, G. A front-tracking method for computation of interfacial flows with soluble surfactants. Journal of Computational Physics 2008, 227, (4), 2238-2262. (24) Hua, J.; Stene, J. F.; Lin, P. Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method. Journal of Computational Physics 2008, 227, (6), 3358-3382. (25) Sussman, M.; Smereka, P.; Osher, S. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics 1994, 114, (1), 146-159. (26) Croce, R.; Griebel, M.; Schweitzer, M. A. Numerical simulation of bubble and droplet deformation by a level set approach with surface tension in three dimensions. International Journal for Numerical Methods in Fluids 2010, 62, (9), 963-993. (27) Chakraborty, I.; Biswas, G.; Ghoshdastidar, P. A coupled level-set and volume-of-fluid method for the buoyant rise of gas bubbles in liquids. International Journal of Heat and Mass Transfer 2013, 58, (1), 240-259. (28) Hirt, C. W.; Nichols, B. D. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 1981, 39, (1), 201-225. (29) Chen, R. H.; Tian, W. X.; Su, G. H.; Qiu, S. Z.; Ishiwatari, Y.; Oka, Y. Numerical investigation on coalescence of bubble pairs rising in a stagnant liquid. Chemical Engineering Science 2011, 66, (21), 5055-5063. (30) Zhang, A.; Sun, P.; Ming, F. An SPH modeling of bubble rising and coalescing in three dimensions. Computer Methods in Applied Mechanics and Engineering 2015, 294, 189-209. (31) Brackbill, J.; Kothe, D. B.; Zemach, C. A continuum method for modeling surface tension. Journal of Computational Physics 1992, 100, (2), 335-354. (32) Clift, R.; Grace, J. R.; Weber, M. E. Wall effects. In Bubbles, Drops, and Particles, Academic Press: New York, 2005; pp 221-241. (33) Ata, S. Coalescence of bubbles covered by particles. Langmuir 2008, 24, (12), 60856091. (34) Bournival, G.; Ata, S. Packing of particles on the surface of bubbles. Minerals Engineering 2010, 23, (2), 111-116. (35) Bournival, G.; de Oliveira e Souza, L.; Ata, S.; Wanless, E. J. Effect of alcohol frothing agents on the coalescence of bubbles coated with hydrophobized silica particles. Chemical Engineering Science 2015, 131, 1-11. (36) Bournival, G.; Ata, S.; Wanless, E. J. The behavior of bubble interfaces stabilized by particles of different densities. Langmuir 2016, 32, (25), 6226-6238. (37) Yeoh, G. H.; Tu, J. Solution methods for multi-phase flows. In Computational Techniques for Multiphase Flows, Elsevier: Oxford 2010; pp 95-242. (38) Tu, J.; Yeoh, G. H.; Liu, C. Some advanced topics in CFD. In Computational Fluid Dynamics: A Practical Approach, Elsevier: Oxford, 2012; pp 349-394. (39) Haynes, W. M. CRC handbook of chemistry and physics. In 97th Edition ed.; CRC press 2017. (40) Ozkan, F.; Worner, M.; Wenka, A.; Soyhan, H. S. Critical evaluation of CFD codes for interfacial simulation of bubble-train flow in a narrow channel. International Journal for Numerical Methods in Fluids 2007, 55, (6), 537-564.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 39

38

606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637

(41) Youngs, D. L. Time-dependent multi-material flow with large fluid distortion. Numerical Methods for Fluid Dynamics 1982, 24, (2), 273-285. (42) Kirkpatrick, R.; Lockett, M. The influence of approach velocity on bubble coalescence. Chemical Engineering Science 1974, 29, (12), 2363-2373. (43) Kim, J. W.; Lee, W. K. Coalescence behavior of two bubbles in stagnant liquids. Journal of Chemical Engineering of Japan 1987, 20, (5), 448-453. (44) Eggers, J.; Lister, J. R.; Stone, H. A. Coalescence of liquid drops. Journal of Fluid Mechanics 1999, 401, 293-310. (45) Menchaca-Rocha, A.; Martínez-Dávalos, A.; Nunez, R.; Popinet, S.; Zaleski, S. Coalescence of liquid drops by surface tension. Physical Review E 2001, 63, (4), 046309-1. (46) Duchemin, L.; Eggers, J.; Josserand, C. Inviscid coalescence of drops. Journal of Fluid Mechanics 2003, 487, 167-178. (47) Egan, E.; Tobias, C. Measurement of interfacial re-equilibration during hydrogen bubble coalescence. Journal of the Electrochemical Society 1994, 141, (5), 1118-1126. (48) Rider, W. J.; Kothe, D. B. Reconstructing volume tracking. Journal of Computational Physics 1998, 141, (2), 112-152. (49) Sussman, M.; Puckett, E. G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. Journal of Computational Physics 2000, 162, (2), 301-337. (50) Gerlach, D.; Tomar, G.; Biswas, G.; Durst, F. Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows. International Journal of Heat and Mass Transfer 2006, 49, (3), 740-754. (51) Liu, M.; Beattie, J. K.; Gray-Weale, A. The surface relaxation of water. The Journal of Physical Chemistry B 2012, 116, (30), 8981-8988. (52) Schulze, H. J. Analysis of the elementary stages of the flotation process. In In Physicochemical Elementary Processes in Flotation : An Analysis from the Point of View of Colloid Science including Process Engineering Considerations, Elsevier: Amsterdam, 1984; pp 53234. (53) Bournival, G.; Pugh, R. J.; Ata, S. Examination of NaCl and MIBC as bubble coalescence inhibitor in relation to froth flotation. Minerals Engineering 2012, 25, (1), 47-53. (54) Ohnishi, M.; Azuma, H.; Straub, J. Study on secondary bubble creation induced by bubble coalescence. Advances in Space Research 1999, 24, (10), 1331-1336.

638


Page 39 of 39


39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Analysis of Bubble Coalescence Dynamics and Postrupture

Recommend Documents