Bubble Formation and Motion in Liquids—A Review

Abstract

In flotation, a bubble acts as a carrier for attached particles. The properties of the gas–liquid interface of the bubble are one of the main factors determining the bubble motion and flotation efficiency. Monitoring of the bubble motion may deliver interesting information about the state of the gas–liquid interface. In the case of pure liquids, a bubble surface is fully mobile, while the presence of surface-active substances (e.g., surfactants) causes diminishing bubble velocity due to the retardation of the interface fluidity. The theoretical prediction of the terminal velocity value for the bubble has been investigated for over a century, delivering a number of various models describing bubble motion in a liquid. This narrative review is devoted to the motion of the bubble in stagnant liquids and is divided into three main sections describing: (i) experimental techniques for tracking bubble motion, (ii) bubble motion and shape deformation in clean water, and (iii) bubble motion in solutions of surface-active substances.

1. Introduction

Froth flotation is an industrial process, where gas bubbles act as carriers of hydrophobic substances (e.g., minerals, liquid drops) to the froth layer [1,2]. Air bubble–particle interactions during the collision and ‘collection’ of mineral particles by the rising bubbles are an essence of the elementary flotation act. The formation of a stable bubble–grain aggregate can be considered as a result of three sub–processes: (i) collision (encounter), (ii) attachment, and (iii) detachment. For process efficiency, all of those steps need to happen in a very short (millisecond) time frame [3].

Flotation efficiency can be affected by hydrodynamics conditions (i.e., shear forces, turbulence), the presence of substances affecting the stability of the liquid films formed, a degree of dispersion of the gas phase (size and shape of the bubbles), properties of the gas–liquid interface, etc. [3,4]. Despite the complexity of the process, it is possible to focus on the key events and elementary acts of flotation, such as: (i) the dynamics of the bubble collision with various interfaces (i.e., liquid–gas and liquid–solid); (ii) factors affecting the stability of the thin liquid film formed during the bubble collisions (dynamic conditions); (iii) the timescale of the bubble rupture and attachment to the interface; (iv) bubble velocity; and (v) the state of the dynamic adsorption layer over the bubble surface [5,6,7,8,9].

There are three main groups of reagents employed in flotation: (i) collectors, (ii) frothers, and (iii) modifiers (such as activators, depressants, dispersants, and pH regulators) [1,2,10]. Collectors and frothers are surface-active substances (SAS), which are responsible for the modification of the liquid–solid and liquid–gas interfaces, respectively. Collectors are added to increase the hydrophobicity of particles [1,2]. Frothers are responsible for the higher dispersion of the air phase, preventing bubble coalescence, increasing the probability of the formation of bubble–particle aggregates by lowering the rising velocity of bubbles [11,12,13]. Thus, bubble velocity might be of great importance for process efficiency.

In pure liquids, the velocity of rising bubbles depends on the viscosity and density of the continuous phase and the bubble diameter, while in solutions of surface-active substances the bubble rising velocity is lowered as a result of the formation of the motion-induced dynamic architecture of the adsorption layer (DAL) over the bubble surface, which retards the fluidity of the interface. Moreover, the kinetics of the thin liquid film rupture in SAS solutions of different concentrations are changed not only due to the lowering of the bubble impact velocity, but also as a result of (i) the retardation of the film thinning velocity and (ii) the increased stability of the films formed against external disturbances [14,15].

From the experimental point of view, tracking the bubble motion in liquids does not require sophisticated techniques. From the other side, the theoretical prediction of the terminal velocity value for a bubble is a rather complicated task and it has been investigated for over a century. The amount of work conducted on this topic can be realized when reading the literature reviews delivered by Clift et al. [16], Kulkarni and Joshi [17], and Chen et al. [18]. Following recent developments in modeling bubble motion, this paper presents a short review of analytical and semi-analytical models for predicting bubble terminal velocity in ultra-pure water and solutions of surface-active substances. Besides presenting various concepts derived for both ultra-pure and ‘contaminated’ (SAS) systems, two models are emphasized. The first one, was presented by Legendre et al. [19], together with the Moore theory [20,21], gives great agreement with experimental data for ultra-pure water. The other model was proposed by Kowalczuk et al. [22] and allows us to predict bubble terminal velocity in ‘contaminated’ systems for a whole SAS concentration range.

2. Experimental Methods to Study Bubble Motion

The investigation of bubble motion in liquids does not require the implementation of sophisticated experimental methods. Bubble velocity is defined as:

$$U=\frac{\sqrt{\left(x_{i+1}-x_i\right)+\left(y_{i+1}-y_i\right)+\left(z_{i+1}-z_i\right)}}{\mathsf{\Delta }t}$$

(1)

where (x_i+1, y_i+1, z_i+1) and (x_i, y_i, z_i) are the coordinates of the subsequent positions (in three axes) of the bubble and ∆t is the time interval; thus, visual observation is a commonly used technique.

The simplest method of measuring bubble velocity is following a bubble movement between two marks by an observer equipped with a stopwatch [23] (which can be replaced by photodiodes). Probably the only advantage of this method is its simplicity, while a number of disadvantages are presented. To mention few, (i) it is almost impossible to measure local velocity with sufficient accuracy, (ii) to determine terminal velocity it is necessary to place the first (starting) marker above the acceleration/deceleration stages, (iii) to increase measurement accuracy the bubble pathway needs to be elongated, (iv) no data on bubble shape are delivered.

The most common method used for studying bubble rise in liquids is visual observation with a camera. A camera equipped with a proper objective allows us to measure with high accuracy bubble local velocity and shape deformation. An experimental set-up needs to be equipped with components allowing for fast enough image acquisition. These issues can be solved in two ways: (i) a standard camera with acquisition of 10–30 fps and a stroboscopic lamp with a frequency of at least 100 Hz or (ii) a speed–camera with acquisition of >100 fps. Of course, the aforementioned parameters depend on the type of bubble–liquid system studied and should be individually adjusted. A schematic presentation of a typical set-up for tracking bubble rise is shown in Figure 1.

The main weakness of this method is that with one camera it is possible to track the bubble motion in only two axes. In ‘contaminated’ (SAS) liquids, the bubble often rises with either a zig-zag or a spiral path [25,26,27]; thus, an additional camera, placed 90° horizontally to the first camera, might be useful for improved tracking of bubble movement [28]. The smart budget solution allowing us to avoid using the additional camera is the integration of a vision system with a set of mirrors [29,30,31], as presented in Figure 2.

Tracking bubble movement with cameras is the most popular and reliable method. This method allows us to gain data on the bubble velocity, deformation, and path, giving comprehensive information about the air–liquid interface. The main disadvantage is that the collected images need to be processed with computer software to extract data on the bubble movement. The other limitation of the visual method is related to the visibility of the tracked object. Thus, it is necessary that liquid should show a transparency good enough for monitoring the bubble motion.

The determination of bubble velocity on-line and in a non-transparent medium is possible with an ultrasound technique, as presented by Borkowski and Zawala [32]. In this method, the ultrasonic 5 MHz transmitter and receiver are mounted on the bottom of the liquid column. The bubble rising velocity is determined by analyzing the variations in the temporal evolution of a position of the registered signal formed as a result of ultrasonic waves reflected from the rising bubble surface. The authors concluded that the ultrasonic method of determining the rising velocity of a single bubble is reliable, but not as accurate as a visual method. Observed deviation is probably related to the chosen speed of sound in a studied liquid, which is used during the on-line conversion of the signal into the bubble velocity. The main disadvantage is a lack of information on the bubble shape, which is significant for a better description of the state of the dynamic adsorption layer.

The formation of gas bubbles in liquids can be realized either as a result of nucleation in oversaturated liquids or, more commonly, by gas dispersion [33,34]. In dispersion methods, the bubbles are generated as a result of mixing liquid and gas phases with an energy input. Bubbles are commonly produced by sparging, that is, pumping gas through a capillary or a frit into the liquid. In this process, multi-body interactions between the bubbles cause the generated bubbles to be of various diameters. The addition of the surface-active substances can prevent the coalescence of bubbles, which leads to the low scatter of bubble diameter in dispersion. By increasing the SAS concentration, the degree of the bubble coalescence decreases, and at a particular concentration (critical coalescence concentration (CCC)), the coalescence of the bubbles is almost completely prevented [11,34,35,36]. However, in the case of slow bubble formation at a single capillary orifice, where no multi-body interactions are observed, the diameter of formed bubbles (d_b) can be well controlled and is described by Tate’s law [37]:

$$d_b={\left(\frac{6d_c\sigma cos\theta }{\mathsf{\Delta }\rho g}\right)}^{1/3}$$

(2)

where d_c is the inner diameter of the capillary, σ is the surface tension, Δρ is the density difference between liquid and gas, g is the gravitational acceleration, and $\theta$ is the contact angle between gas and capillary material. As the bubble motion causes the deformation of its shape, due to the resistance of the continuous (viscous) phase, the equivalent diameter (d_e) is often used in size descriptions of the rising bubbles. The bubble equivalent diameter is defined as a diameter of the sphere with the same volume as the rising bubble. The real shape of the rising bubble can be described with a good approximation as the oblate spheroid of horizontal (d_h) and vertical (d_v) diameters, which gives a formula for the equivalent diameter in the form:

$$d_e={\left(d_h^2{d}_v\right)}^{1/3}$$

(3)

while the bubble shape can be described using the bubble deformation degree (χ), defined as:

$$\chi =d_h/d_v$$

(4)

Figure 3 presents a comparison of the bubble diameter for experimental results and theoretical values calculated using Equation (2) as a function of surface-active substance concentration. For all the presented SAS (n-octanol, α-terpineol, and n-cetyltrimethylammonium bromide (CTAB)), a good agreement between the experimental and calculated diameters is observed, even for high solution concentrations.

From a practical point of view, the easiest way of delivering gas to the system is through connection to the gas tank, where the gas flow rate is regulated by controlling the pressure in the line with a pressure regulator. Another popular method is the use of a syringe pump, where the gas flow is regulated by a pumping speed. In both methods, bubbles are produced continuously at equal time intervals on the capillary orifice, which defines the bubble size, as described by Equation (1). In 1993, Oguz and Prosperetti [39] proposed the possibility of controlling bubble size by injecting air into flowing liquid. Nowadays, the proposed idea is realized with the use of microfluidic chips, giving the opportunity to form bubbles in a very wide range of diameters [40,41,42,43,44].

Vejrazka et al. [45] designed a generator for the “on-demand” production of bubbles. The device is equipped with a movable needle, through which air is injected. The moment of bubble detachment is controlled by a rapid needle movement. In order for the instrument to work properly, a few conditions need to be fulfilled: (i) the instantaneous flow rate should not exceed the critical value, (ii) the waiting period between bubbles should not be long compared to the bubbling period, (iii) the needle acceleration should be high enough. The last condition might have a substantial impact on bubble velocity just after detachment. The authors summarized that a relatively small change in the needle acceleration results in an important change in the bubble velocity, and, in the case of high needle acceleration, a bubble can be launched at a velocity which is higher than the terminal one.

Sanada and Abe [46] proposed the application of inducing acoustic pressure waves of different amplitudes into a deformable silicone tube with a cut slitting as a bubble generation orifice to generate bubbles. This technique allows bubbles to be produced with a wide distribution of diameters, from ca. 1 to 10 mm. However, the bubble generation process, i.e., its growth and detachment, at the slitting is rather turbulent, resulting in significant changes in the bubble shape, affecting the adsorption process on the growing bubble. Thus, the authors concluded that this method could be used for studying the motion of bubbles in ultra-pure water, while in solutions of SAS, such initial bubble area oscillations are disadvantageous for the kinetics of DAL formation.

Najafi et al. [47] used a pressure pulse technique to produce a single bubble. In this method, changing the pressure impulse (maximum bubble pressure), taper length, and inclination angle of the micropipette allows bubbles of diameters ranging between 400 and 1200 µm to be generated with very good reproducibility.

A slightly different approach to bubble formation using a ‘bubble-on-demand’ (BoD) generator is presented in Zawala and Niecikowska [48]. In this generator, the main emphasis is placed on the adjustment of bubble detachment frequency, which is controlled by an applied pressure pulse, while the bubble diameter depends on the inner diameter of the capillary used. The need to design such an instrument is related to the other part of the whole set-up, i.e., a bubble trap for controlling the SAS adsorption time over the bubble surface (Figure 4). The duration of the residue of the bubble in such a trap could be adjusted. After the desired adsorption time, the trap is opened, releasing the bubble saturated with the surfactant molecules to a desired degree. Using such an approach, the degree of initial bubble coverage by SAS is precisely controlled [49].

3. Bubble Motion in Pure Liquids

The motion (rising or falling) of a body in liquid with its terminal velocity (U_T) is described by the balance between the buoyancy (F_B) and the drag force (F_D):

$$F_B=\mathsf{\Delta }\rho V_{b}g$$

(5)

$$F_D=0.5C_D{\rho }_L{U}_T^2{S}_b$$

(6)

where C_D is the drag coefficient, ρ_L is the liquid density, and V_b and S_b are the volume and surface of projection on a horizontal plane of the body, respectively. Then, for a spherical body of the diameter d_b, the equation for terminal velocity can be expressed as:

$$U_T={\left(\frac{4d_b\mathsf{\Delta }\rho g}{3C_D{\rho }_L}\right)}^{1/2}$$

(7)

The main problem in the analytical solving of Equation (7) is the determination of the drag coefficient. The C_D depends on the conditions of motion, often expressed as a function of the dimensionless Reynolds number (Re) [16]:

$$Re=\frac{d_b{U}_T{\rho }_L}{{\eta }_L}$$

(8)

where η_L is the liquid viscosity. In fluid mechanics, the Reynolds number gives a measure of the ratio of inertial forces to viscous drag and is often used for the description of flow conditions. It is assumed that a bubble moves under creeping flow conditions when Re ≪ 1, a spherical regime for Re < 20, an ellipsoidal regime for 20 < Re < 4700, and a spherical or ellipsoidal cap regime for Re > 4700 [50].

The Reynolds number, together with the Morton (Mo) and Eovots (Eo) numbers are used for the characterization of the hydrodynamic conditions of the body motion in fluid. The Morton and Eovots (also called Bond) numbers present the shape of bubbles or drops moving in a surrounding fluid [16]:

$$Mo=\frac{{\eta }_L^{4}g\mathsf{\Delta }\rho }{{\rho }_L^2{\sigma }^3}$$

(9)

$$Eo=\frac{d_b^{2}g\mathsf{\Delta }\rho }{\sigma }$$

(10)

The other, often used, dimensionless numbers are the Weber (We), Archimedes (Ar), and Lyshchenko (Ly) numbers [16]:

$$We=\frac{d_b{U}_T^2\mathsf{\Delta }\rho }{\sigma }$$

(11)

$$Ar=\frac{d_b^{3}g\mathsf{\Delta }\rho {\rho }_L}{{\eta }_L^2}$$

(12)

$$Ly=\frac{U_T^3{\rho }_L^2}{g\mathsf{\Delta }\rho {\eta }_L}$$

(13)

In the case of a single bubble formed on a capillary orifice, two stages at the profile of the local velocity can be distinguished: (i) an acceleration stage, when the local velocity increases monotonically, and (ii) a terminal stage, when a constant value of the velocity is achieved. Figure 5 presents the local velocity profiles for bubbles of various diameters as a function of the distance from the capillary orifice. It can be observed that the bubble accelerates rapidly and, after ca. 30–40 mm (depending on the bubble size), the velocity of the bubble starts to be constant. The establishment of the terminal velocity is a consequence of the balance between the drag force and buoyancy.

One of the very first models of the C_D, known as Stokes’ law, was obtained by solving the Navier–Stokes equation for a symmetrical solid sphere falling in liquid under creeping flow conditions [52]:

$$C_D=\frac{24}{\mathit{Re}}$$

(14)

which, together with Equations (7) and (8), gives:

$$U_T=\frac{d_b^{2}g\mathsf{\Delta }\rho }{18{\eta }_L}$$

(15)

In the early 20th century, Hadamard [53] and Rybczynski [54] independently, showed that the terminal velocity of a fluid sphere is up to 50% higher than the velocity of a rigid sphere of the same size and density. This effect is related to the smaller viscous drag exerted by the liquid phase at the liquid–gas interface and the internal gas circulation induced inside a bubble, as presented schematically in Figure 6. Thus, the equation for terminal velocity takes the form [53,54]:

$$U_T=\frac{d_b^{2}g\mathsf{\Delta }\rho }{6{\eta }_L}\cdot \frac{{\eta }_L+{\eta }_G}{2{\eta }_L+3{\eta }_G}\approx \frac{d_b^{2}g\mathsf{\Delta }\rho }{12{\eta }_L}$$

(16)

where η_G is the gas viscosity.

It is interesting to mention that Pawliszak et al. [56] showed that the terminal velocities of bubbles of 215–550 μm in diameter deviate from the theoretical ones predicted for a fully mobile water—air interface, while for bubbles of below 180 μm in diameter good compliance with the Stokes’ model is observed. This shift from full mobility to an immobile water—air interface is associated with the sensitivity of tiny (d_b < 500 μm) bubbles to any surface-active impurities in water. Vakarelski et al. [57] confirmed, by carrying out similar experiments in perfluorocarbon liquids, that for the clean continuous phase the bubble is at the mobile liquid—air interface, giving higher terminal velocity than that predicted by the Stokes’ law, which follows Equation (16).

Levich [58] and Ackeret [59] independently elaborated a model for the potential flow of spherical drops and bubbles in clean water, where the expression for the C_D for Re < 50 and d_b < 0.5 mm is given as:

$$C_D=\frac{48}{\mathit{Re}}$$

(17)

Bubbles of large diameters have non-spherical shapes, which influences the drag coefficient. For such an oblate spheroid, Moore [20,21] obtained, considering the dissipation of energy on the boundary layer, the following relation of the C_D:

$$C_D=\frac{48}{Re}G\left(\chi \right)\left[1+\frac{H\left(\chi \right)}{\sqrt{Re}}+O\left(\frac{1}{\sqrt{Re}}\right)\right]$$

(18)

where G(χ) and H(χ) are functions of the bubble deformation degree (χ):

$$G\left(\chi \right)=\frac{1}{3}{\chi }^{4/3}{\left({\chi }^2-1\right)}^{3/2}\frac{\sqrt{{\chi }^2-1}-\left(2-{\chi }^2\right){\mathrm{s}}^{-1}\chi }{{\left({\chi }^2{\mathrm{s}}^{-1}\chi -\sqrt{{\chi }^2-1}\right)}^2}$$

(19)

$$H\left(\chi \right)=0.0195{\chi }^4-0.2134{\chi }^3+1.7026\chi -1.5732$$

(20)

Both functions, G(χ) and H(χ), were simplified for χ < 2 by Loth [50]:

$$G\left(\chi \right)=0.1287+0.4256\chi +0.4466{\chi }^2$$

(21)

$$H\left(\chi \right)=0.8886+0.5693\chi -0.4563{\chi }^2$$

(22)

The deformation can be found from the relationship between the Weber number and χ. For small distortions, i.e., χ < 1.1, this is [21]:

$$\chi =1+\frac{9}{64}We+O\left(W{e}^2\right)$$

(23)

while, for moderate and great deformations (χ > 1.1):

$$We=4{\chi }^{-4/3}\left({\chi }^3+\chi -2\right){\left({\chi }^2{\mathrm{s}}^{-1}\chi -\sqrt{{\chi }^2-1}\right)}^2{\left({\chi }^2-1\right)}^{-3}$$

(24)

For χ < 2, Equations (23) and (24) can be approximated by a simple relationship:

$$\chi =1+\frac{9}{64}We-0.0089W{e}^2+0.0287W{e}^3$$

(25)

The Moore model [20,21] is said to be applicable for a wide range of the Reynolds number, 100 < Re < 10,000, but a comparison of the prediction with experimental data [16,60] showed that this model underestimates values of the bubble terminal velocities for bubbles of d_b > 1.2 mm. Thus, Duineveld [60] proposed the following formula:

$$\chi ={\left(\frac{a_1}{W{e}_0-We}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$

(26)

where We₀ and a₁ are 4.412 and 4.39, respectively. The main limitation of Equation (26) is that for We → We₀, the χ is increasing indefinitely.

Based on data for ultra-pure water over a range of Reynolds numbers from 0.2 to 5000, Loth [50] proposed an approximate correlation for the aspect ratio of ‘uncontaminated’ bubbles as:

$$\chi ={\left[1-\left(1-E_{\mathit{min}}\right)\mathit{tanh}\left(c_{E}We\right)\right]}^{-1}$$

(27)

where:

$$E_{\mathit{min}}=0.250+0.55\mathit{exp}\left(-0.09\mathit{Re}\right)$$

(28)

$$c_E=0.165+0.55\mathit{exp}\left(-0.30\mathit{Re}\right)$$

(29)

Legendre et al. [19] proposed a simple formula closely describing the bubble deformation in water for χ < 3:

$$\chi ={\left(1-\frac{9}{64}We\right)}^{-1}$$

(30)

In general, the Weber number, which is defined as the ratio of continuous-fluid stresses (causing deformation) to surface tension stresses (resisting deformation), describes bubble geometry qualitatively by the following relationships [50]:

• We ≪ 1—bubble shape tends to a spherical geometry;
• We ~ 1—moderate deviations, oblate spheroid is observed;
• We ≫ 1—large bubble shape distortion occurs; spherical cap and oblate ellipsoidal cap are observed.

In their monograph, Clift et al. [16] reviewed most of the known models and compared them with experimental data. They emphasized the influence of bubble deformation on the drag coefficient for bubbles rising in clean water, which resulted in the presentation of two equations, for two separate regimes:

• for Re < 150

$$C_D=14.9R{e}^{-0.78}$$

(31)

• for Re > 565:

$$U_T={\left(\frac{2.14\sigma }{d_b{\rho }_L}+0.505d_{b}g\right)}^{1/2}$$

(32)

The other semi-analytical formula for the C_D, which fits with experimental data for Re < 130 in clean water, was proposed by Masliyah et al. [61]:

$$C_D=\frac{16}{Re}\left(1+0.077R{e}^{0.65}\right)$$

(33)

After analyzing the available experimental data and the correlations for the motion of rising gas bubbles, Karamanev [62,63] proposed a semi-analytical equation linking the bubble rising velocity and its geometry. Moreover, the drag coefficient was described in a more convenient way—in terms of the Archimedes number (Ar). Using these considerations, the terminal velocity of the bubble can be calculated as follows:

$$U_T=\frac{d_e}{d_h}\sqrt{\frac{4g{d}_b}{3C_D}}$$

(34)

and the drag coefficient is expressed as:

$$C_D=\frac{432}{Ar}\left(1+0.047A{r}^{2/3}\right)+\frac{0.517}{1+154A{r}^{-1/3}}$$

(35)

for Ar < 13,000, and C_D = 0.95 for Ar > 13,000. The term d_e/d_h is related to the changes in the bubble area projected on a horizontal plane, and can be calculated using the correlation proposed by Clift et al. [16]:

$$\frac{d_e}{d_h}={\left(1+0.163E{o}^{0.757}\right)}^{-1/3}$$

(36)

for Eo < 40, and d_e/d_h = 0.62 for Eo > 40. This model does not take into account the fluidity of the liquid–air interface, giving rather poor predictions when compared with experimental data, as seen in Figure 7.

Subsequently, Mei et al. [64,65] proposed a theoretical prediction of the terminal velocity based on the numerical solution of their model for the unsteady motions of particles (including droplets and bubbles) in a carrier fluid. For Re < 50, the drag coefficient is given as [65]:

$$C_D=\frac{16}{\mathit{Re}}\left\{1+{\left[\frac{8}{\mathit{Re}}+0.5\left(1+\frac{3.315}{{\mathit{Re}}^{1/2}}\right)\right]}^{-1}\right\}$$

(37)

Rodrigue [66,67] proposed another generalized correlation, which uses the flow number (Fl) and the velocity number (Ve) introduced earlier by Abou-El-Hassan [68]:

$$Fl=g{\left(\frac{{\rho }_L^5{d}_b^8}{\sigma {\eta }_L^4}\right)}^{1/3}$$

(38)

$$Ve=U_T{\left(\frac{{\rho }_L^2{d}_b^2}{\sigma {\eta }_L}\right)}^{1/3}$$

(39)

Based on the experimental data available, Rodrigue found the following relationship [69]:

$$Ve=\frac{Fl}{12}\left[\frac{{\left(1+1.31\cdot {10}^{-5}\times M{o}^{11/20}F{l}^{73/33}\right)}^{21/176}}{{\left(1+0.020F{l}^{10/11}\right)}^{10/11}}\right]$$

(40)

Figure 7. Terminal velocity of an air bubble rising in ultra-pure water as a function of the bubble diameter (points—experimental results [48,56], lines—theoretical models [16,19,20,21,50,55,60,61,63,64,65,66,67,70,71,72]).

Figure 7. Terminal velocity of an air bubble rising in ultra-pure water as a function of the bubble diameter (points—experimental results [48,56], lines—theoretical models [16,19,20,21,50,55,60,61,63,64,65,66,67,70,71,72]).

Figure 7 presents a comparison between experimental results and most of the described models within this review for the terminal velocity of bubbles rising in ultra-pure water. As can be seen, most of the models overestimate bubble velocities for d_b < 1.0 mm, while underestimating bubble velocities for d_b > 1.0 mm. A good agreement between the experimental data and the Moore model combined with the relationship between the Weber number and χ proposed by Legendre et al. [19] is observed.

4. Bubble Motion in the Solution of a Surface-Active Substance

The presence of surface-active substances (SAS) and their adsorption at the bubble surface affect (retard) the fluidity of the surface, which leads to a lowering of the bubble velocity [49,73,74,75,76,77,78,79,80]. The mechanism responsible for this effect was described for the first time by Frumkin and Levich [55,81]. They postulated that a bubble detaching from a capillary has uniform adsorption coverage of the SAS molecules over its surface but viscous drag exerted by the liquid on the rising bubble surface induces the uneven distribution of the SAS molecules. Depletion at the upstream part of the bubble and accumulation of the SAS molecules in the rear part of the bubble results in an inducement of the surface tension gradients and a tangential ‘Marangoni’ stress opposing the flow shear stress (Figure 8) [81]. Thus, the formation of this dynamic adsorption layer over the bubble interface is the reason for the retardation of the bubble surface fluidity and internal circulation, which leads to an increase in the drag force towards that of a rigid sphere [82,83,84,85,86].

As seen in Figure 9, in a solution of SAS, up to four bubble motion stages can be determined [49,74,76,79,80]: (i) acceleration, (ii) maximum velocity, (iii) deceleration, and (iv) terminal velocity. All physicochemical parameters (e.g., the value of the initial acceleration, the height and width of the velocity maximum, the period of bubble deceleration and the distance (time) needed for terminal velocity establishment, and bubble shape pulsations) describing the motion of the bubble for each mentioned stage can be used to track the DAL development and the moment of the formation of its full architecture, causing the immobilization of the rising bubble (liquid–gas) interface.

The bubble, detached from the capillary orifice, accelerates and attains its terminal velocity at a distance, which is highly dependent on the type and concentration of SAS (Figure 9). Moreover, the existence of the acceleration/deceleration stage for some concentrations of SAS, before the establishment of the terminal velocity, can be addressed by the formation of the DAL architecture during the initial stage of the bubble motion [24,75,76].

There is also a minimum degree of adsorption coverage (concentration), different for various SAS, which is sufficient for the full immobilization of the bubble surface [24,75,76]. A further increase in the SAS coverage, above this ‘threshold’ coverage, does not practically affect the bubble terminal velocity, as seen in Figure 10. The full retardation of the bubble surface mobility means that its motion is similar to that of a solid sphere of identical dimensions and density.

The terminal velocity of the bubble in solutions of SAS, often called ‘contaminated’ liquids, was the topic of many research studies, mainly focused on bubbles with fully immobilized surfaces, whose motion showed some similarities to the falling of rigid spheres. In the case of a rise in the laminar flow conditions, the drag coefficient can be described using the Stokes equation (Equation (14)), but even for such small velocities, deviations, first noted by Oseen [87], were observed. Ossen, after simplifying the Navier–Stokes equation by linearization, showed that the drag coefficient becomes larger compared to Equation (14) and for Re < 0.1 it can be calculated from [87]:

$$C_D=\frac{24}{Re}\left(1+\frac{3}{16}Re\right)$$

(41)

In the case of a higher Reynolds number (for Re < 800), Schiller and Naumann proposed the following formula [88]:

$$C_D=\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right)$$

(42)

Clift et al. [16] collected and analyzed data for ‘contaminated’ liquids, similarly as for bubbles’ motion in pure liquids, and presented the following relationship for bubble terminal velocities:

$$U_T=\frac{{\eta }_L}{{\rho }_L{d}_b}M{o}^{-0.149}\left(J_C-0.857\right)$$

(43)

where:

$$J_C=0.94H_C{}^{0.757}$$

(44)

for 2 < H_C < 59.3, and

$$J_C=3.42H_C^{0.441}$$

(45)

for H_C > 59.3, and H_C is described as:

$$H_C=\frac{4}{3}EoM{o}^{-0.149}$$

(46)

Karamanev [63] assumed that Equation (35), used for the clean water system, can be used for an estimation of the C_D in the case of a free rising sphere in ‘contaminated’ liquids. If the bubble in ‘contaminated’ liquid has a spherical shape, then there is no need to add a geometrical term in the formula for bubble terminal velocity. Thus, the terminal velocity of a free rising sphere should be calculated from Equations (7) and (34) for Ar < 1.18·10⁶d_b², and C_D = 0.95 for Ar > 1.18·10⁶d_b².

Nguyen et al. [89] also showed that the terminal velocity of a settling solid sphere can be directly predicted on the basis of the Archimedes and Lyaschenko numbers, but in a more straightforward way than that proposed by Karamanev [63]. Later, using a similar derivation [90], the formulas for bubble terminal velocity in the ‘contaminated’ system were presented. Nguyen [90] assumed that for Re < 130, the bubble’s shape is spherical and the bubble drag coefficient is equal to that of a solid particle. The following correlation for bubble terminal velocity was proposed:

$$U_T=\frac{d_b^{2}g\mathsf{\Delta }\rho }{18{\eta }_L}{\left[1+\frac{Ar/96}{{\left(1+0.079A{r}^{0.749}\right)}^{0.755}}\right]}^{-1}$$

(47)

When Ar is small, the second part of Equation (48) converges to unity, transforming the equation to the Stokes law.

For Re > 130, Nguyen [90] proposed, similarly to Karamanev, that C_D = 0.95, but in their formula for the bubble terminal velocity there are additional numerical parameters (a,b) describing the shape of rising bubble:

$$U_T=\sqrt[3]{\frac{g{\eta }_L}{{\rho }_L}}{\left(\frac{4a^{2}M{o}^{0.46b}}{2.85}\right)}^{1/\left(2-2b\right)}A{r}^{\frac{2b+1}{6-6b}}$$

(48)

The values of a and b depend on the Archimedes number and they are collected in the work of Nguyen [90].

Ng et al. [91,92] predicted the C_D with the use of the solution of the Oseen law, and for 0.2 < Re < 20,000 the following formula was derived:

$$C_D=\frac{24}{Re}\left(\frac{7}{6}R{e}^{0.15}+0.02Re\right)$$

(49)

Figure 11 presents comparisons of experimental data and the velocity values predicted from the models for bubble motion in clean and ‘contaminated’ water (i.e., in the presence of various SAS (n-pentanol, n-hexanol, n-octanol, n-nonanol, n-octyltrimethylammonium bromide (OTAB), and n-octylsulphate sodium salt (OctylS-Na)). The existing models describe relatively well the boundary cases (i.e., ultra-pure and fully ‘contaminated’ systems), while proper predictions for the U_T for intermediate states (i.e., the case of the partial immobilization of the bubble surface) are not covered by the models proposed. This issue was addressed in the work of Davis and Acrivos [93]. They proposed a model for the creeping flow of a ‘contaminated’ bubble that postulates the existence of a stagnant cap over the rear of the bubble, the size of which influences the drag force acting on the rising bubble. Later, Sadhal and Johnson [94] derived an analytical solution for the spherical stagnant cap. In their model, the bottom part of the rising bubble is completely covered with SAS, resulting in a tangentially immobile air–liquid interface, while the top part is free of SAS, giving a fully mobile boundary condition. They proposed a simple formula for the drag force of a ‘contaminated’ bubble:

$$F_D=2\pi d_b{\eta }_L\left[\frac{1}{4\pi }\left(2\psi +\mathit{sin}\psi -\mathit{sin}2\psi -\frac{1}{3}\mathit{sin}3\psi \right)+1\right]$$

(50)

where ψ is the stagnant cap angle (see Figure 8). For the limiting case of ultra-pure water, i.e., ψ = 0, Equation (50) gives the Hadamard–Rybczynski equation. Similarly, for a fully ‘contaminated’ bubble surface, i.e., ψ = π, Equation (50) transforms into the Stokes equation.

Figure 12 presents bubble terminal velocity calculated using the Sadhal and Johnson model [94] for an air bubble of a diameter equal to 0.1 mm. Assuming that the stagnant cap angle reflects the adsorption coverage, which is related to SAS concentration, then qualitative similarity to data, presented in Figure 10, especially theoretical lines, is seen. Unfortunately, a simple quantitative relationship between SAS adsorption coverage and the size of the stagnant cap is not available in the form of a direct analytical equation and needs to be calculated using numerical procedures [83,85].

Tomiyama et al. [70,71,72] proposed a different approach to the problem of the ‘contaminated’ system. Based on the previous works of Peebles and Garber [95] and Ishii and Chawla [96], they proposed formulas of the C_D for three ‘contamination’ regimes:

(A)

for a clean system

$$C_D=max\left\{min\left[\frac{16}{Re}\left(1+0.15R{e}^{0.687}\right),\frac{48}{Re}\right],\frac{8}{3}\frac{Eo}{Eo+4}\right\}$$

(51)

(B)

for a slightly ‘contaminated’ system

$$C_D=max\left\{min\left[\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right),\frac{72}{Re}\right],\frac{8}{3}\frac{Eo}{Eo+4}\right\}$$

(52)

(C)

for a fully ‘contaminated’ system

$$C_D=max\left\{\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right),\frac{8}{3}\frac{Eo}{Eo+4}\right\}$$

(53)

However, the conditions for a partially ‘contaminated’ system are not clearly described in the paper of Tomiyama et al. [70,71,72]. It is specified that a fully ‘contaminated’ bubble reaches its terminal velocities after covering a distance less than 17 mm from the point of bubble formation [71].

Yan et al. [97] made an attempt to deliver a single equation covering a wide range of fluid parameters (10⁻³ < Re < 10⁵, 10⁻² < Eo < 10³, and 10⁻¹⁴ < Mo < 10⁷) for a bubble rising in both pure and ‘contaminated’ liquids. Using nonlinear fitting to the experimental data, they proposed the following correlation:

$$C_D=max\left\{\begin{array}{l}\mathit{min}\left[\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right),\frac{72}{\mathit{Re}}\right],\\ \frac{24}{Re}\left(1+0.15R{e}^{0.687}\right)\frac{{\mathit{Re}}^{0.55}E{o}^{0.95}W{e}^{-1.10}}{12.6}{\left(\frac{Mo}{M{o}^{*}}\right)}^{0.048}\end{array}\right\}$$

(54)

where Mo* is the Mo of water at 293 K (2.6 × 10⁻¹¹).

Based on an analysis of the literature and experimental data for surfactants used as flotation frothers, Kowalczuk et al. [22] introduced the determination of a concentration at the minimum bubble velocity (CMV). Their empirical model allows us to predict the bubble terminal velocity as a function of the surfactant concentration (c) for a wide range of bubble diameters:

$$U_T=U_{\mathit{min}}+\left(U_{\mathit{max}}-U_{\mathit{min}}\right)e^{-3{\left(c/CMV\right)}^2}$$

(55)

In this model, it is necessary to know the maximum (U_max) and minimum (U_min) terminal velocities, as well as the value of the CMV. For the most popular surface active substances, the CMV data are collected in the paper [22]. The values of U_max and U_min can also be estimated using models for bubbles with fully mobile (pure liquid) and fully immobile (rigid sphere) interfaces, respectively.

5. Limitations and Future Directions

Most of the presented models are either empirical, semi-analytical, or analytical solutions derived with many simplifying assumptions. More sophisticated models of bubble motion in solutions of surface-active substance require the numerical solution of the Navier–Stokes equation [82,98], mainly with the implementation of the rear stagnant cap model [73,94,99,100,101] and a finite-rate mass exchange between the interface and bulk fluid [82,83,102,103,104]. Similarly, models describing the motion of a bubble in a stagnant liquid have a limited application in technological processes, where a gas–liquid flow is turbulent [18,105,106]. In such conditions, bubbles tend towards coalescences, breakups, and deformations [16,107,108,109]. Moreover, the motion of a bubble is characterized by its instability with continuous acceleration and deceleration stages [110,111,112]. The prediction of bubble motion becomes even more complicated for bubble swarms (groups of bubbles dispersed in a liquid phase), which are typical for industrial processes [113,114,115]. In this situation, advanced numerical methods, i.e., computational fluid dynamics (CFD), need to be implemented [116]. Particular attention should be paid to two CFD methods: direct numerical solutions (DNS) [112,117] and discrete element methods (DEM) [118]. The former technique allows us to resolve the turbulent fluid motion [119,120], but the price is high: most of the time it demands a computing power which can be assured by supercomputers only [120]. The discrete element method, which is less demanding in terms of computing resources, is implemented in the study of bubble–solid particles interactions, especially in terms of flotation processes [121,122,123,124]. The further development of DNS and DEM is rather inevitable, as both methods deliver significant and often satisfactory outcomes.

6. Conclusions

• Experimental methods to study bubble motion: visual observation (with use of cameras) is still the most reliable method for tracking bubbles. This method delivers the most comprehensive information about bubble motion, i.e., velocity, deformation, and path.
• Bubble motion in water: it was shown that predictions of the model proposed by Moore [20,21], supplemented by recent semi-empirical formulas (Legendre et al. [19]) describing the geometrical parameters of the rising bubble (deformation ratio), agree almost perfectly with experimental data for pure liquids.
• Bubble motion in liquid in the presence of surface-active substances (SAS): for SAS solutions, existing models describe relatively well the boundary case of fully ‘contaminated’ systems. It was shown that the semi-empirical model proposed by Kowalczuk et al. [22] is a very convenient tool for the description and prediction of bubble terminal velocities as a function of surfactant concentration for a wide range of bubble diameters.