Numerical simulation of continuum scale electrochemical hydrogen bubble evolution

ContentslistsavailableatScienceDirect

Applied Mathematical Modelling

journalhomepage:www.elsevier.com/locate/apm

Numerical simulation of continuum scale electrochemical hydrogen bubble evolution

Kurian J. Vachaparambil, Kristian Etienne Einarsrud

^∗

Department of Materials Science and Engineering, Norwegian University of Science and Technology (NTNU), Trondheim 7491, Norway

a rt i c l e i nf o

Article history:

Received 18 August 2020 Revised 18 April 2021 Accepted 13 May 2021 Available online 27 May 2021 Keywords:

Volume of Fluid Multiphysics simulations Electrochemical bubble evolution

a b s t ra c t

Oneoftheimportantaspectsinimprovingtheefficiencyofelectrochemicalprocesses,such aswaterelectrolysis,istheefficientremovalofbubbleswhichevolvefromtheelectrodes.

NumericalmodellingbasedonComputationalFluidDynamics(CFD)candescribethepro- cess, provideinsights intoits complexity, elucidate the underlyingmechanisms ofhow bubblesevolveandtheireffectaswellas aidindevelopingstrategiestoreducetheim- pactofthebubble.

Inthispaper, aVolumeofFluid(VOF) basedsimulationframework tostudythe evolu- tion ofhydrogenbubbles inthe orderof fewhundred micrometers,refered toas con- tinuumscalebubbles,isproposed.Theframeworkaccountsforthemultiphasenatureof the process,electrochemical reactions,dissolvedgas transport,chargetransport,interfa- cialmasstransferandassociatedbubblegrowth.Theproposedsolverisverified,fortwo- dimensionalcases,bycomparisontoanalyticalsolutionofbubblegrowthinsupersaturated solutions,stationarybubble,risingbubblesandqualitativeanalysisbasedonexperimental observationsofthevariationsincurrentbasedonstaticsimulations.Theproposedsolver isused tosimulatetheevolutionofasinglebubbleundervariouswetting conditionsof theelectrodeaswellasthecoalescencedrivenevolutionoftwobubbles.Theresultsshow thatasthebubblesdetach,itssurfaceoscillatesandtheshapeoftherisingbubbleisdeter- minedbythebalancebetweendragforceandsurfacetension.Thesesurfaceoscillations, whichcausesthe bubbletogetflattenedand elongated,resultsintemporalvariationof theelectricalcurrent.Thereductionofcurrentduetobubblegrowthisvisibleonlywhen thesesurfaceoscillationshavereduced.Thesimulationsalsoshowthecurrentasafunc- tionofthepositionofthebubbleintheinterelectrodegap.Theframeworkalsopredicts the increaseincurrentas aresult ofbubblesleavingthe surfacewhichislargerwhen theprocessiscoalescencedriven.Thesimulationsindicatethatbubblecoalescenceisthe underlyingmechanismforcontinuumscalebubbledetachment.

© 2021TheAuthors.PublishedbyElsevierInc.

ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/)

∗Corresponding author.:.

E-mail address: [email protected] (K.E. Einarsrud).

https://doi.org/10.1016/j.apm.2021.05.007

0307-904X/© 2021 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

1. Introduction

Oneofthepotentialwaystoaddresstheintermittenciesinenergyproductionviarenewablesourcesistoconvertsurplus energy into hydrogen using waterelectrolysis. This hydrogen can be used as an energy vector which would reduce the dependenceonfossilfuels,reducecarbonfootprintandfoster themove totheenvironmentally-benignhydrogeneconomy [1].Inordertoenablethistransition,thecostofhydrogenproductionfromwaterelectrolysismustbereducedusingcheaper electricityfromrenewablesourcesonenergy(aselectricitycanaccountupto70%ofthecosts[1])andmoreefficientwater electrolyzers. Ongoing research intomakingefficient waterelectrolyzershasbeen focusedprimarily on thedevelopment ofactive anddurableelectrocatalystsforthewatersplittingreactions,see[2–4].Anotheraspect toimprovetheefficiency oftheelectrolyzeristoremove bubbles,whichreducethearea oftheelectrodeincontactwithelectrolyte aswell asthe effectiveelectrolyte conductivity,see[5,6].Additionally,convectionassociatedwithbubbleevolution hasbeenreportedto increasemasstransferratesintheelectrochemicalsystems[6].Asefficientremovalofbubblesduringwaterelectrolysiscan resultinsaving10–25%oftheenergysupplied[5],understanding thedynamicsofbubbleevolutioncanaidindeveloping strategiestoimprovetheefficiencyofthewaterelectrolyzers.

Althoughexperimentalworkshaveprovidedsubstantialinsightsintothephysicsunderlyingtheelectrochemicalgasevo- lution, [7–10], numericalsimulations can provide a fundamental understanding of the coupled nature ofthe process as well asthe temporalandspatialvariationsofthe flowparameters asbubbles evolve.As electrochemical gasevolution is a multiphysics-multiscale process, the choiceof numericalsimulations employed dependson the phenomenaofinterest, see[11].Broadlyspeaking,atomisticprocessliketheelectrochemicalreactionsandthebubblenucleation[12]aretypically studiedusingmoleculardynamicsimulations[13]whereascontinuumscaleprocess(whicharefewhundredmicrometeror larger)likebubblegrowthanddetachmentcanbestudiedusingComputationalFluidDynamics(CFD)[14]andtheinterme- diate scales relevantforionmigrationandcontinuousbubbleevolution canbe investigatedusingmeso-scalemodelslike LatticeBoltzmannmethod[15,16].Thispaperdelvesintotheevolutionofcontinuumscalehydrogenbubbles,whicharefew hundredmicrometersinradius,observedinexperimentalworkslike[7–9],whicharestudiedusingmultiphase-multiphysics CFDapproaches.

In order to study the dynamics of continuum scale electrochemical bubble evolution it is necessary to simulatethe relevant multiphysics, summarized in Fig. 1. The numerical modelling ofmultiphase flows can be divided into dispersed phase andinterface resolved modelling [17]. Dispersed phase modelling (which includes Euler-Euler, Mixture andEuler- Lagrange approaches) requiresclosure models to describe themomentum transfer between thephases asthe individual bubbles arenotresolved.Thesemomentum closuretermsreliesontheassumptionofabubblesizewhichistypically set based on experimentally observed detachment diameter, see [18]. The dispersed phase modelling approaches, which are usedinmajorityofthenumericalsimulationoftheelectrochemicalgasevolution,seethereviewby [18,19]orworkslike [19,20],are typically employed to studythedynamics ofindustrial scale electrolyzers wherethe largerflow features are ofinterest. Ontheotherhandinterface resolved modellingapproaches, liketheVolume ofFluid(VOF), aretypicallyused tostudythedetailedbehaviour oftheinterfacewithouthavingtouseapproximationsinordertocapturethemomentum transfer betweenthe phaseslikeindispersed phasemodelling.Asa resultofresolving thebubble,VOFprovidesan ideal frameworktostudythedetailsofevolutionofbubbleswhichareofinterestinthiswork.TheVOFmodelusesascalar,the

Dissolved gas

Flow fields

Multiphase flow Charge transport

T ransp ort Electro

chemical

reactions In terfacial

mass transfer

Properties Prop

erties Transp

ort

Fig. 1. Illustration of the coupled nature of electrochemical gas evolution.

volumefractionofliquid(

α

1),toidentifythebubble,liquidandinterfacewhichcorrespondsto

α

1equaltozero,unityand (0,1)respectively[17].Theinterfaceiscapturedbyadvecting

α

1andthesharpnessoftheinterfaceisensuredbyusingeither the computationallycheaperalgebraicVOF,see[21],orthe sharperbutcomputationallymore demandinggeometric VOF, see[22].Boththesemethodsgeneratenumericalartifactsknownas’spuriousvelocities’duetoinaccuraciesincalculating the interfacialcurvature andinconsistent descritization, seediscussion in [21].It should be noted that thesub-cell level reconstructionofinterfaceinthegeometricVOFmethodgeneratessmallerspuriousvelocitieswhencomparedtoalgebraic VOFmethods,butthisgeometricreconstructionofinterfaceismorecomputationallydemandingandhassome limitations initsapplication,see[21–23].Thesespuriousvelocities,whicharedominantwhensimulatingsurfacetensiondrivenflows, are very well studiedin literature[21,22,24–26].Readers interestedinthe sourceof thesespurious velocitiesandreview ofon-goingresearch toaddressthemarereferredtoworkslike[24].Duetothecomputationalcostassociatedwithusing geometricVOFandchallengesrelatedtoitsimplementation,algebraicVOFmethodsareverycommonlyusedtostudytwo phaseflowscenariosincludingsurfacetensiondominantflows[21,25–31].

Anotherrelevantfeatureinelectrochemicalgasevolutionofcontinuumscalebubbleisthetransportofthedissolvedgas which requirestreatmentofinterfacialjumpconditionsforconcentration acrosstheinterface andinterfacialtransmission condition(whichensuresthattheinterfacedoesnotaccumulatethespecies),see[28,32].Twowaystosimulatethespecies transportaresingleandtwofieldapproaches.Thetwofieldapproach,describedby[32,33],usesindividualtransportequa- tions forspecies ineach phase and the interfacialjump conditionsare treated likeboundary conditions foreach phase.

The single field approach, describedby [28,29], uses a unifiedgoverning equation that accountsforthe interfacialjump conditionsto describethetransportofspeciesinbothphases. Althoughdissolvedgas,whichispresentonlyintheliquid phase,transportcanbesimulatedbyboththeseapproachesbylimitingmasstransferacrosstheinterface(likedescribedin [27]),twoandsinglefieldapproachesarerecommendedtobeusedwithgeometricandalgebraicVOFmethodsrespectively, see[28].Oncethedissolvedgasdistributionisknown,theinterfacialmass transferandtheassociatedbubblegrowthcanbe computed basedonFick’s 1stlaw(see [27,34,35]) orwithflow scenariobasedSherwoodnumber basedcorrelations(see [36]).AsSherwoodnumbercorrelationsareapplicableforspecificflowscenarios[34],bubbleevolution,whichisassociated withcomplexflowpatternsaswellasinterfacedeformation,cannotbeaccuratelydescribedbyasingleSherwoodnumber correlationtothebestoftheauthors’knowledge.

Anotheraspectofsimulatingelectrochemicalgasevolutionisthechargetransport,whichdrivestheelectrochemicalre- actions asitisproportional tothe currentdensity,basedonFaraday’slawofelectrolysis.Due totheinsulatingnature of bubbles,theirpresenceresultsinaredistributionofcurrent,whichincreasestheelectricalresistance[37,38].Thecommonly usedapproachtoensurechargeconservationisbasedonGauss’s law,oneofthefourMaxwell’sequationsforelectromag- netism,andcurrentdensityexpressedbasedontheOhm’slawwhichhasbeenusedinworkslike[39,40].Abenefitofusing theinterfaceresolvedapproach,likeVOF,isthattheeffectiveconductivityoftheheterogenousmedia,i.e.bubble-electrolyte, can becomputedwithouthavingtousesemi-empiricalrelations, liketheBruggemancorrelation[6],whichisrequiredin dispersedphasemodellingapproaches[19,20]asshownby[41].

Due tothecoupled multiphysicsrelevantin electrochemicalgasevolution, literaturewhich employsinterface resolved simulationstostudythedynamicsofelectrochemicallygeneratedbubblesarequitelimitedtothebestoftheknowledgeof theauthors.TheseworkscanbedividedintopureandhybridVOFmethods.In’pure’VOFmethods,like[36],theevolution ofasinglehydrogenbubblewasstudiedbysimulatingthegrowthdriven byinterfacialmasstransfer(basedonSherwood number)withoutaccountingforthechargetransport.Inourpreviouswork[41],wedescribedthevariousmodulesneeded to simulatecontinuum scale electrochemicalgasevolution, seeFig. 1,ina decoupled mannerinaddition tobriefly illus- trating thepotentialoffullycoupledVOFsolverviasimulationofasingleslidingbubbleona verticalelectrode.Although notbasedonVOFapproach,therecentworkby[42]usedanotherinterfaceresolvedapproach,basedontheFrontTracking method,tostudytheimpactofbubbleevolutioninducedconvectiononionictransportinwaterelectrolysis.In’hybrid’VOF approaches,like[39,40,43],VOFwascoupledwithsub-gribbubblestreatedviadispersedmodellingapproachesandthere- solvedbubblewereassumedtogrowonlyviacoalescence.ThesehybridVOF,typicallyusedtosimulatethecarbondioxide evolutionattheanodeduringaluminiumproduction,isusedtosimulatebubbleevolutiononanwholeelectrode[39,40].re- portedthetransientevolutionofvoltage(underconstantcurrentcondition)asresultofbubbleevolutionfromtheelectrode whereas[43] didnotsimulatethecurrentdistributionnorthetransportofdissolvedgasgeneratedfromtheelectrochemi- calreactions.Apartfromthesestudies,othersimplifiedstudieshaveemployedinterfaceresolvedsimulationstoinvestigate the behaviour ofthe interface without considering any multiphysics effects like [44].Although the previous works have providedsubstantialknowledgeintothemodellingofevolutionofelectrochemicallygeneratedbubbles,thereisstillalack of computationalmodels that treatthecomplexityassociated withthemultiphysics andmultiscalenature oftheprocess ashighlightedintherecentreviewby[11].Additionally,bubbleevolutionduringwaterelectrolysisduetocoalescence,ob- served inexperimentalstudieslike[7,8,10],hasnotbeeninvestigatedextensivelyusingnumericalsimulationsasdiscussed inthereviewby[45].

Inthispaper,whichisbasedonourrecentworkwhichwaspresentedatCFD2020[41],weaddressthislackofknowl- edge bydevelopinga coupledmultiphysicssolverthat canhandlethecontinuum scalehydrogen bubbleevolutionduring waterelectrolysis.TheproposedsolverisbasedonthealgebraicVOFframeworkavailableinOpenFOAM® 6whichismod- ifiedtoaccountfortransportofdissolvedgas[29] alongwiththeassociatedsupersaturationdriven bubblegrowth(based on[27])andchargetransport(basedonGauss’slawandOhm’slaw,see[39]) whicharecoupledbasedonelectrochemical

reactiondefinedbasedonFaraday’slawofelectrolysis.Inordertoreducethecomplexity,theproposedsolverisdeveloped basedonthefollowingsimplifications/assumptions:

• The bubbleevolutionandtheelectrochemical reactions occursonly atthecathodeandtheanode/counter-electrodeis assumedtonotaffecttheprocess,

• Theproposed modeltreats theevolutionofonlycontinuum scalebubbles,whichareintheorderoffew hundredmi- crometers,

• The liquid and the continuum scale bubbles treated in the simulations are assumedto have a constant density and viscosity,

• Theflowisassumedtobeisothermalandlaminar,

• Theinterfaceisassumedtobealwayssaturatedandsaturationconcentrationisassumedtobeconstantwiththevaria- tionofhydrostaticpressureandLaplacepressureinsub-millimetersizedbubblesasitevolves,

• The proposed solver assumes a constant value of surfacetension in the simulations. This assumption results in not treatinggradientsinsurfacetensionwhichhasbeenobservedtogenerateMarangoniconvectionduringelectrochemical hydrogenevolution[9,46],

• TheinterfacialmasstransferistreatedusingtheFick’s1stlawasthemasstransfercoefficientsareunknownapriori,

• The spatial andtemporal gradient of ions in the computational domain is assumedto be negligible due to the high concentrationionspresenttypicallyintheelectrolyteusedforwaterelectrolysis[47],

• Thesystemisassumedtobeunderconstantpotentialdifference,sothebubbleevolutionleadstochangesincurrent.The proposedsolveraccountsforthecurrentvariationasaresultoftheohmiccontributionofthebubblespresentinthebulk andonthesurface. Thechange incurrentdueto thecontributionsfromthe surfaceandconcentration overpotentials, see[48],asaresultofbubbleevolutionisneglectedinthispaper.

Theabilityoftheproposedsolverisshowcasedvia2Dsimulationsoftheevolutionofthesinglehydrogenbubblefrom cathode forvarious wetting conditionatthe electrode. In addition to singlebubblesimulations, the evolution of bubble drivenbycoalescenceisalsoinvestigated.Thetemporalchangeinthecurrentasbubbleevolvesisanalyzedforthesecases.

The noveltyoftheproposed solver, whencomparedto previousworksthat employ multiphysicsbasedinterface resolved simulationslike[36,42,44],isitsabilitytosimulatecoupledphenomenarelatedtodissolvedhydrogentransport,interfacial masstransferduetosupersaturationandcurrentvariationasaresultofbubbleevolutiondrivenbycoalescenceandbubble growth.Inthe interestofknowledge disseminationandtoencouragefurtherscientific advances,theproposed solverwill bereleasedatGitHub[49].

2. Proposedsolver 2.1. Governingequations

In the Volume of Fluid(VOF) method, see [21,50], the interface and the two phases are representedusing a volume averagedscalarcalledvolumefractionofliquid(

α

1),definedas

α

¹

(

^x,t

)

=

_{0(withinPhase2orgas)}

0<

α

^{1<1(attheinterface)}

1(withinPhase1orliquid).

(1)

Basedonvalueof

α

1,thevolumefractionofgas(

α

2)iscalculatedas1−

α

1.Asinglefielddescriptionofthegoverningequa- tions forthesimulation ofmultiphaseflows canbe developedby conditional volume-averagingtechniquewhichinvolves conditioningtheconservationequationvalidinliquidandgasregions forphasediscrimination anditssubsequentvolume averaging,see[28,51]fordetails.Thisconditionalvolume-averagingtechniqueisusedtoobtainthesinglefieldequationsfor continuity,momentumandvolumefractiontransportusedbyVOFapproach,thisderivationisdiscussedin[28,52,53].For thesake ofsimplicity,onlythefinalformsoftheseequationsusedintheVOFapproach,afterthenotationsofconditional volume-averaging aredropped, aredescribedinthiswork whichissimilarto representationusedcommonlyinliterature [21,25–27,29,50,51].

Theinterfacedynamicsarecapturedusingatransportequationofliquidvolumefraction:

∂α

¹

∂

^{t +}

∇

^·

(

^U1

α

¹

)

=S˙_α, (2)

whereU₁ isthe velocity oftheliquidphase andS˙_α isa sourcetermthat accountsforbubblegrowth[50].Asinglefield descriptionofthevelocity(U)canbedescribedas

U =

α

1U1+

α

2U2, (3)

whereU₁ andU₂ are the velocity fieldsin liquidandgas phase respectively[50]. Therelative velocity between thetwo phases(Ur)iscomputedas

Ur=U1−U2. (4)

NowU₁isdefinedintermsofthesinglefieldformulationofvelocity(U)andrelativevelocity(Ur)as

α

1U₁=

α

1U+

α

1(¹−

α

1)^Ur,from

α

1×Eq.3and

α

1(¹−

α

1)×Eq.4(see[50]),whichcanbeusedtorewriteEq.2as

∂α

¹

∂

^{t +}

∇

^·

( α

^1U

)

+

∇

^·

( α

¹

(

¹−

α

¹

)

^Ur

)

=S˙_α. (5) The third termin Eq.5, which is non-zero only atthe interface due to the product

α

1(¹−

α

1)^{, is used to sharpen the} interface.Eq.5belongstoAlgebraicVOFmethod,see[21],andthecalculationofUrisestimatedbasedonEq.21[25–27,54]. ThesourcetermusedinEq.5totreatthebubblegrowthiscomputedbasedonworkslike[27,34]as

S˙_α=

α

1

∇

^{·U. (6)}

The fluid properties (

χ

^{) like density(}

ρ

^{),viscosity (}

ν

^{) and electrical} conductivity(^{) are calculated based onvolume} fractionweightedaverageas

χ

=

α

1

χ

1+

α

2

χ

2. (7)

Thecontinuityequationiswrittenas

∇

·U =m˙

ρ

^{, (8)}

wheretermontherightsideontheaboveequationaccountsforthebubblegrowthduetosupersaturationandm˙ iscom- putedbasedonEq.17.

Themomentumequationis

∂ρ

^U

∂

^{t +}

∇

^·

( ρ

^UU

)

⁼⁻

∇

^p+

∇

^·

μ ( ∇

^U+

∇

^UT

)

+

ρ

^{g+FST, (9)}

where F_ST accountsforsurfacetension,gis theacceleration duetogravity and

∇

·

μ

(

∇

^U+

∇

^UT)

is theviscous term.

Definingp_rghasequaltop−

ρ

^g·x,−

∇

^p+

ρ

^{gcanbewrittenas}−

∇

^prgh−g·x

∇ ρ

^{[21].Theviscoustermsinthemomentum}

equation,

∇

·

μ

(

∇

^U+

∇

^UT)

canbewrittenas

∇

·(

μ∇

^U)+

∇

^U·

∇μ

^{,see[21,28].Thesurfacetensiontermiscomputed}

usingtheSharpSurfaceForce(SSF)model,describedin[25,26],as

FST=

σ κ

f inal

∇α

sh, (10)

where

κ

f inal isthe interface curvature (whose calculation is describedin Appendix A) and

α

sh is the sharpenedvolume fractionofphase1whichisdeterminedas

α

sh= 1 1−C_sh

min

max

α

^1,C₂^sh

,1−C_sh 2

−C_sh 2

, (11)

whereC_sh isauserdefinedsharpeningcoefficientC_sh∈[0,1)^{,see[26] fordiscussion onitsinfluence.Additional noteson} thesurfacetensionmodellingisavailableinAppendixA.

Nowsubstitutingthedefinitionsofp_rgh,viscoustermsandsurfacetensiontermusingSSFmodelinEq.9toget

∂ρ

^U

∂

^{t +}

∇

^·

( ρ

^UU

)

=−

∇

^prgh+

∇

^·

( μ∇

^U

)

+

∇

^U·

∇μ

^−g·x

∇ρ

⁺

σκ

f inal

∇α

sh. (12)

Thetransportofspeciesneedstoaccountforconcentrationjumpacrosstheinterfaceinadditiontointerfacialtransmis- sionconditionapartfromspeciestransportineachphase,see[28].Asinglefieldformulationofthespeciestransportcanbe obtainedby conditionalvolume-averagingofspeciestransportequation ineachphase alongwithincorporatingthe afore- mentionedinterface jumpconditions,see[28]forderivation.Avariantofthesinglefield formulationforspeciestransport istheCompressiveContinuousSpeciesTransfer(CCST)model,whichwasproposedby[29],:

∂

^Ci

∂

^{t +}

∇

^·

(

^UCi

)

=

∇

^·

Dˆi

∇

^Ci−DˆiBCi

∇ α

^1−B

α

¹

α

^2UrCi

+Si, (13)

whereS_iisasinkterm(computedinEq.18),Bisequalto(¹−He_i)/(

α

1+

α

2He_i)^,Heidescribesconcentrationjumpcondi- tion(ratioofthespeciesconcentrationingasphasetoliquidphaseattheinterface),U_risthecompressivevelocitydefined in Eq. 21, and Dˆ_i is the harmonic average of the diffusion coefficients of the phases. It should be pointed out that the conditional volume-averaged notationsofthe derived singlefieldformulationof speciestransportequation is dropped in Eq. 13for simplicity,similar to [29]. Using thesinglefield formulation ofconcentration in CCST, dissolved hydrogencan be simulated bylimiting theinterfacialmass transfer,that is inherentlytreatedusing theinterfacialjump conditions,by settingHe_itoasmallvalue[27,34,41].AnadditionalconsequenceofusingasmallvalueofHe_iisthattheC_iattheinterface becomesclosetozerowhichcanbeinterpretedastheconcentrationfieldobtainedfromCCSTwhichisoverthesaturation condition,see[27,34].

As a result of simulating the dissolved gas, the treatment of interfacial mass transfer andassociated bubble growth requiredeterminationofthesourceterminEq.8andsinktermatinterfacefordissolvedgasinEq.13.Inthispaper,these

source andsinkterms arecalculated basedon theapproach initiallydeveloped forevaporation[55,56]andsubsequently extendedforsupersaturationdrivenbubblegrowth[27,34].Thedrivingforceforinterfacialmasstransferiscomputedbased onFick’s1stlawas

j=MiDi,1

|∇

^Ci

|

, (14)

whereM_iisthemolarmassofhydrogen(i)andD_i,1isthediffusioncoefficientofthehydrogenintheliquid[27,34].Itshould be notedthatEq.14isbasedontheassumptionthattheconcentration gradientinthetangentialdirectiontointerface is negligible andthisapproximationisconsistent withthederivationofthesinglefieldformulationofspeciestransport,i.e.

Eq.13,see[28].Theaboveformulationofdrivingforceforinterfacialmasstransferimplicitlyassumesthatthespeciesare dilute,whilethisisstrictlynotthecaseinageneralsetting,weusethisformulationasabasisofourframework.Basedon Eq.14,thelocalmasstransferrateiscomputedattheliquidsideoftheinterfaceas

ψ

0=N j

α

1

|∇α

1

|

, (15)

where N is normalization factor computed as

|∇α

1

|

^dV/

α

1

|∇ α

1

|

^{dV, see [27,55,56].}

ψ

0 is smeared, asproposed by [27,55,56],using:

D

^t

∇

²

ψ

=

ψ

−

ψ

0, (16)

where D^{t isthe userdefinedvaluewhich controlstheextentofsmearing. Using}

ψ

^{,thesource termforthe continuity}

equationwhichaccountsforthebubblegrowth,isredistributedintheregionwhere

α

1<0.001as

˙

m=A

α

2

ψ

, (17)

where A is a normalizationfactor computedas

ψ

0dV/

α

2

ψ

^{dV, see [27,55,56]. Thisprocedure to define m˙, which is}

described further in thework by [55,56], enablesthe interface to get advected just by the velocity field, without being influenced by the sourceterm. Thesink termforthe dissolved hydrogento account forthe interfacialmass transport is computedattheliquidsideoftheinterfaceas

S_i=−N

α

¹

(

^j

|∇α

¹

| )

M_i , (18)

whereNisthenormalizationfactorusedinEq.15,see[27,34].Generationofdissolvedhydrogencausingsupersaturationin thevicinity oftheelectrode duetoelectrochemicalreactions, whichare proportionaltothelocalvalueofcurrentdensity, areaddedasboundaryconditions(furtherdescribedinSection.3.1).ThechargeconservationissolvedbasedonGauss’slaw as

∇

^{·i=0, (19)}

whereiisthelocalcurrentdensitywhichiscalculatedbasedonOhm’slawas

i=−

∇

⁽²⁰⁾

where ^{isthevolumefractionweightedaverageelectrical}conductivity,basedonEq.7,and^{istheelectricalpotential.}

It should be noted that Eq.20 is obtained from Nernst-Planck equation under the assumption of negligible spatial and temporalgradientsofionsinthecomputationaldomaini.e.electroneutrality,see[57].Thisapproximationcanbephysically motivatedby highconcentrationsofionspresentinthe electrolytewhichistypicallyused forwaterelectrolysis,see[47]. Eq. 19 andEq. 20has been previously used to simulatecharge transport inelectrochemical gas evolution in works like [39–41].

2.2. Solutionalgorithmandnumericaldetails

The proposedsolverisdevelopedbasedontheVOFsolverinOpenFOAM® 6,interFoam,whichsupportsalgebraicVOF basedsimulationoftwoimcompressibleandisothermalflowwithoutanyphasechange.Theproposedsolveraddsmodules, illustrated inFig. 1,to interFoam in orderto enabletreatment of the multiphysicscoupling relevant incontinuum scale electrochemicalhydrogenevolution.Theoverallsolutionalgorithmcanbesummarizedas:

1. Compute

α

1bysolvingEq.5usingthesemi-implicitMultidimensionalUniversal LimiterwithExplicitSolution(MULES) methodwhichusesaimplicitpredictorandexplicitcorrectorstepstoensuretheboundednessof

α

1betweenzeroand unity [21,58].Once

α

1 iscomputed, the volume fractionof phase 2 orbubble (

α

2) iscomputed as

α

2=1−

α

1. In a singlefieldformulationofvelocity,Ur whichisusedinthesurfacecompressionterminthevolumefractionadvection equation,Eq.5,isunknownandestimated,basedon[25–27,54],as

Ur=C_α

φ

|

^Sf

|

^{n, (21)}

wherenistheunitnormaltotheinterface,S_f istheareavectorofthecellface,

φ

^{isthevolumeflux,andC}α istheuser definedcompressionfactorwhichisusuallysetbetweenvaluesofzeroandfour[59,21,25].

Table 1

Discretisation schemes.

Modelling term Keyword Scheme

Time derivatives ddtSchemes Euler [61]

Divergence term

∇^·(ρ^{U U} ) vanLeerV [61]

∇^·(U α1), ∇^·(U C i) vanLeer [61]

∇·(U rα1(1 −α1)) interfaceCompression [21,61]

∇·(D ^ˆ iBC i∇α1) vanLeer [29,61]

∇·(Bα1α2U rC i) vanLeer [29,61]

Gradient term gradSchemes linear [61]

Laplacian term laplacianSchemes linear corrected [61]

Other snGradSchemes corrected [61]

interpolationSchemes linear [61]

Table 2

Solvers used for the discretised governing equations, see [61] for further details.

Flow

variable Linear solver

Smoother/

preconditioner Tolerance

p rgh PCG GAMG 10 ⁻¹⁰

U smoothSolver symGaussSeidel 10 ⁻¹⁰

α1 smoothSolver symGaussSeidel 10 ⁻¹⁰

GAMG DICGaussSeidel 10 ⁻¹⁰

C i PBiCGStab diagonal 10 ⁻¹⁰

ψ GAMG DICGaussSeidel 10 ⁻¹⁰

2. Oncetheinterfaceisknownatthecurrenttimestep,fluidproperties(i.e.densityandviscosity)iscomputedbasedon Eq.7.

3. The drivingforce forinterfacialmass transferisthen computedbased onthe Fick’s1stlaw,i.e.Eq.14, whichisthen usedtocomputethelocalmasstransferrate(

ψ

0)usingEq.15.

4. Inordertoimprovenumericalstability,

ψ

0 issmearedtoobtain

ψ

^{basedonEq.16.}

5. Nowsolvecontinuityandmomentumequations,i.e.Eq.8andEq.12,aresolvedusingPressureImplicitwithSplittingof Operators(PISO)algorithm[60].Thiscanbebriefly summarizedasiterativeprocedureinwhichpredictedvelocityfield iscorrectedbasedonpressurecorrectionequation,see[21]fordetailsoftheimplementedalgorithm.

6. ComputetheelectricalpotentialandcurrentdensityusingEq.19andEq.20. 7. SolveforthedissolvedhydrogentransportbasedonEq.13.

8. Advancetothenext timestep.Adaptivetimestepping isusedinthesimulationstoallowforcomputingthetime step suchthatmaximumCourantnumber,whichisdefinedbytheuser,isalwayssatisfied.Ifamaximumtimestepisdefined, thesolvertakesthe minimumofthe userdefinedmaximumtime stepandthetime stepassociatedwiththeCourant number[61].

All thevariables arestoredatthe cellcenter likeina collocatedarrangement butcomputationsare performedat the cell faces(which wasinitiallyproposed by[62])toprevent checkerboardtypeerrorsaswell asinconsistentdescritization (discussedin[21,50]).TheFinite VolumeMethodisemployedto solvethegoverningequationswithspatialandtemporal termsdiscretizedbasedonTable1.Thediscretizedgoverningequationsaresolvedwithiterativesolverswiththerelevant preconditioners/smoothers tabulated in Table2.The maximum numberofiterationsis setsuch that thesolution offlow variableconvergesandsatisfiesthetolerancecriteriondefinedinTable2.Asummary oftheparametersusedinthesolu- tionalgorithm isprovidedinTable3.Unlessotherwisenoted, thetimesteptakenbytheproposed solverisbasedonthe maximumCourantnumberwhichissetequalto0.05.

3. Definitionoftestcases

3.1. Computationaldomainandboundaryconditions

Thetwo-dimensionalcomputationaldomainusedinthesimulationsisarectangle,oflength(Le)andheight(Lac)equalto 10mmand5mmrespectively,asillustratedinFig.2.Althoughthesimulationsdiscussedinthepaperaretwo-dimensional, OpenFOAM® requires afinite thicknesswitha singlecellresolution inthethird direction,representedby h,whichis set to1micrometer.These’additional’boundariesaresetto’empty’boundaryconditionstoperforma2Dsimulation.Thearea of each electrode, A, is computedas L_e×h and the volume of the entire domain is equal to A×L_ac. The left andright boundariesareassignedzerogradientforvolumefractionofliquid(

α

1),dissolvedgasconcentration(C_i),localmasstransfer rate(

ψ

^{) andelectrical potential() whereassinglefield velocity(U) andmodifiedpressure(p}rgh)are setasinletOutlet (zerogradientforoutflowbutbackflowintothedomainisrestrictedbysettingzerovelocity)andfixedvalue(equalto0Pa)

Table 3

Summary of the other settings used in the solution algorithm.

Parameter Value Notes

nAlphaCorr 2 Number of α1correction [27] .

nAlphaSubCycles 1 Number of sub-cycles of α1equation within each time step [61] . cAlpha ( C α) 1 Used for interface compression in Eq. 21 [21,61] .

MULESCorr yes Switches on semi-implicit MULES [58] .

nLimiterIter 3 Number of MULES iterations over the limiter [58] . momentumPredictor no Controls solving of the momentum predictor [61] .

nOuterCorrectors 1 PISO algorithm is selected by setting this parameter equal to unity in PIMPLE algorithm [61] .

nCorrectors 3 The number of times the PISO algorithm solves the pressure and momentum equation in each time step [61] .

nNonOrthogonalCorrectors 0 Used when meshes are non-orthogonal [61] .

relaxationFactors 1 Specifies the under-relaxation factors used for fields and equations, see [25] . C sh 0.3 Sharpening coefficient used in Eq. 11 : set to enable simulation of a sub-millimeter

bubble, see [26,27] .

D t 10 ⁻⁶ Smearing of the local mass transfer rate (based on sensitivity studies performed in [27] ).

He i 10 ⁻⁴ Concentration jump across interface, used in Eq. 13 : set to simulate dissolved gas, see parametric study in [27] .

Initialized bubble Bottom/Electrode

Top

Left Right

δ

L

_ac

L

_e

g C

_i

based on Eq.22

α

₁

= 1

C

_i

= 0 α

₁

= 0

Fig. 2. Illustration of computational domain used in the simulations and the initial conditions used in the simulations.

respectively.Thetopboundaryisassignedzerogradientfor

α

1,fixedvalue(equalto0mol/m³)forC_i,fixedvalue(of0.2V) forelectrical potential (^), fixedFluxPressure for p_rgh,zerogradientfor

ψ

^{andno-slipforvelocity.The bottomboundary,}

referred toaselectrode,isassignedzerogradientfor

α

1alongwiththestaticcontactangle(

θ

^{) conditionattheelectrode,}

fixed gradient (which is computedbased on Faraday’slawof electrolysis¹) forthe dissolved gasconcentration (C_i), fixed value (of0.0V)for^,fixedFluxPressure formodifiedpressure,zerogradientfor

ψ

^{andno-slipforvelocity.Theboundary}

conditionsofelectricalpotential(^{)arechosen suchthatthecurrentdensityobtainedfromthesystemwhenbubblesare} notpresentinthedomainiscalculatedfromOhm’slaw,as

κ

1(^top−electrode)/Lac,tobeequalto1210.8A/m²whichisin thepracticalrange(between1000–3000A/m²)usedincommerciallyusedalkalineelectrolyzersystems,see[63].

AsC_iindicates theconcentration ofthedissolved gaswhichexceedsthe saturationcondition,which forhydrogendis- solvedinwaterunder1atmand25^◦Cis0.79mol/m³ [64],theconcentrationboundarylayerattheelectrodeisinitialized basedontheverticaldistancefromtheelectrode(y),asshowninFig.2,suchthat

C_i= C_i,y=δ(wheny>

δ )

^,

Ci,y=0+

(

^Ci,y=δ−Ci,y=0

)

y/

δ

^(wheny≤

δ )

, (22)

where

δ

^istheconcentrationboundarylayerthicknesssetto0.5mm(seeAppendixE),C_i,y=0andC_i,y=δistheconcentration atthe electrode anddistanceof

δ

^{(or larger)which isset equalto125.66 mol/m3 (based on}supersaturation reportedin [65]) and 0 mol/m³ (which corresponds to saturation)respectively. It should be noted that the initialized concentration boundarylayerdoesnotaccountforthedepletedconcentrationofthedissolvedgasnearthebubbleduetointerfacialmass

1The boundary condition for C iat the electrode, which is computed based on Faraday’s law of electrolysis, is ∂^nC i= |ⁱ |α1/ (nF D _i,1), where n is the number of electrons transferred for the electrochemical reaction to produce hydrogen (equal to 2), F is the Faraday’s constant (equal to 96,485 As/mol) and D _i,1is the diffusion coefficient of dissolved hydrogen in the liquid.

Table 4

The fluid properties, based on 1atm and 25 ^◦C, for Phase 1 (liq- uid/electrolyte) and Phase 2 (hydrogen bubble) used in the simulations.

Properties Dimensions Phase 1 Phase 2

Density ( ρ^{) kg/m 3 1075.05 0.082}

Viscosity ( ν) m ²/s 9.89 ×10 ⁻⁷ 0.00011 Diffusion coefficient ( D i) m ²/s 4.80 ×10 ⁻⁹ 1 ×10 ⁻⁵ Electrical conductivity ( ) S/m 30.27 1 ×10 ⁻¹³ Molar mass ( M i) kg/mol 2 ×10 ⁻³

Surface tension ( σ) N/m 0.072

transferacrosstheinterface.Thesinglebubbleinthedomain,whichisinitializedasinFig.2,islocatedsuchthatitscenter is5mmfromtheleftboundaryand0mmfromthebottomboundary(yc).Inthecaseoftwobubbles,usedincoalescence studies,withradiiequaltoR_2aandR_2b,thebubblesareinitializedassemicirclessuchthatthecentersofthebubblesarea distanceof5×10⁻³−R_2aand5×10⁻³+R_2bfromtheleftboundary.Astheelectrodesarehorizontalwithrespecttogravity, whichissetas

|

^g

|

=9.81m²/s, theevolutionofthebubblewouldleadtodetachmentandsubsequentverticalrisedueto buoyancy.The fluidpropertiesusedinthe simulationsisdescribedinTable4.Thecomputational geometryismeshedby 1600×800hexahedralcellsbasedonthemeshconvergencestudyindescribedinAppendixC.

3.2. Definitionofparameters

Althoughvisual comparisonisaveryusefultoolstointerprettheresults,tocomparevarioussimulations andquantify its results,itisimportantandeasierto usestandardizedparameters. Thestandardized parameterswhich areusedinthe paperforanalyzingtheresultsaredefinedinthissection.

Normalizedbubblevolume,Thevolumeofthebubbleatanytime(V)iscalculatedas

α

2dVandthecorrespondingarea ofthebubbleisequaltoV/h,wherehisunit cellthickness.Thenormalizedbubblevolumeiscalculatedasthefractionof thecomputationaldomainoccupiedbythebubble(f)whichisdeterminedas

f= V A×Lac

(23)

whereAistheareaoftheelectrodeandL_acistheinter-electrodedistance.ForthegeometryusedinthesimulationA×L_ac isequalto5×10⁻¹¹m⁻³,seeFig.2.

Risevelocity.Themeanvelocitywithwhichabubblerisesiscomputedas Urise=

v α

2dV

α

^{2dV (24)}

where

v

istheverticalcomponentofthevelocityvector.Therisevelocity,computedbasedonEq.24,hasbeenbeenprevi- ouslyusedinworkslike[66]and[25]tostudythedynamicsofsinglerisingbubble.

Bubblecoverageofelectrode.Thefractionoftheareaoftheelectrodecoveredbythebubbleiscomputedas

=

α

2

|

^S

|

|

^S

|

^{, (25)}

whichiscomputedattheelectrodeboundarywith

|

^S

|

representingthemagnitudeofthesurfaceareaoftheindividualmesh cellattheboundaryand

|

^S

|

^{istheareaoftheelectrode(A),whichisequalto10−8m2.Thefractionoftheelectrodearea}

incontactwiththeelectrolyteisequalto1−^{.Whenthebubblehasdetachedorwhenthewholeelectrodeisincontact} withtheelectrolyte,^{reducestozero.}

Bubbledetachment.Basedonthetemporalvariationof1−^{,thebubbleisconsideredtohavedetachedwhen1}−>

0.999andthecorrespondingtime,whichisindicatedby

τ

d,isconsideredasthedetachmenttimeofthebubble.

Average current.As theproposed solver computeslocalcurrentdensity(i),thecurrentobtainedfromthe system(I) is equalto

I=

i·S, (26)

whichisthe sumofdotproduct oflocalvalue ofcurrentdensityandcellarea vectorofindividualmesh attheelectrode boundary.Forsimplicity,anormalizedaveragecurrent(I/I₀)isusedinthispaperwhichisdeterminedastheratioofItoI₀ whichisthecurrentobtainedwhenbubblesareabsentinthesystem.²

Bubbledeformation. Asthebubbledetaches fromthesurface, itundergoes deformationwhichismeasured horizontally (^{x)andvertically(y).Thehorizontal}deformationofthebubbleatanytimestepiscalculatedasthedifferencebetween

2Theoretically the current which can be obtained for the applied potential difference is when the system does not have any bubbles can be calculated by substituting relevant values in I 0= 1A (topboundary−electrode)/L acwhich givesI 0equal to 1.2108×10 ⁻⁵A.

Fig. 3. Temporal change in α1as initialized square bubble reaches the equilibrium shape i.e. a stationary circular bubble..

maximum andminimumxcoordinates oftheinterface (

α

1=0.5). Similarly^{y is thedifference betweenmaximumand} minimumycoordinatesoftheinterface.

4. Verificationoftheproposedsolver

Inthiswork,theproposedsolverisverifiedinamodularizedmannerwhichallowsforsequentialverificationofpartsof theproposed solverwithexistinganalyticalandcomputationalbenchmarks. Thebenchmarksusedinthisworkarebased on static 2D bubble, dynamic caseof rising 2D bubbles, growth of 2D bubble driven by supersaturation andqualitative assessment ofthevariation ofcurrent. Apartfromtheseverificationcases,the relevantimbalance (massof liquid/bubble

Fig. 4. Temporal change in the error in calculating the Laplace pressure, calculated as ^pc−_pc^pwhere p cis determined based on Eq. 27 , p = ^α_α²₂^pdV_dV −p 0

and p 0is the operating pressure used in the simulation (equal to 0 Pa)..

Fig. 5. Evolution of spurious velocities, calculated as max( |U |), during stationary bubble simulation..

Fig. 6. Comparison of bubble morphology at t = 3 s predicted by the solver and results reported in literature [25,66] : a) TC1 b) TC2. The data used for comparison are obtained from benchmark simulations of [66] and our previous work [25] .

Fig. 7. Comparison of the prediction of the proposed solver with the Extended Scriven model to predict the bubble growth from a pre-existing bubble in a uniformly supersaturated solution: a) Bubble radius b) Growth rate.

Fig. 8. The distribution of |i | for static bubble simulations: a) TC1a b) TC2a and c) TC2d. The white lines represents the current lines in the liquid phase whereas the black line indicates the interface (i.e. α1= 0 . 5 )..

0.00450 0.005 0.0055

0.0002 0.0004 0.0006 0.0008 0.001

x m)

y (m)

a

0.00450 0.005 0.0055

0.0002 0.0004 0.0006 0.0008 0.001

x m)

y (m)

b

0.00450 0.005 0.0055

0.0002 0.0004 0.0006 0.0008 0.001

x m)

y (m)

c

Fig. 9. The evolution of the bubble (isosurface of α1= 0.5) at 0s ( ), 5 ×10 ⁻⁴s( ), 1 ×10 ⁻³s( ), 1.5 ×10 ⁻³s( ) and 2 ×10 ⁻³s( ) for the various wetting condition at the electrode: (a) 15 ^◦in SCT0, (b) 30 ^◦in SCT1, and (c) 45 ^◦in SCT2.

and dissolved gas), mesh andtemporal convergence studies have also been performedin Appendix B, Appendix C and AppendixDrespectively.

4.1. Stationarybubble

Thefirststepintheverificationoftheproposedsolverisbasedontheabilitytosimulateevolutionofthe2Dsquarebub- ble,initializedinthecomputationaldomain,whichduetosurfacetension(andabsenceofgravity)reachestheequilibrium shape i.e.astationarycircularbubble.Atequilibrium,theLaplacepressureofthe2Dbubblecanbecalculatedanalytically usingtheYoung-Laplaceequationas

^pc=

σ

R, (27)

where

σ

^{is the surfacetension andR isthe 2D bubble radius.} Additionally,atequilibrium, thevelocities in the domain shouldtheoreticallyreducetozero.

Fig. 10. Temporal change in the fraction of area of the electrode in contact with the electrolyte ( 1 −) as the bubble evolves in SCT0 ( θ= 15 ^◦), SCT1 ( θ= 30 ^◦) and SCT2 ( θ= 45 ^◦). The timesteps at which the isosurface of α1= 0.5 is extracted to plot Fig. 9 is illustrated using vertical dashed lines.

Fig. 11. The temporal variation of various relevant parameters during bubble detachment in SCT0 ( θ= 15 ^◦): a) normalized current ( I/I 0), b) horizontal deformation ( x ), c) vertical deformation ( y ), d) normalized volume of bubble ( f). The vertical black line indicates the detachment time, τd, and the vertical red lines are equidistant grid line at every 0.001s to enable comparison between the plots.

Fig. 12. The temporal variation of various relevant parameters during bubble detachment in SCT1 ( θ= 30 ^◦): a) normalized current ( I/I 0), b) horizontal deformation ( x ), c) vertical deformation ( y ), d) normalized volume of bubble ( f). The vertical black line indicates the detachment time, τd, and the vertical red lines are equidistant grid line at every 0.001s to enable comparison between the plots.

Inorderto simulatetheabove mentionedflowscenario,gasisinitialized asa squareregion(equal to0.6×0.6mm²) inthecenterofthecomputationaldomain(ofdimensions1.8×1.8mm²).Thefluidpropertiesusedforthesimulationsis describedinTable4withtheexception ofgravitywhichissettozero.Theboundaryconditionsusedfor

α

¹,

ψ

^{andU are}

zerogradientwhere p_rgh usesa fixedvalue(equalto0Pa)onall fourboundaries.C_iusesfixedvalue(equal to0mol/m³) onthetopboundaryandzerogradientontheremainingthreeboundaries.For^{,theleftandrightboundariesareassigned} fixed value (equal to 0 V) andother boundaries are set to be zerogradient. The computational domain is meshed with 120×120cells.Theoretically,theradiusofcircle(whichhasthesameareaastheinitializedsquare)is338.6 micrometers andthecorrespondingLaplacepressure,calculatedusingEq.27,isequalto212.64Pa.

TheevolutionofthebubblefromtheinitializedsquaretotheequilibriumcircularshapeisshowninFig.3.Asthebubble reachesequilibrium, theLaplacepressure obtainedfromthe simulationalsoreachesreadyaconstant value,seeFig.4.At t=0.02 s, the Laplace pressure obtained fromthe simulations has an absolute errorof around 9% with respect to the correspondinganalyticalsolution,whichissimilartoerrorsreportedinworkslike[25,26].Althoughthebubblereachedthe equilibriumshape,theexistenceofspurious velocities,whicharise fromerrorsinestimatingtheinterfacialcurvature(see AppendixA),causesanon-zerovelocitiesinthecomputationaldomainwhichhasbeenshowninFig.5.

4.2. Risingbubble

Thenextverificationcaseusedinthispaperisbasedontheworkby[66],whoproposedcomputationalbenchmarkfor two2D risingbubbleswhichhasbeenextensivelyusedinliterature,like[25,30,31].Theworkby[66]reportedthebubble morphology,risevelocitiesandcircularityfortworisingbubblescenarios(whichdifferbasedontheCapillarynumber)using codeslikeTP2D,FreeLIFEandMoonNMD.Inthispaper,verificationoftheproposedsolverisbasedonbubblemorphology reported bythe originalwork of[66]usingFreeLIFEcode and[25] whoperformedtheserising bubblesimulations using VOFsolverinOpenFOAM® withsurfacetensionmodelledusingtheSharpSurfaceForcemodel.

The simulations use a computationaldomain of dimensions1× 2 m² wherea bubble ofdiameterequal to 0.5m is initialized such that its center ishalf a meterfromthe bottom boundaryand equidistantfrom theside walls. The com- putational domainis meshed using160 × 320cells, see[25,30,31,66]. It shouldbe pointedout that fluid properties like density,viscosity,accelerationduetogravityandsurfacetensionaresetbasedon[25,66]buttheremainingparametersare

Fig. 13. The change in the dissolved hydrogen concentration (mol/m ³) as the bubble evolves in the case of SCT0. The white line represents the interface which is plotted at α1= 0 . 5 . See supplementary materials for more information.

setequaltovaluesdescribedinTable4.The tworisingbubblecasescanbedifferentiated basedontherelevantCapillary number(Ca)

Ca=

ρ

1

ν

1Ug

σ

^{, (28)}

whereUg is equalto

|

^g

|

^{LandL isthe} characteristiclength (equal to thediameter ofthe initialized bubble)[25,31,66]. The twocases,representedby TC1andTC2,hasan associatedCawhichisequalto0.268and3.571respectively. Thefour boundaries inthe domainare set to zerogradient for

ψ

,

α

1 whereas C_i is setto fixed value (equal to 0mol/m³) at top boundary andremaining boundariesare set to zero gradient. For p_rgh, the top boundary is set to fixed value (equal to 0mol/m³)andremainingboundariesaresettozerogradient.For^{,thetopandbottomboundariesaresettozerogradient} whereas leftandrightboundariesareassignedtobe fixedvalue (equalto0V). TheboundaryconditionsforU aresetas no-slipattopandbottomboundarieswhereasleftandrightboundariesaresetasslipwalls.

Thecase,TC1,whichcorrespondstoflowwithCa=0.268,asthesurfacetensioneffectsaredominanttherisingbubble deforms toan ellipsoidalshape,seeFig.6a. Ontheotherhand,TC2whichhasa higherCapillarynumber,i.e.Ca=3.571, (when compared toTC1) resultsinforming askirted bubblewithfilamentsalong thesides, seeFig.6b.The final bubble shapeinboththesecasesareobtainedatt=3sanditcomparedagainstbubblemorphologyreportedby[66]and[25]in Fig.6.Thepredictedbubblemorphologyforboth2Drisingbubblesimulationsagreereliablywiththeresultsfrom[25,66]. The deviationsbetweenthepredictions fromthe proposedsolver andthe datareportedinliterature canbe attributedto spuriousvelocities,seediscussionsin[25,30,31].Itisalsoworthpointingoutthatthedeviationbetweentheproposedsolver anddatareportedin[25],usingtheSSFmodel,canalsobe attributedtospurious velocitiesasitsevolutionisdependent on the maximum time step used, which in the former is based on Courant number (equal to 0.05) whereas the latter used a time step constraintproposed by [21,67]. The effect oftime step used in the simulations is further discussed in AppendixD.Anotherpossiblereasonforvariationinthespuriousvelocitiesbetweenthepresentandpreviousworkisthe differentvaluesofsharpeningcoefficient whichissetequalto0.3and0.5respectively(Eq.11),see[26] fordiscussion on theinfluenceofC_shonthesenumericalartifacts.

Fig. 14. The change in the dissolved hydrogen concentration (mol/m ³) as the bubble evolves in the case of SCT1. The white line represents the interface which is plotted at α1= 0 . 5 . See supplementary materials for more information.

0 0.005 0.01 0.015 0.02

0 0.05 0.1 0.15 0.2 0.25

Time (s) U rise (m/s)

SCT0 SCT1

Fig. 15. The temporal variation of the rise velocity for SCT0 ( θ= 15 ^◦) and SCT1 ( θ= 30 ^◦) with the detachment time ( τd) represented using dashed lines (blue for SCT0 and red for SCT1).