# Releasing long bubbles trapped in thin capillaries via tube centrifugation and inclination

Alice Marcotte, Pier Giuseppe Ledda, Valentin Buriasco, Paul Dené, François Gallaire, Ludovic Keiser

► To cite this version:

Alice Marcotte, Pier Giuseppe Ledda, Valentin Buriasco, Paul Dené, François Gallaire, et al. Releasing long bubbles trapped in thin capillaries via tube centrifugation and inclination. Journal of Fluid Mechanics, 2024, 999, pp.A9. 10.1017/jfm.2024.746. hal-05226690

HAL Id: hal-05226690

https://hal.science/hal-05226690v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Releasing long bubbles trapped in thin capillaries via tube centrifugation and inclination

Alice Marcotte<fnref n="1" />,†, Pier Giuseppe Ledda<fnref n="2" />, Valentin Buriasco<fnref n="3" />, Paul Dené<fnref n="3" />, François Gallaire<fnref n="3" /> and Ludovic Keiser<fnref n="4" />

(Received 26 April 2024; revised 15 July 2024; accepted 18 July 2024)

In confined systems, the entrapment of a gas volume with an equivalent spherical diameter greater than the dimension of the channel can form extended bubbles that obstruct fluid circuits and compromise performance. Notably, in sealed vertical tubes, buoyant long bubbles cannot rise if the inner tube radius is below a critical value near the capillary length. This critical threshold for steady ascent is determined by geometric constraints related to the matching of the upper cap shape with the lubricating film surrounding the elongated part of the bubble. Developing strategies to overcome this threshold and release stuck bubbles is essential for applications involving narrow liquid channels. Effective strategies involve modifying the matching conditions with an external force field to facilitate bubble ascent. However, it is unclear how changes in acceleration conditions affect the motion onset of buoyancy-driven long bubbles. This study investigates the mobility of elongated bubbles in sealed tubes with an inner radius near the critical value inhibiting bubble motion in a vertical setting. Two strategies are explored to tune bubble motion, leveraging variations in axial and transversal accelerations: tube rotation around its axis and tube inclination relative to gravity. By revising the geometrical constraints of the simple vertical setting, the study predicts new thresholds based on rotational speed and tilt angle, respectively, providing forecasts for the bubble rising velocity under modified apparent gravity. Experimental measurements of motion threshold and rising velocity compare well with theoretical developments, thus suggesting practical approaches to control and tune bubble motion in confined environments.

**Key words:** thin films, bubble dynamics, core-annular flow

†Email address for correspondence: alice.marcotte@sorbonne-universite.fr

<footnote n="1">¹Institut Jean le Rond d'Alembert, Sorbonne Université, Paris, France</footnote>

<footnote n="2">²Dipartimento di Ingegneria Civile, Ambientale e Architettura, Università degli Studi di Cagliari, Cagliari, Italy</footnote>

<footnote n="3">³Laboratory of Fluid Mechanics and Instabilities, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland</footnote>

<footnote n="4">⁴Institut de Physique de Nice, Université Côte d'Azur, Nice, France</footnote>
# 1. Introduction

Air entrapment into liquid-filled channels is encountered in a broad range of applications, from simple hydraulic systems for intravenous filling (Groell, Schaffler & Rienmueller 1997) to CO₂ sequestration in depleted geological oil reservoirs (Corapcioglu, Cihan & Drazenovic 2004; Oldenburg & Lewicki 2006; Wang & Clarens 2012), embolism in circulatory biological systems (Brodribb, Bienaimé & Marmottant 2016; Li et al. 2021) and multiphase microfluidics flows (Baroud, Galllaire & Dangla 2010). In miniaturized fluid systems, air bubbles can be exploited for transport of particles or for mixing processes (Stone, Stroock & Ajdari 2004; Baroud, de Saint Vincent & Delville 2007; Baroud et al. 2010). Conversely, long gas bubbles may represent a challenging issue, since they can occlude the entire cross-section of a channel and reduce the performance of a fluid circuit (Jensen, Goranović & Bruus 2004; van Steijn, Kreutzer & Kleijn 2008; Brodribb et al. 2016).

In application fields involving narrow liquid channels, long bubbles may be challenging to eliminate, while disrupting fluid flow, causing pressure fluctuations and affecting mixing processes. In perfusion systems for cell cultures, these bubbles can have several detrimental effects on cell health and experimental outcomes, such as localized nutrient deprivation, altered pH levels and accumulation of waste products, all of which can negatively impact cell viability and function (Sung & Schuler 2009). As an additional example, in fuel cells, the oxidation of methanol leads to the formation of CO₂ bubbles that reduce a cell's efficiency (Litterst et al. 2006). Thus, a considerable effort has been dedicated to the removal of bubbles in these circuits (see e.g. Sung & Schuler 2009; Cheng & Lu 2014; Guo, Liu & Ran 2022).

Conversely, transport of long bubbles in microfluidic channels can be cleverly exploited, for instance for particle sieving. Since the speed of a bubble is intrinsically linked to the thickness of its surrounding lubricating film (Fairbrother & Stubbs 1935), tuning the velocity of the bubble may be used to separate particles based on their size (Yu, Khodaparast & Stone 2018): monitoring the speed of the bubble may prevent particles larger than the film thickness from reaching the fluid region past the bubble. The bubble thus acts as an active filter that has the high advantage of preventing clogging. Thus, enabling and controlling the motion of elongated bubbles in capillaries can enhance the efficiency of these microfluidic systems.

Many hydraulic and microfluidic systems rely on vertical settings (see e.g. Kaigala *et al.* 2011), thereby calling for a better understanding of how bubble transport is influenced by gravity forces. In a vertical configuration, a gas volume in a liquid-filled channel is expected to rise owing to buoyancy. The more specific case of buoyant ascent in a vertical tube of a gas volume with an equivalent spherical diameter larger than the tube inner radius has been investigated by Dumitrescu (1943) and Davies & Taylor (1950), which provided a prediction for the rising velocity of long bubbles in tubes, subsequently termed Taylor bubbles. However, it was observed over a century ago (Gibson 1913) that long bubbles within a sealed vertical tube with a sufficiently narrow diameter exhibit an interesting behaviour: they cease to rise and appear to be stuck. This observation was puzzling considering the existence of a thin lubricating film surrounding long bubbles in tubes, allowing for the drainage of the fluid and thus for the rising of the bubble. Bretherton (1961) showed that if the tube inner radius was smaller than a critical value $R_c$ close to the capillary length $\ell_c$ of the liquid (more precisely, $R_c \approx 0.918\ell_c$), no valid bubble shape was compatible with a steady rising motion. This threshold stems from the asymptotic matching between the upper cap profile, which results from the equilibrium between surface tension and gravity, and the thin film surrounding the elongated part
Releasing long bubbles trapped in thin capillaries

<figure><img src="image_3.png" /><figcaption>Figure 1. (a) Schematics of a long bubble immersed in a viscous liquid inside a sealed capillary. The top part of the bubble can be divided into an upper cap and an elongated part surrounded by a thin film. For a buoyant bubble to rise, mass conservation requires the fluid displaced by the tip of the bubble to drain through the thin film. (b) Sketch of the upper cap profile of a long air bubble within a sealed tube of radius $R$, in the vertical setting studied by Bretherton (1961). The profile exhibits an inflection point denoted by $I$. For $R > R_c$, the matching with the thin-film region at the inflection point is possible. (c,d) Sketch of the configurations investigated in this study. In the first case (c), the tube is held vertically and rotates around its symmetry axis at angular frequency $\omega$. In the second case (d), the tube is tilted with respect to gravity and makes an angle $\alpha$ with the horizontal plane.</figcaption></figure>

of the bubble, where viscous, surface tension and gravity forces are balanced. These two regions are depicted in figure 1(a). The condition for the onset of motion can be summarized as a geometrical constraint which imposes, for the existence of a steadily ascending bubble, that the upper cap profile exhibits an inflection point with negative slope (for an upward-oriented vertical axis); see figure 1(b). At the critical condition $R = R_c$, both the slope and curvature vanish at the solid wall. In addition, for $R$ slightly larger than $R_c$, Bretherton (1961) predicted the bubble rising velocity, by exploiting mass conservation through the thin film and the variation of the slope at the inflection point with the tube's radius.

Below the threshold $R < R_c$, Lamstaes & Eggers (2017) studied the unsteady bubble motion and predicted the occurrence of a self-similar pinch-off singularity of the thin lubricating film around the bubble, thus hindering any further flow and eventually stopping the progression of the bubble, found to travel a finite distance over infinite time. That prediction is supported by recent interference microscopy experiments that have demonstrated that the bubble is apparently stuck by an infinitely slow flow taking place in the surrounding thin liquid film whose nanometric thickness results from an equilibrium between capillary stress and disjoining pressure (Dhaouadi & Kolinski 2019).

For bubble ascent in sealed tubes, it is necessary for the fluid displaced by the tip of the bubble to drain through the thin film. Thus, enabling the motion of the bubble in sealed tubes with inner radii smaller than the critical value $R_c$ requires the development of some strategies that would act on the thickness of the surrounding lubricating film. Zhou & Prosperetti (2021) showed numerically that ‘encaging’ the bubble by means of thin vertical rods regularly arranged on a circle coaxial with the tube could effectively expand the gap between the air-liquid interface and the inner solid wall, thus facilitating the downward flow of the liquid and increasing the rising velocity of the bubble. Bi & Zhao (2001) and Bico & Quéré (2002) demonstrated that using angular tubes could effectively promote the rising of the bubble even under strong confinement, owing to the presence of corners that allow for a more efficient drainage of the liquid around the bubble (Funada et al. 2005). In the same spirit, another strategy consists of using textured inner walls: because of the imbibition of the roughness, the effective thickness of the lubricating film is actually larger than on a smooth surface (Bico, Tordeux & Quéré 2001).
A. *Marcotte and others*

However, in some applications where the geometry of the tube cannot be modified adequately, the film thickness could be varied by adjusting the pressure distribution in the surrounding liquid by means of an external force field, which could be easily tuned so as to precisely control the ascent velocity of the bubble. In this context, it has been shown that imposing a liquid flow in the tube effectively thickens the lubricating film around the bubble (Yu *et al.* 2021). In particular, Magnini *et al.* (2019) demonstrated that when the external flow is oriented in the same (upward) direction as buoyancy, it can enable the rise of bubbles in tubes with radius $R < R_c$. Kubie (2000) documented a significant increase in the ascent velocity of a Taylor bubble enclosed in a vertical tube subjected to horizontal oscillations. In the case of a vertically oscillated tube, Brannock & Kubie (1996) reported instead experimental evidences of the slowing down of the bubble, while Madani, Caballina & Souhar (2009, 2012) observed a more nuanced behaviour: as the acceleration of the oscillations is gradually increased, the rising of the bubble initially slows down, but then increases at larger accelerations. More recently, Zhou & Prosperetti (2024) studied the rising behaviour of a Taylor bubble exhibiting volume oscillations imposed either by forcing the liquid column above the bubble to oscillate or by imposing a pulsating pressure field at the top liquid surface. Their numerical simulations evidence that the gas volume oscillations result in the thinning of the lubricating film around the bubble and thus in the decrease of the drainage flow and rising velocity.

Here, we focus on the transport of long bubbles in sealed tubes filled with a viscous liquid, with an inner radius close to the critical value below which the bubble is stopped in a vertical configuration. We investigate two different strategies to enable bubble motion and tune its velocity, namely rotating the tube around its symmetry axis and inclining it with respect to gravity (figure 1c,d). In both cases, we leverage theoretical developments to predict the new threshold for the onset of motion that depends on the rotational speed and on the tilt angle, respectively. We also provide a prediction for the rising velocity of the bubble as a function of the liquid properties, the tube geometry and the (modified) gravity field. Our theoretical findings are then compared with the outcomes of dedicated experimental campaigns.

The paper is organized in two parts. In the first part (§2), we report our investigation on bubble motion in rotating tubes. Section 2.1 develops the theoretical prediction for the cap profile of the bubble and the matching conditions between this cap and the flat-film region, from which we derive the theoretical threshold for the onset of motion and the prediction of the bubble velocity, in terms of the rotational speed. Section 2.2 presents the experimental set-up and a comparison of the results against the theoretical findings. The second part (§3) presents the same structure as the previous part, but investigates the effect of tube inclination, with theoretical predictions and comparison with experimental measurements of bubble transport in tilted tubes.

## 2. Effect of centrifugation

Fluid centrifugation pertains to extensive applications, ranging from the segregation of complex or biological fluids (Svedberg & Fåhraeus 1926) to numerous industrial processes, such as wastewater treatment (Turano *et al.* 2002) or crude oil refining (Gary *et al.* 2007). In interfacial flows, spinning rods (Than *et al.* 1988) and spin-coating (Emslie, Bonner & Peck 1958) are used to deposit uniform thin films onto diverse substrates such as optical lenses for anti-reflective properties (Krogman, Druffel & Sunkara 2005), or silicon wafers for organic semiconductor fabrication (Yuan *et al.* 2014). This method precisely controls the film thickness through the modulation of the angular velocity, essential for achieving high-quality coatings. Additionally, centrifugation can be employed in the
generation of surface roughness in curing polymer melts (Marthelot, Strong & Brun 2018; Jambon-Puillet, Royer Piéchaud & Brun 2021), where centrifugal instabilities (Rietz et al. 2017) are harnessed to facilitate the formation of periodic patterns.

In sealed tubes, the effect of centrifugation on the shape of capillary interfaces has been exploited in spinning drop experiments (Vonnegut 1942; Rosenthal 1962; Princen, Zia & Mason 1967; Torza 1975). These experiments can measure very low interfacial tensions (Drelich, Fang & White 2002) by rotating a horizontal tube containing a drop of lower-density liquid within a higher-density fluid. For high enough rotation rate, the (transverse) gravity acceleration can be neglected, and the equilibrium shape of the drop results from the balance of the centrifugal force, which tends to elongate the drop along its axis (and thus to thicken the surrounding liquid film), with surface tension, which promotes a spherical shape.

In this context, Manning, Collicott & Finn (2011) studied the case of a tube of inner radius R partially filled with a liquid of density ρ and surface tension γ, rotated around its symmetry axis at angular velocity ω under weightlessness. They derived a criterion for the occlusion of the tube by a static meniscus spanning the cross-section of the channel, with a contact angle φ, and computed a critical angular velocity ω₀:

$$ \frac{\rho\omega_0^2 R^3}{\gamma} = 32 \sin^3 \left( \frac{\pi + 2\phi}{6} \right), \qquad (2.1) $$

such that the tube cannot occlude if ω > ω₀. In the case of an occluding meniscus forming a bubble, the gas-liquid interface meets the solid wall tangentially, i.e. φ = 0. Thus, (2.1) indicates that under weightlessness, the bubble cannot occlude the channel if ρω²R³/γ > 4.

To the best of our knowledge though, the combined effect of axial gravity and transverse centrifugal force on the rising motion of long bubbles in a vertical setting has not been studied yet. Building on the demonstrated ability of centrifugation to elongate light drops or bubbles in tubes, and thereby to thicken their lubricating film, we now study how centrifugation can facilitate the release of long bubbles that are trapped in sealed capillaries due to surface tension.

We consider a long bubble of length L immersed in a viscous fluid of dynamic viscosity μ, density ρ and surface tension γ, both contained in a vertically oriented circular tube of radius R ≪ L, sealed at both ends. The bubble ascends along the vertical axis at a constant velocity U_b under the influence of gravity. The tube's radius is assumed to be of the order of the capillary length ℓ_c = √γ/ρg, where g is the acceleration due to gravity, so that the Reynolds number Re = ρU_bR/μ is sufficiently small to neglect any inertial effects.

Bretherton's solution describing the bubble's ascent at a constant velocity is valid only if the tube radius exceeds a critical value R_c ≈ 0.918ℓ_c. As the tube radius R approaches this critical value, the bubble's ascending speed diminishes, eventually reaching zero. This phenomenon can be explained through a simple mass conservation consideration: the sealed tube requires the rising bubble to displace the liquid above, creating drainage through its peripheral lubricating film. However, for R < R_c = 0.918ℓ_c, surface tension becomes dominant, causing the bubble to expand and occupy the entire tube cross-section, preventing liquid drainage.

We now examine the scenario where the vertical tube undergoes constant rotation around its symmetry axis with an angular velocity ω. We can readily anticipate that the centrifugal force will push liquid towards the solid tube wall, thickening the fluid film around the bubble and facilitating its ascent. Therefore, a steady rising motion of the bubble may be achievable even in tubes with R < R_c, provided the angular velocity ω
is sufficiently high. In the subsequent section, we revisit Bretherton's theory (Bretherton 1961) to predict the new threshold $R_c(\omega)$ and the steady rising velocity $U_b$ as a function of the rotational speed.

## 2.1. Theoretical prediction for the threshold and rising velocity

Since our focus lies in describing motion near the threshold characterized by a vanishing velocity $U_b$, we preliminarily assume a small capillary number $Ca = \mu U_b / \gamma \ll 1$. Viscous stresses at the gas-liquid interface thus play a significant role only in regions where the fluid is strongly confined, i.e. where the interface is very close to the solid wall. Consequently, the upper part of the bubble's profile can be divided into two regions; see figure 1(a). The outer region corresponds to the top of the bubble (cap) where viscous effects are negligible: the equilibrium is controlled by an interplay among surface tension, gravity and centrifugal forces. Conversely, a thin liquid film resulting from a balance among viscous forces, surface tension, gravity and centrifugal forces defines an inner region of small axial curvature. We now derive the bubble's profiles in these two regions.

### 2.1.1. The static cap

In the outer region, the fluid around the cap of the bubble can be considered at rest (Bretherton 1961; Lamstaes & Eggers 2017) so that in the cylindrical reference frame co-rotating with the tube and translating with the bubble, the pressure $P$ in the surrounding fluid satisfies

$$0 = -\nabla P + \rho\omega^2 re_r - \rho g e_z. \qquad (2.2)$$

By integrating the radial component of (2.2), we obtain

$$\left.
\begin{aligned}
P(r, z) &= \frac{1}{2} \rho \omega^2 (r^2 - r_1(z)^2) + \gamma \kappa + P_{\text{air}}, \\
\kappa &= -\frac{1}{r_1(z)(1 + r'_1(z)^2)^{1/2}} + \frac{r''_1(z)}{(1 + r'_1(z)^2)^{3/2}},
\end{aligned}
\right\}
\qquad (2.3)$$

where $r_1(z)$ and $\kappa$ denote the location of the air-liquid interface measured from the central axis (oriented upwards) and its curvature, respectively.

From the axial component of the momentum conservation equation (2.2), it follows that

$$\gamma\kappa - \frac{1}{2}\rho\omega^2 r_1(z)^2 + \rho gz = \text{cst.} \qquad (2.4)$$

By denoting as $s$ the arc-length of the interface profile measured from the tip of the bubble and $\theta$ its tangent angle with respect to the horizontal (see figure 2), the static interface profile is given by

$$-\gamma \left[ \frac{d\theta}{ds} + \frac{\sin\theta}{r_1(s)} \right] - \frac{1}{2} \rho \omega^2 r_1(s)^2 + \rho gz(s) = \text{cst}, \qquad (2.5)$$

where the coordinates $(r_1(s), z(s))$ locate the position of the gas-liquid interface at $s$ in the $(r, z)$ plane. Using that $dr_1/ds = \cos(\theta)$ and $dz/ds = -\sin(\theta)$, differentiating with respect to the curvilinear coordinate $s$ gives

$$\gamma \left[ \frac{d^2\theta}{ds^2} + \frac{\cos\theta}{r_1(s)} \frac{d\theta}{ds} - \frac{\cos\theta \sin\theta}{r_1(s)^2} \right] = -\rho\omega^2 r_1(s) \cos\theta - \rho g \sin\theta. \qquad (2.6)$$
<figure><img src="image_4.png" /><figcaption>Figure 2. Sketch of a bubble in a vertical tube that rotates around its central axis with angular velocity $\omega$. Here $s$ is the arc-length of the interface measured from the tip of the bubble at $r = 0$. In the static cap region, the air-liquid interface is located at a distance $r_1(s)$ from the central axis, and the angle its tangent makes with the horizontal axis is denoted as $\theta(s)$. In the inner region, where a two-dimensional Cartesian system $(x, y)$ is used, the interface is located instead by its distance from the solid wall $y_1(x)$.</figcaption></figure>

Finally, with the dimensionless variables $\bar{r}_1 = r_1/R$ and $\bar{s} = s/R$, (2.6) becomes

$$ \frac{d^2\theta}{ds^2} + \frac{\cos(\theta)}{\bar{r}_1} \frac{d\theta}{ds} - \frac{\cos(\theta)\sin(\theta)}{\bar{r}_1^2} = -Bo \sin(\theta) - \bar{r}_1 Ce \cos(\theta), \quad (2.7) $$

where the Bond number $Bo = \rho g R^2 / \gamma = (R / \ell_c)^2$ is introduced as the square of the ratio between the tube radius and the capillary length. The centrifugal number $Ce = \rho \omega^2 R^3 / \gamma$ can be seen as a rotational Bond number where the centrifugal acceleration $R\omega^2$ plays the role of the gravitational acceleration.

For a given set of parameters $(Bo, Ce)$, two boundary conditions are required to integrate (2.7) from the position $(r_1(s=0) = 0, z(s=0) = 0)$. A first condition is provided by the symmetry of the problem, which imposes $\theta(0) = 0$ at the top of the static cap. The second boundary condition will be determined upon matching of this static profile with that of the inner region.

### 2.1.2. The thin-film region

In the inner region where the bubble is surrounded by a thin lubricating film, the film's thickness is extremely small compared with the tube's radius. Following Bretherton (1961), we thus neglect the azimuthal curvature of the air-liquid interface and consider the thin-film region as planar instead of annular.

Under these assumptions, we introduce the two-dimensional, stationary, Cartesian coordinate system $(x, y)$, where $x = z - U_b t$ opposes gravity, and $y = R - r$ represents the distance to the solid wall. In the framework of the lubrication approximation, the viscous flow in the thin film is driven by a pressure gradient resulting from a combination of gravity, capillarity and centrifugal force. The axial velocity accordingly is written as (see Appendix A.1 for a detailed derivation)

$$ u(x, y) = \frac{\gamma}{2\mu} \left( -y_1''' + \frac{\rho\omega^2 R}{\gamma} y_1' + \frac{\rho g}{\gamma} \right) (y^2 - 2y_1 y) - U_b, \quad (2.8) $$
where $y_1(x)$ denotes the distance of the air-liquid interface to the solid wall of the tube and a prime denotes the derivative with respect to $x$. In (2.8), the first term of the right-hand side stems from surface tension effects, the second from the centrifugal force and the third from gravity. Upon integration within the thin film, the volume flux reads

$$Q = -2\pi R U_b y_1 - 2\pi R \frac{\gamma}{3\mu} \left( -y_1''' + \frac{\rho\omega^2 R}{\gamma} y_1' + \frac{\rho g}{\gamma} \right) y_1^3. \quad (2.9)$$

This flux must equate the volume of fluid displaced per unit time by the top of the bubble, which is equal to $\pi R^2 U_b$. Since $y_1/R \ll 1$, the $-2\pi R U_b y_1$ term in the expression of the flow rate is a negligible correction. Finally, by imposing flux continuity with the region far away from the tip, where the film thickness can be considered as uniform and equal to a constant $b$, we obtain the following thin-film equation:

$$y_1''' = \frac{\rho g}{\gamma} \left( 1 - \frac{b^3}{y_1^3} \right) + \frac{\rho \omega^2 R}{\gamma} y_1'. \quad (2.10)$$

Since $b$ is the length scale governing the flow in the inner region, we non-dimensionalize, as in Bretherton (1961), with

$$y_1 = \eta b, \quad x = \zeta b (\rho g b^2 / \gamma)^{-1/3}. \quad (2.11a,b)$$

This leads to the following ordinary differential equation:

$$\eta''' = \frac{\eta^3 - 1}{\eta^3} + a\eta', \quad (2.12)$$

where $a = (Ce/Bo^{2/3})(b/R)^{2/3}$. In (2.12), the left-hand side represents the surface tension term, while on the right-hand side, the first term accounts for gravity, the second for viscous dissipation and the third for the effect of centrifugation.

### 2.1.3. Matching

We aim at matching the inner solution, which is described by (2.12), with the static cap solution provided by (2.6). Given that $b \ll R$, this requires taking the limit $\eta \to \infty$ in (2.12) for the inner solution. In this limit, the equation behaves as

$$\eta''' = 1 + a\eta', \quad (2.13)$$

whose general solution reads

$$\eta = c_1 e^{\sqrt{a\zeta}} + c_2 e^{-\sqrt{a\zeta}} + c_3 - \frac{\zeta}{a}. \quad (2.14)$$

The values of $c_1$, $c_2$ and $c_3$ can be found through interpolation using the numerical solution of the complete equation (2.12), whose initial conditions are obtained from the uniform
film solution. Indeed, when $\eta \to 1$, (2.12) becomes

$$ \eta''' = 3(\eta - 1) + a\eta', \qquad (2.15) $$

for which the only non-oscillating solution is $\eta_0 = 1 + C \exp(\mathcal{F}\zeta)$, where $C$ is an integration constant and

$$ \mathcal{F} = \frac{\left(\frac{2}{3}\right)^{1/3} a}{(27 + \sqrt{3}\sqrt{243 - 4a^3})^{1/3}} + \frac{(27 + \sqrt{3}\sqrt{243 - 4a^3})^{1/3}}{2^{1/3} 3^{2/3}}. \qquad (2.16) $$

Since the value of $C$ can be adjusted by shifting the origin of $\zeta$, we can set $C = 1$ and a large, negative initial value $\zeta_0$ to initialize the integration. This procedure yields initial conditions for the full nonlinear equation (2.12), i.e.

$$ \eta(\zeta_0) = 1 + \exp(\mathcal{F}\zeta_0), \quad \eta'(\zeta_0) = \mathcal{F}\exp(\mathcal{F}\zeta_0), \quad \eta''(\zeta_0) = \mathcal{F}^2 \exp(\mathcal{F}\zeta_0). \quad (2.17a-c) $$

We solve (2.12) (using the built-in MATLAB ODE solver ode45) to obtain the inner solution $\eta$ for various values of $a$, and fit the outer profile $\eta = c_1 e^{\sqrt{a}\zeta} + c_2 e^{-\sqrt{a}\zeta} + c_3 - \zeta/a$ in the region where $\eta \gg 1$. This allows us to retrieve the coefficients $c_1(a), c_2(a)$ and $c_3(a)$. From figure 3(a), it is evident that the outer profile $\eta$ of the inner solution exhibits an inflection point. We thus translate the origin to the position where $\eta'' = 0$, located at the coordinate

$$ \zeta^* = \frac{1}{2\sqrt{a}} \log \left( -\frac{c_2}{c_1} \right). \qquad (2.18) $$

Using the shifted variable $\chi = \zeta - \zeta^*$, we can now define $\eta(\chi) = c_1^* e^{\sqrt{a}\chi} + c_2^* e^{-\sqrt{a}\chi} + c_3^* - \chi/a$, where the new coefficients $c_1^*, c_2^*$ and $c_3^*$ are expressed as

$$ c_1^* = c_1 e^{\sqrt{a}\zeta^*} = c_1 \sqrt{-\frac{c_2}{c_1}} = \operatorname{sgn}(c_1)\sqrt{-c_1c_2}, \qquad (2.19a) $$

$$ c_2^* = c_2 e^{-\sqrt{a}\zeta^*} = c_2 \sqrt{-\frac{c_1}{c_2}} = \operatorname{sgn}(c_2)\sqrt{-c_1c_2}, \qquad (2.19b) $$

$$ c_3^* = c_3 - \frac{\zeta^*}{a} = c_3 - \frac{1}{2a^{3/2}} \log \left( -\frac{c_2}{c_1} \right) = \eta(\chi = 0). \qquad (2.19c) $$

The new coefficients are well fitted (see figure 3b) by

$$ c_1^* = -c_2^* \approx 0.500a^{-3/2} + 0.286a^{-1/2}, \qquad (2.20a) $$

$$ c_3^* \approx 1.10, \text{ independently of the value of } a. \qquad (2.20b) $$

Thus, the distance from the wall at which the outer profile exhibits an inflection point can be evaluated as

$$ y_1(0) = \eta(0)b = (c_1^* + c_2^* + c_3^*)b \approx 1.10b \ll R. \qquad (2.21) $$

Furthermore, the slope of the profile at the inflection point is given by

$$ \begin{aligned} y_1'(0) &= \eta'(0)(\rho gb^2/\gamma)^{1/3} = \left(c_1^* \sqrt{a} - c_2^* \sqrt{a} - \frac{1}{a}\right) (\rho gb^2/\gamma)^{1/3} \\ &= 0.572(\rho gb^2/\gamma)^{1/3} > 0. \end{aligned} \qquad (2.22) $$

Remarkably, the conditions on film thickness (2.21) and slope (2.22) at the inflection point are the same as described in Bretherton (1961). Therefore, the centrifugal force alters
<figure><img src="image_5.png" /><figcaption>Figure 3. (a) The inner region profile $\eta$ as a function of the dimensionless height $\zeta$. The circles represent the solution $\eta$ of the full equation (2.12) for various values of $a$, while the black solid lines represent the outer profiles of the inner region $\eta = c_1 e^{\sqrt{a}\zeta} + c_2 e^{-\sqrt{a}\zeta} + c_3 - \zeta/a$, where $c_1(a)$, $c_2(a)$ and $c_3(a)$ are obtained by fitting with the full inner solution, in the $\eta \gg 1$ region. The outer profiles clearly exhibit an inflection point, at a distance referred to as $\zeta^*(a)$. For each value of $a$, the origin is then shifted so that $\eta''(0) = 0$. (b) Shifted coefficient $c_1^* = -c_2^*$ as a function of $a$ (red circles). The black solid line corresponds to $c_1^* = 0.500a^{-3/2} + 0.286a^{-1/2}$. Inset: shifted coefficient $c_3^* = \eta(0) \equiv \eta^*$ as a function of $a$. This coefficient does not vary significantly with $a$.</figcaption></figure>

the inner region solution and the static cap profile, but the matching conditions (and thus the boundary conditions for the static cap solution) remain surprisingly unchanged from those of Bretherton (1961).

The above analysis provides the missing information required to solve the static cap profile described by (2.7). Specifically, the static profile for a given set of parameters (Bo, Ce) is obtained through a shooting method, searching for the first derivative $\dot{\theta}_0$ at the tip of the cap, which is such that the integration of (2.7) from initial conditions ($\theta_0 = 0, \dot{\theta}_0$) results in a profile where the inflection point $\ddot{r}_1(z) = 0$ is reached for $r_1(z) = R$. The numerical integration of (2.7) is performed using the MATLAB built-in ODE solver ode23t, with a spacing along the curvilinear coordinate of 0.01R, while the shooting method is implemented by means of the nonlinear MATLAB system solver fsolve. The resulting slope $r'_1(z)|_{r=R}$ at the inflection point is then computed from the generated profile.

For fixed $Bo < Bo_{c,0} = (R_c/\ell_c)^2$ and varying $Ce$, it appears that some profiles are unphysical: below a critical value $Ce_c$ that depends on the Bond number $Bo$, the static cap shape exhibits a positive slope at the inflection point at the solid wall, causing the upper profile to extend beyond the fluid domain $r < R$, as illustrated in figure 4, which is not compatible with the matching condition $0 < y'_1(0) = -r'_1(z)|_{r=R}$. By progressively increasing the centrifugal number $Ce$, the slope at the inflection point at the wall decreases, leading to a reduction of the bulge outside the fluid domain, as shown in figure 4(c). Ultimately, for $Ce > Ce_c$, the slope becomes negative, causing the entire static cap to reside within the fluid domain; see figure 4(b,c).

In the following, we denote as $\phi$ the resulting angle between the liquid-air interface and the vertical axis at the inflection point, i.e. $\phi = \tan^{-1}(-r'_1(z)|_{r=R}) = \tan^{-1}(y'_1(0))$. A closer inspection reveals that within a small range around $Ce_c$, i.e. for $|Ce - Ce_c(Bo)| < 0.2$, $\phi$ varies linearly with $Ce$, as shown in figure 5(a). Within this range, the angle $\phi$ (in radians) is well approximated by

$$ \phi(Ce, Bo) \approx 0.144(Ce - Ce_c(Bo)), \qquad (2.23) $$
Releasing long bubbles trapped in thin capillaries

<figure><img src="image_6.png" /><figcaption>Figure 4. Static cap profile for a Bond number $Bo = 0.55$ and a centrifugal number $Ce$ that is (a) below and (b) above the threshold $Ce_c(Bo)$. Below the threshold, the slope $r'_1(z)|_{r=R}$ is positive at the inflection point, causing the upper profile to escape the fluid domain $r < R$. Angle $\phi$ is the angle between the tangent to the static cap profile at the wall and the vertical axis: $\phi = \tan^{-1}(-r'_1(z)|_{r=R})$ and is thus negative in (a) and positive in (b). (c) Evolution of the static cap profile at fixed $Bo = 0.55$, when increasing $Ce$ from a value below the threshold $Ce_c \approx 1$ to a value slightly above threshold.</figcaption></figure>

where the factor 0.144 is independent of $Bo$ (up to variations less than 0.001 radians). Interestingly, this prediction is consistent with the occlusion criterion derived by Manning et al. (2011) under weightlessness (i.e. $Bo = 0$), reported in the introduction of this section (equation (2.1)). Indeed, $\phi$ represents the slope of the gas-liquid interface at the inflection point, located at a distance 1.10b from the solid wall of the tube (equation (2.21)). In this sense, it can be seen as the contact angle of a meniscus occluding a virtual channel of radius ($R - 1.10b$), rotating at angular velocity $\omega$. Thus, in the absence of gravity, and in the limit of small contact angle $\phi$, the occlusion criterion (2.1) becomes

$$Ce(\phi, Bo = 0) = \frac{\rho\omega^2(R - 1.10b)^3}{\gamma} = 32 \sin^3\left(\frac{\pi + 2\phi}{6}\right) \approx 4(1 + \sqrt{3}\phi), \quad (2.24)$$

which can be recast as $\phi(Ce, Bo = 0) \approx 0.144(Ce - Ce_{c,0})$, where $Ce_{c,0} \equiv Ce_c(Bo = 0) = 4$.

Thus, at fixed $Bo$, the critical centrifugal number $Ce_c(Bo)$ for vanishing angle $\phi$ (or equivalently for vanishing slope) is retrieved as the value of $Ce$ at which the best linear fit of $\phi(Ce, Bo)$ cancels out. Its dependency on the Bond number is depicted in figure 5(b). Note that we performed a convergence analysis and observed no further variations of $Ce_c$, within a tolerance of 0.07%, when increasing 10 times the resolution on the spacing along the curvilinear coordinate used to integrate the static cap profiles. For $Ce < Ce_c$, the geometrical constraint $\phi = \tan^{-1}(y_1(0)) \approx y'_1(0) > 0$ cannot be satisfied, so that this value corresponds to the threshold for the onset of motion.

A second-order polynomial approximation of $(Ce_c(Bo), Bo)$ is $Bo \approx 0.842 - 0.295 Ce_c(Bo) + 0.020 Ce_c^2(Bo)$, whose solution, shown in figure 5(b), reads

$$Ce_c = \frac{0.295 - \sqrt{0.295^2 - 0.080(Bo_{c,0} - Bo)}}{0.040}, \quad (2.25)$$
A. Marcotte and others

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) Angle $\phi$ of the static cap profile between the vertical axis and the tangent to the static cap profile (obtained by integrating (2.7) while requiring that the interface reaches the solid wall with an inflection point) as a function of the centrifugal number $Ce$, for various Bond numbers $Bo$. The circles are the values computed from the shape of the interface, while the solid lines are the best linear fit. For each $Bo$, the critical centrifugal number $Ce_c$ is defined as the value of $Ce$ for which the best linear fit cancels out. (b) The critical centrifugal number $Ce_c$ as a function of the Bond number $Bo$. At a given value of $Bo$, the matching with the inner region is only possible if $Ce > Ce_c(Bo)$. The circles are the values of $Ce$ that cancel the linear approximation of $\phi(Ce, Bo)$ for each Bond number, while the back line represents the approximation (2.25). The red circle with coordinates $(Bo_{c,0} = 0.842, Ce_c = 0)$ locates the threshold in the absence of centrifugation.</figcaption></figure>

where $Bo_{c,0} = 0.842$ is the critical Bond number in the absence of centrifugation ($Ce = 0$). Note that in the limit $Bo \to 0$, the approximation (2.25) yields $Ce_c(Bo = 0) \approx 3.9$, which is close to the threshold $Ce_{c,0} = 4$ computed by Manning et al. (2011) under weightlessness. Conversely, in the limit ($Ce_c \to 0, Bo \to Bo_{c,0}$), (2.25) simplifies into

$$Ce_c \approx \frac{1}{0.295}(Bo_{c,0} - Bo). \qquad (2.26)$$

By injecting (2.26) into (2.23), we obtain, in the limit $Ce \to 0$:

$$\phi \approx 0.49(Bo - (Bo_{c,0} - 0.295Ce)), \qquad (2.27)$$

reminiscent of the expression derived by Bretherton for a non-rotating capillary tube ($\phi = 0.49(Bo - Bo_{c,0})$). The similarity of these expressions highlights the role of the centrifugation as a downward shift in the critical Bond number for the onset of motion.

In addition, this analysis provides a prediction for the rising velocity of the bubble for $Ce > Ce_c(Bo)$. Indeed, for $Ce$ close enough to the threshold $Ce_c(Bo)$, the slope of the inner solution at the inflection point should verify, using (2.21), (2.22) and (2.23):

$$\begin{aligned} y'_{1}(0) &= 0.572(\rho gb^2 / \gamma)^{1/3} \\ &= \tan(\phi) \approx \phi \approx 0.144 \left( \frac{\rho \omega^2 (R - 1.10b)^3}{\gamma} - Ce_c(Bo) \right). \end{aligned} \qquad (2.28)$$

Since the volume of fluid displaced per unit time by the tip of the bubble $\pi R^2 U_b$ should be equal to the volume flux in the uniform film region, we can relate the thickness $b$ to the inner radius $R$ and the velocity $U_b$ through $\rho gb^3/3\mu U_b = R/2$. Expanding (2.28) up to first order in $(b/R)$ yields the following expression of $Ca$ as an implicit function of $Ce$
and Bo:

$$Ce - Ce_c = 3.78Ce \left(\frac{Ca}{Bo}\right)^{1/3} + 4.35Bo^{1/3} \left(\frac{Ca}{Bo}\right)^{2/9}, \quad (2.29)$$

where $Ce_c$ is the function of $Bo$ described above. The ratio $Ca/Bo$ compares the bubble velocity $U_b$ with the settling velocity $U^* = \rho g R^2 / \mu$. Note that (2.29) admits an analytical solution: by setting the unknown $x = (Ca/Bo)^{1/9}$, (2.29) can indeed be recast as an equation of the type $x^3 + a_1 x^2 + a_2 x + a_3 = 0$.

## 2.2. Experiments on centrifugated bubbles

In this section, we outline our experimental set-up and procedure and compare the results against the above theoretical developments.

### 2.2.1. Experimental set-up and procedure

Cylindrical borosilicate capillary tubes (Hilgenberg GmbH) are partially filled with silicone oil (Sigma Aldrich, $\rho = 964 \text{ kg m}^{-3}$, $\mu = 9.64 \times 10^{-2} \text{ Pa s}$, $\gamma = 2.09 \times 10^{-2} \text{ N m}^{-1}$), leaving an air bubble with a length $L$ greater than 10 times the radius $R$ of the tube. Both ends of the tubes are sealed with epoxy resin. The inner radii of the tubes used in our experimental campaign vary between 0.8 and 1.3 mm, corresponding to Bond numbers $Bo \equiv \rho g R^2 / \gamma$ in the range [0.29, 0.76]. Note that the uncertainty in the inner diameters of the tubes is of 0.05 mm.

The tube attachment system, presented in figure 6(a), consists of two circular mounts made of PETG, rigidly connected together via two vertical steel rods, and linked by bearings to a fixed aluminium frame (not represented in the figure). On each mount, a central, threaded circular mouthpiece accommodates a hollow cylinder whose inner radius matches the outer radius of the capillary tube. The extremities of the tube are then inserted into these cylinders, and securely clamped to the mounts using a clamping chunk. By this means, the tube can be easily replaced by a capillary of a different size with minimal adjustments of the set-up. The lower mount is also connected to the shaft of a DC motor that imposes the rotation of the tube attachment system. While the tube, the mounts and the rods rotate collectively, the verticality and stability of the whole set-up are ensured by the fixed aluminium frame.

The motor is voltage-controlled to achieve the desired rotation speed, measured with less than 1% error using a tachometer. To enhance visualization, a light-emitting-diode panel is positioned behind the tube attachment system.

A Basler camera records the evolution of the bubble in the tube, and the velocities of the upper and lower caps of the bubble are obtained through image post-processing performed via a custom MATLAB script. Specifically, a column of pixels aligned with the tube that crosses the upper and bottom profiles of the bubble is extracted from each frame. These slices are then juxtaposed to each other into an image where the horizontal axis represents time. The displacement of the bubble extremities with time are clearly visible on the resulting image, as shown in figure 6(c).

Once the motor is switched on, a transient regime occurs where the upper cap of the bubble rises while the bottom cap remains immobile, resulting in bubble elongation, reminiscent of spinning bubble experiments (Vonnegut 1942). This is accompanied by the progressive thickening of the surrounding film that propagates from the top to the lower cap of the bubble, as seen in figure 6(b). Once the propagation front reaches the bottom extremity, the lower cap starts its ascent at the same (constant) velocity as the upper cap;
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

Figure 6. Experimental set-up and post-processing for rotating bubbles. (a) Tube attachment system. The capillary tube is clamped on both extremities to mounts connected together by two metal rods. The bottom mount is linked to the shaft of a DC voltage-controlled motor that imposes the rotation of the system around its central, vertical axis. (b) Photographs of a long bubble inside a tube filled with silicone oil at different and equally spaced time steps within the transient regime. In the red frame, the motor has been switched on and the upper cap starts rising while the bottom cap remains still. Along with the resulting bubble elongation, the surrounding liquid film gets progressively thicker from the top to the bottom part of the bubble. The dotted line roughly locates the position of the propagation front. Once the front has reached the lower cap, it starts rising. (c) Intensity profile as a function of time along the tube axis. To produce this image, a column of pixels aligned with the central axis of the tube is extracted from each frame of the movie. The columns are then juxtaposed to each other. The locations of the upper and lower cap as a function of time are easily identified as the two roughly parallel black curves limiting a darker domain that corresponds to the position of the bubble itself. At time $t_1$, the motor is switched on. At $t_2$, the bubble dynamics reaches a stationary state: the upper and lower caps rise at same constant velocity, as highlighted by the parallel blue solid lines that are superimposed on the position of the caps as a function of time. For (b,c), $R = 1.2$ mm and $Ce = 1.97$. The transient duration is approximately equal to $t_2 - t_1 \approx 130$ s and the capillary number computed from the steady state is $Ca \approx 3.03 \times 10^{-4}$.

see figure 6(c). The rising regime is assumed to be stationary if the difference between the caps' velocities is less than 5%. The bubble velocity is computed as the mean velocity between the upper and the bottom cap velocities.

For a given inner radius $R$ and a fixed rotational speed $\omega$, both $Bo$ and $Ce$ are fixed. The capillary number $Ca = \mu U_b / \gamma$ is then derived from the measurement of the bubble velocity at steady state $U_b$. For the experiments reported in this section, we specify that the Bond number remains smaller than the threshold $Bo_c = 0.842$. Thus, the bubble does not move at all if the rotational speed is zero. To avoid excessively long working time for the motor, we did not operate it more than 8 hours consecutively. Considering that
Releasing long bubbles trapped in thin capillaries

<figure><img src="image_11.png" /></figure>

Figure 7. (a) Diagram (Bo, Ce), where each circle corresponds to a measurement of Ca > 0 for a given set of parameters (Bo, Ce). The red crosses indicate the couples (Bo, Ce = Ce<sub>c,exp</sub>(Bo)) for which the bubble displacement fell below our detection limit. The black circles indicate the theoretical threshold for the onset of motion and the black solid line represents the approximation (2.25) of Ce<sub>c</sub>(Bo). Below this line, the grey area indicates the region of parameters where the steady rising of a bubble is not possible according to our theoretical analysis. Finally, the shaded area corresponds to the region of parameters that is not accessible with our set-up, due to the constraints on the maximal angular velocity provided by the motor. (b) Capillary number Ca as a function of the centrifugal number Ce measured for various Bond numbers. The circles are the experimental points, and for each Bond number Bo, the solid line is the theoretical prediction (2.29), computed using the corresponding experimental value of Bo indicated in the legend. The dotted lines also represent the prediction (2.29), but for Bo ± ΔBo, where ΔBo accounts for the ±0.05 mm uncertainty on the tube inner diameters. The errorbars represent the measurement uncertainty on the bubble velocity.

with our experimental set-up we cannot precisely detect a motion smaller than 1 mm
between the start and the end of an experiment, the smallest capillary number that is
experimentally measurable is Ca<sub>min</sub> = 1.6 × 10<sup>-7</sup>. A value inferior to this limit will be
accordingly set equal to zero. The maximal rotational speed achieved by the DC motor
is ω<sub>max</sub> = 400 rad s<sup>-1</sup>. For a given tube inner radius, this sets a limit on the maximal
centrifugal number Ce<sub>max</sub> that is experimentally reachable.

2.2.2. *Experimental results and comparison with the theoretical threshold*

For comparison with the theoretical threshold for the onset of motion, we present
our experimental findings in the (Bo, Ce) diagram featured in figure 7(a). Overall, the
theoretical prediction is quantitatively consistent with the experimental results, which
reveal a rapid decay in the bubble velocity as the centrifugal number Ce approaches
the theoretical threshold Ce_c(Bo). As it was challenging to precisely determine the
experimental threshold, we endeavoured to establish a narrow range by identifying the
highest Ce ≡ Ce_c,exp for which the bubble displacement fell below our detection limit.
This lower bound is denoted by red crosses in figure 7(a), and closely aligns with the
theoretical threshold. However, the prediction is less precise for the smallest values of Bo:
the experimental threshold is downward-shifted with respect to the theoretical prediction.

Figure 7(b) reports measurements of the bubble velocity as a function of the rotational speed. The trend is satisfyingly captured by (2.29). We note, however, that close to the threshold, the measured velocities are in general higher than predicted, consistent with the downward shift of the experimental threshold mentioned above. We believe that this can be ascribed to horizontal vibrations of the tube attachment system observed while
A. *Marcotte and others*

operating the motor. As observed by Kubie (2000), the rising velocity of a Taylor bubble within a vertical tube is indeed larger when the tube is oscillated in the horizontal plane, and increases with the oscillation acceleration. This hypothesis is backed up with the photographs of rotating bubbles along their ascent, that show some asymmetry of the bubble profile with respect to the tube axis, as can be seen for instance in figure 6(b). This is compatible with the observations of Kubie (2000) under horizontal oscillations: the relative position of the bubble moves periodically from one side of the tube to the other, which thickens the lubricating film on one or the other side of the bubble, resulting in a more efficient drainage and thus in faster bubble ascent.

Despite these discrepancies, our theoretical analysis seems to provide a good estimation of the threshold for the onset of motion and a satisfying prediction for the general trend of the rising velocity as a function of the rotational speed.

In summary, rotation reduces the critical tube radius for the onset of motion and facilitates bubble ascent. From a theoretical point of view, the most appreciable effect of centrifugation is the modification of the static cap profile, while geometrical constraints stemming from the matching with the thin-film region appear remarkably unchanged from the classical case without rotation. At the same time, experiments demonstrate the thickening of the thin film surrounding the elongated part of the bubble, for increasing rotational speed. This thickening is caused by the centrifugal acceleration which induces a radial, 'gyrostatic' pressure gradient that pushes liquid towards the solid wall. We can thus interpret centrifugation as a means to tune the thickness, and thus the flow rate within the gap between the tube wall and the bubble, resulting in the lowering of the critical Bond number for the onset of motion.

An alternative and simple strategy to modify the hydrostatic pressure gradient is to tilt the tube with respect to gravity, whose effect is investigated in the next section.

3. Effect of tilt

The influence of inclination angle on the mobility of elongated bubbles was first observed by White & Beardmore (1962), who pointed out the necessity of careful positioning of pipes for precise measurement of the rising velocity. Since then, many studies have been dedicated to the motion of long bubbles in inclined pipes (e.g. Zukoski 1966; Maneri & Zuber 1974; Bendiksen 1984; Weber, Alarie & Ryan 1986; Couët & Strumolo 1987; Shosho & Rya 2001; Boucher *et al.* 2023). All studies reported a non-monotonic dependency of the rising velocity on the tilt angle: starting from a horizontal position, the velocity of elongated bubbles increases with the inclination of the pipe, reaching a maximum value. Subsequently, the velocity decreases until the vertical position is attained. These observations are reminiscent of the so-called Boycott effect (Boycott 1920; Acrivos & Herbolzheimer 1979) in the case of settling suspensions in sealed tubes, as seminally observed by Boycott (1920) with blood corpuscles sedimenting in serum, that demonstrated a several-fold increase in their sedimentation rate when the tube was inclined.

Most of these analyses are interested in the inertial regime, with large Bond numbers, and only a few studies have been dedicated to the regime close to the onset of motion, which is dominated by surface tension (low Bond number). Zukoski (1966) conducted an extensive series of experiments focusing on the velocity of elongated bubbles in tubes within a large range of Bond numbers, delving into the impact of liquid viscosity and surface tension on bubble velocity. For low Bond number (Bo = 0.870), elongated bubbles exhibited no detectable movement in horizontal or vertical positions but could rise in inclined tubes with angles ranging from 20° to 80° with the horizontal, with a maximum
Releasing long bubbles trapped in thin capillaries

<figure><img src="image_12.png" /><figcaption>Figure 8. (a) Sketch of the static cap in a tube tilted with angle $\alpha$ with respect to the horizontal plane. A Cartesian coordinate system $(x, y, z)$ is used, where $z$ is the direction aligned with the central axis of the tube. The height of the liquid-air interface is denoted as $h(x, y)$ and the gas-liquid interface meets the wall with an angle $\phi$, in all radial directions. (b) Sketch of the thin-film region, assumed to be axisymmetric. A two-dimensional Cartesian coordinate system $(\tilde{x}, \tilde{y})$ is used, where $\tilde{x}$ is the direction aligned with the central axis of the tube. The distance of the liquid-air interface from the solid wall is denoted by $y_1(\tilde{x})$. (c) Comparison between the solution of (3.1) (red solid line) and the axisymmetric solution of (2.6) with no rotation ($\omega = 0$) (pink dotted line), for tilt angle $\alpha = 90^\circ$ (vertical tube), $\phi = 0.50^\circ$ and Bond number $Bo = 0.86$.</figcaption></figure>

velocity reached for about $50^\circ$. This observation suggests that tilting the tube with respect to gravity may enable the motion of long bubbles that are stuck in a vertical configuration owing to surface tension.

In a similar context, Collicott & Manning (2014) studied the stability of a liquid mass in a tube above a capillary interface spanning the cross-section of the channel, for various contact angles and tube inclinations with respect to gravity. At fixed contact angle, they computed the critical Bond number as the threshold above which the SURFACE EVOLVER simulations do not converge to a solution of finite axial extent, thereby identifying the critical Bond number as a stability threshold for the capillary interface. Within this approach, the critical Bond number for a $0^\circ$ contact angle should correspond to the stability threshold of a long static bubble, expanding over the entire cross-section of the tube. However, difficulties of modelling perfectly wetting conditions prevented the authors from computing the critical Bond number as a function of inclination in this case.

To the best of our knowledge, there is currently no predictive analysis of the mobility enhancement of long bubbles due to tilted gravity in very narrow capillaries. In this study, we investigate how the direction of gravity affects the mobility of long bubbles in the low-$Bo$ regime, focusing on the angle-dependent threshold for the initiation of motion.

## 3.1. Theoretical prediction for the threshold and rising velocity

### 3.1.1. The three-dimensional static cap

We first introduce the equilibrium equation for the static three-dimensional shape of the upper cap of the bubble. We define a Cartesian coordinate system $(x, y, z)$, where $z$ is the direction aligned with the central axis of the tube. The gravity vector reads $g = (g \cos(\alpha), 0, -g \sin(\alpha))$, where $\alpha$ is the tilt angle of the tube ($\alpha = 90^\circ$ corresponds to a vertical tube); see figure 8(a).

The evaluation of the static interface profiles of the upper cap is based on the two-dimensional Young-Laplace equation, where length scales are non-dimensionalized with the tube radius (Manning et al. 2011; Rascón, Parry & Aarts 2017):

$$ \nabla \cdot \left( \frac{\nabla \bar{h}}{\sqrt{1 + (\nabla \bar{h})^2}} \right) = Bo(\cos(\alpha)\bar{x} - \bar{h}\sin(\alpha)), \quad (3.1) $$
where $\bar{h}(\tilde{x}, \tilde{y})$ denotes the height of the static cap; the quantity on the left-hand side is the curvature $\kappa$ of the liquid-gas interface while the term on the right-hand side corresponds to the hydrostatic contribution. Equation (3.1) is complemented with the boundary condition at the solid wall: $(\nabla h/\sqrt{1+(\nabla h)^2}) \cdot \mathbf{n} = -\cos(\phi)$, where $\mathbf{n}$ is the outwards-oriented vector normal to the tube, and the angle $\phi$ is defined similarly to that in the first part of this study, as the angle between $e_z$ and the tangent at the wall to the intersection of the liquid-gas interface with the plane $(n, e_z)$. We thus require the interface to reach the solid wall with a specified slope, which is assumed to be the same for all directions $\mathbf{n}$, for consistency with the matching with the thin-film profile, assumed in the following to be radially symmetric. Note that this differs from the first part of this study, where we imposed a vanishing radial curvature at the wall and computed the angle $\phi$ a posteriori. Here, the slope at the wall is specified as a boundary condition, with no requirement on the curvature. We acknowledge that requiring the gas-liquid interface to reach the wall for all directions $\mathbf{n}$ is somehow counterintuitive, given that for small tilt angles $\alpha$ (i.e. for a strongly inclined tube with respect to gravity), we expect the film surrounding the bubble to be thicker in the direction $x > 0$, and the liquid-gas interface to be relatively far from the solid wall in this region. However, we focus here on the vicinity of the threshold for the onset of motion, where surface tension is dominant and causes the bubble to expand in the entire fluid domain $\sqrt{x^2 + y^2} = r < R$ in all directions, as experimentally observed (see for instance figure 11d). The derivation of (3.1) can be found in Appendix B.

### 3.1.2. The thin-film region and matching

Here, we opt for a simplified description of the inner region, which is assumed to be radially symmetric. Neglecting the azimuthal curvature allows us to describe the bottom thin film with a two-dimensional, stationary Cartesian reference frame $(\tilde{x} = z - U_b t, \tilde{y})$, where $e_{\tilde{x}}$ is aligned with the tube central axis and points upwards (such that $e_{\tilde{x}} \cdot g = -g \sin(\alpha)$), and $e_{\tilde{y}}$ is the inward vector normal to the inner solid wall, such that $e_{\tilde{y}} \cdot g = -g \cos(\alpha)$; see figure 8(b).

Within the lubrication framework, the viscous flow in the thin film is driven by Laplace and hydrostatic pressure gradients. The axial velocity in the bottom thin film is accordingly written as (see Appendix A.2 for a detailed derivation)

$$u(\tilde{x}, \tilde{y}) = \frac{\gamma}{2\mu} \left[ -y_1''' + \frac{\rho g \cos(\alpha)}{\gamma} y_1' + \frac{\rho g \sin(\alpha)}{\gamma} \right] (\tilde{y}^2 - 2y_1 \tilde{y}) - U_b, \quad (3.2)$$

where $y_1$ denotes the distance of the air-liquid interface from the solid wall of the tube. Upon integration within the thin film and neglecting azimuthal variations of curvature and film thickness, the volume flux reads

$$Q \approx -2\pi R U_b y_1 - 2\pi R \frac{\gamma}{3\mu} \left( -y_1''' + \frac{\rho g \cos(\alpha)}{\gamma} y_1' + \frac{\rho g \sin(\alpha)}{\gamma} \right) y_1^3. \quad (3.3)$$

Owing to mass conservation, this flux must be equal to the volume of fluid displaced per unit time by the top of the bubble, which is $\pi R^2 U_b$. Given that $y_1/R \ll 1$, the term $-2\pi R U_b y_1$ in the flow rate expression is a negligible correction. Finally, by enforcing flux continuity with the region far from the tip, where the film thickness is uniformly constant
and equal to $b$, we derive the following thin-film equation:

$$y_1''' = \frac{\rho g \sin(\alpha)}{\gamma} \left(1 - \frac{b^3}{y_1^3}\right) + \frac{\rho g \cos(\alpha)}{\gamma} y_1'. \quad (3.4)$$

We non-dimensionalize with

$$y_1 = \eta b, \quad \tilde{x} = \zeta b (\rho g b^2 \sin(\alpha) / \gamma)^{-1/3}, \quad (3.5a,b)$$

which leads to the following ordinary differential equation:

$$\eta''' = \frac{\eta^3 - 1}{\eta^3} + a\eta', \quad (3.6)$$

where $a = \cos(\alpha) \sin(\alpha)^{-2/3} Bo^{1/3} (b/R)^{2/3}$.

Upon introduction of the parameter $a$, this equation is exactly the same as (2.12) describing the inner region in a centrifugated tube, which has been solved in § 2.1. Thus, shifting the origin to the position where $\eta'' = 0$, the distance from the wall at which the inner solution exhibits an inflection point is

$$y_1(0) = \eta(0)b = 1.10b \ll R, \quad (3.7)$$

and the slope of the inner region profile at the inflection point is given by

$$y_1'(0) = \eta(0) (\rho g b^2 \sin(\alpha) / \gamma)^{1/3} = 0.572 Bo^{1/3} \sin(\alpha)^{1/3} \left(\frac{b}{R}\right)^{2/3} > 0. \quad (3.8)$$

For the matching of the thin-film region with a two-dimensional cap, we would need to determine the value of (positive) angle $\phi$ which leads to zero curvature at the wall. For the matching with the previously introduced three-dimensional shape of the static cap, we extend this analysis by searching for the angle $\phi$ that gives rise to zero radial curvature in at least one point of the matching boundary. Note that although we imposed as a boundary condition a constant angle $\phi$ at the interface when reaching the wall, the height of the interface at the wall, and so the curvature, varies along the azimuthal direction. Since the angle $\phi$ at the point of vanishing curvature should be positive according to the matching condition (3.8), we identify the critical Bond number as the $Bo$ value for which $\phi = 0$, and the radial curvature at the wall vanishes in at least one point. For smaller Bond numbers, the geometrical constraint on the slope cannot be satisfied in at least one point of the domain.

Equation (3.1), together with its boundary conditions, is implemented in the finite-element solver COMSOL Multiphysics. We exploit fourth-order Lagrangian shape functions, solving for the height $h$ and the mean curvature $\kappa$ in a grid composed of quadrangular elements, with 10 boundary layers of 1.3 stretching factor to properly capture the curvature at the boundary. For each tilt angle $\alpha$, solutions are obtained for different values of the angle $\phi$ and Bond number $Bo$ using the built-in Newton algorithm, initialized with the zero solution. We then perform a continuation study by gradually decreasing the angle $\phi$ from $90^\circ$. Note that the boundary condition $\phi = 0^\circ$ cannot be imposed in this framework, as it implies infinite directional derivatives for the thickness. We thus study solutions in the close vicinity of $\phi = 0^\circ$ and extrapolate the retrieved behaviour for $\phi \to 0$; however, this limitation will not significantly affect the evaluation of the threshold for the bubble rise. For fixed Bond number, a convergence analysis from a characteristic size of
<figure><img src="image_13.png" /><figcaption>Figure 9. (a) Static cap profiles computed as solutions of (3.1) for various tilt angles $\alpha$, with $\phi = 0.5^\circ$ and $Bo \gtrsim Bo_c(\alpha)$. The colourbar represents the radial curvature $\kappa_r$ computed along the height profiles in the plane $y=0$. The heights of the profiles for various $\alpha$ have been translated for visualization purposes. (b) Three-dimensional static cap shape close to critical conditions $\phi = 0.5^\circ$, for various tilt angles. The colourbar represents the profile height $h - h_{max}$.</figcaption></figure>

$0.05R$ to $0.01R$ (i.e. from 3200 to 33 012 elements) showed variations of $\sim 10^{-4}$ rad in the value of $\phi$ resulting in a zero radial curvature. The numerical code for $\alpha = 90^\circ$ (i.e. for a vertical tube), $\phi = 0.5^\circ$ and $Bo = 0.86$ is compared against the axisymmetric solution of (2.6) with no rotation ($\omega = 0$). The result of the comparison is reported in figure 8(c) and a good agreement is observed.

Figure 9 shows the static cap of the bubble for different inclination angles and Bond numbers, for same angle $\phi = 0.5^\circ$. For $\alpha < 90^\circ$, the static cap is not axisymmetric: as the tilt angle increases, the apex of the cap moves towards negative $x$. Conversely, an elongated region (tongue) develops in the vicinity of the $x$ axis, in the direction $x > 0$ (i.e. in the direction of positive gravitational acceleration) and becomes longer as the tilt angle increases. The elongated region presents abnormal values of mean curvature with respect to the rest of the cap. The highest (negative) curvature is observed to be localized at the tip point of this tongue.

To obtain the critical conditions, we fix the tilt angle and the Bond number and progressively decrease the angle $\phi$. A preliminary analysis showed that the highest radial curvature is obtained at the tip point of the tongue (of coordinates $(\bar{x} = 1, \bar{y} = 0)$), in agreement with the above observations, and increases with decreasing $\phi$. For each angle $\phi$, we thus compute the radial curvature $\kappa_r = (\partial^2\bar{h}/\partial\bar{x}^2)/(1 + (\partial\bar{h}/\partial\bar{x})^2)^{3/2}$ at the extremity of the tongue. Note that $(\partial\bar{h}/\partial\bar{y})(\bar{x}, \bar{y} = 0) = 0$ because of symmetry. The angle $\phi$ is then decreased until $\kappa_r$ vanishes. This limit value (which is obtained through linear interpolation, when a change of sign is detected, of the values of curvature between two successive values of $\phi$, with a step of $9 \times 10^{-5}$ rad) is denoted $\phi_{lim}(\alpha, Bo)$. For the set of parameters $(Bo, \alpha, \phi = \phi_{lim}(\alpha, Bo))$, the liquid-air interface exhibits then an inflection point at the wall, and its tangent plane makes an angle $\phi_{lim}(\alpha, Bo)$ with the $z$ direction.
Releasing long bubbles trapped in thin capillaries

<figure><img src="image_14.png" /><figcaption>Figure 10. (a) Angle $\phi_{lim}$ for which the radial curvature of the liquid-air interface vanishes at one point at the wall as a function of $Bo$ for various values of $\alpha$ (coloured circles). For each panel, the black solid line is the linear fit $\phi_{lim}(\alpha, Bo) = \beta(\alpha)[Bo - Bo_c(\alpha)]$ performed in order to retrieve the threshold for the onset of motion $Bo_c(\alpha)$, that corresponds to $\phi_{lim} = 0$. (b) Threshold $Bo_c$ as a function of the tilt angle $\alpha$ (filled circles). The maximum extrapolation error is of the order of 0.001 and is smaller than the marker size. The black solid line represents the polynomial approximation (3.10). The open circles represent instead the threshold $Bo_c^{2D}(\alpha)$ retrieved from the matching of the thin film with a two-dimensional static cap profile. The black dotted line is a guide for the eyes. Inset: coefficient $\beta$ as a function of the tilt angle $\alpha$. The errorbars represent the 95% confidence interval. The black solid line represents the polynomial approximation (3.11).</figcaption></figure>

Repeating the same procedure varying the Bond number while fixing the tilt angle $\alpha$, we can retrieve $\phi_{lim}$ as a function of $Bo$. In the range $\phi_{lim} \in [0.5^\circ, 2^\circ]$, $\phi_{lim}(\alpha, Bo)$ varies linearly with the Bond number $Bo$, as shown in figure 10(a). For each tilt angle $\alpha$ we interpret the Bond number value at which $\phi_{lim}(\alpha, Bo_c) = 0$ as the threshold $Bo_c(\alpha)$ for the onset of motion. We retrieve this value by performing for each tilt angle $\alpha$ a linear fit of $\phi_{lim}(\alpha, Bo)$ for $\phi_{lim}$ varying between $0.5^\circ$ and $2^\circ$, and by extrapolating the value $Bo_c(\alpha)$ that corresponds to $\phi_{lim} = 0^\circ$. The fit is performed by considering at least eight points within the declared range. We verified that the threshold and slope do not vary appreciably by decreasing the number of points while keeping a constant distance between the remaining points. The slope of the fit is also a function of $\alpha$, so that overall, $\phi_{lim}(\alpha, Bo)$ is approximated by

$$ \phi_{lim}(\alpha, Bo) = \beta(\alpha)(Bo - Bo_c(\alpha)). \quad (3.9) $$

The result of this procedure is displayed in figure 10. Our study clearly indicates that the threshold for the onset of motion is lowered by tilting the tube, with a minimum that is reached for a tilt angle of $45^\circ < \alpha_{opt} < 50^\circ$. Overall, the critical Bond number and the slope $\beta$ as a function of the tilt angle are well approximated by

$$ Bo_c(\alpha) \approx 0.54 \left[ 1 + 0.5 \left( \frac{\pi}{180} \right)^2 (\alpha - 48^\circ)^2 + \left( \frac{\pi}{180} \right)^4 (\alpha - 48^\circ)^4 \right], \quad (3.10) $$

$$ \beta(\alpha) \approx -\frac{2}{3} \left(\frac{\pi}{180}\right)^2 (\alpha - 51^\circ)^2 + 0.85 \quad [\text{rad}], \quad (3.11) $$

as shown in figure 10(b). Note that in the vertical case, we retrieved values for the critical Bond number $Bo_c(\alpha = 90^\circ)$ and for the slope $\beta(\alpha = 90^\circ)$ that match Bretherton's values
### A. Marcotte and others

$(Bo_c = 0.842, \beta = 0.49 \text{ rad})$ within 0.03% and 2% of relative error, respectively, thus validating further the procedure.

We now aim at providing a prediction for the ascent velocity. The matching of the two-dimensional thin-film region with the static cap shape at the inflection point amounts to enforcing

$$0.572 Bo^{1/3} \sin(\alpha)^{1/3} \left(\frac{b}{R}\right)^{2/3} = \phi_{lim}(\alpha, Bo) = \beta(\alpha) \left[ \frac{\rho g}{\gamma} (R - 1.10b)^2 - Bo_c(\alpha) \right]. \quad (3.12)$$

The volume of fluid displaced per unit time by the tip of the bubble $\pi R^2 U_b$ being equal to the volume flux in the uniform film region, $(b/R) = (3Ca/2Bo \sin(\alpha))^{1/3}$. Expanding (3.12) up to first order in $(b/R)$ finally yields the following implicit function for the bubble velocity:

$$(Bo - Bo_c(\alpha)) \sin(\alpha) = 2.52(Bo \sin(\alpha))^{2/3} Ca^{1/3} + \frac{0.63 \sin(\alpha)}{\beta(\alpha)} (Bo \sin(\alpha))^{1/9} Ca^{2/9}. \quad (3.13)$$

As for the centrifugated case, (3.13) can be recast as an equation of the type $x^3 + a_1 x^2 + a_2 x + a_3 = 0$, by setting $x = Ca^{1/9}$.

## 3.2. Experiments on tilted bubbles

### 3.2.1. Experimental set-up and procedure

The same silicone oil used in § 2.2 is employed to partially fill capillary tubes, that are then sealed on both ends, trapping a long air bubble inside. The inner radii of the tubes vary between 1.08 and 1.96 mm, corresponding to Bond numbers in the range $Bo \in [0.53, 1.73]$. As in the previous section, the uncertainty in the tube inner diameters is 0.05 mm.

The experimental set-up is depicted in figure 11(a). The tube is attached on an aluminium arm that can be tilted by an angle $\alpha \in [0^\circ, 180^\circ]$ with respect to the horizontal plane. A light-emitting-diode panel is positioned behind the set-up for visualization purposes. Once the tilt angle is fixed, a camera records the rising motion of the bubble along the central axis of the tube. From the recorded footage, we can then retrieve the bubble velocity as previously described in § 2.2, and as illustrated in figure 11(b,c). For the narrowest tubes where the bubble velocity is the smallest (if not zero), we use time lapses instead of movies.

Unlike the case of motor-driven rotating tubes, there is in principle no limitation on the observation time, allowing the detection of much slower bubble displacements. In practice, we consider the bubble velocity to be zero if the displacement of the bubble over a week is smaller than our resolution limit of 1 mm. This implies that the smallest experimentally measurable capillary number is $Ca_{min} = 7.6 \times 10^{-9}$.

### 3.2.2. Experimental results and comparison with the theoretical prediction

Our experimental findings are summarized and compared with our theoretical predictions in figure 12. Firstly, we observe that the bubble velocity strongly depends on the tilt angle $\alpha$
<figure><img src="image_15.png" /><figcaption>Releasing long bubbles trapped in thin capillaries</figcaption></figure>

Figure 11. (a) Sketch of the experimental set-up. (b) Photographs of a long bubble inside a tube filled with silicone oil, at different and equally spaced time steps. Here, $Bo = 1$ and the tube is tilted by $α = 35°$ with respect to the horizontal axis. (c) Intensity profile as a function of time along the tube axis. To produce this image, a column of pixels aligned with the central axis of the tube is extracted from each frame of the movie. The columns are then juxtaposed to each other. The locations of the upper and lower cap as a function of time are easily identified as the two roughly parallel black curves limiting a slightly darker domain that corresponds to the position of the bubble itself. The rising velocity is given by the slope of these black lines. As in (b), $Bo = 1$ and $α = 35°$. (d) Photograph of a bubble in a tube tilted by $α = 50°$, with $Bo = 0.7$. For these parameters, the system is close to critical conditions for the onset of motion.

and reaches its maximum at $α ≈ 50°$, a value independent of $Bo$ within the range of Bond numbers investigated here, as shown in figure 12(b). Furthermore, for $Bo_c(α = 50°) = 0.5401 < Bo < 0.842 = Bo_c(α = 90°)$, tilting the tube by the appropriate angle actually enables the motion of a bubble that would otherwise be stuck in a vertical configuration, as illustrated for instance by the cases $Bo = 0.71$ and $Bo = 0.65$ reported in figure 12(b). No motion at all is observed below the threshold $Bo_c(α ≈ 50°)$. Those observations align well with our theoretical analysis.

Finally, the bubble velocity as a function of the tilt angle $α$ at low Bond numbers seems to be well described by (3.13), without any fitting parameter; see figure 12(b). We note that the agreement with the theoretical prediction appears to slightly deteriorate at larger Bond numbers. Indeed, several assumptions made in the theoretical analysis, reasonable in the vicinity of the threshold, are expected to fail in the large-$Bo$ regime. Notably, Atasi et al. (2017) measured the top and bottom film thickness around long bubbles in horizontal tubes and observed that the film asymmetry grows with the Bond number. Thus, in large capillaries, the thin-film thickness cannot be considered as uniform along the azimuthal direction: the lubricating film is indeed much thicker in the direction $x > 0$ (Zukoski 1966). Similarly, requiring the static cap profile to expand in the entire fluid domain $r < R$ is likely to become inadequate as $Bo$ increases. All together, however, the comparison tends to validate the relevance of a two-dimensional analysis to describe the thin lubricating film surrounding the bubble, even in a tilted configuration, in the low-$Bo$ regime.

It is worth mentioning that as a first attempt to describe the phenomenon, we opted for a fully two-dimensional description of the air-liquid interface. By matching the two-dimensional static cap with the thin-film profile, we obtained the threshold $Bo_c^{2D}(α)$ reported in figure 10(b). This threshold exhibits the same non-monotonic trend as a function of the tilt angle, with a minimum reached for $α_{opt} ≈ 50°$. However, it is downward-shifted with respect to the critical Bond number relying on a three-dimensional description of the static cap, which provides a much better agreement with experimental measurements; see figure 12(a). From this comparison, we conclude that while a
A. Marcotte and others

<figure><img src="image_16.png" /><figcaption>Figure 12. (a) Diagram (Bo, α) where each circle corresponds to a measurement of Ca > 0 for a given set of parameters (Bo, α). The red crosses indicate the couples (Bo, α) for which the bubble displacement fell below our detection limit. The black circles indicate the theoretical threshold for the onset of motion and the black solid line represents the polynomial approximation (3.10). Below this line, the grey area corresponds to the region of parameters where the steady rising of a bubble is not possible according to our theoretical analysis. (b) Velocity of the bubble as a function of the tilt angle α for various Bond numbers. The circles are the experimental points while for each Bond number, the solid line is described by (3.13), using the corresponding experimental value of Bo reported above each panel. No fit parameter is used here: the values of β(α) and Bo_c(α) in (3.13) are those displayed in figure 10. The dotted lines also represent the prediction (3.13), but for Bo ± ΔBo, where ΔBo accounts for the ±0.05 mm uncertainty on the tube inner diameters. Here, the marker size represents the maximal measurement uncertainty on the bubble velocity.</figcaption></figure>

simplified, two-dimensional description of the thin-film region is acceptable, a proper characterization of the phenomenon requires accounting for the three-dimensional shape of the static cap.

We can now rationalize the theoretical and experimental results: the increase in transversal acceleration due to the tilt angle tends to increase the film thickness at the tongue of the static cap, enabling higher velocities within the tube for the same axial gravity. However, tilting the tube decreases the driving buoyancy force, which in turn reduces the bubble velocity. The interplay between these two effects leads to the observed non-monotonic dependency of the rising velocity on the tilt angle. In the limit case α → 0° (horizontal tube), there is no motion within the tube since the driving force disappears.

4. Conclusion

In this study, we investigated theoretically and experimentally two different strategies aimed at enabling the motion of long air bubbles trapped in narrow, sealed capillaries partially filled with a viscous liquid. Both strategies, namely centrifugating the tube or tilting it with respect to its central axis, amount to modifying the pressure distribution in the film surrounding the bubble by means of an external force field (centrifugal force or tilted gravity). This impacts both the shape of the static cap of the bubble, and the profile of the liquid–air interface in the thin-film region. In particular, the resulting pressure gradients lead in both cases to the thickening of the lubricating film, thus enabling bubble ascent. The threshold for the onset of motion and the rising velocity above threshold as functions of the rotational speed and of tilt angle, respectively, are retrieved by the matching of the static cap and thin-film profiles, which conditions the steady ascent of the bubble. Remarkably, the matching conditions in terms of film thickness and slope at
the inflection point are in both cases the same as described in Bretherton (1961) for the classical vertical setting (without rotation). However, both centrifugation and inclination alter the inner region solution and the static cap profile, making this matching possible for smaller Bond numbers. Thus, tunable parameters such as the rotational speed or the tilt angle can effectively lower the threshold for the onset of motion, thus allowing the transport of bubbles even in very narrow capillaries.

The first part of this study was dedicated to the case of a vertical tube in rotation around its central symmetry axis. We extended Bretherton's analysis (Bretherton 1961) to account for the radial pressure gradient resulting from the tube centrifugation. By computing the shape of the tip of the bubble and solving the lubrication equation describing the thin-film region, we could derive a matching condition yielding a theoretical prediction for the ascent velocity of the bubble, together with a new threshold for the onset of motion. Our theoretical findings highlight that centrifugating the tube acts as a downward shift on the critical bubble confinement. Our experimental campaign corroborated this analysis and confirmed the relevance of this strategy for releasing bubbles trapped in very narrow capillaries.

In the second part, we explored how tilting the tube with respect to gravity could influence the transport of the bubble trapped inside. The three-dimensional static cap shape of the bubble was computed numerically, while the thin-film region was assumed to be axisymmetric. By matching these profiles at the point of vanishing radial curvature, we could derive a prediction for the steady velocity of the bubble, which can only hold if the inner radius is larger than an angle-dependent critical value. This threshold varies non-monotonically with the tilt angle, with a minimum reached about $\alpha_{opt} \approx 48^{\circ}$. Those predictions, although relying on a simplified description of the thin-film region, align well with our experimental findings.

Overall, these strategies seem well suited to many microfluidics applications where it is instrumental to get rid of trapped bubbles, without compromising the integrity of the capillary. The use of a tunable external force field provides a practical way to precisely monitor the motion of long bubbles. For further practical uses, we recall here the approximated expressions of the thresholds derived during this study:

$$ 
\left. 
\begin{array}{l}
Bo_c(Ce) \approx 0.842 - 0.295 Ce + 0.020 Ce^2 \quad \text{for centrifugated tubes and} \\
Bo_c(\alpha) \approx 0.54 \left[ 1 + 0.5 \left( \frac{\pi}{180} \right)^2 (\alpha - 48^{\circ})^2 + \left( \frac{\pi}{180} \right)^4 (\alpha - 48^{\circ})^4 \right] \quad \text{for tilted tubes.}
\end{array}
\right\} 
(4.1)
$$

At a more fundamental level, these strategies provide an interesting framework to examine the infinitely slow dynamics of pinch-off, a phenomenon explored theoretically by Lamstaes & Eggers (2017) and experimentally investigated by Dhaouadi & Kolinski (2019) in capillary tubes with inner radii $R < R_c$. For instance, starting from a moving bubble within a rotating capillary and subsequently halting the rotation offers a practical means to establish a precisely defined initial condition, from which the pinching process starts.

**Acknowledgements.** The authors thank Hervé Elettro and Gaétan Lerisson for their initial contribution to this project.

**Funding.** We acknowledge funding by the Swiss National Science Foundation under grant 200341.

**Declaration of interests.** The authors report no conflict of interest.
Author ORCIDs.

● Alice Marcotte https://orcid.org/0009-0001-1769-0621;

● Pier Giuseppe Ledda https://orcid.org/0000-0003-4435-8613;

● François Gallaire https://orcid.org/0000-0002-3029-1457;

● Ludovic Keiser https://orcid.org/0000-0002-1681-9606.

Appendix A. Derivation of the thin-film velocity profile for the centrifugal and tilted cases

A.1. Thin-film velocity profile in a centrifugated tube

Here, we derive the velocity profile in the thin-film region when the (vertical) tube is rotated around its central axis at angular frequency $\omega$. In the rotating reference frame, translating with the bubble at steady velocity $U_b e_z$, the stationary Navier–Stokes equations are written as

$$ \nabla \cdot \mathbf{u} = 0, \quad \rho[(\mathbf{u} \cdot \nabla)\mathbf{u} + 2\mathbf{\Omega} \times \mathbf{u}] = -\nabla p - \rho\mathbf{\Omega} \times (\mathbf{\Omega} \times \mathbf{r}) + \mu\Delta\mathbf{u} + \rho g, \quad (A1) $$

where $\mathbf{\Omega} = \omega e_z$ is the rotation vector.

We introduce the dimensionless variables $\bar{u}$, $\bar{z}$, $\bar{r}$ and $\bar{p}$, such that $\mathbf{u} = U_b \bar{u}$, $z = R\bar{z}$, $r = R\bar{r}$ and $p = \pi\bar{p}$, where $\pi = \rho\omega R U_b$. The dimensionless Navier–Stokes equations read

$$ \nabla \cdot \bar{\mathbf{u}} = 0, \quad \mathrm{Ro}(\bar{\mathbf{u}} \cdot \nabla)\bar{\mathbf{u}} + 2\bar{\mathbf{e}}_z \times \bar{\mathbf{u}} = -\nabla \bar{p} - \frac{1}{\mathrm{Ro}}\bar{\mathbf{e}}_z \times (\bar{\mathbf{e}}_z \times \bar{\mathbf{r}}) + E\Delta\bar{\mathbf{u}} - \frac{g}{\omega U_b}\bar{\mathbf{e}}_z, \quad (A2) $$

where $\mathrm{Ro} = U_b/\omega R$ is the Rossby number and $E = \nu/\omega R^2$ is the Ekman number. In the experiments presented in this study, $\mathrm{Ro} \ll 1$ and the nonlinear terms of the Navier–Stokes equation can therefore be neglected. Under these assumptions, and enforcing axisymmetry, the stationary Navier–Stokes equations in cylindrical coordinates are

$$ 0 = \frac{1}{r} \frac{\partial}{\partial r} (r u_r) + \frac{\partial u_z}{\partial z}, \quad (A3) $$

$$ -2\rho\omega u_{\theta} = -\frac{\partial p}{\partial r} + \rho\omega^2 r + \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_r) \right) + \frac{\partial^2 u_r}{\partial z^2} \right], \quad (A4) $$

$$ 2\rho\omega u_r = \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_{\theta}) \right) + \frac{\partial^2 u_{\theta}}{\partial z^2} \right], \quad (A5) $$

$$ 0 = -\frac{\partial p}{\partial z} + \mu \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) + \frac{\partial^2 u_z}{\partial z^2} \right] - \rho g. \quad (A6) $$

Furthermore, in cylindrical coordinates, the stress tensor is written as

$$ \sigma = \begin{pmatrix} -p + 2\mu \frac{\partial u_r}{\partial r} & \mu r \frac{\partial}{\partial r} \left( \frac{u_{\theta}}{r} \right) & \mu \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \right) \\ \mu r \frac{\partial}{\partial r} \left( \frac{u_{\theta}}{r} \right) & -p + \frac{2\mu}{r} u_r & \mu \frac{\partial u_{\theta}}{\partial z} \\ \mu \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \right) & \mu \frac{\partial u_{\theta}}{\partial z} & -p + 2\mu \frac{\partial u_z}{\partial z} \end{pmatrix}, \quad (A7) $$

and the vector normal to the interface is $\mathbf{n} = (1/\sqrt{1+r'_1(z)^2})(-1, 0, r'_1(z))^T$.
<figure><img src="image_17.png" /><figcaption>Figure 13. Sketch of the upper part of a long bubble in a sealed vertical tube, rotating around its symmetry axis at angular velocity $\omega$. The flow around the bubble is first described in cylindrical coordinates $(r, \theta, z)$, where $e_z$ is aligned with the tube axis. We focus on the thin-film region, which can be considered as planar instead of annular, and describe then the lubricating film in Cartesian coordinates $(x, y = R - r, z)$.</figcaption></figure>

Therefore, the dynamic and kinematic boundary conditions at the fluid–air interface are

$$ \gamma\kappa = (p - p_{\text{air}})(1 + r'_{1}(z)^{2}) - 2\mu \left[ \frac{\partial u_{r}}{\partial r} + r'_{1}(z)^{2} \frac{\partial u_{z}}{\partial z} \right] + 2\mu r'_{1}(z) \left[ \frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r} \right], \quad (\text{A8}) $$

$$ 0 = 2r'_{1}(z) \left[ \frac{\partial u_{r}}{\partial r} - \frac{\partial u_{z}}{\partial z} \right] + (1 - r'_{1}(z)^{2}) \left[ \frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r} \right], \quad (\text{A9}) $$

where $\kappa = -1/(r_1(z)\sqrt{1+r'_1(z)^2}) + r''_1(z)/(1+r'_1(z)^2)^{3/2}$ is the curvature. Finally, the no-slip boundary condition at the solid wall implies $u_z(r=R) = -U_b$.

### A.1.1. Change of coordinates

In the thin-film region, the film thickness is very small compared with the radius, so that the flow can be treated as if the region were planar, instead of annular. Accordingly, we describe the flow in the Cartesian coordinate system $(x, y = R - r, z)$; see figure 13. Furthermore, we introduce the modified pressure field: $P = p + \rho\omega^2 R y + \rho g z$. In this system of coordinates, the Navier–Stokes equation become

$$ 0 = \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} - \frac{u_y}{R - y}, \quad (\text{A10}) $$

$$ -2\rho\omega u_x = \frac{\partial P}{\partial y} - \rho\omega^2 y + \mu \left[ \frac{1}{R-y} \frac{\partial u_y}{\partial y} + \frac{u_y}{(R-y)^2} - \frac{\partial^2 u_y}{\partial y^2} - \frac{\partial^2 u_y}{\partial z^2} \right], \quad (\text{A11}) $$

$$ -2\rho\omega u_y = \mu \left[ -\frac{1}{R-y} \frac{\partial u_x}{\partial y} - \frac{u_x}{(R-y)^2} + \frac{\partial^2 u_x}{\partial y^2} + \frac{\partial^2 u_x}{\partial z^2} \right], \quad (\text{A12}) $$

$$ 0 = -\frac{\partial P}{\partial z} + \mu \left[ -\frac{1}{R-y} \frac{\partial u_z}{\partial y} + \frac{\partial^2 u_z}{\partial y^2} + \frac{\partial^2 u_z}{\partial z^2} \right]. \quad (\text{A13}) $$
Likewise, the dynamic and kinematic boundary conditions become

$$ \gamma\kappa = (P - P_{\text{air}})(1 + y'_{1}(z)^{2}) - 2\mu \left[ \frac{\partial u_{y}}{\partial y} + y'_{1}(z)^{2} \frac{\partial u_{z}}{\partial z} \right] + 2\mu y'_{1}(z) \left[ \frac{\partial u_{y}}{\partial z} + \frac{\partial u_{z}}{\partial y} \right], \quad (\text{A14}) $$

$$ 0 = 2y'_{1}(z) \left[ \frac{\partial u_{y}}{\partial y} - \frac{\partial u_{z}}{\partial z} \right] + (1 - y'_{1}(z)^{2}) \left[ \frac{\partial u_{y}}{\partial z} + \frac{\partial u_{z}}{\partial y} \right], \quad (\text{A15}) $$

$$ \text{where } \kappa = -1/((R - y_1(z))\sqrt{1 + y_1(z)^2}) - y_1''(z)/(1 + y_1(z)^2)^{3/2}. $$

### A.1.2. Lubrication approximation

We non-dimensionalize as follows: $u_z = U_b \bar{u}_z$, $u_y = U_y \bar{u}_y$, $u_x = U_x \bar{u}_x$, $P = P_0 \bar{P}$, $y = b\bar{y}$, $y_1 = b\bar{y}_1$, $z = R\bar{z}$, where $\epsilon = b/R \ll 1$. According to the least degeneracy principle applied to the mass conservation equation, $U_y = \epsilon U_b$ and the mass conservation equation becomes $\partial \bar{u}_y/\partial \bar{y} + \partial \bar{u}_z/\partial \bar{z} = 0$.

Furthermore, upon introduction of the dimensionless fields and variables, the momentum conservation equations are written as

$$ -2\rho\omega U_x \bar{u}_x = \frac{P_0}{b} \frac{\partial \bar{P}}{\partial \bar{y}} - \rho\omega^2 b \bar{y} + \epsilon \frac{\mu U_b}{b^2} \left[ \frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_y}{\partial \bar{y}} + \epsilon^2 \frac{\bar{u}_y}{(1 - \epsilon \bar{y})^2} - \frac{\partial^2 \bar{u}_y}{\partial \bar{z}^2} - \epsilon^2 \frac{\partial^2 \bar{u}_y}{\partial \bar{z}^2} \right], \quad (\text{A16}) $$

$$ -2\rho\omega\epsilon U_b \bar{u}_b = \frac{\mu U_x}{b^2} \left[ -\frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_x}{\partial \bar{y}} - \epsilon^2 \frac{\bar{u}_x}{(1 - \epsilon \bar{y})^2} + \frac{\partial^2 \bar{u}_x}{\partial \bar{z}^2} + \epsilon^2 \frac{\partial^2 \bar{u}_x}{\partial \bar{z}^2} \right], \quad (\text{A17}) $$

$$ 0 = -\frac{P_0}{R} \frac{\partial \bar{P}}{\partial \bar{z}} + \frac{\mu U_b}{b^2} \left[ -\frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_z}{\partial \bar{y}} + \frac{\partial^2 \bar{u}_z}{\partial \bar{y}^2} + \epsilon^2 \frac{\partial^2 \bar{u}_z}{\partial \bar{z}^2} \right]. \quad (\text{A18}) $$

The least degeneracy principle applied to the momentum conservation equations along the x and z axes implies that $U_x = 2\epsilon\rho\omega U_b b^2/\mu$ and $P_0 = \mu U_b/b\epsilon$.

Thus, the momentum conservation equations along the y axis become

$$ -4\epsilon Re \bar{u}_x = \frac{Ca}{\epsilon^4 Ce} \frac{\partial \bar{P}}{\partial \bar{y}} - \bar{y} + \frac{Ca}{\epsilon^2 Ce} \left[ \frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_y}{\partial \bar{y}} + \epsilon^2 \frac{\bar{u}_y}{(1 - \epsilon \bar{y})^2} - \frac{\partial^2 \bar{u}_y}{\partial \bar{z}^2} - \epsilon^2 \frac{\partial^2 \bar{u}_y}{\partial \bar{z}^2} \right], \quad (\text{A19}) $$

with $Re = \rho U_b b/\mu$, $Ca = \mu U_b/\gamma$ and $Ce = \rho\omega^2 R^3/\gamma$. From the volume conservation constraint $\rho g b^3/3\mu U_b = R/2$, we deduce that $Ca/Ce \sim \epsilon^3 (2g/3\omega^2 R) = O(\epsilon^3)$. Thus, the system of equations reduces to

$$ 0 = \frac{\partial \bar{u}_y}{\partial \bar{y}} + \frac{\partial \bar{u}_z}{\partial \bar{z}}, \quad 0 = \frac{\partial \bar{P}}{\partial \bar{y}}, \quad 0 = -\frac{\partial \bar{P}}{\partial \bar{z}} + \frac{\partial^2 \bar{u}_z}{\partial \bar{y}^2}. \quad (\text{A20a-c}) $$
The dynamic and kinematic boundary conditions are

$$ \gamma\kappa = \frac{\mu U_b}{b\epsilon}(\bar{P} - \bar{P}_{air})(1 + \epsilon^2 \overline{y_1'}(\bar{z})^2) - 2\epsilon \frac{\mu U_b}{b} \left[ \frac{\partial \overline{u_y}}{\partial \bar{y}} + \epsilon^2 \overline{y_1'}(\bar{z})^2 \frac{\partial \overline{u_z}}{\partial \bar{z}} \right] \\ + 2\epsilon \frac{\mu U_b}{b} \overline{y_1'}(\bar{z}) \left[ \epsilon^2 \frac{\partial \overline{u_y}}{\partial \bar{z}} + \frac{\partial \overline{u_z}}{\partial \bar{y}} \right], \qquad (A21) $$

$$ 0 = 2\epsilon^2 \overline{y_1'}(\bar{z}) \left[ \frac{\partial \overline{u_y}}{\partial \bar{y}} - \frac{\partial \overline{u_z}}{\partial \bar{z}} \right] + (1 - \epsilon^2 \overline{y_1'}(\bar{z})^2) \left[ \epsilon^2 \frac{\partial \overline{u_y}}{\partial \bar{z}} + \frac{\partial \overline{u_z}}{\partial \bar{y}} \right]. \qquad (A22) $$

Thus, including the no-slip boundary condition at the solid wall, the boundary conditions at leading order are

$$ \bar{P} - \bar{P}_{air} = \gamma\kappa \frac{\epsilon b}{\mu U_b}, \quad \frac{\partial \overline{u_z}}{\partial \bar{y}} = 0, \quad \overline{u_z}(\bar{y} = 0) = -1. \qquad (A23a-c) $$

Finally, going back to the dimensional form, and reintroducing the original pressure field $p = P - \rho\omega^2 R y - \rho g z$, the full problem is written as

$$ 0 = \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z}, \quad \frac{\partial p}{\partial y} = -\rho\omega^2 R, \quad \frac{\partial p}{\partial z} = \mu \frac{\partial^2 u_z}{\partial y^2} - \rho g, \qquad (A24a-c) $$

and is complemented by the following boundary conditions:

$$ p(y = y_1, z) - p_{air} = \gamma\kappa, \quad \left. \frac{\partial u_z}{\partial y} \right|_{y=y_1} = 0, \quad u_z(y = 0, z) = -U_b. \qquad (A25a-c) $$

The integration of the pressure field is straightforward and leads to

$$ p(y, z) = p_{air} + \gamma\kappa + \rho\omega^2 R(y_1(z) - y). \qquad (A26) $$

Finally, by injecting this pressure field in the axial component of the momentum equation, we can derive the following equation for the velocity in the thin film:

$$ \mu \frac{\partial^2 u_z}{\partial y^2} = \gamma \kappa' + \rho \omega^2 R y_1' + \rho g, \qquad (A27) $$

which results in

$$ u_z(y, z) = -U_b + \frac{\gamma}{2\mu} \left( \kappa' + \frac{\rho\omega^2 R}{\gamma} y_1' + \frac{\rho g}{\gamma} \right) (y^2 - 2y_1 y). \qquad (A28) $$

### A.2. Thin-film velocity profile in a tilted tube

We now aim at deriving the velocity profile in the thin-film region when the tube is tilted by an angle $\alpha$. Since the Reynolds number $Re = \rho U_b b / \mu$ characterizing the flow in the thin-film region is very small, we can safely neglect the effect of inertia. In the stationary cylindrical system of coordinates $(r, \theta, z)$, translating with the bubble at steady velocity
<figure><img src="image_18.png" /><figcaption>Figure 14. (a) Sketch of the upper part of a long bubble in a sealed tube tilted by an angle $\alpha$ with respect to the horizontal plane. The flow in the thin lubricating film is described in cylindrical coordinates $(r, \theta, z)$, where $e_z$ is the direction aligned with the tube axis, such that $e_z \cdot g = -g \sin(\alpha)$. The origin of $\theta$ is chosen such that $e_r(\theta = 0) \cdot g = g \cos(\alpha)$. (b) Sketch of the cross-section of the channel and of the bubble. We focus on the region in the vicinity of the plane $\theta = 0$, described in Cartesian coordinates $(x, y = R - r, z)$, where $e_x \cdot g = 0$, $e_y \cdot g = -g \cos(\alpha)$ and $e_z \cdot g = -g \sin(\alpha)$.</figcaption></figure>

$U_b e_z$, where $z$ is aligned with the central axis of the (tilted) tube (see figure 14a), the Navier-Stokes equations read

$$0 = \frac{1}{r} \frac{\partial}{\partial r} (r u_r) + \frac{1}{r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial u_z}{\partial z}, \qquad (A29)$$

$$0 = -\frac{\partial p}{\partial r} + \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_r) \right) + \frac{1}{r^2} \frac{\partial^2 u_r}{\partial \theta^2} - \frac{2}{r^2} \frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial^2 u_r}{\partial z^2} \right] + \rho g \cos(\theta) \cos(\alpha), \qquad (A30)$$

$$0 = -\frac{\partial p}{\partial \theta} + \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_{\theta}) \right) + \frac{1}{r^2} \frac{\partial^2 u_{\theta}}{\partial \theta^2} + \frac{2}{r^2} \frac{\partial u_r}{\partial \theta} + \frac{\partial^2 u_{\theta}}{\partial z^2} \right] - \rho g \sin(\theta) \cos(\alpha), \qquad (A31)$$

$$0 = -\frac{\partial p}{\partial z} + \mu \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 u_z}{\partial \theta^2} + \frac{\partial^2 u_z}{\partial z^2} \right] - \rho g \sin(\alpha). \qquad (A32)$$

We know from our analysis of the three-dimensional cap profile that the matching with the thin-film region profile should be imposed at the tip of the tongue exhibited by the static cap, i.e. at the point of coordinates $(r = R, \theta = 0, h(R, 0))$. We thus restrict the study of the thin-film solution to the plane ($\theta = 0$); see figure 14(b). By assuming a vanishing azimuthal curvature $\sim 1/R$, the region of size $\sim R d\theta$ in the close vicinity of $\theta = 0$ can be considered infinite. Yet, in the vicinity of $\theta = 0$, the derivative with respect to $\theta$ should vanish by symmetry. Therefore, in this region, the Navier-Stokes equations reduce to

$$0 = \frac{1}{r} \frac{\partial}{\partial r} (r u_r) + \frac{\partial u_z}{\partial z}, \qquad (A33)$$

$$0 = -\frac{\partial p}{\partial r} + \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_r) \right) + \frac{\partial^2 u_r}{\partial z^2} \right] + \rho g \cos(\alpha), \qquad (A34)$$

$$0 = \mu \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r u_{\theta}) \right) + \frac{\partial^2 u_{\theta}}{\partial z^2} \right], \qquad (A35)$$

$$0 = -\frac{\partial p}{\partial z} + \mu \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) + \frac{\partial^2 u_z}{\partial z^2} \right] - \rho g \sin(\alpha). \qquad (A36)$$
These equations are complemented by the following dynamic and kinematic boundary conditions:

$$ \gamma\kappa = (p - p_{air})(1 + r'_{1}(z)^{2}) - 2\mu \left[ \frac{\partial u_{r}}{\partial r} + r'_{1}(z)^{2} \frac{\partial u_{z}}{\partial z} \right] + 2\mu r'_{1}(z) \left[ \frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r} \right], \quad (A37) $$

$$ 0 = 2r'_{1}(z) \left[ \frac{\partial u_{r}}{\partial r} - \frac{\partial u_{z}}{\partial z} \right] + (1 - r'_{1}(z)^{2}) \left[ \frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r} \right], \quad (A38) $$

where $\kappa = -1/(r_1(z)\sqrt{1+r_1'(z)^2}) + r_1''(z)/(1+r_1'(z)^2)^{3/2}$ is the curvature, and, by the no-slip boundary condition at the solid wall, $u_z(r=R) = -U_b$.

### A.2.1. Change of coordinates

As before, we neglect curvature in the azimuthal direction, and describe the thin-film region within the Cartesian coordinate system $(x, y = R - r, z)$; see figure 14(b). Furthermore, we introduce the modified pressure field: $P = p + \rho g \cos(\alpha)y + \rho g \sin(\alpha)z$. In this system of coordinates, the mass conservation and momentum conservation equations along the y and z directions become

$$ 0 = \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} - \frac{u_y}{R - y}, \quad (A39) $$

$$ 0 = \frac{\partial P}{\partial y} + \mu \left[ \frac{1}{R-y} \frac{\partial u_y}{\partial y} + \frac{u_y}{(R-y)^2} - \frac{\partial^2 u_y}{\partial y^2} - \frac{\partial^2 u_y}{\partial z^2} \right], \quad (A40) $$

$$ 0 = -\frac{\partial P}{\partial z} + \mu \left[ -\frac{1}{R-y} \frac{\partial u_z}{\partial y} + \frac{\partial^2 u_z}{\partial y^2} + \frac{\partial^2 u_z}{\partial z^2} \right]. \quad (A41) $$

Likewise, the dynamic and kinematic boundary conditions become

$$ \gamma\kappa = (P - P_{air})(1 + y'_{1}(z)^{2}) - 2\mu \left[ \frac{\partial u_y}{\partial y} + y'_{1}(z)^{2} \frac{\partial u_z}{\partial z} \right] + 2\mu y'_{1}(z) \left[ \frac{\partial u_y}{\partial z} + \frac{\partial u_z}{\partial y} \right], \quad (A42) $$

$$ 0 = 2y'_{1}(z) \left[ \frac{\partial u_y}{\partial y} - \frac{\partial u_z}{\partial z} \right] + (1 - y'_{1}(z)^2) \left[ \frac{\partial u_y}{\partial z} + \frac{\partial u_z}{\partial y} \right], \quad (A43) $$

$$ \text{where } \kappa = -1/((R - y_1(z))\sqrt{1 + y_1'(z)^2}) - \frac{y_1''(z)}{(1 + y_1'(z)^2)^{3/2}}. $$

### A.2.2. Lubrication approximation

We non-dimensionalize as follows: $u_z = U_b \bar{u}_z$, $u_y = U_y \bar{u}_y$, $P = P_0 \bar{P}$, $y = b \bar{y}$, $y_1 = b \bar{y}_1$, $z = R \bar{z}$, where $\epsilon = b/R \ll 1$. According to the least degeneracy principle applied to the mass conservation equation, $U_y = \epsilon U_b$ and the mass conservation equation becomes $0 = \partial \bar{u}_y / \partial \bar{y} + \partial \bar{u}_z / \partial \bar{z}$.
Furthermore, upon introduction of the dimensionless fields and variables, the momentum conservation equations along the y and z directions are written as

$$0 = \frac{P_0}{b} \frac{\partial \bar{P}}{\partial \bar{y}} + \frac{\epsilon \mu U_b}{b^2} \left[ \frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_y}{\partial \bar{y}} + \frac{\epsilon^2}{(1 - \epsilon \bar{y})^2} \bar{u}_y - \frac{\partial^2 \bar{u}_y}{\partial \bar{y}^2} - \epsilon^2 \frac{\partial^2 \bar{u}_y}{\partial \bar{z}^2} \right], \quad (A44)$$

$$0 = -\frac{P_0}{R} \frac{\partial \bar{P}}{\partial \bar{z}} + \frac{\mu U_b}{b^2} \left[ -\frac{\epsilon}{1 - \epsilon \bar{y}} \frac{\partial \bar{u}_z}{\partial \bar{y}} + \frac{\partial^2 \bar{u}_z}{\partial \bar{y}^2} + \epsilon^2 \frac{\partial^2 \bar{u}_z}{\partial \bar{z}^2} \right]. \quad (A45)$$

The least degeneracy principle applied to the momentum conservation equation along the z axis implies that $P_0 = \mu U_b/b\epsilon$. At leading order, the problem reduces then to

$$0 = \frac{\partial \bar{u}_y}{\partial \bar{y}} + \frac{\partial \bar{u}_z}{\partial \bar{z}}, \quad 0 = \frac{\partial \bar{P}}{\partial \bar{y}}, \quad 0 = -\frac{\partial \bar{P}}{\partial \bar{z}} + \frac{\partial^2 \bar{u}_z}{\partial \bar{y}^2}. \quad (A46a-c)$$

The dynamic and kinematic boundary conditions are

$$\gamma \kappa = \frac{\mu U_b}{b\epsilon} (\bar{P} - \bar{P}_{air}) (1 + \epsilon^2 \bar{y}_1'(\bar{z})^2) - 2\epsilon \frac{\mu U_b}{b} \left[ \frac{\partial \bar{u}_y}{\partial \bar{y}} + \epsilon^2 \bar{y}_1'(\bar{z})^2 \frac{\partial \bar{u}_z}{\partial \bar{z}} \right] \\ + 2\epsilon \frac{\mu U_b}{b} \bar{y}_1'(\bar{z}) \left[ \epsilon^2 \frac{\partial \bar{u}_y}{\partial \bar{z}} + \frac{\partial \bar{u}_z}{\partial \bar{y}} \right], \quad (A47)$$

$$0 = 2\epsilon^2 \bar{y}_1'(\bar{z}) \left[ \frac{\partial \bar{u}_y}{\partial \bar{y}} - \frac{\partial \bar{u}_z}{\partial \bar{z}} \right] + (1 - \epsilon^2 \bar{y}_1'(\bar{z})^2) \left[ \epsilon^2 \frac{\partial \bar{u}_y}{\partial \bar{z}} + \frac{\partial \bar{u}_z}{\partial \bar{y}} \right]. \quad (A48)$$

Therefore at leading order, and including the no-slip boundary condition at the solid wall, the boundary conditions are written as

$$\bar{P} - \bar{P}_{air} = \gamma \kappa \frac{\epsilon b}{\mu U_b}, \quad \frac{\partial \bar{u}_z}{\partial \bar{y}} = 0, \quad \bar{u}_z(\bar{y} = 0) = -1. \quad (A49a-c)$$

Finally, going back to the dimensional form, and reintroducing the original pressure field $p = P - \rho g \cos(\alpha)y - \rho g \sin(\alpha)z$, the full problem reduces to

$$0 = \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z}, \quad \frac{\partial p}{\partial y} = -\rho g \cos(\alpha), \quad \frac{\partial p}{\partial z} = \mu \frac{\partial^2 u_z}{\partial y^2} - \rho g \sin(\alpha), \quad (A50a-c)$$

and is complemented by the following boundary conditions:

$$p(y = y_1, z) - p_{air} = \gamma \kappa, \quad \left. \frac{\partial u_z}{\partial y} \right|_{y=y_1} = 0, \quad u_z(y = 0, z) = -U_b. \quad (A51a-c)$$

The pressure field integrates straightforwardly into

$$p(y, z) = p_{air} + \gamma \kappa + \rho g \cos(\alpha)(y_1(z) - y). \quad (A52)$$

Finally, by injecting this pressure field in the axial component of the momentum equation, we can derive the following equation for the velocity in the thin film:

$$\mu \frac{\partial^2 u_z}{\partial y^2} = \gamma \kappa' + \rho g \cos(\alpha) y_1' + \rho g \sin(\alpha), \quad (A53)$$

which leads to the velocity profile:

$$u_z(y, z) = -U_b + \frac{\gamma}{2\mu} \left( \kappa' + \frac{\rho g \cos(\alpha)}{\gamma} y_1' + \frac{\rho g \sin(\alpha)}{\gamma} \right) (y^2 - 2y_1 y). \quad (A54)$$
Appendix B. Derivation of the equilibrium equation for the static cap

Following previous works (see e.g. Manning *et al.* 2011; Lubbers *et al.* 2014; Rascón *et al.* 2017), we derive from energy principles the three-dimensional equilibrium equation for the equilibrium of the static cap. We introduce the coordinate system (x, y, z), with the z axis aligned with the central tube axis and the x axis aligned with the gravity component normal to the z axis (so that $e_x \cdot g = g \cos(\alpha)$; see figure 8a). The location of the air-liquid interface is denoted by $h(x, y)$.

The Gibbs free energy associated with the cap interface can be written as

$$E(h) = \gamma A + G, \qquad (B1)$$

where the first term on the right-hand side is the surface energy. As in Rascón *et al.* (2017), the surface A is computed as

$$A = \int_{\Omega} \sqrt{1 + (\nabla h)^2} \, dx \, dy, \qquad (B2)$$

with Ω the cross-section of the capillary. The term G represents in turn the gravitational potential energy, which in general form reads (Pitts 1973)

$$G = - \int_V \rho g \cdot r \, dx \, dy \, dz, \quad r = (x, y, z), \quad g = g(\cos \alpha, 0, -\sin \alpha), \qquad (B3)$$

where the volume V is given by

$$V = \int_{\Omega} dx \, dy h(x, y). \qquad (B4)$$

Upon introduction of the Lagrange multiplier λ to ensure volume conservation, the functional to be minimized to obtain equilibrium reads

$$F(h) = \gamma \int_{\Omega} \sqrt{1 + (\nabla h)^2} \, dx \, dy + \rho g \int_{V} (-x \cos \alpha + z \sin \alpha) \, dx \, dy \, dz - \lambda \int_{V} dxdy \, dz. \qquad (B5)$$

Upon integration along the z direction between 0 and h:

$$F(h) = \gamma \int_{\Omega} \sqrt{1 + (\nabla h)^2} \, dx \, dy + \rho g \int_{\Omega} \left( -x \cos \alpha + \frac{1}{2} h \sin \alpha \right) h \, dx \, dy - \lambda \int_{\Omega} h \, dx \, dy. \qquad (B6)$$

Formal minimization of the functional $F(h)$ with respect to $h$ leads to the following partial differential equation:

$$\gamma \nabla \cdot \left( \frac{\nabla h}{\sqrt{1 + (\nabla h)^2}} \right) = \rho g(\cos(\alpha)x - h \sin(\alpha)) + \lambda, \qquad (B7)$$

with the constant slope condition at the wall:

$$\frac{\nabla h}{\sqrt{1 + (\nabla h)^2}} \cdot \mathbf{n} = - \cos(\phi), \qquad (B8)$$

where n is the unit exterior normal to the tube wall. The value of λ can be set by integrating the resulting equilibrium equation within the whole domain, leading to the following expression:

$$\lambda = -\frac{2\gamma}{R} \cos(\phi) + \rho g h_0 \sin \alpha, \qquad (B9)$$

where $h_0 = V/\Omega$ is the reference average value of the static cap height. By imposing $\lambda = 0$, the reference height reads $h_0 = 2\ell_c^2 \cos(\phi)/R \sin \alpha$, reminiscent of the well-known
A. Marcotte and others

Jurin height. Upon non-dimensionalization with the tube radius *R*, one obtains the equilibrium equation (3.1) reported in the main text.

REFERENCES

ACRIVOS, A. & HERBOLZHEIMER, E. 1979 Enhanced sedimentation in settling tanks with inclined walls. *J. Fluid Mech.* **92**, 435–457.

ATASI, O., KHODAPARAST, S., SCHEID, B. & STONE, H.A. 2017 Effect of buoyancy on the motion of long bubbles in horizontal tubes. *Phys. Rev. Fluids* **2**, 094304.

BAROUD, C.N., GALLLAIRE, F. & DANGLA, R. 2010 Dynamics of microfluidic droplets. *Lab on a Chip* **10**, 2032–2045.

BAROUD, C.N., DE SAINT VINCENT, M.R. & DELVILLE, J.-P. 2007 An optical toolbox for total control of droplet microfluidics. *Lab on a Chip* **7**, 1029–1033.

BENDIKSEN, K.H. 1984 An experimental investigation of the motion of long bubbles in inclined tubes. *Intl J. Multiphase Flow* **10**, 467–483.

BI, Q.C. & ZHAO, T.S. 2001 Taylor bubbles in miniaturized circular and noncircular channels. *Intl J. Multiphase Flow* **27**, 561–570.

BICO, J. & QUÉRÉ, D. 2002 Rise of liquids and bubbles in angular capillary tubes. *J. Colloid Interface Sci.* **247**, 162–166.

BICO, J., TORDEUX, C. & QUÉRÉ, D. 2001 Rough wetting. *EPL* **55**, 214.

BOUCHER, A., KARP, J., BELT, R. & LINÉ, A. 2023 Dynamics of elongated bubbles in slightly inclined pipes with viscous fluids. *J. Fluid Mech.* **977**, A40.

BOYCOTT, A.E. 1920 Sedimentation of blood corpuscles. *Nature* **104**, 532–532.

BRANNOCK, D. & KUBIE, J. 1996 Velocity of long bubbles in oscillating vertical pipes. *Intl J. Multiphase Flow* **22**, 1031–1034.

BRETHERTON, F.P. 1961 The motion of long bubbles in tubes. *J. Fluid Mech.* **10**, 166–188.

BRODRIBB, T.J., BIENAIMÉ, D. & MARMOTTANT, P. 2016 Revealing catastrophic failure of leaf networks under stress. *Proc. Natl Acad. Sci. USA* **113**, 4865–4869.

CHENG, H.B. & LU, Y.W. 2014 Applications of textured surfaces on bubble trapping and degassing for microfluidic devices. *Microfluid Nanofluid* **15**, 855–862.

COLLICOTT, S.H. & MANNING, R.E. 2014 Capillary stability in a tilted circular cylinder. *Microgravity Sci. Technol.* **25**, 319–326.

CORAPCIOGLU, M.Y., CIHAN, A. & DRAZENOVIC, M. 2004 Rise velocity of an air bubble in porous media: theoretical studies. *Water Resour. Res.* **40**, W04214.

COUÈT, B. & STRUMOLO, G.S. 1987 The effects of surface tension and tube inclination on a two-dimensional rising bubble. *J. Fluid Mech.* **184**, 1–14.

DAVIES, R.M. & TAYLOR, G.I. 1950 The mechanics of large bubbles rising through extended liquids and through liquids in tubes. *Proc. R. Soc. Lond. A* **200**, 375–390.

DHAOUADI, W. & KOLINSKI, J.M. 2019 Bretherton's buoyant bubble. *Phys. Rev. Fluids* **4**, 123601.

DRELICH, J., FANG, C. & WHITE, C.L. 2002 Measurement of interfacial tension in fluid-fluid systems. *Encycl. Surf. Colloid Sci.* **3**, 3158–3163.

DUMITRESCU, D.T. 1943 Strömung an einer Luftblase im senkrechten Rohr. *Z. Angew. Math. Mech.* **23** (3), 139–149.

EMSLIE, A.G., BONNER, F.T. & PECK, L.G. 1958 Flow of a viscous liquid on a rotating disk. *J. Appl. Phys.* **29** (5), 858–862.

FAIRBROTHER, F. & STUBBS, A.E. 1935 Studies in electro-endosmosis. Part VI. The 'bubble-tube' method of measurement. *J. Chem. Soc.* **1**, 527–529.

FUNADA, T., JOSEPH, D.D., MAEHARA, T. & YAMASHITA, S. 2005 Ellipsoidal model of the rise of a Taylor bubble in a round tube. *Intl J. Multiphase Flow* **31**, 473–491.

GARY, J.H., HANDWERK, J.H., KAISER, M.J. & GEDDES, D. 2007 *Petroleum Refining: Technology and Economics*. CRC Press.

GIBSON, A.H. 1913 On the motion of long air-bubbles in a vertical tube. *Lond. Edinb. Phil. Mag.* **26**, 952–965.

GROELL, R., SCHAFFLER, G.J. & RIENMUELLER, R. 1997 The peripheral intravenous cannula: a cause of venous air embolism. *Am. J. Med. Sci.* **314**, 300–302.

GUO, L., LIU, Y., RAN, P., WANG, G., SHAN, J., LI, X., LIU, C. & LI, J. 2022 A bioinspired bubble removal method in microchannels based on angiosperm xylem embolism repair. *Microsyst. Nanoeng.* **8**, 34.

JAMBON-PUILLET, E., ROYER PIÉCHAUD, M. & BRUN, P.-T. 2021 Elastic amplification of the Rayleigh–Taylor instability in solidifying melts. *Proc. Natl Acad. Sci.* **118** (10), e2020701118.
JENSEN, M.J., GORANOVIĆ, G. & BRUUS, H. 2004 The clogging pressure of bubbles in hydrophilic microchannel contractions. *J. Micromech. Microeng.* **14**, 876.

KAIGALA, G.V., LOVCHIK, R.D., DRECHSLER, U. & DELAMARCHE, E. 2011 A vertical microfluidic probe. *Langmuir* **27**, 5686–5693.

KROGMAN, K.C., DRUFFEL, T. & SUNKARA, M.K. 2005 Ultra-high mobility transparent organic thin film transistors grown by an off-centre spin-coating method. *Nanotechnology* **16** (7), S338.

KUBIE, J. 2000 Velocity of long bubbles in horizontally oscillating vertical pipes. *Intl J. Multiphase Flow* **26**, 339–349.

LAMSTAES, C. & EGGERS, J. 2017 Arrested bubble rise in a narrow tube. *J. Stat. Phys.* **167**, 656–682.

LI, Z., LI, G., LI, Y., CHEN, Y., LI, J. & CHEN, H. 2021 Flow field around bubbles on formation of air embolism in small vessels. *Proc. Natl Acad. Sci. USA* **18**, e2025406118.

LITTERST, C., ECCARIUS, S., HEBLING, C., ZENGERLE, R. & KOLTAY, P. 2006 Increasing µDMFC efficiency by passive CO₂ bubble removal and discontinuous operation. *J. Micromech. Microeng.* **16**, 248–253.

LUBBERS, L.A., WEIJS, J.H., BOTTO, L., DAS, S., ANDREOTTI, B. & SNOEIJER, J.H. 2014 Drops on soft solids: free energy and double transition of contact angles. *J. Fluid Mech.* **747**, R1.

MADANI, S., CABALLINA, O. & SOUHAR, M. 2009 Unsteady dynamics of Taylor bubble rising in vertical oscillating tubes. *Intl J. Multiphase Flow* **35**, 363–375.

MADANI, S., CABALLINA, O. & SOUHAR, M. 2012 Some investigations on the mean and fluctuating velocities of an oscillating Taylor bubble. *Nucl. Engng Des.* **252**, 135–143.

MAGNINI, M., KHODAPARAST, S., MATAR, O.K., STONE, H.A. & THOME, J.R. 2019 Dynamics of long gas bubbles rising in a vertical tube in a cocurrent liquid flow. *Phys. Rev. Fluids* **4**, 023601.

MANERI, C.C. & ZUBER, N. 1974 An experimental study of plane bubbles rising at inclination. *Intl J. Multiph. Flow* **1**, 623–645.

MANNING, R., COLLICOTT, S. & FINN, R. 2011 Occlusion criteria in tubes under transverse body forces. *J. Fluid Mech.* **682**, 397–414.

MARTHELOT, J., STRONG, E.F. & BRUN, P.-T. 2018 Designing soft materials with interfacial instabilities in liquid films. *Nat. Commun.* **9** (1), 4477.

OLDENBURG, C.M. & LEWICKI, J.L. 2006 On leakage and seepage of co₂ from geologic storage sites into surface water. *Environ. Geol.* **50**, 691–705.

PITTS, E. 1973 The stability of pendent liquid drops. Part 1. Drops formed in a narrow gap. *J. Fluid Mech.* **59**, 753.

PRINCEN, H.M., ZIA, I.Y.Z. & MASON, S.G. 1967 Measurement of interfacial tension from the shape of a rotating drop. *J. Colloid Interface Sci.* **23** (1), 99–107.

RASCÓN, C., PARRY, A. & AARTS, D. 2017 Correction for Rascón *et al.*, geometry-induced capillary emptying. *Proc. Natl Acad. Sci. USA* **114**, E3162–E3162.

RIETZ, M., SCHEID, B., GALLAIRE, F., KOFMAN, N., KNEER, R. & ROHLFS, W. 2017 Dynamics of falling films on the outside of a vertical rotating cylinder: waves, rivulets and dripping transitions. *J. Fluid Mech.* **832**, 189–211.

ROSENTHAL, D.K. 1962 The shape and stability of a bubble at the axis of a rotating liquid. *J. Fluid Mech.* **12**, 358–366.

SHOSHO, C.E. & RYA, M.E. 2001 An experimental study of the motion of long bubbles in inclined tubes. *Chem. Engng Sci.* **56**, 2191–2204.

VAN STEIJN, V., KREUTZER, M.T. & KLEIJN, C.R. 2008 Velocity fluctuations of segmented flow in microchannels. *Chem. Engng J.* **135**, 159–165.

STONE, H.A., STROOCK, A.D. & AJDARI, A. 2004 Engineering flows in small devices: microfluidics toward a lab-on-a-chip. *Annu. Rev. Fluid Mech.* **36**, 381–411.

SUNG, J.H. & SCHULER, M.L. 2009 Prevention of air bubble formation in a microfluidic perfusion cell culture system using a microscale bubble trap. *Biomed. Microdevices* **11**, 731–738.

SVEDBERG, T. & FÄHRAEUS, R. 1926 A new method for the determination of the molecular weight of the proteins. *J. Am. Chem. Soc.* **48** (2), 430–438.

THAN, P., PREZIOSI, L., JOSEPH, D.D. & ARNEY, M. 1988 Measurement of interfacial tension between immiscible liquids with the spinning road tensiometer. *J. Colloid Interface Sci.* **124** (2), 552–559.

TORZA, S. 1975 The rotating-bubble apparatus. *Rev. Sci. Instrum.* **46**, 778–783.

TURANO, E., CURCIO, S., DE PAOLA, M.G., CALABRÒ, V. & IORIO, G. 2002 Ultra-high mobility transparent organic thin film transistors grown by an off-centre spin-coating method. *J. Membr. Sci.* **209** (2), 519–531.

VONNEGUT, B. 1942 Rotating bubble method for the determination of surface and interfacial tensions. *Rev. Sci. Instrum.* **13** (1), 6–9.
WANG, S. & CLARENS, A.F. 2012 The effects of CO₂-brine rheology on leakage processes in geologic carbon sequestration. *Water Resour. Res.* **48**, W08518.

WEBER, M.E., ALARIE, A. & RYAN, M.E. 1986 Velocities of extended bubbles in inclined tubes. *Chem. Engng Sci.* **41**, 2235-2240.

WHITE, E.T. & BEARDMORE, R.H. 1962 The velocity of rise of single cylindrical air bubbles through liquids contained in vertical tubes. *Chem. Engng Sci.* **17**, 351-361.

YU, Y.E., KHODAPARAST, S. & STONE, H.A. 2018 Separation of particles by size from a suspension using the motion of a confined bubble. *Appl. Phys. Lett.* **112**, 181604.

YU, Y.E., MAGNINI, M., ZHU, L., SHIM, S. & STONE, H.A. 2021 Non-unique bubble dynamics in a vertical capillary with an external flow. *J. Fluid Mech.* **911**, A34.

YUAN, Y., Giri, G., AYZNER, A.L., ZOOMBELT, A.P., MANNSFELD, S.C.B., CHEN, J., NORDLUND, D, TONEY, M.F., HUANG, J & BAO, Z. 2014 Ultra-high mobility transparent organic thin film transistors grown by an off-centre spin-coating method. *Nat. Commun.* **5** (1), 3005.

ZHOU, G. & PROSPERETTI, A. 2021 Faster Taylor bubbles. *J. Fluid Mech.* **920**, R2.

ZHOU, G. & PROSPERETTI, A. 2024 Volume oscillations slow down a rising Taylor bubble. *J. Fluid Mech.* **978**, A13.

ZUKOSKI, E.E. 1966 Influence of viscosity, surface tension, and inclination angle on motion of long bubbles in closed tubes. *J. Fluid Mech.* **25**, 821-837.