On the effect of the approach velocity on the coalescence of fluid particles

ContentslistsavailableatScienceDirect

International Journal of Multiphase Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

On the effect of the approach velocity on the coalescence of fluid particles

Suat Canberk Ozan

^∗

, Hugo Atle Jakobsen

Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), Trondheim N- 7491, Norway

a rt i c l e i n f o

Article history:

Received 15 January 2019 Revised 16 May 2019 Accepted 24 July 2019 Available online 25 July 2019 Keywords:

Film drainage Coalescence time Approach velocity Boundary Integral Method

a b s t r a c t

Theeffectoftheconstantrelativeapproach velocityoftwoNewtonian fluidparticlesontheircoales- cenceisinvestigatedviaafilmdrainagemodel.Theinterfacesaredeformableandallowedtohaveany degreeofmobility.Thethinningequationisderivedbasedonthelubricationtheory,andthetangential velocityoftheinterfaceisdeterminedviatheBoundaryIntegralMethod(BIM).Thedetailsforthetreat- mentofthe inherentsingularityinBIMareprovided.Astheapproachvelocityincreases,regardlessof thestrengthofthevanderWaalsforcesandthevalueoftheviscosityratio(orthedegreeoftheinterfa- cialmobility),threetypesofbehaviorofthecoalescencetimeareidentified:thelinearslow,thedimpled andthemultiple-rimdrainageregimes.Inthefirsttworegionsthecoalescencetimedecreaseswiththe approachvelocity,eventuallypassesthroughaminimumandstartstoincreaseinthethirdregion.The slopeinthelinearregionandthecriticalvelocitiesseparatingtheregimesareusedtoquantifythebe- havior,andgivenasfunctionsoftheHamakerconstantandtheviscosityratio.Theresultsareshownto beingoodagreementwithseveralrecentexperimentalstudiesintheliterature.

© 2019TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

(http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Coalescence offluid particles plays an importantrole in wide variety of phenomenaincluding naturalones such asthe forma- tion oftheraindroplets,aswell asseveralprocessesinchemical andbiochemicalindustriessuch asfoodandbeverageproduction orpetroleum refining,where multiphaseflows arefrequently en- countered.Intheseindustries,chemicalandbiochemicalreactors, fermentationunits,boilingandcondensationequipment,andsep- arators are widely employed; theperformance and theefficiency of which strongly depend on thecharacteristics of the dispersed flowwithinthem,i.e.,thespatialdistributionofthefluidparticles in thecontinuous phase andtheir size. Therefore,it is crucialto understandthenatureoffluidparticlecoalescence(andbreakage) onthescaleoftheequipmentofinterest,inwhichalargenumber of fluid particles interact at the same time enabling coalescence andbreakageto occursimultaneouslywithin theunit.That,how- ever,firstrequirestheunderstandingofthephenomenaonamore fundamental level: the single events of coalescence of two fluid particles,andthebreakageofasingleparticle.Inthiswork,wefo- cus onlyonthecoalescence oftwofluid particles,particularlyon

∗ Corresponding author.

E-mail address: [email protected] (S.C. Ozan).

theeffectofthe relativeapproachvelocityoftheparticles onthe coalescencetime.

There has been huge interest and effort in gaining better comprehension of coalescence and its effects on multiphase en- gineering units for decades, which seemingly resulted in some conflicting observations and conclusions. One earlier example of these contradictions might be seen as in between Shinnar and Church (1960), Shinnar (1961) and Howarth (1964). The former claimsthatimmediatecoalescenceoftwofluidparticlesrarelyoc- curs,butinsteadtheparticlesrathercoherefirstentrappingathin filmofcontinuousphaseinbetween.Whereas,thelattersupports the idea of a critical approach velocity, after which immediate coalescence prevails and before which the probability of coales- cenceisverylow.AsreviewedbyLiaoandLucas(2010),thesetwo proposalslaythefoundationfortwodistinctmodelingapproaches, thefilmdrainagemodelsandtheenergymodels,respectively.On theother hand,another criticalapproach velocity is observedby Lehr etal. (2002), afterwhich the collision of thefluid particles resultin bouncing instead of coalescence and before which coa- lescence is rapid.LiaoandLucas (2010) adds Lehr etal. (2002)’s model next to the film drainage and energy models as a third branch.Theoreticalestimatesthat areon thesameorderofmag- nitudeasLehretal.(2002)’smeasurementforthecriticalvelocity, aregiveninthetheoreticalworkofChestersandHofman(1982)by https://doi.org/10.1016/j.ijmultiphaseflow.2019.07.016

0301-9322/© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license. ( http://creativecommons.org/licenses/by-nc-nd/4.0/ )

Nomenclature

Vapp Relativeapproachvelocityofthefluidparticles tc Coalescencetime

λ

c Coalescenceefficiency t_drainage Drainagetime tcontact Contacttime

λ

^{Dispersedtocontinuousphaseviscosityratio}

a Reducedparticleradius a₁,a₂ Fluidparticleradii

h₀ Minimuminitialfilmthickness

μ

d Dispersedphaseviscosity

μ

^{c Continuousphaseviscosity}

R1,R2 Fluidparticleradii

V₁,V₂ Approachvelocitiesoftheparticles r Radialdistanceinthefilm z Axialdistanceinthefilm

t Time

h Thicknessofthethinfilm

¯

a Widthofthethinfilm R_p Equivalentparticleradius

h¯ Characteristicmeasureoffilmthickness

¯

r Characteristicmeasureoffilmwidth

ε

1,

ε

2,

^{Lengthscaleratios}

p Excesspressureinthefilm

v_r r-componentofthecontinuousphasevelocity vz z-componentofthecontinuousphasevelocity v_d Dispersedphasevelocityfield

p_d Dispersedphasepressure

U Tangentialvelocityoftheinterface

τ

d Particlesidetangentialstressevaluatedatthein- terface

A Hamakerconstant

σ

^{Surfacetension}

r_∞ Alargeradialdistance

I BoundaryIntegralMethodintegrand

ρ

^{Integralvariableinradialaxis}

θ

^Azimuthangle

λ

^{∗ Dispersed to continuous phase viscosity ratio}

multipliedby

A^∗ DimensionlessHamakerconstant Ca Capillarynumber

Uˆ Tangentialvelocityoftheinterfacemultipliedby

λ

^∗

t¯ Characteristictimescale

Y₁−Y₄ Transformed variablesemployed inspectral ele- mentbasedsolvers

δ

^{Asmallarbitrarynumber}

R Radial distance in the polar coordinate system usedforsingularitytreatment

ϕ

^{Azimuth angle in the polar coordinate system}

usedforsingularitytreatment

[A] IntegrationmatrixusedfortheBoundaryIntegral equation

V_dimp Critical approach velocity separating the linear andthedimpleddrainageregimes

V_mult Criticalapproachvelocityseparatingthedimpled andthemultiple-rimdrainageregimes

Iav AveragedBoundaryIntegralMethodintegrand E TrapezoidalareaunderZ

τ

dbetween[r−

δ

,r+

δ

^]

Z Integrationkernel v Arbitraryscalarquantity

v

ˆ vinChebyshevpseudospectrum

N Numberofcollocationpointsminusone T Matrixrelatingaquantityanditscounterparton

theChebyshevpseudospectrumbasedonCheby- shevpolynomials

b Numericalintegrationvector [B] Numericalintegrationmatrix

[C₁]−[C₃] Integrationmatricesemployedtocompute[A]

considering the changes in the kinetic and the surface energies of the particles, and by Kirkpatrick and Lockett (1974) via a film drainage model. Although these studies aim to explain the phenomenon based ondifferent arguments,they are not entirely conflicting, but rather might represent different regions of the particle approachvelocity spectrum. Yaminskyet al.(2010) iden- tifies three distinct regimes of coalescence in their experiments withairbubblesinwater:atverylowapproachvelocities(smaller than 1 μm/s), the thin film in between bubbles appears to be stable in finite time, at intermediate velocities (upto 150 μm/s) coalescenceoccursintherangeof10–100swithvisibledimpling of the interface, andfinally atlarge velocities (much larger than 150 μm/s) almost immediate coalescence is observed with very little hydrodynamic resistance from the film and no observable dimpling. Theyassociate these three regions with DLVO stability meaning that thestabilitystems fromthe electrical forcesacting between charged interfaces, viscous film drainage in which the filmflattensanddimples,andinertialfilmdrainagewherecoales- cenceoccursbefore thedeformationofthe particles,respectively.

Horn et al. (2011) discusses several works from the literature (Chesters and Hofman, 1982; Klaseboer et al., 2000; Yaminsky et al., 2010; Del Castillo et al., 2011) to end up with a single chart ofdifferentpossibleregimesofcoalescenceasafunctionof the bubble approach speed and the NaCl concentration in their Figs.2and3.Fromtheirmap,itisevidentthatinthelimitofzero NaClconcentration,astheapproachspeed(orthekineticcollision energy)increasessequentialregimesof

1. DLVOstabilityforverylowvelocitycollisions,wherenocoales- cenceispossibleinfinitetime;

2. Film drainage governed by viscous effects for collisions with higher velocity, where time of coalescence might be as large as100s;

3. Inertia governed rapid film drainage for energetic collisions, wherecoalescenceisvirtuallyimmediate;

4. Bouncingregime forevenmoreenergeticcollisions, whereco- alescencedoesnotoccur;

are expected. However, they do not provide information on how the coalescence time changes within these individual regimes.

Orvalho etal.(2015) observesa power-lawtyperelationbetween the coalescence time, t_c,and therelative approach velocity, V_app, in their experiments with air bubbles and liquids of different viscosities. Theyproposethat tc is alinearfunction ofV_app^−0.85.On the other hand,Del Castilloet al. (2011)observes similar trends in their experiments with air bubbles in pure water, with the additional observation of a minimum in Vapp versus tc as the approachvelocityfurtherincreases.

Inthefilmdrainagemodels,oncethefluidparticlesarebrought intocontactbytheexternalflow, threeconsecutivestepsarepro- posed(ShinnarandChurch,1960)

• Athin film of the continuousphase is entrapped inbetween twocollidingparticles,

• Thefilmdrainsuntilitsthicknessreachesacriticalvalue,

• The film ruptures and coalescence occurs upon reaching the criticalthickness.

Iftheparticlesareallowedtobeincontactbytheexternalflow forsufficientlylongtime,suchthat thetimerequiredforthefilm todrainuntilthecriticalthicknessisreached,issmallerthanthe duration oftheinteraction betweentheparticles,coalescence oc- curs.Basedonthefilmdrainageapproach,Coulaloglou(1975)gives a statistical formula forthe coalescence efficiency,

λ

^{c, using two}

characteristictimescales,thedrainageandthecontacttimes,as

λ

^c=exp

−tdrainage

tcontact

(1)

whichisfrequentlyemployedtoestimate coalescencebehaviorin the scale of the multi-phase equipment such as chemical reac- torsandseparators.Themodelisrestrictedto’gentle’collisions,in whichtheparticlesaresignificantlylargerthantheradiusaffected bythecollision(orinotherwords,thewidthoftheemergingthin film)(Chesters,1991).Theenergymodels,ontheother hand,con- sidermore’energetic’collisionsandrelates

λ

^{c tothekineticcolli-}

sionenergyandtheinterfacialenergyofthefluidparticleinstead oft_drainageandtcontact.

To determinethe drainagetime, t_drainage,in Eq.(1),numerous modelswithvaryingdegree ofcomplexityhavebeenproposedin the literature.These modelsare classifiedbased onthe deforma- bilityandthetangentialmobilityoftheinterfaces(LeeandHodg- son,1968; Chesters,1991;LiaoandLucas,2010). Theformerclas- sification differentiates between the models for the rigid spheri- cal fluid particles andtheones fortheparticles withdeformable interfaces,which are capableofmodeling the experimentally ob- served dimple formationduring thefilm drainage. Inaddition to theformationofthedimple,thepimpling,wimplingandrippling of the interface are also possible (Chanet al., 2011). The pimple and the wimple refer to the emergence of an additional rim at the interface, andof a local maximumbetween the rimandthe centeroftheinterface, respectively.Whereas, thedevelopmentof multipleadditionalmaximaandminimaattheinterfaceistermed as rippling. Next,by following Chesters (1991) and Liaoand Lu- cas (2010),thedeformableclasscanbe furtherdividedintothree subgroupsbasedontheinterfacialmobility,asimmobile,partially mobile andfully mobile models. Immobilemodels correspond to theinterfaceswithzerotangentialvelocity,whichmightstemfrom veryhighdispersedphaseviscosityand/or interfacialtensiongra- dientsemergingduetopresenceofimpuritiesandsurfactants(Lee andHodgson, 1968; Chesters,1991). Inthat case, thedrainageof thefilmisgovernedbyviscousforceswithinthefilm(resultingin a parabolic velocity profile) andnot coupledto the velocity field of the dispersed phase (Liao andLucas, 2010; Chanet al., 2011).

Partially mobile and fully mobile names, on the other hand, are usedforthecaseswherethecontributionofthetangentialveloc- ityof theinterfaceson thefilm drainageare non-negligible.This non-negligibletangentialvelocityofthemobileinterfacesyieldsa plug flow-likecontribution to the velocity ofthe film. According to Chesters (1991),for the fully mobile interfaces, the tangential stress atthe interface iseffectivelyzeroandthe drainageisgov- erned eitherby thedeformation oracceleration,respectively cor- responding to theviscous andthe inertial filmdrainage regimes.

Whereasinthepartiallymobileinterfacecasethetangentialstress attheinterface isnon-zero andthe drainageiscontrolled bythe dispersedphaseviscosity.Basedonscalingargumentsonthethin filminbetweenthefluidparticles,Davisetal.(1989)arguesthat if the dispersed to continuous phase viscosity ratio,

λ

^{, is much}

greater than

a/h₀, wherea= a^a₁¹+^aa²₂ is the reduced particle ra- dius, a₁ and a₂ are particle radii, and h₀ is the minimum initial thickness of the film, the fluid particles act as rigid particles, if

λ

<<

a/h₀ the interfaces of the particles are fully mobile, and if

λ

^and

a/h₀ are oncomparablemagnitudes,theinterfacesare partiallymobile.

Modeling of the interface as a mobile one, regardless of its deformability, requires the coupling of the flow fields inside the fluid particle and within the film. Davis et al. (1989) employs the Boundary Integral theory to determine the tangential veloc- ityofthe interface,which, then,enablesthe couplingoftwo ve- locityfields through the no-slipcondition atthe interface, with- out requiring the computation of the velocity profile inside the fluid particle. This approximation lowers the computational ef- forts significantly and later adapted by many others (Yiantsios and Davis, 1990; 1991; Abid and Chesters, 1994; Saboni et al., 1995;Klaseboeretal.,2000;Bazhlekovetal.,2000).Yiantsiosand Davis (1990) concludes that a dimple is always formed during bouyancy-driven interactions between a fluid particle and a sur- face regardless of the viscosity ratio. In their following work (YiantsiosandDavis,1991),theyextendtheirworktotheinterac- tions betweentwofluid particles,in thelimitsofvery smalland very large interfacial mobilities, and reveal that without the at- tractivevan der Waals forcesthe coalescence is not possible for deformableinterfacesinfinitetime. Theyalsoidentifytwo differ- entrupturingtypes,thenoseandtherimruptures,whichareob- served at strong andweak van der Waals forces, respectively. In thenoserupture,thevanderWaalsforcesbecomesignificantbe- foredimplingoccursresultinginruptureatthecenteroftheinter- face,whereastherelativelyweakervanderWaalsforcesallowthe capillaryforcestoactfirst,resultinginformationofadimpleand a rimaway fromthe axisof symmetry. Inthat case, the rupture startsfromthe rimearning the namethe rim rupture.Abid and Chesters(1994)andSabonietal.(1995)studycenterlinecollisions of droplets in the presence of van der Waals forces, assuming a constantapproachvelocityandaconstantdropletinteractionforce, respectively.Theydeterminethe criticalfilmrupturethicknessas afunctionoftheHamakerconstant.Intheirmodels,thedropvis- cosity determines therateofdrainage, meaning that thevelocity profilesareapproximatedasplug flows,byneglectingtheviscous effectswithin the film.Klaseboer etal.(2000) comparetheir ex- perimentalresultsobtained forthecoalescence ofvarious liquid- liquidcouplestooutcomesoftwotheoreticalmodels,inwhichthe interfaces are deformable andthe droplets havea constant rela- tiveapproachvelocity.Theyconcludethat thetangentiallyimmo- bile interface model,where the velocity in the film is parabolic, fitstheir experimentaldatamoreaccurately incomparisonto the mobileinterfacemodel,wheretheyassumeaplug-likeflowwithin thefilm.Theysuspectthatthepossibilityofverysmallamountof surfactantsbeingpresentintheexperimentsmightberesponsible fortheimmobilizationoftheinterface.Bazhlekovetal.(2000)in- vestigatethe effectof theviscosity ratioon thefilm drainageby consideringeitherconstant approachvelocityorconstant interac- tionforcecollisionsoftwofluid particles.Theirmodeltakesboth theparabolicandplug-likecontributions tothefilmflowintoac- count,thereby closingthegapbetweenpreexisting immobileand mobilemodels inthe literature. Theypresentestimations forthe minimumfilmthickness asafunction oftime andviscosityratio inacompactform.However, intheir expressionstheeffectofthe vanderWaalsforcesisnotconsidered.

Our models follow those of Bazhlekov et al. (2000)’s and Klaseboeretal. (2000)’s,forthe mobile andimmobileinterfaces, respectively, withthe addition of van der Waals effects. The ad- dition ofthe van der Waals forcesto the models makes the es- timation of the coalescence time possible. Furthermore, we ren- derourequationsdimensionlessfollowinga differentroute toset therelativeapproachvelocity oftheparticlesasamodelparame- ter,in orderto investigate its effecton the coalescence time ex- plicitly, together with the effects of the disperse to continuous phase viscosity ratio andHamaker constant. Although our mod- elsaresimilartothepreexisting onesintheliterature,theability toshow theindividual effect oftheapproach velocityon theco-

Fig. 1. Physical system.

alescencetime directly,enables usto compareourtheoretical re- sultstotheobservationsfromnumerousexperimentalworkssuch asthose of Yaminsky et al. (2010), Del Castillo et al.(2011) and Orvalhoetal.(2015),inanefficientandclearmanner,therebyal- lowingustoilluminatesomeofthemissinglinksbetweenthethe- oreticalandexperimentalworksintheliterature.

Thearticleisoutlinedasfollows:Thephysicalconfigurationof interest and the mathematical model is presented in Section 2. Thenumericalprocedure includingthetreatmentofthe singular- ityinherenttotheBoundaryIntegralMethodisgiveninSection3. Section4presentsthecoalescencetime andthetimeevolutionof thefilmthicknessatdifferentvaluesofthephysicalparametersof interest,andthecomparisonofthemodels’outcometotheexperi- mentalobservationsintheliterature.Finally,theconclusionsofthe workaresummarizedinSection5.

2. Physicalsystemandmathematicalmodel

This work considers the axisymmetrical interactions between twoNewtonianfluidparticles(eitherdropletsorbubbles)ofsame fluid approaching each other at a constant relative approach ve- locity,Vapp,alongtheir centerlines,throughacontinuousmedium ofa Newtonianfluid.Thedepictionofthephysicalsystemispre- sentedinFig.1.Theviscositiesofthecontinuousandthedispersed phasesaredenoted by

μ

^{c and}

μ

d,respectively,andthefluidpar- ticlesareallowedtobeofdifferentradii,R₁andR₂.Theinterfaces betweenthedispersedandthecontinuousphasesareimmiscible, deformableandcharacterized by a constantvalue of surfaceten- sion,

σ

^{. Asthe fluid particlesapproach each other,they entrapa}

thinfilmofthe continuousphasein betweenthem,which even- tuallystartsto drain.Thethicknessof thisemergentthinfilm,h, isafunctionofthe radialpositionrandthetime t.It ispossible toidentifythreedistinctlengthsgoverningthefilmdrainagephe- nomenon:theparticleradiiR₁andR₂,thethicknessofthefilmh(r, t) andthewidthofthefilm a¯. Althoughthecharacteristiclength scalesforthe firsttwo are straightforward,the widthofthefilm isnot easytomeasure.Thereforeitsmagnitudewillberelatedto thelengthscales fortheparticlesizeandthefilmthicknesslater onduringthederivationofthemodel.Providedthatthecollision isa gentle one,meaning that both particleradii are muchlarger than the characteristic film width, the equivalent radius, defined by

1 Rp =1

2

₁

R1 + 1 R2

(2)

canbeusedasthelengthscaleforbothparticles(Chesters,1991).

Duetothisapproximation,thecollisionofunequal-sizedparticles canbemodeledasthecollisionofequal-sizedones,meaningthat symmetry aroundr axisisalso introducedto themodel inaddi- tion to theaxisymmetry. This newsymmetry creates fourequiv- alentquadrants, enablingusto seekforthe solution onlyin one ofthem.Here,onlythesolutioninr≥0,z≥0quadrantissought, where the interface position is given by z=h(^r,t)/2. The gentle collision resultinginemergence ofathinfilm suggestsarelation betweenthethreelengthscales:

h¯ <<a¯<<Rp (3)

whereh¯ anda¯aremeasuresofthethicknessandthewidthofthe film,respectively.Inthelightofthisargumentonthelengthscales, twodimensionlessratiosaredefinedas

h¯

¯

a=

ε

1<<1, a¯

Rp =

ε

2<<1 (4)

whichfurthersuggeststhatasmallnumber

^{canbedefinedas}

_¯

h Rp

1/2

=√

ε

1

ε

2=

⁽⁵⁾

Then,

^{canbeutilizedtorelatethethreedistinctlengthscalesin}

themodelas h¯ =

^2Rp, r¯=

ε

1

ε

2

¯

a=

^{Rp (6)}

where,r¯isemployedinordertounifytheorderofmagnitudedif- ference betweenthe length scales andhereafter, adopted as the r-directionlengthscaleinsteadofa¯.

2.1. Governingequationsandboundaryconditions

Based on the arguments in Section 2 (h¯<<r¯) the lubrication theory is applicable within the film. Following the analysis of YiantsiosandDavis(1990),as

=h¯/r¯<<1,theflowofthefilmis governed bythe simplifiedversionsofthe Navier–Stokesandthe continuityequations

∂

^p

∂

^{r =}

μ

^c

∂

²

v

r

∂

^z2,

∂

^p

∂

^{z =0 (7)}

and

∂v

z

∂

^{z +1r}

∂

∂

^r

⁽

^r

v

r

)

^{=0 (8)}

Fig. 2. Sketch of the coordinate transformation from ( ρ, θ) to ( R , ϕ) and the integration domain. Labels 1, 2 and 3 stand for the integration subdomains that have to be handled separately.

Fig. 3. Fig. 4 of Klaseboer et al. (20 0 0) reproduced for validation of the immobile solver. Thin film thickness as a function of r and t as λ^∗→ ∞ . Here, to match Klaseboer et al. (20 0 0) , t is shifted by −17 to set the onset of dimpling as t = 0 . Profiles are obtained with A ^∗= 0 , V app = 1 , h 00= 10 and r ∞= 15 .

wherepistheexcesspressure inthefilm,and,vz andvrarethe z and r components of the continuous phase velocity. The flow withintheparticles,ontheotherhand,isapproximatedbythein- compressibleStokesequations

μ

d

∇

^2vd=

∇

^pd (9)

and

∇

^·vd=0 (10)

where p_d and v_d are the dispersed phase pressure and velocity fields,respectively.Theflowsinbothphasesarecoupledviaaset ofboundaryconditionsvalidattheinterface,z=h(^r,t)/2:theno- slipcondition,thekinematiccondition,andthetangentialandthe normalcomponentsofthestressbalancegivenby

v

r

|

z=h/2=U (11)

1 2

∂

^h

∂

^{t =}

v

z

|

z=h/2−1 2

∂

^h

∂

^r

v

r

|

z=h/2 (12)

−

μ

^c

∂ v

r

∂

^z

z=h/2

=

τ

d (13)

p=2

σ

Rp −

σ

₂¹_r

∂

∂

^r

r

∂

^h

∂

^r

+ A

6

π

^{h3 (14)}

respectively.Thekinematiccondition,Eq.(12),appearsduetothe continuitybetweenthe normal componentof thevelocity of the thinfilm andthe normal speed of the interface (given by ¹₂^∂_∂^h_t).

Here,Uisthetangentialvelocityoftheinterface,

τ

distheparticle- sidetangentialstressevaluatedattheinterface,A istheHamaker constant,and thedeviations inthe particle-side pressureare ne- glectedinEq.(14)duetothegentlecollisionassumption.Further- more,itisassumedthat ata largeenough radialdistance,atr_∞, theshapeoftheinterfaceandtheapproachvelocityareunaffected

bythecollision.Then,theconditions, p

|

r=r∞=0 and

∂

^h

∂

^t

r=r∞

=−Vapp (15)

hold. Eq. (15) form the set of boundary conditions for p and h together with the symmetry conditions at r=0. The initial film thicknessistakenas

h=h00+ r² Rp

(16)

to resemble the distance between two spherical fluid parti- cles. To determine U, the boundary integral form of the Stokes flow is used (reader may refer to Davis et al. (1989) and Ladyzhenskaya(1969)forfurtherinformationonthemethod.) U= 1

μ

d r_∞ 0

I

( ρ

^,

θ ) τ

dd

ρ

⁽¹⁷⁾

where I

( ρ

,

θ )

=

ρ

2

π

π 0

cos

θ

r²+

ρ

^2−2r

ρ

^cos

θ

^d

θ

⁽¹⁸⁾

2.2.Dimensionlessequations

Byapplyingthetransformations h˜ = h

^{2Rp, ˜r=}

r

^Rp,

v

^˜r=

v

r

μ

^c

³

σ

^,

v

^˜z=

v

z

μ

^c

⁴

σ

^,

V˜app=Vapp

μ

^c

⁴

σ

^{, p˜=}

pRp

σ

^,

τ

^˜d=

τ

dRp

σ

^{, t˜=}

t

σ μ

dRp

(19)

theEqs.(7)–(8)and(11)–(16)arerendereddimensionlessas

∂

^p˜

∂

^r˜=

∂

²

v

^˜r

∂

^{z˜2 ,}

∂

^p˜

∂

^{z˜=0 (20)}

∂ v

^˜z

∂

^{z˜ +}

1

˜ r

∂ ∂

^r˜

⁽

^r˜

v

^˜r

)

^{=0 (21)}

v

˜r

|

z˜=˜h/2=U˜ (22)

1 2

λ

^∗

∂

^h˜

∂

^{t˜ =}

v

^˜z

|

˜z=˜h/2−1 2

∂

^h˜

∂

^r

v

^˜r

|

z˜=h˜/2 (23)

−

∂ v

^˜r

∂

^z˜

˜ z=h˜/2

=

τ

˜d (24)

˜

p=2− 1 2r˜

∂ ∂

^r˜

˜ r

∂

^h˜

∂

^r˜

+A^∗

h˜³ (25)

˜

p

|

r˜=r˜_∞=0,

∂

^h˜

∂

^t˜

˜ r=r˜∞

=−

λ

^{∗V˜app (26)}

and

˜h=˜h00+r˜² (27)

whereA^∗= _6π^A_6R2

pσ and

λ

^∗=^μ_μ^d_c.The non-dimensionalized nor- malstress balancealsorequiresarestrictiononthemagnitudeof thecapillarynumber,Ca= ^μcV_σ^app ≈

⁴, asdiscussed forthinfilm analysesin generalbyOron etal.(1997).Onthe other hand,the boundaryintegralequation,Eq.(17),becomes

U˜= 1

λ

^∗

˜ r∞ 0

I˜

( ρ

^˜,

θ ) τ

^˜dd

ρ

˜= 1

λ

^{∗Uˆ (28)}

andsubstitutingr=˜rand

ρ

=

ρ

˜ ^{intoEq.(18)givesI˜}(

ρ

^˜,

θ

)^.Here- after, all equations are given in dimensionless form; therefore,

tildes are omittedin the transformed variables. Solving Eq.(20), and,applyingEq.(22)atz=h/2togetherwiththesymmetrycon- ditionatraxis,yields

v

r=1 2

∂

^p

∂

^r

z²−

h 2

2

+U (29)

clearlyindicatingtwo distinct contributionstothe radialvelocity, theparabolicflowdriven bythepressuregradientinside thefilm andtheinterfacially-drivenplugflow.Then,thetangentialcompo- nentofthestressbalance,Eq.(24),becomes

−h 2

∂

^p

∂

^{r =}

τ

d (30)

Oncevrisknown,thez componentofthefilmvelocity, vz,isde- terminedvia Eq.(21).Finally,by substituting vz andvr (given by Eq.(29)intoEq.(23),thethinningequationisfound:

∂

^h

∂

^{t =}

1 r

λ

^∗

12

∂ ∂

^r

r

∂

^p

∂

^rh3

−

∂

∂

^r

rUˆh

(31)

Despitethe differences inthe characteristic scales employed, the equationisexactlythesameastheBazhlekovetal.(2000)’sthin- ningequation,andrevealsthattheeffectsoftheparabolicandthe plug flows inthefilm are weightedbythe dispersedto thecon- tinuousphaseviscosityratiovia

λ

^{∗.Thatbringsouttheclassifica-}

tionoftheinterfacesintermsoftheirmobility.Atveryhighvalues of

λ

^{∗, theinterface attains}infinitesimalvaluesofUˆ andbecomes tangentiallyimmobile,whereasatlow

λ

^{∗thethinningequationis}

dominatedbytheplugflow’scontribution,earningthetitleoffully mobileinterfacefollowingDavisetal.(1989)’sdescription.Inthese extremecases, Eq.(31) can be simplified by onlytaking intoac- count the dominantterms forthe respectivecaseand neglecting theother.Then,fortheimmobilecase

∂

^h

∂

^{t =}

1 12r

∂ ∂

^r

r

∂

^p

∂

^rh3

(32)

andforthefullymobilecase

∂

^h

∂

^{t =−}

1 r

∂ ∂

^r

rUˆh

(33)

canbeusedasthethinningequation.Noticethat,intheimmobile case,

λ

^{∗ is included inthe time scale by setting t¯}=^μc2^Rσ^p, which further impliesthatthe righthand sideof thevelocity boundary condition in Eq. (26) becomes −V_app. The thinning equation for thefully mobilecase, Eq.(33), canbe obtainedby setting

λ

^∗=0 in Eq.(31). Therefore, thefully mobile case isreferred to asthe

λ

^∗=0 case, hereafter. However, this only implies that the first termontherighthandsideofEq.(31)isnegligible,anddoesnot suggest thephysical value of

λ

^{∗ to beactually zero.}Additionally, even though the name ‘the

λ

^∗=0 case’ used for the fully mo- bile case evokes the opposite, the

λ

^{∗ appearing in the boundary}

condition, Eq. (26), is not zero, but it is rather a small physical value. Onthe other hand, atmoderatevalues of

λ

^∗ (corresponds toDavisetal.(1989)’spartiallymobile interface),contributionsof bothflows havenon-negligibleeffectsonthinning andEq.(31)is directlyemployedwithoutanyfurthersimplifications.Inallcases, thethinningequationiscoupledwithEq.(25),andexcepttheim- mobilecase,Eqs.(28)and(30)areusedtodetermineUˆ.

3. Numericalprocedure

Foreach oneofthethree mobilitycases(the partiallymobile, theimmobileandthefullymobileinterfaces)giveninSection2.2, adifferentsolverisemployed.However, thesesolvers sharesome characteristics: In all cases, the corresponding thinning equation (Eqs.(31),(32)or(33))issolvedsimultaneouslywithEq.(25).The

boundaryconditionsgiveninEq.(26)holdtogetherwiththesym- metry conditionsforbothhandpatr=0.Thespatialdiscretiza- tion schemeiseithera spectral onebased onChebyshevpolyno- mials,orconsists ofspectral elements, whichare againbasedon Chebyshev polynomials. The elementsare used in some cases to obtain higherspatialresolutionwithout drasticallyincreasing the computational efforts.However, whenspectralelementsareused, since Chebyshev polynomials are only c₀ continuous on the ele- ment boundaries (Patera, 1984; Karniadakis and Sherwin, 2013), the solver fails to give accurate results for differential equations withorderhigherthan one.Therefore,thesecond orderdifferen- tial equations, the thinning equation (one of Eqs. (31)–(33)) and Eq. (25), are converted into a set of four first order differential equations in the spectral element solvers. The transformation is done by setting Y₁=h, Y₂=^∂_∂^h_r, Y₃=p, Y₄= ^∂_∂^p_r. In all solvers, thesecondorderbackwarddifferentiationisusedtodiscretizethe timederivatives.

In mobile cases,computation ofUˆ isrequired via the Bound- aryIntegralequation,Eq.(28).However,theequationhasaninher- entsingularityat

ρ

=rand

θ

=0ascanbeseenfromEqs.(17)to (18).Tothebestofourknowledge,theworksintheliteraturethat simulate similar thinning equations to ours, donot explicitlyre- veal their methods of coping withthe singularity. Therefore, we attempttogive someinsightonourmethodinthissection.First, startingwithEq.(28),theintegraldefiningUˆisdividedintothree partstoisolatethesingularity:

Uˆ = ^r∞

0

I

( ρ

,

θ ) τ

dd

ρ

= ^r−δ

0

I

( ρ

,

θ ) τ

dd

ρ

+ ^r+δ

r−δ I

( ρ

,

θ ) τ

dd

ρ

+ ^r∞

r+δI

( ρ

,

θ ) τ

dd

ρ

⁽³⁴⁾

where

δ

^{isa smallarbitrarynumber. Here,the} singularityis iso- lated inthe second term on theright hand side ofEq. (34), en- ablingefficientandaccurate computationofthe firstandthelast terms. The kernels in theseterms, I(

ρ

^,

θ

^{), are merely geometric}

functions,andasaresult,computedonlyonceatthebeginningof thecomputation.Thesingularintegral,however,istreatedfurther.

Following Section8.2ofScottetal.(2013),acoordinate transfor- mation(sketchedinFig.2)isappliedtoshiftthesingularpointto theoriginofapolarcoordinatesystem,from(

ρ

^,

θ

^)to(R,

ϕ

^)via

ρ

^−r=Rsin

ϕ

^,

θ

^=Rcos

ϕ

⁽³⁵⁾

As can be seenfromFig. 2,the transformationrequiresthe inte- grationdomaintobehandledinthreeseparateparts.Thesingular integral,then,iswrittenintermsofthesummationofthesethree integralsas

r+δ

r−δ I

( ρ

^,

θ ) τ

d

( ρ )

^d

ρ

= 3

k=1

_R

k

0 ϕk+1

ϕk

I

(

^R,

ϕ ) τ

d

(

^R,

ϕ )

^Rd

ϕ

^dR

(36)

where R_k=

−

δ

sin

ϕ

^,

π

cos

ϕ

^,

δ

sin

ϕ

,

ϕ

k=

−

π

2 ,tan⁻¹

−

δ π

, tan⁻¹

δ

π

,

π

2

and

I

(

^R,

ϕ )

⁼

(

^r+Rsin

ϕ )

^cos

(

^Rcos

ϕ )

r²+

(

^r+Rsin

ϕ )

^2−2r

(

^r+Rsin

ϕ )

^cos

(

^Rcos

ϕ )

⁽³⁷⁾

Next, by assuming

τ

d is constant in the

δ

^{vicinity of r}=

ρ

, Eq.(34)turnsouttobe

Uˆ= ^r−δ

0

I

( ρ

^,

θ ) τ

dd

ρ

⁺

τ

d

3

k=1

_R

k

0 ϕk+1

ϕk

I

(

^R,

ϕ )

^Rd

ϕ

^dR

+ ^r∞

r+δI

( ρ

^,

θ ) τ

dd

ρ

⁽³⁸⁾

in which, now, the singular part of the integral is also given as a purely geometrical function, that is computed only once throughoutthesimulation.Finally,alltheintegrals,inEq.(38),are transformedinto matrixform usingquadraturerules, so that the BoundaryIntegralequationcanbewrittenas

Uˆ=[A]

τ

d (39)

where[A]istheintegrationmatrix.Somekeydetailsinderivation of[A]isprovidedintheAppendixA.Theaccuracyoftheapproxi- mationmightdependheavilyonthevalueofthesmallparameter

δ

^{whenitisselectedpoorly.Inour}simulationsthisvalue istaken aslowas10⁻⁴-10⁻⁵,aroundwhichnofurthersignificantdepen- denceofUˆonthevalueof

δ

^{isobserved.Formobilecases(except}

for

λ

^∗=0),first,theEqs.(25)and(31)aresolved simultaneously atagiventimestep;then,throughEqs.(30)and(39),

τ

dandUˆare determined;andfinally,Uˆ isfedtothenexttimestepasaninput.

Whereas,inthefullymobilecase(

λ

^∗=0),Eqs.(25),(30),(33)and (39)aresolved simultaneouslyto obtaina morestablenumerical scheme,whichisfoundtobeamust,especiallyintheearlystages ofthesimulation.

4. Resultsanddiscussion

Althoughthemodelsinthisworkarenotcapableofrendering theruptureoftheinterface,itispossibletoforeseeitbyconsider- ingnon-zeroA^∗values.Atacriticalfilmthicknessthatdependson themagnitude ofA^∗,the attractivevan derWaals forces become muchmoreinfluentialthantheresistancewithinthefilm.Thisre- sultsinveryrapidthinningofthefilm,afterwhichtheruptureof theinterfaceisextremelylikely.Sincethetimescaleoftherupture ismuchsmallerthanthetimescaleofthefilmdrainageuntilthe criticalrupturethicknessisreached,thecoalescencetimeisdeter- minedby approximatingitasthedrainagetime,withoutactually modelingthe rupturingphenomenon itself.Dueto thisapproach, the value given asthe dimensionless coalescence time, tc, is the simulationtime between thebeginning of thenumerical simula- tion and the time right before the rupture in the corresponding simulation.

4.1. Immobileinterfaces

Intheir Fig.4,Klaseboer etal.(2000) presentsthe filmthick- nessprofilesatdifferenttimesobtainedfortheirimmobilemodel withthe constant approach velocity boundary condition. Tovali- datetheimmobile solverdescribedinSection 3,their results are reproducedbysettingA^∗=0andVapp=1,andgiveninFig.3.The results of our immobile solver seem to be in perfect agreement withthose ofKlaseboeretal.(2000)’s.Fort>0,theformationof thedimpleis clearlyseenin Fig.3.However, duetothe absence ofthe disjoiningpressure (A^∗=0),the ruptureof thefilm is not estimated.

Fig. 4 presents the coalescence time, tc, as a function of the relative approach velocity for different non-zero values of A^∗. For all the values of A^∗ considered here, three distinct types of drainage/coalescence behavior are observed. First, at very low approachvelocities,tc decreases withVapp followinga powerlaw type relation,i.e., alinear relationbetweenlog(^tc)^andlog(^Vapp)^. Anexampleofthe‘lowvelocitydrainage’canbeseeninFig.5(a),

Fig. 4. Coalescence time for immobile interfaces as a function of the relative ap- proach velocity. Profiles are obtained with h 00= 10 and r _∞= 15 .

inwhichthetimeevolutionofthefilmthicknessforVapp=0.003 andA^∗=10⁻⁴ is shown. In that type of drainage, the attractive vanderWaalsforcesbecomesignificantbeforethecapillaryforces act substantially, resulting in rupture at r=0, the center of the fluid particle (nose rupture). The low velocity drainage behavior continues until a critical velocity is reached. After this critical velocity,thecapillaryforcesbecomeinfluentialbeforetherupture

occurs, and as a result, a dimple shape emerges. Hereafter, this critical velocity is referred to as the dimpling velocity, V_dimp. The results indicate that V_dimp increases as the attractive forces gets stronger meaning that stronger capillaryforces are required for the dimple formation as A^∗ increases. Fig. 5(b) reveals the typical drainageprocess (obtainedforVapp=0.01 andA^∗=10⁻⁴) inthisregime.Although,firsttwoprofiles aresimilartotheones observed in the low velocity drainage regime, in later stages of the drainage(aftert>1040) the dimplebecomesvisible,andthe rupture occurs on the rim instead of the center (rim rupture).

In thisregime, tc continues todecrease withincreasing Vapp,but the decrease becomes lessand lessdramatical as Vapp increases, and eventually, tc passes through a minimum. Right after the minimum, further increase in the velocity brings out the third regime,wheremultiplerim-like structures emergebeforerupture andtcincreaseswithVapp.Thissecondcriticalvelocityisdenoted as V_mult. Two examples of the time evolution of the film thick- ness in the ‘multiple-rim’ regime are provided in Fig. 5(c) and (d), which fall into the category of pimpling andrippling of the interface,respectively.Here,thepimplereferstotheemergenceof anadditionallocalminimumsuchthatthemainrimispositioned inbetweenthecenterofthefluidparticleandthenewlocalmin- imum,whereas thetermrippleis usedwhen multipleadditional minimaandmaximaareencountered(ReadermayrefertoSection 3.1.2ofChanetal.(2011)formoredetailandvisualizationofthe interfaceshapes).Theslope oflog(^Vapp)^versuslog(^tc)^inthefirst region,andthe criticalvelocities,V_dimp andV_mult,can beused to identifyeach curve inFig.4. Thecorresponding valuesare deter- mined as: d(^log(^tc))/d(^log(^Vapp))=[−0.994,−0.992,−0.988], V_dimp=[0.2,0.05,0.01] and V_mult=[0.9,0.3,0.09], for A^∗=

Fig. 5. Time evolution of the film thickness in different regimes: (a) the low velocity regime, the nose rupture, (b) the dimpled drainage regime, the rim rupture, (c) and (d) the multiple rim regime (pimpling and rippling behaviors, respectively). Profiles are obtained with A ^∗= 10 ⁻⁴, h 00= 10 and r _∞= 15 .

Fig. 6. Fig. 4 of Bazhlekov et al. (20 0 0) reproduced for validation of the mobile solvers. Thin film thickness as a function of r and t for (a) λ^∗= 0 , (b) λ^∗= 10 , (c) λ^∗= 100 and (d) λ^{∗= 10 0 0 . In all cases A ∗= 0 ,} λ^∗V app= 1 , h 00= 2 and r _∞= 30 .

[10⁻²,10⁻³,10⁻⁴], implying that, just as Vapp does, V_mult also increases with increasing A^∗, whereas the slope decreases. The analysis iscarried out by assuming that the interface immobility is due to the very high dispersed to continuous phase viscosity ratio, however the applicability of the observations should be independent ofthesource oftheinterfacialimmobilization,asin alltypesofimmobilizationthetangentialvelocityoftheinterface, Uˆ,iszero.

4.2. Mobileinterfaces

To validate the mobile solvers describedin Section 3,the re- sults of Bazhlekov et al. (2000) for the constant approach ve- locity boundary condition are reproduced by setting A^∗=0 and

λ

^∗Vapp=1.Fig.6correspondstotheirFig.4andpresentstheevo- lutionofthefilm thicknessasafunction oftimeandradial posi- tionatdifferent

λ

^{∗.Forallvaluesof}

λ

^{∗presentedintheFig.6dim-}

pleformationattheinterface isobserved.Althoughitisnotpos- sible to compare the numerical values for the film thickness di- rectly, excellent visual agreement is observed for partially mo- bilecases

λ

^∗=10,

λ

^∗=100 and

λ

^∗=1000.For

λ

^∗=0,however, our solver doesnot produce the very sharp gradients shown by Bazhlekovetal.(2000)att=96,althoughattheearliertimesthe results seem to be in agreement. This discrepancymight appear duetothedifferentnumericaltechniquesemployed,i.e.,finitedif- ference in Bazhlekov et al. (2000) and spectral methods in this work, ordueto thedifferencesin treatmentof theboundaryin- tegral. Here, we cannot give anyfurther discussion on thelatter

possibility,sinceBazhlekovetal.(2000)doesnotprovideanyde- tailedinformationontheirtreatment.

As mentioned earlier in the end of Section 2.2, the

λ

^∗=0 caserefers to the fullymobile interfaces. In this case, the value of

λ

^{∗ is sufficiently small (yet non-zero) enough to render the}

first termon the right hand side of Eq.(31) negligible, enabling usto approximate the drainagevia Eq. (33)by setting

λ

^∗=0 in Eq.(31). However, in the boundary condition, Eq.(26),

λ

^{∗ isnot}

taken as zero, but rather as the small physical value itself, im- plyingthat

λ

^∗Vapp=0.Then,aviscosityratiosweepiscarriedout todetermine the

λ

^{∗ value, belowwhich the drainagecan beap-}

proximatedwiththe

λ

^∗=0 case. Tocarry out the viscosityratio sweep,first, the results forthe very low viscosity ratiolimit are obtained by solving the

λ

^∗=0 case for each value of the non- dimensionalized Hamaker constant, A^∗=[10⁻²,10⁻³,10⁻⁴]. Then, startingwith

λ

^∗=100,theviscosityratioisloweredtoonetenth ofitspreviousvalue, untiltheresultsvisuallyconvergetothere- sultsof

λ

^∗=0,i.e.,thelowviscosityratiolimit.ForA^∗=10⁻²and A^∗=10⁻³,thisvalueisfound as

λ

^∗=0.1,whereas forA^∗=10⁻⁴,

λ

^∗=1seemssufficientlyclose.Foranyphyiscalvalueof

λ

^∗smaller

thaneachcorrespondinglimit,theresultsfromthe

λ

^∗=0casecan beusedtoestimatethedrainage/coalescencebehavior.

TheleftcolumnofFig.7revealsthebehavioroftcasafunction of

λ

^∗Vapp,andtheactualVappdependenceoftcispresentedonthe rightcolumnbydividingeach datasetby thecorrespondingvalue of

λ

^{∗.Ontherightcolumn, the}

λ

^∗=0resultsarenotgiven, since it implies a division by zero. However, the curves given on the left columnfor

λ

^∗=0can be adjustedby dividing

λ

^{∗Vapp by the}

smallphysicalvalueof

λ

^∗.Theonlyrestrictionisthatthisphysical

Fig. 7. Coalescence time for mobile interfaces as a function of λ^∗V appon the left, (a) - (c), and as a function of V appon right, (d) - (f). Each curve corresponds to a different λ^∗value: solid line ( λ^{∗= 0 ), stars (} λ^{∗= 0 . 1} ), dashed line ( λ^{∗= 1 ), dotted line (} λ^{∗= 10} ), dashed-dotted line ( λ^{∗= 100 ). All} results are obtained with h 00= 2 and r _∞= 30 .

valueshouldbesmallerthan0.1forA^∗=10⁻² andA^∗=10⁻³,and smallerthan 1 forA^∗=10⁻⁴ based onthe insight obtainedfrom theviscosity ratio sweeps mentioned earlierin this section. This impliesthat thechangesin

λ

^{∗ inthatlimit donotactually affect}

thebehavioroftcagainstVapp,butratheronlyshiftsthevaluesof Vapp.

Similar tothemobile interfaces,threedistinct typesofbehav- iorare alsoobservedfortheimmobileinterfaces:the linearslow drainageregionforlowvelocities,thedimpleddrainageregionaf- terV_dimp isachieved, andafter passinga minimumthe multiple- rimregionwheretcbeginstoincreasewithVapp.Thelinearregime has also been observed experimentally by Orvalho et al. (2015). Furthermore,the regime impliesthat asVapp approacheszero,tc

attains verylarge values.This observationcoincides withthe ex-

perimentalresultsofYaminskyetal.(2010),wheretheyobservea stablefilminfinitetimeintheverysmallVapplimit.Boththefirst andsecondregimes,aswellastheminimumpoint,havebeenob- servedexperimentallybyDelCastilloetal.(2011).However,tothe bestof ourknowledge,an experimental work thatcould be used tovalidatethemultiple-rimregiondoesnotexistintheliterature.

TherightcolumnofFig.7showsthatatagivenvalueofVapp,tcis largerforsmaller

λ

^{∗aslongasVappfallsintothelinearorthedim-}

pleddrainageregionsforbothvaluesof

λ

^{∗compared,orinother}

wordsaslongasmultiple-rimsarenotpresentinneitherofthe

λ

^∗

values.Ontheotherhand,ascanbeconcludedfromthepresence ofintersecting curvesinthemultiple-rimregime (e.g.

λ

^∗=1and

λ

^∗=10curvesinFig.7(f)),tccanbelargerforhigher

λ

^∗.Asanin-

creasein

λ

^{∗canbeanalyzedasadecreaseinthecontinuousphase}