Undulations on the surface of elongated bubbles in conﬁned gas-liquid ﬂows

A systematic analysis is presented of the undulations appearing on the surface of long bubbles in conﬁned gas-liquid ﬂows. CFD simulations of the ﬂow are performed with a self-improved version of the open-source solver ESI OpenFOAM (release 2.3.1), for Ca = 0 . 002 − 0 . 1 and Re = 0 . 1 − 1000, where Ca = µU/σ and Re = 2 ρUR/µ , with µ and ρ being, respectively, the viscosity and density of the liquid, σ the surface tension, U the bubble velocity and R the tube radius. A model, based on an extension of the classical axisymmetric Bretherton theory, accounting for inertia and for the curvature of the tube’s wall, is adopted to better understand the CFD results. The thickness of the liquid ﬁlm, and the wavelength and decay rate of the undulations extracted from the CFD simulations, agree well with those obtained with the theoretical model. Inertial eﬀects appear when the Weber number of the ﬂow We = Ca Re = O (10 − 1 ) and are manifest by a larger number of undulation crests that become evident on the surface of the rear meniscus of the bubble. This study demonstrates that the necessary bubble length for a ﬂat liquid ﬁlm region to exist between the rear and front menisci rapidly increases above 10 R when Ca > 0 . 01 and the value of the Reynolds number approaches 1000.


I. INTRODUCTION
The steady propagation of a long gas bubble through a viscous liquid within a circular tube is a classical problem in fluid mechanics, since the pioneering studies of Bretherton [1] and Taylor [2], and its interest has been motivated by applications such as flow in porous media, coating processes, and enhanced oil recovery. Under the dominant effect of capillary forces, the elongated bubble traps a thin film of liquid against the channel wall, and the front and rear menisci are separated by a uniform film thickness region; see the schematic in Fig. 1(a). Bretherton [1] utilized a lubrication theory to analyze the flow in the liquid film and concluded that, in the limit that Ca = µU/σ 1 and Re = 2ρU R/µ 1 (with µ and ρ being, respectively, the viscosity and density of the liquid, σ the surface tension, U the bubble velocity and R the tube radius), the thickness of the liquid film h 0 scales as h 0 /R ≈ Ca 2/3 . Taylor [2] performed experimental liquid film thickness measurements for very viscous fluids and showed that h 0 /R reached an asymptotic value of about 1/3 when Ca → 2. Ever since, the dynamics of the elongated bubble, in particular of its front meniscus, have been investigated extensively, e.g. effects of inertia [3][4][5][6][7], profiles of the bubble [6,8,9], streamlines ahead of the bubble tip [4,10] and pressure drop [4,8,11,12]. A summary of the studies that are relevant to the present work (Re 1) is provided in Table I.
While the literature concerned with the fluid dynamics of the front meniscus of a bubble in confined geometries is relatively vast, the rear meniscus has received far less attention; see Table I. The liquid film thickness decreases monotonically from the bubble nose to the uniform film region, whilst the gas-liquid interface presents stationary undulations near the rear meniscus of the bubble (see Fig. 1(a)). In the limit that Ca 1 and Re 1, the theoretical analysis of Bretherton [1] indicated that oscillations induce a minimum film thickness of 0.716h 0 , while their wavelength scales as h 0 Ca −1/3 [13]. However, the results of numerical simulations in the Re 1 regime showed significant deviations from the theoretical predictions of minimum film thickness and wavelength already at Ca ∼ 10 −3 − 10 −2 [3,14]. Furthermore, both experimental [9] and numerical simulations [3,9,14] emphasized that, when Re 1, the amplitude of the undulations increases, such that they become more apparent, and their wavelength decreases (for fixed Ca).
In summary, it has been observed that the characteristic features of the undulations on the surface of elongated, confined bubbles depend on both Ca and Re, but the literature still lacks Ca ≤ 2 Lubrication theory including inertia Data of [2] and [5] Re 10 3 applied to the flow in the liquid film well predicted by the model  [3] Ca ≤ 0.06 Numerical simulations The measured wavelength is lower than Re 10 3 of the entire bubble predicted by [13]; undulations become more apparent as Re increases Giavedoni and Saita [14] Ca 1 Numerical simulations Same as [3]; at fixed Ca, 0 < Re ≤ 70 of the rear meniscus wavelength decreases as Re increases a systematic analysis under a wide range of capillary and Reynolds numbers. Furthermore, no theoretical model is yet available to satisfactorily predict their geometrical parameters when Re 1. This topic is now of interest to a series of novel engineering applications and scientific problems where the topology of the ripples and their minimum film thickness are more important than the thickness of the uniform film region. For instance, in the removal of colloids from a microchannel by a translating gas-liquid interface, the local thickness of the liquid layer determines the surface cleaning efficiency [15]. In microchannel flow boiling, peaks of the heat transfer performance are observed upon the crests of these undulations [16] but, on the other hand, large amplitude oscillations may trigger local liquid film dryout [17,18], which has a negative impact on the performance and life of a micro-evaporator.
The objective of the present work is to perform a systematic analysis of the undulations appearing on the surface of long bubbles in confined gas-liquid flows, in the range Ca = 0.002 − 0.1 and Re = 0.1 − 1000. CFD simulations based on a Volume Of Fluid (VOF) method [19], and a theoretical model based on the work of de Ryck [6], are utilized to generate numerical and theoretical bubble profiles. The parameters characterizing the undulations in both the CFD simulations and the analytical model results are compared, and their trends versus the governing nondimensional groups are investigated. In particular, we show that in the visco-inertial regime the Weber number of the flow, We = Ca Re = 2ρU 2 R/σ, which represents the ratio of inertial stresses to capillary stresses, organizes the basic features of the bubble shape and thin film profile. This article is organized as follows: in Section II, a description of the problem is presented; the CFD model is introduced in Section III; the theoretical model is described in Section IV; finally, the results of the analysis are discussed in Section V.

II. DESCRIPTION OF THE PROBLEM
An elongated gas bubble translates at a velocity U in a channel of circular cross-section (radius R) filled with liquid. Figure 1(a) shows a sketch of the flow configuration under analysis, in a reference frame attached to the gas bubble. The flow field and the bubble profile are assumed axisymmetric. In order to describe the thin liquid film, and by analogy with the notation of Bretherton [1], the origin of the vertical axis y is on the tube's wall and its direction is radially inward. The horizontal coordinate x is directed downstream, and the reference x = 0 will be changed, as convenient, during this work. In these coordinates, the axial and radial velocities are defined as u and v, respectively, with v being positive when directed towards the axis of the channel. The walls of the tube have velocity u = −U . The thickness of the liquid film surrounding the bubble is h(x).
Following the notation and description of Bretherton, we assume that there exists a region of uniform film thickness, region CD in Fig. 1(a), where the liquid is moving with a uniform velocity −U and the film has thickness h 0 ; this feature will be confirmed by numerical simulations. The thin film region is axially bounded by two dynamic menisci Figure 1. (a) Sketch of a confined elongated bubble flowing within a tube and notation used in this work. Regions AB and EF represent the front and rear static menisci, respectively. Following  Fig. 1(a).
We assume that the flow is laminar, steady and incompressible, and that the fluid is Newtonian. The flow in the liquid film is governed by the steady-state Navier-Stokes equations, here written in nondimensional form: where the hats indicate dimensionless entities. The characteristic length scale chosen in Eqs. (1) and (2) At the gas-liquid interface, the boundary conditions are obtained from the dynamic balance of the stresses, here expressed by neglecting the viscous stresses within the gas phase: where n and κ indicate the interface normal vector and its dimensionless curvature, respectively. The capillary number is defined as Ca = µU/σ, where σ indicates the surface tension. Figure 1(b) illustrates the procedure adopted to characterize the decay and wavelength of the ripple appearing at the rear meniscus of the bubble. These details will be described in Section V B. In the figure, the liquid film thickness is rescaled by the uniform liquid film thickness, Y = h/h 0 , and the horizontal coordinate is rescaled by a length . The length is on the order of the length of the dynamic menisci, BC and DE in Fig. 1(a), and will be identified below when the lubrication analysis is presented.

III. CFD SIMULATIONS
Numerical simulations of the flow configuration reported in Fig. 1 [20]. In order to improve the accuracy of the surface tension representation with respect to the original interFoam solver, the interface normal vector and curvature are estimated from a smoothed VOF function field [21] when Ca > 0.005 and from a reconstructed level set function [22] when Ca ≤ 0.005. Details of the flow equations solved by interFoam and of the numerical methods adopted can be found in Ref. [9]. The computational mesh is made of uniform square cells, which are gradually refined at the wall boundary in order to capture the fluid mechanics within the thin liquid film surrounding the bubble. To identify parameters for the simulations a grid independence test was run by utilizing the same mesh arrangements described by Khodaparast et al. [9]. It was observed that the main parameters characterizing the bubble profile (uniform liquid film thickness, wavelength of the undulations and minimum film thickness) for a computational mesh including 40 square cells of uniform size along the radial direction, followed by 10 gradually refined mesh elements in the near wall region, deviated by less than 2 % from those measured with more refined grids. This is the computational grid utilized in the present work (see [9] for more details).
The present numerical model was originally validated by Khodaparast et al.  Re=1000, and different bubble lengths L b . x F locates the intersection between the tube axis and the rear meniscus of the bubble (see Fig. 1(a)).
where Y = h/h 0 and the subscript indicates a third-order derivative with respect to the nondimensional variable X. The length scale = h 0 Ca −1/3 has been used to define X = x/ . Figure 2(a) shows that for Re = 1 the profiles at the bubble back obtained with the CFD simulations agree reasonably well with Bretherton's theory, and the flat film region is recovered after two undulation crests. A deviation of the uniform film thickness of 9.2 % with Bretherton's theory is measured, with h 0 /R = 0.0195 in the simulations versus h 0 /R = 1.34Ca 2/3 = 0.0213. This difference can be ascribed to the limited accuracy of Bretherton's law when Ca > 0.001. It is observed that, for Ca = 0.002, the empirical fit of Aussillous and Quéré [5] to Taylor's data predicts a film thickness of h 0 /R = 1.34Ca 2/3 /(1 + 3.35Ca 2/3 ) = 0.0202, with only a 3.6 % deviation from the present CFD data. When Re = 1000 ( Fig. 2(b)), the liquid film thickness increases as reported by many authors [5,6,8], the transition region for the rear meniscus becomes longer, and the undulations become more apparent. Their amplitude increases or, as will be discussed below, their decay rate decreases, and three peaks become clearly visible in Fig. 2 there is no appreciable effect of Re and the bubble profiles collapse onto one shape. In these conditions, the shapes of the rear menisci agree well with that calculated according to Bretherton's low Ca and low Re theory. As the Reynolds number is increased above 100, the profiles are modified. The minimum film thickness (measured at the most upstream crest) stays approximately constant, while the amplitude of the undulations increases and more crests appear. Note also that the wavelength, estimated as the horizontal distance between two consecutive peaks, decreases slightly. All of these parameters characterizing the undulations will be studied in detail in Section V by systematically varying the capillary and Reynolds numbers. x F locates the intersection between the tube axis and the rear meniscus of the bubble (see Fig. 1(a)).

ING INERTIAL EFFECTS
The model used in this work to derive theoretical profiles of the rear meniscus of the bubble was originally presented by de Ryck [6], who applied it to investigate the front meniscus of the bubble. Before introducing the complete model, simpler analytical models obtained by increasing levels of complexity starting from Bretherton's theory are briefly reviewed below.
The equations governing the fluid flow in the dynamic menisci regions for the axisymmetric flow configuration depicted in Fig. 1 with u y = 0, and p = −σκ, where the pressure inside the bubble is assumed constant and is set to zero. Under a small slope approximation, |h x | 1, By setting = h 0 Ca −1/3 , Bretherton obtained a third-order differential equation for the nondimensional profile Y (X), as shown in Eq. (5), which holds as long as δ = h 0 / = Ca 1/3 1 and Re 1.
The effects of inertia in the dynamic menisci regions (BC and DE in Fig. 1(a)) can be introduced in the equations above by retaining the convective term in the x−momentum The continuity and momentum equations along x are now replaced by equations integrated along the y coordinate from y = 0 to y = h(x): where the boundary condition at the free surface Eq. (7) has been used to write p x = −σκ x .
These equations are still difficult to solve because u and v depend on both x and y. A method for going further has been proposed by Shkadov [24] and has been afterwards adopted by Esmail and Hummel [25] and Koulago et al. [26] to solve fluid mechanics problems associated with thin liquid film flows in the presence of inertia. The velocity profile is supposed to remain parabolic, as in the purely viscous case, but the effective pressure gradient that induces the flow is taken to be unknown. The velocity profile must satisfy the no-slip condition at the wall and Eq. (7) at the free surface, so that u(x, y) has the form: where the term between square brackets comes from the exact solution in the purely viscous case. Equation (11) which yields: These steps lead to a nonlinear differential equation for h(x) which, by adopting the dimensionless coordinates Y = h/h 0 and X = x/(h 0 Ca −1/3 ), takes the nondimensional form where H = h 0 /R is the nondimensional uniform film thickness. The Reynolds number is still defined as Re = 2ρU R/µ. The second term on the right-hand side of Eq. (15) where Bretherton's law H ≈ Ca 2/3 has been used to express H.
This theoretical model is regarded as approximate, because inertia in the x−momentum equation is the only first order term added to the leading order Eqs. (6) and (7), while other first order terms, that will be included in the next models developed below, are here with u y = 0, and p = −σκ + 2µv y , where the µu xx term is included, although being O(δ 2 ), in order to avoid reduction of the order of the boundary condition at the free surface. Integrating Eq. (17) where the last term on the right-hand side arises from the more accurate boundary condition at the free surface. Details of the derivation are given in the Supplemental Material. A similar equation was derived by Esmail and Hummel [25] to describe the film entrained from a bath of liquid by the vertical withdrawal of a sheet at high velocity.
Due to the small slope assumption utilized to express the interface curvature, Eq. (19) is expected to work well as long as |h x | 1. Such a hypothesis can be relaxed by considering the complete expression for the interface curvature, When introducing Eq. (20) into the free surface boundary condition (18), the final differential equation for the film thickness becomes where and the term between braces coincides with the right-hand side of Eq. (19). The first and second terms on the right-hand side of Eq. (21) implement the effects of the slope of the interface, and hence Eq. (21) reduces to Eq. (19) when |Y X | 1.
Finally, de Ryck [6,27] introduced the curvature of the tube's wall in the formulation of the problem by retaining the related O(δ) term within the x−momentum equation, The differential equation governing the dynamic meniscus profile then becomes with Note that the terms between braces are formally similar to those included in Eq.
This algebraic equation has one real and two complex conjugate roots, and therefore the general solution of the linear problem is where β, α and λ derive from the solution of Eq. (26) and are functions of both Ca and Re.
The procedure for the numerical integration of Eq. (24) at the front and rear menisci is illustrated in the following sections.

A. Numerical integration at the front meniscus
The integration for the front meniscus starts from close to the uniform film region (point C in Fig. 1(a)) and proceeds towards X → +∞. The linearized initial conditions are: where in our calculation we take = 10 −4 .
Compared to the Bretherton problem governed by Eq. (5), the additional terms appearing in the theoretical model used here introduce two difficulties. First, the differential equation for the liquid film profile now depends on the nondimensional uniform film thickness H.
Second, integration of Eq. (24) shows that Y X and Y XX quickly diverge towards very large values as Y increases. Therefore, the asymptotic matching adopted by Bretherton [1] cannot be used. The procedure to match the film profile in the dynamic region to a static meniscus of spherical shape (point B in Fig. 1(a)) is as follows. For a given set of Ca and Re, there exists a range of values of H, 0 < H ≤ H * , that allow matching the meniscus to a sphere.
The critical value H * is the largest H allowing the matching. As an example, Fig. 4 Figure 4(b) shows the nondimensional curvature of these interface profiles, κ = κ/(h 0 / 2 ), calculated as and that of the matching spheres, κ sph , calculated as It can be seen that, for each value of H < H * = 0.0642, there is always one point where κ = κ sph . This represents the curvature matching point, point B in Fig. 1(a). As H is  Fig. 1(a)).   bubble length in all of the simulations is about 20R. It can be seen that the values of the maximum nondimensional film thickness H that allow matching the bubble nose to a sphere extracted from the theoretical model are in very good agreement with the CFD results, with an average deviation of 3.9 %. Figure 5 also shows the predictions given by the empirical correlations of Han and Shikazono [7] and Langewisch and Buongiorno [11], which exhibit deviations from the present CFD data of, respectively, 4.2 % and 8.5 %.
Although these correlations were obtained by fitting large experimental [7] and numerical [11] databases, they do not perform better than the present theoretical model when applied to predict the CFD data. The three simplified theoretical models based on Eqs. (15), (19) and (21)  It is worth noting that the theoretical model, the CFD results and the empirical correlations all show that, for a fixed Ca, the uniform film thickness exhibits a decreasing trend with Re for Re 10 2 , followed by a monotonically increasing trend at larger Re numbers (Fig. 5). This behavior is observed also for the simpler theoretical models, e.g. Eq. (15) which includes only the viscous and inertial terms, and it originates from the expression of the inertial term. For Re 1, the inertial term gives a negligible contribution to Y XXX . As Re is increased above 1, the inertial term gives an increasing negative contribution to Y XXX for Y < √ 6 (see Eq. (15)), and an increasing positive contribution to Y XXX for Y > √ 6.
The overall behavior of Y XXX from Y = 1 to Y = H at each value of Re changes in such a way that, as Re is increased in the Re This mixed effect of inertia on the film thickness can be also derived by a scaling analysis of the forces acting on the front meniscus. Aussillous and Quéré [5] expressed the curvature of the bubble nose as where the −h 0 in the denominator of the second term accounts for the presence of the liquid film which increases the curvature of the bubble nose. However, inertia further increases the bubble nose curvature as it makes the bubble nose more slender. This effect can be accounted for by modifying Eq. (32) as where f κ > 0 increases with both Ca and Re [3,8,9]. Balancing the viscous force with the inertial force and the pressure gradient in the equation of motion along the dynamic meniscus yields while the matching between both static and dynamic menisci suggests Using Eq. (35) to replace in Eq. (34), a scaling law for the liquid film thickness is obtained: where we choose γ = 10 −6 . The nondimensional uniform film thickness H is set to the value extracted from the numerical integration at the front meniscus. The solution now depends on the parameter ϕ, which shifts horizontally the profile characterizing the initial condition.
Similarly to what was done for the front meniscus, here ϕ is set to the value that allows matching the profile of the bubble back to a sphere (point E in Fig. 1(a)) at the farthest distance from the tube's wall. It will be shown in the next sections that this criterion yields theoretical profiles of the rear meniscus that compare quite well with the CFD results. Next, a low capillary (Ca = 0.002) and high Reynolds number (Re = 1000) case is shown in Fig. 6(b). The presence of two peaks and two valleys is now apparent. The two methodologies give very similar features of the undulations.
A relatively larger capillary (Ca = 0.05) and Weber (We = 5) number case is presented in Fig. 6(c). The case discussed has Ca 2/3 ≈ 0.14 (model assumptions are valid as long as  Fig. 1(a)).
number, Ca = 0.005, as illustrated in Fig. 6(d). The main characteristics of the ripples when increasing Re are all well captured by the model. Similarly to Fig. 3, the liquid film profiles begin to change when Re > 100, i.e. We 0.5. The liquid film becomes thicker, and more undulation crests become apparent, while the minimum film thickness stays almost constant.

B. Systematic analysis of the undulations
The profiles of the rear meniscus displayed in Fig. 6 indicate that, when moving towards X → +∞, the undulation maintains a sinusoidal shape while its amplitude is gradually reduced. As such, its profile can be well approximated as a cosine function multiplied by an exponential damping term, analogous to the linear solution of Eq. (24) at the bubble back, where Y = h/h 0 and X = x/(h 0 Ca −1/3 ). Here, λ identifies the wavelength of the undulation, while α represents its decay rate, both in nondimensional units; α > 0 in order for the amplitude to decrease in the direction of increasing X. Estimates of λ and α as a function of the capillary and Reynolds numbers are already provided as roots of Eq. (26).
The wavelength and decay rate in actual rear meniscus shapes predicted by the CFD simulations and by the numerical integration of Eq. (24) are extracted as illustrated in Fig. 1(b). The ripple on the bubble surface is assumed to be well represented by Eq. (38).
(X 1 , Y 1 ) denote the nondimensional coordinates of the most upstream crest, which is also the largest amplitude. Also, (X 2 , Y 2 ) indicate the coordinates of the second crest along the positive X direction. The wavelength is calculated as the nondimensional distance between these two most upstream crests Note that in the dimensional units based on the lubrication description, this distance is where λ is the nondimensional wavelength given by Eq. (39). According to Eq. (38), the decay rate of the undulation is calculated by considering the decrease of its nondimensional amplitude between the first and second most upstream crests where basic results of the Bretherton problem have been utilized; the capillary effects enter in the expression above because they determine the thickness of the liquid layer, i.e. the characteristic length along which viscous effects act. It is common to interpret We as the ratio of inertial stresses to capillary stresses. Nevertheless, as Eq. (41) suggests, for thin film flow where the local pressure gradient is set by capillary effects, as in the Bretherton problem, then the ratio of inertial effects to viscous effects is also set by the Weber number.
Now we report the wavelength of the ripples as a function of the Weber number in  Fig. 7(a)). It is interesting to note that the wavelength obtained from the linear model is always very close to the values extracted from the complete theoretical model, despite being slightly higher when We < 1.
As general trends, when We < 10 −1 the results collapse onto a single curve for Ca up to 0.05 and are independent of the Weber number. This suggests that inertial effects are negligible for this parameter range. Within this range, the nondimensional wavelength has the value λ = 4.7 − 4.8. In the asymptotic limit of Ca 1/3 1, Bretherton [1] obtained the value of λ = 5.03 from the linear solution of Eq. (5), and a value of λ = 4.82 is extracted from the profile obtained from the numerical integration of Eq. (5).
When We > 10 −1 the inertial effects manifest with a decrease of the wavelength as the Weber number increases, which is in line with the results obtained by Giavedoni and Saita [14] and Edvinsson and Irandoust [3] by means of numerical simulations, while the capillary number has a negligible effect. de Ryck and Quéré [28] observed that in fibre coating the thickness of the liquid film began to deviate from the low Reynolds number theory when We ≈ 1, which is the value typically assumed in the practice to judge the importance of inertial effects. The present results suggest that inertial effects begin to appear when We > 10 −1 , and become evident when We 1.
When We > 10, time-dependent patterns appear in the CFD simulations analogously to what was observed both experimentally and numerically by Khodaparast et al. [9]. The wavelength becomes a function of time and the values reported in Fig. 7(a) have to be interpreted as time-averaged quantities. We ≈ 10 represent a limiting condition also for the theoretical model. Above this critical value (actually for values of We that slightly decrease when increasing Ca), the numerical integration of Eq. (24) towards the bubble back yields oscillations that grow unbounded as X → −∞, thus making it impossible to match the profile to a sphere.
The decay rate of the ripples as a function of the Weber number is presented in Fig. 7(b), for different values of Ca. It can be seen that the capillary number is now an influential parameter, as larger values increase the decay thus making the amplitude of the ripples decrease rapidly. Therefore, the oscillation tends to disappear after the first crest and valley, and the existence of a smooth and long liquid film region around the bubble is promoted. This trend is in agreement with results reported by Giavedoni and Saita [14]. Inertial effects make the trends dependent on the Weber number when We > 10 −1 . The spatial decay of the undulations decreases, and therefore more crests and valleys become visible, as was already observed in previous works [9,29] and in agreement with the trends in Fig. 3.
When time-dependent patterns occur in the CFD simulations (We ≥ 10), the timeaveraged values of α reported in Fig. 7(b) are not really meaningful because they are characterized by a very large standard deviation, which is on the order of the mean value itself.
Furthermore, the amplitude of the second crest can become larger than that of the first crest at certain time instants, thus making α negative. The analysis of this time-dependent regime for the flow of the bubbles is outside of the scope of the present paper.
The average deviation between the CFD results and theoretical model based on Eq. (24) for all of the data points included in Fig. 7 (20), yields a significant improvement of the description of the liquid film profile.

C. Minimum bubble length
It has been observed both in Fig. 3 and Fig. 6 The inset in Fig. 8(a) shows a sketch of the procedure for the calculation of the length at the front meniscus. This length is evaluated as the axial distance between the tip of the bubble nose (point A) and the location where the liquid film thickness h decreases below 1.01 times the uniform film thickness h 0 . The inset in Fig. 8(b) illustrates the procedure for the bubble back. The length of the dynamic region of the rear meniscus is taken as the axial distance between the tip of the bubble back (point F) and the point where the liquid film thickness deviates from h 0 by less than 1 %. The 1 % threshold for both nose and tail is a small value that was chosen arbitrarily.
The plots of the front meniscus and of the entire bubble (front plus rear) length are given in Fig. 8. The curves are presented versus the Reynolds number because this is still more frequently used than the Weber number when reporting results for confined bubble flows in narrow channels. When inertial effects are absent, the transition length is a function of the capillary number only and it increases with Ca. As Ca → 0, L f /R → 1 and L b /R → 2 (L b = L f + L r ), i.e. the bubble exhibits spherical caps. In these conditions, the transition region at the bubble back is slightly longer than that at the front due to the presence of an undulation at the matching point between the dynamic and static menisci. The minimum bubble length in this visco-capillary regime is always relatively short, below 4R up to Ca = 0.1.
Inertial effects appear at decreasing values of the Reynolds number when increasing Ca (at We ≈ 10 −1 ). As is evident in Fig. 8(b), the necessary length increases steeply with Re, in particular at high Ca, and it grows above 10R when the Reynolds number approaches 1000.
Note that curves for Ca ≥ 0.02 are interrupted when the matching of the rear meniscus to a sphere is no longer possible, which also corresponds to the onset of time-dependent patterns in the CFD simulations. Air bubbles flowing in water within a tube of diameter of 0.5 mm at Ca = 0.02 have Re of about 900. Figure 8(b) demonstrates that bubbles with length 6R, such as those used in the work of Khodaparast et al. [9], are not sufficient to recover a flat film region under such conditions.
The steps in the curves in Fig. 8(b) at high Re reveal the appearance of new crests on the surface of the bubble at the rear meniscus, i.e. oscillations whose amplitude has grown above 0.01h 0 . It can be seen that a second crest appears at Re = 35 − 60, which corresponds to Weber numbers that range from 0.12 (at Ca = 0.002) to 4 (Ca = 0.1). The appearance of the third crest is anticipated as the capillary number is increased, which happens for We ranging from 1.42 (Ca = 0.002) to 4.5 (Ca = 0.05).

VI. CONCLUSIONS
In this study we present a systematic analysis of the undulations appearing on the surface The equations governing the flow in the liquid film region surrounding the bubble, in the cylindrical axisymmetric coordinates introduced in Fig. 1(a), are written: ρ (uu x + vu y ) = −p x + µ u xx + u yy − u y R − y , (A2) When neglecting the viscous stress within the gas phase, and setting the pressure within the bubble to a zero reference value, the boundary conditions at the free interface are and p + σκ + µh x (u y + v x ) − 2µv y = 0, at y = h(x).
The normal vector at the interface has components and the interface curvature is expressed as shown in Eq. (20).