The impact and bounce of air bubbles at a flat fluid interface

Abstract

The rise and impact of bubbles at an initially flat but deformable liquid–air interface in ultraclean liquid systems are modelled by taking into account the buoyancy force, hydrodynamic drag, inertial added mass effect and drainage of the thin film between the bubble and the interface. The bubble–surface interaction is analyzed using lubrication theory that allows for both bubble and surface deformation under a balance of normal stresses and surface tension as well as the long-range nature of deformation along the interface. The quantitative result for collision and bounce is sensitive to the impact velocity of the rising bubble. This velocity is controlled by the combined effects of interfacial tension via the Young–Laplace equation and hydrodynamic stress on the surface, which determine the deformation of the bubble. The drag force that arises from the hydrodynamic stress in turn depends on the hydrodynamic boundary conditions on the bubble surface and its shape. These interrelated factors are accounted for in a consistent manner. The model can predict the rise velocity and shape of millimeter-size bubbles in ultra-clean water, in two silicone oils of different densities and viscosities and in ethanol without any adjustable parameters. The collision and bounce of such bubbles with a flat water/air, silicone oil/air and ethanol/air interface can then be predicted with excellent agreement when compared to experimental observations.

---

Introduction

When a bubble approaches an initially flat free surface with a large enough velocity, a film of liquid is trapped, creating a pressure build-up between the bubble and the free surface. Lubrication and deformation forces can cause a rebound of the bubble prior to film rupture. Advances in high-speed photography allow for very precise measurements of the rise velocity, impact and bounce of bubbles against solid surfaces,^1,2 soft deformable surfaces^3–6 and also from compound films.^7,8 The liquid film eventually breaks and the bubble bursts through the free surface causing small droplets of the liquid phase to be propelled into the air and smaller bubbles can also form in the liquid phase.⁹ This bursting phenomenon occurs on a milliseconds timescale whereas the bounce and drainage processes usually take place over much longer time scales of up to seconds. Thus the rebound and drainage of the film between the bubble and the free surface are the important rate determining steps in the dynamics of a bubble-interface encounter.

Modelling the impact of bubbles against deformable surfaces is challenging because it requires tracking the deformation of the bubble and the deformable surface simultaneously. It also requires modelling the detailed thin film drainage that occurs at separations on micron scale as well as the motion of the bubble of millimetre size to over say a centimetre length scale.

The velocity of the bubble rising in a straight path is predicted using a force balance method where drag, buoyancy, added mass and film forces are taken into account. This approach is an extension of a theory that has been applied to model the experimental data^10,11 on the collision between bubbles and solid surfaces in ultraclean¹² and contaminated^13,14 systems. Advantages of this model include its efficiency and ease of implementation. In contrast, numerical solutions of the full Navier–Stokes equations using, for example, the volume of fluid method¹⁵ require refined grids and are more expensive computationally, even though such approaches can provide more details about the complete flow field.

We use lubrication theory to treat the small length scale film drainage stage in which the interaction between the colliding interfaces cause them to deform. Such an approach was proven to be accurate when the separation between the interfaces became much smaller than the interaction region as is the case for interactions involving bubbles and drops^16,17 approaching at low speeds (low Reynolds numbers) and also bubbles approaching a solid surface at high speeds.^18–20

A quantitative account of the impact and bounce of a bubble at an initially flat but deformable interface requires precise knowledge of the shape and velocity of the bubble as it rises in bulk as an initial condition. In itself, the topic of terminal velocity and shape of rising bubbles in clean liquid systems has been an important area of research both experimentally^21–23 and theoretically.^24–26 Recently, research is focused on the onset of instabilities, when the bubble rise trajectory changes from a rectilinear to a zig-zag or spiral path.^27,28 In this article, we will only be concerned with rectilinear trajectory regimes that are axisymmetric during bubble rise and its subsequent impact with the surface. We consider clean systems in which the boundary condition at the bubble surface is mobile, that is, the tangential stress on the bubble vanishes.

However, it has been observed experimentally that the initial few bounces of a rising bubble in water at a solid hydrophilic (water-wet) surface such as a smooth glass plate,¹ at a solid hydrophobic (non-water-wet) surface such as a slightly rough Teflon plate² or at a deformable water/air interface, are all very similar^3,4 even though the final states of the bubble are obviously quite different. This leads to the conclusion that the hydrodynamic boundary conditions at all three surfaces must be very similar, if not identical. For our modelling we therefore consider the fluid boundary condition of the free surface to be tangentially immobile, that is, the fluid velocity at the free air/water interface follows the velocity of the moving and deforming interface rather than obeying the condition of zero tangential stress.

The manuscript is organized as follows. A prerequisite for a quantitative description of a bubble bouncing at a solid or deformable surface is the correct prediction of the approach velocity. If this velocity is too low or too high, the extent of the rebound is under- or over-estimated. Therefore we will first present a model to predict this approach velocity. This model also needs to be able to predict the correct deformation of the shape of the rising bubble in a self-consistent way. Such preliminary considerations form the starting point for constructing a predictive model of the interaction between a bubble and a free surface. We then make detailed comparisons with results of bubble bounce experiments in water, in ethanol and in two silicone oils of different densities and viscosities at the free air/liquid interface for a range of bubble radii. We also use our model to predict characteristic physical quantities that are yet to be measured experimentally.

Stokes–Reynolds–Young–Laplace equations

We present the equations used to describe the deformation of the initially flat interface and the bubble. Let us consider a bubble of initial radius R that rises with a time dependent velocity V(t) in a liquid with density ρ and viscosity μ. The axisymmetric bubble impacts at a deformable air–liquid surface z(r,t) as shown in Fig. 1 where t indicates time and r the radial coordinate. A thin liquid film h(r,t) forms between the bubble and the surface. The initially flat undeformed liquid surface defines the plane z = 0 of the coordinate system and the axis of symmetry corresponds to r = 0.

Fig. 1 Schematic of a bubble with equivalent radius R rising with the velocity of the centre of mass V(t) impacting on a deformable surface. The axisymmetric shape of the surface z(r,t) and the bubble z_b(r,t) as well as the film thickness or separation h(r,t) between the bubble and the surface are indicated. The plane z = 0 locates the undeformed surface prior to bubble impact.

The density of air is negligible compared to the density of the liquid ρ and will not be considered. The axisymmetric shape, z(r,t) of the initially flat interface obeys the following equation, with the assumption that the slope of the deformation is small: |dz/dr| ≪ 1

where the left hand side represents surface tension σ times the curvature and on the right hand side, g = 9.82 m s⁻² is the acceleration due to gravity, p is the pressure due to fluid motion and Π(h) is the disjoining pressure due to surface forces that are functions of the film thickness, h between the bubble and the surface.

The shape of the top of the bubble z_b(r,t) is given by

with the assumption |dz_b/dr| ≪ 1. The term 2σ/R represents the Laplace pressure of the bubble. Eqn (1) and (2) describe the shape of the free surface and the bubble. The equation for the film thickness h = z − z_b is found by combining eqn (1) and (2)Eqn (3) is critical for our theory since surface tension is the only element capable of storing energy. During rebound, the stored energy is given back. Energy loss occurs due to viscous losses in the film and due to the drag force. As we shall see during the initial approach of the bubble towards the surface, the film is relatively thick so the disjoining pressure Π is small and can be neglected during the initial approach.

The Reynolds number within the film is small and the drainage therefore obeys Stokes flow. The velocity profile in this thin film region is parabolic and the drainage process is calculated from the classical Stokes–Reynolds equation. Assuming that the immobile boundary condition holds at the deformable flat surface whereas the bubble surface has a mobile boundary condition we have the following film thinning equation^12,17

If both surfaces are immobile factor 3 should be replaced by 12. When both surfaces are mobile, the flow in the film is no longer parabolic, but instead becomes a plug flow. In that case we cannot use lubrication theory because the flow in the film is inertia controlled. Chesters and Hofman²⁹ proposed a model to deal with this situation. Their model has the disadvantage that the pressure does not decay fast enough, resulting in a non-converging film force.¹⁷

Analytical shape of the surface under an applied force

To obtain a scale of the deformation of the free surface, we first consider the equilibrium deformation of the air–liquid interface due to the presence of a bubble of radius R pushing the surface upward from below as shown in Fig. 1. The equilibrium state is established by the balance between the buoyancy force and the surface force due to deformation of the free surface. A repulsive disjoining pressure between the free surface and the bubble mediates this interaction. For the purpose of the following derivation, the detailed form of the disjoining pressure, Π(h), is unimportant as long as it is repulsive. At equilibrium the hydrodynamic pressure vanishes, p = 0. The analysis is based on matching the deformation of the initially flat surface in the outer region far from the bubble (r ≫ R), with the solution in the interaction region (r ≤ R).

Outer region of the free surface

In the outer region r ≫ R, h is large, so Π(h) can be neglected and eqn (1) becomes a Bessel equationwhich has an analytical solution, in terms of the modified Bessel function of the second kind of order zero³⁰where the capillary length, λ, is defined as The constant A will be found by matching with the inner solution. To do so, we note that the asymptotic form of eqn (6) when r ≪ λ iswhere γ_E = 0.57721566 is the Euler constant. More details on the asymptotic behaviour of K₀ are given in the Appendix.

To obtain the constant A in eqn (6), the outer solution has to be matched with the inner solution for the film where the disjoining pressure, Π is important but the term ρgz can be neglected whereby eqn (1) becomes

The first integration of eqn (9) yieldssince axial symmetry requires dz/dr = 0 at r = 0. The second integration yields in the limit r → ∞whereis the force between the bubble and the deformed surface. Matching the coefficients of the logarithmic terms in eqn (8) and (11) gives A = F/(2πσ) so that the shape of the deformable surface outside the interaction zone between the bubble and the surface is given by This is the result that will be used as a boundary condition for the numerical calculations in which the shape of the horizontal surface is calculated using eqn (1) and (4) with boundary conditions dz/dr = 0 at r = 0 and z(r_m) = (F/2πσ)K₀(r_m/λ), with r_m being some large radial value at which the disjoining pressure has become negligibly small, but still satisfies r_m < λ, as illustrated in Fig. 2.

Fig. 2 The analytical solution, eqn (16)–(18), for the equilibrium deformation of the free interface (solid line) due to buoyancy force on a bubble resting beneath it in (a) water (R = 0.74 mm) and in (b) ethanol (R = 0.81 mm). The bubble profile (dashed line) was calculated numerically. Note that the axis on the left and right side of the figures are different in order to show the different length scales more clearly; i.e. bubble radius R, capillary length λ and the extent of the film region r₀ (eqn (17)).

Inner region of the surface

We can obtain a complete approximate analytical solution by deriving an expression for the inner region where the bubble and the surface interact. We assume that a repulsive disjoining pressure maintains a thin equilibrium film between the bubble and the deformed free surface that has a near uniform thickness, (h ∼ constant). Then from eqn (3), we can make the approximation: Π ∼ σ/R when ρgz is small, so that eqn (1) becomes Integrating this twice together with the axisymmetric condition dz/dr = 0 at r = 0 yields the inner solutionwhere z₀ = z(r = 0).

Finally, the complete approximate analytical solution for the equilibrium shape of the interface is the combination of eqn (13) and (15)

To find the constant r₀ where the two functional forms change over, we equate the absolute value of the force F with the force due to the approximate disjoining pressure that has the form, Π = σ/R, 0 < r < r₀, and Π = 0, r > r₀, acting over the film of area πr₀². This gives The constant z₀ can be found by equating the two forms of the solution in eqn (16) at r₀ to give This is just the maximum equilibrium central deformation of the free surface due to the presence of the bubble beneath it being pushed up by buoyancy. Furthermore, it can be shown that the derivative of z(r) is also continuous at r₀.

In Fig. 2a and b, we show the approximate analytic solution given by eqn (16)–(18) for bubbles of radius R = 0.74 mm in water and R = 0.81 mm in ethanol. These two radii correspond to the size of bouncing bubbles that will be investigated later. We take F = 4πR³ρg/3 as the absolute value of the buoyancy force on a spherical bubble of radius R. It is assumed that a stable film is maintained between the bubble and the free surface by a repulsive disjoining pressure and a final equilibrium configuration is reached when the surface tension force of the deformation of the free surface and the buoyancy force of the bubble balance each other. The radius of curvature at the apex of the free surface, see eqn (15), is exactly twice the original radius the bubble.

In this approximate analysis for each liquid, the disjoining pressure in the film is taken to be Π ≈ σ/R for 0 < r < r₀ and it is zero for r > r₀.

Prediction of the approach velocity

To model the bouncing behaviour accurately, it is essential to start with a model that can predict the approach speed correctly. In this section we describe the equation of motion of a rising bubble with surface tension σ, which started from rest to reach a constant terminal velocity V_T in a medium with viscosity μ and density ρ. As the velocity of the bubble increases, it deforms from a sphere into an approximate oblate ellipsoid due to inertial effects. In the Appendix we derive an approximate analytical relation for the aspect ratio of the deformed bubble rising in ultraclean systems due to the balance of inertia and surface tension as characterised by the Weber number, We = 2RρV_T²/σ. In the inset of Fig. 3b we show the ellipsoidal bubble as well as the spherical bubble with equivalent radius R, which is defined to have the same volume as the deformed bubble through R³ = R_h²R_v, where R_h and R_v are the horizontal and vertical radii of the deformed bubble.

Fig. 3 Comparison between experiments of Zawala et al.,^3,5 Suñol et al.,⁶ Duineveld²¹ and Wu and Gharib²³ (symbols) with the theoretical prediction (lines) using Moore's theory, eqn (19)–(24), for (a) terminal velocity, V_T of bubbles rising in silicone oil, ethanol and ultrapure water. The dashed line is the variation for water assuming the bubble remained spherical. (b) Vertical-to-horizontal aspect ratio of bubbles in water and ethanol as a function of the Weber number, We = 2RρV_T²/σ according to eqn (24). The inset shows definitions of vertical and horizontal radii.

Bubble approach velocity: buoyancy vs. drag

Buoyancy force will cause a bubble in a liquid to rise. This force points in the direction opposite to the acceleration vector due to gravity and is given by The bubble attains a constant approach velocity when the buoyancy force balances the hydrodynamic drag force²⁴ The drag force is characterised in terms of the drag coefficient, C_D and the instantaneous Reynolds number, Re = 2Rρ|V|/μ. The product C_D Re is given in terms of the aspect ratio χ of the oblate ellipsoid that approximates the shape of the deformed bubble according to the theory of Moore²⁴and The function K in eqn (21) was tabulated by Moore,²⁴ but will be approximated by the following polynomial;¹³

K(χ) = 0.0195χ⁴ − 0.2134χ³ + 1.7026χ² − 2.1461χ − 1.5732 (23)

A relation between the aspect ratio χ and the Weber number, We = 2RρV_T²/σ, is derived in the Appendix as This completes the approximate theory for the drag force on rising bubbles in bulk liquid. Eqn (19)–(24) give the hydrodynamic drag force, F_D on a bubble as a function of its velocity, V that accounts consistently for the deformation of the spherical bubble with the initial radius, R into an oblate ellipsoid with the aspect ratio χ due to inertial forces. This derivation assumes that the shear stress vanishes on the surface of the bubble and is therefore applicable for bubbles in ultra-clean liquids.

Comparison with experiments: terminal velocity and shape

We consider the experimental data for the terminal velocity of bubbles rising in bulk in ultraclean water, silicone oil and ethanol in which bouncing bubble experiments have been performed. The physical parameters for these systems are presented in Table 1.

Table 1 Physical parameters for water,⁴ silicone oil⁵ and ethanol⁶ used in the numerical calculations

Parameter	Water	Oil A	Oil B	Ethanol
Density (kg m^−3) ρ	1000	750	850	789
Viscosity (mPa s) μ	1.0	0.5	1.3	1.2
Interfacial tension (mN m^−1) σ	72	16	17	22

---

As shown in Fig. 3a, the terminal velocity as a function of bubble size predicted using eqn (19)–(24) compares well with the experimental data from the literature.^3,6,21,23 The terminal velocity of bubbles of the same radius in ethanol is considerably lower than that in water due to smaller buoyancy, larger viscosity and especially larger deformation caused by the lower interfacial tension (see Table 1). Results for silicone oil B (not shown in Fig. 3a) overlap those for ethanol due to their very similar liquid properties.

In Fig. 3b, we show the experimental data from various sources^6,21,22 for the variation of the aspect ratio, χ, of the deformed rising bubble with bubble size expressed in terms of the Weber number, We. Comparison with eqn (24) shows a very good agreement for a wide range of Weber numbers, even though the derivation of this result assumes small deformations (see appendix). There is a critical inverse aspect ratio of about 1/χ ∼ 0.5, corresponding to We ∼ 3, below which bubbles no longer rise along a straight path, but follow a zig-zag or spiral path.^21,27 This transition is evident from the scatter of data points at around We ∼ 3.

In the next section, we use the bubble velocity calculated from the balance of drag and buoyancy forces as the initial condition to the rising bubble before it impacts the free surface.

Bubble interacting with a soft deformable surface

When a bubble decelerates (here due to collision with a soft deformable surface) the surrounding fluid must also be decelerated. This will give rise to an added mass force:¹⁴where C_m = 0.5 is the added mass coefficient for a spherical bubble in bulk liquid. Miloh³¹ has shown that the added mass coefficient for a spherical bubble touching a flat free surface is 0.4198. However, in the current situation, the free surface also moves upwards and it seems to be justified to assume that the added mass coefficient has a value very close to C_m = 0.5.

When the bubble approaches the surface, lubrication theory provides a relationship between the pressure, p in between the two surfaces and the separation h according to eqn (4). The lubrication pressure, p, in the thin liquid film builds up and generates a film force that can be found by integrating p over the axisymmetric film region:

where r_m is the domain size for the numerical computation at which p has essentially decreased to zero.

The velocity V of the centre of mass is obtained by equating all forces acting on the bubble, which results in the following equation of motion for the bubble as it rises and impacts the surface

F_A + F_B + F_D + F_F = 0. (27)

Using eqn (19), (20), (25) and (26) we obtain a point force model for the centre of mass of the bubbleEqn (3) (with Π = 0) and (4) constitute the partial differential system to be solved numerically for the film thickness h(r,t). Apart from using the terminal velocity as the initial velocity, the initial film thickness is given bywhere H₀₀ is the initial distance between the top of the bubble and the z = 0 plane that defines the undeformed free surface.

To complete the formulation we need to provide four boundary conditions. The outer solution for the free surface given by eqn (6): z(r) = AK₀(r/λ) provides an analytical expression z(r,t) that is valid for large r. The numerical solution for the inner solution of the shape of the deformable surface will be matched to this analytical solution. To derive a boundary condition at the large radial position r_m, we use eqn (13) and (29) to write the separation h at r_m

where the force acting on the surface is now the film force, that is F = F_F, given by eqn (26). Taking the time derivative of eqn (30) and assuming dH₀/dt = −V, the boundary condition at the outside border r_m is given by Furthermore, we take p = 0 at r = r_m. At the axis of symmetry (r = 0), dp/dr = 0 and dh/dr = 0. The force F_F on the bubble is computed from eqn (26) using Simpson's rule. It is essential to apply this correct boundary condition at r_m to obtain results that are independent of the domain size. However, the computational domain must be large enough to be able to describe the drainage process completely. Here we take r_m = 1.2R. The system of equations with the above boundary conditions is solved using a standard differential algebraic solver in Matlab. Note that the constructed model is free of any fitting parameters.

Comparison with bouncing experiments

The experimental data chosen for comparison were performed for bubbles in ultraclean water,^3,4 silicone oil⁵ and ethanol.⁶ The characteristic bouncing behaviour depends strongly on the size and the approach velocity of the bubble. Each bubble was taken to be ellipsoidal during rise then its shape was changed to a sphere after the first impact when the velocity becomes zero for the first time, that is, setting the aspect ratio χ = 1 after the first bounce. We will compare predictions of our model over a wide range of radii and approach velocities obtained experimentally by releasing the bubble at different initial separations from the surface.

In Fig. 4a, two bubbles with different radii are released sufficiently far away from the surface to reach a constant terminal velocity before impact. In this case, we assume that the bubbles become spherical after the first impact around 10 ms as observed experimentally. The theory agrees well with the experimental data of Zawala et al.³ until film rupture occurred in the experiments at roughly 32 ms and 80 ms for these bubbles. We make no attempt to predict the coalescence time as we do not have detailed information about the surface chemistry that is responsible for the development of the disjoining pressure.

Fig. 4 Comparison between theory (lines) and experiment (symbols) for the rise and impact of bubbles in ultraclean water. (a) Experiments of Zawala et al.,³ in which two bubbles with different radii (R = 0.74 and 0.50 mm) are released 250 mm away from the surface and impact with an approach speed of 34.5 and 27.8 cm s⁻¹. (b) Experiments of Kosior et al.⁴ in which a bubble (R = 0.74 mm) is released from the tip of a syringe placed 3 mm away from the deformable surface. (c) Comparison between the experimental data from (a) and (b) and the data from (b) shifted in time to match at the first bounce.

In Fig. 4b we show a comparison with experimental bubble velocity obtained by Kosior et al.⁴ for a bubble in ultraclean water (R = 0.74 mm) released from the tip of a syringe placed 3 mm away from the air–water interface, therefore the initial separation between the top of the bubble and the surface H₀₀ ∼ 1.52 mm is very close to the value H₀₀ = 1.38 mm used in the model. Thus the result in Fig. 4b is for a bubble of the same size as in Fig. 4a but released close to the interface and therefore had not yet attained terminal velocity. The agreement is impressive. In this experiment, the film was observed to rupture after about 55 ms from release. However, note that no experimental values were reported for t < 10 ms. The drag was calculated assuming the bubble remained spherical (χ = 1) for the entire collision and bounce process. The initial acceleration of a bubble in bulk released with zero velocity will be 2g as a consequence of the added mass coefficient C_m = 0.5 in eqn (25).

In Fig. 4c, we overlay the experimental results for R = 0.74 mm of Fig. 4a (bubble approaching with interface at terminal velocity) and the results of Fig. 4b, but shifted by 20 ms so that the first peaks of the two data sets coincide. Since the theoretical curve in Fig. 4a now matches almost exactly both sets of data, this suggests that any possible flow disturbance caused by the first bounce has little influence on subsequent bounces.

Comparison with experiments of Zawala et al.⁵ for air bubbles impacting an air–oil interface in two different silicone oils (see Table 1 for properties) is shown in Fig. 5. In Fig. 5a, two bubbles of R = 0.53 and 0.32 mm are released in silicone oil A (see Table 1) while in Fig. 5b two bubbles of similar radii of R = 0.58 and 0.33 mm are released in silicone oil B. For this system, our model also predicts the dynamic behaviour of the bubble free surface system very well. Both the amplitudes of the velocity fluctuations and their timing are predicted correctly in our model without fitting parameters. Only two representative experimental bubble sizes of Zawala et al.⁵ are shown here, the largest and the smallest bubble radii. The cases not shown also show excellent agreement.

Fig. 5 Comparison between theory (lines) and experimental data (symbols) of Zawala et al.⁵ for the impact and bounce or air bubbles in two different silicone oils.

As a last experimental system we investigate bubbles in ethanol. Comparison between our theory and experiments of Suñol et al.⁶ for bubbles with different radii rising in ethanol is shown in Fig. 6. An arrow indicates the time when each film ruptured and the bubble coalesced with the air above the free interface.

Fig. 6 Comparison between theory (lines) and experiment (symbols) from Suñol et al.⁶ for bubbles in ethanol. (a) Position of the centre of mass of the bubble relative to the originally flat free surface. Note that the curve for R = 0.21 mm has been shifted to the right for the figure not to be too cluttered. (b) Velocity of the centre of mass of bubbles with different radii rising in ethanol and impacting with a deformable air–ethanol surface. The arrows indicate the time when film rupture is observed experimentally.

In Fig. 6a we show the position of the centre of mass of the bubble as a function of time. The centre of mass travels above the z = 0 line for the largest bubble (R = 0.81 mm) at t = 22 ms. The maximum deformation of the surface decreases as the bubble size decreases and similarly for the amplitudes of the bounces, which become smaller with smaller bubble size. No bounce was observed for the smallest bubble with radius R = 0.21 mm where immediate rupture occurred on first contact with the flat ethanol/air interface. The film thickness at rupture was estimated to be about 10 μm and is around 100 times larger than many other estimates of the film thickness at coalescence triggered by attractive surface forces.¹⁶ We can speculate that impurities or small dust particles at the top pool surface might be responsible for causing film rupture at such large thicknesses. However, we do not have any direct physical evidence for this so it is perhaps one interesting aspect that more detailed future experimental work can uncover.

It is also interesting to note from Fig. 6a that the effect of the free surface on the trajectory of the rising bubble only comes into effect when the two surfaces are less than a radius apart. This is a result that holds at a high Reynolds number where the influence of the free surface is not important until the bubble is very close (for the R = 0.21 mm case, Re ∼ 20). In contrast, under Stokes flow at Re ≪ 1, the hydrodynamic interaction between the bubble and the surface would start to affect the trajectory of the bubble when they are over 10 radii apart.³² In Fig. 6b we compare the velocity of the centre of mass for two selected cases from Fig. 6a.

On comparing the point of film rupture across water, silicone oil and ethanol in Fig. 4–6 it is observed that film rupture occurs consistently just after the moment the velocity of the centre of mass of the bubble reaches a local maximum whereas the number of bounces before rupture depends on the bubble size and approach velocity. It is not clear at this stage, if this is always the case. Perhaps experimental data could further confirm if this is generally true or a mere coincidence. The authors have no explanation why this should be the case from a theoretical point of view.

After showing that the model provides very good agreement for bubble collision and bounce results in different liquids and at different bubble sizes, we now use the model to predict features that were not yet measured experimentally. In Fig. 7 we show the film thickness h(r,t) at the maximum pool surface deformation during consecutive impacts using the same parameters as in Fig. 4a for a bubble with R = 0.74 mm rising in water. When the bubble impacts the surface a pressure builds up and inverts the curvature of the film thickness in the interaction region so that a dimple forms. As a consequence, the minimum thickness h_m(t) is no longer at the centre, r = 0, but at some larger radial position as shown in Fig. 7. The times corresponding to letters A, B, C and D will also be indicated in Fig. 8a.

Fig. 7 Film thickness at the maximum deformation of each bounce during multiple impacts of a bubble in water with R = 0.74 mm. For subsequent bounces the deformation becomes less pronounced. The central h₀ and minimum h_m film thicknesses are also defined.

Fig. 8 Time evolution of key physical parameters for a bubble with R = 0.74 mm in water colliding with a free surface corresponding to the case shown in Fig. 4a. (a) Time variation of the central deformation, z₀ of the free surface. The full numerical result (solid red) is almost identical to that given by the expression in eqn (34) using the numerical force F_F. Letters A, B, C and D correspond to the times in Fig. 7. (b) Time variation of different forces acting on the bubble. (c) Film thickness at centre h₀(t) and minimum film thickness h_m(t) at the rim of the dimple.

Surface deformation, forces and minimum film thickness h_m(t) and central film thickness h₀(t) as defined in Fig. 7 as a function of time are shown in Fig. 8. The maximum deformation of the free surface at r = 0 is shown in Fig. 8a. During the first impact, the free surface deforms by over 0.6 mm, a magnitude comparable to the radius of the bubble (R = 0.74 mm). The numerical solution is compared with the analytical solution by assuming that the force F in the equilibrium shape in eqn (16) is equal to the film force F_F calculated using eqn (26). However, note that the result in eqn (16) was derived for a system in equilibrium in which the force F was always positive. However, in a dynamic situation the film force F_F can become negative during retraction. Therefore, we need to be cautious in deriving an estimate of the radial extent of the interaction region, r₀ by

When the film force, F_F is negative the analytical shape is written as The maximum deformation z₀ can be found by matching both solutions in eqn (33) at r₀. The expression for z₀ is identical to the one previously derived in eqn (18) but is now a function of timeand is valid for positive or negative film force F_F that is calculated numerically using eqn (26). Here we have used eqn (8) for the shape of the free surface in the inner region.

In Fig. 8b we show the evolution of forces during consecutive bounces of the bubble. We can see that most of the bouncing behaviour is dominated by the balance between added mass and film force whereas buoyancy and drag play a greater role during the bubble rise stage when they balance each other. Just prior to bubble bounce we observe that the film force changes sign, meaning a ‘suction’ effect corresponds to the free surface deformation becoming negative. The amplitude of the interaction force becomes smaller at each impact but the duration of each impact remains mostly the same (∼10 ms).

In Fig. 8c, we show the film thickness at the centre h₀(t) and the minimum film thickness, h_m(t) at the rim of the dimple region. The results show very similar qualitative behaviour for each bounce that occurs at successively smaller average film thickness. Dimple formation is observed at each impact (i.e. when h_m ≠ h₀) and the minimum thickness becomes thicker before it becomes much thinner for a very short time just before the bubble fully separates from the free surface. Experimentally, coalescence was observed at 80 ms in this case when the minimum film thickness estimated from our model is around 10 μm. This value is quite large indicating that rupture was possibly not caused by surface forces that have a typical range of less than 1 μm.

In Fig. 9a we show in detail the close agreement between the numerical solution (solid red) and the result predicted by eqn (34) (dashed green) for the central deformation z₀ of the free surface. Unfortunately, no experimental data are available for comparison. The interaction region r₀ from eqn (32) is also shown for comparison (dotted blue). The kink at time around 13.5 ms represents the moment where the force becomes negative. We see from this plot that the maximum of the interaction region has almost the same size as the radius of the bubble. This justifies using a domain size that is larger than the radius of the bubble at r_m = 1.2R in our model.

Fig. 9 (a) A close up view of the time variation of the central deformation, z₀ (solid and dashed lines) of the free surface and the magnitude of the horizontal interaction region r₀ (dotted line) around the first impact of the bubble shown in Fig. 8a. (b) Shape of the deformable free surface during approach (left – dashed red) and rebound (right – solid blue) at different times during the first impact of the bubble. The dotted line indicates the free surface initial position.

In Fig. 9b we show the shape of the deformed free surface at selected times. During approach (red curves on the left) the surface deforms upward until a maximum deformation of about 0.6 mm is reached at around t = 10 ms. After that the surface starts its downward movement and overshoots becoming negative at around 13.5 ms and reaching the maximum free surface depression (or depth) at 14.2 ms before the surface returns to its original position and the bubble completely detaches from the surface. Later it will approach a second time and similar behaviour is observed but with lower amplitudes due to a lower impact velocity.

Conclusions

The model proposed for the rise and bounce of bubbles against soft deformable interfaces performs surprisingly well without any fitting parameters through a combination of lubrication, interfacial deformation and a force balance for the centre of mass of the bubble. The main novelty of the current approach is the incorporation of the deformation of the horizontal surface through a boundary condition that takes into account the behaviour of the surface at large radial distances, which is given analytically by a modified Bessel function of the second kind A range of bubble radii and approach velocities, as well as release distances were investigated and all gave excellent agreement with the experimental data. The model performed equally well for different liquid systems as long as the liquid is clean so that the appropriate boundary condition at the bubble surface is that the tangential stress vanishes. However, the time scale of the bubble surface collision is consistent with the assumption that the liquid at the boundary of the free surface exhibits tangentially immobile behaviour. It appears that the free surface, which is exposed to the laboratory atmosphere has accumulated sufficient trace impurities whereby the zero tangential stress condition no longer applies.³³

Appendix

Bessel function

The free surface obeys eqn (5) and in the main text it was shown that the solution of eqn (5) is a Bessel function of the second kind of order zero, K₀(r/λ). In the limit of very small and very large r/λ, K₀ can be approximated by a logarithmic and an exponential function respectively. When r ≪ λwhere the Euler constant γ_E = 0.57721566. On the other hand, when r ≫ λ The function K₀ and its asymptotic forms in eqn (A1) and (A2) are shown in Fig. 10. This justifies using the logarithmic approximation in the inner region for the theory.

Fig. 10 Bessel function K₀ of eqn (6) and its asymptotic behaviour for small and large radial distance according to eqn (A1) and (A2).

Vertical-to-horizontal aspect ratio vs. Weber number

To derive a relation between the bubble aspect ratio and the Weber number, We, we assume that an ellipsoidal bubble is rising with velocity V_T in a liquid as shown in Fig. 11.

Fig. 11 Schematic of an ellipsoidal bubble with equivalent radius R rising under buoyancy with velocity V_T in a bulk liquid.

Here we follow a derivation similar to Moore.²⁴ We start with the equation for an elliptic representation of the bubble

where the curvatures at the top and side of the bubble are given by The pressure P along a sphere with radius R according to potential flow iswhere P₀ is the ambient pressure and θ is the angle starting from the top. The pressure jump across the bubble surface is given by surface tension times the curvature. Using the definitions in Fig. 11, and assuming that the pressure inside the bubble P_in remains constant, we obtain for the top and side of the bubble:

P_in − P_top = σk_top (A6)

P_in − P_side = σk_side (A7)

The pressure at the top is obtained using eqn (A5)and on the side where θ = 90° so sin θ = 1.Using eqn (A6)–(A9)and eqn (A4):where we define the horizontal-to-vertical aspect ratio χ = R_h/R_v. Using eqn (A10) and R³ = R_h²R_v = χ²R_v³, allows us to write, R = χ^2/3R_v Using the definition of the Weber number (We = 2RρV_T²/σ) we can rewrite eqn (A12) as It is observed experimentally that the relation between the Weber number and the inverse of the aspect ratio appears to be linear in the range of current interest (see Fig. 3b). Using this physical insight and performing a Taylor expansion with respect to 1/χ neglecting higher order terms results in the following relation:Eqn (A14) has also been proposed by Legendre et al.²⁵ It gives very good agreement with experimental observations of bubbles rising in all liquids considered here: water,⁴ silicone oil⁵ and ethanol.⁶

Acknowledgements

This work was supported in part by the Australian Research Council through a Discovery Project Grant to DYCC.

References

1. K. Malysa, M. Krasowska and M. Krzan, Adv. Colloid Interface Sci., 2005, 114–115, 205–225.
2. D. Kosior, J. Zawala, M. Krasowska and K. Malysa, Colloids Surf., A, 2014, 441, 788.
3. J. Zawala and K. Malysa, Langmuir, 2011, 27, 2250.
4. D. Kosior, J. Zawala, R. Todorov, D. Exerowa and K. Malysa, Colloids Surf., A, 2014, 460, 391.
5. J. Zawala, S. Dorbolo, D. Terwagne, N. Vandewalle and K. Malysa, Soft Matter, 2011, 7, 6719.
6. F. Suñol and R. González-Cinca, Colloids Surf., A, 2010, 365, 36.
7. J. Feng, M. Roché, D. Vigolo, L. N. Arnaudov, S. D. Stoyanov, T. D. Gurkov, G. G. Tsutsumanova and H. A. Stone, Nat. Phys., 2014, 10, 606.
8. P. S. Stewart, J. Feng, L. S. Kimpton, I. M. Griffiths and H. A. Stone, J. Fluid Mech., 2015, 777, 27.
9. J. C. Bird, R. de Ruiter, L. Courbin and H. A. Stone, Nature, 2010, 465, 759.
10. H.-K. Tsao and D. L. Koch, Phys. Fluids, 1997, 9, 44.
11. M. Fujasová-Zedníková, L. Vobecká and J. Vejrazka, Can. J. Chem. Eng., 2010, 88, 473.
12. R. Manica, E. Klaseboer and D. Y. C. Chan, Langmuir, 2015, 31, 6763.
13. E. Klaseboer, J.-Ph. Chevaillier, A. Maté, O. Masbernat and C. Gourdon, Phys. Fluids, 2001, 13, 45.
14. E. Klaseboer, R. Manica, M. H. W. Hendrix, C.-D. Ohl and D. Y. C. Chan, Phys. Fluids, 2014, 26, 092101.
15. A. Albadawi, D. B. Donoghue, A. J. Robinson, D. B. Murray and Y. M. C. Delauré, Int. J. Multiphase Flow, 2014, 65, 82.
16. D. Y. C. Chan, E. Klaseboer and R. Manica, Soft Matter, 2011, 7, 2235.
17. D. Y. C. Chan, E. Klaseboer and R. Manica, Adv. Colloid Interface Sci., 2011, 165, 70.
18. R. Manica, M. H. W. Hendrix, R. Gupta, E. Klaseboer, C.-D. Ohl and D. Y. C. Chan, Soft Matter, 2013, 9, 9755.
19. R. Manica, M. H. W. Hendrix, R. Gupta, E. Klaseboer, C.-D. Ohl and D. Y. C. Chan, Appl. Math. Model., 2014, 38, 4249.
20. M. H. W. Hendrix, R. Manica, E. Klaseboer, D. Y. C. Chan and C.-D. Ohl, Phys. Rev. Lett., 2012, 108, 247803.
21. P. C. Duineveld, J. Fluid Mech., 1995, 292, 325.
22. F. Raymond and J.-M. Rosant, Chem. Eng. Sci., 2000, 55, 943.
23. M. Wu and M. Gharib, Phys. Fluids, 2002, 14, 49.
24. D. W. Moore, J. Fluid Mech., 1965, 23, 749.
25. D. Legendre, R. Zenit and J. R. Velez-Cordero, Phys. Fluids, 2012, 24, 043303.
26. J. Magnaudet and I. Eames, Annu. Rev. Fluid Mech., 2000, 32, 659.
27. G. Mougin and J. Magnaudet, Phys. Rev. Lett., 2002, 88, 014502.
28. R. Zenit and J. Magnaudet, Phys. Fluids, 2008, 20, 061702.
29. A. K. Chesters and G. Hofman, Appl. Sci. Res., 1982, 38, 353.
30. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, 1972.
31. T. Miloh, J. Eng. Math., 1977, 11, 349.
32. R. Manica, L. Parkinson, J. Ralston and D. Y. C. Chan, J. Phys. Chem. C, 2010, 114, 1942.
33. G. H. Kelsall, S. Tang, A. L. Smith and S. Yurdakul, Faraday Trans., 1996, 92, 3879.

---

Footnote

† Electronic supplementary information (ESI) available: Video highlighting numerical calculations of a bubble (R = 0.74 mm) impacting a fluid interface corresponding to the data of Fig. 4a. See DOI: 10.1039/c5sm03151f

---