## Abstract

The steady axisymmetric behaviour of a relatively small bubble moving with a flowing liquid in a straight round tube is studied by computationally solving the nonlinear Navier–Stokes equations, using a Galerkin finite-element method with boundary-fitted mesh, for wide ranges of capillary number *C**a* and Reynolds number *R**e*. Here a bubble is considered relatively small when its volume-equivalent radius is less than that of the tube. At small values of *R**e*, the velocity of a bubble increases with bubble size for large values of *C**a* but decreases with bubble size for small values of *C**a*. At large values of *R**e*, however, a bubble of large size appears to move at a slower velocity for any given value of *C**a*. When *R**e* is large (e.g. *R**e* = 100) and *C**a* > 0.1, a bubble of radius greater than half of the tube radius moves at a velocity that seems to be independent of bubble size. The strong inertial effect at large *R**e* makes a small bubble of radius greater than a quarter of the tube radius to deform into a noticeable oblate shape as *C**a* increases from very small value, and then to be elongated into a bullet shape with further increasing *C**a* after *C**a* reaches an intermediate value. Even very small bubbles (e.g. of radius equal to one-tenth of the tube radius) can still be significantly deformed provided that the value of *C**a* is adequately large. Despite significant shape deformations that may still occur, bubbles of radius less than a quarter of that of the tube almost always move at the same velocity as that of the local liquid flow at the tube centreline (i.e. twice that of the average liquid velocity), regardless the values of *R**e* and *C**a*. This fact suggests that very small bubbles are basically carried by the local liquid flow.

## 1. Introduction

Amidst a flowing liquid in a tube, a gas bubble is subjected to hydrodynamic stresses that tend to deform it from the spherical shape otherwise maintained by surface tension and to make it move relative to the tube wall faster than the average liquid flow speed. Besides the interest in advancing our fundamental understanding of two-phase flow phenomena, the subject of bubble behaviour in a tube with flowing liquid can be brought to bear on important applications in the micro-gravity environment. When studying a gas bubble in a tube moving with flowing liquid (with negligible buoyancy effect), however, most authors have restricted their attentions to long bubbles whose sizes are much larger than the tube size (e.g. Fairbrother & Stubbs 1935; Bretherton 1961; Taylor 1961; Cox 1962; Park & Homsy 1984; Reinelt & Saffman 1985; Schwartz *et al.* 1986; Martinez & Udell 1989; Giavedoni & Saita 1997; Feng 2009).

The purpose of the present work is to explore the behaviour of relatively small bubbles moving in a tube with flowing liquid (which seems to be lacking in the published literature1), as a natural extension to the recent computational analysis of Feng (2009) for the situation of long bubbles. Unlike long bubbles whose basic fluid mechanic behaviour is independent of the bubble size, small bubbles in relatively large tubes are expected to exhibit significant behavioural change as their (relative) size varies.

## 2. Mathematical problem description

Considered here is an isolated gas bubble moving at a speed *U* (relative to the tube wall) with a liquid of constant density *ρ*, viscosity *μ* and surface tension *γ* in a straight round tube of radius *R*. As in Feng (2009), a reference frame moving with the bubble is used with the coordinate origin fixed at the bubble’s centroid. With the hydrodynamic stresses due to the flow of gas inside the bubble being ignored, the axisymmetric, laminar liquid flow around the bubble is governed by the steady incompressible Navier–Stokes equation system
2.1
2.2
where *R**e* denotes the Reynolds number defined as 2*ρ**U**R*/*μ*, ** I** the identity tensor, and superscript ‘T’ stands for the transpose. Here, the variables are non-dimensionalized with length measured in units of the tube radius

*R*, velocity

**in units of bubble’s velocity**

*v**U*relative to the tube wall, and pressure

*p*in units of

*μ*

*U*/

*R*.

The mathematical problem description is completed with the boundary conditions (cf. Feng 2009)
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
where *C**a* ≡ *μ**U*/*γ* is the capillary number, the local unit normal vector ** n** at the free surface points from the liquid into gas, the local unit tangent vector

**points in the direction of increasing**

*t**s*along the boundary and relates to

**in such a way that**

*n***×**

*n***=**

*T*

*e*_{θ}in the right-handed coordinate system (

*z*,

*r*,

*θ*). Here

*R*

_{b}denotes the volume-equivalent radius of the bubble, and

*e*_{z}and

*e*_{r}denote the unit vectors in the

*z*- and

*r*-directions, respectively. A cylindrical (

*z*,

*r*)-coordinate system is used with the

*z*-axis coinciding with the axis of symmetry (i.e. the centreline of the tube) and pointing in the same direction as the wall velocity relative to the bubble.

With the present mathematical formulation, the Reynolds number *R**e* and capillary number *C**a* are the two independent parameters that can be conveniently specified. (With given values of *R**e* and *C**a*, the value of Weber number *W**e* can be readily obtained as *R**e* *C**a*.) The value of *p*_{i} in equation (2.9) is solved as part of the solution, becoming another independent parameter associated with the mathematical system. But the value of *p*_{i} depends on the locations of the inlet and outlet boundaries; therefore, it cannot be easily connected to a practically measurable quantity. As a direct consequence of *p*_{i}, the average liquid flow velocity (relative to the tube wall),
2.11
relates to the (constant) liquid flow rate in the tube independent of the location of inlet and outlet boundaries. Hence as in Feng (2009), instead of *p*_{i} is used in the present work as an independent parameter for describing the driving force for the bubble motion. Increasing the value of corresponds to increasing the ‘extra pressure drop’ as referred to by previous authors (e.g. Ho & Leal 1975; Martinez & Udell 1989, 1990) when considering drops and bubbles moving in a tube with flowing liquids.

## 3. Computational method

As described by Feng (2007) and citations therewith, solutions of the present problem can be computed by discretizing the partial differential equation system (2.1)–(2.10) with the Galerkin method of weighted residuals using finite-element basis functions (Strang & Fix 1973; Kistler & Scriven 1983).2 In doing so, the problem domain is divided into a set of quadrilateral elements (see figure 1), with biquadratic basis functions being used for expanding the velocity field and linear discontinuous basis function for pressure. The positions of finite-element mesh nodes around the deformable bubble surface are determined by a pair of elliptic partial differential equations that are also discretized by the Galerkin finite-element method (Christodoulou & Scriven 1992; de Santos 1991).

Then, the set of nonlinear algebraic equations of Galerkin’s weighted residuals is simultaneously solved by Newton’s method of iterations (Ortega & Rheinoldt 1970). At each Newton iteration, the Jacobian matrix of sensitivities of residuals to unknowns is evaluated with the values of unknowns determined in the previous iteration. The resulting linear algebra system is then solved by direct factorization of the Jacobian matrix with a modified version of Hood’s frontal solver (Hood 1976). The iteration is continued until the *L*_{2} norm of residual vector becomes less than 10^{−8}.

Critical to Newton’s iteration is a sufficiently accurate initial iterate (or initial ‘guess’). For solving the Navier–Stokes equation system, a convenient initial iterate can be obtained by computing a case for the Stokes flow (*R**e* = 0) at a small value of capillary number, e.g. *C**a* = 0.01, which usually corresponds to a bubble shape well approximated by a sphere. Thereafter, steady solution families can be computed by simple continuation schemes such as zeroth or first-order continuation (Riks 1972; Keller 1977) either in *R**e* while holding *C**a* fixed, or in *C**a* at fixed *R**e* (cf. Feng 2007).

## 4. Results and discussion

As a convenient starting iterate, solutions at *R**e* = 0 can be computed relatively easily due to reduced nonlinearities. With increasing surface tension, i.e. reducing the value of *C**a*, bubbles of radius *R*_{b} < 1 approach spherical shape. By continuation in *C**a* towards larger values, significant bubble deformations can be obtained, especially for *C**a* > 1. Figure 2 shows streamlines around deformed bubbles of *R*_{b} = 1,0.75,0.5 and 0.25 for *R**e* = 0 and *C**a* = 1.

For *R*_{b} = 0.75, the bubble shape appears quite similar to a case presented by Martinez & Udell (1990) for a neutrally buoyant liquid drop of a radius 0.726 (in units of that of the tube) and viscosity 0.19 (in units of that of the surrounding liquid) at *C**a* ∼ 0.5. However, the present result at *C**a* = 1 does not show the dimple (‘negative curvature’) on the bubble tail like that by Martinez & Udell (1990) at the same value of *C**a*, although liquid density inside the drop is not expected to play a role here at *R**e* = 0. But, the similar dimple on the bubble tail can be observed in the case of *R*_{b} = 0.5 and *C**a* = 5 at *R**e* = 0 (cf. figure 4 in this work). It should be noted that Martinez & Udell (1990) had a ‘viscosity ratio’ of 0.19, whereas that in the present work equals zero.

Unlike the case of buoyancy-driven bubble motion (Feng 2008), the flow field relative to the bubble is intrinsically fore–aft asymmetric (having a parabolic profile in the far-field) even at *R**e* = 0. At a given value of *C**a* (e.g. *C**a* = 1 as in figure 2), the significance of bubble surface deformation decreases with the bubble (relative) size; a bubble of *R*_{b} = 0.25 exhibits a nearly spherical shape though with a fore–aft asymmetric streamline distribution. In fact, in the small region adjacent to the tube centreline the liquid flow is almost uniform with diminishing local velocity relative to the bubble (as indicated by ). Without significant relative flow around the bubble, there cannot be substantial hydrodynamic stresses to cause free surface deformation. But as a bubble becomes large enough, i.e. *R*_{b} > 0.75, its shape exhibits characteristics quite similar to that of a typical long bubble (cf. Feng 2009), especially for *C**a* ≥ 1, as also found with buoyancy-driven bubbles in tubes (Feng 2008).

Table 1 shows the computed values of , (the minimum *z*-coordinate on bubble surface), (the maximum *z*-coordinate on bubble surface), (the maximum *r*-coordinate on bubble surface) and (the maximum tangential velocity on bubble surface) for bubbles of various radius *R*_{b} at *R**e* = 0 and *C**a* = 1. For a bubble of *R*_{b} = 0.25 with negligible surface deformation, the computed value of becomes basically the same as the value 0.5004 as predicted by the analytical formula of Haberman & Sayre (1958) for a spherical bubble of *R*_{b} = 0.25 in a tube at *R**e* = 0. In a coordinate system moving with the bubble, the liquid flow velocity at the tube centreline (due to a pressure gradient) becomes zero when approaches 0.5. Thus, a bubble of sufficiently small *R*_{b} is basically carried by the flowing liquid at the same velocity as the local liquid flow along the tube centreline. In other words, the intensity of liquid flow relative to the bubble becomes negligible for very small bubbles moving along the tube centreline. This fact is consistent with the diminishing bubble deformation with shrinking bubble size, as seen in figure 2, despite that the value of *C**a* is not small. However, for a bubble of *R*_{b} = 0.25, there is a measurable free surface deformation in terms of the axis ratio
4.1
which takes a value of 0.9956. This indicates a prolate deformation (when *α* < 1). If *R*_{b} is increased to 0.5, the value of *α* becomes 0.9012 indicating a quite noticeable free-surface deformation that leads to a computed noticeably less than 0.5129 as predicted by Haberman & Sayre (1958) for a corresponding spherical bubble. As *C**a* is reduced to 0.1, a bubble of *R*_{b} = 0.5 approaches the spherical shape and the computed becomes 0.5150. Because the values of in table 1 are basically ≤0.5, no obvious recirculating flow appears in figure 2 ahead of the bubble nose, according to the reasoning of Taylor (1961) based on a kinematic consideration.

Noteworthy in table 1 is that at *C**a* = 1 and *R**e* = 0, bubbles of larger *R*_{b} moves at faster velocities (i.e. smaller ) than those of smaller *R*_{b}. The computed results over a wide range of *C**a* show that such a trend holds true for *C**a* ≥ 1 but reverses for *C**a* ≤ 0.7. As the value of *C**a* is reduced, the surface tension effect becomes stronger and the cross-sectional area of a relatively larger bubble increases. Hence, the bubble tends to become less streamline (or bullet) shaped, and therefore move slower. For *R*_{b} ≤ 0.5 at *R**e* = 0, the value of remains within 10 per cent of 0.5 as long as *C**a* ≤ 5.3

At *R**e* = 100 and *C**a* = 1, much more severe bubble deformations are shown in figure 3, especially for bubbles of *R*_{b} > 0.5. Recirculating eddies appear behind the tails of bubbles of *R*_{b} = 0.75 and 1 due to large curvatures formed at the bubble tail rims. Again, the significance of bubble deformation decreases with bubble size; a bubble of *R*_{b} = 0.25 exhibits a nearly spherical shape, as the local liquid flow relative to the bubble diminishes.

As in tables 1 and 2 shows the computed values of , , , , and for bubbles of various radius *R*_{b} at *R**e* = 100 and *C**a* = 1. Unlike the case for *R**e* = 0, larger bubbles at *R**e* = 100 and *C**a* = 1 move slower than the smaller ones. This trend seems to be maintained throughout all values of *C**a* according to computed results over a wide range of *C**a*. Not surprisingly, recirculating flows ahead of the bubble nose appear in figure 3 corresponding to the fact that all the cases in table 2 have . In contrast to the cases at *R**e* = 0 with prolate free-surface deformations, bubbles of *R*_{b} = 0.25 and 0.5 at *R**e* = 100 and *C**a* = 1 are deformed to slightly oblate shapes with *α* = 1.008 and 1.039, respectively.

For a small bubble of *R*_{b} = 0.5 at *C**a* = 5, the prolate free-surface deformations become quite noticeable even at *R**e* = 0 and are further enhanced with increasing *R**e* (cf. figure 4). At *R**e* ≥ 200, namely the bubble-size-based Reynolds number *R**e*_{b} ≡ *R*_{b} *R**e* ≥ 100, the free surface takes a shape quite similar to that of long bubbles at *R**e* ≥ 100 with a skirt and sizable recirculating eddies forming at the bubble tail (as in Feng 2009). The computational result for a large range of *C**a* and *R**e* indicate that the average liquid flow velocity for a bubble of *R*_{b} = 0.5 varies within ±10 per cent around 0.5, as long as *C**a* is not too large (e.g. →10 and greater).

If *R*_{b} is reduced to 0.25, the computed *α* would become 0.7935, 0.7971, 0.8373 and 0.8211, with computed as 0.4910, 0.4917, 0.4962 and 0.4993, respectively, at *R**e* = 0, 40, 400 and 800 when *C**a* = 10. The average liquid flow velocity for a bubble of *R*_{b} = 0.25 varies no more than ±2 per cent around 0.5. Increasing *C**a* from 10 to 20 at *R**e* = 800 leads to a much more enhanced bubble deformation with *α* = 0.7169; yet the value of (as part of the numerical solution) is still 0.4945—within 1 per cent from 0.5. This points to the fact that a small enough bubble is basically carried by the flowing liquid, moving at the same speed as the local liquid flow. The deformed bubble surface has little influence on the bubble’s moving speed along the tube centreline, unlike the buoyancy-driven small bubble motion as represented by the value of Froude number (cf. Feng 2008).

For a very small bubble of *R*_{b} = 0.1, noticeable free-surface deformations can only be observed at very large values of *C**a* as expected in view of the diminishing relative liquid flow around the bubble. For example, at *C**a* = 100, figure 5 shows quite significant prolate bubble deformations for various values of *R**e*. The elongation of a bubble of *R*_{b} = 0.1 even at *R**e* = 0 must be a consequence of the wall-induced shear component in the background liquid flow.4 Again for *R**e*_{b} ≥ 100, corresponding to *R**e* ≥ 1000, the free surface takes a shape quite similar to that of long bubbles at *R**e* ≥ 100 (when *C**a* is not small), with skirt and sizable recirculating eddies forming at the bubble tail. Not surprisingly, the average liquid flow velocity for a bubble of *R*_{b} = 0.1 only varies within ±1 per cent around 0.5, despite the variety of shapes that a bubble may exhibit.

As in the case of buoyancy-driven bubble motion in tube (cf. Feng 2008), a bubble of *R*_{b} = 1 behaves just like a typical long bubble. For example, the computed values of for *R*_{b} = 1 at *R**e* = 0 (100) are 0.9731 (0.9726), 0.9007 (0.9068), 0.6972 (0.7236) and 0.4754 (0.5401), for *C**a* = 0.001, 0.01, 0.1 and 1, as comparable to 0.9753 (0.9734), 0.8967 (0.8958), 0.6945 (0.7287) and 0.4754 (0.5402) computed for typical long bubbles with flowing liquid in tube (Feng 2009). Thus, bubbles of *R*_{b} ≥ 1 (corresponding to bubble volume ≥4*π*/3) in tubes can be considered as long bubbles whose basic behaviour is independent of bubble size. For long bubbles moving with flowing liquid in tubes, Feng (2009) found that the (non-dimensionalized) average liquid flow velocity can be described as a function of *C**a* and *R**e* in the logistic dose–response form, , with , and being functions of *R**e*.

But for small bubbles, especially for those with 0.25 < *R*_{b} < 0.95, their moving velocities with flowing liquid in tubes can be much more complicated as the values of *R**e* and *C**a* vary. As shown in figure 6 for *R**e* = 100, the value of typically increases with *C**a*, as the bubble becomes more oblate (*α* > 1), until a maximum is reached and then decreases thereafter with further increasing *C**a*, as the bubble becomes more elongated (*α* < 1). This type of bubble deformation behaviour is most significant when *R*_{b}∼0.75 and diminishes as the bubble size is reduced to *R*_{b}∼0.25 or enlarged to *R*_{b}∼0.95, and only appears when *R**e* is large. At *R**e* = 0, however, the axis ratio *α* monotonically decreases with *C**a* from unity, and so is the value of from the value for spherical bubble as predicted by Haberman & Sayre (1958) which can be fitted in the logistic dose–response function similar to that in Feng (2009).

Figure 7 shows the variation of bubble shape with streamlines for *R*_{b} = 0.75 at *R**e* = 100 as *C**a* increases from 0.001 to 0.1, 0.25 and 1. Although both the bubbles of *C**a* = 0.001 and 0.25 have the axis ratio *α* about unity (1.0025 and 1.0191), their actual shapes differ considerably. The one of *C**a* = 0.001 has a plain spherical shape, whereas that of *C**a* = 0.25 takes a bullet-like shape which becomes firmly established at *C**a* = 1. For the bubbles of *C**a* = 0.001 and 0.25, the values of also differ somewhat: one is 0.6107 and the other 0.6594 suggesting that the bubble of *C**a* = 0.25 moves relatively slower than that of *C**a* = 0.001 despite that they both have *α*∼1. At *C**a* = 0.1, the bubble takes about the most deformed oblate shape with *α* = 1.1596 and .

## 5. Concluding remarks

The general behaviour of a short bubble (*R*_{b} ≤ 1) moving in a tube with flowing liquid is analysed computationally using a Galerkin finite-element method with boundary-fitted mesh as an extension to the work by Feng (2009) for long bubbles. This computational approach enables an in-depth study of (small, short) axisymmetric bubbles of various sizes in steady states for wide ranges of *C**a* and *R**e*. Unlike its long bubble counterpart, a relatively small, short bubble exhibits significant behavioural dependence upon its size (or volume) for a given set of *R**e* and *C**a*.

At small values of

*R**e*, the velocity of a bubble increases with bubble size for large values of*C**a*but decreases with bubble size for small values of*C**a*. For example, at the limit of creeping flow*R**e*= 0 as shown in table 1, decreases as*R*_{b}increases for a given value of*C**a*≥ 1 but increases with*R*_{b}for*C**a*< 0.7.At large values of

*R**e*(e.g. at*R**e*= 100 as shown in table 2), a bubble of larger size appears to move at a slower velocity for a given value of*C**a*.At large values of

*R**e*(e.g. at*R**e*= 100 as shown in figure 6), a bubble of*R*_{b}≥ 0.75 seems to have independent of bubble size*R*_{b}when*C**a*> 0.1.At large values of

*R**e*(e.g. at*R**e*= 100 as shown in figures 6 and 7), a bubble of 0.25 <*R*_{b}< 0.95 can be deformed into an oblate shape (with axis ratio*α*> 1) at some intermediate values of*C**a*with a corresponding increase in the value of (namely, slower bubble moving velocity for larger cross-sectional area). But for bubbles of*R*_{b}< 0.25 and >0.95, their axis ratios are almost always less than unity indicating an elongated shape in the liquid flow direction.At large values of

*R**e*(e.g. at*R**e*= 100 as shown in figure 7), a bubble of 0.25 <*R*_{b}< 0.95 takes a nearly spherical shape at very small*C**a*with*α*∼1, and then can be deformed into a nearly bullet shape again with*α*∼1 after the oblate deformation diminishes with increasing*C**a*. Here the value of axis ratio*α*does not convey much detailed information about the bubble shape (as also reflected by the difference in the value of ).Very small bubbles (e.g.

*R*_{b}= 0.1) can still be significantly deformed even at*R**e*= 0 (cf. figure 5) provided that the value of*C**a*is adequately large, despite that the relative local liquid flow around the bubble is diminishing.Despite significant shape deformations that may still occur (e.g. as shown in figure 5 at very large values of

*R**e*and*C**a*), bubbles of*R*_{b}< 0.25 (namely, the volume-equivalent radius less than a quarter of that of the tube) almost always move at the same velocity as that of the local liquid flow at the tube centreline (i.e. twice that of the average liquid velocity corresponding to ), regardless the values of*R**e*and*C**a*. This fact suggests that very small bubbles are basically carried by the local liquid flow with the relative liquid flow around the bubble diminishing.

## Footnotes

↵1 However, there were many publications on motion of liquid drops in a tube of comparable size with a flowing immiscible liquid, probably for the ease of minimizing the fluid density difference in an earthbound laboratory (e.g. Ho & Leal 1975; Olbricht & Leal 1982; Martinez & Udell 1990; Olbricht & Kung 1992; Borhan & Pallinti 1998, as well as citations therewith), often focusing on the low Reynolds number cases.

↵2 The actual computational code is called FECAW—finite-element computational analysis widget—as was used for computing the results in Feng (2007, 2008, 2009) and can be acquired from the website http://james.q.feng.googlepages.com/FECAWwelcome.html.

↵3 But at

*C**a*= 10 and 20 for*R*_{b}= 0.5, decreases to 0.4440 and 0.4384, respectively. In computing the creeping motion of neutrally buoyant drops through circular tubes, however, Martinez & Udell (1990) stated that ‘drops with radius ratios (corresponding to*R*_{b}here) less than 0.7 are insensitive to substantial variation in capillary number and viscosity ratio, etc.’ without clearly mentioning the range of parameter values.↵4 According to the theoretical analysis of (Saito 1913), a deformable bubble (or drop) in a uniform background flow field (in absence of wall effect) should retain a spherical shape at

*R**e*= 0 for arbitrary value of*C**a*.- Received May 28, 2009.
- Accepted October 2, 2009.

- © 2009 The Royal Society