Experimental and numerical study of the periodic bubbling regime in planar co-flowing air–water sheets

International Journal of Multiphase Flow 50 (2013) 106–119

Contents lists available at SciVerse ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

Experimental and numerical study of the periodic bubbling regime in planar co-flowing air–water sheets C. Gutiérrez-Montes a, R. Bolaños-Jiménez a, A. Sevilla b, C. Martínez-Bazán a,⇑ a b

Área de Mecánica de Fluidos, Departamento de Ingenierı´a, Mecánica y Minera, Universidad de Jaén, Campus de las Lagunillas, 23071 Jaén, Spain Área de Mecánica de Fluidos, Departamento de Ingenierı´a, Térmica y de Fluidos, Universidad Carlos III de Madrid, 28911 Leganés, Spain

a r t i c l e

i n f o

Article history: Received 21 August 2012 Received in revised form 29 October 2012 Accepted 11 November 2012 Available online 28 November 2012 Keywords: Bubble formation Bubbling frequency Planar co-flowing jets

a b s t r a c t We present an experimental and numerical study of the periodic breakup of an air sheet surrounded by two co-flowing water streams. For a fixed liquid-to-gas thickness ratio, h = hw/ha ’ 5.27, we explore the dependence of the bubbling regime on the two main governing parameters, namely the Weber number, We ¼ qw u2w ha =r, and the velocity ratio, K = uw/ua, where uw and ua are the mean water and air velocities at the exit slit, respectively, and ha and hw are the half-thicknesses of the air and water sheets at the outlet. The bubble formation frequencies obtained experimentally are in good quantitative agreement with the results of numerical computations performed with the volume-of-fluid technique. The analysis of the results obtained allows a detailed characterization of the dynamics of the bubble formation process in planar co-flows. In particular, using the pressure fluctuations in the gas stream and the time evolution of the bubble interface, both obtained from the numerical simulations, the bubbling process is described through a two-stage mechanism governed by distinct physical phenomena: the first stage is characterized by the formation of a neck that moves downstream at the liquid velocity, and it collapses towards the symmetry plane during the second stage. A simple scaling law is proposed for the bubbling frequency, that shows a good agreement with our experimental and numerical results.  2012 Elsevier Ltd. All rights reserved.

1. Introduction Bubble formation is a phenomenon which embraces numerous natural phenomena and an increasing number of industrial applications, such as in chemical, pharmaceutical or food industries, among many others. One of the most traditional ways to generate bubbles is by injecting a gas flow through a cylindrical needle which discharges in a still liquid. This configuration has been ~uz and Prosperetti, 1993; Longuet-Higgins widely studied (see Og et al., 1991; Davidson and Schuler, 1960; Kumar and Kuloor, 1976, among others), and it is particularly useful under conditions of low gas flow rate since a well-controlled periodic generation of bubbles can be achieved. Under these conditions, the minimum size of the bubble generated is limited by the Fritz Volume, given by a balance between surface tension and buoyancy forces. However, for flow rates larger than a critical value, the periodic behavior is no longer observed, bubble coalescence occurs and bubbles of different sizes are produced. Such limitations make this method inappropriate in many emerging engineering applications, where monodisperse micron-sized bubbles are necessary, such as ⇑ Corresponding author. E-mail addresses: [email protected] (C. Gutiérrez-Montes), [email protected] (R. Bolaños-Jiménez), [email protected] (A. Sevilla), [email protected] (C. Martínez-Bazán). 0301-9322/$ - see front matter  2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2012.11.003

those related to medical and material industries. Therefore, new techniques are demanded to decrease the bubble size and, consequently, a deeper understanding of the physical mechanisms controlling the bubble formation processes is required. A widely extended method to reduce the bubble size, avoid bubble coalescence and increase the bubbling frequency, is the so-called co-flow configuration, which consists of injecting the gas flow inside a liquid stream flowing in the same direction. Due to its importance, the cylindrical co-flow configuration has been intensively investigated and it is fairly well understood nowadays (Chuang and Goldschmidt, 1970; Sevilla et al., 2002; Sevilla et al., 2005a,b; Gordillo et al., 2007). However, the planar, co-flow configuration, where an air film discharges between two parallel liquid sheets, has not received as much attention as the cylindrical one, although it could represent an alternative way to produce bubbles of a well controlled size in a massive way. This type of geometry has been investigated in the opposite problem of atomization of droplets in air-blasted liquid sheets (Lozano and Barreras, 2001; Lozano et al., 2001; Lozano et al., 2005; Hauke et al., 2001; López-Pagés et al., 2004; Park et al., 2004). Regarding the break-up of air sheets in liquids, only a few theoretical studies can be found in the literature for the particular case of planar air sheets without a liquid co-flow (Li and Bhunia, 1996; Li and Bhunia, 1997; Lezzi and Prosperetti, 1991). Considering the co-flowing configuration, in their previous study, Bolaños-Jiménez

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

et al. (2011) identified two different flow regimes in planar water– air–water sheets: a bubbling regime, characterized by the periodic breakup of the air sheet into individual elongated bubbles close to the exit slit, and a jetting regime, where the air sheet does not break-up in the near field. They defined the transition from the bubbling to the jetting regime in terms of the Weber number, We, and the liquid-to-gas velocity ratio, K, for a given value of the initial liquid-to-gas thickness ratio, h ’ 5.27, showing that, unlike in the cylindrical configuration, surface tension has an stabilizing effect in the planar case. Once the jetting-bubbling transition in plane air–water sheets has been understood and characterized, the purpose of the present work is to analyze and describe the bubbling regime with the aim of modeling the breakup process and, thus, controlling the bubbling frequency and the bubble size. A complete experimental study of this problem is nearly impossible due to the limitations of the optical access, difficulties in performing an exhaustive parametric study and the disparity of length and time scales taking place in the bubbling process. Thus, detailed numerical simulations have been combined with experiments to tackle the problem. Different kinds of two-phase numerical techniques can be applied to track the interface, such as front-tracking (Unverdi and Tryggvason, 1992), Volume of Fluid (VoF) method (Hirt and Nichols, 1981; Scardovelli and Zaleski, 1999; Gueyffier et al., 1999), Level Set method (LS) (Sussman et al., 1994) or multiple hybrid approaches (Sussman and Puckett, 2000; Enright et al., 2002; Menard et al., 2007; Herrmann, 2008; Chakraborty et al., 2011), among others. In the present work, the VoF method has been used, a technique which has been successfully applied in the study of different problems, such as liquid atomization (Siamas and Jiang, 2007; Fuster et al., 2009), liquid sheets between co-flowing gas streams (López-Pagés et al., 2004), bubble generation (Pianet et al., 2010; Arias et al., 2011), drop impact onto a liquid layer (Berberovic´ et al., 2009), liquid bridges in a coaxial gas flow (Herrada et al., 2011) or Taylor-bubble generation (Zheng et al., 2007), among others.

(a)

107

Thus, in this paper we focus on characterizing the periodic generation of bubbles in an innovative planar co-flow configuration by means of experiments and numerical simulations, providing a suitable scaling law for the bubbling time as a function of the control parameters. The work is organized as follows. The experimental set-up, control parameters and experimental techniques are described in Section 2.1. Section 2.2 is devoted to the definition of the numerical method. In Section 3.1, the bubbling process is described in detail. The results corresponding to the bubble characteristics are reported in Section 3.2, a scaling law for the bubbling frequency is proposed in Section 3.3 and a description of different mechanisms leading to suction in the air stream is given in Section 3.4. Finally, Section 4 is devoted to conclusions. 2. Experimental and numerical techniques 2.1. Description of the experiments The co-flowing air–water sheets were created with the same experimental set-up used in Bolaños-Jiménez et al. (2011), displayed in Fig. 1. A detailed view of the nozzle employed to generate the sheets is sketched in Fig. 1b. It consists of a planar nozzle with an inner channel to inject the air film, and two outer channels located at both sides of the air slit to produce a pair of symmetric coflowing water sheets. The central gap was formed by two parallel brass plates separated 2Hi = 0.47 mm through which the air flow was injected. The walls of both plates were carefully machined at their sides in contact with the water stream, obtaining a sharp edge at its tip of thickness Ho  Hi = 0.22 mm. The water nozzle, made of Plexiglas, was designed to provide two parallel water streams of width Hw  Ho = 1.945 mm (Hw = 2.4 mm) with a uniform velocity profile at the exit. The length of the air–water nozzle was 150 mm, long enough to make sure that a fully developed, parabolic profile was provided in the air flow. Moreover, with the aim of obtaining a quasi-two-dimensional configuration, the spanwise length was b = 41.75 mm long, i.e. two orders of magnitude larger

(b)

Fig. 1. (a) Sketch of the experimental facility. (b) Detail of the nozzle used to generate the co-flowing sheets, together with the main geometrical parameters: Hi = 0.235 mm denotes half the gap between the central plates confining the air film, Ho = 0.455 mm is the half-distance between their outer sides at the tip, Hw = 2.4 mm denotes half the distance of the exterior walls at the outlet, and b = 41.75 mm is the spanwise length.

108

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

Fig. 2. Detail of the air and water sheets at the exit of the supplying nozzle, together with the dimensional parameters involved in the physical problem.

than the thickness of the air stream. Finally, all the pieces were joined by two lateral Plexiglas walls which sealed the device, letting us visualize the flow inside. As can be observed in Fig. 1a, both the air and water flows were supplied through a cylindrical reservoir placed over the nozzle. To obtain a constant water flow rate during each experiment, the water stream was supplied from a constant-pressure tank, and injected into the distribution vessel through a pair of symmetric inlets placed at its bottom. A piece of foam was set at the entrance to attenuate the turbulence generated by the discharging jet and, thus, smooth the flow. In addition, the reservoir had an inner wall to guide the water stream towards the inlet of the nozzle, as well as to divide the water flow equally between the two independent channels. The water flow rate was varied from Qw ’ 10 l/min to Qw ’ 30 l/min, indicating that the water velocity at the exit, defined as uw = Qw/[(Hw  Ho)2b], changed from uw ’ 1 to uw ’ 3 m/s. Regarding the air flow, it was supplied from a pressurized bottle, controlled through a precision valve and measured with a digital mass flow meter. The air was injected into the central channel of the nozzle through several holes distributed along the top of the nozzle slit. The gas flow rate, Qa, was varied from Qa ’ 10 l/min to Qa ’ 30 l/min, and thus, the mean air velocity at the exit, ua = Qa/(2Hib), changed from ua ’ 8 m/s to ua ’ 25 m/s. Fig. 2 shows a sketch of the physical problem indicating the dimensional parameters where ha and hw denote, respectively, the half-thicknesses of the air and water sheets at the exit of the nozzle, ua and uw are the mean velocities of both sheets at the outlet, r is the air–water surface tension coefficient, qa and qw are the air and water densities, la and lw their dynamic viscosities, and g represents the gravitational acceleration. The experiments and numerical simulations performed revealed that the inner air– water interface at the exit was always pinned at the outer wall of the splitter plate, thus giving ha = Ho. On the other hand, the outer water–air interface at the exit was pinned at the inner edge of the water nozzle, so that hw = Hw. Regarding the dimensionless parameters governing the problem, the air-to-water density ratio, S = qa/qw = 0.0012, the initial water-to-air thickness ratio, h = hw/ ha = Hw/Ho ’ 5.27, and the dimensionless wall thickness, b = Hi/ Ho ’ 0.52, were kept constant in the present work. In addition, since the range of the water Reynolds numbers was 2000 [ Rew = qwuw (Hw  Ho)/lw [ 6000, viscous effects in the water sheet were considered negligible. Similarly, the air Reynolds number, Rea = qauaHo/la, varied in the range 245 < Rea < 770, so that the flow in the air sheet was considered to be a Hagen-Poiseuille flow. In addition, the effect of gravity was also negligible in the bubble formation process, since the Froude number based on the breakup distance of the bubbles was always very large, Fr ¼ u2w =ðgli Þ  1, where li is a typical value of the breakup length. Consequently, the two main control parameters governing the

present configuration are the Weber number, We ¼ qw u2w Ho =r, and the water-to-air velocity ratio, K = uw/ua, which in the present work varied in the ranges 6 6 We 6 50 and 0.04 6 K 6 0.4. Experimental measurements were performed by recording high-speed movies of the sheets in the spanwise (x, z) plane by means of a Photron camera together with a Sigma 105 mm micro-lens. The videos were taken at frame rates between 10,000 s1 and 20,000 s1, with acquisition boxes going from 256  592 to 256  352 pixels, and spatial resolutions of around 100 lm/pixel. An halogen lamp was used as a light source, and a diffuser was placed to illuminate the measuring region with a backlighting technique (see Fig. 1a). 2.2. Numerical method Numerical simulations of the problem were conducted with the free, open source Computational Fluid Dynamics software package OpenFOAM.1 To perform the simulations, an unsteady, laminar and incompressible, two-dimensional flow was assumed, neglecting, therefore, the spanwise variations, and without imposing the symmetry condition. Under these assumptions, the solver interFOAM was used. This solver is specific for two incompressible, isothermal, immiscible fluids, implementing a Volume of Fluid (VoF) interfacecapturing approach to describe the time evolution of the air–water interface. Thus, a single set of continuity and momentum equations was solved for both phases together with a transport equation for an indicator function, i.e. the volume fraction, c. In this case, the governing equations take the form,

r  u ¼ 0;

ð1Þ

@ðquÞ þ r  ðquuÞ  r  ðlruÞ  ðruÞ  rl ¼ rp þ qg þ rjrc; @t ð2Þ @c þ r  ðucÞ þ r  ½ur cð1  cÞ ¼ 0; @t

ð3Þ

where u is the velocity vector and c the volume fraction of the liquid phase (water in our case). The surface tension forces at the interface were evaluated by means of the continuum-surface-force (CSF) model (Brackbill et al., 1992) as,

fr ¼ rjrc;

ð4Þ

where r is the water-to-air surface tension coefficient and j is the mean curvature of the free surface. Eq. (3) includes an additional term, r  [urc(1  c)], where ur is the liquid–gas relative velocity, which ‘‘compresses’’ the free surface, improving the interface 1

http://www.openfoam.com/.

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

resolution (see Rusche, 2002; Berberovic´ et al., 2009 for further details). In the VoF simulations, the physical properties at any location of the domain correspond to those of the fluid occupying the position, varying only at the interface according to,

q ¼ ql c þ qg ð1  cÞ; l ¼ ll c þ lg ð1  cÞ;

ð5Þ ð6Þ

where ql and qg, represent the liquid and gas densities, respectively, and ll and lg, the corresponding viscosities. Fig. 3a shows a sketch of the computational domain and its main dimensions. It includes the whole injection nozzle, to simulate the air and water flows inside the nozzle and thus determine the actual velocity profile at the exit, and part of the quiescent discharging atmosphere of length 426Hi and width 102Hi. The entire computational domain consisted of a structured, non-uniform mesh (see detail in Fig. 3b) with overall dimensions of 1064Hi (575Ho) long and 220Hi (119Ho) wide, containing nearly 3.8  105 cells. The domain was delimited by four different boundaries, namely the air inlet area, Xi,a, the water inlet area, Xi,w, the solid walls, Xs, and the outer discharging atmosphere, Xo, where the following boundary conditions were imposed,

n  rp ¼ 0; u ¼ ð0; ua ; 0Þ; c ¼ 0; for x 2 Xi;a ; t P 0;

ð7Þ

n  rp ¼ 0; u ¼ ð0; uw ; 0Þ; c ¼ 1; for x 2 Xi;w ; t P 0;

ð8Þ

n  rp ¼ n  rc ¼ 0; u ¼ ð0; 0; 0Þ; for x 2 Xs ; t P 0; p ¼ 0; n  ru ¼ n  rc ¼ 0; for x 2 Xo ; t P 0;

ð9Þ ð10Þ

In Eqs. (7)–(10), n represents the normal surface vector directed outwards from the computational domain. In particular, Eqs. (7) and (8) indicate that uniform velocity profiles were set for the air and the water streams at the entrance of the injection nozzle. The non-slip condition, given by Eq. (9), was imposed at the solid surfaces, and outflow conditions, Eq. (10), were specified at the far field of the discharging domain, Xo. To discretise the equations, the Gauss Linear, central differencing, discretisation scheme was used for the spatial terms (Lax and Wendroff, 1960; Jasak, 1996), except

109

for the convective terms, which were discretised using a van Leer total variation diminishing (TVD) scheme (Van Leer, 1974; Jasak, 1996). A Crank-Nicolson scheme (Crank and Nicolson, 1947) with an off-centering coefficient of 0.5 was used for the temporal terms and the PISO algorithm was used for the pressure–velocity coupling (Issa, 1986). Finally, the time step was controlled limiting the maximum Courant-Friedrichs-Lewy (CFL) number to 0.5, in order to ensure the stability of the solution procedure. A grid sensitivity study was performed to guarantee that the results did not depend on the grid density. In this sensitivity study, the entire air–water injection nozzle was not included and only the discharging domain and short inlet channels for the air and the water streams were considered. Thus, at the entrance channels, close to the discharging domain, x = 22Hi, a parabolic profile was imposed for the air stream and uniform flow was specified for the water flow. Table 1 shows the bubbling frequency obtained for We = 38.8 and K = 0.147. Note that, increasing the number of cells from 3.6  105 in mesh 4 to 3.8  106 in mesh 6, shows differences on the bubbling frequency lower than 1.5%, increasing the computational time by a factor of 87 and, thus, mesh 4 was adopted in our simulations. In addition, we studied the influence of the size of the discharging domain on the bubbling frequency performing two different tests: (1) the length of the discharging atmosphere was varied from 213Hi to 852Hi keeping its width equal to 102Hi and (2) the width of the discharging domain was varied from 102Hi to 153Hi maintaining its length equal to 426Hi. The results obtained indicated that the bubbling frequency did not vary with the size of the domain for domains larger than 102Hi  213Hi, thus mesh 4 in Table 1 seemed adequate to perform the simulations reported in this work.

3. Analysis of results In this section, we will present the experimental and numerical results on the bubbling process. Based on the gathered evidences, a

Fig. 3. (a) Scketch and main dimensions of the computational domain and (b) detail of a mesh near the exit of the injection nozzle.

110

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

Table 1 Bubbling frequency, fb, obtained from the numerical simulations using different grid densities, for We = 38.8 and K = 0.147. The CPU time indicates the computational hours per simulated second. The relative error is computed with respect to finest grid (sixth row). The computational domain has width 102Hi and length 426Hi. In this table, Mesh number 4, highlighted in bold, was used in our simulations. Mesh

Cells

CPU time (h)

fb (Hz)

Relative error

1 2 3 4 5 6

1.1e5 1.5e5 2.2e5 3.6e5 1.2e6 3.8e6

36.78 56.65 89.58 151.90 1812.67 13145.64

157 156 154 149 150 151

0.0397 0.0331 0.0199 0.0132 0.0066 –

physical explanation of the periodic bubbling process will be provided, allowing us to propose a simple scaling law for the bubble formation frequency. It is important to point out that the experiments were performed under constant air flow rate conditions, corresponding to the fact that the pressure drop along the air channel was always much larger than the typical pressure fluctuations inside the forming bubble (see Gordillo et al., 2007 for details). Accordingly, the simulations were performed at constant air flow rate conditions as well. 3.1. Description of the bubbling process Fig. 4 illustrates the two different flow regimes that exist depending on the water and the air velocities, as reported in detail by Bolaños-Jiménez et al. (2011). For a fixed value of the water velocity, a jetting regime is observed if the air velocity is sufficiently small (Fig. 4a and e). However, for air velocities larger than a certain critical value, that depends on the water velocity, the flow changes to a periodic bubbling regime (Fig. 4b–d and f–h), whose detailed characterization is the main purpose of the present work. To this end, it is useful to consider the side view of the inner air– water interface, i.e. a cut through the (x, y) plane shown in Figs. 2 and 4. Unfortunately, the experiments performed did not allow us to obtain clear views of the flow in this lateral plane. In fact, only the recordings taken in the spanwise direction, or (x, z) plane, like

those shown in Fig. 4a–d, were used to obtain all the experimental results presented in Section 3.2. Thus, to deduce the physical mechanisms underlying the periodic formation of bubbles, we have made use of numerical evidences. Note that this procedure is justified since, as corroborated below, the numerical simulations provide a good overall quantitative agreement with the experiments as far as the global flow features, like the bubbling frequency or the bubble geometry, are concerned. The main events associated to the formation of a bubble in the periodic bubbling regime may be briefly summarized as follows. After the breakup of a bubble, a lump of air of length li, hereafter denoted intact ligament, remains attached to the injection nozzle (see sketch in Fig. 2). The intact ligament is slightly thicker than Ho and, therefore, a small step is formed at its joint with the injection nozzle. This step will be the origin of a neck that will eventually collapse generating a new bubble. While the air is injected, both the neck and the tip of the ligament propagate downstream at a velocity close to the water velocity, during the formation process of the next bubble. As the neck moves downstream, it collapses towards the symmetry plane, following an accelerated closure process that is completed at a distance li from the outlet, thereby providing the initial conditions for the next bubble formation cycle. A detailed picture of the bubbling process is provided by Fig. 5, extracted from a reference numerical simulation performed for We = 38 and K = 0.147. Specifically, Fig. 5a shows the temporal evolution of the gauge pressure, in the air stream, at the center of the exit slit, pe(s)  p0, nondimensionalized with twice the air dynamic pressure, qa u2a , as a function of the dimensionless time, s = t uw/Ho, while panels Fig. 5b–h display seven snapshots of the dimensionless inner air–water interface at different instants of the periodic process, indicated by square symbols in Fig. 5a. It should be emphasized that pe  p0 actually represents the difference between the air and water pressures at the outlet, since the interface remains nearly parallel to the exit walls in the region close to the outlet, as can be seen in Fig. 5b–h. The analysis of Fig. 5 indicates that, similarly to the cylindrical case reported by Sevilla et al. (2005a), a bubbling event can be described by a

(a)

(b)

(c)

(e)

(f)

(g)

(d)

(h)

Fig. 4. Top row: sequence of experimental images in the (x, z) plane corresponding to We = 26.7 at different water-to-air velocity ratios: (a) K = 0.167, (b) K = 0.151, (c) K = 0.134, (d) K = 0.110. Bottom row: corresponding two-dimensional numerical simulations in the (x, y) plane for We = 25.6, and (e) K = 0.169, (f) K = 0.145, (g) K = 0.127 and (h) K = 0.112. The black regions in (e)–(h) represent the water stream.

111

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

  Fig. 5. (a) Dimensionless gauge pressure inside de air stream at the center of the exit slit, ðpe  p0 Þ= qa u2a , as a function of the dimensionless time, s = t uw/Ho. Panels b–h show the inner air–water interface in dimensionless form, for several time instants indicated with squares in the main plot. All the information corresponds to a numerical simulation performed for We = 38 and K = 0.147.

two-stage process: a neck formation stage, followed by a neck collapse stage. Nevertheless, the case at hand presents specific features worth being analyzed in detail. Fig. 5b), when an air ligament, referred to as intact ligament above, remains attached to the outer wall of the exit slit through a small protuberance or perturbation. The latter protuberance appears due to the fact that the initial thickness of the intact ligament, hl(0), is of the order of the distance between the outer sides of the air injector walls, but is slightly inflated when the previous bubble detaches, hl(0) P Ho, as explained in more detail below. Since the gas ligament is exposed to a negative gauge pressure, as deduced from Fig. 5a, the perturbation generated at the joint between the intact ligament and the injection nozzle evolves forming an incipient neck shortly thereafter (Fig. 5c). Both the incipient neck, as well as the tip of the air ligament, move downwards with a velocity very close to the water velocity, as can be qualitatively observed in Fig. 5b–h, and quantitatively confirmed in Fig. 6a, which shows that the dimensionless streamwise position of the neck, xn/Ho, increases linearly with s, with a slope of approximately equal to one. Moreover, both Fig. 5b–h and Fig. 6b show that, for values of s J 7, the semi-thickness of the neck, yn/Ho, decreases continuously with time while it moves downstream. Therefore, an increasing pressure drop takes place across the neck, as a consequence of the kinetic energy dissipation associated to the air jet discharging into the forming bubble. Notice that the collapsing neck can also be interpreted as a closing valve. Thus, in order to supply a constant air flow rate through the injection nozzle, as imposed both in the experiments and in the numerical simulations, the pressure at the exit of the gas injector must increase, as indeed observed in Fig. 5a for s > 12. As time progresses, the collapse process becomes more violent: while the neck continues to accelerate inwards, the gas pressure at the exit increases very quickly (from s  25) and, consequently, the air sheet begins to inflate upstream of the neck prior to the detachment of the bubble forming downstream. This inflation mechanism takes most of the flow rate supplied through the injection nozzle, decreasing the gas flow rate feeding the forming bubble through

the neck, Qn, (Fig. 6c). During the last instants, the ligament inflation process induces a fast decrease in the exit pressure (from s  29), that eventually becomes negative at the break-up time. Note that the air sheet expansion process can be observed only during the latest instants of the bubbling period in Fig. 5b–h. Since the process is periodic, Fig. 5b represents both the end of the bubble formation period, as well as the beginning of the process leading to the next bubble. Notice also from Fig. 6c that the dimensionless flow rate through the neck, Qn(s)/Qa, takes an approximately constant value Qn/Qa ’ 0.7 for s [ 25. Both the constancy of Qn, as well as its value, can be easily explained by the fact that the ligament thickness is hl ’ Ho during the neck formation stage, while the neck moves dowsntream at a velocity equal to uw. Thus, global volume conservation requires that Qn = Qa  2b uwHo and, consequently, Qn/Qa = 1  K/b ’ 0.71 for the particular conditions of Fig. 6. Another important aspect of the bubble formation process revealed by Fig. 5a is the fact that the negative pressure at the exit of the injection nozzle persists during most of the bubble formation time. Thus, when the neck starts to form (Fig. 5b), the gas sheet does not inflate significantly upstream of the neck, as happens in the cylindrical configuration (Sevilla et al., 2005a). However, the portion of the air film located downstream of the neck, where the pressure is slightly larger than the ambient pressure, starts to expand. Therefore, during this phase, the gas domain observed in Fig. 5b–h can be divided into two different spatial regions: (1) an initial part, upstream of the neck, which stretches downstream at the water velocity without inflating transversely and where, correspondingly, the interface remains nearly parallel to the air injector walls and (2) a second region, located downstream of the neck, that inflates transversely while moving downstream. Note that this second part of the gas stream will form the new bubble after the collapse of the neck. For the parameter values corresponding to Figs. 5 and 6, namely We = 38 and K = 0.147, the region downstream of the neck has an approximately constant length, lb, corresponding to the fact that the tip of the bubble moves at the water velocity in this particular case. However, a similar analysis performed for other parameter

112

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

1

(a)

1

15

0.4

5

0.2 10

15

20

25

30

0.6

0.6

10

5

(c)

0.8

0.8

20

yn/Ho

xn/Ho

25

0

(b)

Qn/Qa

30

0

0.4 0.2 5

10

τ

15

τ

20

25

30

0

5

10

15

20

25

30

τ

Fig. 6. (a) Dimensionless downwards position of the neck, xn/Ho, (b) dimensionless transversal position of the neck yn/Ho, and (c) dimensionless air flow through the neck, Qn/ Qa, as functions of s, for the same conditions as in Fig. 5.

  Fig. 7. Time evolution of the contours of dimensionless gauge pressure, ðp  p0 Þ= qa u2a , inside the air stream, for the same conditions as in Fig. 5.

combinations revealed that the latter conclusion is not generally valid. Indeed, although the downwards velocity of the neck is always very close to uw, the bubble tip velocity increases as K decreases giving, for instance a value of approximately 1.2uw for K = 0.1. Therefore, for small enough values of K, the bubble length lb increases with time during the bubble formation process. The time evolution of the contours of dimensionless gauge pres  sure, ðp  p0 Þ= qa u2a , is displayed in Fig. 7, for the same reference case as in Figs. 5 and 6. It can be deduced from Fig. 7 that, for s > 25, the collapse of the neck produces an increase in pressure upstream of the neck, and a corresponding decrease inside the forming bubble. The latter pressure decrease is produced due to the fact that, since the flow rate filling the bubble, Qn, decreases with time (see Fig. 6c), the growth rate of the forming bubble is faster than the gas flow rate feeding it through the collapsing neck. Consequently, the pressure inside the forming bubble must decrease to slow down the rate of volume increase. Finally, since Qn suddenly decreases during the final instants, most of the air flow rate supplied through the injection nozzle is used to fill the gas volume upstream of the neck, what produces a rapid increase of the pressure and, thus, a fast transverse acceleration of the liquid stream surrounding the gas stem. However, since the growth rate of volume induced by the liquid transverse acceleration is faster

than the air flow rate supplied, Qa, the gas pressure must suddenly decrease when the neck collapses, providing the negative pressure observed in Fig. 5. A quantitatively detailed picture of the pressure field in the air   stream is provided by Fig. 8a, where the value of ðp  p0 Þ= qa u2a upstream of the neck is represented as function of x/Ho for several values of s. Fig. 8a reveals that, prior to the ligament inflation phenomenon discussed above, the air pressure increases both with time and with downstream position. Nevertheless, notice that from the point of view of the bubbling frequency, the most relevant information is provided by the gauge pressure at the position of the neck, p[x = xn(s), y = 0]  p0 = pn(s)  p0, represented in Fig. 8b for two different cases, namely We, K = 38.8, 0.147 and 24.9, 0.107, respectively. Indeed, the fact that pn  p0 < 0 "s explains the transverse acceleration of the liquid towards the symmetry plane in the neck region, responsible for the formation of the neck and its subsequent collapse. In addition, it is clear that the characteristic bubbling time may be obtained through the y-component of the liquid momentum equation, once the characteristic negative pressures shown in Fig. 8b are estimated as a function of the governing parameters. In turn, to estimate the characteristic pressures, the dominant suction mechanism must be identified, as discussed in detail in Sections 3.3 and 3.4.

113

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

0

-0.5

0

-1 τ = 6.41 11.07 15.73 20.84 25.05 27.38

We=38.8, Λ=0.147

-0.5

0

5

10

15

20

We=38.8, Λ=0.147 We=24.9, Λ=0.107

25

30

0

5

10

15

-1.5

20

τ

x/Ho

2

0.5

2

(p-p0)/(ρaua )

(b)

(pn-p0)/(ρaua )

(a)

1

25

30

    Fig. 8. (a) Downstream evolution of ðp  p0 Þ= qa u2a at the symmetry plane, y = 0, up to the neck position, for several values of s. (b) Temporal evolution of ðp  p0 Þ= qa u2a , evaluated at the instantaneous position of the neck, (x = xn(s), y = 0).

(b)

uw=1.6 m/s 1.8 2.0 2.3 2.5 2.7 2.9

300

400 350 300

250

250

200

200

150

150

100

1.6

1.8

2

2.2

2.4

2.6

2.8

3

u w (m/s)

15

17

19

21

23

25

f b (Hz)

350

f b (Hz)

(a)

ua=15.3 m/s 17.0 18.7 20.4 22.1 23.8

400

100

u a (m/s)

Fig. 9. (a) Bubble break-up frequency, fb, as a function of the water velocity, uw, for several values of the mean air velocity, ua and (b) dependence of fb on ua for several values of uw.

3.2. Bubble characteristics The bubbling frequency, fb, was experimentally determined, for several values of the air and water velocities, by processing the images recorded with a high-speed digital camera. The reproducibility of the results was guaranteed by performing the experiments several times on different days. The frame rate was varied between 10,000 s1 and 20,000 s1, thus ensuring a good time resolution in the determination of the bubble formation frequency, which was found to be of the order of a few hundred bubbles per second. Fig. 9a shows the dependence of the experimental bubbling frequency on the water velocity, uw, for several values of the mean air velocity, ua. Similarly, Fig. 9b displays fb as a function of ua for different values of uw. Alternatively, since the volume of the bubble can be obtained from the bubbling frequency as Vb = Qa/ fb, similar figures could be shown to analyze the dependence of Vb on uw and ua. As can be observed in Fig. 9a, fb increases monotonically with uw for a fixed value of ua. Moreover, Fig. 9b reveals that the same conclusion holds for the dependence of fb on ua when uw is kept constant. The bubbling frequency was also determined from the numerical simulations. In this case, the temporal evolution of the pressure inside the air stream, at the exit of the nozzle, (x = y = 0), was extracted during the numerical bubbling process, for different values of the air and the water velocities, and a spectral analysis was applied to the signals to determine the bubble formation frequencies. Fig. 10a shows an example of the time series of the pressure obtained for

the reference case, where We = 38 and K = 0.147. The associated power spectral density is displayed in Fig. 10b, showing that the bubbling frequency is fb = 196 Hz in this particular case. The numerically determined frequencies are compared with the experimental ones in Fig. 11 for four intermediate values of the water velocity. It can be observed in this figure that the experimental and numerical bubbling frequencies show an excellent agreement, following the same increasing trend with the air velocity. Consequently, since the numerical simulations are two-dimensional, it can be concluded that the periodic bubbling phenomenon can be described without performing more complex and time consuming three-dimensional simulations. Another relevant quantity obtained from the experiments was the intact ligament length, li, defined above (see Fig. 2). Figs. 12a and b show the dependence of li on the velocities of the water and air streams, respectively. While it can be deduced from Fig. 12a that the intact length barely depends on the water velocity, Fig. 12b clearly reveals a decreasing trend of li with the air velocity. This result is consistent with the evolution of the bubbling frequencies with ua and uw observed in Fig. 11. The intact ligament corresponds to the position where the neck collapses and, since it propagates downstream at the water velocity, li ’ uw/fb. However, given that the bubbling frequency increases almost linearly with uw (see Fig. 11), the dependence of li with uw is expected to be very weak, in agreement with the results of Fig. 12a. On the other hand, since fb increases linearly with ua, li ’ uw/fb / 1/ua when uw is kept constant, in agreement with Fig. 12b. From these

114

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

(b)

PSD

pressure signal

(a)

0

1.14 1.15 1.16 1.17 1.18 1.19

1.2

1.21

0

500

time (s)

1000

2000

1500

frequency (Hz)

Fig. 10. (a) Gauge pressure at the outlet as a function of time, extracted from a numerical simulation performed for We = 38 and K = 0.147. (b) Power spectral density of the pressure signal represented in (a), from which the numerical bubbling frequency can be easily extracted.

ua (m/s) 19

17

ua (m/s) 21

23

25

13

15

19

17

21

23

25 300

(a)

(b)

uw = 1.8 m/s

uw = 2.0 m/s

250

200

200

150

150

experiments

experiments

simulations

simulations

100

100

300

300

(d)

(c) fb (Hz)

250

fb (Hz)

fb (Hz)

250

15

uw = 2.5 m/s

uw = 2.3 m/s

250

200

200

150

fb (Hz)

300

13

150

experiments

experiments

simulations

simulations

100

100 13

15

17

19

21

23

25

ua (m/s)

13

15

17

19

21

23

25

ua (m/s)

Fig. 11. Comparison of the experimental and numerical bubbling frequencies as functions of the air velocity, corresponding to four different values of the water velocity.

considerations, although not presented here to avoid showing redundant data, it is also clear that the volume of the bubbles, given by Vb = Qa/fb, decreases with the water velocity but it is barely affected by the air velocity, in agreement with our experimental and numerical results. In addition, Fig. 13 shows a comparison of the numerical and experimental intact lengths obtained for different water and air velocities. A fairly good quantitative agreement is found in Fig. 13, specially at lower water velocities, for which the maximum relative differences are about 10%. Note that the same general trend is observed in numerics and experiments, the intact length

being monotonically reduced, in both cases, as the air velocity increases. 3.3. Scaling of the bubbling frequency A simple scaling law for the bubbling frequency will be proposed in this section, based on the evidences provided in Section 3.1. As has already been pointed out in Section 3.1, the y-component of the liquid momentum equation can be used to provide a simple estimate of the bubbling frequency, once the dominant mechanism

115

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

(a)

18

li (mm)

(b)

ua=15.3 m/s 17 18.7 20.4 22.1 23.8

uw=1.6 m/s 1.8 2.0 2.3 2.5 2.7 2.9

22

18

14

14

10

10

6

6 1.6

1.8

2

2.2

2.4

2.6

2.8

14

16

18

uw (m/s)

20

22

li (mm)

22

24

ua (m/s)

Fig. 12. Dependence of the experimental intact length, li, on (a) the water velocity, uw and (b) the air velocity, ua.

ua (m/s) 15

(a)

20

30

20

(b)

uw = 1.8 m/s

20

25

30

uw = 2.0 m/s

experiments

experiments

simulations

simulations

20

15

15

10

10

5

5

25

25

(c)

uw = 2.3 m/s

20

l i (mm)

(d)

uw = 2.5 m/s

experiments

experiments

simulations

simulations

20

15

15

10

10

5

5

10

15

20

25

30

10

15

ua (m/s)

20

25

li (mm)

l i (mm)

25

15

li (mm)

10

ua (m/s) 10

30

ua (m/s)

Fig. 13. Dependence of li on ua, as obtained from the numerical simulations (solid symbols), and measured from the experiments (open symbols), for different values of uw.

governing the generation of negative pressures in the air stream has been identified. To that end, following Sevilla et al. (2005a) and Gordillo et al. (2007), we will only consider the transverse dynamics of the water stream, in response to the pressure fluctuations induced by the air stream. Specifically, let us assume that, in a frame of reference moving at the water velocity, uw, the liquid velocity field is purely transverse, vw = v ey. Note that this hypothesis is equivalent to assume long wavelength perturbations of the air–water interface. Then, from the liquid continuity equation, it is immediately deduced that v is only a function of time, v = v(t), and is thus related to the instantaneous semi-thickness of the neck, hn(t), through v ¼ h_ n ðtÞ, where a dot will hereafter denote time

derivative. With these assumptions, the y-component of the momentum equation in the water domain writes

qw h€n ¼ 

@pw ; @y

ð11Þ

which can be integrated from y = hn to y = hw to give,

qw ðHw  Ho Þh€n ¼ pðhn Þ  p0 ;

ð12Þ

where p(hn) is the pressure on the liquid side of the interface at the neck position, p0 is the atmospheric pressure, and it has been assumed that the thickness of the water stream is constant within

116

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

our long-wavelength approximation, that is, hw  hn = Hw  Ho. Indeed, since the liquid velocity is independent of y, both the inner and outer interfaces are simultaneously displaced with the same transverse velocity. From Eq. (12) it is deduced that the neck formation and its subsequent collapse require negative values of p(hn)  p0. Indeed, neglecting surface tension effects, p(hn) = pn, and thus the results of Fig. 8b indicate that the negative air pressures in the neck region are responsible for the inwards acceleration of the neck through Eq. (12), leading to the bubble formation process. A detailed discussion on the role of surface tension forces in our configuration will be provided below in Section 3.4, where it is argued that, although surface tension effects can play an important role during the collapse stage of the bubble, their role is subdominant during the neck formation stage. Thus, since our theoretical objective in the present work is just to provide a simple estimate of the characteristic bubbling time, tc, and keeping in mind that the neck formation time is considerably longer than the corresponding collapse time (see Fig. 5), the characteristic bubbling time, tc, may be obtained from Eq. (12) as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qw ðh  1ÞH2o ; tc ¼ Dpc

ð13Þ

where h = Hw/Ho ’ 5.27 in our experiments and numerical simulations, and Dpc is the characteristic value of the negative gauge pressures in the air stream during the bubbling process. It is important to emphasize that Eq. (13) is valid irrespective of the detailed processes leading to suction inside the air stream. However, to close the model for the bubbling frequency, an explicit expression for Dpc as a function of the governing parameters must be specified. In fact, although we have identified several physical mechanisms that may contribute to Dpc, as discussed in Section 3.4, in this section we will focus on the only mechanism able to explain the negative pressures observed in Fig. 5a, namely the sudden expansion of the air stream from the inner wall thickness, Hi, to the outer one, Ho. Indeed, as mentioned above, the lateral expansion of the air stream is due to the fact that the contact line is pinned at the outer edge of the wall separating the air and water streams. Thus, considering that the liquid pressure in the neck region is approximately atmospheric during the neck formation stage, the air pressure at the exit must be negative due to the pressure loss associated to the planar sudden expansion,

pe  p0 ¼ qa u2a bð1  bÞ ) Dpc ¼ qa u2a bð1  bÞ:

ð14Þ

Eq. (14) predicts that the characteristic negative pressures at the exit, scaled with twice the air dynamic pressure, qa u2a , should depend only on the wall thickness through the value of the expansion parameter b = Hi/Ho, which takes the constant value b ’ 0.52 in all the experiments and numerical simulations performed in this work.   Specifically, Eq. (14) gives a value of ðpe  p0 Þ= qa u2a ’ 0:25, a prediction which is in remarkable agreement with the numerical results represented in Figs. 5a and 8. Substituting Eq. (14) in Eq. (13), the following expression for the characteristic time is obtained,

Ho tc ¼ ua

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qw h  1 : qa bð1  bÞ

ð15Þ

Taking into account the fact that the only pressure mechanism considered in Eq. (15) is the sudden expansion of the air stream, it is clear that it leaves out many relevant aspects that also contribute to the time evolution of the pressure at the propagating neck. In particular, Eq. (15) disregards the Venturi effect due to the quasisteady gas acceleration at the neck, i.e. the Kelvin–Helmholtz instability due to the existence of shear at the air–water interface (Sevilla et al., 2005a; Gordillo et al., 2005), as well as the suction associated to the bubble growth at a constant flow rate (Gordillo

et al., 2007), that will be discussed in Section 3.4. Thus, the time given by Eq. (15) may be interpreted as the characteristic neck formation time. Nevertheless, we shall use it below to scale the total bubbling time, on the basis that the collapse stage is considerably shorter. It is important to emphasize that Eq. (15) is, by construction, independent of the liquid velocity, in contradiction with the results shown in Fig. 9. Although, as discussed in Section 3.4, this fundamental limitation can be overcome by considering additional physical effects, it would lead to scaling laws considerably more complex than Eq. (15), thus strongly compromising the simplicity of the approach. Fig. 14 displays the nondimensional experimental and numerical frequencies, fbtc, as a function of K, for several values of We, with tc given by Eq. (15). It can be deduced from Fig. 14 that 0.5 [ fbtc [ 0.8 indicating that, indeed, the characteristic time given by Eq. (15) roughly captures the dependence of the bubble break-up frequency in a wide range of values of the two main control parameters, namely 0.05 [ K [ 0.2 and 10 [ We [ 60. Since fbtc < 1, the value of t1 overestimates the bubbling frequency, in c correspondence with the fact that the values of Dpc, given by Eq. (14), overestimate the air suction at the neck position during most of the bubbling time, as clearly seen in Fig. 8. Notice also that the scatter observed in Fig. 14 is not random, the data being organized in approximately constant series corresponding to the different values of We or, equivalently, to different values of the liquid velocity, such that fbtc increases monotonically with We. The latter behaviour can be simply understood as a consequence of the fact that fb increases with uw, as shown in Fig. 9. Therefore, improving the scaling law (15) would require a proper modeling of the dependence of fbtc on uw, by taking into account the contributions of the Venturi and suction-by-growth mechanisms discussed in the next section, to Dpc. Nevertheless, the equation

ua fb ¼ 0:65 Ho

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qa bð1  bÞ ; qw h  1

ð16Þ

may be used for practical engineering or design purposes. Let us also point out that, in Eq. (16), the relative wall thickness parameter b plays an essential role, since the only physical mechanism behind Dpc in Eq. (14) is the sudden expansion of the air stream from Hi to Ho. Consequently, in designs where the relative wall thickness is either very large, b  1, or very small, 1  b  1, alternative physical mechanisms, like those discussed in the next section, would become an essential ingredient of the theoretical description.

2

1.5

fbtc 1

We = 21.6 (experiment) 26.6 32.2 38.3 52.2 21.6 (simulation) 26.0 32.0 39.5 56.9

0.5

0 0.05

0.1

Λ

0.15

0.2

Fig. 14. Dependence of the dimensionless experimental and numerical bubbling frequencies, fbtc, on the velocity ratio, K, for several values of We. The characteristic time tc is given by Eq. (15).

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

3.4. Other mechanisms leading to suction in the air stream It is important to point out that, besides the sudden expansion mechanism used in the previous section to scale the bubbling frequency observed in our experimental set-up, other physical effects would also contribute to the bubble formation process, whose relative importance will depend on the planar or cylindrical geometry of the bubble, as well as on details of the injection system, like the thickness of the walls separating the air and water streams at the nozzle exit. Let us first consider the influence of surface tension on the bubble formation process. In the two-dimensional geometry under consideration, the only contribution of surface tension forces comes from the longitudinal curvature, so that the pressure jump across the inner air–water interface at the position of the neck is 00 00 given by pðhn Þ  pn ¼ rhn , where hn is the interfacial curvature evaluated at the bubble neck. Thus, including surface tension forces, the air suction pressure responsible for the generation and subsequent collapse of the neck, namely p0  p(hn), can be esti00 00 mated as ðp0  pn Þ  rhn . Taking into account that hn > 0, surface tension effects have a stabilizing influence in the formation of planar bubbles (Bolaños-Jiménez et al., 2011), in contrast with their destabilizing effect in the cylindrical configuration, where the radial curvature at the cylindrical neck (Sevilla et al., 2005a) accelerates the collapse process. Note that, this line of reasoning is consistent with the results shown in Fig. 14, since the dimensionless frequencies increase with increasing values of We. This fact indicates that surface tension forces stabilize the interface, delaying the bubble formation process. Another physical mechanism of relevance in the case of axisymmetric bubbles is the Bernoulli suction associated to the quasisteady air acceleration through the neck (Sevilla et al., 2005a; Gordillo et al., 2005), also present in the planar configuration under study. Nevertheless, as revealed by Figs. 6 and 8b, this suction mechanism is not active until the last stages of the process, and cannot explain the negative pressures prevailing near the outlet, that are responsible for the formation of the neck (Gordillo et al., 2007). Moreover, it can be assumed that, as revealed by Figs. 6 and 8, the neck collapse stage is considerably shorter than the neck formation stage. Therefore, in cases where the wall thickness is very small, 1  b  1, for which the sudden expansion mechanism described in Section 3.3 does not play a role, alternative phenomena leading to negative gauge pressures in the air stream must be identified to properly estimate the bubbling time. In the axisymmetric case, a mechanism alternative to the Bernoulli suction at the neck was identified in Gordillo et al. (2007), where it was argued that filling a bubble at constant flow rate implies a radial deceleration, leading to negative pressures inside the forming bubble. Although the latter conclusion is valid for cylindrical bubbles of constant downstream extent, in the case of planar bubbles this mechanism requires that the length of the bubble increases with time. Indeed, in the case at hand the volume conservation equation reads

Qa d ¼ qa ¼ ½hl ðtÞxn ðtÞ þ hb ðtÞlb ðtÞ; dt 2b

ð17Þ

where, for convenience, we have defined qa = Qa/2b as half the air flow rate per unit spanwise length. During the neck formation stage, Eq. (17) reduces to

qa ¼ uw Ho þ h_ b ðtÞlb ðtÞ þ hb ðtÞ_lb ðtÞ;

ð18Þ

where it has been taken into account that x_ n ¼ uw and hl ’ Ho. The € ðtÞ, can be obtained transverse acceleration of the forming bubble, h b by imposing the constancy of air flow rate, q_ a ¼ 0, to yield

€ _ _ € ¼  2hb lb þ hb lb : h b lb

117

ð19Þ

Substituting h_ b from Eq. (18) into Eq. (19), one gets

€ _ € ¼  2lb ðq  uw Ho  h _l Þ  hb lb : h b b b a 2 lb lb

ð20Þ

Eq. (20) clearly shows that, in the planar geometry under consideration, the transverse acceleration of the bubble, and thus the presence of pressure disturbances in the neck region due to the growth of the bubble at constant flow rate, requires the length of the forming bubble, lb, to change with time. However, as revealed by Fig. 5, the value of lb is in fact approximately constant with time for the reference case, where the governing parameters are We = 38 and K = 0.147. As already mentioned in Section 3.1, lb(t) becomes a weakly increasing function of time for smaller values of K. Nevertheless, it is clear that this transient mechanism alone cannot explain the observed bubbling behaviour and, in particular, it cannot account for the persistent negative exit pressures represented in Fig. 5a, nor for those observed in Fig. 8a in the vicinity the outlet. However, these negative pressures are responsible for the formation of the neck from the initial interfacial disturbance pinned to the injector wall. Consequently, although it would be possible to obtain an expression for Dpc through the combination of Eq. (20) with the y-component of the liquid momentum equation outside the forming bubble, this will not be attempted in the present work. Finally, let us point out that in the thin-wall limit, 1  b  1, for which the sudden expansion mechanism does not apply, the neck formation process would be governed by the suction pressures induced by the transverse decelerations of the interface, described by Eq. (20). Indeed, prior to the appearance of the neck, the intact ligament would increase in length according to _lb ’ uw , and the neg€ stemming from this stretching effect through Eq. ative values of h b (20) would lead to a suction pressure in the air ligament, causing the generation of a neck that would propagate downstream at the liquid velocity and produce a new bubble after collapsing. 4. Conclusions We have studied the dynamics of a plane air film surrounded by two co-flowing water sheets discharging in stagnant air by means of experiments and numerical simulations. Bolaños-Jiménez et al. (2011) showed that this configuration presents two different flow regimes, a jetting and a bubbling regime, depending on the Weber number, We, and the velocity ratio, K. However, in the present work we have mainly focused on the bubbling regime with the aim of describing the bubble formation process. Thus, for a fixed value of the initial liquid-to-gas thickness ratio, h = hw/ha ’ 5.27, we have systematically varied the air and water flow rates, inducing the periodic formation of bubbles as a result of the two-dimensional breakup of the air sheet into individual elongated bubbles. The bubble formation event has been described as a two-stage process: a neck formation stage and a neck collapse stage. The periodic bubbling process starts with the pinch-off of the previous bubble, that leaves an initial lump of air, of length li and slightly thicker than Ho, attached to the outer wall of the air nozzle through a small protuberance. Thus, the gas stream suffers a sudden expansion from the inner thickness of the air nozzle, Hi, to the outer one, Ho, generating a negative gauge pressure inside the air lump. Consequently, as time evolves, the initial protuberance propagates downstream at the water velocity and begins to develop towards the symmetry plane, forming the bubble neck. The persistent negative pressure makes the neck to start closing while the newly forming bubble grows at a constant rate downstream of the neck. When the neck begins to collapse, the flow rate feeding the

118

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119

forming bubble decreases, increasing the pressure upstream of the neck. This increase of pressure produces the slight inflation of the air lump that remains attached to the air nozzle when the neck collapses, and the forming bubble pinches-off. On the other hand, just before the bubble pinches-off, since the flow rate filling the forming bubble through the collapsing neck decreases, the pressure inside the bubble decreases to slow down its volume growth rate. Regarding the bubbling frequency, it has been experimentally measured by means of image processing, finding a monotonic increase with both, the water and the air velocities. This result indicates that the bubble volume, obtained from the bubbling frequency as Vb = Qa/fb, decreases as the water velocity increases but does not vary with the air velocity. The dependence of the bubbling frequency on ua and uw has also been obtained numerically, showing a remarkable agreement with the experimental results. Furthermore, the numerical simulations have been used to describe the time evolution of the pressure inside the air stream during the bubble formation process. In fact, as mentioned above, the numerics corroborate the negative values of the gauge pressure at the exit of the nozzle during the neck formation stage due to the sudden expansion that the air stream suffers when it discharges from the air nozzle of thickness Hi to the forming bubble, attached to the outer wall, of thickness Ho. In addition, the length of the intact ligament, li, has been also determined numerically, finding that, although it barely depends on water velocity, it decreases almost linearly with the air velocity. The numerical predictions for the intact length agree well with the experiments, showing the same trend. Such dependence on the air velocity is consistent with the linear variation of the bubbling frequency with ua at constant uw. Finally, we have proposed a simple scaling law for the bubbling frequency taking into account the effect of the wall thickness of the air nozzle on the pressure decrease in the air stream, when it discharges from a channel of thickness Hi to the forming bubble of thickness of the order of Ho, at the nozzle exit. The model indicates that the dimensionless frequency does not depend on K as corroborated from the experimental and numerical results, which showed that the bubbling frequency increases linearly with both, the air and the water velocities. The dimensionless frequencies given by the model compare fairly well with those obtained from the experiments and the numerics, and can be used for engineering and design purposes. However, it fails to capture the dependence of the bubbling frequency on the Weber number. Thus, a more complex model, taking into account the pressure decrease provided by other mechanisms, should be elaborated to obtain more accurate results, especially in the cases where the relative wall thickness is either very large, b  1, or very small, 1  b  1. Nevertheless, a quick inspection of Fig. 14 indicates that the dimensionless frequency, fbtc, slightly increases with the Weber number and, although not based on first principles, rescaling fbtc with the squared root of We collapses the data on a single curve given by fbtc  0.11 We1/2, with tc provided by Eq. (15). Acknowledgments This work has been supported by the Spanish MINECO, Subdirección General de Gestión de Ayudas a la Investigación, and European Funds under projects number DPI2011-28356-C03-02 and DPI2011-28356-C03-03. Financial support from the University of Jaén under project UJA2011/12/50 is also acknowledged. References Arias, S., Legendre, D., González-Cinca, R., 2011. Numerical simulation of bubble generation in a T-junction. Comput. Fluids 56, 49–60. Berberovic´, E., van Hinsberg, N., Jakirlic´, S., Roisman, I., Tropea, C., 2009. Drop impact onto a liquid layer of finite thickness: dynamics of the cavity evolution. Phys. Rev. E 79, 036306.

Bolaños-Jiménez, R., Sevilla, A., Gutiérrez-Montes, C., Sanmiguel-Rojas, E., MartínezBazán, C., 2011. Bubbling and jetting regimes in planar co-flowing air–water sheets. J. Fluid Mech. 682, 519–542. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modelling surface tension. J. Comput. Phys. 100, 335–354. Chakraborty, I., Biswas, G., Ghoshdastidar, P., 2011. Bubble generation in quiescent and co-flowing liquids. Int. J. Heat Mass Trans. 54, 4673–4688. Chuang, S.C., Goldschmidt, V.W., 1970. Bubble formation due to a submerged capillary tube in quiescent and co-flowing streams. Trans. ASME D: J. Basic Eng. 92, 705–711. Crank, J., Nicolson, E., 1947. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Proc. Camb. Philos. Soc. 43, 50–67. Davidson, J.F., Schuler, B.O.G., 1960. Bubble formation at an orifice in an inviscid liquid. Trans. Inst. Chem. Engrs. 38, 335–342. Enright, D., Fedkiw, R., Ferziger, J., Mitchell, I., 2002. A hybrid particle level set method for improved interface capturing. J. Comput. Phys. 183, 83–116. Fuster, D., Bagué, A., Boeck, T., Le Moyne, L., Leboissetier, A., Popinet, S., Ray, P., Scardovelli, R., Zaleski, S., 2009. Simulation of primary atomization with an octree adaptive mesh refinement and VOF method. Int. J. Multiphase Flow 35, 550–565. Gordillo, J.M., Sevilla, A., Martínez-Bazán, C., 2007. Bubbling in a co-flow at high Reynolds numbers. Phys. Fluids 19, 077102. Gordillo, J.M., Sevilla, A., Rodríguez-Rodríguez, J., Martínez-Bazán, C., 2005. Axisymmetric bubble pinch-off at high Reynolds numbers. Phys. Rev. Lett. 95, 194501. Gueyffier, D., Li, J., Nadim, A., Scardovelli, S., Zaleski, S., 1999. Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows. J. Comput. Phys. 152, 423–456. Hauke, G., Dopazo, C., Lozano, A., Barreras, F., Hernández, A.H., 2001. Linear stability analysis of a viscous liquid sheet in a high-speed viscous gas. Flow Turbul. Combust. 67, 235. Herrada, M., Lopez-Herrera, J., Vega, E., Montanero, J., 2011. Numerical simulation of a liquid bridge in a coaxial gas flow. Phys. Fluids 23, 012101. Herrmann, M., 2008. A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids. J. Comput. Phys. 227, 2674–2706. Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201. Issa, R., 1986. Solution of the implicitly discretized fluid flow equations by operatorsplitting. J. Comput. Phys. 62, 40–65. Jasak, H., 1996. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. PhD thesis, Imperial College of Science, Technology and Medicine, London. Kumar, R., Kuloor, N.R., 1976. The formation of bubbles and drops. Chem. Eng. Sci. 31, 453. Lax, P., Wendroff, B., 1960. Systems of conservation laws. Commun. Pure Appl. Math. 13, 217–237. Lezzi, A.M., Prosperetti, A., 1991. Stability of an air film in a liquid flow. J. Fluid Mech. 226, 319–347. Li, X., Bhunia, A., 1996. Temporal instability of plane gas sheets in a viscous liquid medium. Phys. Fluids 8, 103–111. Li, X., Bhunia, A., 1997. Instability of plane compressible gas sheets. Acta Mech. 123, 117–133. Longuet-Higgins, M.S., Kerman, B.R., Lunde, K., 1991. The release of air bubbles from an underwater nozzle. J. Fluid Mech. 230, 365–390. López-Pagés, E., Dopazo, C., Fueyo, N., 2004. Very-near-field dynamics in the injection of two-dimensional gas jets and thin liquid sheets between two parallel high-speed gas streams. J. Fluid Mech. 515, 1–31. Lozano, A., Barreras, F., 2001. Experimental study of the gas flow in a air-blasted liquid sheet. Exp. Fluids 31, 367–376. Lozano, A., Barreras, F., Hauke, G., Dopazo, C., 2001. Longitudinal instabilities in an air-blasted liquid sheet. J. Fluid Mech. 437, 143–173. Lozano, A., Barreras, F., Siegler, C., Löw, D., 2005. The effects of sheet thickness on the oscillation of an air-blasted liquid sheet. Exp. Fluids 39, 127–139. Menard, T., Tanguy, S., Berlemont, A., 2007. Coupling level set/VOF/ghost fluid methods: validation and application to 3D simulation of the primary break-up of a liquid jet. Int. J. Multiphase Flow 33, 510–524. ~uz, H.N., Prosperetti, A., 1993. Dynamics of bubble growth and detachment from Og a needle. J. Fluid Mech. 257, 111–145. Park, J., Huh, K., Li, X., Renksizbulut, M., 2004. Experimental investigation on cellular breakup of a planar liquid sheet from an air-blast nozzle. Phys. Fluids 16, 625. Pianet, G., Vincent, S., Leboi, J., Caltagirone, J., Anderhuber, M., 2010. Simulating compressible gas bubbles with a smooth volume tracking 1-Fluid method. Int. J. Multiphase Flow 36, 273–283. Rusche, H., 2002. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. PhD thesis, Imperial College of Science, Technology and Medicine, London. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Ann. Rev. Fluid Mech. 31, 567–603. Sevilla, A., Gordillo, J.M., Martínez-Bazán, C., 2002. The effect of the diameter ratio on the absolute and convective instability of free co-flowing jets. Phys. Fluids 14, 3028–3038. Sevilla, A., Gordillo, J.M., Martı´nez-Bazán, C., 2005a. Bubble formation in a coflowing air–water stream. J. Fluid Mech. 530, 181–195. Sevilla, A., Gordillo, J.M., Martı´nez-Bazán, C., 2005b. Transition from bubbling to jetting in a coaxial air–water jet. Phys. Fluids 17, 018105.

C. Gutiérrez-Montes et al. / International Journal of Multiphase Flow 50 (2013) 106–119 Siamas, G., Jiang, X., 2007. Direct numerical simulation of a liquid sheet in a compressible gas stream in axisymmetric and planar configurations. Theor. Comput. Fluid Dyn. 21, 447–471. Sussman, M., Puckett, E., 2000. A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows. J. Comput. Phys. 162, 301–337. Sussman, M., Smereka, P., Osher, S., 1994. A level set method for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146–159.

119

Unverdi, S., Tryggvason, G., 1992. A front-tracking method for viscous incompressible multi-fluid flows. J. Comput. Phys. 100, 25–37. Van Leer, B., 1974. Towards the ultimate conservative differencing scheme – II. Monotonicity and conservation combined in a second-order scheme.. J. Comput. Phys. 14, 361–370. Zheng, D., He, X., Che, D., 2007. CFD simulations of hydrodynamic characteristics in a gas–liquid vertical upward slug flow. Int. J. Heat Mass Transfer 50, 4151– 4165.