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

PHYSICAL REVIEW FLUIDS 10, 053603 (2025)
Two-dimensional global s abili y analysis o elonga ed bubbles
mo ing in a ho izon al ube
Mi co Magnini 1and Miguel A. He ada 2
1Depa men o Mechanical, Ma e ials and Manu ac u ing Enginee ing,
Uni e si y o No ingham, NG7 2RD, Uni ed Kingdom
2Dep . Ing. Ae ospacial y Mecánica de Fluidos, Uni e sidad de Se illa,
Camino de los Descub imien os s/n, 41092 Se illa, Spain
(Recei ed 2 Decembe 2024; accep ed 28 Ma ch 2025; published 7 May 2025)
The linea s abili y o an elonga ed axisymme ic gas bubble anspo ed by a liquid in a
capilla y ube is analyzed h ough he use o nume ical simula ions. The s udy ocuses on
he in luence o ine ia, cha ac e ized by he Reynolds numbe (Re) and he imposed low
a e, cha ac e ized by a capilla y numbe (Cal), on he s abili y o he bubble ail, which
exhibi s ipples a high Re. The nume ical app oach u ilizes a bounda y- i ed me hod in
a e e ence ame ancho ed o he bubble, combined wi h a mixed spa ial disc e iza ion
based on spec al colloca ion in he adial di ec ion and ini e di e ences in he axial
di ec ion. This amewo k enables he compu a ion o s eady nonlinea solu ions h ough
a New on i e a ion scheme and acili a es he linea s abili y analysis o hese solu ions.
Di ec nume ical simula ions using he olume o luid me hod in OpenFOAM a e also
pe o med o co obo a e he esul s o he s abili y analysis. We pe o m sys ema ic
simula ions o Cal=0.005 −0.04 and obse e ha he sys em becomes uns able, wi h he
eme gence o oscilla ions a he ea o he bubble, when he Reynolds numbe g ows abo e
a c i ical alue, designa ed as Re∗; his c i ical alue is dependen on he capilla y numbe .
The ins abili y is due o he inc eased ine ia o he eci cula ing low in he liquid behind
he bubble, which impinges i s ea meniscus. A modi ied Webe numbe Wep, based on
he ela i e eloci y be ween he ex e nal low and he bubble, is in oduced o desc ibe
he compe i ion be ween he des abilizing p essu e o ce ac ing on he bubble ea and
su ace ension. Ou esul s show ha he bubble dynamics become uns able o a c i ical
alue, We∗
p≈3.65, which emains qui e uni o m ac oss he ange o capilla y numbe s
es ed, and di ides he Cal−Re diag am in o s able and uns able egimes. Ou indings
o e insigh s in o he beha io o bubbles in mic o luidic applica ions, wi h implica ions
o hea ans e , mass ans e , and cleaning p ocesses in mic ochannels.
DOI: 10.1103/PhysRe Fluids.10.053603
I. INTRODUCTION
When a long gas bubble p opaga es wi hin a ci cula capilla y ube illed wi h a we ing liquid,
i aps a hin liquid ilm agains he channel wall which, su icien ly a om he on and ea
menisci, [1] akes a cons an hickness. The hickness o his ilm depends on he capilla y numbe ,
[2,3]Ca=μlU/γ , wi h μlbeing he liquid iscosi y, U he bubble o liquid a e age speed, and γ
he su ace ension, which quan i ies he impo ance o iscous e ec s o e su ace ension, and he
Published by he Ame ican Physical Socie y unde he e ms o he C ea i e Commons A ibu ion 4.0
In e na ional license. Fu he dis ibu ion o his wo k mus main ain a ibu ion o he au ho (s) and he
published a icle’s i le, jou nal ci a ion, and DOI.
2469-990X/2025/10(5)/053603(20) 053603-1 Published by he Ame ican Physical Socie y
MIRCO MAGNINI AND MIGUEL A. HERRADA
(a)
(d) (e) ( ) (g)
(b) (c)
FIG. 1. (a)–(c) Expe imen al [7] and (d)–(g) compu a ional [9] snapsho s o he ea o a long gas bubble
a eling in a liquid- illed capilla y o Ca =0.023 and Re =1065. Flow is om le o igh . In he expe imen ,
ai and wa e we e used as luids, he channel was ho izon al wi h a diame e o 0.5 mm, and he bubbles we e
imaged om he op o he channel wi h a ame a e o 104s−1. Th ee-dimensional simula ions we e un wi h
a coupled olume o luid and le el se sol e [9] in OpenFOAM.
Reynolds numbe , [4,5]Re=2ρlUR/μl, wi h ρlbeing he liquid densi y and R he channel adius,
which desc ibes he compe i ion o ine ia and iscous o ces. In he limi ha Ca 1 and Re 1,
he bubble a els s eadily, he shape o bo h he on and ea menisci o he bubble app oxima e
well ha o a sphe e o adius R, and an undula ion appea s on he liquid ilm a he ma ching
poin wi h he ea meniscus [2]. As bo h capilla y and Reynolds numbe s inc ease, he cu a u e
o he bubble nose inc eases and ha o he bubble ail dec eases, hence he bubble ail la ens
[6,7]. When he Webe numbe o he low, We =Ca Re =2ρlU2R/γ , g ows o abo e 0.1, ine ial
e ec s become impo an and mani es wi h he appea ance o mul iple (s a iona y) undula ions
on he bubble su ace, a he ma ching poin be ween he ea meniscus and he hin liquid ilm
[1]. These undula ions a e al eady p esen in he B e he on’s solu ion, howe e , o We <0.1 hei
ampli ude decays quickly such ha only one c es emains isible. As he Webe numbe is inc eased
abo e 0.1, he decay a e o he undula ions dec eases such ha mo e c es s become isible and hei
ampli ude inc eases [8]. Magnini e al. [1] epo ed ha upon a u he inc ease o he Webe numbe
abo e a alue o abou 10, ime-dependen pa e ns appea ed on he ea meniscus, which s a ed
o oscilla e, whe eas he bubble nose con inued a eling a a cons an speed. A simila obse a ion
was epo ed by Khodapa as e al. [7], who pe o med expe imen al measu emen s o bubble shape,
speed, ilm hickness, and on and ea menisci cu a u es, o long ai bubbles anspo ed in wa e
in a 0.5 mm ho izon al ci cula channel. When he capilla y numbe was below 0.01 he bubbles
a eled s eadily, bu o highe capilla y numbe s (and Reynolds numbe s app oaching 103) he
ea meniscus o he bubble was obse ed o lap signi ican ly as he bubble p og essed along he
channel, see snapsho s in Fig. 1, wi h a consequen impac on he su ace undula ions.
Besides he quali a i e obse a ions ou lined abo e, no sys ema ic s udy has been ye a emp ed
o cha ac e ize hese uns eady pa e ns o explain hei appea ance. The low o elonga ed bubbles
in mic ochannels a high Webe numbe is now ecei ing inc easing a en ion due o eme ging
applica ions u ilizing low- iscosi y luids such as wo-phase cooling in mic o-e apo a o s [10,11],
mic o-chemical eac o s [12], cleaning o bac e ial cells om medical su aces [13], geological
CO2seques a ion [14], o en i onmen al p oblems such as de achmen and anspo o pollu an s
053603-2
TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF …
FIG. 2. Ske ch o he low geome y conside ed in his s udy: an axisymme ic bubble in a ho izon al ube
wi h an applied axial low.
in unsa u a ed soil [15]. Uns eady pa e ns nea he bubble ail may ha e a de imen al e ec on he
applica ions in oduced abo e, because hey may p omo e d aining and dewe ing o he lub ica ing
ilm, comp omising he wall- luid exchanges such as hea /mass ans e and cleaning e iciency.
Fo example, Bo hani, Agos ini, and Thome [16] u ilized high-speed imaging o he wo-phase low
wi hin a mic ochannel e apo a o and epo ed ha liquid ilm d you was ini ia ed om d y pa ches
o ming in he p oximi y o he bubble ail, due o a combined e ec o hyd odynamic ins abili ies
and liquid e apo a ion.
In his wo k, he linea s abili y o elonga ed bubbles in ho izon al ubes will be analyzed using
a a ian o he nume ical echnique de eloped by He ada and Mon ane o [17] o explain he
p e iously men ioned obse a ions and iden i y he ansi ion be ween s eady and uns eady low
egimes. The esul s o he linea s abili y analysis will also be compa ed o hose ob ained wi h
nons a iona y nume ical simula ions using OpenFOAM.
The emainde o he pape is s uc u ed as ollows. In Sec. II, we p esen he go e ning equa-
ions and he nume ical echniques employed o pe o m he compu a ions o s eady axisymme ic
lows and s abili y analysis. In Sec. III, we p esen he de ails o he nonlinea uns eady simula ions
pe o med wi h OpenFOAM. In Sec. IV, we discuss he base s a e esul s and p esen he linea
s abili y analyses, in addi ion o he compa isons wi h he esul s ob ained wi h OpenFOAM. Finally,
in Sec. V, we p o ide concluding ema ks.
II. MODEL FORMULATION
We conside he con igu a ion ske ched in Fig. 2, whe e an elonga ed bubble o a gas o densi y
ρg, iscosi y μg, and olume Vis mo ing in a cylind ical ho izon al capilla y o adius Rcon aining
a liquid o densi y ρland iscosi y μl. We assume ha he low is axisymme ic and hus we
adop a wo-dimensional axisymme ic model, he alidi y o which was e i ied by means o
h ee-dimensional simula ions un wi h OpenFOAM, as documen ed in he Appendix. To model
he mo emen o he bubble unde he ac ion o an imposed low a e h ough he ube Q, he sys em
is sol ed in a ame mo ing wi h he bubble whe e a nonine ial axisymme ic cylind ical coo dina e
sys em ( ,z) aligned wi h he ube axis and ancho ed o he bubble is selec ed. This sys em is, in u n,
linked o an ine ial cylind ical coo dina e sys em ( o,zo) ancho ed o he ube walls and also aligned
wi h he ube axis. To model he appa en o ces in he momen um equa ions in his nonine ial
coo dina e sys em, we ha e o include all o he accele a ions o he sys em. In ou case, since bo h
053603-3
MIRCO MAGNINI AND MIGUEL A. HERRADA
coo dina e sys ems a e aligned, we will only ha e o conside he axial accele a ion a=dU
d eZo he
e e ence ame, whe e U( ) is he axial eloci y o he igh poin o he liquid-gas in e ace (g een
a ow in Fig. 2) and Ul=Q/(πR2) is he mean eloci y o he liquid lowing h ough he ube (bo h
eloci ies measu ed in he ine ial coo dina e sys em). Figu e 2also indica es he compu a ional
domain used in he simula ions (a cylinde o adius Rand leng h H) and he bounda y condi ions
applied a he ou e bounda ies o his domain.
A. Go e ning equa ions
The conse a ion o mass and a balance o linea momen um in he liquid (i=l) and gas (i=g)
subdomains is gi en by
∇· i=0,(i=l,g),(1a)
ρi∂ i
∂ +( i·∇) i=∇·σi,(i=l,g),(1b)
whe e i=wiez+uie is he eloci y ield and σiis he s ess enso o ma e ial i(i=g,l). We
conside in bo h egions he incomp essible low o a New onian luid, whe e he s ess enso akes
he o m
σi=−piI+μi[∇ i+(∇ i)T],(i=l,g) (1c)
whe e pi(i=l,g) is he educed p essu e, pi=p∗
i+ρidU
d z,p∗is he s a ic p essu e, and he o he
e m a ises om he po en ial associa ed wi h he linea accele a ion o he ame o e e ence.
A he le bounda y, z=−2H/3, (see Fig. 2) we assume a pa abolic axial eloci y p o ile,
wl=−U( )+2Ul[1 −( /R)2],ul=0.(1d)
A he igh bounda y, z=H/3, he s a ic p essu e is se o ze o,
p∗
l=0.(1e)
A he wall, =R, no-slip bounda y condi ions a e applied,
wl=−U( ),ul=0.(1 )
Ac oss he gas-liquid in e ace pa ame ized in e ms o a me idional a c leng h s(0 ⩽s⩽1)
I=F(s, ) and zI=G(s, ), we impose ha he eloci y ield mus be con inuous, i.e.,
wl=wg,ul=ug,(1g)
and impose a balance o no mal and angen ial s esses be ween he liquid and he gas in he o m
n·(σl−σg)·n=γκ +(ρl−ρg)dU
d zI,(1h)
·(σl−σg)·n=0,(1i)
whe e
n=Gse −Fsez
G2
s+F2
s1/2, =Gsez+Fse
F2
s+G2
s1/2,(1j)
a e no mal (n) and angen ial ec o s ( ) o he su ace, he subsc ip s s ep esen s de i a i es wi h
espec o s,κ=∇·nis wice he mean cu a u e, and γis he su ace ension. In addi ion, he
kinema ic bounda y condi ion on he in e ace is w i en
ul−∂F
∂ ∂G
∂s−wl−∂G
∂ ∂F
∂s=0.(1k)
053603-4
TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF …
To ensu e a uni o m dis ibu ion o poin s along he a c leng h s, he ollowing equa ion is
applied:
∂F
∂s
∂2F
∂s2+∂G
∂s
∂2G
∂s2=0.(1l)
A he axis, =0, egula i y condi ions a e conside ed.
Finally, o close he p oblem we need an addi ional condi ion o he bubble eloci y U( ), which
is unknown because he liquid low a e Qis he con olled pa ame e in his wo k (and, usually, in
expe imen s) and hus he bubble speed is a esul o he simula ion. This is achie ed by imposing
ha a pa icula poin o he in e ace has a ze o axial eloci y in he nonine ial ame o e e ence.
The e o e, we choose ha he igh poin o he in e ace c ossing he axis e i ies
wl=0,zI=G=0, I=F=0,a s =1.(1m)
The go e ning equa ions a e in eg a ed wi h a a ian o he nume ical me hod p oposed by He ada
and Mon ane o [17] and He ada [18]. The spa ial domain occupied by he liquid and he gas a e
mapped in o wo ec angula domains [0 ⩽sl⩽1] ×[0 ⩽ηl⩽1] and [0 ⩽sg⩽1] ×[0 ⩽ηg⩽
1] using quasi-ellip ic ans o ma ions [19]. These mappings a e applied o he go e ning equa-
ions (1) and he esul ing equa ions a e disc e ized in he η-di ec ion wi h nηland nηgChebyshe
spec al colloca ion poin s in he liquid and gas domains, espec i ely. Con e sely, in he s-di ec ion,
we use second-o de ini e di e ences wi h nsl and nsequally spaced poin s in he liquid and gas
domains, espec i ely. The esul s p esen ed in his wo k we e ob ained using ns=851, nsl =1353,
nηg=10, and nηl=11.
B. Axisymme ic s eady solu ions
The s eady-s a e solu ions o he nonlinea equa ions (1) wi h all a iables independen o ime
a e ob ained by sol ing all disc e ized equa ions simul aneously (a so-called monoli hic scheme)
using a New on echnique.
We deno e he s eady axisymme ic solu ion o he sys em wi h he subsc ip b. The s eady bubble
p o ile hb(s), is de ined as
hb(s)=R−Fb(s),0⩽s⩽1.(2)
As desc ibed in [20], hb o a long bubble wi hou ine ia (Re =0) is cha ac e ized by h ee
dis inc egions: (I) he bubble “nose” hb= c e; (II), uni o m ilm egion hb(s)=c e =b; and (III),
he bubble “ ail” hb= c e (see Fig. 4). We ace he s eady solu ions as a unc ion o he model
pa ame e s and quan i y hem using he s eady axial bubble leng h, Lzb, de ined as
Lzb =Gb(s=1) −Gb(s=0).(3)
C. Small ampli ude 2D pe u ba ions
To es he s abili y o a gi en s eady axisymme ic s a e we calcula e he linea global modes by
assuming he empo al dependence
( ,z; )=b( ,z)+δ( ,z)e−iω (1),(4)
whe e ( ,z; ) ep esen s any dependen a iable while b( ,z) and δ( ,z) deno e he base
(s eady) solu ion and he spa ial dependence o he eigenmode o ha a iable, espec i ely, while
ω=ω +iωiis he equency (an eigen alue). Bo h he eigen equencies and he co esponding
eigenmodes a e calcula ed as a unc ion o he go e ning pa ame e s. The dominan eigenmode is
ha wi h he la ges g ow h ac o ωi, so i i is posi i e, he base low is asymp o ically uns able.
The nume ical p ocedu e used o sol e he s eady p oblem can be easily adap ed o sol e he
eigen alue p oblem ha de e mines he linea global modes o he sys em. In his case, he empo al
053603-5

MIRCO MAGNINI AND MIGUEL A. HERRADA
de i a i e is compu ed assuming he dependence (4). The spa ial dependence o he linea pe u ba-
ion δ(q)is he solu ion o he gene alized eigen alue p oblem J(p,q)
bδ(q)=iωQ(p)
bδ(q), whe e
J(p,q)
bis he Jacobian o he sys em e alua ed wi h he basic solu ion (q)
band Q(p,q)
baccoun s o
he empo al dependence o he p oblem. This gene alized eigen alue p oblem is sol ed using he
MATLAB eigs unc ion.
D. Con ol pa ame e s
To nondimensionalize he sys em we use he ube adius R, he su ace ension γ, and he liquid
iscosi y μl. The esul ing p oblem is go e ned by i e dimensionless pa ame e s,
Re =ρlUl2R
μl
,Cal=μlUl
γ,ˆρ=ρg
ρl
,ˆμ=μg
μl
,ˆ
V=V
(4/3)πR3.(5)
He e, Re is a Reynolds numbe he based on he mean liquid speed Ul,Ca
lis he capilla y numbe ,
ˆ
V he dimensionless bubble olume, and ˆρand ˆμa e he densi y and iscosi y a ios, espec i ely.
In his wo k, we ix densi y and iscosi y a io, ˆρ=0.001 and ˆμ=0.01. In he simula ions, he
leng h o he compu a ional domain will be kep ixed, ˆ
H=H
R=15. We will also conside he
ixed bubble olume, ˆ
V=3.3. This alue is la ge enough o ensu e ha he bubble leng h is la ge
e sus he ube adius.
III. NONLINEAR UNSTEADY OPENFOAM SIMULATIONS
To e i y bo h he s eady s able solu ions and he esul s o he linea s abili y analysis, non-
linea simula ions ha e been pe o med wi h a di e en nume ical me hod. Simula ions we e un
using he TwoPhaseFlow lib a y [21] in ESI OpenFOAM 2106. The lib a y sol es he uns eady
Na ie -S okes equa ions o wo immiscible, New onian luids in incomp essible low, using a
geome ic olume o luid (VOF) me hod [22,23] o cap u e he liquid-gas in e ace dynamics. The
compu a ional model is based on he solu ion o con inui y and momen um equa ions, o mula ed
as
∇· =0,(6)
∂(ρ )
∂ +∇·(ρ )=−∇p+∇·[μ∇ +μ(∇ )T]+Fγ,(7)
whe e indica es he luid eloci y, ρ he mix u e luid densi y, he ime, p he p essu e, μ he
mix u e luid dynamic iscosi y, and Fγ he su ace ension o ce ec o . The su ace ension o ce
is calcula ed using he con inuum su ace o ce me hod [24]asFγ=γκ∇α, whe e he in e ace
cu a u e κis es ima ed using de i a i es o a econs uc ed dis ance unc ion om he in e ace
[21]. The scala αdeno es he liquid olume ac ion, which is upda ed as ime elapses by he
solu ion o he anspo equa ion
∂α
∂ + ·∇α=0.(8)
The mix u e luid p ope ies a e e alua ed as weigh ed a e ages o he liquid-only and gas-only
p ope ies, o example, he densi y ρ=αρl+(1 −α)ρg.
The go e ning equa ions a e in eg a ed in ime wi h a i s -o de implici scheme, excep o
he olume ac ion equa ion which is in eg a ed explici ly using he isoAd ec o lib a y [22].
All spa ial de i a i es a e disc e ized using second-o de me hods. The p essu e- eloci y coupling
is handled using OpenFOAM’s buil -in PISO (p essu e implici spli ing o ope a o s) algo i hm
[25]. The ime-s ep o he simula ions is adap i e and limi ed by a maximum allowed Cou an
numbe o 0.1.
The low equa ions a e disc e ized and sol ed on a wo-dimensional axisymme ic domain, which
ep oduces a ci cula channel o adius Rand leng h L=80R. A he inle bounda y, liquid en e s he
053603-6
TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF …
FIG. 3. Bubble p o ile, con ou s o eloci y ield, and s eamlines ob ained using he TwoPhaseFlow [21]
lib a y in OpenFOAM, o Cal=0.005 and Re =1800. The eloci y con ou s a e escaled by he mean liquid
speed in he channel Ul. The s eamlines a e calcula ed on a ame o e e ence ancho ed o he bubble,
acco ding o he eloci y ield −Ubez. The axial coo dina e is ˆz=z/R.
low domain wi h a ully-de eloped pa abolic eloci y p o ile, o mean eloci y Ul. A he channel
wall, a no-slip condi ion is imposed. A he ou le bounda y, he p essu e is se o a ze o e e ence
alue while he eloci y g adien along he s eam di ec ion is se o ze o. As an ini ial condi ion o
he olume ac ion ield, an elonga ed bubble is pa ched nea he inle bounda y o he domain. The
bubble olume and luid p ope ies a e se o ma ch he alues chosen in he linea s abili y model
and he co esponding capilla y and Reynolds numbe s. Simula ions a e un o e ime un il he low
eaches a s eady s a e (whe e his exis s) o a s eady-pe iodic egime. To check his, he bubble’s
axial leng h, Lz( ), de ined as he axial dis ance be ween he ail end and he head end o he bubble,
is moni o ed as a unc ion o ime.
The compu a ional mesh is a s uc u ed o hogonal g id made o squa e cells, which a e g adually
e ined in a nea -wall laye ; he hinnes nea -wall compu a ional cell has a hickness o 0.0007R
and his ine nea -wall mesh ensu es ha he dynamics o liquid ilms o hickness down o alues
o 0.005Ra e always well esol ed. This nume ical se up was ecen ly alida ed by El Mellas e al.
[26] upon a compa ison o he liquid ilm hickness measu ed a ound he long bubbles wi h he
empi ical co ela ion o Aussillous and Qué é [4] wi hin he ange o Cal=0.002 o 0.5, wi h
de ia ions always wi hin 5%. Addi ional alida ion es s agains expe imen al da a and empi ical
co ela ions o lows in he isco-ine ial egime we e p esen ed by Khodapa as e al. [7] and
Magnini e al. [1].
The bubble p o ile, con ou s o eloci y ield, and s eamlines ob ained using he TwoPhaseFlow
[21] lib a y in OpenFOAM o a ep esen a i e case un wi h Cal=0.005 and Re =1800 a e
shown in Fig. 3. This case is selec ed because, ha ing he smalles capilla y numbe wi hin he
ange o in e es o his s udy, i is mo e p one o be a ec ed by spu ious eloci ies [27]. The
TwoPhaseFlow lib a y is e y e ec i e in minimizing spu ious eloci ies, which a e negligible in
Fig. 3. The s eamlines o he bubble- ela i e eloci y ield in he liquid nea he on and ea
menisci o he bubble exhibi he ypical o oidal o ices, wi h he liquid impinging he bubble ea
a he channel axis and eci cula ing om he wall owa d he channel cen e ahead o he bubble.
IV. RESULTS
The p esen sec ion is de o ed o an analysis o he e ec o inc easing he ine ia, cha ac e ized
by he Reynolds numbe Re, on he low s uc u e and s abili y o he low o a ious alues o he
capilla y numbe Cal. We p edic he bubble capilla y numbe , Cab=μlUb/γ , whe e Ubis bubble
speed, he s eady dimensionless axial bubble leng h ˆ
Lzb =Lzb/R, and he global linea s abili y
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. S eady bubble p o iles ˆ
hb, ob ained wi h he nume ical me hod desc ibed in Sec. II, as a unc ion o
he axial coo dina e ˆz=z/R, o di e en Reynolds numbe s keeping capilla y numbe ixed: (a) Cal=0.05,
(b) Cal=0.01, (c) Cal=0.02, and (d) Cal=0.04.
o he axisymme ic base as a unc ion o Re o elonga ed bubbles. Unless o he wise speci ied,
he esul s desc ibed below we e ob ained wi h he nume ical me hod in oduced in Sec. II, while
he esul s o he nonlinea simula ions un wi h OpenFOAM se e as a alida ion o he model’s
p edic ions.
A. Axisymme ic s eady solu ions
The axisymme ic, s eady p o iles o he bubble o ou alues o he liquid capilla y numbe ,
Cal=0.005,0.01,0.02, and 0.04, and di e en alues o he Reynolds numbe , a e p esen ed in
Fig. 4. Fo each alue o he capilla y numbe , each sub igu e epo s he bubble p o ile in he
absence o ine ial o ces (Re =0), he p o ile a a Reynolds numbe nea he c i ical alue o
ansi ion o an uns eady egime, and he p o ile o an in e media e alue o Re. As p e ious s udies
ha e shown [19,20], when he e is no ine ia (Re =0), he dimensionless uni o m ilm wid h, ˆ
b=
b/R, o egion II inc eases as Calinc eases as shown in Fig. 4. In his igu e, he bubble p o ile, hb,
is plo ed as a unc ion o he axial coo dina e z o di e en Reynolds numbe s wi h ixed capilla y
numbe s. As he Reynolds numbe inc eases, undula ions appea on he su ace and become mo e
p onounced wi h inc easing ine ia [1].
Ano he no ewo hy obse a ion om Fig. 4is ha he inc ease in ilm hickness wi h inc easing
ine ia is accompanied by an elonga ion o he bubble as seen in Fig. 5(a), which shows ˆ
Lzb as a
unc ion o Reynolds numbe o di e en capilla y numbe s. Figu e 5(b) illus a es he a ia ion
in bubble capilla y numbe , Cab, wi h bo h he Reynolds numbe and he capilla y numbe . Cal
ep esen s he dimensionless a e age eloci y imposed on he ube ( low a e) while Cabis he
dimensionless eloci y o he bubble. The esul s in his igu e show ha excep o he case o la ge
Cal, an inc ease in ine ia does no signi ican ly modi y he bubble eloci y.
In o de o gain a comp ehensi e unde s anding o he in luence o ine ia on he low dynamics
a ound he bubble, Fig. 6depic s he con ou s o he axial eloci y and he bubble shape (pink line)
o a ixed capilla y numbe and h ee dis inc Reynolds numbe s. I can be obse ed ha he bubble
unde goes a leng hening and na owing p ocess, accompanied by he eme gence o ipples a i s
pos e io . The eme gence o ipples and he na owing o he bubble ail sugges ha ine ial o ces
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) S eady dimensionless bubble leng h and (b) s eady bubble capilla y numbe as a unc ion o Re
o a ious alues o Cal.
each a magni ude compa able o su ace ension o ces, enabling he modi ica ion o he bubble
shape.
Fu he insigh in o he low and p essu e ield nea he bubble ea is p o ided in Fig. 7,
whe e wo cases wi hou (Cal=0.04, Re =0) and wi h (Cal=0.04, Re =350) undula ions a e
compa ed. The s eamlines in Figs. 7(a) and 7(c) a e calcula ed om he liquid eloci y ield ela i e
o he bubble, l−Ubez. Fi s , i can be obse ed ha he low ield behind he bubble exhibi s a
eci cula ion egion a he channel cen e , whe e he liquid a he axis impinges he bubble ea wi h
ela i e speed 2Ul−Ub, de ia es adially ou wa d along he ea meniscus o he bubble, and hen
lea es he bubble a abou ˆy=0.7, a eling ups eam. He e, i joins wi h he liquid ha bypasses
FIG. 6. Axial eloci y con ou s o Cal=0.04 and h ee di e en Reynolds numbe s. He e, ˆy=ˆ cos(θ)
being θ he azimu hal 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
s eady
uns eady
Wep=3.75
Wep=3.5
Wep=3 65 Wep=3.8
Wep*=3.65
FIG. 15. Diag am summa izing s eady and uns eady low egimes depending on he Caland Re numbe s.
The symbols iden i y he (Cal,Re) combina ions s udied, wi h emp y ( ull) symbols used o indica e s able
(uns able) condi ions. The c i ical Webe numbe s iden i ied o each Cala e epo ed. The ed dashed line
illus a es he We∗
p=3.65 h eshold, plo ed by sol ing he equa ion We(2 −Ub/Ul)2/2=3.65, wi h Ub/Ul
es ima ed as a unc ion o Caland Re using Han and Shikazono [5] co ela ion, o Cal=0.004 −0.05.
iscous o ces o he squa e oo o he p oduc o capilla y and ine ial o ces, he small alues
o Oh sugges ha iscous o ces a e ela i ely weak in he cases unde conside a ion, which
allows he ins abili y o be cha ac e ized using a single c i ical pa ame e ha emains cons an ,
namely We∗
p.
V. CONCLUSIONS
We ha e in es iga ed he linea s abili y o elonga ed axisymme ic gas bubbles anspo ed
by a liquid low in a ho izon al cylind ical capilla y ube. Ou app oach combined global linea
s abili y analysis wi h nonlinea uns eady simula ions pe o med wi h a olume o luid me hod in
OpenFOAM. The in es iga ion co e ed a ange o liquid capilla y numbe s Cal=0.005 −0.04 and
e ealed ha he bubble ea becomes uns able when he Reynolds numbe g ows abo e a c i ical
alue, which dec eases om abou Re∗=2000 o Cal=0.005, o Re∗=400 when Cal=0.04;
he co esponding c i ical Webe numbe We∗=CalRe∗is also no uni o m and a ies in he ange
o We∗=10 −16 o he capilla y numbe s es ed. The ins abili y is asc ibed o he inc eased
ine ia o he eci cula ing low ha ails he bubble and impinges i s ea meniscus. As he
Reynolds numbe inc eases o a ixed alue o he capilla y numbe , his impinging low gene a es a
low-p essu e egion along he ea meniscus o he bubble, gi ing ise o a nonmono onic in e acial
p essu e p o ile. Abo e a ce ain h eshold in he Reynolds numbe , su ace ension can no longe
oppose he dynamic p essu e o he impinging low and he bubble ea becomes uns able h ough
a Hop - ype bi u ca ion, exhibi ing pe iodic oscilla ions wi h well-de ined pe iod and ampli ude.
We ha e ound ha a modi ied e sion o he Webe numbe , Wep, calcula ed by using he ela i e
speed o he low impinging he bubble ea as he ele an eloci y scale desc ibing he ine ial
e ec s, is able o p edic ha he ins abili y appea s o a c i ical alue o We∗
p, which emains
qui e uni o m, We∗
p≈3.65, ac oss he ange o capilla y numbe s es ed. The cu e iden i ying he
s abili y h eshold in a Cal−Re diag am di ides i in o wo egions, one whe e he bubble dynamics
a e s eady (Wep<3.65) and one whe e he bubble ea becomes uns eady (Wep>3.65).
053603-16

TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF …
ACKNOWLEDGMENTS
The au ho s wish o hank P o . H. A. S one (P ince on Uni e si y) and P o . J. Egge s
(Uni e si y o B is ol) o ui ul discussions du ing he de elopmen o his wo k, and D . S.
Khodapa as (Uni e si y o Leeds) o kindly p o iding he expe imen al images shown in Fig. 1.
M.M. acknowledges suppo om he U.K. Enginee ing & Physical Sciences Resea ch Council
(EPSRC), unde he BONSAI (EP/T033398/1) g an . M.A.H. acknowledges unding om he Span-
ish Minis y o Economy, Indus y, and Compe i i eness unde G an No. PID2022-14095OB-C21.
The OpenFOAM simula ions we e pe o med using he Sulis Tie -2 HPC pla o m hos ed by he
Scien i ic Compu ing Resea ch Technology Pla o m a he Uni e si y o Wa wick; Sulis is unded
by EPSRC G an No. EP/T022108/1 and he HPC Midlands+conso ium. We a e g a e ul o access
o he Uni e si y o No ingham’s Ada HPC se ice.
APPENDIX: INVESTIGATION OF THREE-DIMENSIONAL EFFECTS USING OPENFOAM
The expe imen al images depic ed in Figs. 1(a)–1(c) show a nonaxisymme ic bubble ea
as ime-dependen pa e ns de elop. Howe e , Fe a i, Magnini, and Thome [9] pe o med 3D
simula ions a he same condi ions and obse ed ha he bubble ea oscilla ed main aining ax-
isymme ic p o iles, see Figs. 1(d)–1(g). I is hus likely ha some mino expe imen al e o s such
as impe ec ions in he channel shape, nonaxisymme ic elease o bubbles om he T-junc ion, and
especially he di icul y o gene a ing su icien ly spaced bubbles a high low a es, we e he cause
o he nonaxisymme ic bubble p o iles in he expe imen pe o med wi h high Reynolds numbe s.
To e i y he alidi y o he esul s ob ained in his wo k, 3D simula ions we e pe o med wi h
OpenFOAM o Cal=0.04 and Re =350, and Cal=0.04, and Re =450, which co esponded
o s eady and uns eady bubble dynamics in 2D axisymme ic simula ions. The mesh u ilized is a
block-s uc u ed mesh wi h nea -wall e ined cells o cap u e he lub ica ing ilm, which was al eady
adop ed and alida ed in he wo k o Yu e al. [20].
We s a i s by showing a compa ison o 2D and 3D bubble p o iles o a s eady case in
Fig. 16(a). He e, he 2D p o ile is epo ed only on he ˆy>0 side o he g aph. I can be seen
ha 2D and 3D p o iles ma ch e y well, especially along he ea meniscus o he bubble, whe e
he ins abili y igge ing ime-dependen pa e ns akes place. I is hus expec ed ha he c i ical
Reynolds numbe o he onse o an uns eady low is he same. This is e i ied by he plo in
Fig. 16(b), which shows he axial eloci y o he bubble ea o e ime in he s eady (Re =350)
and uns eady (Re =450) cases. Hence, he 3D simula ions con i m ha he low o Cal=0.04
becomes uns eady when Re ≈400, in ag eemen wi h he 2D axisymme ic esul s. Las , he
ime-dependen p o iles o he bubble ea o Cal=0.04 and Re =450 a e p esen ed in Fig. 16(c),
emphasizing ha he bubble in e ace emains axisymme ic e en a e he onse o he oscilla ions.
Al hough i is possible ha he sys em becomes nonaxisymme ic a la ge Reynolds numbe s
(as obse ed in he expe imen s), a condi ions li le abo e c i ical he bubble p o iles emain
axisymme ic.
Al hough no epo ed he e, u he es s we e conduc ed o assess he e ec o buoyancy ac ing
in he c oss-s eam di ec ion, in he limi o he Bond numbe pe inen o he expe imen o
Khodapa as e al. [7]. Fo a channel adius o R=0.25 mm, he Bond numbe in he expe imen
was Bo =ρlgR2/σ =0.009, and hus e y small. We epea ed he 3D simula ions o Cal=0.04
wi h Re =350 and 450, and obse ed ha buoyancy had a negligible impac on he bubble p o iles,
causing only a 1.5% de ia ion on he op and bo om liquid ilm hicknesses compa ed o he Bo =0
case. Hence, we can sa ely conclude ha o channel diame e s on he o de o he millime e and
below, he c i ical Reynolds numbe s iden i ied by he p esen 2D axisymme ic model emain alid.
This is also consis en wi h he wo k o Mo an e al. [34], who epo ed ha buoyancy e ec s on
he dynamics o elonga ed bubbles a eling in ho izon al ci cula channels become signi ican only
when he Bond numbe app oaches uni y.
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 axisymme ic
OpenFOAM, 3D
-7.2 -7 -6.8 -6.6 -6.4
-1
-0.5
0
0.5
1
/T=0
/T=0.2
/T=0.4
/T=0.6
/T=0.8
/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. Resul s o 3D simula ions un wi h OpenFOAM. (a) Compa ison o bubble p o iles o 2D
axisymme ic and 3D simula ions wi h Cal=0.04 and Re =350; he p o ile on a x=0 slice is epo ed o
he 3D case. (b) Axial eloci y o he ea meniscus o he bubble o e ime o 3D simula ions wi h Cal=0.04
and Re =350 (s eady), and Re =450 (uns eady) cases. (c) Time-dependen p o iles o he ea meniscus o
he bubble o e a pe iod o oscilla ion o a 3D simula ion wi h Cal=0.04 and Re =450.
[1] M. Magnini, A. Fe a i, J. R. Thome, and H. A. S one, Undula ions on he su ace o elonga ed bubbles
in con ined gas-liquid lows, Phys. Re . Fluids 2, 084001 (2017).
[2] F. P. B e he on, The mo ion o long bubbles in ubes, J. Fluid Mech. 10, 166 (1961).
[3] G. I. Taylo , Deposi ion o a iscous luid on he wall o a ube, J. Fluid Mech. 10, 161 (1961).
[4] P. Aussillous and D. Qué é, Quick deposi ion o a luid on he wall o a ube, Phys. Fluids 12, 2367 (2000).
[5] Y. Han and N. Shikazono, Measu emen o he liquid ilm hickness in mic o ube slug low, In . J. Hea
Fluid Flow 30, 842 (2009).
[6] M. Gia edoni and F. Sai a, The ea meniscus o a long bubble s eadily displacing a New onian liquid in
a capilla y ube, Phys. Fluids 11, 786 (1999).
[7] S. Khodapa as , M. Magnini, N. Bo hani, and J. R. Thome, Dynamics o isola ed con ined ai bubbles in
liquid lows h ough ci cula mic ochannels: an expe imen al and nume ical s udy, Mic o luid Nano luid
19, 209 (2015).
053603-18
TWO-DIMENSIONAL GLOBAL STABILITY ANALYSIS OF …
[8] M. Magnini, A. Beisel, A. Fe a i, and J. R. Thome, Po e-scale analysis o he minimum liquid ilm
hickness a ound elonga ed bubbles in con ined gas-liquid lows, Ad . Wa e Resou . 109, 84 (2017).
[9] A. Fe a i, M. Magnini, and J. R. Thome, A lexible coupled le el se and olume o luid ( lexCLV)
me hod o simula e mic oscale wo-phase low in non-uni o m and uns uc u ed meshes, In . J. Mul iphase
Flow 91, 276 (2017).
[10] F. Municchi, I. El Mellas, O. K. Ma a , and M. Magnini, Conjuga e hea ans e e ec s on low boiling
in mic ochannels, In . J. Hea Mass T ans . 195, 123166 (2022).
[11] F. Municchi, C. N. Ma kides, O. K. Ma a , and M. Magnini, Compu a ional s udy o bubble, hin- ilm
dynamics and hea ans e du ing low boiling in non-ci cula mic ochannels, Appl. The m. Eng. 238,
122039 (2024).
[12] P. Angeli and A. Ga iilidis, Hyd odynamics o aylo low in small channels: A e iew, P oc. Ins . Mech.
Eng. C 222, 737 (2008).
[13] S. Khodapa as , M. K. Kim, J. Silpe, and H. A. S one, Bubble-d i en de achmen o bac e ia om con ined
mic o-geome ies, En i on. Sci. Technol. 51, 1340 (2017).
[14] D. Picchi and P. Poesio, Dispe sion o a passi e scala a ound a aylo bubble, J. Fluid Mech. 951,A22
(2022).
[15] Q. Zhang, S. M. Hassanizadeh, B. Liu, J. F. Schij en, and N. K. Ka adimi iou, E ec o hyd ophobici y
on colloid anspo du ing wo-phase low in a mic omodel, Wa e Resou . Res. 50, 7677 (2014).
[16] N. Bo hani, B. Agos ini, and J. R. Thome, A no el ime s ip low isualiza ion echnique o in es iga ion
o in e mi en dewe ing and d you in elonga ed bubble low in a mic ochannel e apo a o , In . J. Hea
Mass T ans . 53, 4809 (2010).
[17] M. A. He ada and J. M. Mon ane o, A nume ical me hod o s udy he dynamics o capilla y luid sys ems,
J. Compu . Phys. 306, 137 (2016).
[18] M. A. He ada, Jacobian Analy ical Me hod (JAM), h ps://miguelhe ada.gi hub.io/JAM.
[19] M. A. He ada, Y. E. Yu, and H. A. S one, Global s abili y analysis o bubbles ising in a e ical capilla y
wi h an ex e nal low, J. Fluid Mech. 958, A45 (2023).
[20] Y. E. Yu, M. Magnini, L. Zhu, S. Shim, and H. A. S one, Non-unique bubble dynamics in a e ical
capilla y wi h an ex e nal low, J. Fluid Mech. 911, A34 (2021).
[21] H. Scheu le and J. Roenby, TwoPhaseFlow: A amewo k o de eloping wo phase low sol e s in
open oam, OpenFOAM J. 3, 200 (2023).
[22] J. Roenby, H. B edmose, and H. Jasak, A compu a ional me hod o sha p in e ace ad ec ion, R. Soc.
Open Sci. 3, 160405 (2016).
[23] H. Scheu le and J. Roenby, Accu a e and e icien su ace econs uc ion om olume ac ion da a on
gene al meshes, J. Compu . Phys. 383, 1 (2019).
[24] J. U. B ackbill, D. B. Ko he, and C. Zemach, A con inuum me hod o modeling su ace ension,
J. Compu . Phys. 100, 335 (1992).
[25] R. I. Issa, Solu ion o he implici ly disc e ized luid low equa ions by ope a o -spli ing, J. Compu . Phys.
62, 40 (1986).
[26] I. El Mellas, F. Municchi, M. Ica di, and M. Magnini, Dynamics o long bubbles p opaga ing h ough
cylind ical mic o-pin in a ays, In . J. Mul iphase Flow 163, 104443 (2023).
[27] T. Abadie, J. Aubin, and D. Legend e, On he combined e ec s o su ace ension o ce calcula ion and
in e ace ad ec ion on spu ious cu en s wi hin olume o luid and le el se amewo ks, J. Compu . Phys.
297, 611 (2015).
[28] M. T. K eu ze , F. Kap eijn, J. A. Moulijn, C. R. Kleijn, and J. J. Heiszwol , Ine ial and in e acial e ec s
on p essu e d op o Taylo low in capilla ies, AIChE J. 51, 2428 (2005).
[29] D. W. Moo e, The eloci y o ise o dis o ed gas bubbles in a liquid o small iscosi y, J. Fluid Mech.
23, 749 (1965).
[30] P. C. Duine eld, The ise eloci y and shape o bubbles in pu e wa e a high Reynolds numbe , J. Fluid
Mech. 292, 325 (1995).
[31] M. A. He ada and J. G. Egge s, Pa h ins abili y o an ai bubble ising in wa e , P oc. Na l. Acad. Sci.
USA 120, e2216830120 (2023).
053603-19
MIRCO MAGNINI AND MIGUEL A. HERRADA
[32] P. Bonne is, J. Sie a-Ausin, D. Fab e, and J. Magnaude , Pa h ins abili y o de o mable bubbles ising in
New onian liquids: a linea s udy, J. Fluid Mech. 980, A19 (2024).
[33] A. Cioncolini and M. Magnini, Soli a y bubbles ising in quiescen liquids: A c i ical assessmen o expe -
imen al da a and high- ideli y nume ical simula ions, and pe o mance e alua ion o selec ed p edic ion
me hods, Phys. Fluids 37, 021304 (2025).
[34] H. Mo an, M. Magnini, C. N. Ma kides, and O. K. Ma a , Ine ial and buoyancy e ec s on he low o
elonga ed bubbles in ho izon al channels, In . J. Mul iphase Flow 135, 103468 (2021).
053603-20