Introduction

The concurrent existence of various fluid phases can be found in daily life and in a variety of chemical processes such as distillation, absorption, flotation, or in multiphase reactors (Sommerfeld and Horender 2012). Multiphase processes in reactors and columns use continuous contact of two or more phases, which is achieved mainly by a gaseous phase bubbling into a liquid phase. Knowledge of the physicochemical properties of both phases, i.e., density, viscosity, or surface tension, is crucial to design such processes, as they influence the hydrodynamic behavior of the whole system and bubble size distribution in the particular process. The main hydrodynamic properties of aerated apparatuses are mostly dependent on the bubble rise velocity, which defines (i) the bubble residence time, that is how long the bubble rises through the bulk, (ii) how long the bubble stays in contact with the liquid phase, which determines the mass transport, usually of oxygen (Kulkarni and Joshi 2005), and (iii) the gas hold-up in the column (Joshi et al. 1998). An accurate description of the bubble behavior is necessary for a proper estimation of the system hydrodynamics. To this day, numerous works exist in which the authors propose simplified theoretical models of the bubble motion in a liquid. Despite their effort, there is no universal model that could describe all bubble motion regimes in different liquids (Tripathi et al. 2015). In aerated columns apparatuses, it is not possible to characterize bubbles individually as they are rising in clusters; therefore, collisions and coalescence of bubbles can appear simultaneously. For a theoretical description and CFD modeling of this phenomenon, knowledge of small bubble clusters (in order of units of bubbles) and single bubble behavior is essential (Sanada et al. 2009).

With increasing computational power, progress has been made in modeling using CFD tools for direct numerical solution of Navier–Stokes equations. To resolve the motion of moving phase interface, there are two types of approaches—interface-capturing and interface-tracking methods (Tryggvason et al. 2001). The former method employs several tracking models: Volume of fluid (VOF), Level set (LS), Phase field (PF), and Constrained interpolation profile (CIP). The latter approach is used in Front-tracking method and Marker-and-cell method (Mirjalili et al. 2017). There are many solvers for fluid dynamics that are either commercial or are made as in-house code. COMSOL Multiphysics is a commercial CFD solver that uses the finite element method (FEM) discretization (COMSOL Multiphysics Reference Manual 2017). COMSOL allows a wide range of multiphase phenomena to be solved using Eulerian, Arbitrary Lagrangian–Eulerian, and Euler–Euler approaches. Two built-in interface-capturing models, solving the motion of the interface on an Eulerian fixed grid, can be found in COMSOL, namely the Level set and Phase field. The last built-in method uses the Arbitrary Lagrangian–Eulerian formulation (ALE), the Moving mesh method. These methods differ greatly in their physical background, resolution of the captured interface, and required computational resources (COMSOL Multiphysics Reference Manual 2017).

Many studies have dealt with the description and simulations of multiphase flow phenomena by benchmarking and testing various CFD solvers (both commercial and in-house codes). These studies can be divided into four main groups according to (i) used solver (commercial or in-house code), (ii) used method for interface computation, (iii) dimensions of the computational domain (2D axisymmetrical domain or fully 3D domain), and (iv) whether single or multiple bubbles were studied. The two-dimensional benchmark cases done by Hysing (2007) are still used to this day. In his work, comparison of their code TP2D with other in-house codes (FreeLIFE, MooNMD) and commercial codes (COMSOL Multiphysics 3.3a, Ansys Fluent 6.3) was done. Their code uses the FEM discretization method, and the bubble interface is resolved by the Level set method. The bubble shape deformations and their terminal rise velocity were studied. Klostermann et al. (2012) extended the work of Hysing by computing the benchmark cases using the VOF method in the open-source CFD solver OpenFOAM and proved its capabilities. Eiswirth et al. (2011) used COMSOL Multiphysics 3.3a to compute the rise characteristics and shape deformations of toluene droplets in water using the Level set method. Their results obtained from COMSOL were validated by experimental data and showed great agreement. Nichita et al. (2010) coupled the built-in Volume of fluid method with Level set method, using UDFs (User-defined functions), in the commercial solver Ansys Fluent and performed simulations for static and dynamic cases in 2D and 3D domains. Their results showed that Coupled Level set and Volume of fluid method perform better, compared to experimental results, than using the Volume of fluid itself. Tryggvason et al. (2001) used the front-tracking method to compute various multiphase flows phenomena, including homogenous bubbly flows, influence of surfactants on bubble dynamics, solidification, and other phenomena. Considerable work on simulation of fully three-dimensional bubble dynamics has been done by Tripathi et al. (2015). They managed to perform around 130 simulations of rising bubbles of various diameters using the Volume of fluid method in open-source code Gerris. Mass transfer during the rise of a single bubble in a fully three-dimensional domain was studied by Özkan (2016). Balcázar et al. (2015) simulated dynamics of fully 3D single bubble, multiple bubbles, and bubble pair interactions using Conservative Level set method in an in-house code called ThermoFluids. Their predicted values agree well with the experimental and numerical results from the literature.

From the above overview, it is obvious that in-house and open-source CFD codes are more used in the field of bubble motion modeling, as they offer more options in terms of code optimalization and adjustment. However, commercial codes are convenient to use and easy to implement because of their problem-solving capabilities in process optimization. They also offer a support from the developers and thus, are used by many companies, as they depend on the developers guarantee. Similarly, as in-house and open-source codes, commercial codes need to be verified on simple multiphase flow problems before simulating more complex phenomena.

Industrial liquid batches in aerated columns or reactors differ in physicochemical properties, especially in density, viscosity, and surface tension. However, most of the published studies focus on water, slightly viscous batches or idealized systems. This study deals with simulation of an isolated bubble in three model liquids: water (low viscosity, high surface tension), 1-propanol (low viscosity, low surface tension), and glycerol (high viscosity, high surface tension). The presence of contaminant, i.e., surface-active agents, is not considered. The study is focused on the description of spherical and slightly deformed bubbles that move in a straight line. For the simulation, the conservative Level set method in COMSOL Multiphysics 5.4 has been selected. This approach has been used successfully several times to describe various multiphase flow phenomena, but various fluids have not been tested (Crha et al. 2021; Gollakota and Kishore 2018; Hysing 2007). The numerical simulations of bubble dynamics have been conducted in the fully 3D domain, which follows on from our previous study done in 2D axisymmetric setup (Crha et al. 2021). The steady-state values of the terminal velocities and bubble deformations were compared with original experimental data as well with theoretical models, which exist for the single bubble case. Dynamic behavior of the rising bubble was also observed; however, it was not in the main scope of this work, and only the averaged values of velocities and the shape deformations in the steady-state region are mentioned.

Theoretical description

The steady-state bubble rise velocity (terminal velocity) has a crucial influence on the hydrodynamics of multiphase apparatuses. A simple expression of the bubble terminal velocity can be derived from a force balance, which are acting upon the bubble during the steady-state rise, where only gravitational, buoyant and drag force are considered (Tomiyama et al. 1998). This expression is usually written as follows:

$${U}_{b}=\sqrt{\frac{4\left({\rho }_{l}-{\rho }_{g}\right)g{D}_{b}}{3{C}_{D}{\rho }_{l}}} ,$$
(1)

where ρl stands for the density of the liquid, ρg denotes the density of the gaseous phase, g is the gravitational constant, Db is the diameter of the bubble, and CD is the drag coefficient. A wide range of expressions for the calculation of the drag coefficient exists. They differ in their range of applicability, which can be classified according to the Reynolds number of the bubble (Reb = Ubρl Db/μl), liquid contamination, and the shape of the bubble (Kulkarni and Joshi 2005).

The theoretical derivation of the drag coefficient of spherical and non-spherical bubbles exists only for limiting cases of flow, where surface tension and/or inertial effects can be neglected. For flow past a spherical bubble at low Reynolds numbers (Re < 1), the flow field stays symmetrical and there is no wake behind the bubble. In pure systems, the liquid slips along the bubble surface, as it rises through the liquid bulk and the no-shear stress condition can be imposed on the surface. The following expression was derived by Hadamard and Rybczynski (the detailed derivation can be found in (Clift et al. 1978)) on the basis of the assumptions mentioned above:

$${C}_{D, H-R}=\frac{16}{\mathrm{Re}} .$$
(2)

However, this expression is not useful for many applications, as bubbles in most systems do not satisfy the restriction of a small Reynolds number value (Re < 1).

Mei et al. (1994) proposed a drag law based on numerical simulations for spherical bubbles applicable for a wide range of Reynolds numbers. The Mei expression for the drag coefficient can be found in this form:

$${C}_{\mathrm{D},\mathrm{Mei}}=\frac{16}{\mathrm{Re}}\left\{1+{\left[\frac{8}{\mathrm{Re}}+\frac{1}{2}\left(1+3.315{\mathrm{Re}}^{-1/2}\right)\right]}^{-1}\right\} .$$
(3)

This function is valid for all values of Reynolds numbers; however, it neglects bubble deformations, which arise at higher Re values. In the air–water system, shape deformations start to appear for bubbles with diameter above Db ≥ 0.83 mm (Legendre et al. 2012). The bubble deformation is usually defined by the bubble aspect ratio (χ), which is a ratio of the semi-major to semi-minor axis in a symmetrical ellipsoid, and it acquires values larger than 1 (value of 1 is for a sphere). Moore (1965) studied terminal velocities of small-distorted bubbles rising in a stagnant liquid of small viscosity. He extended his theory from previous work for the boundary layer on a spherical bubble and found a relationship between the aspect ratio and the Weber number of the bubble (Web = DbρUb2/σ). This expression is written in this way:

$$\chi =1+\frac{9}{64}\mathrm{We}+O\left({\mathrm{We}}^{2}\right) ,$$
(4)

where O(We2) is the correction factor. Rastello et al. (2011) expanded this expression by taking the results of Re-χ dependency from Blanco and Magnaudet (1995) and fitted them with a least-square fit. The resulting correlation is written in this way:

$$\upchi =1+\frac{9}{64}\mathrm{We}+\frac{3}{250}{\mathrm{We}}^{2}+O\left({\mathrm{We}}^{3}\right) .$$
(5)

The Moore’s drag coefficient for deformed bubbles can be calculated using this expression:

$${C}_{D,\mathrm{Moore}}=\frac{48}{\mathrm{R}e}G\left(\upchi \right)\left[1-\frac{2.21H\left(\upchi \right)}{{\mathrm{Re}}^{1/2}}\right] ,$$
(6)

where G(χ) and H(χ) are geometrical factors. Their complete form can be found in (Moore 1965), but more often the polynomial approximation by Loth (2008) is used:

$$G\left(\upchi \right)\simeq 0.1287+0.4256\,\upchi +0.4466\,{\upchi }^{2} ,$$
(7)
$$H\left(\upchi \right)\simeq 0.8886+0.5693\,\upchi -0.4563\,{\upchi }^{2}.$$
(8)

This form of geometrical factors is more straightforward to use in calculations, and the accuracy is guaranteed. Moore’s drag coefficient is applicable in the range of Re > 100, We < 4, and χ < 2 for bubbles with unseparated wake. For the higher values of the Weber number, the values of the CD are underestimated, which is primarily due to the onset of wake separation (Loth 2008). Rastello et al. (2011) modified Mei’s expression to account for the bubble’s deformation and it can be expressed in this form:

$${C}_{D,\mathrm{Rastello}}=\frac{16}{\mathrm{Re}}\left[\frac{1+\frac{8}{15}\left(\upchi -1\right)+0.015\left(3G\left(\upchi \right)-2\right)Re}{1+0.015\mathrm{Re}}+{\left[\frac{8}{\mathrm{Re}}+\frac{1}{2}\left(1+\frac{3.315H\left(\upchi \right)G\left(\upchi \right)}{{\mathrm{Re}}^\frac{1}{2}}\right)\right]}^{-1}\right]$$
(9)

In industrial applications, most liquids are not completely pure and can contain surface-active impurities (surfactants) or they are added into the liquid batch on purpose—i.e., surfactant-like behavior of diluted aqueous solutions of low-carbon alcohols, which are used for surface tension reduction (Basařová et al. 2022). In these systems, surfactant molecules preferably adsorb on the surface of the bubble and are advected by the flow of the surrounding fluid to the rear of the bubble and partially immobilize the surface. In certain cases, the immobilization can affect the whole surface and the bubble acts as a rigid sphere (Cuenot et al. 1997). In that case, a no-slip condition should be applied to the bubble surface for modeling purposes. Assuming a rigid particle/bubble, settling, or rising in a stagnant viscous liquid, moving at low velocity (Re <  < 1) Stokes derived the following drag law:

$${C}_{D,\mathrm{Stokes}}=\frac{24}{Re} .$$
(10)

Another expression for the rise of an immobile bubble was found out by Schiller and Naumann, which is mentioned in many studies, i.e. (Manica et al. 2014). This empirical correlation is valid for a large range of the Reynolds number of the bubble (Re < 800) and is written in this form:

$${C}_{D,S-N}=\frac{24}{Re}\left(1+0.15{\mathrm{Re}}^{0.687}\right) .$$
(11)

The drag coefficients mentioned above are only a fragment of a wide range of existing correlations. These were chosen and described because they are widely used and part of them was also used to calculate theoretical velocities of rising bubbles to validate results from COMSOL Multiphysics.

Experimental

The experiments were carried out in a glass vessel having dimensions of 30 cm in height, 8 cm in width and 6 cm in depth. The scheme of the apparatus is visualized in Fig. 1. Pressurized air (1) comes from an air compressor to the air chamber (2), right below the glass cell (3). The glass column has dimensions of 30 cm in height, 8 cm in width, and 6 cm in depth. The bubble is generated on an exchangeable capillary having the inner diameter 10 μm, which is connected to the bubble generator. The diameter of the bubble is controlled by the filling time of the air supplied into a chamber beneath the cell. After a specified period of time, the bubble is released by a quick pulldown of the capillary. The delay before bubble release is to ensure that the liquid is stagnant and free of any motion that could be caused by previous bubble motion. The detailed description of the experimental apparatus’s functionality can be found in (Basařová et al. 2018). The freely rising bubble was captured by a high-speed camera (4) Photron SA1.1 at 1000 fps with a Navitar macro-objective. The size of the captured images was 1024 × 1024 px with an average resolution of 2–3 px/μm.

Fig. 1
figure 1

Schematic of experimental apparatus (1) pressurized air inlet (2) air chamber (3) glass cell (4) high-speed camera (5) cold light source

Full size image

The use of single camera setup allows obtaining only orthographic projection of the bubble in the plane of focus; therefore, it is not suitable to gather full information about the trajectory of the rising bubble. In our experiments, all captured bubbles stayed sharp, meaning bubbles were rising in the plane of focus. This and other preliminary experiments confirm that bubbles had a rectilinear trajectory.

The desired quantities, such as the coordinates of the bubble center, aspect ratio, and bubble diameter, were obtained for each image in the image sequence, using image analysis software NIS-Elements. The velocity of the bubble was determined from this expression:

$${U}_{b}=\frac{\sqrt{{\left({x}_{2}-{x}_{1}\right)}^{2}+{\left({y}_{2}-{y}_{1}\right)}^{2}}}{{t}_{2}-{t}_{1}} ,$$
(12)

where x and y are horizontal and vertical positions of the bubble center, t is time and indexes 1 and 2 denote image’s order in the sequence. The velocity, obtained from an image sequence, was calculated as an average value of at least 20 consecutive images.

The bubble terminal velocity was measured in three different liquids, namely distilled water, 1-propanol and glycerol. The distilled water was purified by the “ULTRAPURE” system produced by Millipore to provide the highest possible purity. 1-Propanol of > 99.5% purity from Penta and pharmaceutical grade glycerol with declared purity of > 99% were used without further purification. Density and surface tension were measured using the Kruss K11 tensiometer, using the Du Noüy ring method; dynamic viscosity was measured using a Stabinger Viscometer produced by Anton Paar (SVM3000). The accuracy of the measured physical quantities is 0.1%. The measured values of the physicochemical properties are described in Table 1.

Table 1 Values of experimental physical quantities used as input in the simulation
Full size table

The dynamic viscosity and density of pure glycerol were calculated using the empirical model found in (Volk and Kähler 2018) for a temperature of 22.5 °C. The value of the dynamic viscosity should be 1127 mPa s. However, the experimentally obtained value was 575 mPa s, which corresponds to the viscosity of the 96% aqueous solution of glycerol, which means that glycerol absorbed air moisture during the experiment.

Computational methods

The bubble motion was simulated in COMSOL Multiphysics 5.4, which is a commercial CFD solver based on the Finite element method (FEM). The fluid flow for both phases is resolved by the Navier–Stokes equations (Eq. 13) together with the continuity equation (Eq. 14):

$$\rho \left(\frac{\partial {\varvec{u}}}{\partial t}+{\varvec{u}}\cdot \nabla {\varvec{u}}\right)=\nabla \cdot \left[-p{\varvec{I}}+\mu \left(\nabla {\varvec{u}}+{\left(\nabla {\varvec{u}}\right)}^{T}\right)\right]+{\varvec{F}},$$
(13)
$$\uprho \nabla \cdot \left({\varvec{u}}\right)=0,$$
(14)

where ρ is the fluid density, u is the velocity vector, p is the pressure, μ is the dynamic viscosity, I is the identity matrix, T stands for matrix transpose, and F is other external forces, such as surface tension and gravity. These equations were solved together with initial conditions of zero velocity field in whole domain, initial position of the bubble, and boundary conditions, namely no-slip walls and pressure outlet at the top of the domain. The motion of the bubble interface was resolved by the Level set method.

Level set method

The Level set method was introduced by Osher and Sethian (1988). The interface is represented by an implicit function; the so-called level set function, Φ, defined as a signed distance function with a range of values from 0 to 1 across the phases. The interface itself is defined as a 0.5 contour of the level set function. The values of density and viscosity greatly differ across the interface; therefore, the smeared-out Heaviside function H(Φ) is introduced:

$$H\left(\Phi \right)=\left\{\begin{array}{c}0\\ \frac{1}{2}+\frac{\Phi }{2\varepsilon }+\frac{1}{2\pi }sin\left(\frac{\pi\Phi }{\varepsilon }\right)\\ 1 \end{array}\right.\begin{array}{c} for\,\Phi <-\varepsilon \\ for -\varepsilon \le\Phi \le \varepsilon \\ for\,\Phi >\varepsilon \end{array} .$$
(15)

Here ε is the parameter that controls the thickness of the interface, set by default as half of the largest mesh element that could be encountered by the moving interface. The main weakness of the Level set method is the fact that it is not naturally conservative as the VOF method. To overcome the issue arising from the loss of mass, the conservative level set method was introduced by Olsson and Kreiss (2005). This formulation was selected in COMSOL Multiphysics to ensure that the mass of the bubble would be conserved.

The motion of the level set function is maintained by an advection equation, which in conservative form reads as follows:

$$\frac{\partial\Phi }{\partial t}+\nabla \cdot \left({\varvec{u}}\Phi \right)=\upgamma \nabla \cdot \left(\upvarepsilon \nabla\Phi -\Phi \left(1-\Phi \right)\frac{\nabla\Phi }{\left|\nabla\Phi \right|}\right),$$
(16)

where γ is the reinitialization parameter defining the interval of the level set function reestablishment. The level set function is also used to define the physicochemical properties (density and viscosity) of each phase:

$$\rho ={\rho }_{1}+\left({\rho }_{2}-{\rho }_{1}\right)\Phi ,$$
(17)
$$\mu ={\mu }_{1}+\left({\mu }_{2}-{\mu }_{1}\right)\Phi ,$$
(18)

in which index 1 belongs to the part of the domain, where the level set function acquires a value less than 0.5 and index 2 marks the region, where the function’s value is larger than 0.5.

The surface tension forces in COMSOL are implemented in the source term F in Eq. 13. The used model is the continuous surface force (CSF) model, introduced by Brackbill et al. (1992) and is given by:

$${F}_{\mathrm{ST}}=\mathrm{\sigma \kappa \delta }n$$
(19)

where σ is the surface tension, κ is the surface curvature, δ is the Dirac function, and n is the unit normal vector to the surface. By definition, surface tension is a surface force. However, during the calculation of the force, it is converted to a volume force. This conversion can be written in the following form:

$${F}_{\mathrm{ST}}=\mathrm{\sigma \kappa \delta }n=\mathrm{\sigma \kappa }\left(\upphi \right)\nabla\upphi$$
(20)

in which κ and n are given by

$${\varvec{n}}=\frac{\nabla\upphi }{\left|\left|\nabla\upphi \right|\right|},\upkappa \left(\upphi \right)=-\nabla \cdot {\varvec{n}}$$
(21)

Simulation setup

The rise of the bubble was computed non-stationary in a fully three-dimensional computational domain, built as a square rectangular cuboid with lengths of the square base of 6 mm and height of at least 20 mm, depending on the studied bubble size and liquid to ensure that the bubble will reach its terminal velocity and will not be influenced by the walls. Krishna et al. (1999) studied this phenomena experimentally, and their results showed that the ratio of bubble diameter to the diameter (width) of the domain should be below 0.125. However, their study dealt with large bubbles (3–80 mm) rising in water, which rise non-rectilinearly, usually at zigzag or helical path. The ratio of 0.125 is significantly overestimated for our arrangement with the assumed rectilinear movement of bubbles. Mukundakrishnan et al. (2007) studied the wall effects numerically using a front-tracking method together with a level contour reconstruction procedure. The numerical simulations were performed for different initial bubble diameters, rising in various liquids and in variously wide domains. Their results showed that the computational domain should be at least three times wider and eight times higher than the bubble diameter, which our domain fulfills.

The initial bubble diameters in our simulations were in the range of 1.0–1.6 mm, where different values of steady-state aspect ratios are expected. The domain was divided into two parts and discretized with different sizes of mesh elements. The subdomain surrounding the rising bubble was constructed as a cylinder with a diameter at least twice larger than the diameter of the bubble. The free tetrahedral computational mesh was used for both domains, and the size of elements was selected according to results from mesh-dependency study, which has been done in (Crha et al. 2021) for single rising bubble in 2D axisymmetrical domain. The size of the inner mesh elements ranged from 0.1 mm to 0.2 mm, depending on the diameter of the bubble. The outer domain was meshed using elements whose size varied spatially, from smaller at the boundary of the inner domain to maximum size at the walls of the outer domain. The whole domain consisted of approximately 4 million elements.

The iterative segregated approach was selected instead of the direct fully coupled approach to avoid a large consumption of the system’s RAM. The Algebraic multigrid solver method, using the GMRES (generalized minimal residual method) solver, was used, and the absolute tolerance was set to 0.001. The time stepping method Generalized alpha was used with free time step choice. The fluid flow was discretized using P3 + P2 elements, meaning third-order elements for the velocity components and second-order elements for the pressure field. The level set variable was discretized using quadratic elements. The results for postprocessing were saved every 0.001 s.

The bubble rising velocity was evaluated in each time step, using the following expression in the two-dimensional form (Hysing 2007):

$${U}_{\mathrm{b}} =\frac{{\int }_{{\Omega }_{2}}\mathrm{w\,dz}}{{\int }_{{\Omega }_{2}}1\,\mathrm{dz}} ,$$
(22)

in which w is the vertical component of the velocity field, and Ω2 denotes the region of the domain occupied by the bubble. The terminal velocity is an average value of at least twenty values of instantaneous velocities. The bubble shape deformation χ and the equivalent diameter Db were determined by an image analysis of the bubble image sequence in steady state, which was acquired from an orthographic projection of the yz-plane. The steady-state value of the bubble aspect ratio is also an average of ten values, and the value of the equivalent diameter served to choose the corresponding experimental data. The trajectory of the bubble center of mass was also evaluated to analyze whether the motion is rectilinear:

$$\mathrm{Bubble\,\,center}=\frac{{\int }_{{\Omega }_{2}}\mathrm{i\,di}}{{\int }_{{\Omega }_{2}}\mathrm{di}} ,$$
(23)

where i is the i-th component of the coordinate system. The degree of the non-rectilinearity was defined as a standard deviation from the ideal motion of the bubble center. Standard deviations for all bubbles rising in all liquids were negligible. In the case of the bubbles in 1-propanol and glycerol, the deviation of the bubble center in thy xy-plane was below the size of one mesh element. In the case of bubbles in water, the deviations were higher. However, the bubble center did not travel further than 2.25 times the size of the mesh elements.

Additionally, the values obtained from our simulations of Galileo and Eötvös numbers were compared against the chart of bubble rise regimes in (Tripathi et al. 2015). Bubbles rising in 1-propanol and glycerol clearly belong to the axisymmetric regime. Values of Ga and Eo of larger bubbles (1.4 and 1.5 mm) rising in water are on the overlap of the axisymmetric and asymmetric regions. For further analysis of the trajectory of the bubble in water, the fluctuations of instantaneous velocities obtained from COMSOL were also studied, and they did not show any mentionable deviations; thus, the trajectories of the bubbles in water can be categorized as rectilinear. The Galileo number Ga = ρlg1/2Rb3/2/μl expresses the ratio of the gravity force to the fluid dynamic viscosity. The Eötvös number Eo = ρlgRb2/s reflects the influence of surface tension s. The smaller the Eötvös number, the stronger the influence of surface tension is. Finally, the Morton number Mo = Eo3/Ga4 characterizes the multiphase fluid properties. Rb is the radius of the bubble.

The reinitialization parameter γ in Eq. 16 was tuned according to the expected velocity magnitude in the domain; for water and 1-propanol, it was selected as a 1 m s−1 and for glycerol 0.002 m s−1. In COMSOL, the bubble terminal velocities were calculated using Eq. 22, which was implemented by applying the integration operator intop in the whole computational domain. The time-dependent values were obtained in the postprocessing section as a global evaluation, and the final value of the terminal velocity was taken as an average value from at least twenty time steps in the region of steady-state rise. The average calculation time of one bubble rise was approximately in the range of 80–180 h (Fig. 2).

Fig. 2
figure 2

a Domain with the bubble in its initial position b Top view of the meshed domain

Full size image

Results

The results are given in Table 2 and in Figs. 3, 4, 5, and 6. Here, the data from COMSOL simulations are compared with both experimental data and results calculated according to theoretical models. The input values of the physicochemical properties are summarized in Table 1. In Table 2, the values of commonly used dimensionless quantities such as Reynolds, Weber, and Eötvös numbers are displayed. Initial diameters of bubbles from 1.0 to 1.6 mm were chosen to reach different values of aspect ratio and to achieve different regimes of the bubble deformation (loss of fore- and aft-symmetry). The results are divided into separate paragraphs according to the studied liquid. At first, the terminal velocities are evaluated and compared with the experimental data and theoretically calculated velocities. In the following paragraphs, the shape deformations are assessed. Finally, our results are compared to existing data from the literature, in cases where they are available.

Table 2 Experimental, theoretical, and calculated values of aspect ratio c, terminal bubble velocity U, Weber We, Eötvös Eo, Reynolds number Re
Full size table
Fig. 3
figure 3

Steady-state velocities and aspect ratios obtained from COMSOL (green triangles, rhombi) for different bubble sizes, rising in 1-propanol, compared to experimental values (red circles, squares) and theoretical values (solid line, dashed line) (calculated using Eq. 9)

Full size image
Fig. 4
figure 4

Steady-state velocities and aspect ratios obtained from COMSOL (green triangles, rhombi) for different bubble sizes, rising in water, compared to experimental velocities (red circles) and theoretical velocities (solid line) (calculated using Eq. 9, with experimental values of aspect ratio)

Full size image
Fig. 5
figure 5

Steady-state velocities obtained from COMSOL for different bubble sizes, rising in glycerol, compared to experimental velocities and theoretical velocities (calculated using Eq. 10)

Full size image
Fig. 6
figure 6

Graphical comparison of the instantaneous values of the aspect ratio. Experiment: a 1.6 mm bubble rising in 1-propanol b 1.5 mm bubble rising in water c 1.2 mm bubble rising in glycerol COMSOL: d 1.6 mm bubble rising in 1-propanol, e 1.5 mm bubble rising in water, f 1.2 mm bubble rising in glycerol

Full size image

The simulation setup differed in the time length of the bubble rise to achieve steady-state velocity. In COMSOL, the bubble terminal velocities were calculated using Eq. 22, which was implemented by applying the integration operator intop in the whole domain. The time-dependent values were obtained in the postprocessing section as a global evaluation, and the final value of the terminal velocity was taken as an average value from at least twenty time steps in the region of steady state.

The results of the simulations were compared with the experimental data and theoretical models. However, it is necessary to realize how individual quantities depend on each other in the calculation. When calculating the velocity of a bubble of a given size (according to Eq. 1), we must know the drag coefficient. If we use Eq. 6 to calculate it, it requires the deformation of the bubble. The aspect ratio can be estimated using the Weber number, which is a function of bubble velocity. Therefore, bubble velocity calculations according to theoretical models are always iterative. In Table 2, the method of the velocity calculation Utheor is specified for each liquid further in the text. The accuracy of the relationships (Eqs. 4 and 5) allowing the estimation of aspect ratio on the Weber criterion has not been tested on our experimental fluids. The estimation error would then significantly affect the results of the comparison of the individual models.

1-propanol

1-Propanol is a liquid with low surface tension and low viscosity. The experimental data are in excellent agreement with the theoretical model. The theoretical velocity calculation was performed iteratively. For 1-propanol, the mobile bubble surface was assumed. The aspect ratio value χtheor was calculated according to Eq. 5, the drag coefficient according to Eqs. 7–9, and the terminal velocity according to Eq. 1. In Fig. 3, the data for the terminal velocities from the experiment are depicted as circles, the theoretical data are marked with a solid line, and the CFD velocities are displayed as triangles. The resulting aspect ratios obtained from COMSOL, experiment, and Eq. 5 are also shown in Fig. 3.

The simulations of bubble motion were done for initial diameters of 1, 1.2, 1.4, 1.5, 1.6 mm, and the results are outlined in Fig. 3 as full squares. The CFD velocity (139 mm s−1) of the 1 mm bubble is slightly higher than the measured velocity (132 mm s−1). This agrees well with the smaller bubble deformation obtained from COMSOL than the observed value (χCFD = 1.05 vs. χexp = 1.07). The explanation is simple. The drag force acts on the area (bubble projection). The lower value of the aspect ratio means a smaller bubble projected area; thus, the drag coefficient is lower, and the velocity is higher. This trend of higher velocity as a result of the lower bubble deformation can also be observed for the cases of bubbles of 1.2 and 1.4 mm. The computed velocity of the 1.2 mm bubble is higher than the experimental value (168 mm s−1 vs 162 mm s−1), and the corresponding of aspect ratios is χCFD = 1.121 vs. χExp = 1.15. The resulting velocities from COMSOL for the bubbles of initial diameter 1.5 and 1.6 mm are slightly lower than experimental velocities (199, 205 mm s−1 vs. 203, 214 mm s−1), and the difference is below 5%. The aspect ratio values of both bubbles agree with the experimental data. In general, the resulting velocities for 1-propanol agreed excellently with the experimental and theoretically calculated velocities. The velocities obtained from COMSOL differ from experimental values, averagely around 5% in absolute values, and this difference is also valid for the comparison with theoretically calculated velocities (Eq. 9). The aspect ratio values from COMSOL do not largely differ from experimental and theoretical values.

The simulation results were not validated against other numerical studies, as there was not found any, which was dealing with the single rising bubble in 1-propanol. The total time of one simulation was in the range of 75–98 h, depending on the diameter of the bubble. The simulation converged well and did not exhibit any complications during the calculation. The mass conservation of the bubble was also monitored, and the mass of the bubble stayed constant during the rise.

Water

Water is a liquid with high surface tension and low viscosity. Bubble velocities in water are often greatly affected by the presence of surfactants, even in small amounts, and therefore, great emphasis needs to be placed on the accuracy of the experimental data. Whereas in glycerol and 1-propanol, the degree of surfactant effects is lowered due to the high viscosity and low surface tension, respectively. To confirm that the water is pure of any surface-active agents, it is important to compare the experimental velocities, which must match the estimates for the mobile bubble surface and measure the surface tension value, which should match the values available in the literature. Any deviations from the reference values can indicate the presence of surfactants.

The velocity calculation was performed according to Eq. 1, where the Rastello relation (Eqs. 7–9) along with the assumption of a mobile bubble surface was used for the drag coefficient. The experimental aspect ratio was used for the calculation. Estimated aspect ratio from Eq. 5 is not accurate, and the error for a bubble with a radius of 1.5 mm was 9%. The Legendre estimation (2012) is also recommended for water and gives better results, as the error was 3%. Even this error, in the calculation of the aspect ratio, could lead to incorrect comparisons for velocities, and thus, the experimentally obtained values of the aspect ratio were used for the calculation of theoretical velocity.

Although we have experimental data for the entire bubble size range, the simulations of bubble motion in water converged for bubbles of initial sizes 1, 1.4, and 1.5 mm. Simulations of other bubble sizes did not converge, which will be discussed later. The results are outlined in Fig. 4. The computed velocities from COMSOL (depicted as green triangles in Fig. 4) did not agree that greatly for all bubble sizes as in the case of 1-propanol. The best agreement with the experimental data was for the case of a 1 mm bubble, where the difference between the CFD velocity and the experimental velocity was below 3.5% (269 mm s−1 vs. 279 mm s−1), although the aspect ratio value obtained from COMSOL was almost 10% lower than the experimental value (χCFD = 1.07 vs. χexp = 1.18). This trend of underestimation of bubble deformation is also seen in the case of the 1.4 and 1.5 mm bubbles, where the aspect ratio was 19% and 20% lower, respectively, than the experimental value (χCFD = 1.22 vs. χexp = 1.51 and χCFD = 1.28 vs. χexp = 1.61). The CFD computed terminal velocities are 13% lower in both cases. The velocities of the 1.4 mm and 1.5 mm bubbles are 305 and 313 mm s−1, while the corresponding experimental values are 352 and 359 mm s−1. The experimental values of terminal velocity values were validated by results available in the literature (Duineveld 1995) and are in excellent agreement with them. The resulting velocities obtained from COMSOL were compared with previous research (Crha et al. 2021), which dealt with bubbles rising in the 2D axisymmetrical domain. In 2D axisymmetrical simulation, the 1.5 mm bubble reached a terminal velocity of 250 mm s−1, which is significantly lower. Therefore, the 2D axisymmetrical setup is not suitable for larger deformed bubbles. The terminal aspect ratio values obtained for all bubbles from COMSOL are underestimated when compared to those obtained by the experiment. The comparison with the calculated values of the aspect ratios using Eq. 5 is missing, as this expression was not used for aspect ratio estimation, and the theoretical velocities were calculated with inputted values of aspect ratio obtained from experiment.

In general, the simulations converged well, with the simulation time of a single bubble rise around 180 h. The velocities and shape deformations of bubbles with initial diameters of 1.4 and 1.5 mm do not agree with the experimental values. The bubble deformation clearly corresponds (see Eqs. 6, 9) to the drag coefficient value and therefore, to the terminal velocity; the higher the aspect ratio, the lower terminal velocity of a bubble with constant equivalent diameter. However, the bubble terminal velocities are in clear contradiction with obtained values of aspect ratio when compared to experiment. Lower values of the aspect ratios obtained from COMSOL should lead to a higher value of the terminal velocity. This could be due to numerical errors during the calculation of the surface tension force, which is crucial for correct calculation of the bubble shape deformation and thus, the correct bubble terminal velocity.

Glycerol

Glycerol is a liquid with high surface tension and high viscosity. Unfortunately, glycerol is also highly hygroscopic and can absorb air moisture when the experimental apparatus is opened to the surrounding atmosphere. As our glycerol, used during the experiments, contained almost 4% water, there was a significant decrease in viscosity compared to pure glycerol (see Table 1). Figure 5 illustrates the resulting velocities. The experimental velocities are indicated by red circles in Fig. 5. Theoretical velocities, calculated for the dynamic viscosity of 575 mPa s, are indicated by a dark gray solid line. Due to the dominant effect of viscosity and surface tension and consequently low velocity of the bubble (low inertia effect leads to low value of Reynolds and Weber number—Re, We <  < 1), the bubbles have a completely spherical shape (Loth 2008), and no bubble shape deformations were considered for glycerol. Therefore, the Mei et al. (1994) relation for the drag coefficient (Eq. 3) was used. The experimental velocities are lower than theoretical, computed using the Mei expression (Eq. 3). This indicates the presence of impurities that adsorb to the phase interface and reduce the mobility of the bubble surface. For a bubble of 1.0 mm diameter, the experimental velocity is 1.508 mm s−1 and the theoretical velocity, assuming full surface mobility, is 1.78 mm s−1. If we assumed an immobile bubble surface, i.e., calculation of the drag coefficient according to Eq. 8, the velocity of the rising bubble would be equal to 1.19 mm s−1. Similarly for the 1.4 mm diameter bubble, the experimental velocity is 3.24 mm s−1, the theoretical velocity assuming full surface mobility is 3.48 mm s−1, and the theoretical velocity assuming full surface immobility is 2.33 mm s−1. From the velocity comparison, the surfaces of the bubbles are evidently partially immobilized. However, the degree of immobilization cannot be directly determined only from the velocity comparison. The experimental data could not be used to validate the simulation, as data show that the bubble behavior is influenced by the surfactant presence. The simulation of the bubble motion with the viscosity value of 575 mPa s was done only for the initial diameter of 1.4 mm. The final velocity was 1.99 mm s−1, which is slightly lower than the value for the immobile surface of the bubble (2.31 mm s−1).

Further simulations were performed for the case of pure glycerol with a dynamic viscosity of 1127 mPa s. The initial bubble diameters were 1, 1.2, and 1.4 mm. Due to the low value of the Reynolds number, it is possible to consider the creeping flow, where two limit analytical solutions exist—i) Stokes’ drag law for the immobile bubbles (Eq. 10) and rigid spheres or ii) the Hadamard–Rybczynski solution (Eq. 2) for the mobile bubbles (Manica et al. 2016). In the region of low Reynolds number (Re < 0.1), where the difference between H-R (Eq. 2) solution and expression of Mei (Eq. 3) is below 1%, the drag coefficient obtained from Mei (Eq. 3) equals to the H-R solution and the velocity can be calculated using Eq. 3. The resulting velocities obtained from COMSOL (green triangles) and calculated using Eq. 3 (gray dashed line) and Eq. 10 (black dashed line) are shown in Fig. 5. The terminal velocities computed by COMSOL agree with the values calculated using the Stokes drag coefficient, which is in contradiction with expected results—that the bubbles will rise with a velocity predicted by the theoretical results considering mobility of the surface (Eq. 3).

This shows that COMSOL largely overestimates the resulting drag of the bubble in the case of highly viscous liquids and the terminal velocity matches the velocity of an immobile bubble. This applies for both cases of tested glycerol viscosity values (1127 and 575 mPa s). However, from the resulting velocity field, it is clear that the recirculation of the bubble is not reduced, and the surface is still mobile. The exact cause of the computed low velocity is not clear, and the dependence of the simulation settings (reinitialization parameter, absolute tolerance, etc.) on the bubble velocity should be further studied.

Shape deformations

The shape deformations of bubbles obtained from CFD simulations were compared to the experimental data and, for bubbles in 1-propanol and glycerol, to theoretically calculated values using Eq. 5. Both the simulation and experiment values were obtained by image analysis of the bubble rise sequence using the NIS-Elements image analysis software. The steady-state value for each bubble is an average value of at least 20 images. The visual comparison is outlined in Fig. 5. Bubble sizes are approximate because different pixel calibrations were used in the experiments.

The correctly computed shape deformations are essential for the accurate description of the fluid flow around the bubble and hence, for the correct bubble terminal velocity. In Fig. 5, the steady-state values of the aspect ratios, obtained from COMSOL, are compared to the deformations observed during the experiment. From a visual point of view, the degree of shape deformation seems highly similar. However, even small divergence (in units of percent) of the aspect ratio value can cause a large difference in the terminal velocity. In the case of 1-propanol and glycerol, the values of averaged aspect ratios computed by COMSOL for all bubble diameters agreed with the experimental results and with values calculated from Eq. 5. In the case of bubbles rising in the water, the CFD simulation fails and does not compute the bubble deformation correctly. The reason is not exactly clear. Substituting into Eqs. 7 and 8 (for G(χ) and H(χ)), Eq. 9 (for CD) and Eq. 1 (for Ub), it can be easily proved that as the aspect ratio χ increases, the drag coefficient increases and the velocity decreases. A detailed analysis of the pressure and velocity fields around the rising bubble would probably provide an explanation. Several studies (Chakraborty et al. 2013; Harvie et al. 2006; Klostermann et al. 2012) also mention parasitic or spurious currents. These nonphysical errors arise due to inaccurate calculation of the surface tension force, which results in incorrect shape deformations of the bubble. However, a detailed analysis of spurious currents and resulting velocity and pressure fields was not the main objective of this work.

Comparison with the literature

It was intended to directly compare the computed velocities and other quantities with the data available in the literature. Therefore, all dimensionless quantities (Re, We, Eo, Mo, and Ga) used in the cited studies were calculated. The Reynolds, Weber, and Eötvös numbers are shown in Table 2. The Morton number is constant for each liquid and the values were (1.75·10–8; 1.64·10–11; 50.18) for 1-propanol, pure water, and glycerol, respectively. The corresponding ranges of the Galileo number for the smallest and largest bubbles were (13.7; 27.8) for 1-propanol, (39.8; 79.5) for water, and (0.08; 0.16) for glycerol. Most of the numerical studies, mentioned in the theoretical part, deal with simulation of larger bubbles (i.e., (Balcázar et al. 2015), (Gumulya et al. 2016)) and idealized two-phase systems (low viscosity or density ratios, low value of surface tension—i.e. (Klostermann et al. 2012)). Other studies assume that the system can be modeled using 2D planar domain, which is applicable only for a narrow range of systems (i.e., Chakraborty et al. 2013; Islam et al. 2015; Krishna and Van Baten 1999)). For example, a bubble motion in a narrow channel can be modeled by this approach (Šimčík 2008).

Comparable results were found in the publications of Sharaf et al. (2017) and Tripathi et al. (2015). In the first publication, bubbles rising in water, aqueous solutions of glycerol and pure glycerol are observed experimentally and the CFD results from the latter publication were used for comparison. The velocities and shape deformations of bubbles rising in pure 1-propanol could not be directly compared, as they are not present in any published study. Nevertheless, the bubbles with low values of Galilei and Eötvös number (Ga < 30, Eo < 10) should be categorized in axisymmetric regime, meaning rectilinear rise without any path oscillations (Tripathi et al. 2015). In our simulations, bubbles rising in 1-propanol showed negligible variations from rectilinear rise, which is consistent with results of above-mentioned study.

The region of higher Galilei and low Eötvös number (approximately 30 < Ga < 160, Eo < 10) is called asymmetric (Tripathi et al. 2015), characterized by non-rectilinear rise and asymmetric wake of the bubble. This is where our results of rising bubbles in the water belong. However, the sharp boundaries in the phase plot do not reflect the reality, as the bubble behavior in water strongly depends on the purity of used water and the bubble release mechanism. Duineveld (1995) conducted experiments in hyperclean water and determined the boundary of the onset of trajectory instabilities to be at Reb = 660 and We = 3.3. Both our computed and experimental velocities of the 1.5 mm bubble (see Table 2) are below these limit values mentioned by Duineveld (1995). Pesci et al. (2018) mentioned that the lateral movement and path instability of the bubble rising in pure water may be caused by the unstructured computational mesh. Additionally, slight differences in the bubble behavior can be expected, when compared to results of Tripathi et al. (2015), as they were computed with lower viscosity ratio of the air–water system. At approximately 27 °C, their 1 mm diameter bubble reached Ga ~ 50 and Eo ~ 0.15. The degree of deformation is not directly stated in the study. Borkowski and Zawala (2021) recently published a detailed study focused on bubble motion in pure water. They compared the experimental data with the simulation using open software Gerris. Our experimental data, namely the aspect ratio and velocity, match perfectly. The bubble velocity calculated by the Gerris simulation was always slightly (units of percent) lower than experimental. Our results from COMSOL were also underestimated.

Both our experimental results and CFD results of bubbles rising in glycerol were compared with the data from the literature. Sharaf et al. (2017) published a detailed study focused on bubble motion in water-glycerol mixtures; unfortunately, they do not state exact velocity or aspect ratio values. For pure glycerol, they state the viscosity of 1657 mPa s at 30 °C and Morton number of 230. Our value of the Morton number is 50.2, which is due to the lower viscosity (1127 mPa s). Our results for bubbles in pure glycerol with ranges of Eötvös and Galilei numbers of (0.2–0.4) and (0.08–0.13), respectively, belong to the axisymmetric region according to their phase plot. Bubbles stay spherical and show no distinctive variations in the rising path. The agreement between the two studies is clear.

Conclusions

The transient motion of bubbles of different initial diameters rising through three quiescent liquids was studied directly in a fully 3D approach using Level set method implemented in the commercial CFD solver COMSOL Multiphysics. Averaged values of velocities and aspect ratios in steady state were obtained for validation purposes. The aim of the study was to assess the capabilities of a commercial solver in the bubble rising problem. The motion of the bubbles was simulated in a bulk of liquid with low surface tension and low viscosity (1-propanol), high surface tension and low viscosity (pure water), and high surface tension and high viscosity (glycerol). The simulations were performed for initially spherical bubbles to study the different bubble shape deformations simultaneously with the bubble terminal velocities. Trajectories were also investigated to verify the presumed type of motion. These results were compared to experimental data, widely used theoretical correlations, and findings in the literature.

The simulations converged successfully, and the use of conservative Level set method overcame known issues with mass conservation. The results obtained for bubbles rising in 1-propanol agreed with the experimental data, theoretically calculated values, and with the values found in the literature (Eo–Ga plot in (Tripathi et al. 2015)). The resulting velocities and shape deformations for bubbles rising in water agreed with experiment and theory only for the case of a 1 mm bubble; the velocities and shape deformations of other simulated bubbles (1.4 and 1.5 mm) were underestimated. In the case of glycerol, bubbles remained spherical, which agrees with the creeping flow condition mentioned in (Loth 2008) and theoretically calculated shape deformations from Eq. 5. However, the velocities of the bubbles agreed with the calculated values using the Stokes drag law (Eq. 10), which is applicable to rigid spherical particles or fully immobilized bubbles. In our case, the fully mobile bubble surface was considered.

Our results confirm that the larger bubbles, which tend to show non-symmetrical shape deformations, have to be modeled in a fully 3D approach, as the variations from rectilinear rise may appear. COMSOL Multiphysics with the conservative Level set method has been tested. We proved its ability to simulate the two-phase flow systems with large differences in physical properties (viscosity and density ratios), which is known as one of the main causes of numerical errors in CFD solvers (Klostermann et al. 2012). Nevertheless, the settings both of the Level set parameters and the solver should be further studied to ensure the best results.