Oscillation dynamics of a bubble rising in viscous liquid

Abstract

In the present work, we study the oscillation dynamics of 3 mm-diameter bubbles generated through an orifice submerged in viscous liquids. The viscosity of those liquids is varied to change the behavior of the rising bubble. The details of the rising motion and shape oscillation of the bubbles are measured using a combination of high speed, high-resolution imaging, and an accurate digital image processing technique. Direct Numerical Simulations that mimic the experimental conditions are also performed using a front-tracking technique, called the Local Front Reconstruction Method. The predictions of the bubble shape and rising velocity obtained by the numerical simulations show good agreement with the experimental results. Our experimental and numerical results show that the oscillation frequency and the damping rate at lower modes can be predicted using available theoretical models found in the literature. However, discrepancies arise between our results with the theoretical predictions at higher order oscillation modes. We conclude that the discrepancies are due to the influence of rising motion and the vortex wave, which is not considered in the theoretical models.

Graphic abstract

1 Introduction

Because of its relevance in many industrial and natural phenomena, the dynamics of a rising gas bubble has been studied for many centuries and continues to be a problem of significant interest nowadays. This includes for example boiling, flotation, and intensification of heat and mass transfer in bubble column reactors. However, the fundamental understanding of bubble dynamics is still incomplete due to the large number of parameters, the non-linearity, and the fully three-dimensional transient nature of the problem.

The process of bubble detachment and rise is illustrated in Fig. 1. In this case, the bubble is formed from an orifice plate. After the detachment of the bubble from the orifice, the buoyant forces drive the rising motion of the bubble in the liquid. While rising, the bubble undergoes shape deformation until it reaches its equilibrium shape and subsequently oscillates around the equilibrium shape. The characteristics of the rising motion (i.e., the shape, rise velocity, and trajectory) strongly depend on the fluid properties as well as the bubble volume, making it very difficult to develop a general understanding for all the cases.

Fig. 1

Bubble rising in a viscous liquid after its detachment from an orifice plate. The dynamics of the bubble consists of two distinct stages: a transition stage and a terminal stage. The bubble undergoes shape deformation (left) until it reaches its equilibrium shape and subsequently oscillates (right) around the equilibrium shape

Interestingly, Wu and Gharib (2002), Tomiyama et al. (2002), and Laqua et al. (2016) have independently reported that the method of bubble release from the orifice, i.e., the initial condition before detachment, affects the terminal rise velocity and the bubble shape. This unexpected experimental finding is attributed by the authors to the initial shape of the bubble at the beginning of its rise, which is in turn influenced by the way that the bubble is formed. Upon release, a bubble created by a small capillary undergoes strong shape oscillations and is found to reach a higher terminal velocity than a bubble of equal volume released from a large capillary, whose detachment is a much gentler process. This behavior is in strong contrast with the more common assumption of a unique terminal rise velocity for a bubble, which is often used in literature (Clift et al. 1978).

The correlation between rising motion and shape oscillations has been investigated by several authors (Meiron 1989; Lunde and Perkins 1997; De Vries et al. 2002; Veldhuis et al. 2008; Lalanne et al. 2013, 2015). Meiron (1989) developed a numerical model to study the steady rise and stability of inviscid bubbles and showed that the interaction of hydrodynamic pressure and surface tension forces does not lead to linear instability of the bubble path. Lunde and Perkins (1997) performed an experimental study on bubbles rising in still tap water. They showed that the low modes of the shape oscillations account for the wobbly and rocking nature of the shape and motion of intermediately sized bubbles. De Vries et al. (2002) performed experiments with 2–4 mm-diameter bubbles interacting with a hot-film anemometer probe in ultra-clean water to study whether bubble shape oscillation affects the rising velocity. Their experiments showed that the oscillating bubbles do not have higher mean velocities than the non-oscillating bubbles, which is in contrast with the results of Wu and Gharib (2002) and Tomiyama et al. (2002). Veldhuis et al. (2008) investigated the surface oscillations on bubbles rising in water. They showed that shape oscillation is linked to the path instability by checking the frequency of the main mode and the vortex shedding. Recently, Lalanne et al. (2015) studied the influence of the rising movement on the shape oscillations of small bubbles using experiments and numerical simulations. They showed that the effect of the rising motion on the oscillations is negligible, provided that the mean shape of the oscillation remains close to a sphere. Their finding indicates that the shape oscillation and the rising motion is one-way coupled. Lalanne et al. (2013) studied the distinct mechanism of the influence of rising motion on the oscillation of drops and bubbles. The oscillation dynamics taken at the transient stage are studied by DNS simulations, in which drops and bubbles are with initially prescribed deformations.

Limited by the requirements of high spatial and temporal resolutions and sophisticated data treatments of such a highly dynamic physics process, there have been only a few detailed studies. This motivates us to investigate larger bubbles, where the interface dynamics becomes even more pronounced. In the present work, we study the rising motion and shape oscillation of 3 mm bubbles generated through an orifice submerged in viscous liquids. The details of rising motion and shape oscillation of the bubbles are measured using a combination of a high-speed, high-resolution camera and an accurate digital image processing technique. Numerical simulations that mimic the experimental results are also performed using a front-tracking model, the Local Front Reconstruction Method (LFRM). Experiments and simulations have different strengths and weaknesses, related to control and accuracy. Therefore, the correspondence of experimental and simulation results makes the presented results more robust. This paper is organized as follows. The description of the experimental setup and measurement techniques is given in Sect. 2. Section 3 provides a detailed description of the numerical model. In Sect. 4, an extensive analysis of oscillation dynamics of a rising bubble in a viscous liquid is given. Finally, a summary of the main conclusions of the present work is provided.

2 Experimental method

2.1 Experimental setup

The experimental setup is schematically shown in Fig. 2. Air bubbles are formed from a submerged orifice plate made of stainless steel attached to the lower wall of the column. The orifice has an inner diameter of 1 mm. A capillary tube with an inner diameter of 0.5 mm and a length of 150 cm is used to generate sufficient pressure drop that ensures a constant flow condition. The volumetric gas flow rate is controlled using a combination of a kdScientific LEGATO 100 syringe pump and a 2.5 mL Hamilton 1000 series GASTIGHT syringe. The experiments are performed using three different glycerol–water mixtures listed in Table 1. The properties of the liquids are determined using a Brookfield DV-E viscometer for measuring the viscosity and K20 EasyDyne digital of Krusse with the Wilhelmy plate method for measuring the surface tension.

Fig. 2

Schematic diagram of the experimental setup of bubble formation from a submerged orifice and the calibration plate

Table 1 Fluid physical properties for air bubble formation in glycerol–water mixture

2.2 Measurement and image processing techniques

The experimental images are captured using a pco.dimax HD high-speed digital video camera with a frame rate of 2000 Hz and a resolution of $1.5 \times 10^{-2}$ mm/pixel. A calibration plate with known size and distance was used to convert pixel distances to mm. On the plate, the center distance of two horizontal/vertical neighboring white dots is 5 mm and the size of the dots is 1.2 mm (Kong et al. 2018). The spatial resolution is high enough to capture the bubble interface without the help of any interpolative reconstruction. In addition, the selected temporal resolution ensures that the bubble moves 3–4 pixels on two consecutive images. The recordings are performed with the help of a back lighting, which is only switched on during the imaging to reduce the heating of the liquid by illumination of the channel. A calibration plate with specific markers is used to obtain the pixel size.

The captured images are then processed using an in-house Matlab code, which utilizes the image processing toolbox. The main image processing steps are the determination of a threshold for the gray-scale image, binarization of the images using this threshold value, and then determination of the bubble volume, center of gravity, and aspect ratio (Fig. 3). The uncertainty of the determined edge is about 0.5% due to the high resolution of camera. The velocity and deformation of a bubble are derived from the centroid and axis lengths of bubbles, respectively. Furthermore, an advanced smooth spline fitting technique based on smoothing parameter determination theories (De Boor et al. 1978; Wahba 1990; Hurvich et al. 1998; Krakauer and Krakauer 2012) is applied to fit the bubble centroid as a function of time to suppress the noise of velocity calculation. The corrected Akaike information criterion (AICC) approach (Hurvich et al. 1998) is chosen to determine the smoothing parameter. In Fig. 3, the black dot, red lines, and blue lines represent the bubble centroid, major axis lines, and minor axis line of a typical bubble, respectively. The deformation is defined by aspect ratio as in the following equation:

$$\begin{aligned} \chi =l_\mathrm{{major}}/l_\mathrm{{minor}}. \end{aligned}$$

(1)

Fig. 3

Raw image (left), processed image (right). Black dot: center, C; blue line: the minor axis; red line: the major axis

In the present work, even though the sizes are relatively large, the surfaces of bubbles are less dynamic due to the higher viscosity. Therefore, the bubbles are generally axisymmetric. This allows the implementation of calculations of center of 3D and bubble volume using Pappus second theorem (Legendre et al. 2012). It should be noted that the 2D centroid, $C_\mathrm{{2D}}$, computed from area averaging of the planar (projected) bubble image, is not the same as the 3D volume averaged centroid, $C_\mathrm{{3D}}$, especially for bubbles that have a strong shape deformation (Fig. 4). In the image processing, $C_\mathrm{{3D}}$ is computed using the assumption of axisymmetry. This is further confirmed by analyzing a series of binary images obtained from the numerical simulation (Sect. 3). In Fig. 5, the rising velocities calculated using the displacement of $C_{\ rm{2D}}$ and $C_\mathrm{{3D}}$ are compared with the result obtained using the numerical simulation. We use $C_\mathrm{{3D}}$ for the comparison of experiments with simulations, because the simulations are fully 3D and thus the 3D centroid is computed.

Fig. 4

2D/3D center shift of a bubble. The black dot is the center of 2D $C_\mathrm{{2D}}$ and the red dot is the center of 3D $C_\mathrm{{3D}}$

Fig. 5

The bubble rising velocity calculated using two different methods (displacement of 2D $C_\mathrm{{2D}}$ and center of 3D $C_\mathrm{{3D}}$) and numerical simulation

3 Numerical model

The numerical model used in the present work is based on the Local Front Reconstruction Method (LFRM), originally developed by Shin et al. (2011) and adjusted by Mirsandi et al. (2018, 2019). In the sections below, the main characteristics of the numerical model are described.

3.1 Governing equations and solution methodology

In the numerical model, both fluids are assumed to be incompressible, immiscible, and Newtonian. A one fluid formulation is used to describe the fluid flow for both phases. The governing mass and momentum conservation equations are expressed as follows:

$$\begin{aligned} \nabla \cdot \mathbf {u}= & {} 0 \end{aligned}$$

(2)

$$\begin{aligned} \rho \frac{\partial \mathbf {u}}{\partial t}= & {} -\nabla p - \rho \nabla \cdot \left( \mathbf {u}\mathbf {u}\right) - \nabla \cdot \varvec{\tau }+\rho \mathbf {g} + \mathbf {F}_{\sigma }, \end{aligned}$$

(3)

where $\mathbf {u}$ is the fluid velocity, p is the pressure, and $\varvec{\tau }$ is the stress tensor given by $- \mu \big [ \nabla \mathbf {u} + \big (\nabla \mathbf {u}\big )^{T} \big ]$. The local averaged density $\rho$ and dynamic viscosity $\mu$ depend on the local fluid phase distribution and, hence, are calculated from the local phase fraction, F, using normal and harmonic averaging, respectively (Prosperetti 2002). The local volumetric force accounting for the effect of surface tension, $\mathbf {F}_{\sigma }$, is obtained by employing the hybrid Lagrangian–Eulerian formulation representation of Shin et al. (2005) to minimize unphysical parasitic currents in the vicinity of the interface using the following:

$$\begin{aligned} \mathbf {F}_{\sigma } = \sigma \kappa _{H} \nabla F, \end{aligned}$$

(4)

where $\sigma$ is the surface tension coefficient and $\kappa _{H}$ is twice the mean interface curvature field calculated on the Eulerian grid using the information from the Lagrangian interface. Once the fluid flow is calculated, the Lagrangian marker points, which are used to track the interface, are moved using a fourth-order Runge–Kutta time stepping scheme with the locally cubic spline interpolated fluid velocities. Finally, the phase fraction in each Eulerian cell is updated using a geometrical method based on the location of the marker elements (Dijkhuizen et al. 2010).

Due to the advection of the interface, the size of the marker elements changes, decreasing the quality of the interface mesh. To remedy this situation, the LFRM procedure is periodically performed to ensure that each part of the interface is represented with a sufficient resolution. The details of this procedure can be found in the work of Mirsandi et al. (2018).

3.2 Computational setup

The schematic of the computational domain is given in Fig. 6. The air bubble is injected through an orifice of radius $R_\mathrm{{o}}$ submerged in an initially quiescent liquid. The circular orifice is located at the bottom center of the numerical domain and represented with a staircase approximation. The gas injection is assumed to be at constant flow rate Q. The flow in the gas inlet is assumed fully developed laminar and a parabolic inflow velocity is imposed. Initially, a hemispherical bubble is positioned above the orifice with a radius equal to $R_\mathrm{{o}}$. The bubble base is pinned to the orifice during the growth process. At the top of the simulation domain, the pressure-prescribed outlet boundary condition is imposed where the outlet pressure is set to atmospheric pressure, while at the side and lower walls, the no-slip boundary condition is used. The computational domain has a width, length, and height of 5, 5, and 10 equivalent bubble diameters, respectively, to ensure that the bubble formation and rise are not significantly influenced by any wall effect.

Fig. 6

Schematic representation of the computational domain for the simulation of bubble formation and rise in quiescent liquid

For all simulations presented in this paper, a grid size $\Delta = 1.25 \times 10^{-4}$ m and a time step of $\Delta t = 4 \times 10^{-5}$ s are used, for which the results are grid independent. The time step $\Delta t$ is chosen, such that it satisfies both Courant–Friedrichs–Lewy (CFL) and capillary time step restrictions as follows (Brackbill et al. 1992):

$$\begin{aligned} \Delta t< & {} \Delta t_{\mathrm {CFL}} = \frac{\Delta }{v_{\max }} \end{aligned}$$

(5)

$$\begin{aligned} \Delta t< & {} \Delta t_{\sigma } = \root \of {\frac{(\rho _\mathrm{{l}}+\rho _\mathrm{{g}}) \Delta ^3}{4 \pi \sigma }}. \end{aligned}$$

(6)

Here, $v_{\max }$ is the maximum fluid velocity in the computational domain.

4 Results and discussion

4.1 Bubble shape and rising velocity

Table 2 Equivalent diameter of the detached bubble and the corresponding dimensionless numbers for different fluids

Figure 7 shows typical photographs of bubble detachment and rise from the orifice for three different liquids (Table 2). It can be seen that once the bubble detaches, the bubble undergoes shape deformation while rising. The shape evolves from a pendant shape at the time of detachment to an oblate spheroid. This shape deformation becomes more pronounced with decreasing liquid viscosity.

The corresponding bubble shapes obtained using the numerical model, which are represented with red lines, are superimposed on top of the photographs. The spatial evolution of the numerically predicted bubble shapes is very close to the corresponding experimental results. The deviations in the detached bubble diameter in all cases are less than 1%, meaning that the total volume detaching from the orifice plate is roughly the same and unaffected by fluid viscosity.

The comparison of experimental and numerical rising velocity in the case of 60 wt% glycerol mixture is shown in Fig. 8. It is clear that the experimental result agrees well with the simulation result for $Q = 1.1$ ml/min, whereas for $Q = 2.2$ ml/min and $Q = 5.1$ ml/min, the rising velocities are higher in the experiments. The discrepancy is probably due to the wake effect from the preceding bubble, since this influence is nonexistent in the simulation. The comparison of experimental and numerical rising velocity and aspect ratio in the case of $Q = 1.1$ ml/min for three different liquids is shown in Fig. 9. The numerically predicted rising velocity and aspect ratio agree well with the experimental results. This indicates that not only the wake effect but also the influence of surface-active impurities is negligible. For this reason, $Q = 1.1$ ml/min is chosen for further detailed analysis. The corresponding bubbling frequency is 1 Hz.

Fig. 7

Bubble detachment and rising in three different liquids: a 60 wt% glycerol, b 40 wt% glycerol, and c 20 wt% glycerol, for gas flow rate ± 1.1 ml/min. The simulation results are indicated with the red lines. The time interval between the figures is 3 ms

Fig. 8

Comparison of experimental and numerical rising velocities for 60 wt% glycerol at different gas flow rates

Fig. 9

Comparison of experimental and numerical rising velocity and aspect ratio for three different liquids: a 60 wt% glycerol, b 40 wt% glycerol, and c 20 wt% glycerol

The shape evolution with the rising motion is a measure of the interaction of deformation and rising motion. The correlation of the Weber number and the aspect ratio is shown in Fig. 10. Both of the Weber number and the aspect ratio are instantaneous values from experiments. The correlation is compared with correlations developed by Benjamin (1989) based on potential theory, and Legendre et al. (2012) based on fitting of a collection of results taking viscosity into account, giving the following:

$$\begin{aligned} \chi = 1 / \left( 1 - 9 / 64 We \left( 1 + 0.2 Mo^{ 0.1 } \right) ^{ - 1 } \right) . \end{aligned}$$

(7)

Note that the velocity used to calculate the Weber number is the instantaneous velocities, whereas in those two approaches, velocities are terminal velocities. Therefore, the meaningful comparisons are on averaged sense.

From Fig. 10, it is also clear that the results match best for more viscous liquids at moderate We numbers, which are surface tension dominated. For higher We numbers, the difference between the two correlations becomes larger, and also for the current experimental data, more dynamic behavior can be observed.

Fig. 10

Instantaneous aspect ratio against instantaneous Weber number and comparisons with models of Benjamin and Legendre

There are several theoretical models available to predict the frequency of shape oscillation. Lamb (1932) gives an analytical expression for the frequency of shape oscillations of spherical bubbles as follows:

$$\begin{aligned} f_n=\frac{1}{2\pi }\sqrt{\frac{(n+1)(n-1)(n+2)\sigma }{\rho _{l} R^3}}, \end{aligned}$$

(8)

where n is the mode number. A spherical bubble is an isotropic body, and hence, all modes of the same order have the same frequency. For mode 2 oscillations, the frequency is given by the following:

$$\begin{aligned} f_2=\frac{1}{2\pi }\sqrt{\frac{12\sigma }{\rho _{l} R^3}}. \end{aligned}$$

(9)

A bubble deformed into a spheroid is no longer isotropic and modes of the same order, but different degree will generally have different frequencies (Lunde and Perkins 1997). In addition, the influence of rising motion is not taken into account in the Lamb expression. A model to predict the shape oscillation frequencies for non-spherical rising bubbles was derived in a study by Meiron using stability analysis based on potential flow theory (Meiron 1989). On the other hand, Lunde and Perkins (1997) proposed a simple model based on an assumption of a plane capillary wave traveling along the bubble surface, where the frequency is determined by the wave speed and the wave length. For mode 2,0 shape oscillation, the frequency is given as follows:

$$\begin{aligned} f_{2,0}=\frac{1}{2\pi }\sqrt{\frac{16\sqrt{2}\chi ^2}{(\chi ^2+1)^{1.5}}}\sqrt{\frac{\sigma }{\rho _\mathrm{{l}} R^3}}, \end{aligned}$$

(10)

where $\chi$ is the aspect ratio.

The oscillation frequencies obtained from the experiments and numerical simulations are compared with these models in Table 3. The frequencies are extracted through the peaks of curves readily. It is clear that the present results agree well with the models proposed by Meiron and Lunde, but disagree with the spherical model of Lamb. This is further confirmed by comparing the normalized frequency against deformation in Fig. 11. Surprisingly, Meiron’s model based on the potential flow theory is still valid, even though it does not take into account the influences of wake. Although the assumption of a plane wave in Lunde and Perkins should fail for non-spherical bubbles as pointed out by van Wijngaarden and Veldhuis (2008), it coincidentally works well. Similarly, Lalanne et al. (2013, 2015) also showed that the predictions obtained using a potential flow theory agree well with the DNS results for a slightly deformed spherical bubble. This indicates that the wake and viscosity effects are minor, which needs to be further investigated.

Table 3 Comparison of oscillation frequencies for three different fluids

Fig. 11

Normalized frequency against bubble aspect ratio

4.2 Analysis and discussion

In this section, the shape oscillations are studied in detail by decomposing the interface obtained from both experimental and numerical results on the basis of the Legendre polynomials given as follows: (Becker et al. 1991):

$$\begin{aligned} R(\theta ,t)/R_0=a_0(t)+\sum _{l=2}^\infty a_\mathrm{l}(t)P_\mathrm{l}(\cos \theta ), \end{aligned}$$

(11)

where $a_{0}$ is unity in general due to the assumption that the bubble volume is constant and $a_\mathrm{{l}}$ is the amplitude of the mode. The series is truncated at $l = 20$, which ensures an accurate description of the present bubble shapes. The evolution of amplitudes for orders $l =$ 2–5 in three different liquids is shown in Fig. 12. A satisfactory agreement between the experimental and numerical results with respect to the evolution of various harmonics is obtained, i.e., the amplitudes of modes 2–5 of bubbles are similar. In general, the amplitudes are much lower for higher modes, which is more pronounced for cases with a higher liquid viscosity. Moreover, Fig. 12 shows that in higher viscosity liquids, the amplitude of the shape oscillations dampens faster. In the low viscosity liquid, however, the damping rate gradually weakens, and for certain higher modes, the amplitude increases. It is important to note that the modes are not rigorous eigenmodes due to non-spherical steady-state shape of the rising bubble, while the spherical harmonic solutions assumes oscillatory perturbation on a spherical surface. The eigenmodes shift with the evolving shape. However, the decomposed modes still can be used to quantify the effect of rising motion on the shape oscillation (Lalanne et al. 2013).

Fig. 12

Time evolution of spherical harmonic amplitudes $a_{l}$ for $l = 2$-5 for three different liquids: a 60 wt% glycerol, b 40 wt% glycerol, and c 20 wt% glycerol

Figure 12 also reveals that mode 2 is the main mode. The $a_2$ is basically the representative of main shape of bubble and a similar to the aspect ratio in that it describes the extent of shape deformation. The value starts from a positive value due to the prolateness of the initial shape, while it decreases gradually to negative when the shape evolves to an oblate shape. For all three set of comparisons, the results of experiments and simulation match very well on mode 2, which was also evident in the qualitative visualization in Fig. 7. The match deteriorates for higher modes and longer time. This might be attributed to the spatial resolution of the simulation and the inherent in-stable nature of the higher modes. Furthermore, clearly, all the curves of the modes show that the deformations consist of a main shape evolution and a dampened oscillation. The later one is of our interest.

To investigate the oscillation energy decay of these modes, the following window-moving average decomposition is applied:

$$\begin{aligned} Y=\bar{Y}+\tilde{Y}, \end{aligned}$$

(12)

where Y is the instantaneous variable, $\bar{Y}$ is the average value, and $\tilde{Y}$ is the fluctuation. Then, the instantaneous frequency can be calculated by the time interval of peaks, $f=1/(t_{i+1}-t_{i})$ and the damping factors are calculated by fitting the peaks with an exponential damping function $y= ae^{-bt}$. Note that the shape information at the instant after the detachment is excluded from calculation. These mode fluctuations are shown in Fig. 13.

Fig. 13

Fluctuations of a mode 2 and b modes 3–5 for three different liquids

Prosperetti (1980) derived a model by an eigenvalue analysis of spectrum of the vortex form of Naiver–Stokes equation neglecting gravity and buoyancy forces. A general equation that describes the damping of oscillation can be reduced to Lamb’s equation if the surrounding liquid is less viscous (${Oh} \ll 0.1$):

$$\begin{aligned} b_n=(n+2)(2n+1)\mu /(\rho _{l} R^2). \end{aligned}$$

(13)

Then, a dimensionless theoretical damping factor for each mode is given by the following:

$$\begin{aligned} b^*_n=(n+2)(2n+1)Oh. \end{aligned}$$

(14)

The dimensionless damping factors extracted from experimental results (Fig. 13) and the theoretical values are shown in Fig. 14 (left). The damping factors generally increase with increasing Oh number. The oscillation energy needs less time to decay in more viscous liquids. Moreover, it shows that for mode 2, the match between our study and the theoretical equation proposed by Prosperetti (Eq. 14) is satisfactory. However, our study shows that higher modes need more time to decay in contrast to the prediction from the theory. This is mostly due to the mode shift. The weight of high modes that contributes to the total fluctuation becomes more significant gradually with the more and more flattened shape. Note that for the case of ${Oh} = 5 \times 10^{-3}$, the damping factors for mode 4 and above are negative, i.e., the energy no longer decays (not shown). Simulation studies of Lalanne et al. (2013) also revealed this departure of higher modes from the theoretical study (Fig. 14). Moreover, they also found a negative damping factor on mode 4 as well.

Fig. 14

The dimensionless damping factor (left) and frequency (right) against Ohnesorge number. Mode 2: red circle; Mode 3: green filled square; Mode 4: blue circle; Mode 5: deep red diamond. Theoretical prediction (Prosperetti 1980) for Mode 2: red solid line; Theoretical prediction for Mode 3: green dash line

Moreover, we observed that most of the surface oscillation energy (mainly mode 2) damps only for a short period of time. In the case of 60% glycerol solution, the oscillations are damped after $t=0.08$ s. In contrast, for bubbles rising in 20% glycerol solution, path instability occurs after $t=0.08$ s, producing an increase of the oscillation amplitude. On the other hand, in the case of 40% glycerol solution, for each mode, the oscillations maintain their amplitude. Here, our results differ from the theoretical results of Prosperetti (1980) and simulation studies of Lalanne et al. (2013). They reported that the oscillation damping of bubbles is close to the theoretical prediction at the initial stages (or transient stages). In other words, the oscillation behavior of the rising bubble at the initial stages is similar to the oscillation of bubble in zero gravity, i.e., all modes dampen gradually. However, in our study, the oscillations are not damped at lower Oh numbers for any of the modes after the damping at the initial stage of bubble rise, which has also been reported by Lunde and Perkins (1997). Similarly, in the study of Gordillo et al. (2012), the agreement between the analytical model and DNS simulations deteriorates for longer times. The authors attribute this to the viscous dissipation including the boundary layer and wake, which are not considered in the analytical model.

Because the initial surface energies are close for all the cases (see Fig. 13), any oscillation energy should eventually dampen if the system is isolated. In the case of low Re number (${Oh}=3.1 \times 10^{-2}$), there is no strong wake effect beneath the bubble and the surface energy is eventually damped. Whereas in bubbles with higher Re numbers (${Oh}=1.1 \times 10^{-2}$, $0.5 \times 10^{-2}$), there are standing vortices and vortex shedding appears in the wakes. In addition, our results also suggest that this energy is mainly added to higher modes instead of mode 2 at the initial stage of rising. Later, mode 2 also gains energy, which makes the oscillation system at all modes a forced oscillation system. Therefore, the oscillation behavior of bubbles at lower Oh numbers are different from bubbles in zero gravity in our results.

The role of rising motion, or more specific the wake, in shape oscillation is still unclear. Nevertheless, our inference is supported by extra measurements of two different size of bubbles rising in 40% glycerol solution (Fig. 15). The diameters of the bubbles are 2.63 mm and 3.22 mm, respectively, and the Oh numbers of these bubbles are almost identical. However, according to the study of Blanco and Magnaudet (1995), the wakes of these bubbles are of two type: one has no vortex wake, whereas the later one has a vortex wake. We have noticed that for the simulations conducted in the study of Lalanne et al. (2013) vortex wakes would develop in none of those cases. The properties of the most oscillating bubble in their studies are close to the bubble of 2.63 mm in the present study. It is clear from the figure that the shape oscillation is damped at the initial stage for the one without vortex wake, while the oscillating deformation continues for the one with vortex wake even after the initial stage of rising. Veldhuis et al. (2008) reported that bubbles with vortex shedding in wakes have a much more dynamic interface, in which the vortex shedding time scale is larger than the typical oscillation damping time scale. This suggests that if the time scale of a sudden move of the bubble or a change of surrounding flow is shorter than the time scale of the damping, the oscillation behavior will be similar to the one in zero gravity. Otherwise, the oscillations gain external energy. This phenomenon could be more pronounced when bubbles interact with turbulence, which includes a range of varied size eddies.

Fig. 15

The rising velocity and aspect ratio of 2.63 mm and 3.22 mm bubbles in 40% glycerol solution

The dimensionless frequency is shown in Fig. 14 (right). It can be seen that the frequency increases with Oh number. The trend of mode 2 is consistent with Fig. 11, which correlates the frequency with the terminal shape of the bubble. Because the aspect ratio is a result of fluid properties and hydrodynamics, it is meaningful to correlate the dimensionless frequency with Oh number.

5 Conclusions

In the present work, we study the oscillation dynamics of 3 mm-diameter bubbles generated through an orifice submerged in viscous liquids by means of experiments and detailed numerical simulations. First, our results confirm that the potential flow theory can still predict the oscillation frequency of large size bubbles, where the dynamics are pronounced. This finding is also in accordance with several experimental studies reported in the literature (Lunde and Perkins 1997; van Wijngaarden and Veldhuis 2008; Gordillo et al. 2012; Lalanne et al. 2015). Second, it is found that the damping factor at mode 2 can still be described by a theoretical model (Prosperetti 1980), which takes into account viscosity effects and neglects the rising motion effect. However, at higher modes, the damping factor cannot be explained by the theoretical model and so far no available theoretical model can be used to interpret the damping behavior. Wake effects seem to play a very important role as an external energy source driving the shape oscillations of rising bubbles, particularly for the higher modes.