Skip to main content Accessibility help
×
Home

Information:

  • Access
  • Open access

Figures:

Actions:

MathJax
MathJax is a JavaScript display engine for mathematics. For more information see http://www.mathjax.org.
      • Send article to Kindle

        To send this article to your Kindle, first ensure no-reply@cambridge.org is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the @free.kindle.com or @kindle.com variations. ‘@free.kindle.com’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘@kindle.com’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        The shape and motion of gas bubbles in a liquid flowing through a thin annulus
        Available formats
        ×

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        The shape and motion of gas bubbles in a liquid flowing through a thin annulus
        Available formats
        ×

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        The shape and motion of gas bubbles in a liquid flowing through a thin annulus
        Available formats
        ×
Export citation

Abstract

We study the shape and motion of gas bubbles in a liquid flowing through a horizontal or slightly inclined thin annulus. Experimental data show that in the horizontal annulus, bubbles develop a unique ‘tadpole-like’ shape with a semi-circular cap and a highly stretched tail. As the annulus is inclined, the bubble tail tends to vanish, resulting in a significant decrease of bubble length. To model the bubble evolution, the thin annulus is conceptualised as a ‘Hele-Shaw’ cell in a curvilinear space. The three-dimensional flow within the cell is represented by a gap-averaged, two-dimensional model, which achieved a close match to the experimental data. The numerical model is further used to investigate the effects of gap thickness and pipe diameter on the bubble behaviour. The mechanism for the semi-circular cap formation is interpreted based on an analogous irrotational flow field around a circular cylinder, based on which a theoretical solution to the bubble velocity is derived. The bubble motion and cap geometry is mainly controlled by the gravitational component perpendicular to the flow direction. The bubble elongation in the horizontal annulus is caused by the buoyancy that moves the bubble to the top of the annulus. However, as the annulus is inclined, the gravitational component parallel to the flow direction becomes important, causing bubble separation at the tail and reduction in bubble length.

1 Introduction

Simultaneous flow of gas and liquid through a conduit, over a wide range of flow rates, exhibits a ‘slug flow’ pattern, which consists of the pseudo-periodic appearance of large bubbles (often called Taylor bubbles) separated by liquid slugs (Fabre & Liné 1992). The shape and motion of Taylor bubbles in tubes have been extensively investigated over the past decades (Davies & Taylor 1950; White & Beardmore 1962; Zukoski 1966; Collins 1967b ; Collins et al. 1978; Reinelt 1987b ; Vanden-Broeck 1984, 1992; Viana et al. 2003; Figueroa-Espinoza & Fabre 2011; Fabre 2016). In contrast, much less effort has been devoted to studying Taylor bubbles in annuli, which are, however, very relevant for many industrial applications such as energy dissipation facilities for chemical/nuclear reactors, casing-tubing annuli for oil production and double-pipe heat exchangers for geothermal systems (Kelessidis & Dukler 1990; Das et al. 1998; Ekberg et al. 1999).

In an annulus, flow is confined between the inner surface of an outer pipe and the outer surface of an inner pipe. If the thickness of the gap between the two surfaces is much smaller than the perimeter of the pipe, the annulus is then similar to a Hele-Shaw cell. The problem of two-phase flow in a Hele-Shaw cell between two parallel flat plates has been the focus of many previous theoretical, experimental and numerical studies (e.g. Taylor & Saffman 1959; Collins 1965a , 1967a ; Maxworthy 1986; Kopf-Sill & Homsy 1988; Roig et al. 2012; Cueto-Felgueroso & Juanes 2014). Large bubbles in a Hele-Shaw cell (the bubble characteristic length, $d$ , is much greater than the gap thickness between the plates, $h$ ) behave very differently from bubbles in an unbounded liquid due to the presence of wall confinement and wettability effects (Thompson 1968; Eck & Siekmann 1978; Saffman & Tanveer 1989; Cueto-Felgueroso & Juanes 2014). The bubble dynamics in a Hele-Shaw cell can be characterised by the bubble Reynolds number $Re=\unicode[STIX]{x1D70C}_{l}u_{\infty }d/\unicode[STIX]{x1D707}_{l}$ ( $u_{\infty }$ is the bubble velocity relative to the ambient liquid having a density $\unicode[STIX]{x1D70C}_{l}$ and a dynamic viscosity $\unicode[STIX]{x1D707}_{l}$ ). In addition, the gap Reynolds number $Re(h/d)^{2}$ is often used to compare the in-plane inertial effect and the cross-gap viscous effect (Thompson 1968; Eck & Siekmann 1978; Bush & Eames 1998; Roig et al. 2012). When $Re(h/d)^{2}\ll 1$ , the flow corresponds to the classical Hele-Shaw flow regime where the inertia is negligible (Saffman & Taylor 1958; Taylor & Saffman 1959; Eck & Siekmann 1978; Maxworthy 1986; Tanveer 1986, 1987; Tanveer & Saffman 1987; Kopf-Sill & Homsy 1988; Meiburg 1989; Saffman & Tanveer 1989; Maruvada & Park 1996; Cueto-Felgueroso & Juanes 2014). When $Re(h/d)^{2}$ is close to or larger than 1, the inertia becomes important, as has been observed for bubbles in inclined or vertical Hele-Shaw cells (Collins 1965a , b ; Grace & Harrison 1967; Lazarek & Littman 1974; Maneri & Zuber 1974; Hills 1975; Vanden-Broeck 1984, 1992; Couet & Strumolo 1987; Bush 1997; Kelley & Wu 1997; Bush & Eames 1998; Huisman, Ern & Roig 2012; Roig et al. 2012; Wang et al. 2014, 2016; Filella, Ern & Roig 2015; Piedra, Ramos & Herrera 2015).

In this paper, we study bubble behaviour in another type of Hele-Shaw cell, i.e. thin annuli. We provide an insight, based on experimental observations, numerical simulations and theoretical analysis, into simultaneous gas–liquid flow through such systems. We focus on the regime where inertia becomes important, i.e. $Re(h/d)^{2}$ is close to or larger than unity. Laboratory experiments are conducted to measure the bubble distribution and evolution in a horizontal or slightly inclined thin annulus (§ 2). A gap-averaged, two-dimensional (2-D) formulation of the Navier–Stokes equations is then derived to represent the three-dimensional (3-D) flow dynamics (§ 3). The numerical model predictions are compared with the experimental data for different inclinations. Further numerical simulations are then performed to elucidate the roles of gap thickness and pipe diameter in bubble formation (§ 4). In addition, a theoretical analysis is presented to interpret the flow mechanisms that govern bubble dynamics in thin annuli (§ 5). Finally, a brief discussion is given and conclusions are drawn (§ 6).

2 Experimental apparatus and procedure

The experimental apparatus consisted of air and water supply systems, an annulus formed between two concentric pipes (comprising three connected subsections and initially filled with water), and a collection tank (figure 1). The parameters of the experimental apparatus and working fluids are as given in table 1. During the experiments, air and dyed water were continuously pumped into a mixer before entering the test section. Visual observation showed the two fluids were well segregated by the time they entered the annulus. Three high-speed cameras were mounted above the annulus subsections for visualisation and recording purposes. The outer pipe was made of Plexiglas to facilitate visualisation, while the inner stainless steel tubing was spray coated to create a white surface which maximised the colour contrast between the air and dyed water. The videos were taken at a frame rate of 17 frames per second with a resolution of $1280\times 720$ pixels, and a colour depth of 24 bits per pixel. The data were imported into MATLAB for post-processing to extract bubble size and velocity data.

Figure 1. Schematic illustration of the experimental apparatus.

3 Mathematical formulation and numerical model

3.1 Governing equations

Figure 2. Schematic of (a) the thin annulus and (b) the equivalent Hele-Shaw cell. (c) An illustration showing how the direction of gravitational acceleration changes along the equivalent Hele-Shaw cell to ensure the flow replicates that seen in the actual annulus.

Table 1. Parameters of the experimental apparatus and fluids.

As shown in figure 2(a), we define a 3-D curvilinear coordinate system ( $x,y,z$ ) based on the mid-plane of the thin annulus (i.e. the distance of this mid-plane from both the outer and inner walls is $h/2$ ), where $x$ is defined to be along the pipe length ( $0\leqslant x\leqslant L$ ), $y$ is around the perimeter ( $-W/2\leqslant y\leqslant W/2$ with $y=0$ and $y=\pm W/2$ corresponding to the top and bottom of the pipe, respectively), and $z$ is the direction across the gap ( $-h/2\leqslant z\leqslant h/2$ with $z=h/2$ and $z=-h/2$ defining the outer and inner channel walls, respectively). The system is therefore equivalent to a Hele-Shaw cell (viz. figure 2 b) under a distorted gravity field (viz. figure 2 c).

We use the approximation that the flow velocity in the $z$ direction is negligible ( $u_{z}=0$ ) given the very high aspect ratio (i.e. $h/W\ll 1$ ). By eliminating the terms involving $u_{z}$ , the governing equations of the incompressible two-phase, air–water flow are obtained, including the mass and momentum conservation equations given respectively by

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}z}\right)+\unicode[STIX]{x1D70C}\boldsymbol{g}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

and an equation for the volume fraction of air given as

(3.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}=0. & & \displaystyle\end{eqnarray}$$

Here, $\boldsymbol{u}=[u_{x},u_{y}]$ is the velocity vector, $\unicode[STIX]{x1D735}=[\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x,\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y]$ is the 2-D gradient, $t$ is the time, $\unicode[STIX]{x1D6FC}$ is the volume fraction of air, the bulk density is $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{a}+(1-\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D70C}_{w}$ , $p$ is the pressure, the bulk dynamic viscosity is $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D707}_{a}+(1-\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D707}_{w}$ , the gravitational acceleration vector is $\boldsymbol{g}=[g\sin \unicode[STIX]{x1D6FD},g\cos \unicode[STIX]{x1D6FD}\sin (2y/d)]$ , $\unicode[STIX]{x1D70E}$ is the surface tension, $\unicode[STIX]{x1D705}$ is the interface curvature, $\unicode[STIX]{x1D6FF}$ is the Dirac delta function and $\boldsymbol{n}$ is the interface outward-pointing unit normal. Surface tension force is treated as a continuous surface force (Brackbill, Kothe & Zemach 1992). Note (3.1)–(3.3) are written in a form following Bush (1997) that eases the derivation of gap-averaged, 2-D equations in the following § 3.2.

3.2 Gap-averaging process

We assume the velocity profile in the annulus to be parabolic in the $z$ -direction (Gondret & Rabaud 1997):

(3.4) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=\frac{3}{2}\left[1-\left(\frac{2z}{h}\right)^{2}\right]\bar{\boldsymbol{u}}, & & \displaystyle\end{eqnarray}$$

where $\bar{\boldsymbol{u}}$ is the gap-averaged velocity. This is consistent with our calculation $\unicode[STIX]{x1D70C}_{w}v_{i}(2h)/\unicode[STIX]{x1D707}_{w}=1134$ , which reveals the flow to be laminar (Beavers, Sparrow & Magnuson 1970).

We derive a coupled set of gap-averaged equations by integrating equations (3.1)–(3.3) along the $z$ direction (Roig et al. 2012):

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\boldsymbol{u}}=0, & \displaystyle\end{eqnarray}$$
(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\bar{\boldsymbol{u}}}{\unicode[STIX]{x2202}t}+\frac{6}{5}\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}[\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}+(\unicode[STIX]{x1D735}\bar{\boldsymbol{u}})^{\text{T}}]-\frac{12\unicode[STIX]{x1D707}}{h^{2}}\bar{\boldsymbol{u}}+\unicode[STIX]{x1D70C}\boldsymbol{g}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}\boldsymbol{n}, & \displaystyle\end{eqnarray}$$
(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x2202}t}+\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D6FC}=0. & \displaystyle\end{eqnarray}$$

We acknowledge the fact that (3.5)–(3.7) may not represent fully all relevant dynamics, e.g. 3-D effects local to the air–water interface (Oliveira & Meiburg 2011). It is, nonetheless, of interest to determine the extent to which our approximate model can capture as many of the phenomena observed experimentally as possible.

3.3 Curvature and contact angle

The shape of the air–water interface is affected by the wettability of the channel walls. To incorporate this effect, the curvature in (3.6) may be replaced by (Thompson, Juel & Hazel 2014)

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{l}+\unicode[STIX]{x1D705}_{t}\approx \unicode[STIX]{x1D705}_{l}+\frac{\cos \unicode[STIX]{x1D703}_{o}+\cos \unicode[STIX]{x1D703}_{i}}{h}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{l}$ is the lateral curvature (on the $x$ $y$ plane) which is estimated using the diffuse-interface approach based on the volume fraction field (Xie et al. 2016), $\unicode[STIX]{x1D705}_{t}$ is the transverse curvature (across the gap) which is estimated from the meniscus profile (figure 3 a), $\unicode[STIX]{x1D703}_{o}$ is the local contact angle at the outer wall and $\unicode[STIX]{x1D703}_{i}$ is that at the inner wall. The local contact angle along the contact line can vary significantly from the advancing value (at the bubble tail) to the receding value (at the bubble front). The experimental measurements by Antonini et al. (2009) showed that the angle can be interpolated using a third-order polynomial function:

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}=2\frac{\unicode[STIX]{x1D703}_{min}-\unicode[STIX]{x1D703}_{max}}{\unicode[STIX]{x03C0}^{3}}\unicode[STIX]{x1D719}^{3}-3\frac{\unicode[STIX]{x1D703}_{min}-\unicode[STIX]{x1D703}_{max}}{\unicode[STIX]{x03C0}^{2}}\unicode[STIX]{x1D719}^{2}+\unicode[STIX]{x1D703}_{min}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ corresponds to the local contact angle at either the outer or inner wall (i.e. $\unicode[STIX]{x1D703}_{o}$ or $\unicode[STIX]{x1D703}_{i}$ , respectively), $\unicode[STIX]{x1D703}_{max}$ and $\unicode[STIX]{x1D703}_{min}$ denote the advancing and receding contact angles of the corresponding surface, respectively, $\unicode[STIX]{x1D719}$ is the angle between the local velocity vector $\bar{\boldsymbol{u}}$ and the interface normal $\boldsymbol{n}$ (figure 3 b). Figure 4 shows the variation of the local contact angle $\unicode[STIX]{x1D703}$ at the outer and inner walls with respect to $\unicode[STIX]{x1D719}$ . Contact angle hysteresis is promoted by roughness of the channel wall surface (Dussan 1979); we also discuss in § 6 the potential effect of roughness on the presence/absence of thin water films on the channel walls.

Figure 3. Schematic of (a) the cross-gap meniscus profile governed by the local contact angles, $\unicode[STIX]{x1D703}_{o}$ and $\unicode[STIX]{x1D703}_{i}$ , at the outer and inner walls, respectively, and (b) the moving contact line in the $x$ $y$ plane with the local advancing/receding state determined by the intersection angle $\unicode[STIX]{x1D719}$ between the velocity vector $\bar{\boldsymbol{u}}$ and the interface normal $\boldsymbol{n}$ .

Figure 4. Interpolation of contact angles around a contact line at the outer or inner wall of the annulus using a third-order polynomial fitting (Antonini et al. 2009);  $\unicode[STIX]{x1D703}$ denotes the local contact angle and $\unicode[STIX]{x1D719}$ denotes the intersection angle between the local velocity vector and the interface normal as shown in figure 3.

3.4 Model set-up and numerical methods

A numerical experiment was designed as shown in figure 5. The length and width of the 2-D domain were the pipe length $L$ and perimeter $W$ , respectively. This domain was originally filled with water that had an initial velocity $v_{i}$ along the $x$ direction. The inlet condition was defined by the inflow velocity $v_{i}$ uniformly imposed at the left boundary, such that the air flowed into the domain through the middle part ( $W/5$ wide), corresponding to the inflow rate condition of $q_{a}/(q_{a}+q_{w})=1/5$ (table 1), while the water was injected from the remaining portion of the left boundary. This is consistent with the experimental observation that air and water were well segregated in the mixer before entering the annulus. A free-slip condition was used for the two longitudinal sides of the domain, considering that the air phase is away from the pipe bottom and the flow across the pipe bottom is negligible. A hydrostatic pressure condition was imposed along the outlet to account for the gravity-induced pressure gradient across the channel.

Figure 5. The 2-D model set-up for numerical simulation.

The numerical model was solved based on a mixed control-volume finite-element method, which has been validated against a series of benchmark cases of single rising bubbles, coalescence of two bubbles and droplet impacts (Xie et al. 2014, 2016, 2017). The computational domain was discretised into an unstructured grid of triangular $\text{P}_{1}\text{DG}\text{-}\text{P}_{2}$ elements (linear discontinuous velocity between elements and quadratic continuous pressure between elements) (Cotter et al. 2009). A finite volume discretisation of the mass conservation equation and a linear discontinuous Galerkin discretisation of the momentum equation were used with an adaptive implicit/explicit time stepping scheme (Xie et al. 2016). Within each time step, the equations were iterated upon using a pressure projection method until all equations are simultaneously balanced (Pavlidis et al. 2014; Xie et al. 2014). Interface dynamics were captured through a compressive advection-based volume-of-fluid approach, which used a novel and mathematically rigorous nonlinear Petrov–Galerkin method for maintaining sharp interfaces (Pavlidis et al. 2016). Surface tension was modelled using a diffused-interface formulation that can accurately estimate the lateral curvature based on the volume fraction field (Xie et al. 2016). Anisotropic mesh adaptation (Pain et al. 2001) was used to place a finer mesh around the interface during its evolution as required to capture the dynamics (Xie et al. 2014, 2016). The criteria for mesh refinement, based on the volume fraction field, included an interpolation error bound of 0.1, a minimum element size of $W/60$ and maximum element size of $L/6$ and $W/20$ in the $x$ and $y$ directions, respectively. The 2-D simulation outputs were projected onto the original 3-D Cartesian coordinate system to facilitate visual comparison with experimental data.

4 Results

4.1 Experimental observation and numerical simulation

Figures 69 show the formation of air bubbles in the water flowing through the annulus at different inclinations. Despite the limitations of the gap-averaged 2-D numerical model in representing some 3-D effects such as the detailed transverse shape of the air–water interface and the cross-gap velocity/gravity component, very similar flow behaviour in the lateral directions (i.e. $x$ and $y$ directions) was observed in the laboratory experiments and numerical simulations.

In the horizontal case (figure 6), air was conveyed by long bubbles having a tadpole-like shape with a semi-circular cap and a highly stretched tail connected via an apparent neck. These bubbles, separated by liquid slugs, were forced to flow at the top of the annulus due to buoyancy. The tail of a bubble was initially connected to the inlet (figure 6 a) and was stretched as the bubble advanced until it finally disconnected from the inlet (figure 6 b). The bubble then contracted as a result of surface tension (figure 6 c) and continued to move with an essentially steady shape through the co-current flowing liquid towards the outlet (figure 6 d). When the bubble was separated from the inlet, another bubble began to form, resulting in the appearance of a sequence of similarly shaped/sized bubbles moving along the top of the annulus.

Figure 6. The formation and evolution of bubbles in the horizontal annulus (top view): (a) generation of the bubble near the inlet (in annulus subsection 1), (b) disconnection of the bubble (in annulus subsection 1), (c) stabilisation of the bubble after slight contraction (in annulus subsection 2) and (d) translation of the steady bubble to the outlet (in annulus subsection 3). Each figure panel includes the experimental observation, the simulation result and the adaptive unstructured mesh used in the simulation.

In the $1.9^{\circ }$ inclined case (figure 7), bubbles still exhibited a tadpole-like shape but with a much shorter tail, whilst the cap shape remained similar to that in the horizontal case. The bubble detached quickly from the inlet after its formation (figure 7 a). Due to the constant inflow of air, another bubble formed, disconnected and migrated following the leading one. In this way, a stream of distinct bubbles developed as in the horizontal case. However, if the spacing between two bubbles was sufficiently small, the trailing bubble might catch up (figure 7 b) and merge with its upstream neighbour (figure 7 c). The larger bubble formed might further coalesce with another bubble upstream or downstream. Finally, the bubbles that formed relatively stable shapes were transported to the outlet by the flowing water (figure 7 d).

Figure 7. The formation and evolution of bubbles in the $1.9^{\circ }$ inclined annulus (top view): (a) generation of bubbles near the inlet (in annulus subsection 1), (b) interaction of bubbles, with the trailing bubble catching up with the leading bubble (in annulus subsection 2), (c) coalescence of the trailing and leading bubbles (in annulus subsection 2) and (d) translation of steady bubbles to the outlet (in annulus subsection 3). Each figure panel includes the experimental observation, the simulation result and the adaptive unstructured mesh used in the simulation.

There appears to be three stages in the coalescence process. (i) Interacting stage (0.0–0.4 s in figure 8). The trailing bubble with a larger velocity caught up with the leading bubble. As the two bubbles became closer, the trailing bubble experienced a marked change in shape, with its cap becoming more pointed with an increased local curvature at the front. The geometry of the leading bubble only changed slightly with the tail becoming shorter and flatter. The changing shape of the two bubbles was caused by the increased pressure gradient between the two bubble interfaces as they got closer together, during which time the liquid between the two bubbles was gradually pushed out. (ii) Merging stage (0.5–0.6 s in figure 8). The cap of the trailing bubble invaded the rear of the leading bubble. The liquid film between the two bubbles was ruptured and drained as the two bubbles united. The merged bubble exhibited a complex, tortuous shape with the cap and tail connected by a narrow neck and a wider body. (iii) Stabilising stage (0.7–1.0 s in figure 8). The merged larger bubble underwent deformation and adjusted itself from a complex, unsteady shape to the stable tadpole-like form.

Figure 8. Simulation results showing the coalescence process for two bubbles in the $1.9^{\circ }$ inclined annulus.

Figure 9. The formation and evolution of bubbles in the $4.6^{\circ }$ inclined annulus (top view): (a) generation and interaction of bubbles close to the inlet (in annulus subsection 1), (b) coalescence of the trailing, middle and leading bubbles (in annulus subsection 1), (c) stabilisation of merged larger bubbles (in annulus subsection 2) and (d) translation of steady bubbles to the outlet (in annulus subsection 3). Each figure panel includes the experimental observation, the simulation result and the adaptive unstructured mesh used in the simulation.

In the $4.6^{\circ }$ inclined case (figure 9), bubbles became more compact and almost lost the tail structure although the cap shape remained similar to that in the horizontal case. Bubbles separated from the inlet quickly and again formed a train (figure 9 a). The bubbles could merge at an earlier stage by virtue of their smaller spacing and simultaneous coalescence of three bubbles could occur (figure 9 a,b). During the coalescence of three bubbles, the leading bubble remained a semi-circular cap, whereas its rear was flattened; the middle bubble was between the leading and trailing bubbles and showed a sharpened front and a flattened rear; the trailing bubble displayed a narrowed front invading into the back of middle bubble. As merging occurred, the liquid films between bubbles broke and dissipated, such that a united larger bubble was formed with a highly tortuous, unsteady shape. The coalesced larger bubble gradually stabilised itself (figure 9 c) and then migrated to the exit (figure 9 d).

The statistics of the normalised length, cap diameter and terminal velocity of stabilised bubbles in the cases of different inclinations are compared in figure 10. Data from ten fully developed bubbles were collected in the region 2–5 m downstream of the inlet, where the bubbles were not influenced by either inlet or outlet effects. As expected from the discussion above, the bubble length reduces significantly as the inclination of the annulus increases (figure 10 a), due to the deterioration of tail structure. However, the cap diameter, measured based on the curvature of the cap front, only shows a slight increase with the inclination of the annulus (figure 10 b). The terminal velocity of the bubbles also increases slightly with inclination (figure 10 c). An excellent quantitative match with the experimental measurements was achieved by the numerical model. The pseudo-periodic nature of the slug flow and variation of bubble properties as observed in the experiment was also well captured in the simulation (figure 10). Although the inlet condition was stationary, the slug flow pattern varied constantly.

Figure 10. Variation of bubble properties with the inclination angle $\unicode[STIX]{x1D6FD}$ : (a) normalised bubble length $l/D$ , (b) normalised bubble cap diameter $d/D$ and (c) normalised bubble terminal velocity $u/v_{i}$ .

4.2 Effects of gap thickness and pipe diameter

We used the numerical model to further elucidate the effects of gap thickness and pipe diameter on bubbles in thin annuli. First, we set the pipe diameter to be $D=0.137$  m, the air-to-water inflow rate ratio to be $q_{a}/q_{w}=1/4$ and the superficial inflow velocity to be $v_{i}=0.162~\text{m}~\text{s}^{-1}$ . We varied the gap thickness $h$ from 0.001 to 0.005 m, for which the flow remained to be laminar (Beavers et al. 1970) and cross-gap gravity effect was not significant. As shown in figure 11(a), when $h=0.001$  m (i.e. $h/D=0.0073$ ), the bubbles in the horizontal and inclined annuli are very similar, having a slender oval shape with a semi-circular front and a gradually narrowed rear (no discernible neck existed). In the horizontal case, as $h$ becomes larger, the tadpole-like shape becomes more pronounced with a much longer tail and a more distinct neck where the tail joins the cap. If the annulus is inclined, as $h$ increases, the tail tends to vanish and the bubble length decreases significantly. In contrast to the marked variation of bubble length with $h$ , the cap diameter and terminal velocity are much less sensitive to $h$ (figure 11 b,c).

Figure 11. Variation of bubble properties with the gap thickness $h$ while the pipe diameter $D$ is fixed to be $0.137$ m: (a) normalised bubble length $l/D$ , (b) normalised bubble cap diameter $d/D$ and (c) normalised bubble terminal velocity $u/v_{i}$ .

Next, we explored the effect of changing pipe diameter $D=0.05$ –0.4 m (the pipe length $L$ also varied with $D$ according to a fixed ratio $L/D=43.8$ ), keeping $h=0.0035$  m and the same inflow conditions as before. Note that, in figure 12, the bubble length and cap diameter in each case are normalised by the corresponding pipe diameter. As shown in figure 12(a), the ratio of bubble length to the pipe diameter generally does not vary significantly with the change of $D$ . However, in the horizontal annulus with $D=0.4$  m (i.e. $h/D=0.0087$ ), the bubble has a very long, extremely stretched, wavy tail, which does not disconnect from the inlet; the bubble length, therefore, corresponds to the whole pipe length. With the decrease of $D$ (i.e. increase of $h/D$ ), the cap diameter relative to the pipe diameter increases while the bubble velocity decreases (figure 12 b,c).

Figure 12. Variation of bubble properties with the pipe diameter $D$ while the gap thickness $h$ is fixed to be 0.0035 m: (a) normalised bubble length $l/D$ , (b) normalised bubble cap diameter $d/D$ and (c) normalised bubble terminal velocity $u/v_{i}$ .

5 Flow mechanisms

We use the numerical model to investigate the pressure and velocity fields around the bubbles, which determine the bubble motion, shape and dynamics, and can thus provide clues for interpreting the mechanisms that control these behaviours. Figure 13(a) shows the typical steady bubbles for three different inclinations, and figures 13(b) and 13(c) give their pressure and velocity fields, respectively. To facilitate comparison, the absolute pressure at the foremost point of each bubble is chosen as the reference pressure in each case, so that the pressure contour gives the differential pressure between the absolute pressure and the reference pressure. It can be seen that inside the bubble, the pressure is generally constant due to the very low air density. The pressure in the liquid behind the bubble is higher than that in front of the bubble, driving the bubble to flow towards the outlet. The velocity of air inside the bubble is significantly higher than that of the surrounding liquid.

Figure 13. Simulation results showing (a) typical steady bubbles in annuli with different inclinations, and their surrounding (b) pressure and (c) velocity fields. To facilitate comparison, the absolute pressure at the frontal point of each bubble is chosen as the reference pressure for that case, so that the pressure contour gives the differential pressure between the absolute pressure and the reference pressure.

Figure 14. (a) The unrolled geometry of a bubble in the $x$ $y$ plane, and the views of (b) a transversal cross-section in the $y$ $z$ plane and (c) a longitudinal cross-section in the $x$ $z$ plane.

We explore the mechanisms controlling the bubble shape and motion through a simple analysis. We define two polar coordinate systems: one based on the unrolled 2-D bubble cap geometry in the $x$ $y$ plane (figure 14 a), and the other based on the transverse cross-section of the pipe in the $y$ $z$ plane (figure 14 b). By applying Bernoulli’s equation to the steady flow around the bubble cap, assumed to be irrotational with the surface tension being neglected, the following relationship is obtained:

(5.1) $$\begin{eqnarray}\displaystyle \frac{u_{\unicode[STIX]{x1D713}}^{2}}{2}=g\left[\frac{D}{2}(1-\cos \unicode[STIX]{x1D6FE})\cos \unicode[STIX]{x1D6FD}+\frac{d}{2}(1-\cos \unicode[STIX]{x1D713})\sin \unicode[STIX]{x1D6FD}\right], & & \displaystyle\end{eqnarray}$$

where $u_{\unicode[STIX]{x1D713}}$ is the relative velocity of the surrounding liquid to the bubble at a point with a polar angle of $\unicode[STIX]{x1D713}$ , and $\unicode[STIX]{x1D6FE}$ is the polar angle in the coordinate system of the transverse cross-section. We assume that flow around the bubble cap is similar to that seen around a circular cylinder to a first approximation (Collins 1967a ). This flow field can be solved using the stream function and velocity potential, with the following relationship obtained (Wilkes 2006):

(5.2) $$\begin{eqnarray}\displaystyle \frac{u_{\unicode[STIX]{x1D713}}^{2}}{u_{\infty }^{2}\sin ^{2}\unicode[STIX]{x1D713}}=4, & & \displaystyle\end{eqnarray}$$

where $u_{\infty }$ is the far-field relative velocity of the flow to the bubble. By combining (5.1), (5.2) and the relation $(d/2)\sin \unicode[STIX]{x1D713}=(D/2)\unicode[STIX]{x1D6FE}$ (see figure 14 a,b), we obtain

(5.3) $$\begin{eqnarray}\displaystyle u_{\infty }^{2}=\left(\frac{1-\cos \unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6FE}^{2}}\frac{d}{D}\cos \unicode[STIX]{x1D6FD}+\frac{1-\cos \unicode[STIX]{x1D713}}{\sin ^{2}\unicode[STIX]{x1D713}}\sin \unicode[STIX]{x1D6FD}\right)\frac{gd}{4}. & & \displaystyle\end{eqnarray}$$

In the limit of small $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D713}$ , equation (5.3) becomes

(5.4) $$\begin{eqnarray}\displaystyle u_{\infty }\approx \frac{1}{2}\sqrt{\left(\frac{d}{D}\cos \unicode[STIX]{x1D6FD}+\sin \unicode[STIX]{x1D6FD}\right)\frac{gd}{2}}. & & \displaystyle\end{eqnarray}$$

If $r\rightarrow \infty$ , $\unicode[STIX]{x1D6FD}=90^{\circ }$ , then solution (5.4) converges to that for plane bubbles rising through a quiescent liquid between vertical, parallel flat plates, i.e. $0.5(gd/2)^{1/2}$ (Davies & Taylor 1950; Collins 1965a ).

For the simultaneous gas–liquid flow here, the bubble terminal velocity $u$ is calculated as (Griffith & Wallis 1961)

(5.5) $$\begin{eqnarray}\displaystyle u=u_{\infty }+v_{i}\approx \frac{1}{2}\sqrt{\left(\frac{d}{D}\cos \unicode[STIX]{x1D6FD}+\sin \unicode[STIX]{x1D6FD}\right)\frac{gd}{2}}+v_{i}, & & \displaystyle\end{eqnarray}$$

where $v_{i}$ is the superficial inflow velocity (table 1). We substitute the measured diameters of bubble cap from the experimental/numerical results into (5.5) to derive the predicted bubble terminal velocity $u_{pred}$ . The predicted values are compared to the actual velocity $u_{meas}$ measured in the experiment and simulation. It can be seen that our simplified analysis gives a good estimation of the bubble velocity in thin annuli for different inclination angles, pipe diameters and gap thicknesses (figure 15).

Figure 15. Comparison of the predicted terminal velocity $u_{pred}$ from (5.5) and the measured terminal velocity $u_{meas}$ from experiments and simulations for different configurations of pipe diameter $D$ , gap thickness $h$ and inclination angle $\unicode[STIX]{x1D6FD}$ . The velocity values are normalised by the superficial inflow velocity $v_{i}$ .

We then identify the following dimensionless groups of the annulus system:

(5.6a-c ) $$\begin{eqnarray}\displaystyle \hat{Fr}_{\Vert }^{-1}=\frac{\sqrt{gD\sin \unicode[STIX]{x1D6FD}}}{v_{i}},\quad \hat{Fr}_{\bot }^{-1}=\frac{\sqrt{gD\cos \unicode[STIX]{x1D6FD}}}{v_{i}},\quad \hat{\unicode[STIX]{x1D716}}=\displaystyle \frac{h}{D}, & & \displaystyle\end{eqnarray}$$

where the two inverse Froude numbers, $\hat{Fr}_{\Vert }^{-1}$ and $\hat{Fr}_{\bot }^{-1}$ , are related to the effects of gravity parallel and perpendicular to the flow direction, respectively, and $\hat{\unicode[STIX]{x1D716}}$ characterises the effects of confinement associated with the annulus. We also characterise the bubble shape and motion based on the bubble aspect ratio $l/d$ , the bubble Reynolds number $Re=\unicode[STIX]{x1D70C}_{w}u_{\infty }d/\unicode[STIX]{x1D707}_{w}$ and the gap Reynolds number $Re(h/d)^{2}$ . Here, $u_{\infty }$ is used because the bubble behaviour is governed by the relative motion between the two fluids.

As shown in figure 16(a), as $\hat{Fr}_{\Vert }^{-1}$ increases, the bubble shape becomes less elongated, i.e. $l/d$ approaches unity, since the gravitational component parallel to the flow direction causes an increased pressure at the bubble rear, which tends to separate the bubble and reduce the bubble length. As $\hat{Fr}_{\bot }^{-1}$ increases, $Re$ increases (figure 16 b), indicating inertia becomes more dominant. This is because, for the annuli with small inclinations to the horizontal, the gravitational component perpendicular to the flow direction mainly controls the flow field around the bubble cap, as can be expected from (5.1). As the annulus confinement decreases (i.e. $\hat{\unicode[STIX]{x1D716}}$ increases), the gap Reynolds number $Re(h/d)^{2}$ increases (figure 16 c), due to the reduced viscous effect across the gap, reaching a plateau with increasing $\hat{\unicode[STIX]{x1D716}}$ .

Figure 16. (a) Variation of the bubble aspect ratio $l/d$ with the inverse Froude number $\hat{Fr}_{\Vert }^{-1}$ . (b) Variation of the bubble Reynolds number $Re$ with the inverse Froude number $\hat{Fr}_{\bot }^{-1}$ . (c) Variation of the gap Reynolds number $Re(h/d)^{2}$ with the confinement ratio $\hat{\unicode[STIX]{x1D716}}$ . Refer to figure 15 for the legend.

6 Discussion and conclusions

In this paper, we presented an experimental and numerical investigation into the shape and motion of gas bubbles in a liquid flowing through a horizontal or slightly inclined thin annulus. We have focused on the regime where inertia is important (the gap Reynolds number is close to or larger than unity) and the bubble behaviour is governed by the complex interplay of inertia, gravity, viscosity and surface tension. Bubbles in the horizontal annulus developed a unique ‘tadpole-like’ shape featuring a semi-circular cap and a highly stretched tail. As the annulus was inclined with respect to the horizontal, the length of the bubble decreased. We developed a gap-averaged, 2-D numerical model to represent the 3-D flow dynamics, which achieved a close match to the experimental data for different small inclinations. The numerical model was used to further elucidate the effects of gap thickness and pipe diameter on the bubble evolution in thin annuli. We found that the bubble velocity is strongly correlated to the cap structure, but is independent of the bubble length, as has also been reported for bubbles in tubes (Davies & Taylor 1950; Zukoski 1966; Fabre & Liné 1992; Viana et al. 2003) and between parallel flat plates (Grace & Harrison 1967; Maneri & Zuber 1974; Hills 1975; Couet & Strumolo 1987). We reported that the elongated bubble shape in horizontal annuli is due to the buoyancy which causes the bubbles to spread along the top of the annulus. The gravitational component along the flow direction, which increases as the annulus is inclined, impinges the liquid slug and causes a reduction of the bubble tail. The gravitational component perpendicular to the flow direction controls the bubble motion and the cap structure. These mechanisms produce the unique tadpole-like shape with a sharp tail tip because of the cross-sectional curvature of the annulus channel, in contrast to the bubble shape with a tongue-like rear seen in flows between parallel flat plates (Maneri & Zuber 1974; Couet & Strumolo 1987).

It is remarkable that the 2-D numerical model well captured the bubble evolution and interaction behaviour as observed in the experimental annulus. It is still worth mentioning that some complex 3-D effects were not fully represented in the gap-averaged formulation, such as the detailed transverse shape of the bubble cap and tail within the gap as well as the influence of the cross-gap velocity/gravity component (Oliveira & Meiburg 2011), which require direct three-dimensional but computationally very expensive numerical simulations. Such limitations of the 2-D model may explain the discrepancy of the bubble geometry, e.g. in the region near the bubble neck (figure 6), between our experimental and numerical results. The gap-averaged 2-D model assumed an absence of liquid films between the bubble and channel walls. We examined different mathematical formulations that account for the presence of liquid films (Pitts 1980; McLean & Saffman 1981; Park & Homsy 1984; Reinelt 1987a ; Kopf-Sill & Homsy 1988), which however resulted in a poor prediction of the experimentally observed bubble patterns (the simulated bubbles were much smaller, probably due to an overestimation of the curvature of the air–water interface across the gap). Saffman & Tanveer (1989) also reported that the thin film hypothesis, compared to the contact angle assumption, gave a less consistent prediction of bubble shapes in Hele-Shaw cells. Furthermore, we did not see any evidence of such films in the experiments. The thickness of thin films (if present under no-gravity conditions) is estimated to be $\unicode[STIX]{x1D6FF}/h\approx 0.67(\unicode[STIX]{x1D707}_{w}u_{\infty }/\unicode[STIX]{x1D70E})^{2/3}$ (Park & Homsy 1984; Reinelt 1987a ; Klaseboer, Gupta & Manica 2014), which may however be reduced further in the top region of the annulus under gravity-induced drainage downwards along the peripheral direction (Tso & Sugawara 1990; Paras & Karabelas 1991). Thus, we expect that $\unicode[STIX]{x1D6FF}$ may be close to the scale of the surface roughness of the Plexiglas and stainless steel walls, i.e. around 0.01 mm (Li et al. 2018). We therefore do not expect the films to persist in a stable way. The observed contact angle hysteresis in our study may also be attributed to the presence of pronounced surface roughness effects (Dussan 1979). The phenomenon that air wets the channel wall has also been found between the air bubble and the upper wall of a horizontal/slightly inclined tube (Fabre & Liné 1992) or Hele-Shaw cell (Maneri & Zuber 1974). We suggest that a good approximation of the gap-averaged 2-D model to 3-D flow holds if the cross-gap effects are minor.

Acknowledgements

We would like to acknowledge Statoil ASA for funding Q.L., M.D.J., C.C.P., O.K.M. and A.H.M., and for granting permission to publish this work. We thank K. Årland from Statoil ASA for project management. O.K.M. and C.C.P. acknowledge the funding provided by the Engineering and Physical Sciences Research Council (EPSRC) through the Programme Grant MEMPHIS (grant EP/K003976/1). D.P. is grateful to the EPSRC (grant EP/M012794/1). P.S. thanks ExxonMobil and the EPSRC (grant EP/R005761/1) for funding. No data were generated in the course of this work. For further information, please contact the corresponding author Q.L. (email: q.lei12@imperial.ac.uk or qinghua.lei@erdw.ethz.ch), the AMCG Group (http://www.imperial.ac.uk/earth-science/research/research-groups/amcg/), the NORMS group (http://www.imperial.ac.uk/earth-science/research/research-groups/norms) or the Matar Fluids Group (http://www.imperial.ac.uk/matar-fluids-group).

References

Antonini, C., Carmona, F. J., Pierce, E., Marengo, M. & Amirfazli, A. 2009 General methodology for evaluating the adhesion force of drops and bubbles on solid surfaces. Langmuir 25 (11), 61436154.
Beavers, G. S., Sparrow, E. M. & Magnuson, R. A. 1970 Experiments on the breakdown of laminar flow in a parallel-plate channel. Intl J. Heat Mass Transfer 13, 809815.
Brackbill, J. U., Kothe, D. B. & Zemach, C. 1992 A continuum method for modeling surface tension. J. Comput. Phys. 100 (2), 335354.
Bush, J. W. M. 1997 The anomalous wake accompanying bubbles rising in a thin gap: a mechanically forced Marangoni flow. J. Fluid Mech. 352, 283303.
Bush, J. W. M. & Eames, I. 1998 Fluid displacement by high Reynolds number bubble motion in a thin gap. Intl J. Multiphase Flow 24 (3), 411430.
Collins, R. 1965a A simple model of the plane gas bubble in a finite liquid. J. Fluid Mech. 22, 763771.
Collins, R. 1965b Structure and behaviour of wakes behind two-dimensional air bubbles in water. Chem. Engng Sci. 20, 851853.
Collins, R. 1967a The cycloidal-cap bubble: a neglected solution in the theory of large plane gas bubbles in liquids. Chem. Engng Sci. 22, 8997.
Collins, R. 1967b The effect of a containing cylindrical boundary on the velocity of a large gas bubble in a liquid. J. Fluid Mech. 28 (1), 97112.
Collins, R., de Moraes, F. F., Davidson, J. F. & Harrison, D. 1978 The motion of a large gas bubble rising through liquid flowing in a tube. J. Fluid Mech. 89, 497514.
Cotter, C. J., Ham, D. A., Pain, C. C. & Reich, S. 2009 LBB stability of a mixed Galerkin finite element pair for fluid flow simulations. J. Comput. Phys. 228 (2), 336348.
Couet, B. & Strumolo, G. S. 1987 The effects of surface tension and tube inclination on a two-dimensional rising bubble. J. Fluid Mech. 184, 114.
Cueto-Felgueroso, L. & Juanes, R. 2014 A phase-field model of two-phase Hele-Shaw flow. J. Fluid Mech. 758, 522552.
Das, G., Das, P. K., Purohit, N. K. & Mitra, A. K. 1998 Rise velocity of a Taylor bubble through concentric annulus. Chem. Engng Sci. 53 (5), 977993.
Davies, R. M. & Taylor, G. I. 1950 The mechanics of large bubbles rising through extended liquids and through liquids in tubes. Proc. R. Soc. Lond. A 200 (1062), 375390.
Dussan, E. B. 1979 On the spreading of liquids on solid surfaces: static and dynamic contact lines. Ann. Rev. Fluid Mech. 11, 371400.
Eck, W. & Siekmann, J. 1978 On bubble motion in a Hele-Shaw cell, a possibility to study two-phase flows under reduced gravity. Ing.-Arch. 47, 153168.
Ekberg, N. P., Ghiaasiaan, S. M., Abdel-Khalik, S. I., Yoda, M. & Jeter, S. M. 1999 Gas–liquid two-phase flow in narrow horizontal annuli. Nucl. Engng Des. 192, 5980.
Fabre, J. 2016 A long bubble rising in still liquid in a vertical channel: a plane inviscid solution. J. Fluid Mech. 797, R4.
Fabre, J. & Liné, A. 1992 Modeling of two-phase slug flow. Annu. Rev. Fluid Mech. 24, 2146.
Figueroa-Espinoza, B. & Fabre, J. 2011 Taylor bubble moving in a flowing liquid in vertical channel: transition from symmetric to asymmetric shape. J. Fluid Mech. 679, 432454.
Filella, A., Ern, P. & Roig, V. 2015 Oscillatory motion and wake of a bubble rising in a thin-gap cell. J. Fluid Mech. 778, 6088.
Gondret, P. & Rabaud, M. 1997 Shear instability of two-fluid parallel flow in a Hele-Shaw cell. Phys. Fluids 9 (11), 32673274.
Grace, J. R. & Harrison, D. 1967 The influence of bubble shape on the rising velocities of large bubbles. Chem. Engng Sci. 22, 13371347.
Griffith, P. & Wallis, G. B. 1961 Two-phase slug flow. Trans. ASME J. Heat Transfer 83, 307318.
Hills, J. H. 1975 The two-dimensional elliptical cap bubble. J. Fluid Mech. 68, 503512.
Huisman, S. G., Ern, P. & Roig, V. 2012 Interactions and coalescence of large bubbles rising in a thin gap. Phys. Rev. E 85, 027302.
Kelessidis, V. C. & Dukler, A. E. 1990 Motion of large gas bubbles through liquids in vertical concentric and eccentric annuli. Intl J. Multiphase Flow 16 (3), 375390.
Kelley, E. & Wu, M. 1997 Path instability of rising air bubbles in a Hele-Shaw cell. Phys. Rev. Lett. 79 (7), 12651268.
Klaseboer, E., Gupta, R. & Manica, R. 2014 An extended Bretherton model for long Taylor bubbles at moderate capillary numbers. Phys. Fluids 26, 032107.
Kopf-Sill, A. R. & Homsy, G. M. 1988 Bubble motion in a Hele-Shaw cell. Phys. Fluids 31, 1826.
Lazarek, G. M. & Littman, H. 1974 The pressure field due to a large circular capped air bubble rising in water. J. Fluid Mech. 66, 673687.
Li, G., Su, H., Kuhn, U., Meusel, H., Ammann, M., Shao, M., Pöschl, U. & Cheng, Y. 2018 Technical note: influence of surface roughness and local turbulence on coated-wall flow tube experiments for gas uptake and kinetic studies. Atmos. Chem. Phys. 18 (4), 26692686.
Maneri, C. C. & Zuber, N. 1974 An experimental study of plane bubbles rising at inclination. Intl J. Multiphase Flow 1, 623645.
Maruvada, S. R. K. & Park, C. W. 1996 Retarded motion of bubbles in Hele-Shaw cells. Phys. Fluids 8 (12), 32293233.
Maxworthy, T. 1986 Bubble formation, motion and interaction in a Hele-Shaw cell. J. Fluid Mech. 173, 95114.
McLean, J. W. & Saffman, P. G. 1981 The effect of surface tension on the shape of fingers in a Hele Shaw cell. J. Fluid Mech. 102, 455469.
Meiburg, E. 1989 Bubbles in a Hele-Shaw cell: numerical simulation of three-dimensional effects. Phys. Fluids A 1 (6), 938946.
Oliveira, R. M. & Meiburg, E. 2011 Miscible displacements in Hele-Shaw cells: three-dimensional Navier–Stokes simulations. J. Fluid Mech. 687, 431460.
Pain, C. C., Umpleby, A. P., de Oliveira, C. R. E. & Goddard, A. J. H. 2001 Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations. Comput. Meth. Appl. Mech. Engng 190 (29–30), 37713796.
Paras, S. V. & Karabelas, A. J. 1991 Properties of the liquid layer in horizontal annular flow. Intl J. Multiphase Flow 17 (4), 439454.
Park, C.-W. & Homsy, G. M. 1984 Two-phase displacement in Hele Shaw cells: theory. J. Fluid Mech. 139, 291308.
Pavlidis, D., Gomes, J. L. M. A., Xie, Z., Percival, J. R., Pain, C. C. & Matar, O. K. 2016 Compressive advection and multi-component methods for interface-capturing. Intl J. Numer. Meth. Fluids 80 (4), 256282.
Pavlidis, D., Xie, Z., Percival, J. R., Gomes, J. L. M. A., Pain, C. C. & Matar, O. K. 2014 Two- and three-phase horizontal slug flow simulations using an interface-capturing compositional approach. Intl J. Multiphase Flow 67, 8591.
Piedra, S., Ramos, E. & Herrera, J. R. 2015 Dynamics of two-dimensional bubbles. Phys. Rev. E 91 (6), 063013.
Pitts, E. 1980 Penetration of fluid into a Hele-Shaw cell: the Saffman–Taylor experiment. J. Fluid Mech. 97, 5364.
Reinelt, D. A. 1987a Interface conditions for two-phase displacement in Hele-Shaw cells. J. Fluid Mech. 183, 219234.
Reinelt, D. A. 1987b The rate at which a long bubble rises in a vertical tube. J. Fluid Mech. 175, 557565.
Roig, V., Roudet, M., Risso, F. & Billet, A.-M. 2012 Dynamics of a high-Reynolds-number bubble rising within a thin gap. J. Fluid Mech. 707, 444466.
Saffman, P. G. & Tanveer, S. 1989 Prediction of bubble velocity in a Hele-Shaw cell: thin film and contact angle effects. Phys. Fluids A 1 (2), 219223.
Saffman, P. G. & Taylor, G. I. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. Lond. A 245 (1242), 312329.
Tanveer, S. 1986 The effect of surface tension on the shape of a Hele-Shaw cell bubble. Phys. Fluids 29 (11), 35373548.
Tanveer, S. 1987 New solutions for steady bubbles in a Hele-Shaw cell. Phys. Fluids 30 (3), 651658.
Tanveer, S. & Saffman, P. G. 1987 Stability of bubbles in a Hele-Shaw cell. Phys. Fluids 30 (9), 26242635.
Taylor, G. I. & Saffman, P. G. 1959 A note on the motion of bubbles in a Hele-Shaw cell and porous medium. Q. J. Mech. Appl. Math. 12, 265279.
Thompson, A. B., Juel, A. & Hazel, A. L. 2014 Multiple finger propagation modes in Hele-Shaw channels of variable depth. J. Fluid Mech. 746, 123164.
Thompson, B. W. 1968 Secondary flow in a Hele-Shaw cell. J. Fluid Mech. 372, 2544.
Tso, C. P. & Sugawara, S. 1990 Film thickness prediction in a horizontal annular two-phase flow. Intl J. Multiphase Flow 16 (5), 867884.
Vanden-Broeck, J.-M. 1984 Rising bubbles in a two-dimensional tube with surface tension. Phys. Fluids 27 (11), 26042607.
Vanden-Broeck, J.-M. 1992 Rising bubbles in a two-dimensional tube: asymptotic behavior for small values of the surface tension. Phys. Fluids A 4 (11), 23322334.
Viana, F., Pardo, R., Yánez, R., Trallero, J. L. & Joseph, D. D. 2003 Universal correlation for the rise velocity of long gas bubbles in round pipes. J. Fluid Mech. 494, 379398.
Wang, X., Klaasen, B., Degrève, J., Blanpain, B. & Verhaeghe, F. 2014 Experimental and numerical study of buoyancy-driven single bubble dynamics in a vertical Hele-Shaw cell. Phys. Fluids 26 (12), 123303.
Wang, X., Klaasen, B., Degrève, J., Mahulkar, A., Heynderickx, G., Reyniers, M.-F., Blanpain, B. & Verhaeghe, F. 2016 Volume-of-fluid simulations of bubble dynamics in a vertical Hele-Shaw cell. Phys. Fluids 28, 053304.
White, E. T. & Beardmore, R. H. 1962 The velocity of rise of single cylindrical air bubbles through liquids contained in vertical tubes. Chem. Engng Sci. 17 (5), 351361.
Wilkes, J. O. 2006 Fluid Mechanics for Chemical Engineers with Microfluidics and CFD, 2nd edn. Prentice Hall.
Xie, Z., Hewitt, G. F., Pavlidis, D., Salinas, P., Pain, C. C. & Matar, O. K. 2017 Numerical study of three-dimensional droplet impact on a flowing liquid film in annular two-phase flow. Chem. Engng Sci. 166, 303312.
Xie, Z., Pavlidis, D., Percival, J. R., Gomes, J. L. M. A., Pain, C. C. & Matar, O. K. 2014 Adaptive unstructured mesh modelling of multiphase flows. Intl J. Multiphase Flow 67, 104110.
Xie, Z., Pavlidis, D., Salinas, P., Percival, J. R., Pain, C. C. & Matar, O. K. 2016 A balanced-force control volume finite element method for interfacial flows with surface tension using adaptive anisotropic unstructured meshes. Comput. Fluids 138, 3855.
Zukoski, E. E. 1966 Influence of viscosity, surface tension, and inclination angle on motion of long bubbles in closed tubes. J. Fluid Mech. 25, 821837.