2.1. ‘Nose-up’ branch: bubble profile with nose pointing upwards

We begin by summarizing the results from Magnini et al. (Reference Magnini, Khodaparast, Matar, Stone and Thome2019), who investigated the dynamics of a confined bubble translating at steady state under an external upward flow. In fact, with a small variation, as we will describe in the following context, it can be extended to describe a more general case: the dynamics of a bubble with its ‘nose’ pointing upwards (figure 1b).

As indicated in figure 1(b), the coordinates are represented by $(x, y)$, with $x$ pointing vertically upwards and $y = R- r$ pointing radially inwards from the tube inner wall. With the assumption that the film thickness is much smaller than the tube radius, $b/R\ll 1$, ($x,y$) can effectively be treated as local Cartesian coordinates, and the lubrication approximation can be applied to simplify the Navier–Stokes equations as

(2.1a)\begin{gather} \text{continuity}\quad u_x +v_y= 0, \end{gather}

(2.1b)\begin{gather}x\text{-direction}\quad u_{yy} = \frac{1}{\mu} \left(p_x + \rho g \right), \end{gather}

(2.1c)\begin{gather}y\text{-direction}\quad p_y = 0, \end{gather}

where subscripts denote derivatives, the velocities in $(x, y)$ are represented by $(u, v)$ and $p$ denotes the fluid pressure.

In a reference frame moving with the bubble at the steady-state speed $U_b$, the velocity profile within the thin film can be obtained by integrating equation (2.1b) with the boundary conditions

(2.2a,b)\begin{equation} u\big\vert_{y = 0} = -U_b \quad \text{and} \quad u_y \big\vert_{y = h(x)} = 0, \end{equation}

which results in

(2.3)\begin{equation} u(y) = \frac{1}{2\mu} (p_x + \rho g ) (y^2 - 2hy ) - U_b. \end{equation}

Note that, due to the gravitational effects, fluid in the thin film is always draining vertically downwards with a parabolic profile. Specifically, the velocity profile within the uniform film region $u_b(y) = \rho g (y^2 - 2by)/(2\mu ) - U_b$ can then be obtained by demanding $h(x) = b$ and $p_x = 0$ from (2.3). An expression for $\kappa _x$, the gradient of curvature along the $x$-direction, can thus be obtained by matching the fluxes from integrating the two velocity profiles, and, consistent with the relative magnitude of terms in the lubrication approximation, expressing the gradient of capillary pressure as $p_x = -\gamma \kappa _x$:

(2.4)\begin{equation} \kappa_x = \frac{3\mu U_b}{\gamma}\frac{(h-b)}{h^3} + \frac{\rho g}{\gamma}\frac{(h^3 - b^3)}{h^3}. \end{equation}

The two terms in (2.4) indicate the contributions from the capillary and buoyancy effects, respectively. This expression for $\kappa _x$, obtained from the lubrication equations, can be equated to its geometrical counterpart, $\kappa _x = (\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {n})_x$, where $\boldsymbol {n}$ denotes the normal vector of the bubble surface pointing inwards to the gas phase (figure 1b), leading to

(2.5)\begin{equation} \frac{3\mu U_b}{\gamma}\frac{(h-b)}{h^3} + \frac{\rho g}{\gamma}\frac{(h^3 - b^3)}{h^3} = \left[ \frac{1}{(h_x^2 + 1)^{1/2}}\frac{1}{R-h} + \frac{h_{xx}}{(h_x^2 + 1)^{3/2}} \right]_x. \end{equation}

As noticed previously (Magnini et al. Reference Magnini, Khodaparast, Matar, Stone and Thome2019), the full expression for the curvature of the film profile is necessary in order for (2.4) to describe the liquid flow in both the thin film and bubble cap regions. One can further proceed with non-dimensionalization by defining $X = x/\ell$, $H = h/b$ and $K = \kappa \ell ^2/b$, with $\ell$ denoting the characteristic length scale in the $x$-direction, which will be determined below. Different from Magnini et al. (Reference Magnini, Khodaparast, Matar, Stone and Thome2019), non-dimensionalizing (2.5) shows that there are two choices for the characteristic scale $\ell$, which lead to two different expressions for $\epsilon \equiv b/\ell$:

(2.6)\begin{equation} \epsilon_1 = {Ca}_b^{1/3} \quad \text{or} \quad \epsilon_2 = ( \alpha^2 {Bo} )^{1/3}, \end{equation}

where $\alpha = b/R$. Note that $\epsilon _1$ and $\epsilon _2$ are consistent with the characteristic scales for the classic Bretherton problems in horizontal (${Bo} = 0$) and vertical orientations (${Ca}_l = 0$), respectively. Furthermore, the ratio between these two choices, ${\lambda \equiv (\epsilon _2/\epsilon _1)^3 = \alpha ^2 {Bo}/{Ca}_b}$, quantifies the relative significance of the buoyancy and capillary effects within the thin film region. While the right-hand side (RHS) of the rescaled equation (2.5) shows

(2.7)\begin{equation} \text{RHS} = \frac{H_{XXX}}{f^3} - \frac{3 \epsilon^2 H_X H_{XX}^2}{f^5} - \frac{\alpha}{1-\alpha H}\frac{H_X H_{XX}}{f^3} + \left( \frac{\alpha}{1-\alpha H}\right)^2\frac{H_X}{\epsilon^2 f}, \end{equation}

with $f = (\epsilon ^2H_X^2 + 1)^{1/2}$, different choices of $\epsilon$ result in a slightly different left-hand side (LHS) of the rescaled equation (2.5):

(2.8a)\begin{gather} \epsilon = \epsilon_1, \quad \text{LHS} = 3 \frac{(H-1)}{H^3} + \lambda \frac{(H^3-1)}{H^3}, \end{gather}

(2.8b)\begin{gather}\epsilon = \epsilon_2, \quad \text{LHS} = \frac{3}{\lambda} \frac{(H-1)}{H^3} + \frac{(H^3-1)}{H^3}. \end{gather}

When the system is provided with an upward external flow, as in Magnini et al. (Reference Magnini, Khodaparast, Matar, Stone and Thome2019), both length-scale choices are well defined and thus perform equivalently. When solving for the film profile with a downward flow, however, choosing $\epsilon _1$ becomes problematic, as ${Ca}_b\to 0$, which leads to an artificial singularity in the term associated with $\lambda$ in (2.8a).

As a result, in order to extend the film profile solution with an upward nose towards the downward flow regime, while avoiding the artificial singularity, the characteristic length scale is chosen interactively during the numerical shooting process (see § 2.3 for more detail): $\epsilon = \epsilon _1$ is chosen when $\lambda < 1$, where capillary effects outweigh buoyancy in the thin film and (2.7) and (2.8a) are solved; otherwise, when $\lambda \geqslant 1$, $\epsilon = \epsilon _2$ is chosen and (2.7) and (2.8b) are solved.

Note that, in the simulations and physical experiments, the strength of the external flow, ${Ca}_l$, is directly controlled rather than the bubble speed, ${Ca}_b$, which enters $\lambda = \alpha ^2{Bo}/{Ca}_b$. Therefore, an additional relationship is needed to link the two capillary numbers, which can be obtained by balancing the flux from the external Poiseuille flow in the far field and the cross-sectional flux in the uniform film region. Integrating the velocity profiles in cylindrical coordinates, we have

(2.9)\begin{equation} {Ca}_b = \frac{{Ca}_l}{(1-\alpha)^2} + {Bo} \left[-\frac{1}{2} + \frac{3 (1-\alpha)^2 }{8} + \frac{1}{8(1-\alpha)^2} -\frac{1}{2} (1-\alpha )^2 \log(1-\alpha) \right]. \end{equation}

With any two of ${Bo}$, ${Ca}_l$, ${Ca}_b$ and $\alpha = b/R$ given, one can solve (2.7)–(2.9) for the other quantities with numerical shooting (see § 2.3). One of the questions of interest is the critical strength of the external downward flow needed, ${Ca}_{l, cr}$, in order to stabilize the bubble, i.e. ${Ca}_b = 0$. In this scenario, the critical magnitude of external flow ${Ca}_{l, cr}$ and the corresponding film thickness $\alpha = b/R$ can be solved based on the two inputs – the Bond number ${Bo}$ and ${Ca}_b = 0$.

Typical results from a numerical solution are shown in figure 2, where the critical downward flow ${Ca}_{l, cr}$ as a function of ${Bo}$ is shown as the black solid curve in figure 2(a). For a fixed ${Bo}$, the bubble rises if the provided external flow is above the curve with ${Ca}_{l}>{Ca}_{l, cr}$, and vice versa. This solution is consistent with the results for a vertically oriented capillary in Bretherton (Reference Bretherton1961), as ${Ca}_{l, cr}$ = 0 remains true for all ${Bo} <{Bo}_{cr}$, and non-zero downward flow is needed to maintain the bubble stationary in the laboratory frame otherwise. The magnitude of ${Ca}_{l, cr}$ remains small for ${Bo} \approx 1$, and rapidly increases as gravitational effects become more dominant.

Figure 2.

Critical condition for a bubble being stabilized in the laboratory frame ($Ca_b = 0$) under an external downward flow. (a) The critical strength of the external flow ${Ca}_{l,cr}$ as a function of ${Bo}$. Consistent with results in the literature, external downward flow is needed to keep the bubble stationary in the laboratory only for ${Bo} > {Bo}_{cr}$, and, as ${Bo}$ increases, a stronger downward flow is required. (b) The uniform film thickness $b/R$ when the critical external flow condition ${Ca}_{l,cr}$ is applied, where $\alpha = b/R$ monotonically increases with ${Ca}_{l,cr}$. In the limit of $\alpha \to 0$, expanding (2.9) shows that $\alpha \sim (-{Ca}_{l,cr}/{Bo})^{1/3}$, and the asymptotic scaling is shown as the red dashed curve.

The corresponding critical film thickness is shown in figure 2(b), where the inset is plotted in log–log scales. Note that, in the limit of $\alpha \to 0$, (2.9) can be expanded. Since both the external flow and buoyancy are significant, taking the leading order of each term and demanding ${Ca}_b = 0$ (a stationary bubble) leads to

(2.10)\begin{equation} \alpha =\left(-\frac{3}{2}\frac{{Ca}_{l,cr}}{{Bo}}\right)^{1/3}. \end{equation}

This asymptotic approximation is shown in figure 2(b) as the red dashed curve, which is in excellent agreement with the numerical shooting solutions for $\alpha = b/R < 0.1$. Furthermore, the rapid increase of ${Ca}_{l, cr}$ with ${Bo}$ can thus be qualitatively explained, since the critical film thickness increases with ${Bo}$, and the scaling shows $|{Ca}_{l, cr}| \sim \alpha ^3{Bo}$.

2.2. ‘Nose-down’ branch: bubble with nose pointing downwards

Equations (2.7)–(2.9) fully describe the dynamics of the bubble under an upward external flow (${Ca}_{l} \geq 0$). For ${Bo} > {Bo}_{cr}$, the solution can be further extended towards the downward flow regime (${Ca}_{l} < 0$), with the bubble nose pointing upwards and ${Ca}_b \to {Ca}_l^+$. However, this solution is only valid for a limited range of ${Ca}_{l} < 0$, beyond which solutions of (2.7)–(2.9) fail to converge within the set tolerance, and thus the branch terminates. In a physical system, the bubble responds to the stronger downward flow by changing its morphology, so that the ‘nose’ points downwards as indicated in figure 1(c), following the flow direction. Both experiments and simulations introduced in the next sections will confirm this response.

The bubble profile in a downward-nose system is similar to the upward-nose case, and therefore, when it comes to solving for the film profile, the new film equations can be obtained by modifying (2.7)–(2.9) simply by redefining the vertical coordinate as $\tilde {x} = -x$ (figure 1c). With this convention, the positive directions of ${Ca}_l$ and ${Ca}_b$ within the model are redefined, while the scalar variables (e.g. $\alpha$ and $H$) remain unchanged. Furthermore, since the gravitational direction remains vertically downwards regardless of coordinates, in the new coordinate system, the terms associated with ${Bo}$ must change sign as well. As a result, the film equations become

(2.11a)\begin{gather} 3 \frac{(H-1)}{H^3} - \lambda \frac{(H^3-1)}{H^3} = \frac{H_{XXX}}{f^3} - \frac{3 \epsilon^2 H_X H_{XX}^2}{f^5} - \frac{\alpha}{1-\alpha H}\frac{H_X H_{XX}}{f^3} + \left( \frac{\alpha}{1-\alpha H}\right)^2\frac{H_X}{\epsilon^2 f}, \end{gather}

(2.11b)\begin{gather}{Ca}_b = \frac{{Ca}_l}{(1-\alpha)^2} - {Bo} \left[-\frac{1}{2} + \frac{3 (1-\alpha)^2 }{8} + \frac{1}{8(1-\alpha)^2} -\frac{1}{2} (1-\alpha )^2 \log(1-\alpha) \right], \end{gather}

where, compared to (2.7)–(2.9), the net effects are sign changes in the terms related to gravity only. In other words, by changing the nose direction from upwards to downwards, the bubble only senses a sign change in the buoyancy force: instead of assisting the bubble to rise, buoyancy effects are now acting as resistance to the bubble motion, which is now mainly driven by the downward external flow.

When the bubble dynamics falls on this solution branch, the bubble always sinks and the film thickness $\alpha$ converges to the origin as ${Ca}_l \to 0^-$ (as we will confirm in § 4). Thus, $\epsilon = \epsilon _1$ is chosen when non-dimensionalizing the equations, resulting in the form of (2.11a) similar to (2.8a).

2.3. Numerical integration

With ${Ca}_l$ and ${Bo}$ given, (2.7)–(2.9) or (2.11) can be solved for ${Ca}_b$, the uniform film thickness $\alpha$ and the film profile $H(X)$, with the boundary conditions

(2.12a–c)\begin{equation} H(0) = 1, \quad H_X(0) = 0 \quad \text{and} \quad H_{XX}(0) = 0, \end{equation}

and $H(X_{nose}) = 1/\alpha$ is demanded as the additional constraint, with $X = 0$ indicating the uniform film region and $X_{nose}$ the location of the bubble nose, which is yet to be determined by the numerical integration scheme. Although the full curvature of the film is used in (2.7)–(2.9) or (2.11), the underlying solution method is effectively the same as the asymptotic expansion as seen in Bretherton (Reference Bretherton1961), which includes the process of solving the film equation for the bubble ‘nose’ and matching to the static spherical cap. The coupled equations are solved by first imposing an initial guess on $\alpha$. With ${Ca}_l$ and ${Bo}$ given, (2.9) or (2.11b) outputs the capillary number of the bubble ${Ca}_b$, which serves as an input for (2.7) and (2.8) or (2.11a). The ordinary differential equation for $H(X)$ is solved by numerical integration with the MATLAB program $ode45$ from the uniform film region towards the front spherical cap in the positive $X$-direction. Numerical integration ceases when $H_X \to \infty$, where the termination location denotes $X_{nose}$ and the value $H(X_{nose})$ is compared with the constraint $H(X_{nose}) = 1/\alpha$. The initial guess of $\alpha$ is then iteratively updated until this additional constraint is met, meaning that the bubble nose is symmetric about the tube centreline.

3.1. Experimental set-up

Experiments are performed in a refractive-index-matching set-up similar to those described in Yu et al. (Reference Yu, Zhu, Shim, Eggers and Stone2018) and Magnini et al. (Reference Magnini, Khodaparast, Matar, Stone and Thome2019), with pure glycerol filling the capillary tube as the continuous phase. The density, viscosity and surface tension of the pure glycerol are measured as $\rho = (1.29 \pm 0.001) \times 10^3\ \textrm {kg}\,\textrm {m}^{-3}$, $\mu = 1.00 \pm 0.04\ \textrm {Pa}\,\textrm {s}$ (Anton Parr, Physica MCG 301) and $\gamma = 65.4 \pm 1.0\ \textrm {mN}\,\textrm {m}^{-1}$ (pendant drop), respectively. Glass capillaries with three different inner diameters ($ID = 4.45, 4.85, 5.65 \pm 0.05\ \textrm {mm}$) are used, yielding the Bond numbers of ${Bo} = 0.97, 1.16, 1.56$, respectively. Each glass capillary is placed within a cuboid glass box, which is also filled with pure glycerol in order to avoid imaging distortions. The top end of the glass capillary is connected to a liquid reservoir by a Teflon hose, and the bottom end is submerged in a bath of pure glycerol. The experimental flow rate is adjusted by controlling the reservoir pressure using the Elveflow$^{\circledR }$ OB1 MK3 pressure and vacuum controller, and calibration between the flow rate and controller pressure is performed for all capillary tubes used. The imaging apparatus is composed of a Nikon D5100 DSLR camera and a Mitutoyo infinity-corrected objective, mounted on an in-house-made tube microscope. The imaging apparatus is aligned horizontally with a collimated LED source, which is located half-way between the two ends of the capillary tube.

In each set of experiments with the same ${Bo}$, the capillary tube is partially filled with pure glycerol, leaving a section of an air column in the bottom end of the capillary. The set-up is then carefully calibrated to be vertical. A single bubble is formed by applying vacuum to the system, providing an upward flow and assisting the formation of the thin liquid film. For experiments corresponding to ${Ca}_b > 0$, experiments start by directly setting the pressure/vacuum to a target value; for experiments with ${Ca}_b < 0$, on the other hand, the target pressure value is set after the bubble rises to the top end of the capillary. As the bubble reaches the region of interest of the imaging apparatus, the dynamics of the bubble are recorded in the bright-field mode at 30 frames per second. After each experiment, the bubble velocity and the film profile are analysed from the image sequence, and the capillary number of the external flow ${Ca}_l$ is calculated based on the calibration between flow rate and pressure control.

3.2. Numerical simulation set-up

Direct numerical simulations of the flow of elongated bubbles in vertical capillaries are performed utilizing the volume-of-fluid (VOF) method (Hirt & Nichols Reference Hirt and Nichols1981) implemented in OpenFOAM. The unsteady mass and momentum equations for an incompressible flow and Newtonian fluid are solved, together with a transport equation for a passive scalar that identifies the gas and liquid phases across the domain. In this formulation, the surface tension force is implemented as a body force according to the continuum surface force method (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992). Both two-dimensional axisymmetric and full three-dimensional simulations are conducted, with the latter enabling us to investigate conditions that may lead to symmetry breaking. The simulation set-up is similar to that adopted in previous works (Magnini et al. Reference Magnini, Ferrari, Thome and Stone2017, Reference Magnini, Khodaparast, Matar, Stone and Thome2019). An elongated bubble, with the shape of a cylinder with spherical ends, is initialized at one end of the flow domain. Bubble lengths of approximately $6D$ are sufficient to achieve a uniform film region between front and rear menisci, under the conditions of interest. At the inlet boundary, a fully developed laminar profile is set for the incoming liquid, while no slip is set at the pipe wall. At the channel outlet, pressure is given a constant value, together with a zero-gradient condition for the velocity. The liquid-to-gas density and viscosity ratios are set to 1000 and 100, respectively. Each simulation is run forwards in time until the bubble translates with a constant speed.

Numerical simulations are performed for the three values of the Bond number tested experimentally and a wide range of ${Ca}_l$, for both upward and downward liquid flow. For each value of ${Bo}$ and each flow direction, a first simulation is run with the largest $|{Ca}_l|$ desired, until steady state. From this steady solution, $|{Ca}_l|$ is reduced and a new simulation is run until a new steady bubble profile and speed are achieved. The procedure continues by stepping towards smaller values of $|{Ca}_l|$ to span the entire range of conditions of interest, each time until steady state.