Next Article in Journal
Stationary Mach Configurations with Pulsed Energy Release on the Normal Shock
Next Article in Special Issue
Header Shape Effect on the Inlet Velocity Distribution in Cross-Flow Double-Layered Microchannel Heat Sinks
Previous Article in Journal
Spin-Up from Rest of a Liquid Metal with Deformable Free Surface in a Cylinder under the Influence of a Uniform Axial Magnetic Field
Previous Article in Special Issue
Computational Study of the Dynamics of the Taylor Bubble
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Shapes and Rise Velocities of Single Bubbles in a Confined Annular Channel: Experiments and Numerical Simulations

by
Andrea Cioncolini
1,* and
Mirco Magnini
2,*
1
Department of Mechanical, Aerospace and Civil Engineering, University of Manchester, Oxford Road, Manchester M13 9PL, UK
2
Department of Mechanical, Materials and Manufacturing Engineering, University of Nottingham, University Park, Nottingham NG7 2RD, UK
*
Authors to whom correspondence should be addressed.
Fluids 2021, 6(12), 437; https://doi.org/10.3390/fluids6120437
Submission received: 11 November 2021 / Revised: 26 November 2021 / Accepted: 29 November 2021 / Published: 2 December 2021

Abstract

:
Shapes and rise velocities of single air bubbles rising through stagnant water confined inside an annular channel were investigated by means of experiments and numerical simulations. Fast video imaging and image processing were used for the experiments, whilst the numerical simulations were carried out using the volume of fluid method and the open-source package OpenFOAM. The confinement of the annular channel did not affect the qualitative behavior of the bubbles, which exhibited a wobbling rise dynamic similar to that observed in bubbles rising through unconfined liquids. The effect of the confinement on the shape and rise velocity was evident; the bubbles were less deformed and rose slower in comparison with bubbles rising through unconfined liquids. The present data and numerical simulations, as well as the data collected from the literature for use here, indicate that the size, shape, and rise velocity of single bubbles are closely linked together, and prediction methods that fail to recognize this perform poorly. This study and the limited evidence documented in the literature indicate that the confinement effects observed in non-circular channels of complex shape are more complicated than those observed with circular tubes, and much less well understood.

1. Introduction

The dynamics of single bubbles rising through stagnant liquids is of pivotal importance for the fundamental understanding of two-phase bubbly flows, which are relevant in a number of applications including gas–liquid reactors, bubbly columns, nuclear reactors, heat exchangers, and environmental flows. Gas–liquid contacting, which is normally achieved by bubbling a gas into a liquid, is in fact one of the most common operations in the process industry in applications such as absorption, distillation, and froth flotation.
Despite the apparent simplicity, the dynamics of single bubbles rising through stagnant liquids is in reality a rather involved free-boundary problem controlled by the interplay among inertia, buoyancy, viscous, and surface tension forces. When bubbles rise in bounded liquids their dynamics are also affected by the walls of the container. Available experimental studies of single bubbles rising through stagnant liquids were performed in containers of finite dimensions, typically vertical tanks of either circular or square cross-sections. Strictly speaking, therefore, wall effects were always present to a greater or lesser extent. However, it is currently accepted that, if the dimensions of the horizontal cross-section of the container are much larger than the size of the bubble (indicatively, 10–20 times or more), then wall effects on the bubble dynamics are small or absent. The bubble can therefore be considered unconfined, and the observed bubble dynamics can be regarded as representative of the free bubble rise through a stagnant unconfined liquid. This is the case for the majority of the experimental studies documented to date [1,2,3].
Wall effects were investigated rather extensively for bubbles rising through circular tubes [4,5,6] and through narrow rectangular channels [7,8,9,10,11,12,13]. For circular tubes, wall effects typically cause the elongation of the bubbles in the vertical direction and the alteration of the wake structure behind the rising bubble, resulting in milder bubble deformation and a reduced rise velocity in comparison with the unconfined case. For narrow rectangular channels, bubbles are typically bigger than the channel gap, so that their dynamic is largely controlled by wall effects.
In contrast, the dynamics of single bubbles rising through confined channels of more complex shape, representative or informative of rod bundles or bubbly columns with internals, were not extensively investigated; the only studies documented in the literature appear to be those by Venkateswararao et al. [14] and by Tomiyama et al. [15]. In particular, Venkateswararao et al. [14] measured the rise velocity of air bubbles rising in stagnant water (tap filtered) at ambient conditions through a rod-bundle column comprising sixteen rods of 12.7 mm diameter, arranged in a square lattice, with pitch of 17.5 mm, and installed inside a 88.9 mm diameter circular pipe (corresponding to a subchannel hydraulic diameter of 18.0 mm). Partial rods were placed along the inner surface of the circular pipe to minimize end effects and make the cross-section of the test piece representative of a larger rod bundle. For bubble sizes in the range of 2–8 mm, the measured rise velocities increased with the bubble diameter: from about 0.21–0.22 m/s for 2–3 mm diameter bubbles to about 0.3 m/s for 7–8 mm diameter bubbles. On the other hand, Tomiyama et al. [15] measured the shapes and rise velocities of air bubbles rising in stagnant water (distilled) at ambient conditions through an inner subchannel comprising four rods of 12.0 mm diameter, arranged in a square lattice, with a pitch of 15.2 mm (corresponding to a subchannel hydraulic diameter of 12.5 mm), and installed inside a square pipe. The rods were in contact with the inner surface of the square pipe. For bubble size in the range of 3–6 mm, the measured rise velocities decreased with bubble diameter: from about 0.21–0.22 m/s for 3 mm diameter bubbles to about 0.17–0.18 m/s for 6 mm diameter bubbles.
Therefore, as can be noted, the documented studies of single bubbles rising through confined channels of complex shape are not only a few in number, they also yield conflicting results regarding the bubble rise velocity, which increases with bubble diameter according to the measurements of Venkateswararao et al. [14], whilst it decreases with bubble diameter according to the measurements of Tomiyama et al. [15]. This clearly indicates that more investigations are needed on single bubble rise in non-circular channels of complex shape, thereby creating the motivation for the present work. In this study, we investigated by means of experiments and numerical simulations the shapes and rise velocities of single air bubbles rising through stagnant water inside an annular channel, a geometry not previously considered in single bubble studies. In particular, we used the numerical simulations to help interpret the data and widen the scope of the experimental observations.

2. Experiments

2.1. Experimental Setup

A schematic representation of the experimental setup is provided in Figure 1a,b. The test piece is a vertical annular channel with a hydraulic diameter of 16.9 ± 0.5 mm, comprising a square cross-section transparent pipe (high-temperature-rated clear acrylic) with nominal inner size of 25 mm (measured actual size of 25.3 ± 0.3 mm), which accommodates an internal stainless-steel rod of 10 mm nominal diameter (measured actual diameter of 9.99 ± 0.01 mm). The present experimental setup was specifically designed to carry out three different types of experimental investigations:
  • Single air bubble dynamics in confined stagnant water;
  • Multiple air bubble dynamics in confined stagnant water with annular channel operating as a small bubble column running in batch mode;
  • Pool boiling in water where bubbling is sustained via Joule heating of the internal rod.
Only the results of the single bubble dynamics experiment from Point 1 are documented herein; the results of the other experiments will be communicated separately (the experiments on multiple air bubble dynamics are scheduled for the first and second quarters of 2022, whilst the experiments on pool boiling are scheduled for the last quarter of 2022). Nonetheless, it is important to clarify beforehand the wider scope of the present experimental setup, whose design is a compromise among the different requirements of the experimental investigations listed above. Clearly, this has informed the experimental procedure adopted here.
Regarding the single bubble dynamics experiments considered here, the testing fluids were micro-filtered and high-purity water (conductivity of 0.005 µS/cm) and filtered air from the mains. The water depth level in the annular channel was of 300 ± 1 mm, and was not changed during the tests. Bubbles were generated by blowing air through a capillary pipe (diameter of 1.6 ± 0.1 mm) protruding at the bottom of the annular channel, as schematically illustrated in Figure 1c. The length of the capillary pipe was of 250 mm (corresponding to a length-to-diameter ratio of 156), which is large enough to ensure that the pressure fluctuations due to the formation of the bubbles have a negligible effect on the air flow rate, which can then be taken as constant during the tests. The air flow rate through the capillary pipe was adjusted with a needle valve to yield a bubbling frequency of 0.45 Hz (corresponding to one bubble generated every 2.2 s), which was not varied during the experiments. This yielded bubble diameters of about 4.3 mm. The residence time of the bubbles in the column is below 2 s, thereby ensuring that, when a bubble is generated at the nozzle, the previously generated bubble has already reached the top of the channel. As confirmed by the CFD simulations discussed later, this provides enough time for the liquid to recover between successive bubble passages, so that each bubble rises through a virtually still body of water. The present set-point with a bubbling frequency of 0.45 Hz, therefore, ensures continuous bubbling at a steady rate with a delay between successive bubbles sufficient to allow liquid recovery and, correspondingly, treat each bubble as a single bubble rising through a stagnant liquid.
In the provision for the experiments on multiple bubble dynamics (Point 2 above), the annular channel in Figure 1 was equipped with four independent and identical air injection lines. However, only one air injection line was used for the single bubble experiments reported herein, whilst the other three air injection lines were inactive: the four nozzles, one active and three inactive, are shown in Figure 1a.
The present bubble generation method differs from that commonly employed in single bubble studies, where bubbles are generated with syringe pumps. These latter are clearly more flexible, and allow a range of bubble sizes to be generated from a nozzle of fixed diameter. In contrast, in the present setup, the bubble size does not vary once the air flow rate is fixed. To overcome this limitation, we used CFD simulations (discussed later) to generalize the experimental observations to bubble diameters in the range of 3–6 mm. The present bubble generation method was chosen because it was more representative of bubble columns which, as explained before, were within the broader scope of the experimental facility. One disadvantage of bubble generation with syringe pumps is that the initial deformation of the bubbles is variable: the larger the bubble generated from a nozzle of given diameter, the larger the initial deformation of the bubble [16,17]. Particularly with air–water systems and for bubbles in the range of about 1–10 mm in diameter, where surface tension plays a dominant role on the bubble dynamics, the initial deformation of the bubble can have a profound influence on the subsequent bubble rise: the bigger the initial deformation, the bigger the subsequent rise velocity [16,17,18,19,20]. The consequence is that bubbles of comparable size generated with syringe pumps from nozzles of different diameters can have different initial deformations, and, subsequently, exhibit different rise velocities. This is in fact one of the reasons behind the large scatter in the documented results for the rise velocity of air bubbles in water [1], the other reason being the presence of contaminants in the water. With the bubble generation method used here, the initial deformation of the bubbles is not variable and the scatter in rise velocity is not observed, as described later.
Measurements were carried out at ambient pressure and room temperature (295 ± 2 K). The relevant properties (liquid and gas densities ρ l and ρ g , liquid and gas viscosities μ l and μ g , and surface tension σ ) of the testing fluids were calculated with NIST-REFPROP [21], and their values are provided in Table 1. As noted previously, during the tests, the water depth level in the annular channel was 300 ± 1 mm. The variation of the air density due to the hydrostatic pressure variation along the channel is within a few percent, and thus can be neglected.
Despite its simplicity, justified by the scope of the present study, which is more of fundamental rather than applied character, the present experimental setup can be informative of more complex configurations of direct industrial relevance, such as rod bundles and bubble columns with internals.

2.2. Measurement Methodology

Rising bubbles were recorded using a digital camera (Bonito CL-400 equipped with a Navitar 25 mm Platinum F8 lens—image resolution: 864 × 1168 pixels; recording frequency: 386 fps) providing a spatial resolution of 35.4 µm/pixel (corresponding to 28.2 pixel/mm), appropriate for resolving the present bubbles, which are a few mm in size. As shown in Figure 1b, the digital camera was positioned midway through the annular channel at an elevation of 150 ± 1 mm above the bottom of the channel, so as to image a portion of channel of 40 mm in length from a vertical elevation of 130 mm up to 170 mm. As confirmed by CFD simulations discussed later, at this elevation, the bubbles reached their terminal rise dynamics. Due to the large difference between the densities of water and air (see Table 1), air bubbles in water have very low inertia. Correspondingly, a rise of just a few centimeters is normally sufficient to extinguish the initial transient and reach the terminal dynamics [20,22,23]. Since only one camera was used, the bubbles were characterized using their planar projections, as seen in individual frames. Even though this is normally considered to be sufficiently accurate [22], multiple cameras could clearly provide a more faithful bubble characterization. The error associated with using only one camera in the present case was estimated with CFD simulations, as discussed later.
As illustrated in Figure 2, the image processing methodology adopted here is articulated in four steps:
  • Discrete points are manually digitized along the bubble border (Figure 2b);
  • Cubic spline interpolation is used to increase the density of the points along the bubble border (Figure 2c);
  • The polygon representing the bubble border is converted to a binary image of the bubble (Figure 2d);
  • Size and shape of the bubble are computed.
Operatively, the manual location of the points along the bubble border was accomplished with the free Web-based tool WebPlotDigitizer [24], whilst the other image processing steps were carried out with the free software GNU Octave version 4.2.2 [25], relying exclusively on built-in features for the sake of reproducibility (particularly the functions poly2mask and regionprops of the GNU Octave package ‘image’). Manual digitizing was performed on the raw image by visually locating the bubble interface (note the rather sharp intensity variation across the bubble border in Figure 2a) and by manually digitizing discrete points (see Figure 2b) along the interface. The number of manually located points necessary for a successful bubble shape reconstruction was determined empirically via trial and error to be in the range of 10–15 (see Figure 2b), whilst an increase via cubic spline interpolation to 100 points (see Figure 2c) was found sufficient to yield a binary image (see Figure 2d) with a reasonably smooth boundary that properly captured the bubble shape.
Following common practice, the size and shape of the bubbles were characterized using the equivalent diameter and the aspect ratio, respectively. In particular, the instantaneous equivalent diameter d e q of the bubble is the diameter of the circle with the same area as the projection of the bubble:
d e q ( t ) = 4   A ( t ) π ,
where A is the area of the bubble projection and t is the time label of the image, i.e., the time the image was recorded.
Different approaches have been used in the literature to compute the aspect ratio of the bubble projection. The method adopted here is based on the equivalent ellipse (built-in within the function regionprops of the GNU Octave package ‘image’), which is the ellipse that has the same normalized second central moments as the binary image of the projected bubble. This follows common practice in image analysis, where moments are used to describe the shape of image features by measuring the distribution of pixels with respect to the horizontal and vertical directions [26,27]. A representative example is provided in Figure 2e, where the equivalent ellipse (in green color) is superimposed onto the raw image of the bubble. The instantaneous aspect ratio E of the equivalent ellipse is simply defined as the ratio of the lengths of the minor ( b ) and major ( a ) axes of the ellipse:
E ( t ) = b ( t ) a ( t ) .  
In addition to the aspect ratio, the inclination of the equivalent ellipse was also computed, which corresponds to the angle between the major axis of the ellipse and the horizontal (the angle is positive if oriented counter-clockwise, negative otherwise). Following common practice, the instantaneous bubble rise velocity V r i s e was computed as the vertical velocity of the bubble centroid as follows:
V r i s e ( t ) = z c ( t + Δ t ) z c ( t ) Δ t ,
where z c is the vertical elevation of the centroid of the bubble and Δ t is the time elapsed between successive frames (note that, strictly speaking, Equation (3) is exact only in the limit of Δ t 0 + ; in the present case, however, Δ t is small enough to make the approximation acceptable).
Measuring errors were on the order of 9–10% for the equivalent diameter, 9–10% for the aspect ratio, and 6–7% for the rise velocity; and were estimated by imaging various calibration standards of known shape and size placed at various positions inside the annular channel.
The relevant dimensionless numbers in single bubble dynamics are the Reynolds number R e , the Eötvös number E o , the Weber number W e , and the Morton number M o ; these were computed as follows:
R e ( t ) = ρ l   V r i s e ( t )   d e q ( t ) μ l ,
E o ( t ) = ( ρ l ρ g )   g   d e q 2 ( t ) σ ,
W e ( t ) = ρ l   V r i s e 2 ( t )   d e q ( t ) σ ,
M o = g   ( ρ l ρ g )   μ l 4 ρ l 2   σ 3 ,
where g is the acceleration of gravity. Note that the dimensionless numbers listed above are all instantaneous, except for the Morton number, which was constant and equal to 2.18 × 10−11. Mean values of the equivalent diameter, the aspect ratio, the rise velocity, and the dimensionless numbers listed above were computed by averaging across all instantaneous values of each individual bubble.
The present image processing methodology was developed specifically for the present experimental apparatus. In particular, as shown in Figure 3, the bubbles sometimes slide in front of the rod. When this happens, back lighting is not sufficient to resolve the bubble in its entirety and front lighting is also required, as schematically illustrated in Figure 1a. When front lighting is used, the bubble casts a shadow onto the rod, as can be noticed in Figure 3a–c,e). Clearly, the presence of a moving shadow makes the image background variable. Consequently, this makes the standard image processing approach, which is based on background subtraction, image binarization, and edge detection, unfeasible because subtracting a variable background is problematic. As is evident in Figure 3, the present image processing methodology is functional despite the variable background. Note that painting the rod to mitigate the bubble shadowing was not feasible because of the pool boiling experiments mentioned in Section 2.1, which required a metallic rod with clean surface. Moreover, the present image processing methodology can easily be applied to multiple bubbles which, as explained in Section 2.1, are within the scope of the present experimental setup, whereas resolving multiple bubbles with the standard image processing approach based on background subtraction, image binarization, and edge detection is not straightforward. Within the limits of this study, the present image processing methodology was found to be robust and accurate. In comparison with the standard image processing approach, which is largely computer-based, the present method is clearly more time-consuming. The present image processing methodology is similar to that used for bubble columns by Besagni and Inzoli [28,29], who identified bubbles by manually locating six points along the bubble border, and then calculated an approximating ellipse by directly interpolating through these six manually located points. Though the starting point is clearly the same, the present approach differs from that of Besagni and Inzoli [28,29] in the number of points used to resolve the border of the bubble, and in the way the bubble is reconstructed and its shape characterized.

3. Numerical Model

3.1. Governing Equations and Discretization Methods

The rise of air bubbles in water was simulated using the open-source package OpenFOAM v.1812 and the Volume of Fluid (VoF) method [30] implemented in the two-phase solver interIsoFoam [31]. According to VoF, liquid and gas are treated as a single-mixture fluid and the volume fraction α identifies the fraction of cell volume occupied by a selected primary phase, so that 0 α 1 . interIsoFoam solves the unsteady volume fraction, continuity and momentum equations for an isothermal, incompressible two-phase flow and Newtonian fluid in the following form:
α t + · ( α u ) = 0 ,
· u = 0 ,
( ρ u ) t + · ( ρ u u ) = p + · [ μ ( u + u T ) ] + ρ g + σ κ α ,
where t is time, u is the fluid velocity vector, ρ and μ are the mixture fluid density and dynamic viscosity, respectively, p is the pressure, and σ κ α introduces the surface tension force estimated via the Continuum Surface Force method [32], with σ being the surface tension coefficient and κ the local interface curvature, here calculated as κ = · ( α / | α | ) .
OpenFOAM discretizes the transport equations above with a finite volume method on a collocated grid arrangement, where all variables are defined at the control volume centres. interIsoFoam is a geometric VoF solver which discretizes Equation (8) according to a two-step procedure. First, an interface reconstruction step finds an approximation of the liquid–gas interface within each cell cut by the interface (where 0 < α < 1 ), by appropriate isosurface calculations. Then, an interface advection step calculates the volume of fluid crossing each control volume boundary and constituting the convective term of Equation (8), under the assumption that the interface translates steadily across the control volume face during the time interval. Details of the algorithm are provided by Roenby et al. [31]. Unlike OpenFOAM’s algebraic VoF solver interFoam, interIsoFoam guarantees a sharp interface representation, without the need of any artificial interface compression strategy. A first-order, bounded, implicit scheme (Euler) is used for the temporal discretization of the flow equations. Second-order schemes are adopted for all spatial derivatives: Gauss limitedLinearV and Gauss vanLeer for the convective terms in the momentum and volume fraction equations, respectively; Gauss linear corrected for all Laplacian schemes and surface normal gradients. OpenFOAM’s PIMPLE algorithm (combination of SIMPLE and PISO) is used for pressure–velocity coupling, with 3 PISO correctors (nCorrectors 3), no momentum prediction (momentumPredictor no), and 2 non-orthogonal correctors (nNonOrthogonalCorrectors 2) to account for the utilized non-orthogonal mesh. The residuals thresholds for the iterative solution of the flow equations are set to 10−7 for the velocity and 10−8 for both volume fraction and pressure. The time step of the simulation is variable and is calculated based on a maximum allowed Courant number of 0.5.
The present study covers values of the bubble Reynolds number, defined in Equation (4), on the order of R e 10 3 , and thus turbulence is likely to be present in the wake of the bubbles. RANS and LES turbulence models utilize empirical constants and wall functions calibrated with single-phase flow data, and therefore their applicability to interface-resolved simulations of partially confined bubbles, rising in an otherwise stagnant liquid, is not guaranteed. As such, the numerical results presented in this work were obtained by solving Equations (8)–(10) without any turbulence model, as done in previous studies with similar Reynolds numbers [33,34,35,36]. The spatial and temporal resolution of the simulation, chosen upon a grid independence analysis whose results are outlined in Section 3.2 below, set the smallest scales of the vortices that can be fully resolved by the numerical model.

3.2. Geometry and Boundary Conditions

The geometry simulated corresponds to the annular channel used for the experiments and was a vertical box of external dimensions 25.3 × 25.3 × 350   mm 3 , with a coaxial cylindrical rod of diameter of 9.99 mm subtracted from it; a sketch of the entire simulated domain is provided in Figure 4a.
In order to describe the flow, we adopt a Cartesian reference frame where z denotes the vertical coordinate, with z = 0 being the bottom surface of the box, while x and y indicate horizontal axes, with x = y = 0 coinciding with the rod axis, and being aligned with the external edges of the box cross-section (see Figure 4b,c). All the domain boundaries except for the top boundary are set as walls, with a no-slip condition for the velocity, a zero-gradient condition for the pressure, and a zero-contact-angle condition for the volume fraction, which models a hydrophilic wall to prevent bubble adhesion. The top wall is set as an outlet boundary, with a constant pressure value and zero-gradient conditions for velocity and volume fraction. The gravitational force is introduced as g = g z ,   with g = 9.81   m / s 2 . The fluids simulated are air and water at 295 K, with properties as shown in Table 1. At time t = 0 , the computational domain is filled with stagnant water and an initially spherical air bubble is patched nearby the channel bottom, centered at z = 0.007   mm and x = y = 0.0075   mm , i.e., along the diagonal of the channel cross-section and about halfway between the rod and external box surface; see illustration in Figure 4b.
Adaptive mesh refinement was utilized in order to provide sufficient resolution to the flow, while maintaining a coarse mesh far from the bubble. The background mesh is composed of cubic hexahedra (see Figure 4b,c) with two recursive levels of refinement nearby the rod surface, where the control volumes are clipped to fit the cylindrical surface. The mesh is dynamically refined at the bubble interface during runtime up to the maximum level of refinements selected, and the same criterion is also applied to the initial mesh at time t = 0 , as can be observed in Figure 4b. The requirements in terms of mesh resolution for the rise of bubbles in stagnant liquid is usually expressed in cells per bubble diameter. Hua et al. [33], Dijkhuizen et al. [34], and Roghair et al. [35] performed direct numerical simulations of air bubbles rising in an infinite pool of stagnant water for Reynolds numbers up to 103 [34,35] and 104 [33], employing 20 cells per bubble diameter. Loisy et al. [37], Balcázar et al. [38], and Esmaeli and Tryggvason [39] studied the rise of bubbles in the spherical and ellipsoidal regime for Reynolds numbers up to 100 with 25–30 cells per bubble diameter. Cano-Lozano et al. [40], Tripathi et al. [41], and Gumulya et al. [42] analyzed bubble shapes and trajectories for Reynolds numbers up to 100 employing adaptive mesh refinement with resolutions of 128 [40], 82 [41] and up to 220 [42] cells per bubble diameter at the liquid–gas interface. Gumulya et al. [36] simulated the rise of air bubbles in stagnant water for Reynolds numbers up to 103 using adaptive mesh refinement with up to 220 cells/diameter at the bubble interface.
In the present work, numerical simulations are run for a range of bubble diameters D b = 3 6   mm , and Reynolds numbers on the order of 103 are expected. We performed a mesh convergence analysis for a reference D b = 4.5   mm case, employing 5 cells/diameter in the bulk flow and 20 (2 recursive refinements), 40 (3), 80 (4) cells/diameter at the bubble interface, and another arrangement with 10 cells/diameter in the bulk flow and 40 (2 refinements) at the interface. The mesh is updated at the end of every time step. The results are reported in Figure 5. For all the configurations studied, the bubble trajectory (see Figure 5a) is approximately zigzag planar during the first 20 D b of the rise, and then develops into a helix of diameter of about 5 mm and pitch 30–35 mm. Throughout the rise, the bubble remains in the quarter of the channel where it was first generated. The velocity of the bubble centroid, calculated as V c = | d x c | / d t with x c = ( x c , y c , z c ) being the centroid position, is reported in Figure 5b. In all cases, the bubble speed oscillates between V c = 0.17   m / s and 0.25   m / s , with local minima occurring when the bubble approaches the walls. The time-average speed, calculated between z c = 0.2 0.3   m , is in the range V c = 0.216 0.22   m / s for all the meshes tested, indicating less than 2% differences. The vertical component of the rise velocity, estimated as V r i s e = d z c / d t , shows time averages in the range V r i s e = 0.199 0.205   m / s for all the meshes tested, which compare well with the experimental value. As such, all the meshes tested exhibited similar results and the configuration with 10–40 cells per bubble diameter in the bulk bubble was adopted for the simulations presented in this work. This corresponded to a total of about 5 million mesh cells for the reference case with D b = 4.5   mm .
A time step dependency analysis was carried out by testing a Courant number of 0.1: no appreciable differences were observed from the 0.5 Courant number case. Note that for a Courant number of 0.5, the simulation time step oscillated around Δ t = 3.5 × 10 5   s (and about Δ t = 0.7 × 10 5   s for a Courant of 0.1), which is well within the capillary time step limit, Δ t = [ ( ρ l + ρ g ) h 3 / ( 4 π σ ) ] 1 / 2 4 × 10 5   s , with h = D b / 40 being the mesh size across the liquid–gas interface.
Further preliminary tests (not documented here) were conducted by starting the simulation with g = 0 , and increasing g gradually as the time elapsed (at a rate of 1   m / s 2 every 0.1   s ) until its actual value was achieved, in order to emulate a gradual initial rise of the bubble with a smaller initial deformation. However, once the nominal value of g was restored, the bubble shape and speed showed no differences with the dynamics obtained by setting g = 9.81   m / s 2 from t = 0 . This showed that the impulsive bubble start adopted in the simulations did not affect the subsequent bubble rise dynamics. If this was the case, then the actual growth and detachment of the bubble would need to be simulated, adding considerable complexity to the numerical model.

3.3. Postprocessing of Numerical Data

From the numerical simulations, bubble interface data are extracted during runtime with regular temporal frequency, typically every 0.005   s , ensuring about 300 frames for each run. Figure 6a shows a snapshot of the bubble during the rise for the D b = 4.5   mm case, run with the coarsest (5–20) mesh. The points representing the liquid–gas interface are identified as the α = 0.5 iso-surface. This list of points is read in Matlab (version R2018b) and the built-in function alphaShape is utilized to create a bounding volume enveloping these points, enabling the calculation of geometrical quantities such as bubble surface area and volume. Further topological queries such as the extraction of centers, nodes and normal vectors of the boundary facets can be addressed with the function triangulation, as illustrated in Figure 6b. In order to characterize the geometry of the bubble, this is modelled using two-dimensional ellipses and three-dimensional ellipsoids. From the list of interface points, the bubble projections on the xz- and yz-planes are first obtained. Then, for each projection the boundary nodes are extracted and utilized to fit a two-dimensional ellipse [43] as shown in Figure 6c,d. From the fit, the ellipse axes lengths and aspect ratio are calculated. A different function is used to fit a three-dimensional ellipsoid [44] to the interface points and calculate axes lengths and aspect ratios. This procedure is repeated for all the saved bubble interface data, to finely resolve the bubble dynamics over time.

3.4. Validation of the Numerical Model

The numerical model was first validated by comparison with the experimental data of Bhaga and Weber [45] for air bubbles rising in stagnant, unconfined, aqueous sugar solutions. Eight different sets of conditions were selected from Bhaga and Weber [45] (see Figure 7) covering values of the Reynolds number from close to unity to the largest value achieved in their work, R e = 259 . The simulations were run in a three-dimensional domain of size 12 D b × 12 D b × 32 D b , with slip boundary conditions applied to all boundaries. Adaptive mesh refinement was employed, with 10 cells/diameter in the bulk flow. For the cases illustrated in Figure 7a–c,e 40 (2 recursive refinements) cells/diameter were employed at the bubble interface, whereas, for the cases illustrated in Figure 7d,f–h, 80 (3 recursive refinements) cells/diameter were used at the bubble interface to better capture the thin bubble skirt. The Courant number of the simulations was limited to 0.025. The bubble was initialized at the domain center, near the bottom wall. The terminal bubble shapes for simulations and experiments are illustrated in Figure 7. The bubble rose following a vertical rectilinear path at all conditions. Six of these cases (Figure 7a–d,f,g) are characterized by very similar values of the Eötvös number, E o = 114 116 , and Morton numbers ranging from M o = 848 to 8 × 10 4 , while three of these (Figure 7e,g,h) have similar Morton numbers, M o 8 × 10 4 , and an Eötvös number ranging from E o = 32 to 237. As the Morton number is decreased or the Eötvös number is increased, the bubble rises faster and shape transitions from a dimpled ellipsoidal-cap to a flattened spherical cap, eventually exhibiting an open wake and an asymmetric shape as the Reynolds number approaches 150. There is qualitatively a good agreement between the numerical and experimental bubble shapes reported in Figure 7.
A quantitative comparison of terminal bubble speed and aspect ratio is offered in Figure 7i–l. The bubble speed predicted by the simulation for R e 100 is always within 5% of the experimental data, whereas there is a systematic tendency to underestimate the experimental data for larger Reynolds numbers, although the maximum deviation remains below 10%. The same tendency of numerical simulations to underestimate the bubble velocities measured by Bhaga and Weber at larger Reynolds numbers was previously reported by the numerical studies of Hua and Lou [46] and Gumulya et al. [42], who used different simulations techniques and observed up to 10% deviations from the experimental data. The aspect ratio of the bubbles depicted in Figure 7l was calculated in the numerical simulations by taking the ratio between the height and width of the projection of the bubble profile on the xz and yz planes, thus disregarding the indentation at the bubble bottom. As the Reynolds number increases, the bubble flattens and the aspect ratio decreases from E 0.7 when R e 2 to E 0.2 when R e > 100 . The agreement between numerical and experimental data is excellent, with deviations below 5%, which is less than the uncertainty of the experimentally measured data for ellipsoidal bubbles indicated by Bhaga and Weber [45].
The second validation test that was performed consists of the simulation of an air bubble rising in an unconfined pool of stagnant water, for a range of bubble diameters D b = 3 5.2   mm , which covers the range analyzed in confined conditions in the Results section. The experimental data for unconfined air–water systems of Tomiyama et al. [17] and Veldhuis [47] are utilized as benchmark for the numerical simulations. Tomiyama et al. [17] studied trajectories, shapes and velocities of air bubbles rising in stagnant water for a range of bubble diameters D b = 0.5 5.5   mm , and reported very different bubble dynamics, depending on the initial shape deformation of the bubble at the instant of release from the injection nozzle. Bubbles with small initial shape deformation rose with a zigzag motion, lower speed and a larger aspect ratio; as opposed to bubbles released with a larger initial shape deformation, which rose with a helical motion, larger speed and lower aspect ratio. Importantly, bubbles released with larger initial deformation exhibited a significant scattering on their terminal speed and shape, as shown in Figure 8. Veldhuis [47] studied the behavior of air bubbles rising in water at similar conditions to Tomiyama et al. [17], with bubbles released from a capillary tube. They observed different bubble trajectories in the range D b = 3 5.2   mm , from zigzag to spiral and, eventually, chaotic as the bubble size increased. Their data for bubble rise velocities and aspect ratio sit at the boundary of the dataset from Tomiyama et al. [17], see Figure 8, with faster and more flattened bubbles. This can be ascribed to the large initial shape deformation, which might be induced by the injection capillary. Figure 8 also includes the rise velocity prediction extracted from an empirical correlation proposed by Park et al. [48], which applies to air bubbles rising in water in the range D b = 0.1 20   mm ; these predictions fall in between Tomiyama et al. [17] and Veldhuis [47] measurements. Numerical simulations were run in a box of size 28 D b × 28 D b × 150 D b . The fluid properties of air and water at 293 K were considered. A spherical bubble was initialized at the domain center, near the bottom wall. The domain was meshed using an adaptive mesh with 40 cells per bubble diameter at the air–water interface and 10 cells/diameter far from it. Meshes with higher refinement at the interface and/or in the bulk flow did not yield appreciable differences in the results. The Courant number of the simulations was set to 0.1. In the range of bubble diameters simulated, the bubbles rose with paths in between zigzag and helical. The terminal bubble rise speed and aspect ratio were calculated as time averages of the instantaneous values extracted for z c > 30 D b , when the rolling mean of the bubble speed became constant. The results from the simulations are presented in Figure 8. Bubble rise speed and aspect ratio sit within the area of the graphs occupied by the dataset of Tomiyama et al. [17] for helically rising bubbles. Although mild, the numerical data confirm the experimentally observed trends that both bubble rise speed and aspect ratio increase as the bubbles size is reduced and the effect of surface tension forces becomes more significant.
In summary, validation simulations of air bubbles rising in a liquid pool were conducted, covering a range of Reynolds numbers R e 2 1100 . Bubble rise velocities and aspect ratios agreed well with the selected experimental benchmark data, thus confirming the validity of the numerical model utilized in this work.

4. Results and Discussion

During the experiments we recorded and analyzed 30 individual bubbles, all generated under seemingly identical operating conditions, corresponding to a total of 597 individual frames, which should be sufficient to provide a reliable bubble shape distribution (approximately 300 or more frames are needed to provide reliable bubble size/shape distributions, as discussed by Besagni and Inzoli [28]). To help interpret the measurements and generalize the experimental observations to bubble diameters in the range of 3–6 mm, we numerically simulated the rise of six individual bubbles with diameters of 3 mm, 3.5 mm, 4 mm, 4.5 mm, 5.2 mm, and 6 mm. The mean values of the main bubbles’ parameters from the present experiments (averaged over all of the 30 bubbles recorded) and numerical simulations are summarized in Table 2.
As explained later in Section 4.1, the bubble mean equivalent diameter deduced from the two-dimensional bubble projection in the experiments is 5–6% smaller than the actual mean equivalent diameter deduced from the bubble volume in the numerical simulations. The mean bubble diameter of 4.3 mm measured in the experiments, therefore, corresponds to an actual bubble diameter of about 4.5 mm, so that the experimental figures summarized in Table 2 can be compared directly with the CFD simulations for the bubble of 4.5 mm. When comparing the present measurements with the numerical simulations, it can be noted that the simulated bubble rises about 10% faster than the measured bubble, and its aspect ratio is about 30% lower. Though not significant, the difference is larger than the present measuring errors, particularly for the aspect ratio (measuring errors are 6–7% for the rise velocity and 9–10% for the aspect ratio, as noted previously). The discrepancy between the present experiments and numerical simulations can be traced back to two main reasons: (1) turbulence in the wake of the bubbles that is not captured in the numerical simulations and (2) a possible contamination in the test water, as explained below.
At the Reynolds numbers covered in the present simulations (678–1194), turbulence is likely to be present in the wake of the bubbles. We performed preliminary tests by using a k-ω SST turbulence model [49]; however, the bubble exhibited a perfectly planar zig-zag trajectory throughout its rise, with an even larger speed. Therefore, the results presented in this work were obtained without any turbulence model, which means that turbulent eddies with spatial and temporal scales smaller than the mesh and time step size of the simulation were not captured. Normally, in external flows, the emergence of turbulence is accompanied by an increase in the drag in comparison with the laminar flow. Therefore, as a consequence of not accounting for turbulence, the numerical model can be expected to somewhat underpredict the drag, and thus overpredict the rise velocity. A higher rise velocity would increase the pressure force acting on the bubble which, in turn, would flatten, thereby yielding a lower aspect ratio.
Aside from turbulence, water contamination can be a cofactor responsible for the mismatch between the present experiments and numerical simulations. As is well known, fully eliminating contamination and surfactants from water is particularly difficult, if at all possible [3]. Therefore, despite the fact that we employed nominally micro-filtered and high-purity water, a minor contamination in the test water cannot be excluded. Clearly, the water in the numerical simulations is perfectly pure. Since bubbles rising through clean liquids rise faster and are more deformed than the corresponding bubbles rising in contaminated systems [1], a minor contamination in the test water could explain the lower rise velocity and the less pronounced bubble deformation observed in the present experiments, in comparison with the numerical simulations.
Despite the limitations in the present experiments and numerical simulations, the present results are nonetheless informative for the single bubble dynamics in the confined annular channel considered here.

4.1. Preliminary CFD Results

Preliminary CFD results that support the design of the present test system and of the associated measuring procedure are presented below.
As explained previously, in the experiments the bubbles were generated by blowing air through a capillary pipe at a frequency of one bubble generated every 2.2 s (bubbling frequency of 0.45 Hz). The CFD results provided in Figure 9 confirm that the delay of 2.2 s is sufficient for the liquid to recover between successive bubble passages. In particular, Figure 9a depicts a control volume of coordinates x > 0 ,   y > 0 ,   0.0115   m < z < 0.016   m , corresponding to a 4.5 mm long slice of a quarter of the annular channel located immediately above the initial position of the bubble, whereas the time variation of the specific kinetic energy of the water contained in this control volume is presented in Figure 9b,c. It is evident that, as soon as the bubble enters the control volume, the kinetic energy of the water increases as a consequence of the agitation induced by the rising bubble. Once the bubble leaves the control volume, the kinetic energy begins to decrease, indicating that viscous dissipation is gradually damping the water motion. By the time the bubble reaches the top of the channel at time 1.5 s, the kinetic energy of the water inside the control volume decreased by about four orders of magnitude, and, correspondingly, the mean velocity (Figure 9d,e) decreased by about two orders of magnitude, as compared to the peak values reached when the bubble rose through the control volume. Therefore, by the time the bubble reaches the top of the channel, the motion of the water in the control volume is practically extinguished. This confirms that the air flow rate set-point used in the experiments ensures continuous bubbling at a steady rate with a delay between successive bubbles sufficient to allow liquid recovery, and thus treat each bubble as a single bubble rising through a practically stagnant liquid.
Predicted bubble rise velocities are presented versus vertical elevation in Figure 10a. As can be noted, the initial acceleration transient is rather fast and ends after a rise of a few centimeters. This confirms that, by the time they reach the vertical elevation where the camera is positioned, the bubbles have reached their terminal dynamics. This latter state is clearly not a steady-state because the predicted rise velocity is fluctuating, as can be noted in Figure 10a and further discussed in the following Section 4.2.
As explained previously, only one camera was used in the experiments and the bubbles were characterized using their planar projections as seen in individual frames. The error associated with using only one camera in the present case is on the order of a few percent. This is illustrated in the CFD results presented in Figure 10b which provides, for all simulated bubbles, the ratio of the mean equivalent diameter deduced from the two-dimensional bubble projection (the diameter of the circle with the same area as the projection of the bubble) to the three-dimensional mean equivalent diameter (the diameter of the sphere with the same volume as the bubble). As noted previously, in the present case the variation of the air density due to the hydrostatic pressure variation along the channel is negligible. Consequently, the volume of the bubbles and the three-dimensional mean equivalent diameter are constants. As can be noted in Figure 10b, the mean equivalent diameter deduced from the two-dimensional bubble projection is 5–6% smaller than the three-dimensional mean equivalent diameter, an error that is rather small and within the present experimental resolution (as previously mentioned, the measuring error for the equivalent diameter is on the order of 9–10%).

4.2. Instantaneous Bubble Dynamics

The experimental instantaneous dynamic of one representative bubble is documented in Figure 11, where the sequence of the bubble projection perimeters extracted from consecutive images is provided together with the corresponding instantaneous variations of the bubble equivalent diameter, rise velocity, aspect ratio and inclination of the equivalent ellipse, as well as Reynolds, Eötvös and Weber number values. The bubble rise trajectory indicates a zigzag motion of small amplitude, and the instantaneous variations of the bubble parameters clearly indicate that the bubble rise dynamic is wobbling, as all bubble parameters show quite pronounced variations. This is also confirmed by the CFD results documented in Figure 10a, which clearly show a non-steady-state terminal dynamic.
The experimental instantaneous dynamic of all the measured bubbles is documented in Figure 12, where the rise trajectories of the centroid of the bubble projections extracted from consecutive images are provided together with the histograms of the corresponding instantaneous variations of the bubble equivalent diameter, rise velocity, aspect ratio and inclination of the equivalent ellipse, and Reynolds, Eötvös and Weber number values. Even though the instantaneous dynamic of each single bubble is wobbling (see Figure 11), the various trajectories appear reasonably similar, suggesting a fairly repeatable bubble rise dynamic. In turn, this suggests that the initial deformation of the bubbles did not vary appreciably during the experiments.
Air bubbles of a few mm in size rising in stagnant unconfined water at ambient conditions have an ellipsoidal shape and exhibit a wobbling rise dynamic [1,50]. The present results show a wobbling rise dynamic, and therefore indicate that the confinement of the annular channel does not change the qualitative character of the bubble dynamic.

4.3. Bubbles Mean Shape

The mean shape results of the bubbles are presented in Figure 13, where the mean aspect ratio from the present measurements and CFD simulations is presented as a function of the mean Eötvös number (panel a), mean Weber number (panel b), mean Tadaki number T a (panel c), and mean Reynolds number (panel d). The Tadaki number is defined as follows:
T a = R e   M o 0.23 ,
where R e and M o are the mean Reynolds and Morton numbers, respectively. The non-circular confined channel data provided by Tomiyama et al. [15] are also included in Figure 13, whereas the non-circular confined channel data provided by Venkateswararao et al. [14] are not included since the authors measured the terminal rise velocity but did not measure the shape of the bubbles.
To better put the non-circular confined channel results into perspective, the unconfined rise data by Tomiyama et al. [17] are also included in Figure 13: these refer to air bubbles at ambient conditions rising through unconfined water, and are therefore informative of the shapes of air bubbles rising through water when confinement effects are not present. Note that all data included in Figure 13 refer to air bubbles rising through clean water at ambient conditions, so that the main difference between individual data subsets is the presence/absence of confinement effects. Moreover, being only a function of the thermo-physical properties of the fluids (see Equation (7)), the Morton number has the same constant value for all data in Figure 13, so that the Tadaki number in the present case is essentially the Reynolds number multiplied by a constant. This is why the present aspect ratio data in Figure 13 show similar clustering when plotted versus the Tadaki number (panel c) and the Reynolds number (panel d).
The aspect ratio prediction methods proposed by Moore [51], Equations (12) and (13), and by Fan and Tsuchiya [52] (as quoted in [20]), for bubbles rising through unconfined clean liquids are also included for comparison in Figure 13 (panes b and c, respectively):
E = ( 1 + 0.1406   W e 0.0089   W e 2 + 0.0287   W e 3 ) 1 ,
E = { 0.77 + 0.24 tanh [ 1.9 ( 0.4 l o g 10 T a ) ] } 2       0.3 < T a < 20 .
In its original formulation, the prediction method proposed by Moore [51] is an implicit relation for the aspect ratio as a function of Weber number. The relation in Equation (12) used here is a convenient explicit approximation developed by Loth [3] for moderate deformations ( E 0.5 ). Finally, the present data are compared, in the parity plot in Figure 14, with the predictions of the aspect ratio prediction method proposed by Loth [3]:
E = 1 ( 1 E m i n )   t a n h ( c E W e )       0.2 < R e < 5000 ,
E m i n = 0.25 + 0.55   e x p ( 0.09   R e )
c E = 0.165 + 0.55   e x p ( 0.3   R e ) .
As can be noted in Figure 13, the Weber number appears to be the dimensionless group that is most effective at clustering the present data. The Tadaki number is also effective, though not as effective as the Weber number, whereas the Eötvös number does not appear to be effective at clustering the present data. Similar findings concerning the clustering effectiveness of the Weber, Tadaki and Eötvös numbers were reported by Liu et al. [20] for air bubbles rising through unconfined water and glycerol aqueous solution. The clustering effectiveness of the Weber number suggests that the shape of the present bubbles is largely controlled by the competition between the fluid-dynamic force acting on the bubble surface, which causes deformation, and the surface tension force, which resists deformation. Viscous dissipation in the liquid also plays a role, as evident from the relatively successful clustering effectiveness of the Tadaki and Reynolds numbers. The usefulness of the Weber number for predicting the mean aspect ratio is also indirectly confirmed in Figure 14, where it can be noted that the Loth [3] prediction method fits the present data fairly well. Even though this prediction method formally depends on both the Weber number and the Reynolds number, the effect of the latter is restricted to low and moderate values of the Reynolds number: indicatively below about 100. Above this limit, which happens to be the case for the present data (as can be noted in Figure 13d), the expressions in Equations (15) and (16) reached their asymptotic limiting values and, correspondingly, the Loth [3] prediction method essentially depends only on the Weber number. As can be noted in Figure 14, this prediction method seems more effective at predicting the confined rise data, as opposed to the unconfined data, part of which are overpredicted. The available data are however too restricted in scope to draw any definite conclusions.
By comparing the aspect ratios in Figure 13 for the data generated with non-circular confined channels (present experiment and Tomiyama et al. [15]) to the unconfined rise data (Tomiyama et al. [17]) and prediction methods it seems clear that, for a given Weber, Tadaki or Reynolds number value, the aspect ratio of the bubbles rising though a confined channel is higher than the aspect ratio of the corresponding bubbles rising through unconfined water. This indicates that the effect of the confinement is to reduce the deformation of the bubbles and, correspondingly, increase the mean aspect ratio.
As can be noted in Table 2, all dimensionless numbers for the simulated bubbles scale in direct proportion to the size of the bubble. It follows that the CFD data points in Figure 13 are oriented left-to-right in proportion to the size of the bubble: the leftmost point refers to the smallest bubble (3 mm), and the rightmost point refers to the largest bubble (6 mm). The present CFD results in Figure 13b,c show a change in trend with increasing Weber or Tadaki number. For small bubble size (3–4 mm), the aspect ratio decreases with the increasing Weber or Tadaki number. This agrees with the trend for an unconfined rise: in fact, the prediction methods in Equations (12) and (13) indicate that, when the bubbles are rising through unconfined liquids, the aspect ratio decreases as the Weber or Tadaki number increase. On the other hand, for larger bubbles (4.5 mm and above) the aspect ratio levels off and then increases with the increasing Weber or Tadaki number, which is the opposite of what happens during unconfined rise. This suggests that, in the present non-circular annular channel, the effect of the confinement on the shape of the bubble scales with the bubble size, which is the same qualitative trend observed with bubbles rising through confined circular tubes (Clift et al. [1], Krishna et al. [6]).

4.4. Bubble Mean Rise Velocity

The mean bubble rise velocities from the present measurements and CFD simulations are presented as a function of the mean bubble diameter in Figure 15a, together with the data for non-circular confined channels provided by Venkateswararao et al. [14] and Tomiyama et al. [15].
To better put the non-circular confined channel results into perspective, the unconfined rise data of Tomiyama at al. [17] are also included in Figure 15a, together with the prediction method developed by Mendelson [53] for bubbles rising through unconfined liquids:
V r i s e = 2   σ ρ l   d e q + 0.5   g   d e q ;
The data are also presented in Figure 15b as the mean drag coefficient versus mean Reynolds number. Following common practice, the mean drag coefficient C D and mean rise velocity were linked through a steady-state force balance between drag and buoyancy:
C D   1 2 ρ l V r i s e 2   π 4 d e q 2 = ( ρ l ρ g )   g   π 6 d e q 3 ;
From Equation (18), the mean drag coefficient can be computed from the mean rise velocity, and vice versa:
C D = 4   ( ρ l ρ g )   g   d e q 3   ρ l V r i s e 2 ,       V r i s e = 4   ( ρ l ρ g )   g   d e q 3   ρ l   C D .
The prediction method proposed by Mei et al. [54] for spherical bubbles rising through clean unconfined liquids, Equation (20), and the prediction method proposed by Loth [3] for spherical-cap bubbles rising through unconfined liquids, Equation (21), are also included in Figure 15b; these provide a lower bound and an upper bound, respectively, for bubbles rising through unconfined liquids:
C D = 16 R e   { 1 + [   8 R e + 1 2 ( 1 + 3.315 R e ) ] 1 } ,
C D = 8 3 + 16 R e ;
For further reference, the drag prediction method proposed by Brown and Lawler [55] for solid spheres is also included in Figure 15b:
C D = 24 R e   ( 1 + 0.150   R e 0.681 ) + 0.407 1 + 8.710 R e ;
Note that all data included in Figure 15 refer to air bubbles rising through clean water at ambient conditions, so that the main difference between individual data subsets is the presence/absence of confinement effects.
The present measurements and CFD simulations and the data for non-circular confined channels provided by Venkateswararao et al. [14] and Tomiyama et al. [15] are compared, in Figure 16, with the modified Mendelson prediction method, Equation (23), which is widely used for bubbles rising through confined circular tubes:
V r i s e = 2   σ ρ l   d e q + 0.5   g   d e q   [ 1 ( d e q D ) 2 ] 3 / 2 ,
where D is the diameter of the confining circular tube. On the right-hand side of Equation (23), the square root term corresponds to the Mendelson [53] prediction method in Equation (17) for bubbles rising through unconfined liquids, whilst the multiplying term is a correction proposed by Clift et al. [1] that accounts for the confinement. In Figure 16, the velocity from measurements or CFD simulations is divided by the unconfined rise velocity predicted from Equation (17), and then the velocity ratio is plotted as a function of the diameter ratio d e q / D . For the present measurements and CFD simulations the hydraulic diameter was used in place of the tube diameter D , whilst for the data of Venkateswararao et al. [14] and of Tomiyama et al. [15] the inner subchannel hydraulic diameter was used (as carried out by Tomiyama et al. [15]). The solid line in Figure 16 corresponds to the above confinement correction by Clift et al. [1]. The data by Krishna et al. [6], which refer to air bubbles at ambient conditions rising through confined water in circular tubes, are also included in Figure 16 for comparison with the non-circular channel data.
Finally, the present measurements and CFD simulations, together with the data of Venkateswararao et al. (1982) [14] and Tomiyama et al. (2003) [15], are compared in the parity plots in Figure 17, with the prediction methods proposed by Tomiyama et al. [56], Equation (24); Tomiyama et al. [17], Equation (25); and Rodrigue [57], Equations (26)–(27):
C D = m a x { m i n [   16 R e ( 1 + 0.15   R e 0.687 ) , 48 R e ] ,   8 3 E o E o + 4 } ,
V r i s e = a s i n 1 E 2 E 1 E 2 1 E 2 8   σ ρ l   d e q E 4 / 3 + ( ρ l ρ g )   g   d e q 2   ρ l E 2 / 3 1 E 2 ,
V r i s e = F 12   ( ρ l 2   d e q 2 μ l   σ ) 1 / 3 ( 1 + 1.31 × 10 5   M o 11 / 20   F 73 / 33 ) 21 / 176 ( 1 + 0.020   F 10 / 11 ) 10 / 11 ,
F = g ( ρ l 5   d e q 8 μ l 4   σ ) 1 / 3 .
The rise velocity corresponding to the drag coefficient from Equation (24) and the drag coefficients corresponding to the rise velocities from Equations (25) and (26) were computed according to Equation (19). The data provided by Venkateswararao et al. [14] are not included in Figure 17c,d because the authors did not measure the shape of the bubbles (the bubble aspect ratio is needed as an input to Equation (25)).
As can be noted in Figure 15a, for small bubble size (below about 3–3.5 mm) the rise velocity data generated with non-circular confined channels (present CFD simulations, Venkateswararao et al. [14], and Tomiyama et al. [15]) agree with each other and do not seem to depart significantly from the unconfined rise data (Tomiyama et al. [17]) included for comparison, albeit these latter data are more widely scattered: possibly the consequence of a more pronounced variation of the bubbles’ initial deformation during the tests. For larger bubble sizes (above about 3.5 mm), the present experiments and CFD simulations, as well as the data of Tomiyama et al. [15], are in fair agreement with each other and indicate that the rise velocity gradually decreases as the bubble size increases, with values of the rise velocity that appear consistently lower than the unconfined rise data and prediction method. On the other hand, the data of Venkateswararao et al. [14] indicate an opposite trend: the rise velocity gradually increases as the bubble size increases, with values of the rise velocity that appear consistently higher than the unconfined rise data and prediction method.
For the present CFD simulations and for the data of Tomiyama et al. [15], therefore, the effect of the confinement is to reduce the deformation of the bubbles, which increases the mean aspect ratio and reduces the mean rise velocity, which is the same qualitative trend observed with bubbles rising through confined circular tubes [1,6]. On the other hand, according to the data of Venkateswararao et al. [14], the effect of the confinement is an increase in the mean rise velocity. Unfortunately, the authors did not measure the shape of the bubbles, and thus it is not possible to provide a complete assessment of their data. However, it is clear that the cross-section of the test section used by the authors was considerably larger than those used here and used by Tomiyama et al. [15], and this could be the reason for the different trend observed. In turn, this would indicate that, with non-circular channels of a complex shape, the size of the channel cross-section may affect the dynamics of the bubbles, and results generated with comparatively small channels may not extrapolate to channel sizes of industrial relevance.
As noted by Loth [3], Equations (20) and (21) provide a lower bound and an upper bound to the drag of bubbles rising through unconfined liquids, respectively. In other words, drag data fall in the area comprised between the two curves. As can be noted in Figure 15b, this is also the case for most of the available confined rise data. Only a few data points by Tomiyama et al. [15] peak above Equation (21), indicating that the reduction in rise velocity due to confinement can yield a substantial increase in the associated drag coefficient.
As can be noted in Figure 16, the agreement between the present data and CFD simulations, and the data of Tomiyama et al. [15] and Krishna et al. [6] with the modified Mendelson prediction method, is fair but the scatter is substantial. This indicates that the confinement effect may also depend on other factors beyond the diameter ratio that essentially capture the geometry of the system, possibly the Reynolds number, as is the case for the solid spheres falling through confined liquids [55]. However, the available data are too restricted in scope to draw any definite conclusions. It is evident in Figure 16 that the trend of the measurements by Venkateswararao et al. [14] is remarkably different from that of the rest of the data. Nonetheless, it is not possible to rule these data out as outliers: note that the present CFD simulation results for the small bubble size (3–3.5 mm) agree with these data. Moreover, as already noted the cross-section of the test section used by Venkateswararao et al. [14] is considerably larger than those used by the others, and this could explain the different trend observed. It may also be the case that the subchannel hydraulic diameter, which is presently used as a representative length scale, does not fully capture the effect of the channel cross-sectional shape on the dynamics of the bubbles. What seems clear is that the available data are too restricted in scope to duly assess confinement effects with channels of a complex shape, and more investigations are clearly required.
From the parity plots in Figure 17, it can be noted that the prediction method of Tomiyama et al. [17] is quite successful at fitting both the unconfined and the confined rise data, whereas the other prediction methods considered here are ineffective. As can be noted from inspecting Equations (24)–(27), the prediction method of Tomiyama et al. [17] is the only one that explicitly incorporates the shape of the bubbles via the aspect ratio. This confirms that the size, shape, and the rise velocity of ellipsoidal bubbles are closely linked together, and should all be incorporated into prediction methods.

5. Concluding Remarks

Using experiments and numerical simulations, we studied the shape and the rise velocity of single air bubbles, measuring 3–6 mm in diameter, rising through water and confined inside a non-circular annular channel. Our main findings are summarized below:
  • The confinement of the present annular channel did not affect the qualitative behavior of the bubbles, which exhibited a wobbling rise dynamic similar to that observed with bubbles rising through unconfined liquids. The effect of the confinement was evident on the shape and rise velocity; the bubbles were less deformed and rose more slowly in comparison with bubbles rising through unconfined liquids;
  • The present results are in fair agreement with previous observations by Tomiyama et al. [15] on the shape and rise velocity of air bubbles rising though water, confined in a small cross-section subchannel, and with available data on confined rise through circular tubes. However, the observations by Venkateswararao et al. [14] on air bubbles rising through a larger cross-section tubular test section show an increase in rise velocity as consequence of the confinement, which is a trend that was not observed in other studies on confined rise through smaller channels. This indicates that confinement effects with non-circular channels of complex shape could be more complicated than those observed with circular tubes, and both the size of the channel cross-section and its shape may affect the dynamics of the bubbles; therefore, results generated with comparatively small channels may not extrapolate to channel sizes of industrial relevance. The available data are, however, too restricted in scope to draw any definite conclusions, and more investigations are clearly needed;
  • The present data and numerical simulations, as well as the other data collected from the literature and used here, indicate that the size, shape, and rise velocity of ellipsoidal bubbles are closely linked together, and this should be considered when designing prediction methods;
  • The image processing methodology developed and used here, based on manually digitizing points along the bubble border, is robust and effective at dealing with the variable image background caused by the bubble shadow. This technique can be extended to multiple bubbles of interest in bubble columns;
  • The synergetic use of experiments and numerical simulations proved to be an effective approach for the study of single bubble rise in confined geometries. In particular, we found the numerical simulations instrumental in providing a better insight into the measurements, and effective for generalizing the experimental observations, thereby compensating for the limitations of the test setup. Still, we found the numerical simulations somewhat hampered by the current limitations in RANS turbulence models for bubbly flows, and further validation work in this respect would be particularly beneficial.

Author Contributions

Conceptualization, A.C. and M.M.; experiments, A.C.; simulations, M.M.; writing—original draft preparation, A.C. and M.M.; writing—review and editing, A.C. and M.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data discussed are directly presented in the paper and are therefore directly accessible.

Acknowledgments

M. Magnini acknowledges the use of Athena at HPC Midlands+, which was funded by the EPSRC grant EP/P020232/1, as part of the HPC Midlands+ consortium. M. Magnini acknowledges the support from EPSRC, through the BONSAI (EP/T033398/1) grant.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Clift, R.; Grace, J.; Weber, M. Bubbles, Drops, and Particles; Dover Publications: New York, NY, USA, 2005; ISBN 0-486-44580-1. [Google Scholar]
  2. Kulkarni, A.A.; Joshi, J.B. Bubble Formation and Bubble Rise Velocity in Gas-Liquid Systems: A Review. Ind. Eng. Chem. Res. 2005, 44, 5873–5931. [Google Scholar] [CrossRef]
  3. Loth, E. Quasi-Steady Shape and Drag of Deformable Bubbles and Drops. Int. J. Multiph. Flow 2008, 34, 523–546. [Google Scholar] [CrossRef]
  4. Uno, S.; Kintner, R.C. Effect of Wall Proximity on the Rate of Rise of Single Air Bubbles in a Quiescent Liquid. AIChE J. 1956, 2, 420–425. [Google Scholar] [CrossRef]
  5. Salami, E.; Vignes, A.; Le Goff, P. Hydrodynamique Des Dispersions. II. Effet de Paroi. Mouvement d’une Goutte Ou d’une Bulle Dans Un Fluide Immobile Contenu Dans Un Tube Vertical de Petit Diamètre. Genie Chim. 1965, 94, 67–77. [Google Scholar]
  6. Krishna, R.; Urseanu, M.I.; van Baten, J.M.; Ellenberger, J. Wall Effects on the Rise of Single Gas Bubbles in Liquids. Int. Commun. Heat Mass Transf. 1999, 26, 781–790. [Google Scholar] [CrossRef]
  7. Collins, R. A Simple Model of the Plane Gas Bubble in a Finite Liquid. J. Fluid Mech. 1965, 22, 763. [Google Scholar] [CrossRef]
  8. Bush, J.W.M.; Eames, I. Fluid Displacement by High Reynolds Number Bubble Motion in a Thin Gap. Int. J. Multiph. Flow 1998, 24, 411–430. [Google Scholar] [CrossRef]
  9. Roig, V.; Roudet, M.; Risso, F.; Billet, A.-M. Dynamics of a High-Reynolds-Number Bubble Rising within a Thin Gap. J. Fluid Mech. 2012, 707, 444–466. [Google Scholar] [CrossRef] [Green Version]
  10. Wang, X.; Klaasen, B.; Degrève, J.; Blanpain, B.; Verhaeghe, F. Experimental and Numerical Study of Buoyancy-Driven Single Bubble Dynamics in a Vertical Hele-Shaw Cell. Phys. Fluids 2014, 26, 123303. [Google Scholar] [CrossRef] [Green Version]
  11. Filella, A.; Ern, P.; Roig, V. Oscillatory Motion and Wake of a Bubble Rising in a Thin-Gap Cell. J. Fluid Mech. 2015, 778, 60–88. [Google Scholar] [CrossRef] [Green Version]
  12. Hashida, M.; Hayashi, K.; Tomiyama, A. Rise Velocities of Single Bubbles in a Narrow Channel between Parallel Flat Plates. Int. J. Multiph. Flow 2019, 111, 285–293. [Google Scholar] [CrossRef]
  13. Hashida, M.; Hayashi, K.; Tomiyama, A. Effects of Fine Particles on Terminal Velocities of Single Bubbles in a Narrow Channel between Parallel Flat Plates. Int. J. Multiph. Flow 2020, 127, 103270. [Google Scholar] [CrossRef]
  14. Venkateswararao, P.; Semiat, R.; Dukler, A.E. Flow Pattern Transition for Gas-Liquid Flow in a Vertical Rod Bundle. Int. J. Multiph. Flow 1982, 8, 509–524. [Google Scholar] [CrossRef]
  15. Tomiyama, A.; Nakahara, Y.; Adachi, Y.; Hosokawa, S. Shapes and Rising Velocities of Single Bubbles Rising through an Inner Subchannel. J. Nucl. Sci. Tech. 2003, 40, 136–142. [Google Scholar] [CrossRef]
  16. Wu, M.; Gharib, M. Experimental Studies on the Shape and Path of Small Air Bubbles Rising in Clean Water. Phys. Fluids 2002, 14, L49. [Google Scholar] [CrossRef] [Green Version]
  17. Tomiyama, A.; Celata, G.P.; Hosokawa, S.; Yoshida, S. Terminal Velocity of Single Bubbles in Surface Tension Force Dominant Regime. Int. J. Multiph. Flow 2002, 28, 1497–1519. [Google Scholar] [CrossRef]
  18. Okawa, T.; Tanaka, T.; Kataoka, I.; Mori, M. Temperature Effect on Single Bubble Rise Characteristics in Stagnant Distilled Water. Int. J. Heat Mass Transf. 2003, 46, 903–913. [Google Scholar] [CrossRef]
  19. Peters, F.; Els, C. An Experimental Study on Slow and Fast Bubbles in Tap Water. Chem. Eng. Sci. 2012, 82, 194–199. [Google Scholar] [CrossRef]
  20. Liu, L.; Yan, H.; Zhao, G. Experimental Studies on the Shape and Motion of Air Bubbles in Viscous Liquids. Exp. Therm. Fluid Sci. 2015, 62, 109–121. [Google Scholar] [CrossRef]
  21. Lemmon, E.W.; Huber, M.L.; McLinden, M.O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport. Properties-REFPROP, Version 9.1; Natl Std. Ref. Data Series (NIST NSRDS); NIST-National Institute of Standards and Technology: Gaithersburg, MD, USA, 2013.
  22. Liu, L.; Yan, H.; Zhao, G.; Zhuang, J. Experimental Studies on the Terminal Velocity of Air Bubbles in Water and Glycerol Aqueous Solution. Exp. Therm. Fluid Sci. 2016, 78, 254–265. [Google Scholar] [CrossRef]
  23. Aoyama, S.; Hayashi, K.; Hosokawa, S.; Tomiyama, A. Shapes of Ellipsoidal Bubbles in Infinite Stagnant Liquids. Int. J. Multiph. Flow 2016, 79, 23–30. [Google Scholar] [CrossRef]
  24. Rohatgi, A. WebPlotDigitizer. 2021. Available online: https://automeris.io/WebPlotDigitizer/ (accessed on 6 July 2021).
  25. Eaton, J.W.; Bateman, D.; Hauberg, S.; Wehbring, R. GNU Octave 4.2 Reference Manual; Samurai Media Limited: Hong Kong, 2017. [Google Scholar]
  26. Solomon, C.; Breckon, T. Fundamentals of Digital Image Processing: A Practical Approach with Examples in Matlab; Wiley-Blackwell: Chichester, UK; Hoboken, NJ, USA, 2011; ISBN 978-0-470-84472-4. [Google Scholar]
  27. Russ, J.C.; Neal, F.B. The Image Processing Handbook, 7th ed.; CRC Press Taylor & Francis Group: Boca Raton, FL, USA; New York, NY, USA; London, UK, 2016; ISBN 978-1-138-74749-4. [Google Scholar]
  28. Besagni, G.; Inzoli, F. Bubble Size Distributions and Shapes in Annular Gap Bubble Column. Exp. Therm. Fluid Sci. 2016, 74, 27–48. [Google Scholar] [CrossRef] [Green Version]
  29. Besagni, G.; Inzoli, F. Comprehensive Experimental Investigation of Counter-Current Bubble Column Hydrodynamics: Holdup, Flow Regime Transition, Bubble Size Distributions and Local Flow Properties. Chem. Eng. Sci. 2016, 146, 259–290. [Google Scholar] [CrossRef] [Green Version]
  30. Hirt, C.W.; Nichols, B.D. Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39, 201–225. [Google Scholar] [CrossRef]
  31. Roenby, J.; Bredmose, H.; Jasak, H. A Computational Method for Sharp Interface Advection. R. Soc. Open Sci. 2016, 3, 160405. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Brackbill, J.U.; Kothe, D.B.; Zemach, C. A Continuum Method for Modeling Surface Tension. J. Comput. Phys. 1992, 100, 335–354. [Google Scholar] [CrossRef]
  33. Hua, J.; Stene, J.F.; Lin, P. Numerical Simulation of 3D Bubbles Rising in Viscous Liquids Using a Front Tracking Method. J. Comput. Phys. 2008, 227, 3358–3382. [Google Scholar] [CrossRef] [Green Version]
  34. Dijkhuizen, W.; Roghair, I.; Van Sint Annaland, M.; Kuipers, J.A.M. DNS of Gas Bubbles Behaviour Using an Improved 3D Front Tracking Model—Drag Force on Isolated Bubbles and Comparison with Experiments. Chem. Eng. Sci. 2010, 65, 1415–1426. [Google Scholar] [CrossRef]
  35. Roghair, I.; Mercado, J.M.; Sint Annaland, M.V.; Kuipers, H.; Sun, C.; Lohse, D. Energy Spectra and Bubble Velocity Distributions in Pseudo-Turbulence: Numerical Simulations vs. Experiments. Int. J. Multiph. Flow 2011, 37, 1093–1098. [Google Scholar] [CrossRef]
  36. Gumulya, M.; Joshi, J.B.; Utikar, R.P.; Evans, G.M.; Pareek, V. Characteristics of Energy Production and Dissipation around a Bubble Rising in Water. Chem. Eng. Sci. 2019, 193, 38–52. [Google Scholar] [CrossRef]
  37. Loisy, A.; Naso, A.; Spelt, P.D.M. Buoyancy-Driven Bubbly Flows: Ordered and Free Rise at Small and Intermediate Volume Fraction. J. Fluid Mech. 2017, 816, 94–141. [Google Scholar] [CrossRef] [Green Version]
  38. Balcázar, N.; Lehmkuhl, O.; Jofre, L.; Oliva, A. Level-Set Simulations of Buoyancy-Driven Motion of Single and Multiple Bubbles. Int. J. Heat Fluid Flow 2015, 56, 91–107. [Google Scholar] [CrossRef] [Green Version]
  39. Esmaeeli, A.; Tryggvason, G. A Direct Numerical Simulation Study of the Buoyant Rise of Bubbles at O(100) Reynolds Number. Phys. Fluids 2005, 17, 093303. [Google Scholar] [CrossRef]
  40. Cano-Lozano, J.C.; Martínez-Bazán, C.; Magnaudet, J.; Tchoufag, J. Paths and Wakes of Deformable Nearly Spheroidal Rising Bubbles Close to the Transition to Path Instability. Phys. Rev. Fluids 2016, 1, 053604. [Google Scholar] [CrossRef] [Green Version]
  41. Tripathi, M.K.; Sahu, K.C.; Govindarajan, R. Dynamics of an Initially Spherical Bubble Rising in Quiescent Liquid. Nat. Commun. 2015, 6, 6268. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Gumulya, M.; Joshi, J.B.; Utikar, R.P.; Evans, G.M.; Pareek, V. Bubbles in Viscous Liquids: Time Dependent Behaviour and Wake Characteristics. Chem. Eng. Sci. 2016, 144, 298–309. [Google Scholar] [CrossRef] [Green Version]
  43. Brown, R. Fitellipse.m. 2021. Available online: https://www.mathworks.com/matlabcentral/fileexchange/15125-fitellipse-m (accessed on 6 July 2021).
  44. Yury. Ellipsoid Fit. 2021. Available online: https://www.mathworks.com/matlabcentral/fileexchange/24693-ellipsoid-fit (accessed on 6 July 2021).
  45. Bhaga, D.; Weber, M.E. Bubbles in Viscous Liquids: Shapes, Wakes and Velocities. J. Fluid Mech. 1981, 105, 61. [Google Scholar] [CrossRef] [Green Version]
  46. Hua, J.; Lou, J. Numerical Simulation of Bubble Rising in Viscous Liquid. J. Comput. Phys. 2007, 222, 769–795. [Google Scholar] [CrossRef]
  47. Veldhuis, C.H.J. Leonardo’s Paradox: Path and Shape Instabilities of Particles and Bubbles. Ph.D. Thesis, University of Twente, Enschede, The Netherlands, 2006. [Google Scholar]
  48. Park, S.H.; Park, C.; Lee, J.; Lee, B. A Simple Parameterization for the Rising Velocity of Bubbles in a Liquid Pool. Nucl. Eng. Technol. 2017, 49, 692–699. [Google Scholar] [CrossRef]
  49. Menter, F.R.; Kuntz, M.; Langtry, R. Ten Years of Industrial Experience with the SST Turbulence Model. In Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass Transfer, Antalya, Turkey, 12–17 October 2003; Begell-House: Antalya, Turkey, 2003. [Google Scholar]
  50. Lunde, K.; Perkins, R.J. Shape Oscillations of Rising Bubbles. Appl. Sci. Res. 1997, 58, 387–408. [Google Scholar] [CrossRef]
  51. Moore, D.W. The Velocity of Rise of Distorted Gas Bubbles in a Liquid of Small Viscosity. J. Fluid Mech. 1965, 23, 749–766. [Google Scholar] [CrossRef]
  52. Fan, L.S.; Tsuchiya, K. Dynamics in Liquids and Liquid-Solid Suspensions; Butterworth-Heinemann: Oxford, UK, 1990. [Google Scholar]
  53. Mendelson, H.D. The Prediction of Bubble Terminal Velocities from Wave Theory. AIChE J. 1967, 13, 250–253. [Google Scholar] [CrossRef]
  54. Mei, R.; Klausner, J.F.; Lawrence, C.J. A Note on the History Force on a Spherical Bubble at Finite Reynolds Number. Phys. Fluids 1994, 6, 418–420. [Google Scholar] [CrossRef]
  55. Brown, P.P.; Lawler, D.F. Sphere Drag and Settling Velocity Revisited. J. Environ. Eng. 2003, 129, 222–231. [Google Scholar] [CrossRef]
  56. Tomiyama, A.; Kataoka, I.; Zun, I.; Sakaguchi, T. Drag Coefficients of Single Bubbles under Normal and Micro Gravity Conditions. JSME Int. J. Ser. B 1998, 41, 472–479. [Google Scholar] [CrossRef]
  57. Rodrigue, D. A General Correlation for the Rise Velocity of Single Gas Bubbles. Can. J. Chem. Eng. 2008, 82, 382–386. [Google Scholar] [CrossRef]
Figure 1. Schematic representation of the experimental setup (vertical dimension not to scale): (a) top view (as explained in Section 2.1, out of the four injection lines present only one was used whilst three were inactive), (b) side view (only the active injection line is shown), (c) detail of the nozzle tip used to generate the bubbles.
Figure 1. Schematic representation of the experimental setup (vertical dimension not to scale): (a) top view (as explained in Section 2.1, out of the four injection lines present only one was used whilst three were inactive), (b) side view (only the active injection line is shown), (c) detail of the nozzle tip used to generate the bubbles.
Fluids 06 00437 g001
Figure 2. Image processing methodology adopted here: (a) raw RGB image (the region delimited with the broken line is highlighted in panels b-e), (b) raw RGB image detail with discrete points (in red color) manually digitized along the bubble border, (c) raw RGB image detail with increased density points (in red color) generated with cubic spline interpolation along the bubble border, (d) binary image of the bubble generated from the polygon representing the bubble border, (e) raw RGB image detail with equivalent ellipse superimposed (in green color).
Figure 2. Image processing methodology adopted here: (a) raw RGB image (the region delimited with the broken line is highlighted in panels b-e), (b) raw RGB image detail with discrete points (in red color) manually digitized along the bubble border, (c) raw RGB image detail with increased density points (in red color) generated with cubic spline interpolation along the bubble border, (d) binary image of the bubble generated from the polygon representing the bubble border, (e) raw RGB image detail with equivalent ellipse superimposed (in green color).
Fluids 06 00437 g002
Figure 3. Representative example of image processing for a bubble sliding in front of the rod: (a) raw RGB image (the region delimited with the broken line is highlighted in panels (be)), (b) raw RGB image detail with discrete points (in red color) manually digitized along the bubble border, (c) raw RGB image detail with increased density points (in red color) generated with cubic spline interpolation along the bubble border, (d) binary image of the bubble generated from the polygon representing the bubble border, (e) raw RGB image detail with equivalent ellipse superimposed (in green color).
Figure 3. Representative example of image processing for a bubble sliding in front of the rod: (a) raw RGB image (the region delimited with the broken line is highlighted in panels (be)), (b) raw RGB image detail with discrete points (in red color) manually digitized along the bubble border, (c) raw RGB image detail with increased density points (in red color) generated with cubic spline interpolation along the bubble border, (d) binary image of the bubble generated from the polygon representing the bubble border, (e) raw RGB image detail with equivalent ellipse superimposed (in green color).
Fluids 06 00437 g003
Figure 4. Sketch of the computational domain: (a) Image of the whole geometry, in slight transparency, with the initial spherical bubble represented in light blue; (b) close-up view of bubble and mesh on selected planes at t = 0; (c) two-dimensional view of the mesh on the channel horizontal cross-section. This setup corresponds to the case with bubble diameter D b = 4.5   mm and mesh with 40 cells per bubble diameter at the liquid–gas interface and 10 in the bulk liquid, which is the mesh used throughout this study.
Figure 4. Sketch of the computational domain: (a) Image of the whole geometry, in slight transparency, with the initial spherical bubble represented in light blue; (b) close-up view of bubble and mesh on selected planes at t = 0; (c) two-dimensional view of the mesh on the channel horizontal cross-section. This setup corresponds to the case with bubble diameter D b = 4.5   mm and mesh with 40 cells per bubble diameter at the liquid–gas interface and 10 in the bulk liquid, which is the mesh used throughout this study.
Fluids 06 00437 g004
Figure 5. Results of grid independence analysis; D b = 4.5   mm . (a) Position and (b) speed of the bubble centroid during the rise. The legend in (a,b) indicates the number of cells per bubble diameter in the bulk and at the bubble interface. The average speed of the bubble, calculated between z c = 0.2 0.3   m , is: 5–20) 0.216 m/s, 5–40) 0.218 m/s, 5–80) 0.22 m/s, 10–40) 0.217 m/s.
Figure 5. Results of grid independence analysis; D b = 4.5   mm . (a) Position and (b) speed of the bubble centroid during the rise. The legend in (a,b) indicates the number of cells per bubble diameter in the bulk and at the bubble interface. The average speed of the bubble, calculated between z c = 0.2 0.3   m , is: 5–20) 0.216 m/s, 5–40) 0.218 m/s, 5–80) 0.22 m/s, 10–40) 0.217 m/s.
Fluids 06 00437 g005
Figure 6. Postprocessing of numerical simulation data; D b = 4.5   mm . (a) Snapshot of the bubble surface (in blue) identified as α = 0.5 iso-surface, with the central rod (dark gray) and the computational mesh over a selected horizontal plane (in slight transparency); (b) Reconstruction of the bubble surface in Matlab; (c,d) Projected bubble on a (c) xz-plane and (d) yz-plane, with identified boundary nodes of the projection and two-dimensional fitted ellipse. For representation purposes, the results shown in this figure were obtained with a coarser mesh, with 20 and 5 cells per bubble diameter at the interface and bulk liquid, respectively.
Figure 6. Postprocessing of numerical simulation data; D b = 4.5   mm . (a) Snapshot of the bubble surface (in blue) identified as α = 0.5 iso-surface, with the central rod (dark gray) and the computational mesh over a selected horizontal plane (in slight transparency); (b) Reconstruction of the bubble surface in Matlab; (c,d) Projected bubble on a (c) xz-plane and (d) yz-plane, with identified boundary nodes of the projection and two-dimensional fitted ellipse. For representation purposes, the results shown in this figure were obtained with a coarser mesh, with 20 and 5 cells per bubble diameter at the interface and bulk liquid, respectively.
Fluids 06 00437 g006
Figure 7. Validation of the numerical model versus the experimental data of Bhaga and Weber [45]. (ah) Experimental and numerical terminal bubble shapes at different conditions, sorted by increasing values of the Reynolds number; (i) Numerical versus experimental bubble Reynolds number and (l) aspect ratio, with the red squares in (l) identifying simulation results for cases (ad,f), and blue circles identifying simulation results for cases (e,g,h).
Figure 7. Validation of the numerical model versus the experimental data of Bhaga and Weber [45]. (ah) Experimental and numerical terminal bubble shapes at different conditions, sorted by increasing values of the Reynolds number; (i) Numerical versus experimental bubble Reynolds number and (l) aspect ratio, with the red squares in (l) identifying simulation results for cases (ad,f), and blue circles identifying simulation results for cases (e,g,h).
Fluids 06 00437 g007aFluids 06 00437 g007b
Figure 8. Validation of the numerical model versus the experimental data of Tomiyama et al. [17] and Veldhuis [47]: (a) bubble rise velocity and (b) aspect ratio versus bubble diameter. The black dashed line in (a) indicates the bubble rise speed predicted using Park et al. [48] correlation.
Figure 8. Validation of the numerical model versus the experimental data of Tomiyama et al. [17] and Veldhuis [47]: (a) bubble rise velocity and (b) aspect ratio versus bubble diameter. The black dashed line in (a) indicates the bubble rise speed predicted using Park et al. [48] correlation.
Fluids 06 00437 g008
Figure 9. (a) View of the control volume (in red) used in the calculation of the kinetic energy of the water during bubble rise ( D b = 4.5   mm ), with the blue sphere depicting the bubble at t = 0 ; time variation of the kinetic energy of the water inside the control volume in linear (b) and logarithmic (c) scale; time variation of the mean velocity (computed as the square root of the kinetic energy) of the water inside the control volume in linear (d) and logarithmic (e) scale. The red vertical lines mark the times the bubble enters (leftmost red line) and leaves (rightmost red line) the control volume.
Figure 9. (a) View of the control volume (in red) used in the calculation of the kinetic energy of the water during bubble rise ( D b = 4.5   mm ), with the blue sphere depicting the bubble at t = 0 ; time variation of the kinetic energy of the water inside the control volume in linear (b) and logarithmic (c) scale; time variation of the mean velocity (computed as the square root of the kinetic energy) of the water inside the control volume in linear (d) and logarithmic (e) scale. The red vertical lines mark the times the bubble enters (leftmost red line) and leaves (rightmost red line) the control volume.
Fluids 06 00437 g009aFluids 06 00437 g009b
Figure 10. Predicted bubble rise velocities plotted versus vertical elevation (a), and ratio of the mean equivalent diameter deduced from the two-dimensional bubble projection to the three-dimensional mean equivalent diameter deduced from the bubble volume (b).
Figure 10. Predicted bubble rise velocities plotted versus vertical elevation (a), and ratio of the mean equivalent diameter deduced from the two-dimensional bubble projection to the three-dimensional mean equivalent diameter deduced from the bubble volume (b).
Fluids 06 00437 g010
Figure 11. Instantaneous dynamic of one representative measured bubble: (a) sequence of the bubble projection perimeters extracted from consecutive images (the time elapsed between successive frames is 10.36 ms) and corresponding instantaneous variations of the bubble equivalent diameter (b), rise velocity (c), aspect ratio (d) and inclination (e) of the equivalent ellipse, and Reynolds (f), Eötvös (g) and Weber (h) number values.
Figure 11. Instantaneous dynamic of one representative measured bubble: (a) sequence of the bubble projection perimeters extracted from consecutive images (the time elapsed between successive frames is 10.36 ms) and corresponding instantaneous variations of the bubble equivalent diameter (b), rise velocity (c), aspect ratio (d) and inclination (e) of the equivalent ellipse, and Reynolds (f), Eötvös (g) and Weber (h) number values.
Fluids 06 00437 g011
Figure 12. Instantaneous dynamic of all measured bubbles: (a) rise trajectories of the centroid of the bubble projections extracted from consecutive images (the vertical dashed lines in blue color denote the position of the rod), and histograms of the corresponding instantaneous variations of the bubble equivalent diameter (b), rise velocity (c), aspect ratio (d) and inclination (e) of the equivalent ellipse, and Reynolds (f), Eötvös (g) and Weber (h) number values.
Figure 12. Instantaneous dynamic of all measured bubbles: (a) rise trajectories of the centroid of the bubble projections extracted from consecutive images (the vertical dashed lines in blue color denote the position of the rod), and histograms of the corresponding instantaneous variations of the bubble equivalent diameter (b), rise velocity (c), aspect ratio (d) and inclination (e) of the equivalent ellipse, and Reynolds (f), Eötvös (g) and Weber (h) number values.
Fluids 06 00437 g012
Figure 13. Mean aspect ratio: measured and computed values versus Eötvös number (a), Weber number (b), Tadaki number (c), and Reynolds number (d) (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15] and this study: confined rise).
Figure 13. Mean aspect ratio: measured and computed values versus Eötvös number (a), Weber number (b), Tadaki number (c), and Reynolds number (d) (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15] and this study: confined rise).
Fluids 06 00437 g013
Figure 14. Mean aspect ratio: present data vs. predictions of the Loth [3] method (the dashed lines are ±20% bounds).
Figure 14. Mean aspect ratio: present data vs. predictions of the Loth [3] method (the dashed lines are ±20% bounds).
Fluids 06 00437 g014
Figure 15. Mean rise velocity versus bubble mean diameter (a) and mean drag coefficient versus mean Reynolds number (b) (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15], Venkateswararao et al. (1982) [14], and this study: confined rise).
Figure 15. Mean rise velocity versus bubble mean diameter (a) and mean drag coefficient versus mean Reynolds number (b) (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15], Venkateswararao et al. (1982) [14], and this study: confined rise).
Fluids 06 00437 g015
Figure 16. Velocity ratio versus diameter ratio.
Figure 16. Velocity ratio versus diameter ratio.
Fluids 06 00437 g016
Figure 17. Mean rise velocity and corresponding drag coefficient: present data vs. predictions of the Tomiyama et al. [56] method (a,b), of the Tomiyama et al. [17] method (c,d), and of the Rodrigue [57] method (e,f); the dashed lines are ± 20% bounds. (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15], Venkateswararao et al. (1982) [14] and this study: confined rise).
Figure 17. Mean rise velocity and corresponding drag coefficient: present data vs. predictions of the Tomiyama et al. [56] method (a,b), of the Tomiyama et al. [17] method (c,d), and of the Rodrigue [57] method (e,f); the dashed lines are ± 20% bounds. (Tomiyama et al. (2002) [17]: unconfined rise; Tomiyama et al. (2003) [15], Venkateswararao et al. (1982) [14] and this study: confined rise).
Fluids 06 00437 g017
Table 1. Relevant properties of water-air at 295 K.
Table 1. Relevant properties of water-air at 295 K.
ρ l (kg/m3) ρ g   ( kg / m 3 ) μ l   ( μ Pa   s ) μ g   ( μ Pa   s ) σ (mN/m)
9981.1895818.372.4
Table 2. Mean values of the main bubbles’ parameters.
Table 2. Mean values of the main bubbles’ parameters.
EXPCFD
d e q (mm) *4.33.03.54.04.55.26.0
V r i s e   (m/s)0.180.2160.2140.2100.2000.1990.190
E 0.810.6790.6260.6020.6050.6080.634
E o 2.491.221.662.172.753.674.88
W e 1.931.942.212.452.502.853.00
R e 80567878288094410831194
(*) In the experiments d e q was deduced from the area of the bubble planar projection, whilst in the numerical simulations it was deduced from the actual bubble volume.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Cioncolini, A.; Magnini, M. Shapes and Rise Velocities of Single Bubbles in a Confined Annular Channel: Experiments and Numerical Simulations. Fluids 2021, 6, 437. https://doi.org/10.3390/fluids6120437

AMA Style

Cioncolini A, Magnini M. Shapes and Rise Velocities of Single Bubbles in a Confined Annular Channel: Experiments and Numerical Simulations. Fluids. 2021; 6(12):437. https://doi.org/10.3390/fluids6120437

Chicago/Turabian Style

Cioncolini, Andrea, and Mirco Magnini. 2021. "Shapes and Rise Velocities of Single Bubbles in a Confined Annular Channel: Experiments and Numerical Simulations" Fluids 6, no. 12: 437. https://doi.org/10.3390/fluids6120437

APA Style

Cioncolini, A., & Magnini, M. (2021). Shapes and Rise Velocities of Single Bubbles in a Confined Annular Channel: Experiments and Numerical Simulations. Fluids, 6(12), 437. https://doi.org/10.3390/fluids6120437

Article Metrics

Back to TopTop