Two-dimensional global stability analysis of elongated bubbles moving in a horizontal tube

PHYSICAL REVIEW FLUIDS 10, 053603 (2025) Two-dimensional global stability analysis of elongated bubbles moving in a horizontal tube Mirco Magnini 1and Miguel A. Herrada 2 1Department of Mechanical, Materials and Manufacturing Engineering, University of Nottingham, NG7 2RD, United Kingdom 2Dept. Ing. Aerospacial y Mecánica de Fluidos, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain (Received 2 December 2024; accepted 28 March 2025; published 7 May 2025) The linear stability of an elongated axisymmetric gas bubble transported by a liquid in a capillary tube is analyzed through the use of numerical simulations. The study focuses on the influence of inertia, characterized by the Reynolds number (Re) and the imposed flow rate, characterized by a capillary number (Cal), on the stability of the bubble tail, which exhibits ripples at high Re. The numerical approach utilizes a boundary-fitted method in a reference frame anchored to the bubble, combined with a mixed spatial discretization based on spectral collocation in the radial direction and finite differences in the axial direction. This framework enables the computation of steady nonlinear solutions through a Newton iteration scheme and facilitates the linear stability analysis of these solutions. Direct numerical simulations using the volume of fluid method in OpenFOAM are also performed to corroborate the results of the stability analysis. We perform systematic simulations for Cal=0.005 −0.04 and observe that the system becomes unstable, with the emergence of oscillations at the rear of the bubble, when the Reynolds number grows above a critical value, designated as Re∗; this critical value is dependent on the capillary number. The instability is due to the increased inertia of the recirculating flow in the liquid behind the bubble, which impinges its rear meniscus. A modified Weber number Wep, based on the relative velocity between the external flow and the bubble, is introduced to describe the competition between the destabilizing pressure force acting on the bubble rear and surface tension. Our results show that the bubble dynamics become unstable for a critical value, We∗ p≈3.65, which remains quite uniform across the range of capillary numbers tested, and divides the Cal−Re diagram into stable and unstable regimes. Our findings offer insights into the behavior of bubbles in microfluidic applications, with implications for heat transfer, mass transfer, and cleaning processes in microchannels. DOI: 10.1103/PhysRevFluids.10.053603 I. INTRODUCTION When a long gas bubble propagates within a circular capillary tube filled with a wetting liquid, it traps a thin liquid film against the channel wall which, sufficiently far from the front and rear menisci, [1] takes a constant thickness. The thickness of this film depends on the capillary number, [2,3]Ca=μlU/γ , with μlbeing the liquid viscosity, Uthe bubble or liquid average speed, and γ the surface tension, which quantifies the importance of viscous effects over surface tension, and the Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. 2469-990X/2025/10(5)/053603(20) 053603-1 Published by the American Physical Society MIRCO MAGNINI AND MIGUEL A. HERRADA (a) (d) (e) (f) (g) (b) (c) FIG. 1. (a)–(c) Experimental [7] and (d)–(g) computational [9] snapshots of the rear of a long gas bubble traveling in a liquid-filled capillary for Ca =0.023 and Re =1065. Flow is from left to right. In the experiment, air and water were used as fluids, the channel was horizontal with a diameter of 0.5 mm, and the bubbles were imaged from the top of the channel with a frame rate of 104s−1. Three-dimensional simulations were run with a coupled volume of fluid and level set solver [9] in OpenFOAM. Reynolds number, [4,5]Re=2ρlUR/μl, with ρlbeing the liquid density and Rthe channel radius, which describes the competition of inertia and viscous forces. In the limit that Ca 1 and Re 1, the bubble travels steadily, the shape of both the front and rear menisci of the bubble approximate well that of a sphere of radius R, and an undulation appears on the liquid film at the matching point with the rear meniscus [2]. As both capillary and Reynolds numbers increase, the curvature of the bubble nose increases and that of the bubble tail decreases, hence the bubble tail flattens [6,7]. When the Weber number of the flow, We =Ca Re =2ρlU2R/γ , grows to above 0.1, inertial effects become important and manifest with the appearance of multiple (stationary) undulations on the bubble surface, at the matching point between the rear meniscus and the thin liquid film [1]. These undulations are already present in the Bretherton’s solution, however, for We <0.1 their amplitude decays quickly such that only one crest remains visible. As the Weber number is increased above 0.1, the decay rate of the undulations decreases such that more crests become visible and their amplitude increases [8]. Magnini et al. [1] reported that upon a further increase of the Weber number above a value of about 10, time-dependent patterns appeared on the rear meniscus, which started to oscillate, whereas the bubble nose continued traveling at a constant speed. A similar observation was reported by Khodaparast et al. [7], who performed experimental measurements of bubble shape, speed, film thickness, and front and rear menisci curvatures, for long air bubbles transported in water in a 0.5 mm horizontal circular channel. When the capillary number was below 0.01 the bubbles traveled steadily, but for higher capillary numbers (and Reynolds numbers approaching 103)the rear meniscus of the bubble was observed to flap significantly as the bubble progressed along the channel, see snapshots in Fig. 1, with a consequent impact on the surface undulations. Besides the qualitative observations outlined above, no systematic study has been yet attempted to characterize these unsteady patterns or explain their appearance. The flow of elongated bubbles in microchannels at high Weber number is now receiving increasing attention due to emerging applications utilizing low-viscosity fluids such as two-phase cooling in micro-evaporators [10,11], micro-chemical reactors [12], cleaning of bacterial cells from medical surfaces [13], geological CO2sequestration [14], or environmental problems such as detachment and transport of pollutants 053603-2 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … FIG. 2. Sketch of the flow geometry considered in this study: an axisymmetric bubble in a horizontal tube with an applied axial flow. in unsaturated soil [15]. Unsteady patterns near the bubble tail may have a detrimental effect on the applications introduced above, because they may promote draining and dewetting of the lubricating film, compromising the wall-fluid exchanges such as heat/mass transfer and cleaning efficiency. For example, Borhani, Agostini, and Thome [16] utilized high-speed imaging of the two-phase flow within a microchannel evaporator and reported that liquid film dryout was initiated from dry patches forming in the proximity of the bubble tail, due to a combined effect of hydrodynamic instabilities and liquid evaporation. In this work, the linear stability of elongated bubbles in horizontal tubes will be analyzed using a variant of the numerical technique developed by Herrada and Montanero [17] to explain the previously mentioned observations and identify the transition between steady and unsteady flow regimes. The results of the linear stability analysis will also be compared to those obtained with nonstationary numerical simulations using OpenFOAM. The remainder of the paper is structured as follows. In Sec. II, we present the governing equations and the numerical techniques employed to perform the computations of steady axisymmetric flows and stability analysis. In Sec. III, we present the details of the nonlinear unsteady simulations performed with OpenFOAM. In Sec. IV, we discuss the base state results and present the linear stability analyses, in addition to the comparisons with the results obtained with OpenFOAM. Finally, in Sec. V, we provide concluding remarks. II. MODEL FORMULATION We consider the configuration sketched in Fig. 2, where an elongated bubble of a gas of density ρg, viscosity μg, and volume Vis moving in a cylindrical horizontal capillary of radius Rcontaining a liquid of density ρland viscosity μl. We assume that the flow is axisymmetric and thus we adopt a two-dimensional axisymmetric model, the validity of which was verified by means of three-dimensional simulations run with OpenFOAM, as documented in the Appendix. To model the movement of the bubble under the action of an imposed flow rate through the tube Q, the system is solved in a frame moving with the bubble where a noninertial axisymmetric cylindrical coordinate system (r,z) aligned with the tube axis and anchored to the bubble is selected. This system is, in turn, linked to an inertial cylindrical coordinate system (ro,zo) anchored to the tube walls and also aligned with the tube axis. To model the apparent forces in the momentum equations in this noninertial coordinate system, we have to include all of the accelerations of the system. In our case, since both 053603-3 MIRCO MAGNINI AND MIGUEL A. HERRADA coordinate systems are aligned, we will only have to consider the axial acceleration a=dU dt eZof the reference frame, where U(t) is the axial velocity of the right point of the liquid-gas interface (green arrow in Fig. 2) and Ul=Q/(πR2) is the mean velocity of the liquid flowing through the tube (both velocities measured in the inertial coordinate system). Figure 2also indicates the computational domain used in the simulations (a cylinder of radius Rand length H) and the boundary conditions applied at the outer boundaries of this domain. A. Governing equations The conservation of mass and a balance of linear momentum in the liquid (i=l) and gas (i=g) subdomains is given by ∇·vi=0,(i=l,g),(1a) ρi∂vi ∂t+(vi·∇)vi=∇·σi,(i=l,g),(1b) where vi=wiez+uieris the velocity field and σiis the stress tensor of material i(i=g,l). We consider in both regions the incompressible flow of a Newtonian fluid, where the stress tensor takes the form σi=−piI+μi[∇vi+(∇vi)T],(i=l,g) (1c) where pi(i=l,g) is the reduced pressure, pi=p∗ i+ρidU dt z,p∗is the static pressure, and the other term arises from the potential associated with the linear acceleration of the frame of reference. At the left boundary, z=−2H/3, (see Fig. 2) we assume a parabolic axial velocity profile, wl=−U(t)+2Ul[1 −(r/R)2],ul=0.(1d) At the right boundary, z=H/3, the static pressure is set to zero, p∗ l=0.(1e) At the wall, r=R, no-slip boundary conditions are applied, wl=−U(t),ul=0.(1f) Across the gas-liquid interface parametrized in terms of a meridional arc length s(0 ⩽s⩽1) rI=F(s,t) and zI=G(s,t), we impose that the velocity field must be continuous, i.e., wl=wg,ul=ug,(1g) and impose a balance of normal and tangential stresses between the liquid and the gas in the form n·(σl−σg)·n=γκ +(ρl−ρg)dU dt zI,(1h) t·(σl−σg)·n=0,(1i) where n=Gser−Fsez G2 s+F2 s1/2,t=Gsez+Fser F2 s+G2 s1/2,(1j) are normal (n) and tangential vectors (t) to the surface, the subscripts srepresents derivatives with respect to s,κ=∇·nis twice the mean curvature, and γis the surface tension. In addition, the kinematic boundary condition on the interface is written ul−∂F ∂t∂G ∂s−wl−∂G ∂t∂F ∂s=0.(1k) 053603-4 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … To ensure a uniform distribution of points along the arc length s, the following equation is applied: ∂F ∂s ∂2F ∂s2+∂G ∂s ∂2G ∂s2=0.(1l) At the axis, r=0, regularity conditions are considered. Finally, to close the problem we need an additional condition for the bubble velocity U(t), which is unknown because the liquid flow rate Qis the controlled parameter in this work (and, usually, in experiments) and thus the bubble speed is a result of the simulation. This is achieved by imposing that a particular point of the interface has a zero axial velocity in the noninertial frame of reference. Therefore, we choose that the right point of the interface crossing the axis verifies wl=0,zI=G=0,rI=F=0,at s =1.(1m) The governing equations are integrated with a variant of the numerical method proposed by Herrada and Montanero [17] and Herrada [18]. The spatial domain occupied by the liquid and the gas are mapped into two rectangular domains [0 ⩽sl⩽1] ×[0 ⩽ηl⩽1] and [0 ⩽sg⩽1] ×[0 ⩽ηg⩽ 1] using quasi-elliptic transformations [19]. These mappings are applied to the governing equations (1) and the resulting equations are discretized in the η-direction with nηland nηgChebyshev spectral collocation points in the liquid and gas domains, respectively. Conversely, in the s-direction, we use second-order finite differences with nsl and nsequally spaced points in the liquid and gas domains, respectively. The results presented in this work were obtained using ns=851, nsl =1353, nηg=10, and nηl=11. B. Axisymmetric steady solutions The steady-state solutions of the nonlinear equations (1) with all variables independent of time are obtained by solving all discretized equations simultaneously (a so-called monolithic scheme) using a Newton technique. We denote the steady axisymmetric solution of the system with the subscript b. The steady bubble profile hb(s), is defined as hb(s)=R−Fb(s),0⩽s⩽1.(2) As described in [20], hbfor a long bubble without inertia (Re =0) is characterized by three distinct regions: (I) the bubble “nose” hb= cte; (II), uniform film region hb(s)=cte =b; and (III), the bubble “tail” hb= cte (see Fig. 4). We trace the steady solutions as a function of the model parameters and quantify them using the steady axial bubble length, Lzb, defined as Lzb =Gb(s=1) −Gb(s=0).(3) C. Small amplitude 2D perturbations To test the stability of a given steady axisymmetric state we calculate the linear global modes by assuming the temporal dependence (r,z;t)=b(r,z)+δ(r,z)e−iωt(1),(4) where (r,z;t) represents any dependent variable while b(r,z) and δ(r,z) denote the base (steady) solution and the spatial dependence of the eigenmode for that variable, respectively, while ω=ωr+iωiis the frequency (an eigenvalue). Both the eigenfrequencies and the corresponding eigenmodes are calculated as a function of the governing parameters. The dominant eigenmode is that with the largest growth factor ωi, so if it is positive, the base flow is asymptotically unstable. The numerical procedure used to solve the steady problem can be easily adapted to solve the eigenvalue problem that determines the linear global modes of the system. In this case, the temporal 053603-5 MIRCO MAGNINI AND MIGUEL A. HERRADA derivative is computed assuming the dependence (4). The spatial dependence of the linear perturbation δ(q)is the solution to the generalized eigenvalue problem J(p,q) bδ(q)=iωQ(p) bδ(q), where J(p,q) bis the Jacobian of the system evaluated with the basic solution (q) band Q(p,q) baccounts for the temporal dependence of the problem. This generalized eigenvalue problem is solved using the MATLAB eigs function. D. Control parameters To nondimensionalize the system we use the tube radius R, the surface tension γ, and the liquid viscosity μl. The resulting problem is governed by five dimensionless parameters, Re =ρlUl2R μl ,Cal=μlUl γ,ˆρ=ρg ρl ,ˆμ=μg μl ,ˆ V=V (4/3)πR3.(5) Here, Re is a Reynolds number the based on the mean liquid speed Ul,Ca lis the capillary number, ˆ Vthe dimensionless bubble volume, and ˆρand ˆμare the density and viscosity ratios, respectively. In this work, we fix density and viscosity ratio, ˆρ=0.001 and ˆμ=0.01. In the simulations, the length of the computational domain will be kept fixed, ˆ H=H R=15. We will also consider the fixed bubble volume, ˆ V=3.3. This value is large enough to ensure that the bubble length is large versus the tube radius. III. NONLINEAR UNSTEADY OPENFOAM SIMULATIONS To verify both the steady stable solutions and the results of the linear stability analysis, nonlinear simulations have been performed with a different numerical method. Simulations were run using the TwoPhaseFlow library [21] in ESI OpenFOAM v2106. The library solves the unsteady Navier-Stokes equations for two immiscible, Newtonian fluids in incompressible flow, using a geometric volume of fluid (VOF) method [22,23] to capture the liquid-gas interface dynamics. The computational model is based on the solution of continuity and momentum equations, formulated as ∇·v=0,(6) ∂(ρv) ∂t+∇·(ρvv)=−∇p+∇·[μ∇v+μ(∇v)T]+Fγ,(7) where vindicates the fluid velocity, ρthe mixture fluid density, tthe time, pthe pressure, μthe mixture fluid dynamic viscosity, and Fγthe surface tension force vector. The surface tension force is calculated using the continuum surface force method [24]asFγ=γκ∇α, where the interface curvature κis estimated using derivatives of a reconstructed distance function from the interface [21]. The scalar αdenotes the liquid volume fraction, which is updated as time elapses by the solution of the transport equation ∂α ∂t+v·∇α=0.(8) The mixture fluid properties are evaluated as weighted averages of the liquid-only and gas-only properties, for example, the density ρ=αρl+(1 −α)ρg. The governing equations are integrated in time with a first-order implicit scheme, except for the volume fraction equation which is integrated explicitly using the isoAdvector library [22]. All spatial derivatives are discretized using second-order methods. The pressure-velocity coupling is handled using OpenFOAM’s built-in PISO (pressure implicit splitting of operators) algorithm [25]. The time-step of the simulations is adaptive and limited by a maximum allowed Courant number of 0.1. The flow equations are discretized and solved on a two-dimensional axisymmetric domain, which reproduces a circular channel of radius Rand length L=80R. At the inlet boundary, liquid enters the 053603-6 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … FIG. 3. Bubble profile, contours of velocity field, and streamlines obtained using the TwoPhaseFlow [21] library in OpenFOAM, for Cal=0.005 and Re =1800. The velocity contours are rescaled by the mean liquid speed in the channel Ul. The streamlines are calculated on a frame of reference anchored to the bubble, according to the velocity field v−Ubez. The axial coordinate is ˆz=z/R. flow domain with a fully-developed parabolic velocity profile, of mean velocity Ul. At the channel wall, a no-slip condition is imposed. At the outlet boundary, the pressure is set to a zero reference value while the velocity gradient along the stream direction is set to zero. As an initial condition for the volume fraction field, an elongated bubble is patched near the inlet boundary of the domain. The bubble volume and fluid properties are set to match the values chosen in the linear stability model and the corresponding capillary and Reynolds numbers. Simulations are run over time until the flow reaches a steady state (where this exists) or a steady-periodic regime. To check this, the bubble’s axial length, Lz(t), defined as the axial distance between the tail end and the head end of the bubble, is monitored as a function of time. The computational mesh is a structured orthogonal grid made of square cells, which are gradually refined in a near-wall layer; the thinnest near-wall computational cell has a thickness of 0.0007R and this fine near-wall mesh ensures that the dynamics of liquid films of thickness down to values of 0.005Rare always well resolved. This numerical setup was recently validated by El Mellas et al. [26] upon a comparison of the liquid film thickness measured around the long bubbles with the empirical correlation of Aussillous and Quéré [4] within the range of Cal=0.002 to 0.5, with deviations always within 5%. Additional validation tests against experimental data and empirical correlations for flows in the visco-inertial regime were presented by Khodaparast et al. [7] and Magnini et al. [1]. The bubble profile, contours of velocity field, and streamlines obtained using the TwoPhaseFlow [21] library in OpenFOAM for a representative case run with Cal=0.005 and Re =1800 are shown in Fig. 3. This case is selected because, having the smallest capillary number within the range of interest for this study, it is more prone to be affected by spurious velocities [27]. The TwoPhaseFlow library is very effective in minimizing spurious velocities, which are negligible in Fig. 3. The streamlines of the bubble-relative velocity field in the liquid near the front and rear menisci of the bubble exhibit the typical toroidal vortices, with the liquid impinging the bubble rear at the channel axis and recirculating from the wall toward the channel center ahead of the bubble. IV. RESULTS The present section is devoted to an analysis of the effect of increasing the inertia, characterized by the Reynolds number Re, on the flow structure and stability of the flow for various values of the capillary number Cal. We predict the bubble capillary number, Cab=μlUb/γ , where Ubis bubble speed, the steady dimensionless axial bubble length ˆ Lzb =Lzb/R, and the global linear stability 053603-7 MIRCO MAGNINI AND MIGUEL A. HERRADA -6 -4 -2 0 0 0.1 0.2 0.3 0.4 Re=0 Re=700 Re=1800 -6 -4 -2 0 0 0.1 0.2 0.3 0.4 Re=0 Re=400 Re=900 -6 -4 -2 0 0 0.1 0.2 0.3 0.4 Re=0 Re=200 Re=500 -6 -4 -2 0 0 0.1 0.2 0.3 0.4 Re=0 Re=150 Re=350 (a) (c) (b) (d) (III) (I) (II) FIG. 4. Steady bubble profiles ˆ hb, obtained with the numerical method described in Sec. II, as a function of the axial coordinate ˆz=z/R, for different Reynolds numbers keeping capillary number fixed: (a) Cal=0.05, (b) Cal=0.01, (c) Cal=0.02, and (d) Cal=0.04. of the axisymmetric base as a function of Re for elongated bubbles. Unless otherwise specified, the results described below were obtained with the numerical method introduced in Sec. II, while the results of the nonlinear simulations run with OpenFOAM serve as a validation of the model’s predictions. A. Axisymmetric steady solutions The axisymmetric, steady profiles of the bubble for four values of the liquid capillary number, Cal=0.005,0.01,0.02, and 0.04, and different values of the Reynolds number, are presented in Fig. 4. For each value of the capillary number, each subfigure reports the bubble profile in the absence of inertial forces (Re =0), the profile at a Reynolds number near the critical value for transition to an unsteady regime, and the profile for an intermediate value of Re. As previous studies have shown [19,20], when there is no inertia (Re =0), the dimensionless uniform film width, ˆ b= b/R, of region II increases as Calincreases as shown in Fig. 4. In this figure, the bubble profile, hb, is plotted as a function of the axial coordinate zfor different Reynolds numbers with fixed capillary numbers. As the Reynolds number increases, undulations appear on the surface and become more pronounced with increasing inertia [1]. Another noteworthy observation from Fig. 4is that the increase in film thickness with increasing inertia is accompanied by an elongation of the bubble as seen in Fig. 5(a), which shows ˆ Lzb as a function of Reynolds number for different capillary numbers. Figure 5(b) illustrates the variation in bubble capillary number, Cab, with both the Reynolds number and the capillary number. Cal represents the dimensionless average velocity imposed on the tube (flow rate) while Cabis the dimensionless velocity of the bubble. The results in this figure show that except for the case of large Cal, an increase in inertia does not significantly modify the bubble velocity. In order to gain a comprehensive understanding of the influence of inertia on the flow dynamics around the bubble, Fig. 6depicts the contours of the axial velocity and the bubble shape (pink line) for a fixed capillary number and three distinct Reynolds numbers. It can be observed that the bubble undergoes a lengthening and narrowing process, accompanied by the emergence of ripples at its posterior. The emergence of ripples and the narrowing of the bubble tail suggest that inertial forces 053603-8 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … 0 500 1000 1500 Re 5 5.5 6 6.5 7 7.5 8 0 500 1000 1500 Re 0 0.01 0.02 0.03 0.04 0.05 0.06 (a) (b) FIG. 5. (a) Steady dimensionless bubble length and (b) steady bubble capillary number as a function of Re for various values of Cal. reach a magnitude comparable to surface tension forces, enabling the modification of the bubble shape. Further insight into the flow and pressure field near the bubble rear is provided in Fig. 7, where two cases without (Cal=0.04, Re =0) and with (Cal=0.04, Re =350) undulations are compared. The streamlines in Figs. 7(a) and 7(c) are calculated from the liquid velocity field relative to the bubble, vl−Ubez. First, it can be observed that the flow field behind the bubble exhibits a recirculation region at the channel center, where the liquid at the axis impinges the bubble rear with relative speed 2Ul−Ub, deviates radially outward along the rear meniscus of the bubble, and then leaves the bubble at about ˆy=0.7, traveling upstream. Here, it joins with the liquid that bypasses FIG. 6. Axial velocity contours for Cal=0.04 and three different Reynolds numbers. Here, ˆy=ˆrcos(θ) being θthe azimuthal angle. 053603-9 MIRCO MAGNINI AND MIGUEL A. HERRADA 0 0.01 0.02 0.03 0.04 0.05 Ca l 0 500 1000 1500 2000 2500 Re steady unsteady Wep=3.75 Wep=3.5 Wep=3 65 Wep=3.8 Wep*=3.65 FIG. 15. Diagram summarizing steady and unsteady flow regimes depending on the Caland Re numbers. The symbols identify the (Cal,Re) combinations studied, with empty (full) symbols used to indicate stable (unstable) conditions. The critical Weber numbers identified for each Calare reported. The red dashed line illustrates the We∗ p=3.65 threshold, plotted by solving the equation We(2 −Ub/Ul)2/2=3.65, with Ub/Ul estimated as a function of Caland Re using Han and Shikazono [5] correlation, for Cal=0.004 −0.05. viscous forces to the square root of the product of capillary and inertial forces, the small values of Oh suggest that viscous forces are relatively weak in the cases under consideration, which allows the instability to be characterized using a single critical parameter that remains constant, namely We∗ p. V. CONCLUSIONS We have investigated the linear stability of elongated axisymmetric gas bubbles transported by a liquid flow in a horizontal cylindrical capillary tube. Our approach combined global linear stability analysis with nonlinear unsteady simulations performed with a volume of fluid method in OpenFOAM. The investigation covered a range of liquid capillary numbers Cal=0.005 −0.04 and revealed that the bubble rear becomes unstable when the Reynolds number grows above a critical value, which decreases from about Re∗=2000 for Cal=0.005, to Re∗=400 when Cal=0.04; the corresponding critical Weber number We∗=CalRe∗is also not uniform and varies in the range of We∗=10 −16 for the capillary numbers tested. The instability is ascribed to the increased inertia of the recirculating flow that trails the bubble and impinges its rear meniscus. As the Reynolds number increases for a fixed value of the capillary number, this impinging flow generates a low-pressure region along the rear meniscus of the bubble, giving rise to a nonmonotonic interfacial pressure profile. Above a certain threshold in the Reynolds number, surface tension can no longer oppose the dynamic pressure of the impinging flow and the bubble rear becomes unstable through a Hopf-type bifurcation, exhibiting periodic oscillations with well-defined period and amplitude. We have found that a modified version of the Weber number, Wep, calculated by using the relative speed of the flow impinging the bubble rear as the relevant velocity scale describing the inertial effects, is able to predict that the instability appears for a critical value of We∗ p, which remains quite uniform, We∗ p≈3.65, across the range of capillary numbers tested. The curve identifying the stability threshold in a Cal−Re diagram divides it into two regions, one where the bubble dynamics are steady (Wep<3.65) and one where the bubble rear becomes unsteady (Wep>3.65). 053603-16 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … ACKNOWLEDGMENTS The authors wish to thank Prof. H. A. Stone (Princeton University) and Prof. J. Eggers (University of Bristol) for fruitful discussions during the development of this work, and Dr. S. Khodaparast (University of Leeds) for kindly providing the experimental images shown in Fig. 1. M.M. acknowledges support from the U.K. Engineering & Physical Sciences Research Council (EPSRC), under the BONSAI (EP/T033398/1) grant. M.A.H. acknowledges funding from the Spanish Ministry of Economy, Industry, and Competitiveness under Grant No. PID2022-14095OB-C21. The OpenFOAM simulations were performed using the Sulis Tier-2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick; Sulis is funded by EPSRC Grant No. EP/T022108/1 and the HPC Midlands+consortium. We are grateful for access to the University of Nottingham’s Ada HPC service. APPENDIX: INVESTIGATION OF THREE-DIMENSIONAL EFFECTS USING OPENFOAM The experimental images depicted in Figs. 1(a)–1(c) show a nonaxisymmetric bubble rear as time-dependent patterns develop. However, Ferrari, Magnini, and Thome [9] performed 3D simulations at the same conditions and observed that the bubble rear oscillated maintaining axisymmetric profiles, see Figs. 1(d)–1(g). It is thus likely that some minor experimental errors such as imperfections in the channel shape, nonaxisymmetric release of bubbles from the T-junction, and especially the difficulty of generating sufficiently spaced bubbles at high flow rates, were the cause of the nonaxisymmetric bubble profiles in the experiment performed with high Reynolds numbers. To verify the validity of the results obtained in this work, 3D simulations were performed with OpenFOAM for Cal=0.04 and Re =350, and Cal=0.04, and Re =450, which corresponded to steady and unsteady bubble dynamics in 2D axisymmetric simulations. The mesh utilized is a block-structured mesh with near-wall refined cells to capture the lubricating film, which was already adopted and validated in the work of Yu et al. [20]. We start first by showing a comparison of 2D and 3D bubble profiles for a steady case in Fig. 16(a). Here, the 2D profile is reported only on the ˆy>0 side of the graph. It can be seen that 2D and 3D profiles match very well, especially along the rear meniscus of the bubble, where the instability triggering time-dependent patterns takes place. It is thus expected that the critical Reynolds number for the onset of an unsteady flow is the same. This is verified by the plot in Fig. 16(b), which shows the axial velocity of the bubble rear over time in the steady (Re =350) and unsteady (Re =450) cases. Hence, the 3D simulations confirm that the flow for Cal=0.04 becomes unsteady when Re ≈400, in agreement with the 2D axisymmetric results. Last, the time-dependent profiles of the bubble rear for Cal=0.04 and Re =450 are presented in Fig. 16(c), emphasizing that the bubble interface remains axisymmetric even after the onset of the oscillations. Although it is possible that the system becomes nonaxisymmetric at larger Reynolds numbers (as observed in the experiments), at conditions little above critical the bubble profiles remain axisymmetric. Although not reported here, further tests were conducted to assess the effect of buoyancy acting in the cross-stream direction, in the limit of the Bond number pertinent to the experiment of Khodaparast et al. [7]. For a channel radius of R=0.25 mm, the Bond number in the experiment was Bo =ρlgR2/σ =0.009, and thus very small. We repeated the 3D simulations for Cal=0.04 with Re =350 and 450, and observed that buoyancy had a negligible impact on the bubble profiles, causing only a 1.5% deviation on the top and bottom liquid film thicknesses compared to the Bo =0 case. Hence, we can safely conclude that for channel diameters on the order of the millimeter and below, the critical Reynolds numbers identified by the present 2D axisymmetric model remain valid. This is also consistent with the work of Moran et al. [34], who reported that buoyancy effects on the dynamics of elongated bubbles traveling in horizontal circular channels become significant only when the Bond number approaches unity. 053603-17 MIRCO MAGNINI AND MIGUEL A. HERRADA -7 -6 -5 -4 -3 -2 -1 0 -1 -0.5 0 0.5 1 OpenFOAM, 2D axisymmetric OpenFOAM, 3D -7.2 -7 -6.8 -6.6 -6.4 -1 -0.5 0 0.5 1 t/T=0 t/T=0.2 t/T=0.4 t/T=0.6 t/T=0.8 t/T=1 0 10203040 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 (a) (b) (c) FIG. 16. Results of 3D simulations run with OpenFOAM. (a) Comparison of bubble profiles for 2D axisymmetric and 3D simulations with Cal=0.04 and Re =350; the profile on a x=0 slice is reported for the 3D case. (b) Axial velocity of the rear meniscus of the bubble over time for 3D simulations with Cal=0.04 and Re =350 (steady), and Re =450 (unsteady) cases. (c) Time-dependent profiles of the rear meniscus of the bubble over a period of oscillation for a 3D simulation with Cal=0.04 and Re =450. [1] M. Magnini, A. Ferrari, J. R. Thome, and H. A. Stone, Undulations on the surface of elongated bubbles in confined gas-liquid flows, Phys. Rev. Fluids 2, 084001 (2017). [2] F. P. Bretherton, The motion of long bubbles in tubes, J. Fluid Mech. 10, 166 (1961). [3] G. I. Taylor, Deposition of a viscous fluid on the wall of a tube, J. Fluid Mech. 10, 161 (1961). [4] P. Aussillous and D. Quéré, Quick deposition of a fluid on the wall of a tube, Phys. Fluids 12, 2367 (2000). [5] Y. Han and N. Shikazono, Measurement of the liquid film thickness in microtube slug flow, Int. J. Heat Fluid Flow 30, 842 (2009). [6] M. Giavedoni and F. Saita, The rear meniscus of a long bubble steadily displacing a Newtonian liquid in a capillary tube, Phys. Fluids 11, 786 (1999). [7] S. Khodaparast, M. Magnini, N. Borhani, and J. R. Thome, Dynamics of isolated confined air bubbles in liquid flows through circular microchannels: an experimental and numerical study, Microfluid Nanofluid 19, 209 (2015). 053603-18 TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF … [8] M. Magnini, A. Beisel, A. Ferrari, and J. R. Thome, Pore-scale analysis of the minimum liquid film thickness around elongated bubbles in confined gas-liquid flows, Adv. Water Resour. 109, 84 (2017). [9] A. Ferrari, M. Magnini, and J. R. Thome, A flexible coupled level set and volume of fluid (flexCLV) method to simulate microscale two-phase flow in non-uniform and unstructured meshes, Int. J. Multiphase Flow 91, 276 (2017). [10] F. Municchi, I. El Mellas, O. K. Matar, and M. Magnini, Conjugate heat transfer effects on flow boiling in microchannels, Int. J. Heat Mass Transf. 195, 123166 (2022). [11] F. Municchi, C. N. Markides, O. K. Matar, and M. Magnini, Computational study of bubble, thin-film dynamics and heat transfer during flow boiling in non-circular microchannels, Appl. Therm. Eng. 238, 122039 (2024). [12] P. Angeli and A. Gavriilidis, Hydrodynamics of taylor flow in small channels: A review, Proc. Inst. Mech. Eng. C 222, 737 (2008). [13] S. Khodaparast, M. K. Kim, J. Silpe, and H. A. Stone, Bubble-driven detachment of bacteria from confined micro-geometries, Environ. Sci. Technol. 51, 1340 (2017). [14] D. Picchi and P. Poesio, Dispersion of a passive scalar around a taylor bubble, J. Fluid Mech. 951,A22 (2022). [15] Q. Zhang, S. M. Hassanizadeh, B. Liu, J. F. Schijven, and N. K. Karadimitriou, Effect of hydrophobicity on colloid transport during two-phase flow in a micromodel, Water Resour. Res. 50, 7677 (2014). [16] N. Borhani, B. Agostini, and J. R. Thome, A novel time strip flow visualization technique for investigation of intermittent dewetting and dryout in elongated bubble flow in a microchannel evaporator, Int. J. Heat Mass Transf. 53, 4809 (2010). [17] M. A. Herrada and J. M. Montanero, A numerical method to study the dynamics of capillary fluid systems, J. Comput. Phys. 306, 137 (2016). [18] M. A. Herrada, Jacobian Analytical Method (JAM), https://miguelherrada.github.io/JAM. [19] M. A. Herrada, Y. E. Yu, and H. A. Stone, Global stability analysis of bubbles rising in a vertical capillary with an external flow, J. Fluid Mech. 958, A45 (2023). [20] Y. E. Yu, M. Magnini, L. Zhu, S. Shim, and H. A. Stone, Non-unique bubble dynamics in a vertical capillary with an external flow, J. Fluid Mech. 911, A34 (2021). [21] H. Scheufler and J. Roenby, TwoPhaseFlow: A framework for developing two phase flow solvers in openfoam, OpenFOAM J. 3, 200 (2023). [22] J. Roenby, H. Bredmose, and H. Jasak, A computational method for sharp interface advection, R. Soc. Open Sci. 3, 160405 (2016). [23] H. Scheufler and J. Roenby, Accurate and efficient surface reconstruction from volume fraction data on general meshes, J. Comput. Phys. 383, 1 (2019). [24] J. U. Brackbill, D. B. Kothe, and C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100, 335 (1992). [25] R. I. Issa, Solution of the implicitly discretized fluid flow equations by operator-splitting, J. Comput. Phys. 62, 40 (1986). [26] I. El Mellas, F. Municchi, M. Icardi, and M. Magnini, Dynamics of long bubbles propagating through cylindrical micro-pin fin arrays, Int. J. Multiphase Flow 163, 104443 (2023). [27] T. Abadie, J. Aubin, and D. Legendre, On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks, J. Comput. Phys. 297, 611 (2015). [28] M. T. Kreutzer, F. Kapteijn, J. A. Moulijn, C. R. Kleijn, and J. J. Heiszwolf, Inertial and interfacial effects on pressure drop of Taylor flow in capillaries, AIChE J. 51, 2428 (2005). [29] D. W. Moore, The velocity of rise of distorted gas bubbles in a liquid of small viscosity, J. Fluid Mech. 23, 749 (1965). [30] P. C. Duineveld, The rise velocity and shape of bubbles in pure water at high Reynolds number, J. Fluid Mech. 292, 325 (1995). [31] M. A. Herrada and J. G. Eggers, Path instability of an air bubble rising in water, Proc. Natl. Acad. Sci. USA 120, e2216830120 (2023). 053603-19 MIRCO MAGNINI AND MIGUEL A. HERRADA [32] P. Bonnefis, J. Sierra-Ausin, D. Fabre, and J. Magnaudet, Path instability of deformable bubbles rising in Newtonian liquids: a linear study, J. Fluid Mech. 980, A19 (2024). [33] A. Cioncolini and M. Magnini, Solitary bubbles rising in quiescent liquids: A critical assessment of experimental data and high-fidelity numerical simulations, and performance evaluation of selected prediction methods, Phys. Fluids 37, 021304 (2025). [34] H. Moran, M. Magnini, C. N. Markides, and O. K. Matar, Inertial and buoyancy effects on the flow of elongated bubbles in horizontal channels, Int. J. Multiphase Flow 135, 103468 (2021). 053603-20