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 Reference Groell, Schaffler and Rienmueller1997) to CO$_2$ sequestration in depleted geological oil reservoirs (Corapcioglu, Cihan & Drazenovic Reference Corapcioglu, Cihan and Drazenovic2004; Oldenburg & Lewicki Reference Oldenburg and Lewicki2006; Wang & Clarens Reference Wang and Clarens2012), embolism in circulatory biological systems (Brodribb, Bienaimé & Marmottant Reference Brodribb, Bienaimé and Marmottant2016; Li et al. Reference Li, Li, Li, Chen, Li and Chen2021) and multiphase microfluidics flows (Baroud, Galllaire & Dangla Reference Baroud, Galllaire and Dangla2010). In miniaturized fluid systems, air bubbles can be exploited for transport of particles or for mixing processes (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Baroud, de Saint Vincent & Delville Reference Baroud, de Saint Vincent and Delville2007; Baroud et al. Reference Baroud, Galllaire and Dangla2010). 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 Reference Jensen, Goranović and Bruus2004; van Steijn, Kreutzer & Kleijn Reference van Steijn, Kreutzer and Kleijn2008; Brodribb et al. Reference Brodribb, Bienaimé and Marmottant2016).
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 Reference Sung and Schuler2009). As an additional example, in fuel cells, the oxidation of methanol leads to the formation of CO$_2$ bubbles that reduce a cell's efficiency (Litterst et al. Reference Litterst, Eccarius, Hebling, Zengerle and Koltay2006). Thus, a considerable effort has been dedicated to the removal of bubbles in these circuits (see e.g. Sung & Schuler Reference Sung and Schuler2009; Cheng & Lu Reference Cheng and Lu2014; Guo, Liu & Ran Reference Guo, Liu, Ran, Wang, Shan, Li, Liu and Li2022).
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 Reference Fairbrother and Stubbs1935), tuning the velocity of the bubble may be used to separate particles based on their size (Yu, Khodaparast & Stone Reference Yu, Khodaparast and Stone2018): 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. Reference Kaigala, Lovchik, Drechsler and Delamarche2011), 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 (Reference Dumitrescu1943) and Davies & Taylor (Reference Davies and Taylor1950), 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 Reference Gibson1913) 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 (Reference Bretherton1961) 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 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 (Reference Bretherton1961) 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 (Reference Lamstaes and Eggers2017) 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 Reference Dhaouadi and Kolinski2019).
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 (Reference Zhou and Prosperetti2021) 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 (Reference Bi and Zhao2001) and Bico & Quéré (Reference Bico and 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. Reference Funada, Joseph, Maehara and Yamashita2005). 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é Reference Bico, Tordeux and Quéré2001).
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. Reference Yu, Magnini, Zhu, Shim and Stone2021). In particular, Magnini et al. (Reference Magnini, Khodaparast, Matar, Stone and Thome2019) 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 (Reference Kubie2000) 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 (Reference Brannock and Kubie1996) reported instead experimental evidences of the slowing down of the bubble, while Madani, Caballina & Souhar (Reference Madani, Caballina and Souhar2009, Reference Madani, Caballina and Souhar2012) 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 (Reference Zhou and Prosperetti2024) 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 Reference Svedberg and Fåhraeus1926) to numerous industrial processes, such as wastewater treatment (Turano et al. Reference Turano, Curcio, De Paola, Calabrò and Iorio2002) or crude oil refining (Gary et al. Reference Gary, Handwerk, Kaiser and Geddes2007). In interfacial flows, spinning rods (Than et al. Reference Than, Preziosi, Josephl and Arney1988) and spin-coating (Emslie, Bonner & Peck Reference Emslie, Bonner and Peck1958) are used to deposit uniform thin films onto diverse substrates such as optical lenses for anti-reflective properties (Krogman, Druffel & Sunkara Reference Krogman, Druffel and Sunkara2005), or silicon wafers for organic semiconductor fabrication (Yuan et al. Reference Yuan, Giri, Ayzner, Zoombelt, Mannsfeld, Chen, Nordlund, Toney, Huang and Bao2014). 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 Reference Marthelot, Strong and Brun2018; Jambon-Puillet, Royer Piéchaud & Brun Reference Jambon-Puillet, Royer Piéchaud and Brun2021), where centrifugal instabilities (Rietz et al. Reference Rietz, Scheid, Gallaire, Kofman, Kneer and Rohlfs2017) 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 Reference Vonnegut1942; Rosenthal Reference Rosenthal1962; Princen, Zia & Mason Reference Princen, Zia and Mason1967; Torza Reference Torza1975). These experiments can measure very low interfacial tensions (Drelich, Fang & White Reference Drelich, Fang and White2002) 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 (Reference Manning, Collicott and Finn2011) studied the case of a tube of inner radius $R$ partially filled with a liquid of density $\rho$ and surface tension $\gamma$, rotated around its symmetry axis at angular velocity $\omega$ 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 $\phi$, and computed a critical angular velocity $\omega _0$:
such that the tube cannot occlude if $\omega > \omega _0$. In the case of an occluding meniscus forming a bubble, the gas–liquid interface meets the solid wall tangentially, i.e. $\phi =0$. Thus, (2.1) indicates that under weightlessness, the bubble cannot occlude the channel if $\rho \omega ^2 R^3/\gamma > 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 $\mu$, density $\rho$ and surface tension $\gamma$, both contained in a vertically oriented circular tube of radius $R \ll 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 $\ell _c=\sqrt {\gamma /\rho g}$, where $g$ is the acceleration due to gravity, so that the Reynolds number $Re=\rho U_b R/\mu$ 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 \approx 0.918 \ell _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 \ell _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 $\omega$. 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 $\omega$ is sufficiently high. In the subsequent section, we revisit Bretherton's theory (Bretherton Reference Bretherton1961) 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 Reference Bretherton1961; Lamstaes & Eggers Reference Lamstaes and Eggers2017) 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
By integrating the radial component of (2.2), we obtain
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
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
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 ${\textrm {d}r_1}/{\textrm {d}s} = \cos (\theta )$ and ${\textrm {d}z}/{\textrm {d}s} = -\sin (\theta )$, differentiating with respect to the curvilinear coordinate $s$ gives
Finally, with the dimensionless variables $\bar {r}_1=r_1/R$ and $\bar {s}=s/R$, (2.6) becomes
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 (Reference Bretherton1961), 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)
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
This flux must equate the volume of fluid displaced per unit time by the top of the bubble, which is equal to ${\rm \pi} R^2 U_b$. Since $y_1/R \ll 1$, the $-2{\rm \pi} RU_by_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:
Since $b$ is the length scale governing the flow in the inner region, we non-dimensionalize, as in Bretherton (Reference Bretherton1961), with
This leads to the following ordinary differential equation:
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 \rightarrow \infty$ in (2.12) for the inner solution. In this limit, the equation behaves as
whose general solution reads
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 \rightarrow 1$, (2.12) becomes
for which the only non-oscillating solution is $\eta _{0}=1 + \mathcal {C}\exp (\mathcal {F}\zeta )$, where $\mathcal {C}$ is an integration constant and
Since the value of $\mathcal {C}$ can be adjusted by shifting the origin of $\zeta$, we can set $\mathcal {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.
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 \,\textrm {e}^{\sqrt {a}\zeta }+c_2\,\textrm {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
Using the shifted variable $\chi =\zeta -\zeta ^*$, we can now define $\eta (\chi )=c_1^*\,\textrm {e}^{\sqrt {a}\chi }+c_2^*\,\textrm {e}^{-\sqrt {a}\chi }+c_3^*-{\chi }/{a}$, where the new coefficients $c_1^*$, $c_2^*$ and $c_3^*$ are expressed as
The new coefficients are well fitted (see figure 3b) by
Thus, the distance from the wall at which the outer profile exhibits an inflection point can be evaluated as
Furthermore, the slope of the profile at the inflection point is given by
Remarkably, the conditions on film thickness (2.21) and slope (2.22) at the inflection point are the same as described in Bretherton (Reference Bretherton1961). Therefore, the centrifugal force alters 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 (Reference Bretherton1961).
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)\vert _{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)\vert _{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)\vert _{r=R})=\tan ^{-1}(y_1'(0))$. A closer inspection reveals that within a small range around $Ce_c$, i.e. for $\vert Ce-Ce_c(Bo)\vert <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
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. (Reference Manning, Collicott and Finn2011) 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
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
where $Bo_{c,0}=0.842$ is the critical Bond number in the absence of centrifugation ($Ce=0$). Note that in the limit $Bo \rightarrow 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. (Reference Manning, Collicott and Finn2011) under weightlessness. Conversely, in the limit $(Ce_c \rightarrow 0$, $Bo\rightarrow Bo_{c,0})$, (2.25) simplifies into
By injecting (2.26) into (2.23), we obtain, in the limit $Ce \rightarrow 0$:
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):
Since the volume of fluid displaced per unit time by the tip of the bubble ${\rm \pi} R^2 U_b$ should be equal to 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 g b^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$:
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$ kg m$^{-3}$, $\mu =9.64\times 10^{-2}$ Pa s, $\gamma =2.09\times 10^{-2}\ {\rm N}\ {\rm 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 Reference Vonnegut1942). 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; 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 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_{min}=1.6 \times 10^{-7}$. A value inferior to this limit will be accordingly set equal to zero. The maximal rotational speed achieved by the DC motor is $\omega _{max}=400$ rad s$^{-1}$. For a given tube inner radius, this sets a limit on the maximal centrifugal number $Ce_{max}$ 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\equiv 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 operating the motor. As observed by Kubie (Reference Kubie2000), 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 (Reference Kubie2000) 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 (Reference White and Beardmore1962), 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 Reference Zukoski1966; Maneri & Zuber Reference Maneri and Zuber1974; Bendiksen Reference Bendiksen1984; Weber, Alarie & Ryan Reference Weber, Alarie and Ryan1986; Couët & Strumolo Reference Couët and Strumolo1987; Shosho & Rya Reference Shosho and Rya2001; Boucher et al. Reference Boucher, Karp, Belt and Liné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 Reference Boycott1920; Acrivos & Herbolzheimer Reference Acrivos and Herbolzheimer1979) in the case of settling suspensions in sealed tubes, as seminally observed by Boycott (Reference Boycott1920) 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 (Reference Zukoski1966) 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^\circ$ to $80^\circ$ with the horizontal, with a maximum 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 (Reference Collicott and Manning2014) 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 $\boldsymbol {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. Reference Manning, Collicott and Finn2011; Rascón, Parry & Aarts Reference Rascón, Parry and Aarts2017):
where $\bar {h}(\bar {x},\bar {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: $({\boldsymbol {\nabla } h}/{\sqrt {1+(\boldsymbol {\nabla } h)^2}})\boldsymbol {\cdot}\boldsymbol {n}=-\cos (\phi )$, where $\boldsymbol {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 $\boldsymbol {e}_z$ and the tangent at the wall to the intersection of the liquid–gas interface with the plane $(\boldsymbol {n}, \boldsymbol {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 $\boldsymbol {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 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 $\boldsymbol {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_bt,\tilde {y})$, where $\boldsymbol {e}_{\tilde {x}}$ is aligned with the tube central axis and points upwards (such that $\boldsymbol {e}_{\tilde {x}}\boldsymbol {\cdot}\boldsymbol {g}=-g\sin (\alpha )$), and $\boldsymbol {e}_{\tilde {y}}$ is the inward vector normal to the inner solid wall, such that $\boldsymbol {e}_{\tilde {y}}\boldsymbol {\cdot}\boldsymbol {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)
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
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 ${\rm \pi} R^2 U_b$. Given that $y_1/R \ll 1$, the term $-2{\rm \pi} RU_by_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:
We non-dimensionalize with
which leads to the following ordinary differential equation:
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
and the slope of the inner region profile at the inflection point is given by
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 \rightarrow 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 0.05$R$ to 0.01$R$ (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.
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
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
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 ($Bo_c=0.842$, $\beta =0.49$ 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
The volume of fluid displaced per unit time by the tip of the bubble ${\rm \pi} R^2U_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:
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$ and reaches its maximum at $\alpha \approx 50^\circ$, a value independent of $Bo$ within the range of Bond numbers investigated here, as shown in figure 12(b). Furthermore, for $Bo_c(\alpha =50^\circ )=0.5401< Bo<0.842=Bo_c(\alpha =90^\circ )$, 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(\alpha \approx 50^\circ )$. Those observations align well with our theoretical analysis.
Finally, the bubble velocity as a function of the tilt angle $\alpha$ 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. (Reference Atasi, Khodaparast, Scheid and Stone2017) 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 Reference Zukoski1966). 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}(\alpha )$ 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 $\alpha _{opt} \lessapprox 50^\circ$. 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 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 $\alpha \rightarrow 0^\circ$ (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 (Reference Bretherton1961) 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 Reference Bretherton1961) 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:
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 (Reference Lamstaes and Eggers2017) and experimentally investigated by Dhaouadi & Kolinski (Reference Dhaouadi and Kolinski2019) 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.
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 \boldsymbol {e}_z$, the stationary Navier–Stokes equations are written as
where $\boldsymbol {\varOmega }=\omega \boldsymbol {e}_z$ is the rotation vector.
We introduce the dimensionless variables $\bar {\boldsymbol {u}}$, $\bar {z}$, $\bar {r}$ and $\bar {p}$, such that $\boldsymbol {u}=U_b\bar {\boldsymbol {u}}$, $z=R\bar {z}$, $r=R\bar {r}$ and $p={\rm \pi} \bar {p}$, where ${\rm \pi} = \rho \omega R U_b$. The dimensionless Navier–Stokes equations read
where $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, $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
Furthermore, in cylindrical coordinates, the stress tensor is written as
and the vector normal to the interface is $\boldsymbol {n}=({1}/{\sqrt {1+r'_1(z)^2}})(-1, 0, r'_1(z))^\textrm {T}$.
Therefore, the dynamic and kinematic boundary conditions at the fluid–air interface are
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 Ry + \rho g z$. In this system of coordinates, the Navier–Stokes equation become
Likewise, the dynamic and kinematic boundary conditions become
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 \overline {u_z}$, $u_y=U_y \overline {u_y}$, $u_x=U_x \overline {u_x}$, $P=P_0 \bar {P}$, $y=b\bar {y}$, $y_1=b\overline {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 \overline {u_y}/\partial \bar {y}+\partial \overline {u_z}/\partial \bar {z}=0$.
Furthermore, upon introduction of the dimensionless fields and variables, the momentum conservation equations are written as
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
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 ^2R}) = O(\epsilon ^3)$. Thus, the system of equations reduces to
The dynamic and kinematic boundary conditions are
Thus, including the no-slip boundary condition at the solid wall, the boundary conditions at leading order are
Finally, going back to the dimensional form, and reintroducing the original pressure field $p=P-\rho \omega ^2 Ry - \rho g z$, the full problem is written as
and is complemented by the following boundary conditions:
The integration of the pressure field is straightforward and leads to
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:
which results in
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 $U_b \boldsymbol {e}_z$, where $z$ is aligned with the central axis of the (tilted) tube (see figure 14a), the Navier–Stokes equations read
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\,\textrm {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
These equations are complemented by the following dynamic and kinematic boundary conditions:
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
Likewise, the dynamic and kinematic boundary conditions become
where $\kappa =-{1}/{((R-y_1(z))\sqrt {1+y_1'(z)^2})}-{y_1''(z)}/{(1+y_1'(z)^2)^{3/2}}$.
A.2.2. Lubrication approximation
We non-dimensionalize as follows: $u_z=U_b \overline {u_z}$, $u_y=U_y \overline {u_y}$, $P=P_0 \bar {P}$, $y=b\bar {y}$, $y_1=b\overline {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 \overline {u_y}/\partial \bar {y}+\partial \overline {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
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
The dynamic and kinematic boundary conditions are
Therefore at leading order, and including the no-slip boundary condition at the solid wall, the boundary conditions are written as
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
and is complemented by the following boundary conditions:
The pressure field integrates straightforwardly into
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:
which leads to the velocity profile:
Appendix B. Derivation of the equilibrium equation for the static cap
Following previous works (see e.g. Manning et al. Reference Manning, Collicott and Finn2011; Lubbers et al. Reference Lubbers, Weijs, Botto, Das, Andreotti and Snoeijer2014; Rascón et al. Reference Rascón, Parry and Aarts2017), 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 $\boldsymbol {e}_x\boldsymbol {\cdot}\boldsymbol {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
where the first term on the right-hand side is the surface energy. As in Rascón et al. (Reference Rascón, Parry and Aarts2017), the surface $\mathcal {A}$ is computed as
with $\varOmega$ the cross-section of the capillary. The term $\mathcal {G}$ represents in turn the gravitational potential energy, which in general form reads (Pitts Reference Pitts1973)
where the volume $V$ is given by
Upon introduction of the Lagrange multiplier $\lambda$ to ensure volume conservation, the functional to be minimized to obtain equilibrium reads
Upon integration along the $z$ direction between 0 and $h$:
Formal minimization of the functional $F(h)$ with respect to $h$ leads to the following partial differential equation:
with the constant slope condition at the wall:
where $\boldsymbol {n}$ is the unit exterior normal to the tube wall. The value of $\lambda$ can be set by integrating the resulting equilibrium equation within the whole domain, leading to the following expression:
where ${h}_0=V/\varOmega$ 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 Jurin height. Upon non-dimensionalization with the tube radius $R$, one obtains the equilibrium equation (3.1) reported in the main text.