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 Reference Brackbill, Kothe and Zemach1992). Note (3.1)–(3.3) are written in a form following Bush (Reference Bush1997) 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 Reference Gondret and Rabaud1997):

(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 Reference Beavers, Sparrow and Magnuson1970).

We derive a coupled set of gap-averaged equations by integrating equations (3.1)–(3.3) along the $z$ direction (Roig et al. Reference Roig, Roudet, Risso and Billet2012):

(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 Reference Oliveira and Meiburg2011). 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 Reference Thompson, Juel and Hazel2014)

(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. Reference Xie, Pavlidis, Salinas, Percival, Pain and Matar2016), $\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. (Reference Antonini, Carmona, Pierce, Marengo and Amirfazli2009) 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 Reference Dussan1979); 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. Reference Antonini, Carmona, Pierce, Marengo and Amirfazli2009);  $\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. Reference Xie, Pavlidis, Percival, Gomes, Pain and Matar2014, Reference Xie, Pavlidis, Salinas, Percival, Pain and Matar2016, Reference Xie, Hewitt, Pavlidis, Salinas, Pain and Matar2017). 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. Reference Cotter, Ham, Pain and Reich2009). 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. Reference Xie, Pavlidis, Salinas, Percival, Pain and Matar2016). Within each time step, the equations were iterated upon using a pressure projection method until all equations are simultaneously balanced (Pavlidis et al. Reference Pavlidis, Xie, Percival, Gomes, Pain and Matar2014; Xie et al. Reference Xie, Pavlidis, Percival, Gomes, Pain and Matar2014). 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. Reference Pavlidis, Gomes, Xie, Percival, Pain and Matar2016). 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. Reference Xie, Pavlidis, Salinas, Percival, Pain and Matar2016). Anisotropic mesh adaptation (Pain et al. Reference Pain, Umpleby, de Oliveira and Goddard2001) was used to place a finer mesh around the interface during its evolution as required to capture the dynamics (Xie et al. Reference Xie, Pavlidis, Percival, Gomes, Pain and Matar2014, Reference Xie, Pavlidis, Salinas, Percival, Pain and Matar2016). 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.

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 Reference Collins1967a ). This flow field can be solved using the stream function and velocity potential, with the following relationship obtained (Wilkes Reference Wilkes2006):

(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 Reference Davies and Taylor1950; Collins Reference Collins1965a ).

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

(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 Reference Davies and Taylor1950; Zukoski Reference Zukoski1966; Fabre & Liné Reference Fabre and Liné1992; Viana et al. Reference Viana, Pardo, Yánez, Trallero and Joseph2003) and between parallel flat plates (Grace & Harrison Reference Grace and Harrison1967; Maneri & Zuber Reference Maneri and Zuber1974; Hills Reference Hills1975; Couet & Strumolo Reference Couet and Strumolo1987). 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 Reference Maneri and Zuber1974; Couet & Strumolo Reference Couet and Strumolo1987).

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 Reference Oliveira and Meiburg2011), 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 Reference Pitts1980; McLean & Saffman Reference McLean and Saffman1981; Park & Homsy Reference Park and Homsy1984; Reinelt Reference Reinelt1987a ; Kopf-Sill & Homsy Reference Kopf-Sill and Homsy1988), 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 (Reference Saffman and Tanveer1989) 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 Reference Park and Homsy1984; Reinelt Reference Reinelt1987a ; Klaseboer, Gupta & Manica Reference Klaseboer, Gupta and Manica2014), which may however be reduced further in the top region of the annulus under gravity-induced drainage downwards along the peripheral direction (Tso & Sugawara Reference Tso and Sugawara1990; Paras & Karabelas Reference Paras and Karabelas1991). 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. Reference Li, Su, Kuhn, Meusel, Ammann, Shao, Pöschl and Cheng2018). 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 Reference Dussan1979). 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é Reference Fabre and Liné1992) or Hele-Shaw cell (Maneri & Zuber Reference Maneri and Zuber1974). 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.