Path instability of deformable bubbles rising in Newtonian liquids: a linear study

Abstract

The first stages of the path instability phenomenon affecting the buoyancy-driven motion of gas bubbles rising in weakly or moderately viscous liquids are examined using a recently developed numerical approach designed to assess the global linear stability of incompressible flows involving freely evolving interfaces. Predictions for the critical bubble size and frequency of the most unstable mode are found to agree well with reference data obtained in ultrapure water and in several silicone oils. By varying the bubble size, stability diagrams are built for several specific fluids, revealing three distinct regimes with different bifurcation sequences. The spatial structure of the unstable modes is analysed, together with the variations of the bubble shape, position and orientation. For this purpose, displacements of the bubble surface are split into rigid-body components and volume-preserving deformations, allowing us to determine how the relative magnitude of the latter varies with the fluid properties and bubble size. Predictions obtained with freely deformable bubbles are compared with those found by maintaining the bubble shape determined in the base state frozen during the stability analysis. This comparison reveals that deformations leave the phenomenology of the first bifurcations unchanged in low-viscosity fluids, especially water. Hence, in such fluids, bubbles behave essentially as freely moving rigid bodies submitted to constant-force and zero-torque constraints, at the surface of which the fluid obeys a shear-free condition. In contrast, deformations change the nature of the primary bifurcation in oils slightly more viscous than water, whereas, somewhat surprisingly, they leave the near-threshold phenomenology unchanged in more viscous oils.

1. Introduction

A great deal of effort has been invested since the middle of the last century to elucidating the origin and predicting the characteristics of the intriguing path instability phenomenon affecting the buoyancy-driven motion of millimetre-sized gas bubbles rising in weakly or moderately viscous liquids. Attempts were mostly experimental during the first decades, before numerical simulation techniques became mature enough to shed some new light on this puzzling phenomenon. Among them, linear stability analyses have kept pace with advances in the computational techniques required to determine efficiently the discretized form of the stationary solutions of the Navier–Stokes equations subjected to appropriate boundary conditions, and to solve subsequently the large-size generalized eigenvalues problems giving access to the possibly linearly unstable eigenmodes.

Early studies considered bubbles with a frozen shape held fixed in a uniform stream. Assuming an oblate spheroidal shape with an arbitrary aspect ratio (the ratio of the major and minor axes lengths, hereinafter denoted as $\chi$), with the short axis of the spheroid aligned with the incoming flow, conditions under which the axisymmetric wake becomes unstable, the nature of the corresponding first bifurcations and the structure of the associated unstable modes could be explored (Yang & Prosperetti 2007; Tchoufag, Magnaudet & Fabre 2013). It was eventually concluded that wake instability occurs within a finite range of Reynolds number when $\chi \gtrsim 2.21$, in agreement with predictions provided by the numerical solution of the full Navier–Stokes equations in the same configuration (Magnaudet & Mougin 2007). The first bifurcation is stationary, i.e. the corresponding eigenvalue is real. Beyond this bifurcation, the stationary wake exhibits a pair of counter-rotating trailing vortices which are responsible for a non-zero lift force. For larger aspect ratios, increasing the Reynolds number beyond the primary threshold while keeping $\chi$ fixed or vice versa, a second bifurcation of Hopf type takes place. The wake keeps the previous planar symmetry unchanged beyond the corresponding threshold but vortices are periodically shed downstream, resulting in periodic oscillations of the lift force about a non-zero mean. This transition scenario, with the axisymmetric wake becoming unstable through a stationary bifurcation and the resulting non-axisymmetric stationary wake undergoing a Hopf bifurcation, is similar to that obeyed by the wake of rigid spheres and discs (Natarajan & Acrivos 1993).

Compared with the actual problem of buoyancy-driven rising bubbles, the configuration contemplated in the above studies is purely academic, as the bubble is held fixed despite the existence of a non-zero lift force. An important step towards physically realistic conditions was achieved during the last decade, starting with global linear stability studies of the flow past freely moving buoyancy/gravity-driven rigid bodies with various shapes. In this situation, the body motion is subjected to two additional constraints, since the hydrodynamic force must balance the net body weight at all times (implying that the horizontal force remains null whatever the body orientation and velocity), and the hydrodynamic torque must also remain zero, provided that the body density is uniform. Therefore, disturbances in the body velocity or orientation modify the velocity and pressure distributions in the surrounding fluid, which in turn induce changes in the force and torque acting on the body that must respect the above constraints. Two-dimensional plates and rods were considered by Fabre, Assemat & Magnaudet (2011) and Assemat, Fabre & Magnaudet (2012). The former study focused on ‘heavy’ bodies with a density much larger than that of the fluid. An asymptotic approach could then be employed and revealed the existence of four ‘aerodynamic’ modes distinct from the ‘fluid’ (or ‘wake’) modes associated with wake instability (qualitatively similar modes are encountered in flight dynamics, and those evidenced in the above studies were coined after them). Two of these aerodynamic modes have real negative eigenvalues, i.e. they are always stable. In contrast, the other two are associated with a pair of complex conjugate eigenvalues, hence with a Hopf bifurcation. The corresponding oscillations are slow compared with those typical of vortex shedding in two-dimensional wakes, and their frequency varies as the inverse square root of the body-to-fluid density ratio, $m^*$, thus becoming vanishingly small for very heavy bodies. Existence of these modes was confirmed for arbitrary $m^*$ in the global linear stability analysis (hereinafter abbreviated as GLSA) performed numerically by Assemat et al. (2012), and the corresponding thresholds were compared with those of the fluid modes for thin plates and square rods. Increasing the Reynolds number, it was observed that, for infinitely thin plates, the primary wake mode always becomes unstable first, no matter what the value of $m^*$ is. The same conclusion was found to hold for square rods as long as the density ratio exceeds $1.22$. In contrast, for lower $m^*$, the aerodynamic oscillatory mode is destabilized first and exhibits a frequency (corresponding to a fluttering of the falling/rising rod about the vertical) smaller than that of the wake mode by approximately a factor of two. Similar conclusions were previously obtained through fully resolved numerical simulations by Alben (2008) with light two-dimensional ellipses, and the critical Reynolds number at which the fluttering motion sets in was found to be approximately half that corresponding to the onset of wake instability. The two-dimensional GLSA approach of Assemat et al. (2012) was extended to axisymmetric bodies by Tchoufag, Fabre & Magnaudet (2014a), who considered the stability of rising/falling discs and short cylinders. Whatever the body aspect ratio, they found the aerodynamic low-frequency mode to be always the most unstable beyond a critical $\chi$-dependent $m^*$ (note the difference with two-dimensional square rods for which the same happens below a critical $m^*$). For lower $m^*$, the most unstable mode is either the one associated with the wake (here characterized by weak deviations of the path in the presence of sustained fluid oscillations at the back of the body), or a stationary mode in which the body follows an inclined path and its wake is made of a pair of steady counter-rotating streamwise vortices. Based on these findings, Tchoufag et al. (2014a) concluded that the first deviation of the path of such freely moving axisymmetric bodies cannot in general be predicted by considering the instability of the sole wake. Rather, the instability of their path results under most conditions from the intrinsic coupling between the body and fluid implied by the constant-force and zero-torque conditions.

Tchoufag, Magnaudet & Fabre (2014b) applied the GLSA approach to freely rising spheroidal bubbles with a prescribed oblateness. They observed that the vertical path becomes unstable within a finite range of Reynolds number when the bubble oblateness exceeds $2.15$, which is slightly lower than the wake instability threshold determined for a fixed bubble (Magnaudet & Mougin 2007; Tchoufag et al. 2013). The unstable Reynolds number range widens as $\chi$ increases. For $2.15<\chi <2.25$, increasing the Reynolds number, i.e. the bubble size, while maintaining the oblateness fixed, the path first becomes unstable within a narrow Reynolds number range through a low-frequency mode similar to the aerodynamic mode encountered with heavy short cylinders and discs. Increasing the Reynolds number further, a stationary mode yielding an inclined path and growing much faster than the previous low-frequency mode becomes dominant. This situation prevails within an intermediate Reynolds number range, beyond which the low-frequency mode becomes dominant again. Last, the system restabilizes beyond a fourth critical Reynolds number. For instance, the path of a bubble with $\chi =2.22$ becomes oscillatory through a low-frequency mode in the two separate ranges $310\lesssim Re\lesssim 340$ and $450\lesssim Re\lesssim 980$, while it takes a non-zero stationary inclination with respect to the vertical in the intermediate range $340\lesssim Re\lesssim 450$, the Reynolds number $Re$ being based on the bubble rise speed and equivalent diameter. Obviously, the description provided by this simplified model suffers from the fact that the bubble shape is not allowed to vary with the Reynolds number, in contrast to what happens under real conditions. Nevertheless, as will become apparent later, the first half of the above scenario, i.e. the succession of a low-frequency mode becoming unstable first and a stationary mode taking over it, is relevant to the case of deformable bubbles rising in low-viscosity high-surface-tension fluids, especially water.

A somewhat more sophisticated track was followed by Cano-Lozano et al. (2016a). Rather than prescribing arbitrarily an oblate spheroidal shape, they computed the actual bubble shape and rise speed corresponding to the stationary axisymmetric base state, using the open software Gerris (Popinet 2003, 2007). Then, they introduced this shape and speed into the GLSA solver of Tchoufag et al. (2014b), keeping the shape frozen despite the possible development of flow and path instabilities. With this approach, they could determine for a wide variety of Newtonian fluids (characterized by their density, viscosity and surface tension), the critical Reynolds number, i.e. the critical bubble size, beyond which the vertical path of a bubble with a nearly realistic shape becomes unstable. They could also compare these predictions with those corresponding to ‘pure’ wake instability, obtained by holding the same bubble fixed in a uniform stream. For low-viscosity high-surface-tension fluids, they found that, in agreement with the conclusions of Tchoufag et al. (2014b), the system first becomes unstable through a Hopf bifurcation leading to low-frequency path oscillations. At the corresponding Reynolds number, the wake of the same bubble held fixed is still stable, which proves that the mechanism responsible for path instability is intimately linked to the constant-force and torque-free conditions that constrain the dynamics of the coupled fluid–body system. For liquids with a viscosity $5$ to $10$ times larger than water and a surface tension 3–4 times smaller, they observed that path instability first arises through an oscillatory mode with a frequency 4–5 times higher than that found in low-viscosity high-surface-tension liquids. However, the picture looked dramatically different in the intermediate range corresponding to ‘moderate’ liquid viscosities and surface tensions. In particular, the instability thresholds of freely moving and fixed bubbles were found to be identical in that range. Hence, the first unstable mode is stationary, similar to that observed for fixed spheroidal bubbles (Magnaudet & Mougin 2007; Yang & Prosperetti 2007; Tchoufag et al. 2013). When compared with reference experiments carried out in ultrapure water (Duineveld 1995; De Vries 2002) or silicone oils (Zenit & Magnaudet 2008; Sato 2009), these predictions experience mixed successes. In low-viscosity high-surface-tension fluids, they properly capture the low-frequency oscillatory nature of the incipient path instability, and slightly over-predict the size of the smallest bubble whose path becomes unstable. The same conclusion holds for the ‘high’-frequency path oscillations observed in oils 5–10 times more viscous than water. In contrast, for oils with an intermediate viscosity, experimental observations provide evidence that the first instability of the system keeps a similar oscillatory nature, at odds with the above predictions.

This blatant disagreement points to the role of transient bubble deformations not accounted for in the above studies. It provided one of the main motivations that decided us to undertake the development of a more general linear stability approach capable of dealing with free gas–liquid boundaries, thanks to which the fate of time-dependent deformations of the bubble–fluid interface may be predicted together with that of flow and path disturbances. The corresponding method was successfully developed and validated by Bonnefis (2019); other groups have pursued similar objectives in parallel (Zhou & Dusek 2017; Herrada & Eggers 2023). These efforts may be considered as the third and final step on the long route to the linear stability analysis of the flow past rising bubbles. To summarize, in the first step, the bubble shape was frozen, its centroid was assumed to remain fixed and only the stability of the wake was assessed (Yang & Prosperetti 2007; Tchoufag et al. 2013). In the next step, the bubble shape was still frozen but the bubble centroid was allowed to move freely under the effect of buoyancy, being only constrained by the constant-force and zero-torque conditions (Tchoufag et al. 2014b; Cano-Lozano et al. 2016a). Most of the investigations performed in these first two steps assumed the bubble to keep an exact oblate spheroidal shape, except the last of them, in which this shape was obtained as part of the stationary solution of the full Navier–Stokes equations. The time-dependent adjustment of the bubble shape to its possibly non-straight unsteady motion was the last ingredient missing to reach a fully realistic description of the bubble–flow interaction. This is what the third step, achieved by Bonnefis (2019), enables. The purpose of this paper is to present and discuss the predictions provided by this approach across a wide range of fluids, from low-viscosity high-surface-tension liquid metals to oils with a viscosity typically ten times larger and a surface tension four times smaller than those of water, and to compare them with reference data when available. Predictions obtained in the specific case of water were already summarized by Bonnefis, Fabre & Magnaudet (2023) and are discussed in more detail here.

The numerical approach, already described by Sierra-Ausin et al. (2022), is summarized in § 2. The characteristics of the base state are presented in § 3, before the neutral curves resulting from the linear stability analysis are discussed in § 4, and the nature and sequencing of the first unstable modes are examined in § 5. Then, the relative magnitude of deformations of the bubble–fluid interface with respect to the lateral drift of the bubble centroid are discussed in § 6. Section 7 summarizes the main findings of the study and discusses the respective roles of wake, bubble deformation and fluid–bubble hydrodynamic couplings in the different regimes encountered by varying the physical properties of the carrying liquid.

2. Problem formulation and numerical approach

We assume that a single gas bubble with volume $\mathcal {V}$ and equivalent diameter $D=(({6}/{{\rm \pi} })\mathcal {V})^{1/3}$ rises through a large expanse of a Newtonian liquid with density $\rho$, viscosity $\mu$ and surface tension $\gamma$ under the effect of gravity $g$. The liquid is at rest at infinity and the disturbance flow induced by the bubble ascent is assumed incompressible. The bubble is initially axisymmetric and is released with its symmetry axis aligned with gravity. Considering that the density and viscosity of the gas enclosed in the bubble are negligibly small compared with those of the surrounding liquid, the problem depends on two dimensionless parameters, namely the Bond number $Bo=\rho g D^2/\gamma$ and the Galilei number $Ga=\rho (gD^3)^{1/2}/\mu$. These characteristic numbers compare the gravitational force $\rho g D^3$ driving the bubble ascent with the capillary force $\gamma D$ and the viscous force $\mu (gD^3)^{1/2}$, respectively. Whatever the bubble size, a given liquid placed in a given gravitational environment is entirely characterized by the value of the so-called Morton number $Mo=Bo^3Ga^{-4}=g\mu ^4/(\rho \gamma ^3)$. By defining the visco-gravitational diffusion length scale $l_\mu =[\mu ^2/(\rho ^2g)]^{1/3}$ (the length over which viscosity diffuses a disturbance within a $[\mu /(\rho g^2)]^{1/3}$-long period of time) and the capillary length scale $l_\gamma =[\gamma /(\rho g)]^{1/2}$, it is seen that $l_\mu /l_\gamma =Mo^{1/6}$. Hence, the Morton number characterizes the ratio of these two length scales. The bubble diameter, $D$, the gravitational time, $(D/g)^{1/2}$, and their ratio, $(gD)^{1/2}$, will be used throughout the paper to make lengths, frequencies and velocities dimensionless. Once the bubble rise speed $u_b$ is known, the non-dimensional rise speed, $U=u_b/(gD)^{1/2}$ allows the bubble Reynolds number, $Re=Ga\,U$, and Weber number, $We=Bo\,U^2$, to be evaluated.

The stability of the coupled fluid–bubble motion is examined with the help of the linearized arbitrary Lagrangian–Eulerian approach, hereinafter abbreviated as L-ALE, developed by Bonnefis (2019) and described by Sierra-Ausin et al. (2022). Here, we only summarize the main steps of this approach, referring readers to the above references for more details. The key idea is to solve the governing equations of the problem, which hold on a time-deforming domain $\varOmega (t)$ since the bubble–fluid interface $\varGamma _b(t)$ experiences time-dependent deformations, on a prescribed reference domain $\varOmega _0$ bounded internally by a fixed surface $\varGamma _0$ representing the reference position of the interface (figure 1). In addition to $\varGamma _0$, $\varOmega _0$ is bounded by a fixed external surface $\varGamma _\infty$ on which suitable far-field conditions are imposed, and $\varGamma _0$ and $\varGamma _\infty$ are connected by the vertical axis, $\varGamma _z$, about which the base state exhibits an axial symmetry. The L-ALE approach is grounded on the explicit building of the bijective mapping connecting the instantaneous position $\boldsymbol {x}$ of a given geometrical point in the deforming domain to its position $\boldsymbol {x}_0$ in $\varOmega _0$. In other words, this approach relies on the determination of the infinitesimal displacement field $\boldsymbol {\xi }(\boldsymbol {x}_0)$ such that $\boldsymbol {x}=\boldsymbol {x}_0+\boldsymbol {\xi }(\boldsymbol {x}_0)$ everywhere in $\varOmega _0\cup \varGamma _0\cup \varGamma _z\cup \varGamma _\infty$. On $\varGamma _0$, the unit normal of which is $\boldsymbol {n}_0$, the normal component of the displacement, $\boldsymbol {\xi }\boldsymbol {\cdot }\boldsymbol {n}_0$, must equal the displacement of the bubble–fluid interface, $\eta$. In contrast, $\boldsymbol {\xi }$ may be arbitrarily chosen anywhere else. To ensure a smooth distribution of the displacement field, we assume that $\boldsymbol {\xi }$ obeys the Cauchy–Navier equation of elastostatics with unit Lamé coefficients in $\varOmega _0$, together with a suitable symmetry condition on $\varGamma _z$ and a Dirichlet condition $\boldsymbol {\xi }={\boldsymbol {0}}$ on $\varGamma _\infty$.

Figure 1. Flow domain and geometrical transformation involved in the L-ALE approach. The black line represents the current bubble–fluid interface $\varGamma _b(t)$ bounding internally the physical fluid domain $\varOmega (t)$, while the blue line corresponds to the initial interface $\varGamma _{0}$ bounding the reference domain $\varOmega _0$. Both domains exhibit a rotational symmetry about the $\varGamma _z$ axis.

Knowing $\boldsymbol {\xi }(\boldsymbol {x}_0)$ and its gradients everywhere, the time and space derivatives involved in the governing equations written in $\varOmega (t)$ may be expressed in $\varOmega _0$ and linearized consistently with respect to $\boldsymbol {\xi }$, together with the variations of the unit normal and local mean curvature of $\varGamma _b(t)$. Governing equations in $\varOmega _0$ were detailed in Sierra-Ausin et al. (2022). They consist of the continuity and Navier–Stokes equations, while the kinematic (no-penetration) and dynamic boundary conditions must be enforced on $\varGamma _0$. Owing to the negligible viscosity of the gas enclosed in the bubble, the tangential component of the dynamic boundary condition reduces to a shear-free condition that must be satisfied on $\varGamma _0$ by the gradients of the liquid velocity field, $\boldsymbol {u}$. Since the gas density is also negligible, the momentum balance implies that the pressure is uniform (but possibly time dependent) within the bubble. Consequently, the normal component of the dynamic boundary condition on $\varGamma _0$ expresses the fact that the difference between the uniform pressure $p_b$ inside the bubble and the local pressure $p$ in the liquid (from which the hydrostatic component has been subtracted) balances the difference between the local capillary pressure and the normal viscous stress in the liquid. Last, two integral conditions have to be added to ensure that the volume of the bubble does not change over time, and that the vertical position of its geometrical centroid does not vary once the problem has been expressed in a reference frame translating with the bubble rise speed, $u_b$.

A key feature of the L-ALE approach is that the entire set of equations governing the evolution of the state vector ${\boldsymbol {q}}=[{\boldsymbol {u}},p,\eta,\boldsymbol {\xi },u_b, p_b]^\text {T}$ ($^\text {T}$ denoting the transpose) is solved simultaneously (such an approach is frequently referred to as ‘monolithic’). Starting from an arbitrary initial solution ${\boldsymbol {q}}_0=[{\boldsymbol {u}}_0,p_0,0,{\boldsymbol {0}},u_{b0},p_{b0}]^\text {T}$ (here assumed to be axisymmetric) and dropping time derivatives in the Navier–Stokes equations and in the no-penetration condition, a Newton algorithm is used to obtain the stationary base state through a series of global iterations in which the state vector and the reference domain are iteratively updated. For this purpose, we express the problem in weak form and make use of the finite element software FreeFem$++$ (Hecht 2012) to build and invert the various matrices involved. The volume fields $({\boldsymbol {u}},p,{\boldsymbol \xi })$ are discretized on the finite element basis suitable for each of them, while the normal displacement of the interface, $\eta$, is discretized on the local Fourier basis. Provided that ${\boldsymbol {q}}_0$ is chosen ‘not too far’ from the stationary solution (i.e. within its basin of attraction), the above algorithm converges quadratically in a few iterations. Moreover, compared with time-marching algorithms, it has the decisive advantage that unstable steady solutions may also be captured, allowing the features of the corresponding branches (if any) of the bifurcation diagram to be studied in detail; see, e.g. Sierra-Ausin et al. (2022).

Once the state vector and the shape of the bubble–fluid interface defining the axisymmetric base state have been obtained, the linear stability of this solution is assessed in the Galilean reference frame translating with the corresponding rise speed. For this purpose, we consider a cylindrical coordinate system $(r,\theta,z)$ whose axis $r=0$ and cross-sectional plane $z=0$ correspond to the symmetry axis and vertical position of the bubble geometrical centroid in the base state, respectively. Then we compute the complex eigenvalues $\lambda =\lambda _r+\text {i}\lambda _i$ associated with disturbances of the form ${\boldsymbol {q}}'(r,\theta,z,t)=[({\hat {\boldsymbol {u}}},\hat {p},\hat {\eta },\hat {\boldsymbol {\xi }}) {\rm e}^{\text {i}m\theta },\hat {p}_b]^\text {T}{\rm e}^{\lambda t}$, with $m$ denoting the azimuthal wavenumber and $\theta$ the azimuthal angle, the hatted amplitudes $({\hat {\boldsymbol {u}}},\hat {p},\hat {\eta },\hat {\boldsymbol {\xi }})$ depending on both $r$ and $z$ (owing to the choice of the reference frame, $\hat {u}_b\equiv 0$, so that the instantaneous position of the bubble centroid no longer necessarily coincides with the origin of the reference frame). In what follows, we shall be primarily interested in the family of non-axisymmetric modes $|m|=1$ which are those associated with the first deviations of the bubble path from the vertical, and in modes $m=0$ and $|m|=2$ associated with the oscillations of the bubble shape. The eigenvector ${\boldsymbol {q}}'$ is obtained via the Newton algorithm already employed to determine the base state. The set of equations to be solved is similar, except that the momentum equation and the no-penetration condition now involve terms proportional to $\lambda$ since the sought solution is time dependent; the reference state vector and domain to be considered are those corresponding to the base state. The complex eigenvalues are finally obtained using a Krylov–Shur projection method available in the SLEPc library.

Obviously, the accuracy with which the base state and the eigenvalues are determined depends critically on the finite element triangulation of $\varOmega _0$, especially along the bubble–fluid interface $\varGamma _0$ and within the boundary layer and wake regions. The position of the fictitious outer boundary $\varGamma _\infty$ with respect to the bubble is also important, the computational domain having to be large enough to avoid spurious confinement effects. It must also be kept in mind that, in low-Morton-number fluids such as water ($Mo\approx 2.54\times 10^{-11}$ under standard conditions) and even more in liquid metals, such as Galinstan ($Mo\approx 1.4\times 10^{-13}$), the onset of path instability takes place in regimes where the bubble Reynolds number $Re=\rho u_bD/\mu$ is of the order of $10^3$ (${\approx }670$ in water, ${\approx }2100$ in Galinstan, see figure 7 below). To ensure that the flow is still fully resolved in the thin boundary layers encountered in such regimes, we design an initial grid in which $18$ nodes are distributed across a $5\times Re^{-1/2}$-thick layer surrounding the bubble. Along the bubble surface, the number of nodes $N_D$ over a $D$-long arc length is approximately $7.2\times Re^{1/2}$, with at least $N_D=64$ when $Re$ becomes less than $80$. Since $u_b$, hence $Re$, is not known a priori, once an approximate steady state has been computed on the initial grid, a new grid is generated following an adaptive mesh refinement procedure. This new grid is designed based on the previous steady-state solution, with some constraints ensuring that the boundary layer remains sufficiently resolved. This procedure is repeated a few times until the solution becomes grid independent. Extensive tests were carried out to evaluate the sensitivity of the base state (appreciated for instance through the variations of $Re$, $\chi$, $\lambda _r$ and $\lambda _i$) to factors such as the domain size, distance from the bubble to the downstream boundary and density of nodes along the bubble surface. For instance, in the case of a bubble with $Bo=0.6$ rising in water, it was found that relative variations of these indicators for the non-axisymmetric modes $m=\pm 1$ are all less than $0.1\,\%$ provided that $N_D\gtrsim 70$ and the distance from the bubble centroid to the downstream end of the domain is at least $30D$, the upstream end and the lateral surface being both located $15D$ away from this centroid (Bonnefis 2019).

3. Base state

3.1. Bubble shapes

Figure 2 shows how the equilibrium shape of the bubble changes with the Bond number in three different fluids. As the Morton number keeps a constant value in each fluid, increasing $Bo$ (or $Ga$) in a given series amounts to increasing the bubble size; since $Bo\propto D^2$, the volume of the bubble grows as $Bo^{3/2}$. In each series, starting from nearly spherical shapes at low Bond number, bubbles take oblate spheroidal shapes when finite-$Bo$ effects become significant. The figure makes it clear that the properties of the carrying fluid have a major influence on the range of sizes in which this geometric approximation is valid: while the shape of a bubble with $Bo=3$ is still close to a spheroid in silicone oil DMS-T05 ($Mo= 6.2\times 10^{-7}$, see table 1 for the physical properties of the various fluids discussed below), a bubble with $Bo=0.5$ is seen to already exhibit a significant fore–aft asymmetry in water. Beyond that range, the front part of the bubble flattens increasingly as its size increases, while its rear part becomes more and more rounded. As indicated by the change in colour of the bubble contour in the figure, the transition between these two ‘regimes’ corresponds approximately to the onset of path instability. Although they do not rise in a straight line, bubbles with a mildly convex front and a pronounced rounded rear are routinely observed in experiments. In contrast, bubbles with an almost flat or even concave front such as those shown in the figure for $Bo\gtrsim 5$ in DMS-T05 or $Bo\gtrsim 0.7$ in water are not. The obvious reason is that such bubbles actually experience large shape oscillations, so that the present stationary solution makes little sense in these regimes.

Figure 2. Equilibrium shapes of bubbles with various Bond numbers in three different fluids; the Bond number is specified within each contour. Contours represent different bubbles; their centroids are arbitrarily shifted in the vertical direction for readability, making the local value of $z$ irrelevant. Blue and red contours refer to stable and unstable configurations, respectively. (a) Silicone oil DMS-T05 ($Mo=6.2\times 10^{-7}$); (b) silicone oil DMS-T02 ($Mo=1.6\times 10^{-8}$); (c) water at $20\,^\circ \text {C}$ ($Mo=2.54\times 10^{-11}$).

Table 1. Physical properties of some specific fluids considered in this study. Note that the surface tension of Galinstan may vary from $535$ to $718\ {\rm mN}\ {\rm m}^{-1}$, depending on the exact alloy composition.

3.2. Rise speed

Variations of the equilibrium rise speed with the bubble size and liquid properties are displayed in figure 3 in the form of a $Re(Bo)$ diagram. Comparison with experimental data reveals an excellent agreement in most fluids, especially in the most challenging case of water. A consistent $10\,\%$ over-estimate may be noticed in the least viscous silicone oil, DMS-T00, presumably because of some variation in the actual viscosity or surface tension of the sample used in the experiments with respect to the reference values provided by the manufacturer (Zenit & Magnaudet 2008). In each series, the location of the incipient path instability detected in the experiments is specified with a bullet in the figure. The threshold Reynolds number is seen to vary by nearly one order of magnitude from water to the most viscous oil, DMS-T11, $9.4$ times more viscous than water. As the three intermediate series make clear, the predicted Reynolds number departs from that determined in the experiments beyond the threshold. This is no surprise, since part of the potential energy of zigzagging or spiralling bubbles is transferred to the radial and azimuthal fluid motions, yielding a decrease in the rise speed compared with the rectilinear regime.

Figure 3. Rise Reynolds number vs the Bond number in several fluids. Solid lines: present predictions; symbols: experimental data, with $\boxdot$ (blue): ultrapure water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$) (Duineveld 1995), ${\times }$ (orange): silicone oil DMS-T00 ($Mo=1.8\times 10^{-10}$), ${\odot }$ (yellow): DMS-T02 ($Mo=1.6\times 10^{-8}$), (green): DMS-T05 ($Mo=6.2\times 10^{-7}$), and (purple): DMS-T11 ($Mo=9.9\times 10^{-6}$) (data for all four series taken from Zenit & Magnaudet 2008). Red bullets mark the onset of path instability detected experimentally in each series.

Data reported in figure 3 may also be used to shed some light on the influence of the bubble shape and properties of the carrying fluid on the bubble rise speed. To this end, the above numerical predictions for the rise Reynolds number are replotted against the Galilei number in figure 4. This figure reveals the succession of two distinct behaviours in each fluid. First, at low enough $Ga$, all data collapse onto a master curve following the power law $Re\propto Ga^{5/3}$. This behaviour corresponds to the regime in which the bubble rise speed does not depend on surface tension, hence on the Morton number. In this regime, the bubble shape remains close to a sphere, i.e. its aspect ratio $\chi$ is close to unity. For low enough $We$, departures from sphericity are known to increase linearly with the Weber number according to the law $\chi =1+\frac {9}{64}We$ (Moore 1965). Since $We=Re^2Mo^{1/3}Ga^{-2/3}$, the lower $Mo$ the wider the $Re$ range pertaining to this first regime, which extends up to $Re\approx 200$ in water. Bubbles become flatter as $We$, hence $Re$, increases, opposing a larger resistance to the fluid. This makes their rise speed grow more slowly with $Ga$ when the Weber number becomes of order unity. Since their frontal area depends on $We$, which itself depends on $Mo$, the Reynolds number is no longer independent of the Morton number in this second regime. To further examine the two behaviours, one has to refer to the force balance on the bubble which, in non-dimensional form, reads $C_D(Re)Re^2=\frac {4}{3}Ga^2$, with $C_D(Re)$ the usual drag coefficient. The drag coefficient of a spherical bubble is known to vary from $16Re^{-1}$ at low Reynolds number to $48Re^{-1}$ in the limit of very large Reynolds number (Batchelor 1967). Hence, for nearly spherical bubbles, the force balance implies $k(Re)Re=\frac {1}{12}Ga^2$, with $k(Re\ll 1)=1$, $k(Re\gg 1)=3$. The slow increase of $k(Re)$ in between these two limits, which follows approximately the power law $k(Re)\sim Re^{1/5}$ up to Reynolds numbers of a few hundreds, is the reason why the Reynolds number grows slightly more slowly than $Ga^2$ in the first regime. Bubbles are significantly distorted in the second regime, which corresponds to the intermediate range $We={O}(1)$ ($We\approx 0.62$ at $Re=200$ in water). Assuming that their frontal area, which is proportional to $\chi ^{2/3}$, grows approximately linearly with the Weber number in that range yields $C_D(Re)\sim k(Re)Re^{-1}We\sim Re^{6/5}Mo^{1/3}Ga^{-2/3}$. The force balance then implies $Re\sim Ga^{5/6}Mo^{-5/48}$, a prediction supported by the numerical data for the various fluids in the $Ga$ range where path instability occurs.

Figure 4. Numerical predictions for the rise Reynolds number, plotted vs the Galilei number. The colour code and the meaning of the red bullets are similar to those of figure 3.

3.3. Flow structure

The flow structure in the base state reveals several interesting features. Figure 5 displays the azimuthal vorticity (left half of each panel) and pressure (right half) and distributions past a bubble rising in three different fluids in the $Bo$-range where the path instability threshold will later be shown to lie. With no surprise, the near-surface pressure distribution reaches its minimum very close to the equatorial plane, defined as the horizontal plane in which the longest axis of the bubble lies. This is a mere consequence of the increase of the fluid velocity along the interface from the apex of the bubble down to its equatorial plane. The norm of the azimuthal vorticity also reaches its maximum in that plane because this quantity is proportional to the product of the interface curvature in the vertical diametrical plane and the relative fluid velocity along the interface (Batchelor 1967), both of which are maximum there. This vorticity results directly from the shear-free condition obeyed by the fluid at the interface and is responsible for the existence of a boundary layer that encapsulates the bubble and turns into a wake downstream of it (Moore 1963, 1965). The wake structure is found to depend dramatically on the Bond number and, for a given $Bo$, on the fluid properties (compare panels (a,f), both for $Bo=3$). In water, all streamlines are open for $Bo=0.4$ (panel g) and they remain so up to $Bo\approx 0.5$. Beyond this point, a standing eddy exists and grows sharply with the Bond number, its length becoming of the same order as that of the long bubble axis for $Bo=0.8$ (panel i). Qualitatively similar trends are observed in the more viscous fluids. However, for a given $Bo$, the normalized rise speed $U=u_b/(gD)^{1/2}$ is significantly smaller, and so is the Weber number $We\equiv \rho u_b^2D/\gamma =Bo\,U^2$. Indeed, keeping $Bo$ unchanged, the Weber numbers in two different fluids, 1 and 2, obey the relation $We_2/We_1=(Re_2/Re_1)^2(Mo_2/Mo_1)^{1/2}$. Hence, in the nearly five times more viscous silicone oil DMS-T05 for instance, the Weber number of a bubble with $Bo=0.5$ is eight times smaller than in water, leaving it nearly spherical. Consequently, compared with water, much larger $Bo$ are required for the tangential fluid velocity and interface curvature to be large enough for a standing eddy to develop at the back of the bubble. In the case of DMS-T05, the required conditions are achieved only for $Bo\gtrsim 2.4$.

Figure 5. Base flow past bubbles rising in three different liquids: (a–c) DMST05 ($Mo=6.2\times 10^{-7}$), with, from left to right, $Bo=3, 5$ and $7$; (d–f) DMST02 ($Mo=1.6\times 10^{-8}$), with $Bo=1$, $2$ and $3$; (g–i) water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$), with $Bo=0.4$, $0.6$ and $0.8$. The left and right halves of each panel display the azimuthal vorticity and pressure distributions, respectively. The thin lines are the streamlines in the reference frame rising with the bubble.

The above information regarding the existence of a standing eddy in the various liquids is summarized in figure 6. This plot shows the critical curve separating the upper region of the $(Bo,Ga)$ plane where a standing eddy exists, from the lower region in which no such structure takes place; the two insets display an example of each situation. Present predictions are seen to be in excellent agreement with those of the fully resolved axisymmetric simulations carried out with Gerris by Cano-Lozano et al. (2013). It will become clear in the next section that, in low-$Mo$ fluids, path instability arises for $(Ga,Bo)$ pairs located below this critical curve. This finding is of special importance in that it establishes that the existence of a standing eddy is not a prerequisite for the path of a freely rising bubble to become unstable, unlike the case of a fixed bubble (Tchoufag et al. 2013). More generally, it is well established that the initial ‘seed’ leading to global wake instability past bluff bodies held fixed in a uniform stream is the growth of disturbances in the core of the standing eddy (Chomaz 2005). That path instability may happen in the absence of such a near-wake structure proves that, in the relevant $Mo$-range, the mechanism that governs this instability is not driven by wake instability per se. This confirms the findings of Cano-Lozano et al. (2016a) who established that, for $Mo\lesssim 10^{-9}$, the path of bubbles with a frozen but realistic shape becomes unstable while their wake is still stable, the low-frequency ‘aerodynamic’ mode being destabilized first.

Figure 6. Critical curve corresponding to the onset of a recirculating region at the back of the bubble in the $(Bo,Ga)$ plane. Yellow line: present results; black dashed line: predictions of Cano-Lozano, Bohorquez & Martinez-Bazán (2013). The left and right insets show some streamlines around a bubble with $Bo=0.3$ in water and a bubble with $Bo=6.5$ in DMS-T05, respectively. The thin dotted lines correspond to constant values of the Morton number, i.e. to a given liquid; $Mo$ values are specified along the upper and right sides of the figure.

4. Neutral curves

4.1. Critical bubble size

Using the procedure described in § 2, we selected a set of $Mo$ values, some of which correspond to specific liquids. For each of these values, we increased the bubble size, i.e. the Bond and Galilei numbers, until the real part of one of the eigenvalues associated with the non-axisymmetric modes $m=\pm 1$ changes sign, which determines the threshold of path instability. The resulting neutral curve, obtained by linking these thresholds, is shown in figure 7. Together with figure 8, this is presumably the most important result of the present study with respect to the original physical problem. For each fluid characterized by a Morton number in the range $10^{-13}\unicode{x2013}10^{-5}$, this curve readily provides the size of the smallest bubble whose vertical path becomes linearly unstable. For instance, the path of a bubble rising in water at a temperature of $20\,^\circ {\rm C}$ becomes unstable at $Bo=0.463$ and $Ga=250$, yielding a critical diameter $D\approx 1.854$ mm (Bonnefis 2019; Bonnefis et al. 2023). Similarly, the critical Bond numbers for Galinstan ($Mo=1.4\times 10^{-13}$) and DMS-T11 ($Mo=9.9\times 10^{-6}$) are $0.19$ and $6.85$, respectively. Therefore, the critical bubble sizes for these two fluids located close to the bounds of the $Mo$-interval considered here are $1.48$ mm and $3.875$ mm, respectively.

Figure 7. Neutral curve in the $(Bo,Ga)$ plane, obtained by connecting the (orange) symbols at which the threshold was determined. The pale and dark grey lines are the two branches of the neutral curve obtained by considering the same base state but keeping the bubble shape frozen in the stability analysis; the thin yellow dashed line is the critical curve of figure 6 beyond which a standing eddy exists at the back of the bubble in the base state. The thin dotted lines are iso-Morton-number lines corresponding to a given liquid; $Mo$ values are specified along the upper and right sides of the figure. In a given fluid, open (respectively closed) symbols indicate stable (respectively unstable) paths observed in the simulations of Cano-Lozano et al. (2016b) (, $\blacklozenge$), in the experiments of De Vries (2002) (${\odot }$, ${\bullet }$, purple) and Duineveld (1995) (${\boxdot }$, ${\blacksquare }$ purple) with ultrapure water (performed at temperatures of $28\,^\circ {\rm C}$ and $20\,^\circ {\rm C}$, respectively), and in those of Zenit & Magnaudet (2008) (,${\blacktriangle }$ purple) and Sato (2009) (,${\blacktriangledown }$ purple) with various silicone oils.

Figure 8. Variation of the reduced frequency $St=\lambda _i D/(2{\rm \pi} u_b)$ at the threshold vs the Morton number. The yellow curve was obtained by connecting the (orange) symbols at which the frequency was determined. The pale and dark grey lines correspond to the predictions along the two branches of the neutral curve predicted by the FSA using the same base state (for readability, only five of the computed points are highlighted with a marker on the dark grey line). The vertical arrow at $Mo=2.7\times 10^{-7}$ indicates the frequency jump associated with the switching from the stationary mode to the oscillatory mode predicted by FSA. $\blacklozenge$ (black): numerical predictions of Cano-Lozano et al. (2016b); ${\blacksquare }$ (purple): experimental data of Duineveld (1994) for $Bo=0.54$ in ultrapure water; ${\blacktriangle }$ (purple): experimental data of Zenit & Magnaudet (2009) for $Bo=3.92$ in DMS-T05.

In figure 7 we also report some experimental data obtained under controlled conditions, ensuring that the carrying fluid obeys a shear-free condition at the gas–liquid interface. The reference data of Duineveld (1995) in ultrapure water at $20\,^\circ {\rm C}$ yield a critical bubble diameter close to $1.82$ mm, which differs by less than $2\,\%$ from the present prediction. De Vries (2002) also determined the onset of path instability in a tank of ultrapure water with an average temperature of $28\,^\circ {\rm C}$ ($Mo=1.13\times 10^{-11}$). However, a $1.1\,^\circ {\rm C}\ {\rm m}^{-1}$ stabilizing temperature gradient was established to visualize the wake thanks to optical index gradients. The influence of this stratification on path instability is unknown and might be at the origin of the $8\,\%$ difference between the experimentally determined critical bubble size and the present prediction. Comparisons with data obtained with four different silicone oils (Zenit & Magnaudet 2008) and with three other silicone oils in which photochromic dye was added to visualize the wake (Sato 2009) prove satisfactory in the sense that the neutral curve lies generally within the ($Bo,Ga$)-interval whose lower (respectively upper) bound corresponds to the largest (respectively smallest) bubble for which a stable (respectively unstable) path could be identified experimentally. The only two exceptions correspond to a probable outlier in Sato's data for $Mo=2.35\times 10^{-7}$, and to DMS-T11 (the most viscous oil) in which the path of a bubble $4\,\%$ smaller than the critical size predicted here was found to be unstable. Numerical data obtained by Cano-Lozano et al. (2016b) with Gerris in the transition region (but not necessarily in the close vicinity of the threshold) are also reported. These data are seen to be in good agreement with present predictions for fluids with $Mo\gtrsim 10^{-8}$. In contrast, as the authors anticipated in their conclusions, these simulations underestimate the critical bubble size in fluids having a lower Morton number. The lower $Mo$ the larger the underestimate, which makes the critical diameter in water under-predicted by 15 %–20 %. Since these simulations made use of an adaptive grid with $128$ cells per bubble diameter on both sides of the interface, this significant underestimate helps appreciate how demanding fully resolved simulations of bubbly flows in low-viscosity fluids are.

In figure 7 we also reported the critical curve of figure 6 corresponding to the onset of a recirculating region at the back of the bubble in the base flow (yellow dashed line). The yellow solid and dashed lines are seen to cross each other at $Mo\approx 10^{-10}$. For lower $Mo$, the critical bubble size at which the path destabilizes first is such that no standing eddy exists in the base state. For the reasons discussed in the previous section, this leads to the conclusion that wake instability cannot be responsible for path instability in such low-$Mo$ fluids. The pale and dark grey lines in the figure also provide interesting insight into the role of time-dependent bubble deformations. To obtain the neutral curve resulting from these two lines, we considered the same base state as above, and imported the corresponding interface shape and rise speed in the stability code used by Tchoufag et al. (2014b) and Cano-Lozano et al. (2016a). In this code, the interface shape is kept frozen, i.e. the bubble can only move as a rigid body. A similar neutral curve was computed by Cano-Lozano et al. (2016a) (black line in their figure 7), using a base state obtained with Gerris, instead of the present global Newton method. Their neutral curve is similar to that based on the two grey lines in figure 7 but the critical $Ga$ they found at a given $Bo$ stands consistently below that determined here, with differences ranging from $5\,\%$ for $Bo\geqslant 5$ to $15\,\%$ for $Bo=0.5$. Again, the reasons for these differences are to be found in the limited spatial resolution of the Gerris-based computations in high-$Ga$ low-$Bo$ configurations.

Since the yellow and grey lines in figure 7 are based on strictly identical base states, time-dependent deformations not accounted for in the stability analysis used to produce the second of them are entirely responsible for the differences observed in the fate of the imposed disturbances. As the comparison of the thresholds predicted by the two approaches in a given fluid evidences, deformations always lower the threshold, i.e. they promote path instability. However, for fluids with Morton numbers less than ${\approx }5\times 10^{-7}$, the lower $Mo$ the weaker this influence. For instance, the ‘frozen-shape’ approximation, hereinafter referred to as FSA, over-predicts the critical bubble size by nearly $23\,\%$ for $Mo=10^{-7}$, but this overestimate is reduced by a factor of three in the case of Galinstan. The FSA actually predicts that the neutral curve involves two distinct modes intersecting at $Mo\approx 2.8\times 10^{-7}$. In the narrow range $4\times 10^{-7}\lesssim Mo\lesssim 6.5\times 10^{-7}$, the mode that becomes most unstable beyond this intersection (dark grey line) exhibits a S shape, implying a destabilization–restabilization scenario. In contrast, the neutral curve resulting from the stability analysis carried out with freely deforming bubbles remains single valued throughout the whole range of Morton number. For larger $Mo$, the FSA prediction based on this second mode is also quite good, with for instance a $16\,\%$ over-prediction of the critical bubble size in DMS-T11.

At this point it is also worth mentioning that Cano-Lozano et al. (2016a) established that the neutral curve predicted by FSA and that corresponding to the onset of ‘pure’ wake instability for a bubble with the same frozen shape held fixed in a uniform stream coincide in the range $10^{-9}\lesssim Mo\lesssim 6\times 10^{-7}$. Outside of this range, the threshold values of the Bond and Galilei numbers corresponding to the onset of wake instability in a given fluid are consistently larger than those at which FSA predicts the occurrence of path instability. For $Mo\lesssim 10^{-9}$, the difference between the two thresholds increases gradually as $Mo$ is decreased. For $Mo\gtrsim 6\times 10^{-7}$, the threshold of wake instability coincides with the upper neutral branch captured by FSA (here the pale grey line). In other words, depending on the Morton number, the threshold corresponding to the onset of wake instability past a fixed bubble is either identical to or larger than that of the path instability predicted by FSA. These findings, combined with the fact that figure 7 reveals that FSA overestimates the threshold of path instability whatever $Mo$, establish that the wake of a fixed bubble with a frozen shape is always stable at the actual onset of path instability. Hence, path instability cannot be explained on the sole basis of the mechanism responsible for wake instability.

4.2. Frequency at threshold

The most spectacular qualitative difference between the predictions of the present approach and those of FSA is seen in figure 8. This figure displays the Strouhal number, or reduced frequency, corresponding to the first unstable mode at the onset of the instability, computed as $St=\lambda _iD/(2{\rm \pi} u_b)$, with $\lambda _i$ the imaginary part of the unstable eigenvalue at the threshold. As the yellow curve evidences, the present approach predicts that the most unstable mode of the system is always oscillatory, in line with experimental observations. The reduced frequency increases continuously with the Morton number, ranging from $St\approx 0.02$ for $Mo=10^{-13}$ to $St\approx 0.12$ for $Mo=10^{-5}$. The FSA predicts that the first non-vertical path of the bubble is associated with ‘low-frequency’ oscillations in low-$Mo$ liquids ($Mo\lesssim 7\times 10^{-11}$) and with ‘high-frequency’ oscillations in liquids with $Mo\gtrsim 3\times 10^{-7}$, in agreement with the conclusions of Cano-Lozano et al. (2016a). Present findings with deformable bubbles are qualitatively consistent with these predictions in both ranges. The reduced frequencies determined in both approaches for low-$Mo$ liquids are even in fairly good quantitative agreement. This confirms that deformations do not bring major changes in the dynamics of isolated bubbles rising in such liquids, which extends to all low-$Mo$ fluids the conclusions drawn by Bonnefis et al. (2023) in the specific case of water. More surprisingly, the agreement is still reasonable in fluids with $Mo\gtrsim 3\times 10^{-7}$, where the FSA neutral curve (now corresponding to the dark grey line) underestimates the reduced frequency by only roughly $20\,\%$. In contrast, the two approaches dramatically disagree in the intermediate $Mo$-range, where FSA predicts that the most unstable mode (pale grey line) is stationary. Hence, one has to conclude that, in the linear framework adopted here, time-dependent deformations promote the instability of the oscillatory mode and succeed in making it more unstable than the stationary mode in this intermediate range, while only the latter mode may become unstable if these deformations are ignored. Some experimental and numerical data are also reported in figure 8. Most of them were obtained significantly above the threshold, which explains why the corresponding frequencies lie above present predictions (for instance, the bubble size is $8\,\%$ beyond the threshold value in the two experimental determinations). Only the two numerical determinations of Cano-Lozano et al. (2016b) in silicone oils DMS-T02 ($Mo=1.6\times 10^{-8}$) and DMS-T11 ($Mo=9.9\times 10^{-6}$) correspond to near-threshold conditions according to figure 7. The corresponding two $St$ values are seen to agree very well with present predictions.

Having determined the normalized rise speed in the base state, we can re-plot the neutral curve of figure 7 in two different forms, to make the variations of the critical Reynolds number, $Re=Ga\,U$, and Weber number, $We=Bo\,U^2$, with the liquid properties apparent. Figure 9(a) shows that the critical Reynolds number decreases sharply as $Mo$ increases, from ${\approx }2070$ for $Mo=10^{-13}$ to $70$ for $Mo=10^{-5}$, via $Re=668$ for water at $20\,^\circ {\rm C}$. Conversely, the critical Weber number (figure 7b) is seen to increase gradually from $3.20$ for $Mo=10^{-13}$ to $6.36$ for $Mo=10^{-5}$, via $We=3.30$ for water. The above values for water are to be compared with those reported by Duineveld (1995), namely $Re=658$, $We=3.275$, from which they only differ by $1.5\,\%$ and $0.75\,\%$, respectively.

Figure 9. Influence of the fluid properties on (a) the critical Reynolds number and (b) the critical Weber number.

5. Unstable modes: order of occurrence and physical nature

In this section, we examine the variations of the most unstable eigenvalues of the coupled bubble–fluid system as the size of the bubble (here measured through the Bond number) is varied in the vicinity of the primary threshold. We also analyse the spatial structure of the associated global modes. It will soon become apparent that the order in which the unstable modes follow one another as the Bond number increases, and even in some cases the way the nature of a given mode varies, depends dramatically on the range of fluid properties under consideration, the complexity of the picture increasing as the Morton number is decreased. For this reason, we discuss these features in descending order of $Mo$.

5.1. Silicone oil DMS-T05 ($Mo=6.2\times 10^{-7}$)

5.1.1. First unstable eigenvalues: thresholds and nature of the associated modes

Figure 10 displays the variations of the first five eigenvalues whose real part becomes positive as the Bond number, hence the size, of a bubble rising in silicone oil DMS-T05 is increased. All but one of these eigenvalues have a non-zero imaginary part, indicating that the corresponding instabilities arise through Hopf bifurcations. The threshold of path instability is encountered at $Bo=3.4$. This first unstable mode, associated with azimuthal wavenumbers $|m|=1$, has a dimensionless frequency $\lambda _i/2{\rm \pi} \approx 0.115$ at the threshold (throughout this section, eigenvalues are normalized using the gravitational time scale $(D/g)^{1/2}$). As panel (c) suggests, this mode is essentially associated with a lateral drift of the bubble centroid accompanied by some rigid-body rotation. Time-dependent deformations are small: it will be shown in § 6 that their maximum magnitude is barely $3\,\%$ of that of the horizontal displacement of the bubble centroid. For this reason, we refer to this mode as ‘rigid’.

Figure 10. Variation of the first five unstable eigenvalues with the Bond number for bubbles rising in DMS-T05 ($Mo=6.2\times 10^{-7}$). (a) Growth rate ($\lambda _r$), normalized by the gravitational time $(D/g)^{1/2}$; (b) radian frequency ($\lambda _i$) normalized similarly; (c) equilibrium shapes at the threshold (labels refer to the transition points marked with bullets in a). In (a,b), dashed (respectively solid) lines represent the stable (respectively unstable) part of the various branches. In (c), the bubble contour in the vertical diametrical plane is shown at the threshold of each of the five successive modes in the base state (black solid line) and, with an arbitrary amplitude of the disturbance, for $t=0$ (red dashed line) and $t={\rm \pi} /(2\lambda _i)$ (blue dashed line); bubble centroids are arbitrarily shifted in the $z$ direction for readability.

We continue to track this mode and other possible ones beyond this primary threshold by computing the base flow, still assumed axisymmetric, at the relevant Bond number. The growth rate of the rigid oscillatory mode continues to increase with $Bo$, and a second mode becomes unstable at $Bo\approx 5.5$ (point $B$ in figure 10a). This mode is axisymmetric and oscillatory, with a dimensionless frequency approximately $5.5$ times higher than that of the rigid mode. Examination of the corresponding successive contours in panel (c) reveals that this mode is associated with shape oscillations known as $l=2,m=0$, or $(2,0)$, in the terminology of spherical harmonics, i.e. the associated interface displacements leave the interface position unchanged at two angular positions located on both sides of the equatorial plane. We shall show in § 6 that these deformations have a magnitude of the same order as the vertical displacement of the bubble centroid with respect to its reference position in the base state. For these reasons, we refer to this mode and those sharing similar properties as ‘shape’ modes. By further increasing $Bo$, a third unstable mode is encountered at point $C$ ($Bo\approx 5.9$). Its growth rate increases strongly beyond the threshold so that this mode would quickly become dominant, were its amplitude similar to that of the primary oscillatory mode. Panel (b) reveals that this mode is stationary and, since its azimuthal wavenumber is $|m|=1$, it is also asymmetric. These two properties imply that this mode induces an inclined path of the bubble centroid. Moreover, panel (c) indicates that the corresponding interface deformations are weak compared with those associated with the previous ‘shape’ mode, so that this stationary mode can be considered ‘rigid’, just like the primary mode. Therefore, two rigid modes coexist in the system beyond $Bo\approx 5.9$, but one is oscillatory while the other is stationary. To the best of our knowledge, steady oblique bubble paths have never been reported in experiments, suggesting that this second rigid mode never becomes dominant. Nonlinearities not being taken into account here, i.e. the fact that the actual base flow at point $C$ is no longer axisymmetric due to the nonlinear corrections brought by the primary oscillatory mode and is presumably even not stationary owing to the shape oscillations associated with the $(2,0)$ shape mode, are likely the reason for this.

A fourth mode becomes unstable at point $D$ ($Bo\approx 6.6$), and restabilizes itself within the narrow interval $8.65\lesssim Bo\lesssim 9.45$ before becoming unstable again at larger $Bo$. This mode has much in common with the $(2,0)$ shape mode, except that it is not axisymmetric. Rather, it is associated with azimuthal wavenumbers $|m|=2$ and, in the $Bo$-range under consideration, its frequency is nearly half that of the $(2,0)$ mode. That the frequencies of two modes with the same $l$ but different $m$ differ in this context is no surprise, since spherical harmonics are no longer the eigenmodes with which a bubble oscillates when its undisturbed shape is not spherical (Meiron 1989). Here, the Bond number is large, and it is only in the small-$Bo$ limit that the two frequencies are expected to get close to each other. Due to the symmetries preserved by angular variations of the form ${\rm e}^{\pm 2{\rm i}\theta }$, this $(2,2)$ mode does not induce any lateral displacement of the bubble. Last, at point $E$ ($Bo\approx 9.45$), a third shape mode with a frequency close to that of the $(2,2)$ mode becomes unstable. Deformations associated with this mode share the same characteristics as those induced by modes $(2,0)$ and $(2,2)$ on both sides of the bubble equatorial plane. In contrast, as panel (c) makes clear, the two halves of the bubble deform in an asymmetric manner at every instant of time in the vertical diametrical plane, one half becoming more pointed while the other becomes more rounded. These characteristics identify this mode as the $(2,1)$ shape mode. Unlike those associated with $m=0$ and $|m|=2$, this mode may contribute to the lateral drift of the bubble at large enough Bond numbers, owing to its azimuthal asymmetry. It must be stressed that the full sequence described above is kept unchanged when the Morton number is increased up to $10^{-5}$.

5.1.2. Spatial structure of the first unstable mode

Figure 11 shows the spatial structure of the first unstable mode slightly above the threshold. The real and imaginary parts of a given mode may be thought of as the associated disturbance at two different instants of time a quarter of a period apart, i.e. shifted by an interval $\tau =({{\rm \pi} }/{2})\lambda _i^{-1}$. Therefore, if the bubble is zigzagging in a vertical plane (as it does if the two disturbances associated with modes $m=+1$ and $m=-1$ have the same amplitude), the real part corresponds to the time instant by which the bubble reaches its maximal lateral excursion (red contours in figures 11(a) and 10(c)), while the imaginary part corresponds to that by which it crosses the midline of its path (blue contours in figures 11(b) and 10(c)). The streamlines pattern evidences the antisymmetric nature of the disturbance flow. An alternation of positive and negative vorticity zones may be noticed in the wake, both in every horizontal plane (varying $r$) and at every radial position from the axis of the base flow (varying $z$). This pattern is indicative of a vortex shedding process. A specific structure in which fluid particles rotate clockwise on both sides of the axis of the base flow is visible near the bottom of panel (b); in a three-dimensional representation, this structure would look like a pair of crescents connected through their horns in the vertical plane $\theta =\pm {\rm \pi}/2$, with fluid particles rotating in the same direction in every azimuthal plane. As panels (c,d) make clear, this structure is actually the first of a series that develops further downstream in the wake, with alternating positive and negative directions of rotation.

Figure 11. Structure of the first unstable mode past a bubble rising in DMS-T05 at $Bo=3.7$. (a,c) Real part of the mode near the bubble and further downstream in the wake, respectively; (b,d) same for the imaginary part. The left and right halves in (a,b) display the pressure and azimuthal vorticity iso-levels, respectively; those in (c,d) display the azimuthal velocity and vorticity iso-levels, respectively; some streamlines defined in the reference frame of the base configuration are also shown. Black, red and blue bubble contours correspond to the base state, and the real and imaginary parts of the interface disturbance ($\hat \eta$), respectively.

5.2. Silicone oil DMS-T02 ($Mo=1.6\times 10^{-8}$)

We now move to silicone oil DMS-T02 which is $2.6$ times less viscous than DMS-T05 and has almost the same surface tension. Figure 12 shows how the first five eigenvalues yielding unstable modes vary with the Bond number in this case. The beginning of the sequence is similar to that described for DMS-T05. That is, path instability first arises through an oscillatory ‘rigid’ mode that becomes unstable at $Bo=1.56$ (point $A$ in panel a). Then the shape mode $(2,0)$ becomes unstable at point $B$ ($Bo=2.2$), followed by the stationary rigid mode at point $C$ ($Bo=2.43$). Further increasing $Bo$, the $(2,1)$ and $(2,2)$ shape modes become unstable at points $E$ ($Bo=3.89$) and $D$ ($Bo=4.12$), respectively. Note that, compared with the case of DMS-T05, these last two modes emerge in reverse order.

Figure 12. Same as figure 10 for bubbles rising in DMS-T02 ($Mo=1.6\times 10^{-8}$).

Nevertheless, the main originality of the present bifurcation diagram compared with that of DMS-T05 is that the growth rate of the primary mode starts to decrease slightly beyond point $B$, and becomes even very weakly negative within a tiny interval around $Bo\approx 2.55$, before turning positive again and growing continuously up to the upper bound of the explored Bond number range. This non-monotonic behaviour is to be related to the fact that DMS-T02 stands in the middle of the $Mo$-range where FSA predicts the emergence of path instability through a stationary bifurcation instead of a Hopf bifurcation, being unable to detect that the oscillatory rigid mode becomes actually unstable first. This failure is presumably closely connected with the fact that, although interface deformations succeed in making the oscillatory mode unstable first, the growth rate of this mode remains quite low throughout the interval bounded by points $A$ and $C$. We shall come back to this point in § 7.

5.3. Water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$)

5.3.1. First unstable eigenvalues: thresholds and nature of the associated modes

We finally consider the case of bubbles rising in pure water maintained at a temperature of $20\,^\circ {\rm C}$. Variations of the eigenvalues leading to the first five unstable modes are shown in figure 13. Again, the sequence begins with the emergence of an oscillatory rigid mode that becomes unstable at $Bo=0.463$ with a dimensionless frequency $\lambda _i/2{\rm \pi} =0.088$, followed by the $(2,0)$ shape mode at $Bo=0.525$. However, the next steps differ deeply from those encountered in the previous two cases. The primary mode grows continuously until $Bo\approx 0.64$. There, this mode splits into two separate branches, both of which are stationary. Therefore, at the point (marked with a square in the figure) where the three branches meet, the pair of complex eigenvalues associated with the oscillatory mode turns into a pair of real eigenvalues. Such a codimension 2 singular point is referred to as ‘exceptional’ in the theory of non-Hermitian Hamiltonian systems (Bergholtz, Budich & Kuns 2021). In the present case, this exceptional point, at which the two eigenvalues change their nature, corresponds to a Takens–Bogdanov bifurcation with ${O}(2)$-symmetry (Dangelmayr & Knobloch 1987). Beyond this point, the upper stationary branch exhibits a rapid and continuous growth. Conversely, the growth rate of the lower branch decreases sharply with the distance to the exceptional point, until it becomes negative at $Bo=0.685$ (point $C$ in the figure). Similar to the cases of DMS-T02 and DMS-T05, and undoubtedly for the same reasons, inclined paths corresponding to the above stationary modes have not been observed in experiments carried out in water, irrespective of the bubble size.

Figure 13. Same as figure 10 for bubbles rising in water at $20\,^\circ {\rm C}$ ($Mo=2.54\times 10^{-11}$). In (a), the square indicates the codimension-two ‘exceptional’ point at which the pair of complex conjugate eigenvalues associated with the first unstable mode turns into a pair of real eigenvalues.

Next, figure 13 indicates that the mode $(2,2)$ becomes unstable at $Bo=0.68$ (point $D$ in the figure) and then grows slowly with the Bond number. Then, a fourth mode becomes unstable at $Bo=0.809$ (point $A'$) with a dimensionless frequency $\lambda _i/2{\rm \pi} = 0.162$. This mode (identified as rigid osc. II in the legend of figure 13a) is associated with an azimuthal wavenumber $|m|=1$, exhibits low-frequency oscillations (panel b), and is not accompanied by significant shape oscillations (panel c). These characteristics point to another rigid mode. Hence, similar to the right part of figures 10 and 12, two rigid modes coexist in the case of water beyond point $A'$. In all three cases, the most amplified of them is stationary (yellow curve in figure 13) while the other (red curve) exhibits slow oscillations. Finally, figure 13 indicates that a fifth eigenvalue crosses the real axis at $Bo\approx 0.88$ (point $E$). As panels (b,c) make clear, this eigenvalue corresponds to the asymmetric shape mode $(2,1)$ already encountered with DMS-T02 and DMS-T05. One may notice that the mode $(2,1)$ becomes unstable at a larger $Bo$ than the mode $(2,2)$ in the cases of water and DMS-T05, while the reverse is true for DMS-T02. It is also worth noting that, among the three shape modes revealed by figure 13, the mode $(2,2)$ is the one with the lowest frequency in the range of Bond numbers relevant here, say $Bo\gtrsim 0.4$, followed by the mode $(2,0)$ and then the mode $(2,1)$. This is in line with the conclusions obtained by Meiron (1989) assuming an inviscid potential flow. In contrast, figures 10 and 12 show that modes $(2,1)$ and $(2,2)$ have nearly equal frequencies, both significantly lower than that of the mode $(2,0)$, in DMS-T05 and DMS-T02. This suggests that vortical effects not accounted for in the above reference significantly modify the dynamics of shape oscillations of strongly deformed bubbles, even in low-viscosity fluids. Indeed, although the viscosity of DMS-T02 is less than twice that of water, its surface tension is four times lower, so that the relevant Bond numbers in DMS-T02 are typically larger than $3.0$, while those in water are less than $1.0$.

To finish with figure 13, it is important to stress that the succession of modes it reveals is the canonical sequence encountered throughout the low-$Mo$-range, say for $Mo\lesssim 10^{-9}$. For instance, the same first five modes become unstable in the same order in the case of DMS-T00 ($Mo=1.8\times 10^{-10}$), with the low-frequency rigid mode splitting into two stationary modes at an exceptional point, just as in figure 13. Only the thresholds and frequencies differ, with points $A$–$B$ then located at $Bo=0.66$ and $0.805$, followed by points $C$ and $D$ almost coincident at $Bo=0.98$, and finally points $A'$ and $E$ located at $Bo=1.135$ and $1.42$, respectively.

5.3.2. Spatial structure of rigid and shape modes

Figure 14 displays the structure of the three successive rigid modes encountered in the sequence analysed above. Panels (a,b) look qualitatively similar to their counterparts in figure 11. However, the ‘double-crescent’ structure identified in the bottom part of the latter is not present here. The reason is that the Strouhal number of the primary low-frequency mode is approximately three times lower in water (see figure 8). For this reason, the vertical distance separating vortices shed in the wake is much larger, so that the double-crescent structure closest to the bubble is located further downstream. Panels (c,e) compare the wake structure of the primary and secondary oscillatory rigid modes. They evidence the higher vortex shedding frequency of the latter, in line with the nearly twice as large dimensionless frequency noticed at point $A'$ in figure 13 compared with that at point $A$. In the case of the stationary mode displayed in panel (d), the azimuthal vorticity and velocity disturbances are seen to keep a constant sign all along the wake (the real part of the azimuthal velocity is null in this case, which is why the imaginary part is shown in the figure). This is the generic hallmark of wakes associated with stationary inclined paths of axisymmetric bodies; e.g. Tchoufag et al. (2014a, b). No vortex shedding takes place in that situation, the inclination of the body path resulting from the constant lift force generated by a pair of semi-infinite counter-rotating streamwise vortices.

Figure 14. Structure of the rigid modes past a bubble rising in water. (a,b) Primary low-frequency mode at $Bo=0.48$ in the bubble vicinity, with the real and imaginary parts shown in (a,b), respectively. The left and right halves of the two panels display the pressure and azimuthal vorticity iso-levels, respectively; some streamlines defined in the reference frame of the base configuration are also shown. Black, red and blue bubble contours correspond to the base state, and the real and imaginary parts of the interface disturbance, respectively. (c–e) Wake structure of the successive rigid modes. (c) Primary low-frequency mode at $Bo=0.48$; (d) most amplified stationary mode at $Bo=0.75$; (e) secondary oscillating mode at $Bo=0.85$. The right half of each panel displays the real part of the azimuthal vorticity iso-levels; the left half of (c,e) (respectively d) displays the real (respectively imaginary) part of the azimuthal velocity iso-levels.

Figure 15 reveals the structure of the $(2,0)$ and $(2,1)$ shape modes at $Bo=0.54$ and $Bo=0.9$, respectively. The axial symmetry of the flow field is obvious in panels (a,b), with a succession of toroidal eddies of alternate signs shed downstream of the bubble. The red and blue contours evidence the fact that the bubble alternatively ‘inflates’ and ‘deflates’ along its axis while simultaneously shortening or lengthening in its equatorial plane. The streamlines in panels (c,d) make the antisymmetry of the mode $(2,1)$ just as obvious. The wake is dominated by a series of double-crescent structures qualitatively similar to those identified in figure 11. The wavelength separating two consecutive structures is slightly shorter than that of the toroidal eddies in panels (a,b), although the dimensionless frequency of the mode $(2,1)$ at $Bo=0.9$ is approximately $30\,\%$ lower than that of the mode $(2,0)$ at $Bo=0.54$ according to figure 13(b). However, these structures are advected downstream by the base flow, so that their wavelength is determined by the Strouhal number $St=\lambda _i/(2{\rm \pi} U)$, not directly by $\lambda _i$. Therefore, what the comparison of the wavelengths for the two modes indicates is that the normalized rise speed $U$ decreases by slightly more than $30\,\%$ from $Bo=0.54$ to $Bo=0.9$, owing to the progressive flattening of the bubble (in the regime where bubbles keep an approximate oblate spheroidal shape, $U$ is known to vary approximately as $Bo^{-1/2}$ for $Bo\ll 1$ and $Bo^{-1}$ for $Bo\gg 1$; see equation (7-3) and the associated figure in Clift, Grace & Weber 1978).

Figure 15. Structure of the first axisymmetric and asymmetric shape modes past a bubble rising in water. (a,b) Axisymmetric $(2,0)$ mode at $Bo=0.54$, with the real and imaginary parts shown in (a,b), respectively. (c,d) Same for the asymmetric $(2,1)$ mode at $Bo=0.90$. For caption see figures 11 and 14.

6. Respective magnitudes of rigid-body motions and time-dependent deformations

6.1. Decomposition technique

In the previous section, we identified ‘rigid’ and ‘shape’ modes qualitatively, stating that the former leave the bubble shape almost unchanged while the latter are associated with significant time-dependent shape oscillations. However, this distinction deserves a more quantitative basis. Moreover, determining how large the deformations at the threshold of path instability are, compared with the lateral drift of the bubble, is key to clarifying the relative role of the various possible causes of instability present in the system. In order to quantify time-dependent deformations, we define a $(x,y,z)$ Cartesian coordinate system, the origin of which stands at the centroid of the bubble in the base state. We assume that the bubble performs planar zigzags in the vertical $(x,z)$ plane corresponding to the plane $\theta =(0,{\rm \pi} )$ in the cylindrical $(r,\theta,z)$ system defined in § 2. At any position ${\boldsymbol {x}}_0$ on the undisturbed bubble–fluid interface, the normal displacement of the interface, $\hat {\eta }({\boldsymbol {x}}_0)$, may be decomposed into a rigid motion inducing a lateral (respectively vertical) displacement $\hat {T}_x$ (respectively $\hat {T}_z$) about the $x$ (respectively $z$) axis and an inclination $\hat {\psi }$ resulting from a rotation about the $y$ axis, augmented by a local volume-preserving deformation $\hat {\zeta }({\boldsymbol {x}}_0)$ in the normal direction. Projecting this decomposition onto the local normal unit vector ${\boldsymbol {n}}_0$, one has

(6.1)\begin{equation} \hat{\eta}({\boldsymbol{x}}_0)=\hat{T}_x({\boldsymbol{e}}_x\boldsymbol{\cdot}{\boldsymbol{n}}_0)+ \hat{T}_z({\boldsymbol{e}}_z\boldsymbol{\cdot}{\boldsymbol{n}}_0)+\hat{\psi}({\boldsymbol{x}}_0\times{\boldsymbol{n}}_0) \boldsymbol{\cdot}{\boldsymbol{e}}_y+\hat{\zeta}({\boldsymbol{x}}_0), \end{equation}

with ${\boldsymbol {e}}_x,{\boldsymbol {e}}_y,{\boldsymbol {e}}_z$ the unit vectors in the $x,y,z$ directions, respectively. Odd modes $m=\pm 1$ do not induce any translation in the vertical direction, so that $\hat {T}_z=0$ for such modes. Conversely, the lateral displacement $\hat {T}_x$ and the inclination $\hat {\psi }$ of the bubble remain null for even modes $m=-2,0,2$. To determine the non-zero displacement and inclination associated with a given mode, we make use of a least-squares fitting technique. That is, we define the local functional $\mathcal {F}_{\hat {T}_x,\hat {T}_z,\hat {\psi }}({\boldsymbol {x}}_0)\equiv \{\hat {\eta }-[\hat {T}_x ({\boldsymbol {e}}_x\boldsymbol {\cdot }{\boldsymbol {n}}_0)+\hat {T}_z({\boldsymbol {e}}_z \boldsymbol {\cdot }{\boldsymbol {n}}_0)+\hat {\psi }({\boldsymbol {x}}_0\times {\boldsymbol {n}}_0) \boldsymbol {\cdot }{\boldsymbol {e}}_y]\}^2$, integrate it over the entire bubble contour and minimize the corresponding residue with respect to the non-zero components of the triplet $(\hat {T}_x,\hat {T}_z,\hat \psi )$. Given the linear framework of the present study, the magnitude of the displacements is arbitrary, and only their relative magnitudes with respect to a common reference makes sense. This is why we normalize the local displacement $\hat {\eta }({\boldsymbol {x}}_0)$ by $\hat {T}_x$ for odd modes and $\hat {T}_z$ for even modes, in order to make the translational displacement of the bubble centroid unity in all cases.

Thanks to this procedure, the evolution of the bubble shape and the relative magnitude of the interface displacements induced by the rigid-body rotation and volume-preserving deformations may be tracked along each branch of the bifurcation diagram. To save space, we do not discuss shape evolutions here (examples are provided by Bonnefis 2019), and focus on the relative variations of the different contributions and their physical origin.

6.2. Variations of the inclination/deformation-induced contributions and path styles

Figure 16(a) considers the case of the most amplified rigid mode for a bubble rising in water at $20\,^\circ {\rm C}$, i.e. the low-frequency oscillating mode up to $Bo\approx 0.64$, followed by the upper branch of the stationary mode for higher $Bo$ (figure 13). It shows how the maximum over the bubble contour of the three contributions to $\hat {\eta }$ varies with the Bond number. In the vicinity of the threshold ($Bo=0.463$), the inclination-induced displacement is seen to be essentially in quadrature with the lateral drift, which is the expected behaviour in a planar zigzagging motion if the bubble axis remains aligned with the path (i.e. its inclination is maximum when the bubble reaches the centreline of its zigzagging path). The relative magnitude of the maximum normal displacement resulting from this inclination is $7\times 10^{-2}$ and experiences little variation up to $Bo\approx 0.55$. At the threshold, the maximum of the volume-preserving deformation is almost in phase with the lateral drift. This is because the rise speed reaches a maximum at each extremity of the zigzag (Mougin & Magnaudet 2006), resulting in some flattening of the bubble. Conversely, the rise speed is minimum when the bubble crosses the centreline of the zigzag, which slightly reduces its oblateness, yielding the small non-zero imaginary component of the deformation (dashed green line in figure 16a). The relative magnitude of the deformation does not exceed $3\times 10^{-3}$ in this $Bo$ range. This provides the confirmation that time-dependent deformations are extremely weak at the onset of path instability for bubbles rising in low-$Mo$ liquids. These findings also validate the ‘rigid’ qualification employed for the primary low-frequency mode throughout § 5. Conversely, the relative magnitude of deformations in modes $(2,0)$ and $(2,1)$ (not shown) is almost unity at their respective thresholds. The real (respectively imaginary) part of the displacements associated with the bubble inclination (respectively deformation) is typically one order of magnitude smaller than its imaginary (respectively real) counterpart. These weak components exhibit abrupt variations close to the threshold (and for one of them around $Bo=0.6$), presumably because such low-magnitude quantities are only determined with a limited accuracy by the least-squares technique.

Figure 16. Decomposition of the normal displacement of the bubble surface into a horizontal uniform translation, an inclination-induced displacement resulting from a rigid-body rotation, and a volume-preserving deformation. (a) Variation with $Bo$ for the most amplified rigid mode $|m|=1$ in the case of a bubble rising in water at $20\,^\circ {\rm C}$; the yellow vertical line indicates the threshold of path instability; (b) variation with the fluid properties at the path instability threshold. Each line figures the maximum of the corresponding component over the bubble contour; for each component, the solid and dashed lines refer to the real and imaginary parts, respectively.

These two weak components increase sharply with $Bo$ and become dominant for $Bo\gtrsim 0.55$, marking a change in the style of the bubble ascent. In particular, the real part of the displacement associated with the bubble inclination exceeds the imaginary part for $Bo\geqslant 0.59$. This behaviour indicates that the bubble now exhibits a non-zero inclination at the extremities of its zigzagging path, so that its minor axis is no longer aligned with its velocity. In between two successive extremities of the zigzag, the bubble now glides along its path rather than facing it. This makes the part of its surface pointing in the direction of the lateral drift slightly more pointed, while the opposite part becomes slightly more rounded. This is the reason why the imaginary part of the deformation is significant in this second style of zigzagging motion. Conversely, owing to the non-zero bubble inclination at the extremities of the zigzag, the rise speed is slightly lower than in the previous zigzag style. This mitigates the flattening of the bubble, thus reducing the real part of the volume-preserving deformation, as observed on the solid green line of figure 16(a) in the range $0.5\lesssim Bo\lesssim 0.6$.

At $Bo=0.64$, the imaginary parts of the rotation- and deformation-induced contributions fall to zero, a consequence of the change in nature of the primary rigid mode at the exceptional point revealed by figure 13. Beyond this point, the rigid mode is stationary, yielding a non-oscillatory inclined path. When the Bond number increases in the range 0.6–0.8, the bubble gradually flattens and the size of the standing eddy at its back grows significantly, as panels (h,i) of figure 5 show. This makes the strength of the trailing vortices that set in at the path instability threshold increase with $Bo$, which translates directly into an increase of the lift force acting on the bubble, hence an increase of its inclination with the Bond number. This scenario is confirmed by the solid red line in figure 16(a), which shows that the maximum displacement associated with the bubble inclination increases by a factor of four from $Bo=0.64$ to $Bo=0.8$. The larger the bubble inclination, the lower its rise speed, which in turn makes its oblateness reduce compared with that in the base state. This is the origin of the sharp increase of the real part of the volume-preserving deformation (solid green line in figure 16a), which gains two orders of magnitude between $Bo=0.6$ and $Bo=0.8$. Nevertheless, deformations remain small compared with the lateral drift, reaching only a relative magnitude of $3\times 10^{-2}$ at $Bo=0.8$.

Figure 16(b) shows how the relative magnitude of the three contributions to $\hat {\eta }$ vary with the Morton number at the path instability threshold, i.e. along the neutral curve of figure 7. Since the corresponding unstable mode is always oscillatory, one expects the relative magnitude of the different contributions to behave qualitatively in the same way as in the case of water in the vicinity of the primary threshold. Indeed, in line with the behaviours noticed above, deformation-induced displacements are essentially in phase with the bubble lateral drift throughout the eight decades of $Mo$ explored here, while those associated with the bubble inclination are almost in quadrature. Again, we attribute the non-smooth variations of the weak real (respectively imaginary) component of the inclination-induced (respectively deformation-induced) contribution observed for $Mo\lesssim 10^{-10}$ to the limitations of the minimization technique. In contrast, the dominant component of these two contributions is seen to increase smoothly with the Morton number. The relative deformation-induced displacements, say $\hat {\eta }^{def}_{max}|_0$, follow the approximate power law $\hat {\eta }^{def}_{max}|_0\sim Mo^{0.23}$, increasing by two orders of magnitude from Galinstan ($\hat {\eta }^{def}_{max}|_0\approx 10^{-3}$) to DMS-T11 ($\hat {\eta }^{def}_{max}|_0\approx 7\times 10^{-2}$). The relative inclination-induced displacements only grow by one order of magnitude, from $4\times 10^{-2}$ to $5\times 10^{-1}$, in between the same two limits. Nevertheless, the fact that these displacements reach $50\,\%$ of the lateral drift for DMS-T11 (to be compared with $7\,\%$ in the case of water) indicates that bubbles wobble significantly when their path becomes unstable in liquids having a Morton number four to five orders of magnitude higher than that of water.

7. Summary and final discussion

7.1. Main findings

Thanks to an innovative numerical approach allowing the linear stability of the viscous flow past freely rising and deforming bubbles to be assessed accurately, even when the rise Reynolds number is large, we explored the conditions under which the path of such bubbles deviates from a straight vertical line and the physical nature of the underlying instability modes over a wide range of properties of the carrying liquids. Present results extend over eight orders of magnitude of the Morton number. Over this range, the critical Bond number increases by a factor of $35$ from Galinstan to silicone oil DMS-T11 (yielding an increase by a factor of $2.6$ of the critical bubble size), whereas the critical Reynolds number decreases by a factor of $30$, from ${\approx}\,2100$ to $70$, in between the same two liquids. We showed that path instability always arises through the destabilization of a non-axisymmetric oscillatory mode associated with weak time-dependent changes of the bubble shape. Increasing the bubble size in a given fluid, the nature of this primary ‘rigid’ mode stays unchanged if the Morton number is large enough, say $Mo\gtrsim 10^{-7}$. Conversely, this mode gives birth to two stationary modes if the Morton number is low enough ($Mo\lesssim 10^{-9}$). In the intermediate range $10^{-9}\lesssim Mo\lesssim 10^{-7}$, this mode remains oscillatory even far from the primary threshold but its growth rate varies non-monotonically with the bubble size, and it may even re-stabilize within a narrow size range. A second rigid mode, which turns to be stationary for intermediate and ‘high’ Morton numbers but oscillatory for low Morton numbers also arises further away from the primary threshold. ‘Shape’ modes distinct from the previous rigid modes in that they exhibit much larger oscillations of the bubble surface also become unstable as the bubble size is increased.

Comparison of present predictions with reference experimental and numerical data revealed a very good agreement, both on the threshold (most often measured through the Bond or Galilei number, or alternatively the rise Reynolds number) and the reduced frequency of path oscillations (measured through a Strouhal number). This agreement fully supports the view that the path instability of real gas bubbles results from a linear mechanism associated with a Hopf bifurcation. In contrast, oblique paths which are the hallmark of stationary modes have not been reported in experiments, nor in fully resolved simulations. We see their absence as an indication that the changes introduced in the base flow by the primary non-axisymmetric oscillatory mode are significant enough to alter the next steps of the actual bifurcation sequence compared with the predictions of the linear approach.

To better quantify time-dependent bubble deformations and their influence on the path instability mechanisms and characteristics, we employed two additional and complementary tools. First, a least-squares minimization technique allowed us to split the displacements of the bubble–fluid interface into contributions associated with a rigid-body motion and volume-preserving deformations. By doing so, we could among other things establish how the amplitude of the deformations varies with the Morton number at the onset of path instability. It turned out that the ratio of the maximum deformation to the maximum lateral drift of the bubble centroid is always much smaller than unity. In particular, this ratio is approximately $10^{-3}$ for Galinstan and $3\times 10^{-3}$ for water. It exceeds $1\,\%$ only for liquids with $Mo\gtrsim 1.5\times 10^{-8}$, reaching $7\,\%$ for DMS-T11 ($Mo=9.9\times 10^{-6}$). Second, we took advantage of the base states computed with the present numerical technique to update the predictions obtained in the framework of the FSA by Cano-Lozano et al. (2016a). With this procedure, any difference in the nature of the primary bifurcation or/and the corresponding threshold may be ascribed to effects of time-dependent bubble deformations. We found that, despite their weak magnitude, these deformations lower significantly the threshold whatever the Morton number. Hence, a general conclusion of the present study is that deformations always promote path instability in the linear framework. Differences between the predicted critical size of deformable and non-deformable bubbles increase gradually from $8\,\%$ at $Mo\approx 10^{-13}$ to $23\,\%$ at $Mo\approx 6\times 10^{-7}$, before decreasing sharply to nearly $15\,\%$ for larger $Mo$. For $Mo\lesssim 7\times 10^{-11}$, FSA correctly predicts that the most unstable mode is a ‘low-frequency’ one with $St={O}(0.01)$ at the threshold. Similarly, it rightly predicts that path instability arises through a ‘high-frequency’ mode with $St={O}(0.1)$ at the threshold for $Mo\gtrsim 3\times 10^{-7}$. The predicted reduced frequencies lie $20$ to $30\,\%$ below those obtained with deformable bubbles, but their variations with the Morton number are qualitatively similar. In contrast, FSA dramatically fails to detect that the most unstable mode is still oscillatory in the intermediate-Morton-number range $7\times 10^{-11}\lesssim Mo\lesssim 3\times 10^{-7}$.

7.2. Discussion

To better understand the above issue and gain some insight into the subtle balances involved in the various regimes encountered by varying the fluid properties, it is useful to put the first stages of the bifurcation diagrams obtained with deformable bubbles in perspective with those determined by Tchoufag et al. (2014b) who assumed a perfectly oblate spheroidal bubble with a frozen shape. By increasing the geometrical aspect ratio $\chi$ in the range $[2.21,2.30]$ they identified three successive phenomenologies summarized in their figures 1 and 2. For $2.21\leqslant \chi \leqslant 2.23$, the system first becomes unstable through a low-frequency mode which, by increasing the bubble size (i.e. the Reynolds number in their case), turns into a pair of stationary modes, one of which is quickly damped as the bubble size is increased further while the other gets strongly amplified (their figure 1). Present findings summarized in figure 13 indicate that this is exactly what happens with deformable bubbles rising in water and more generally in low-$Mo$ liquids. That the scenario relevant to frozen spheroidal bubbles remains unchanged in the case of freely deforming bubbles rising in low-$Mo$ liquids strongly supports the view that time-dependent deformations do not play any causal role in the mechanisms governing the first stages of path instability in such fluids, although they lower the critical bubble size by $18\,\%$ for water and increase the reduced frequency by nearly $30\,\%$, as may be deduced from the two neutral curves in figures 7 and 8 for $Mo=2.54\times 10^{-11}$. The wake alone is not the source of the instability, since it is still stable at the FSA threshold (Cano-Lozano et al. 2016a). Therefore, as already stated by Bonnefis et al. (2023) in the case of water, the key explanation of path instability in low-Morton-number fluids lies in the specific physical processes accounted for in FSA, namely the degrees of freedom allowed by the rigid-body buoyancy-driven motion of the bubble, and the back reaction imposed by the fluid to this motion by the constant-force and zero-torque constraints.

Increasing the bubble aspect ratio to the range $2.23\leqslant \chi \leqslant 2.25$, Tchoufag et al. (2014b) observed that, although it is still present, the low-frequency mode is no longer unstable, its growth rate keeping weak negative values. In the approximation they used, path instability then arises through a stationary bifurcation, similar to what the FSA neutral curves (pale grey lines) in figures 7 and 8 indicate for $7\times 10^{-11}\lesssim Mo\lesssim 3\times 10^{-7}$. The corresponding stationary mode is that associated with the wake instability observed when the bubble is held fixed in a uniform stream. Hence, one has to conclude that the additional degrees of freedom offered by the free motion of a frozen-shape bubble are not sufficient to change qualitatively the overall dynamics of the system in that range. The yellow neutral curve $St=f(Mo)$ in figure 8 proves that, despite their weak relative magnitude (from $0.4\,\%$ for DMS-T00 to $2\,\%$ for $Mo=3\times 10^{-7}$ according to figure 16b), deformations are able to restore the instability of the oscillatory mode and to make it more unstable than the stationary mode, keeping the near-threshold path instability phenomenology unchanged with respect to the one prevailing in the low-$Mo$ range. Nevertheless, expanding the growth rate in the form $\lambda _r\approx \alpha (Mo)(Bo-Bo_c(Mo))$, with $Bo_c$ the threshold Bond number, and considering variations of $\alpha (Mo)$ from water (figure 13a) to DMS-T02 (figure 12a) via DMS-T00 (not shown), one finds $\alpha (Mo)\sim Mo^{-0.15}$. Hence, the larger $Mo$ the slower the increase of the growth rate of the oscillatory mode with the distance to the threshold, despite the fact that the magnitude of surface deformations at the threshold increases as $Mo^{0.23}$ (figure 16b). This confirms that the damping mechanism predicted by FSA, be it with the schematic spheroidal shape prescribed by Tchoufag et al. (2014b) or with the more realistic shape corresponding to the actual base state, is still active in this intermediate $Mo$-range in the case of deformable bubbles. The reasons for this damping are not obvious. For a non-deforming body, the eigenvalues of the system depend directly on the drag and torque associated with a small edgewise translation or a slow rotation about the major axis perpendicular to the plane of motion (Fabre et al. 2011). We believe that, when the bubble geometrical anisotropy is varied, small changes in the structure of the base flow, especially in the near wake, make these loads vary in a non-monotonic manner with $\chi$, hence with $Mo$, which makes eventually the real part of the eigenvalues associated with the oscillatory mode be positive in some ranges of the parameter space but slightly negative in others.

Last, for $\chi \geqslant 2.25$, Tchoufag et al. (2014b) found the first unstable mode to be oscillatory again. Similar to what is noticed on the FSA neutral curve in figure 8, the corresponding oscillations are much faster than those observed for $2.21\leqslant \chi \leqslant 2.23$, the two oscillatory modes predicted by FSA in the two subdomains being distinct (pale and dark grey lines in figures 7 and 8). Deformation effects totally smooth out this abrupt change by allowing the frequency to increase gradually throughout the intermediate $Mo$-range discussed above. Increasing the Reynolds number at a given $\chi$, i.e. the bubble size in a given fluid, Tchoufag et al. (2014b) found that the next mode that becomes unstable is the stationary mode associated with wake instability. The FSA prediction obtained with the actual base flow confirms this conclusion. For instance, the dark grey line in figure 7 shows that, for DMS-T05, the ‘high-frequency’ mode becomes unstable at $Bo\approx 5.0$. Then the stationary mode (pale grey line) becomes unstable in turn at $Bo=5.85$, i.e. for a bubble only $8\,\%$ larger than at the primary threshold. This sequence subsists when the bubble deforms freely, but the gap between the two thresholds increases considerably, with the oscillatory and stationary modes becoming unstable at $Bo=3.4$ and $Bo=5.95$, respectively. Moreover, the influence of deformation manifests itself through the emergence of the $(2,0)$ axisymmetric shape mode that becomes unstable in between these two thresholds. The above comparisons establish that the role of bubble deformations on the most unstable non-axisymmetric mode responsible for path instability deeply differs among the three $Mo$-ranges identified by FSA. This role is limited to mere, albeit significant, quantitative changes in the critical bubble size and oscillations frequency in the low- and ‘high’-$Mo$ regimes, whereas it induces a complete change in the nature of the most unstable mode in the intermediate-$Mo$ regime.

A weakly nonlinear approach inspired by those developed by Fabre, Tchoufag & Magnaudet (2012) and Tchoufag, Fabre & Magnaudet (2015) for buoyancy/gravity-driven rigid bodies is desirable to extend present predictions to slightly supercritical bubble sizes. By incorporating nonlinear corrections to the base state resulting from the primary non-axisymmetric unstable mode, this extension may in particular allow the disappearance of the stationary mode at any Morton number to be rationalized. Such a weakly nonlinear approach would also settle the question of the supercritical or subcritical nature of the primary bifurcation in each Morton-number range, which the present linear approach cannot answer. Another more modest but much cheaper approach worth exploring is the extension to bubbles (and more generally to rigid bodies with a moderate or low body-to-fluid density ratio) of the ‘aerodynamic’ reduced-order model derived by Fabre et al. (2011) for heavy freely falling two-dimensional bodies. This model should be able to predict accurately the most unstable or least stable eigenvalues of the system for bubbles with a frozen shape, provided that the drag and torque coefficients associated with a small edgewise translation and a slow rotation of the bubble are computed in the relevant base flow. This is the route to be followed to quantify the influence of base flow variations resulting from slight changes in the interface geometry on the hydrodynamic reaction experienced by bubbles close to the threshold of path instability.