Adaptive moving mesh algorithm based on local reaction rate Heliyon

(1)

Contents lists available atScienceDirect

Heliyon

journal homepage:www.cell.com/heliyon

Research article

Adaptive moving mesh algorithm based on local reaction rate

Viktória Koncz

^a,∗

, Ferenc Izsák

^b

, Zoltán Noszticzius

^a

, Kristóf Kály-Kullai

^a

aDepartmentofPhysics,BudapestUniversityofTechnologyandEconomics,1521Budapest,Hungary

bDepartmentofAppliedAnalysisandComputationalMathematics,MTAELTENumNetResearchGroup,EötvösUniversity,1518Budapest,Hungary

A R T I C L E I NF O A B S T R A C T

Keywords:

AdaptiveFEM Movingmesh

Reaction-diﬀusionsystems Acid-basediode COMSOLMultiphysics

Anempiricalmeshadaptionalgorithmisintroducedformodelingone-dimensionalreaction-diffusionsystems withlargemovinggradients.Ournewalgorithmisbasedontherevelation,thatinreaction-diffusionsystems thehighmovingconcentrationgradientsappearnearbytotheregionwheretherateofreactionismaximal, thusthelocalreactionratecanbeusedtocontrolthemeshadaption.Wefound,thatthemainadvantage of sucha method isitssimplicityandeasyimplementation. Asanexamplewestudyanacid-base diode, wherelargemovinggradientsappear.ThemathematicalmodelofthediodecontainsseveralparabolicPDEs, coupledwithoneellipticPDE.An𝑟-refinementtechniqueisusedandattachedtothecommercialfiniteelement solverCOMSOL.Weinvestigatedthetime-dependentsalteffectsofthediodewithourdevelopedalgorithm.

Ourmeshadaptionmethodisadvantageousformodelingofanyreaction-diﬀusionsystemswithlocalized high concentrationgradients.

1. Introduction

Acid-basediodes,whicharecomplexreaction-diffusion-ionicmigra- tion systems, couldtheoretically be used asion-sensors. They could evenformpartsofso-calledlabonachipelectroanalyticaldevicesbe- causeoftheirsimplicity.Wewouldliketomaterializesuchsensorsas ourfinalaim.Toachievethisaim,itisimportanttostudyandmodel thetimeevolutionoftheeffects,whichwouldbethebasisoftheirap- plicationasionsensor.

Inanelectrolytediodeahydrogelcylinderconnectsanalkalineand anacidicreservoir.This waydiﬀusion,migrationandreactionofthe componentsareallowed,butadvectionandconvectionarerestrained, unlikeinpurewater.Duetotheappliedelectricpotentialbetweenthe edgesofthegelcylinder,hydrogene(H⁺)andhydroxil(OH⁻)ionscan getinside theconnectingelement, wherewateris formedfrom their reaction.

Itmeans,thatsomewhereinthemodeleddomainthereisareaction zone,wherethischemicalacid-basereactionoccurs.Thisreactionzone ischaracterizedbylargeconcentrationandpotentialgradients.Dynam- icalsimulationsoftransientbehaviorinsuchsystemscanbedifficult duetothelargemovinggradientsofthemodelvariables(concentra- tionandpotential).Themodeledreaction-diffusionsystemisdescribed byasetofparabolicpartialdifferentialequations(PDEs,in thiscase

*

Correspondingauthor.

E-mailaddress:viki.koncz@gmail.com(V. Koncz).

massbalanceequationsforthechemicalcomponents)andoneelliptic partialdiﬀerentialequation(Poisson-equation).

After the spatial discretization of the domain a system of ordi- narydifferentialequationsisobtained,whichisintegratedforwardby timenumerically(methodof lines).Regardlessthespatialdiscretiza- tionmethod(e.g.finitedifferencemethod,finiteelementmethod)the appliedmeshhasspecialimpactonthenumericalsolvability,speedof convergence,andapproximationerror.Thehigherthenumberofgrid points used,the smallertheerror willbe.However, havingadense gridiscomputationallyexpensive,demandingalargechunkofmemory andanextendedruntime(Farmagaetal., 2011,Zimmerman,2006).

Ifthesolutionshaveregions ofhighspatialactivity,astandarduni- formfixed-gridtechniqueiscomputationallyinefficient,sincetoafford anaccuratenumericalapproximationandreachconvergence,itshould contain,ingeneral,averylargenumberofgridpoints.Themeshneeds tobelocallyrefined.Moreover,iftheregionsofhighspatialactivityare movingintime,likesteepmovingfrontsinreaction-diffusionsystems, thenthemeshshouldbeadaptedaccordingly.

Regardingtothespatialgridadaptionstrategies,atleastthreedif- ferentmethodscanbedistinguished:

• 𝑟-reﬁnementorthemovingmesh(relocation)method:inthiscase ﬁxednumber ofgrid points isused, andthe positionsof points move.

https://doi.org/10.1016/j.heliyon.2020.e05842

Received2October2020;Receivedinrevisedform19November2020;Accepted21December2020

2405-8440/©²⁰²⁰^Published^byÊlsevier^Ltd.^Thisîsânôpenâccessârticleûnder^the^CC^BY-NC-ND^license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

(2)

•ℎ-refinementorthelocalrefinementmethod:thenumberofgrid pointsisnotfixed.Pointsareaddedtoorremovedfromthegrid, accordingtothelocalrequirements,theoriginalgridpoints stay inthesamelocation.Practically,thespatialgridstepℎislocally decreasedorincreased.

•𝑝-reﬁnement:inthisapproachaﬁxedmeshisused,andthepoly- nomialdegreeofbasis(𝑝)alters.

Certainly,thecombinationofthesemethodsexists(Schwab,1999).

Although the moving meshmethod has been less popular than the localrefinementmethod,whichisduetothedifficultyofderivingsuit- ablegoverningequationsfortheadaptivemapping,ithassomeuseful featuresinnumericalcomputations.E.g.itis easytoimplement,and coupletoanyfixedmeshsoftware,anddifficultiesfromthechangein numberofgridpointscanbeavoided(Pˇribyletal.,2006).

Usuallythemovingmeshmethodsarebasedontheequidistribution principle(Huangetal.,1994),meaningselectingmeshpointsinsucha way,thatanempiricalerrorisequalizedovereachsubintervalonthe mesh.Theerrordistributionisfollowedbyamonitorfunction,which isderivedfromtheﬁrstandsecondderivativesofdependentvariables.

(Oftenthisdoesnotresembletonumericalerrors.)

Intheliteraturetherearemanyexamplesaboutthesuccessfulnu- mericalmodelingofsimilarPDEsliketheonesdescribingthebehavior ofthediode.Mostofthemethodsrequirethecompletedevelopment andfullimplementationofthespatialdiscretizationandtimeintegra- tionof themodelequations.A movingmeshmethodwas developed (ZegelingandKok,2004) tosolvetheequationsofreaction-diﬀusion systems, namelytheGray-ScottandBrusselator models.Asystemof adaptivemeshPDEsdescribingthemeshmovementissetupandsolved (inthiscasedecoupledfromtheoriginalPDEsystem),infactthesolu- tionistherealizationofacoordinatetransformationbetweenphysical andcomputationalcoordinates. Theadaptive meshPDEsarederived fromtheminimizationofaso-calledmesh-energyfunctional.Thevari- ationalapproachofthemovingmeshmethodsisgeneral(Caoetal., 2003). Usually the original PDEs are coupled with the mesh PDEs (ZegelingandKeppens,2001).Moreexamplesaboutan𝑟-methodare foundin(Coimbraetal.,2004),(ZegelingandBlom,1992).

Anℎ-methodwasdevelopedforsolvingelectrochemistryproblems withlargespatialgradientsapplyingafinitedifferencemethod(Bieni- asz,2000).Laterthismethodwassuccessfullyusedforthesimulationof theNernst-Planck-PoissonequationsandNernst-Planckequationwith theassumption of electroneutrality((Bieniasz,2004a) and(Bieniasz, 2004b)).More examples abouttheℎ-refinement arefoundin (Lang, 1998,Wangetal.,2004,TrompertandVerwer,1991).

Probably themost significant adaptivemesh techniqueregarding ourworkis thealgorithmsuitableformodeling thetransientbehav- iorofthesamesystemwithcompletelydifferentboundaryconditions (Pˇribyletal.,2006).ThismethodisattachedtothecommercialFEM- LAB routines, which arebased on the finiteelement method. Their monitorfunctionisderivedfromthespatialderivativesofhydrogen,hy- droxylionconcentrationsandthepotential.Duetothetime-consuming evaluationofthemonitorfunctionandsuccessivemeshadaptation,the evaluationisnotcarriedoutateachtimeintegrationstep,theadapta- tionandconsecutiveevaluationtimeintervalisdeterminedbasedon theestimationoftransporttimes.

1.1. Aimofthiswork

Inthispaperwewouldlike toreportourempiricalmovingmesh method (𝑟-refinement), which was developed to model the time- dependentsalt-effectsofanelectrolytediode.Thealgorithmisattached tothecommercialCOMSOLequation-basedmodelingsoftware,which providesthespatialandtimediscretizationFEMalgorithmswithwell- definedinterfaces.Themainideaofthismethodis,thatinareaction- diffusionsystem,wherereactionfrontsmove, thelocalreactionrate indicatingthepresenceofanychemicalreactioncanbethebasisofan

Fig. 1.Schematicviewofastrongacid-basediode:Analkalinereservoir(ﬁlled bye.g.0.1MKOHsolution)andanacidicone(ﬁlledbye.g.0.1MHClsolution) areconnectedbyahydrogel.Δ𝜑>0meansreversebiaseddiode,theopposite isforwardbiased.

adaptivemeshreﬁnement.Theproblematicpartofsuchasystemisthe neighborhoodofthelocalchemicalreactions,wherethegradientsand thenumericalerrorofapproximationarethehighest.Thisideacould bethebasisof anyspatialdiscretization method, notonly theﬁnite elementmethod,theadvantageswouldbethesame:easyimplementa- tion,simplemonitorfunction,hencetherequiredmeshdensitywould bedeterminedinanempiricalwaynotdirectlyderivedfromtheerror ofapproximation.

2.Numericalsimulations 2.1. Themodel

Theso-calledacid-basediode(orelectrolytediode)isareaction- diffusion- ion-migrationsystem(NoszticziusandSchubert,1973,Schu- bertandNoszticzius,1977),where,duetotheapplied potentialand concentrationgradients,thechemicalspeciesmove,whileanacid-base reactionoccurs. Previouslysome highlynonlinearphenomenacalled positiveandnegativesalteffectsrelatedtothesalt-contaminationof theacid-basediodewerediscovered,thusitisstraightforwardtouse thenegativesalteffecttosensitivelydetectnonhydrogencations(Chun andChung,2015,Roszoletal.,2010,Kirschneretal.,1998).

Inanacid-basediodeaconnectingelement(eitherahydrogelora membrane)connectsanalkaline(KOH)andanacidic(HCl)reservoir, restrainingconvection,butallowingdiﬀusion,migrationandchemical reaction.(Theschematicviewofanacid-basediodeisshowninFig.1) Iftheacidic reservoirhasa negativepotentialin comparison tothe alkalineone(forwardbiaseddiode),fromthealkalinereservoircations (e.g.K⁺),fromtheacidiconeanions(e.g.Cl⁻)canmigrateintothegel cylinder,wheretheyformawell-conductingsalt(inthiscaseKCl).If theacidicreservoiristhepositiveside(reversebiaseddiode),hydrogen (H⁺)andhydroxil(OH⁻)ionscanmoveintothehydrogel,wherewater isformedfromtheirreaction,resultinginmuchsmallerconductivity.

Mostgelscontain fixedweaklyacidicgroups.Thedissociationof themleadstofixednegativeions,inthepresentworkFA⁻isthefixed charge,andHFA denotesits protonatedform.The current - voltage characteristicofanelectrolytedioderesemblestothecharacteristicofa semiconductorone,meaningthediodecurrentdependsonthepolarity oftheappliedvoltage,asFig.2shows.

Ifthereservoirsofthediodearecontaminatedbyasalt(inthecase ofstrong diodeby KCl),theconcentrationof thecontaminating salt inﬂuencesthecurrentofthereversebiaseddiode:

• If only one reservoir is contaminated (usually the alkaline one duetohighersensitivity),thediodecurrentincreasesinanonlin- earway(Fig.3/a).Thisphenomenoniscalledpositivesalt-eﬀect ((Kirschneretal.,1998,Heged˝usetal.,1999)).

• Ifonesideofthediodeisalreadycontaminated,andtheopposite reservoirofthediodegetscontaminatedaswell,surprisinglythe

(3)

Fig. 2.Simulatedvoltage–currentdensitycharacteristicofastrongacid-base diode.Concentrationsoftheelectrolytes:[KOH]=[HCl]= 0.1M.

Fig. 3.Calculatedcurrentdensityasthefunctionofsaltcontamination.(a):

positivesalteffect,(b):negativesalteffect.[KCl]_ALand[KCl]_ACdenotethesalt concentrationinthealkalineandacidicreservoirs,respectively.Theblackdot insubfigure(a)representsthealkalinesaltconcentrationusedforsubfigure(b).

[KCl]_AL= 60mM,Δ𝜑= 10V.

diode currentdecreases byincreasedcounter saltcontamination (Fig.3/b).Thisphenomenoniscallednegativesalt-eﬀect.

Bothpositiveandnegativesalteﬀectshadalreadybeeninvestigated for stationarystate (Kirschneret al., 1998, Roszoletal., 2010), we startedtofocusonthetime-dependentmodelingofsalt-contaminated diode (Konczetal.,2017),wherewemetunexpected computational complexity.

Atthisstageofinvestigation,themathematicalmodelofthediode isasystemofone-dimensionalpartialdiﬀerentialequations(Roszolet al., 2010, Merkinet al., 2000, Heged˝uset al., 1996, Lindneret al., 2002).Thepreviouslydiscussedsimpliﬁcationsandassumptions((Iván etal.,2002))arekept(onedimensionalmodel,concentrationpolariza- tionandohmicpotentialdroponthereservoirareneglected).

Tocalculatetheconcentrationdistributionofthe𝑖-thcomponentthe mass-balanceequationisused,andtheNernst-Planckequationisused toobtainthecurrentdensity:

𝜕𝑐_𝑖

𝜕𝑡 = −𝜕

𝜕𝑥 (

−𝐷_𝑖𝜕𝑐_𝑖

𝜕𝑥−𝑧_𝑖𝐷_𝑖𝐹 𝑅𝑇 𝑐_𝑖𝜕𝜑

𝜕𝑥 )

+𝜎_𝑖, (1)

where𝑐_𝑖,𝑧_𝑖and𝐷_𝑖denotetheconcentration,chargenumberanddiﬀu- sioncoeﬃcientofthe𝑖-thcomponent(H⁺,OH⁻,K⁺,Cl⁻,FA⁻,HFA), respectively.𝑇 isthetemperature,𝐹 istheFaradayconstant,𝑅isthe molargasconstant,𝜑denotesthespacedependentpotential, 𝑡isthe timecoordinateand𝑥isthespacecoordinatealongtheconnectingele- ment.𝜎_𝑖denotesthereactionrateofthe𝑖-thchemicalcomponent:

𝜎H⁺=𝑘w(𝐾w−𝑐H⁺𝑐OH⁻) +𝑘f(𝐾f𝑐HFA−𝑐H⁺𝑐FA⁻), (2)

𝜎OH⁻=𝑘w(𝐾w−𝑐H⁺𝑐OH⁻), (3)

𝜎K⁺=𝜎Cl⁻= 0, (4)

𝜎FA⁻= −𝜎HFA=𝑘f(𝐾f𝑐HFA−𝑐H⁺𝑐FA⁻), (5) where𝐾w isthe equilibriumconstantof water dissociation,and𝑘w

denotestherateconstantofwaterrecombination.𝐾fisthedissociation constantand𝑘fdenotesthekineticconstantoftheﬁxedcharge.

Poissonequation is used tocalculate thepotential proﬁleof the diode:

𝜕

𝜕𝑥 (

−𝜀𝜕𝜑

𝜕𝑥 )

=𝐹( ∑

𝑖 𝑧_𝑖𝑐_𝑖)

, (6)

where𝜀denotesthepermittivityofwater.

InourmodelcalculationsDirichletboundaryconditionswereap- plied.Donnanequilibriumisusedtocalculatetheconstantboundary concentrationsinsidetheconnectingelementfromtheﬁxedconcentra- tionsinthereservoirs.TheDonnanpotentialwasusedtocorrectthe electricpotential.

Unfortunatelywewerenotabletoreachconvergence,iftheabove discussedfullsystemofdiﬀerentialequationsisapplied (massbalance equationsaresetupforeachspecies).Forstationaryinvestigationsthe concentrationof theﬁxedanions canbecalculatedfromtheequilib- rium,

𝑐FA⁻= 𝑐f

𝑐H⁺∕𝐾f+ 1, (7)

where𝑐fdenotestheestimatedtotalconcentrationoftheﬁxedgroups (Ivánetal.,2004),thereforethemass-balanceequationisusedforthe mobileions(H⁺,OH⁻,K⁺,Cl⁻)only.

This assumption of quasi-stationarity was considered to be valid evenintimedependentcasesTheusageofthissimplification means thattheprotonationanddeprotonationofthefixedchargegroupsare consideredinstantaneous,thusthereactionratebelongingtothefixed groupsisapproximatedbyzero(𝜎FA⁻=𝜎HFA= 0).Theremainingreac- tionrateofH⁺is:𝜎H⁺=𝜎OH⁻=𝑘w(𝐾w−𝑐H⁺𝑐OH⁻).

Withthissimpliﬁcationitbecomespossibletosolvetheremaining PDEswithCOMSOL.

TheappliedvaluesoftheparameterscanbefoundinTable1.

2.2. SolutionmethodwithCOMSOL

ThesimulationswerecarriedoutbyusingthecommercialFEMsoft- warepackageCOMSOLMultiphysics4.3.Themodelwasdeﬁnedfrom thefollowing interfaces: for the Nernst-Planck equation the “Trans- portofDilutedSpecies”and“ClassicalPDEs/Poissonequation”were applied.Afullmodelcalculationcanbesplitintotwomainsteps:af- tercomputingtheboundaryconditions,theinitialconditions(thisisa stationarystateofthediode)arecalculatedbyapreviouslydiscussed parametricstationarysolver(Roszoletal.,2010),whichisfollowedby thetime-dependentstudy.Thesolversettingsareadjustedunderthe

“Study”node,wherethetwostudystepsbelongingtothesolversare deﬁned.Betweenthetwostudystepstheboundaryconditionsmustbe resetforthetransientvalue.(Inthecaseofstep1Computetoselected, inthecaseofstep2Computefromselectedcommandmustbeused.)

Forthenumerical integration of thetime-dependent ODEs(after spatialdiscretization) BDF method with SPOOLES linear solverwas appliedwiththefollowingsettings:automaticscaling,thetime-steps aredetermined freelyby thesolver. Toavoidnumerical oscillations crosswinddiﬀusionstabilizationtechniquewasapplied.Secondorder Lagrangeelementswereused.

ThekeyofFEMsimulationsistheappliedmesh.Thegradientsofthe concentrationsandpotentialarethelargestintheneighborhood ofthe so-calledreactionzone,thisisthehardesttocalculatepartofanacid- basediode.Toachieveconvergence,inthiszoneextremelyﬁnemesh isrequired.Duringthetime-dependentstudythelocationofthereac- tionzoneischanging,henceanadaptivemeshingtechniqueishighly recommendedtoreducecomputationalcosts.

In COMSOL Multiphysics adaptive mesh reﬁnementextension is availableforthetime-dependentstudies.Unfortunatelythisbuilt-intool

(4)

Table 1.Modelparameters(Marcus(1997)).

𝑐KOH(moldm⁻³) Concentration of base 0.1

𝑐HCl(moldm⁻³) Concentration of acid 0.1

𝐷H⁺(m²s⁻¹) Diffusion coefficient of protons 9.31 × 10⁻⁹ 𝐷_OH⁻(m²s⁻¹) Diffusion coefficient of hydroxyl ions 5.28 × 10⁻⁹ 𝐷K⁺(m²s⁻¹) Diffusion coefficient of potassium ions 1.96 × 10⁻⁹ 𝐷Cl⁻(m²s⁻¹) Diffusion coefficient of chloride ions 2.04 × 10⁻⁹ 𝐷_FA⁻(m²s⁻¹) Diffusion coefficient of the fixed charge 0 𝐷HFA(m²s⁻¹) Diffusion coefficient of the protonated fixed group 0

𝜑L(V) Electric potential on the L boundary 0

𝜑_R(V) Electric potential on the R boundary 10 − 20

𝑙(m) Length of the gel 1 × 10⁻³

𝑧_H⁺ Charge number of protons 1

𝑧OH⁻ Charge number of hydroxyl ions −1

𝑧K⁺ Charge number of potassium ions 1

𝑧_Cl⁻ Charge number of chloride ions −1

𝑧FA⁻ Charge number of the ﬁxed charge −1

𝑧HFA Charge number of the protonated ﬁxed charge 0

𝐹(Cmol⁻¹) Faraday constant 9.6487 × 10⁴

𝑅(Jmol⁻¹K⁻¹) Molar gas constant 8.314

𝑇(K) Temperature 298.15

𝜀(AsV⁻¹m⁻¹) Permittivity of water 6.954 × 10⁻¹⁰ 𝑘_w(m³mol⁻¹s⁻¹) Kinetic constant of water recombination 1.3 × 10⁸ 𝐾w(mol²m⁻⁶) Ionic product of water 1 × 10⁻⁸ 𝑐f(molm⁻³) Concentration of the fixed charge 4 𝑘f(m³mol⁻¹s⁻¹) Kinetic constant of the fixed charge 6 × 10⁶ 𝐾f(molm⁻³) Dissociation constant of the fixed charge 1 × 10⁻¹

Fig. 4.The general operation of adaptive meshing algorithms.

andtheCOMSOL’sMovingMeshinterfaceprovedtobeinadequatefor ourpurposes.Thusweattachedanownadaptivealgorithmtooursim- ulationframework.(CarryingoutmodelcalculationsinCOMSOLGUI inseriesiscumbersome,thereforeaneasy-to-useframeworkwasdevel- opedinJavausingtheCOMSOLJavaAPI.)

ForthegeneraloperationofadaptivemeshingalgorithmsseeFig.4.

Thenumericalsolvermustbemonitored,andattheendofeveryiter- ationstepitmustbedecided,whethertheiterationcanbecontinued (meshisstillsatisfying)oranewmeshshouldbegenerated.Incaseofa newmeshthelastsolutionwiththeoldmeshisinterpolatedtothenew one,andthenumericalsolverisrestarted.InCOMSOLtheFEMalgo- rithmsareconsideredtobe“blackboxes”,thatisthenumericalmethod canbereachedonlyviatheinterfacesprovidedbyCOMSOL.Tostop andmonitorthenumericalsolver,astopconditionisformulatedapply- ingtheso-called“DomainPointProbe”.Withapointprobeexpression afunctionofvariablesofthediﬀerentialequationatapre-deﬁnedlo- cationisevaluated.

Atthetwoedgesofthereactionzone,insidethedensemeshzone, twopointprobeexpressionsaredeﬁned.Thevaluesof themarethe absolutevalueoftheacid-basereaction’slocalreactionrateatthegiven location:

𝑝𝑝𝑣=||𝜎H⁺||=||𝜎OH⁻||=||𝑘w(𝐾w−𝑐H⁺𝑐OH⁻)||. (8) Thenthestopconditionisthefollowing:ifoneofthepointprobeex- pressionsdiﬀersigniﬁcantlyfromzero(𝑝𝑝𝑣>1.3⋅10⁻³^mol

dm³s),thesolver

isstopped.This𝑝𝑝𝑣valuemeansthesolverstopswhen𝑐H⁺⋅𝑐OH⁻> 𝐾w, wegetthisthresholdbytrialanderror.Afterthat,fromthebehavior ofthereactionzoneanewmeshisgenerated.WithCOMSOLJavaAPI inthecaseofone dimensionalgeometryonlythevertexcoordinates arecalculated, theedges aretrivial(how thevertexesareconnected toanelement).InCOMSOLthesolutionofthetime-dependentsolver iscopied(Copysolution)toretaintheresults,andthedomainpoint probesarerelocatedbeforethetime-dependentsolverisrestarted.

2.3. Adaptivealgorithmwithgridfunction

Theso-calledgridfunctionwasselectedtodeﬁnethemesh.An𝑥→ 𝑓(𝑥)gridfunctiondeﬁnestherealmesh’scoordinatesasafunctionof theequidistantmesh’scoordinates.Somegridfunctionswithschematic onedimensionalmeshesareshowninFig.5.

Thebasisofourselectedgridfunctionisanequidistantmeshwitha reﬁnedzone(Fig.5/b),wherethevicinityofthereactionzoneiscov- eredbyanextremelyﬁnemesh.Toavoidnumericalproblemsarising fromsuddenchangesinthesizeofneighboringelements,inthevicin- ityofthestraightline’sintersectionsexponentialsmoothingfunctions areapplied.Formoredetailsseethe1stsectionoftheSupplementary Information.

Theadaptivealgorithmbasedonthegridfunctionhasthefollowing parameters:

• 𝑁:totalnumberofmeshpoints

• 𝑦A,𝑦B:inreal-meshcoordinatestheinitiallimitsoftheextremely ﬁnemeshzone

• 𝑅:ratioofthenumberofpointsbetween𝑦Aand𝑦Bandthetotal numberofgridpoints

• 𝛿: widthof theexponentialsmoothing partin equidistant mesh coordinates

• 𝑑: distancebetween 𝑦A or 𝑦B andits correspondingpointprobe expression

• Δ𝑦:howmuchthecontrolpoints(𝑦A,𝑦B)ofthegridfunctionare shiftedintothedirectionofthereaction’smovement,afterthestop conditionisfulﬁlled,andthenextmeshisgenerated

Themainadvantageofapplyingagridfunctionis itseasyimple- mentation, simple mesh generation, and low time-demand of mesh

(5)

Fig. 5.Somegridfunctionswiththecorrespondingschematicmeshillustrations.(a)Equidistantmesh.(b)Equidistantmeshwithareﬁnedzone.(c)Theapplied gridfunctionwithexponentialsmoothing.

construction. Forthe source codesee https://github.com/vikikoncz/ adaptive_FEM_based_reaction_zone.

3. Resultsanddiscussion

Wewouldliketoshowtheabilityoftheadaptivemeshalgorithm tosimulatethetime-dependentbehaviorofanacid-basediode.Inthe presentworkwedemonstrate theoperationofthealgorithmincase of time-dependentpositive salteﬀect withthefollowingsettings. As initialstatethesalt-freestationarydiodewasused.Thecontaminating salt was addedtothealkaline reservoirat 𝑡= 0,[KCl]_AL,t= 60mM.

Approximatelyafter100sthenewsteadystatewasattained,thus150 swassimulatedafterthesystemwasperturbedbyKClcontamination.

Ourworkispresentedinthreesubsections.First,inSection3.1the appearanceofthereactionzoneisinvestigated.Wewouldliketoshow, that thelocation of thechemical reaction (where thereaction ratio diﬀerssigniﬁcantlyfromzero)isveryclosetothelocationwherethe gradientsoftheconcentrationandpotentialarethehighest.

Then,weinvestigatedtheabsoluteandrelativeerrorofapproxima- tiononthemodeleddomainforvalidationpurposes(Section3.2).

Finally,we discussthe eﬀectof diﬀerent meshes onthe error of approximationandsimulationtime(Section3.3).

3.1. Reactionzone

Thelocalreactionrateoftheacid-basereactionis𝜎H₂O= −𝑘w(𝐾w− 𝑐H⁺𝑐OH⁻).Itisthelinearfunctionof[H⁺]⋅[OH⁻],whichinequilibrium is atemperaturedependentconstant.Fig.6shows the localreaction rate(𝜎H₂O)alongthegelcylinderinastationarysalt-freediode.𝜎H₂O

Fig. 6.Thelocalreactionrate(𝜎H2O= −𝑘w(𝐾w−𝑐H⁺𝑐OH⁻))ofwaterrecombina- tionintheuncontaminateddiode.

isdiﬀerent fromzeroonlyat thereaction zone,wheretheacid-base reactionoccurs.

Fig.7showssomecalculated𝜎H₂Oprofilesofatime-dependentpos- itivesalt effect. Atacertaintime𝜎H₂Ois differentfrom zeroinone regiononly,certainlythisisthereactionzone.

Asthealkalinereservoir(lefthandside)isapproachedbythereac- tionzone,theacid-basereactionbecomesfasterandfaster(thepeakof 𝜎H₂Oincreases).Thereactionzonekeepsitspositionafterabout10s, however,thereactionratecontinuesincreasing.

Thelookofthereactionzoneatfourselectedmomentsoftimeis showninFig.8.Whiletheacid-basereactiongetsfaster,theextension ofthereactionzoneshrinks.Inthenewsteadystateitswidthisonly

(6)

Fig. 7.Thelocalreactionrate(𝜎H2O= −𝑘w(𝐾w−𝑐H⁺𝑐OH⁻))ofwaterrecombina- tionatdiﬀerentmomentsoftime.

thethirdoftheuncontaminatedone.Theextensionofthereactionzone iscalculatedfromthehalf-widthofthe𝜎H₂Opeak.

Intheuncontaminatedsteadystatethepositionsoftheextremaof spatialderivativesarevery closetothelocationwheretheacid-base reactionoccurs.Thesederivativesmovetogetherwiththechemicalre- action.Thelocalextremaofconcentrationderivativesareincreasing(in absolutevalue),whilethederivativeofpotentialchangesonlyslightly.

Thesituationisthesameinthecaseofsecondspatialderivatives.

3.2. Validation

Inthecaseofasystemofnonlinearpartialdifferentialequations, itispossible,thatthenumericalsolverfindsaphysicallyinvalidso- lutioneventhoughitsatisfiesthetolerancecriteria(LeVeque(1992)).

Weusedaposteriorierrorindicatorstocharacterizetheprecisityofthe numericalsolutions.

Ameaningfulmeasureofthecomputationalerrorcanbegivenon (𝑥𝑖,𝑥𝑖+1)asfollows:

𝑥_𝑖+1 𝑥∫_𝑖

(̃𝑦(𝑥) −𝑦(𝑥))² 𝑑𝑥=𝑅(𝑥,𝐚), (9)

where𝑦(𝑥)denotestherealsolution ofthediﬀerentialequationand

̃𝑦(𝑥)=∑_𝑛

𝑗=1𝑎_𝑗⋅Φ_𝑗(𝑥)istheﬁniteelementapproximation.Inpractice, therealsolutionisunknown,therefore,anothererrorindicatorhasto befound.

Ifthecoeﬃcients𝑎_𝑗ofthebasisfunctionsareavailable,substituting themintotheoriginalPDE’sleadstoanerrorindicatoroneachelement.

UnfortunatelyduetoCOMSOL’s‘blackbox’approachitisneither possibletogettheLagrangian coeﬃcients,northeerrorcalculatedby COMSOL.

Tofindametricforcomparingdifferentsolutions,anownmethod hadtobedeveloped.ItispossibletoexportfromCOMSOLtheinterpo- latedsolution,andthespatialandtimederivativesofthevariablesof theequationatthepredefinedpointsinagivenmoment.(FromCOM- SOLJavaAPItheso-called‘Interp’functionwasappliedtoretrievethe numericalresults.)Ourmethodisthefollowing:duringtimedevelop- mentateverydiscretemomentoftimetheapproximatedvalueofthe variables,the1stand2ndspatialderivativesandthetime-derivatives areexportedatthemeshpoints.Thesevaluescanbesubstitutedtothe fivedifferentialequations(fourmassbalancesandthePoisson)atevery momentandlocation,leadingtofive(absolute)residualvalues𝑒𝑟𝑟_abs,𝑖. Asanexample,thecalculatedabsoluteresidualvaluesoftheH⁺’s massbalanceequation,andthePoissonequationareshowninFigs.9 and10.

Theabsolute error in itselfdoesnot tellusthe wholepicture,it hastobecomparedtotheabsolutevaluesofthevarioustermsinthe givenequation.Forthisreasonweusetherelativeerrordeﬁnedinthe followingway:

𝑒𝑟𝑟rel,𝑖= 𝑒𝑟𝑟abs,𝑖

|||^𝜕𝑐^𝜕𝑡^𝑖|||+||

||−^𝜕

𝜕𝑥

(

−𝐷_𝑖^𝜕𝑐_𝜕𝑥^𝑖)|

|||+||

||−^𝜕

𝜕𝑥

(

−^𝑧^𝑖^𝐷^𝑖^𝐹

𝑅𝑇 𝑐_𝑖^𝜕𝜑_𝜕𝑥)|

|||+||𝜎_𝑖|| (10) 𝑒𝑟𝑟rel,Poisson= 𝑒𝑟𝑟abs,Poisson

||||^𝜕𝑥^𝜕 (

−𝜀^𝜕𝜑_𝜕𝑥)|

|||+||

||𝐹( ∑

𝑖𝑧_𝑖𝑐_𝑖)|

|||

(11)

Despiteofthelarge absoluteerrorintheneighborhoodofthere- actionzone,therelativeerrors (Figs.11and12) seemtobe mostly acceptable.However,averysharpregionexistsinthemiddleof the reactionzone,wheretheapproximationerrorofthePoissonequation becomes0.65.

3.3. Theeﬀectofthegridfunction’sparametersontheerrorof approximationandonthesimulationtime

Tocomparethediﬀerenttypesofmeshes,thesimulationtimemea- suredinanisolatedvirtualmachine(fordetailsseethe2ndsectionof SupplementaryInformation)wasmonitored,andtheerrorofapproxi- mationwasestimated.Unfortunately,thenumberofiterationsarenot measurableviaCOMSOL,thuswehadtomeasuretherunningtimeof thesimulations.It isalso importantinpractice, butdepends on the hardware.ThebasicstepsofanadaptiveFEMsimulationwiththemea- suredrunningtimevaluesarefoundinthe4thsectionofSupplementary Information.About97%ofthetotalsimulationtimegoestothetime- steppingandmeshinterpolation, thus therunningtimeoftheother stagesofFEMsimulationisnegligible.

Insteadofplottingtheresidualsateverymomentoftime(likein Section3.2),itisbettertouseacumulatederrorindicatortodescribe theapplicabilityofthemesh.Toconstructsuchanaposteriorierror indicator,weconsiderequation (1)asasystemofconservation laws coupledwiththereactivesourceterms(𝜎H⁺,𝜎OH⁻,0,0).Usingthecor- respondingtheoryforthenon-linearconservationlawsinTheorem5.2 in(Cockburn,1999),astraightforwardchoicefortheerrorestimation ofthe𝑗-thcomponentinequation(1)canbecalculatedasfollows:

𝑒𝑟𝑟𝑗=

𝑇

∫

0 𝐿

∫

0

|𝑒𝑟𝑟abs,𝑗|𝑑𝑡𝑑𝑥+

√√

√√max

𝑡∈[0,𝑇]𝑆(𝑐𝑗)(𝑡)⋅

𝑇

∫

0

𝑆(𝑐𝑗)(𝑡)𝑑𝑡⋅√

𝐷𝑗, (12)

where 𝑒𝑟𝑟abs,𝑗 denotes the𝑗-th componentscalculated residual, and 𝑆(𝑐_𝑗)(𝑡)isthesumoftheconcentrationabsolutediﬀerences between thegridpointsontheentiremodeleddomain.Wefound,thattheterm containing𝑆isalmostnegligiblecomparedtothetimeandspatialin- tegraloftheresiduals.Furthermore,itwasfoundthatthistermvaries within1%inthecaseofeverymeshsetting(ifconvergenceisachieved), thusitcanbeneglectedwhencomparingmeshes.

Forasinglecumulativeerrorindicatortheerrorestimationofthe fourcomponentsissummed andnormedby𝑇= 150s (duetothetime integration):

𝑒𝑟𝑟𝑜𝑟=𝑒𝑟𝑟H⁺+𝑒𝑟𝑟OH⁻+𝑒𝑟𝑟K⁺+𝑒𝑟𝑟Cl⁻

𝑇 . (13)

Inthefollowingtablesandﬁguresonlythiscumulativevalueisused tocomparediﬀerentmeshes.

First,theefficiencyof theadaptivemeshing algorithmisshowed withcomparisontotwofixedmeshsimulations.Afterthat,theeffects ofthegridfunctionparametersareexplored,andwetriedtofindan optimalparametercombinationconsideringthesimulationtimeandthe errorofapproximation.(Inthisstudytheparametersarefixedduringa simulation.)Finally,weinvestigatedhowthecontinuousmodification ofthegridfunction’sparametersinfluencesthesimulation.

3.3.1. Fixedmeshvsadaptivemesh

Theperformanceofanequidistantmesh(namedfixed(1)),afixed meshwithadensezone(fixed(2))andthreemovingmesheswiththe proposedgridfunctionwerecompared(Table2).Therunningtimeof

(7)

Fig. 8.Reactionfronts:spatialgradientsandthelocalreactionrateofacid-basereactionindifferentmomentsoftime.Inthefirstcolumnthelocalreactionrateand the1stspatialderivativesareshownatt=0s(a),t=3s(c),t=5s(e)andt=20s(g).Inthesecondcolumnthelocalreactionrateandthe2ndspatialderivatives areshownatt=0s(b),t=3s(d),t=5s(f)andt=20s(h),respectively.Pleasenote:thederivativesandthelocalreactionratearescaled,thusnounitisshown inthisfigure,becauseweareonlyinterestedinthepositionsofpeaks,nottheirmagnitude.

Fig. 9.TheabsoluteerroroftheH⁺ionsmassbalanceequation(𝑒𝑟𝑟abs,H⁺)at 𝑡= 3s.Theparametersoftheappliedmesharethefollowing:𝑁= 10000,𝑦A= 0.192,𝑦B= 0.202,𝑑= 0.002,Δ𝑦= 0.002,𝛿= 0.01,𝑅= 0.8.

Fig. 10.Theabsoluteerrorof thePoissonequation(𝑒𝑟𝑟abs,Poisson)at𝑡= 3s.

Theparametersoftheappliedmesharethefollowing:𝑁= 10000,𝑦A= 0.192, 𝑦_B= 0.202,𝑑= 0.002,Δ𝑦= 0.002,𝛿= 0.01,𝑅= 0.8.

(8)

Fig. 11.TherelativeerroroftheH⁺ionsmassbalanceequation(𝑒𝑟𝑟rel,H⁺)at𝑡= 3s inthereactionzone.Theparametersoftheappliedmesharethefollowing:

𝑁= 10000,𝑦A= 0.192,𝑦B= 0.202,𝑑= 0.002,Δ𝑦= 0.002,𝛿= 0.01,𝑅= 0.8.

Fig. 12.TherelativeerrorofthePoissonequation(𝑒𝑟𝑟rel,Poisson)at𝑡= 3s inthe reactionzone.Theparametersoftheappliedmesharethefollowing:𝑁= 10000, 𝑦A= 0.192,𝑦B= 0.202,𝑑= 0.002,Δ𝑦= 0.002,𝛿= 0.01,𝑅= 0.8.

Table 2. The measured running time and calculated errorofapproximation(error)in the caseoffixedandadaptivemeshes. The fixed(1)meshisanequidistantmesh,fixed(2) isa fixedmesh witha densezone. Param- etersofthegridfunctionarethefollowing:

𝑦A= 0.192,𝑦B= 0.202,𝑑= 0.002,Δ𝑦= 0.002, 𝛿= 0.01,𝑅= 0.8.

Mesh type N Time (min) error

ﬁxed(1) 100000 3648 147

ﬁxed(2) 21000 532 73

adaptive 500 15 183

adaptive 1000 35 97

adaptive 2000 42 49

adaptive 10000 222 10

theequidistantmesh(10000gridpoints)is thebiggest,andapprox- imatelythesameprecision can beachieved usinganadaptive mesh with500- 1000gridpointsdependingonthegridfunctionparameters.

Theﬁxed(2)meshhasthesameresolutionastheequidistantmesh (ﬁxed(1))initsdensepart,whichcoversthatpartofthegel,wherethe reactionzonecanoccur.Theremainingpartofthegeliscoveredbyan 80timeslessdenseequidistantmesh,andthetwozonesareconnected byanexponentialsmoothingpart.

Despitethesameresolutioninthereactionzone,theerrorofapprox- imationisalmosthalvedinthecaseofﬁxed(2),duetotheconvergence

Fig. 13.Themeasuredrunningtime andtheerrorofapproximationasthe functionofthenumberofgridpoints.Parametersofthegridfunctionarethe following:𝑦_A= 0.192,𝑦_B= 0.202,𝑑= 0.002,Δ𝑦= 0.002,𝛿= 0.01,𝑅= 0.8.

criteriaofCOMSOL’snonlinearsolver.Intheaggregatederrorevalu- atedbyCOMSOLaftereveryiterationstep,theerroronanelementis notweightedbytheelement’ssize.ThusinCOMSOL’saggregatederror, inthecaseofequidistantmeshtheelementsbelongingtothereaction zoneareunderweightedcomparedtotheﬁxed(2)mesh.(Furtherdetails aboutCOMSOL’sconvergencecriteriacanbefoundinthe3rdsection oftheSupplementaryInformation.)

Themostaccuratesolution isachievable byamovingmeshwith asmany gridpointsasthedesiredrunningtimeallows.(Inourcase maximum10000meshelementswereused.)Inthecaseofonespeciﬁc adaptivemeshparametersettingthemeasuredrunningtimeandthe errorofapproximationasthefunctionofthenumberofgridpointsare showninFig.13.(Meshtrajectoriesofthismovingmeshcanbefound inthe7thsectionoftheSupplementaryInformation.)

Ifthenumberofgridpointsincreases,theapproximationerrorde- creasesinahiperbolicway.Therunningtimecanshowdiscrepancies, itdoesnotincreasemonotonicallyasafunctionof𝑁.

Theconvergenceratesofthesediscretizationschemes (ﬁxed(2),and adaptive)arecomparedinthe6thsectionoftheSupplementaryInfor- mation.

3.3.2. Theeﬀectofthegridfunctionparameters

Inthissectionwepresentourﬁndingsaboutsystematicinvestiga- tionsoftheeﬀectofgridfunctionparametersontheperformanceof themovingmeshalgorithm.Pleasenote,duringonesimulationthepa- rametersofthegridfunctionwerekeptconstant.Duringthediscovery oftheparameters𝑁= 2000waschosen,becausewefounditagood compromisebetweenrunningtimeanderrorofapproximation(seeTa- ble2).

First,theeffectofΔ𝑦wasinvestigated(Fig.14).Itdoesnotinflu- encetheerrorofapproximationsignificantly,however,themeasured simulationtimedepends onit.ThesmallertheΔ𝑦is,thehigher the runningtimebecomes.Thefrequentmeshadaptationneedstoomuch interpolationbetweenthedifferentmeshes.

Withinawiderangethegridfunctionparametercalled𝑅hasvery similareﬀect ontheerror andsimulation time(Fig. 15).If𝑅 istoo small,thesolverjust cannotconverge(becausetheresolutionis not highenoughinthereactionzone).Certainly,theprecisevaluedepends on thenumberof grid points,and theotherparameters of thegrid function,e.g.if𝑁= 2000, 𝑅mustbe atleast 0.05.However,in this casetheapproximationerrorishuge(2600).(Forthisreasonitisnot showninFig.15.)Theerrorofapproximationhasanoptimumaround 𝑅= 0.94,if𝑅islarger,toofewgridpointsappearoutsidethereaction zone.

(9)

Fig. 14.Themeasuredrunningtimeandtheerrorofapproximationasthe functionoftheΔ𝑦parameter.Parametersofthegridfunctionarethefollowing:

𝑦_A= 0.191,𝑦_B= 0.203,𝑑= 0.003,𝑅= 0.8,𝛿= 0.01,𝑁= 2000.

Fig. 15.Themeasuredrunningtimeandtheerrorofapproximationasthe functionofthe𝑅parameter.Parametersofthegridfunctionarethefollowing:

𝑦A= 0.191,𝑦B= 0.203,𝑑= 0.003,Δ𝑦= 0.002,𝛿= 0.01,𝑁= 2000.

Dependingonthewidthoftheextremelyﬁnemeshzone,andon thegridfunctionparameters,anoptimal𝛿exists,wheretheerrorof approximationisthesmallest(Fig.16).If𝛿increasesfromzero,ﬁrstthe errorofapproximationdecreases,duetosmoothertransitionofelement size.When𝛿becomescomparabletothesizeofthedensemeshzone, thenumberofmeshelementsinthereactionzonestartstodecrease, resultinginincreasederrorofapproximation.

Duetothetimedemandoftoofrequentmeshadaption,anddueto convenience,densemeshwasusuallyusedonalargeenvironmentof thereactionzone.(Itmeans,𝑦B−𝑦Aisatleastﬁvetimesbigger,than thereactionzone’sextension.)Theimpactofthedistancebetween𝑦A

and𝑦B wasinvestigated(Fig. 17) toseewhethercoveringasmaller environmentwithdensemeshstillgivessatisfactoryresults.Theerror of approximation increaseswith thisdistance (because lessgrid ele- mentscoverthereactionzone).Usuallysmallererrorcomeswithlower runningtime,however,heretheoppositehappens:increasederrorde- creasestherunningtime.

Theminimalreasonabledistancebetween𝑦Band𝑦Aistheextension ofthereactionzone:underthislimitconvergenceproblemsappear.In thecaseofminimal𝑦B−𝑦Atheerrorandtherunningtimeasthefunc- tionof𝑁isshowninFig.18.Toreachconvergencewiththesesettings, theminimumnumberofrequiredelementsis𝑁= 150,althoughinthis casetheerrorofapproximationisaround450,andtheproﬁlesarenot

Fig. 16.Themeasuredsimulationtimeandtheerrorofapproximationasthe functionofthe𝛿parameter.Parametersofthegridfunctionarethefollowing:

𝑦_A= 0.191,𝑦_B= 0.203,𝑑= 0.003,Δ𝑦= 0.002,𝑅= 0.8,𝛿= 0.01,𝑁= 2000.

Fig. 17.Themeasuredrunningtime andtheerrorofapproximationasthe functionofthedistancebetween𝑦Aand𝑦B.Parametersofthegridfunctionare thefollowing:𝑑= 0.0005,Δ𝑦= 0.0005,𝑅= 0.8,𝛿= 0.01,𝑁= 2000.

Fig. 18.Themeasuredrunningtime andtheerrorofapproximationasthe functionofthe𝑁.Parametersofthegridfunctionarethefollowing:𝑦A= 0.196, 𝑦B= 0.199,𝑑= 0.0005,Δ𝑦= 0.0005,𝑅= 0.8,𝛿= 0.01.Thedensemeshbarely coversthereactionzone.

smoothenough.If𝑁= 500,thealgorithmhasverysimilarperformance, as𝑁= 2000withthewidelycoveredzone.

(10)

Table 3. The measured running timeandtheerrorofapproxima- tionin thecaseof movingmesh algorithmswithdifferentcontinu- ouslymodifiedgridfunctions.Ini- tialvaluesofthegridfunctionpa- rameters: 𝑁= 10000, 𝑦A= 0.191, 𝑦B= 0.203,𝑑= 0.003,𝑅= 0.8,Δ𝑦= 0.002, 𝛿= 0.01. In the method calledCMG1onlythedistancebe- tween𝑦B−𝑦A,𝑑,andΔ𝑦arede- creased.InthemethodcalledCMG 2𝑅isalsodecreased,furthermore inthemethodcalledCMG3even 𝛿isdecreased.InthecaseofNMG (non-modifiedgridalgorithm)the initialvaluesoftheparametersare the same, and arekept constant duringthesimulation.

Mesh type Time (min) error

CMG 1 233 7

CMG 2 235 14

CMG 3 257 17

NMG 283 11

3.3.3. Continuousmodiﬁcationoftheparameters

Duringthetransitionofthediodetheextensionofthereactionzone shrinks.Inthenewsteadystate,it’sextensionisonethird,thaninthe saltfreediode(section3.1).Isitpossibletoreachbetterperformance byshrinkingtheﬁnemeshzonewiththereactionzone?

Toalterthegridfunctioncontinuously,thepreviouslydiscussedal- gorithmismodified:ifthestopconditiondefinedinCOMSOLisfulfilled andtheiterationstops,thenthedependent variables(and itsspatial derivatives)mustberetrievedatthegridpoints.Thehalfwidthofthe localreaction ratepeakwas usedtoestimate thereactionzone’sextension,fromwhichnotonlythelocationofthefinemeshzone,but otherparametersofthegridfunctioncouldbere-calculated(forfurther detailssee5thsectionofSupplementaryInformation).Unfortunatelyif thetotalnumberofgridpointsistoosmall,theestimationofshrink- ageofthereactionzonebecomesunprecise,leadingtonon-converging mesh.Itmeans,thatinthecaseofsmall𝑁continuousmodificationof theparameterscanbenumericallyunstable.

The performance of algorithms with continuously modiﬁed grid (CMG)andanonmodiﬁedgrid(NMG)issummarized inTable3.The methodcalledCMG1,whereonly𝑦B−𝑦A,𝑑andΔ𝑦areshrinked,gives thesmallesterrorandsimulationtime.

4. Conclusionsandoutlook 4.1. Conclusions

We suggestan empiricalmoving meshalgorithm, which was at- tachedtothecommercialFEMsolverCOMSOL.Realizingthathigher reaction rateand higher concentrationusually appearclose to each other,thelocalreactionratewasusedtocontrolthemeshadaption.The proposed algorithmwas testedfor areaction-diffusion systemcalled acid-basediode,andtheparameterswerefine-tunedforthiscase.Dur- ingthetime-dependentmodelingoftheso-calledsalteffectsinthissys- temweexperiencedhighcomputationalcomplexity,thusanadaptive meshingalgorithmwasnecessarytoreducethetimeofthesimulation.

Wefound,thatusingthelocalreactionrateinsteadofthespatialderiva- tivesmakestheimplementationofthealgorithmveryeasyandsimple.

4.2. Outlook

Themovingmeshalgorithmcanbeextendedtohandletheappear- ance anddisappearance of thereactionzones.Ifthetime-dependent solverisstopped,oraftereveryiterationstep,themodeled domaincan

bescannedwithamonitorfunction,whichevaluatesthelocalreaction rates,andmarksthelocationofthereactionzone.Fromthelocationsa newgridfunctioncanbecreated,inwhichthelesssteepstraightseg- mentsofthepiecewisegridfunctioncorrespondtothereactionzones.

(Theexponentialsmoothing,andhandlingoftooclosereactionzones needtobedealtwith.)

Unfortunately,astraightforwardextensionoftheproposedprocess fortwo-dimensionalgeometricsetupseemstobeimpossible.Forsuch cases, thegrid-transformationbetween consecutivetime stepsis ob- tainedfrom thesolution ofso-called movingmeshPDE (Zegelinget al.,2005).Themeshhastoberegularizedinmanycasesandaninter- polationisalsonecessarybetweentwoconsecutivemeshes.Formany applications,however,one-dimensionalsituationsarealsoofgreatim- portance.

Furthermore,thelocalreactionratecanbeusedasa‘monitorfunc- tion’ not only in 𝑟-refinement, but in ℎ-refinement, 𝑝-refinement or ℎ𝑝-refinementmethodsaswell.Insuchcases,intheneighborhood of thechemicalreactionmeshelementscanbeadded(andlaterremoved), orthepolynomialdegreeof thebasisfunctioncanbe increased(and laterdecreased).Thesetechniquesareeasiertoextendtotwoormore dimensionalproblemssimulatedbyfiniteelementorfinitedifference methods.

Declarations

Authorcontributionstatement

ViktóriaKoncz&KristófKály-Kullai:Conceivedanddesignedthe analysis;Analyzedandinterpretedthedata;Wrotethepaper.Ferenc Izsák:Conceivedanddesignedtheanalysis;Wrote thepaper. Zoltán Noszticzius:Contributedanalysistoolsordata;Wrotethepaper.

Fundingstatement

ThisworkwassupportedbyOTKAGrantK131425andbytheELTE InstitutionalExcellenceProgram(1783-3/2018/FEKUTSRAT)through theHungarianMinistryofHumanCapacities.

Dataavailabilitystatement

Datawillbemadeavailableonrequest.

Declarationofinterestsstatement

Theauthorsdeclarenoconﬂictofinterest.

Additionalinformation

Supplementarycontentrelatedtothisarticlehasbeenpublishedon- lineathttps://doi.org/10.1016/j.heliyon.2020.e05842.

Acknowledgements

TheauthorsthankA.Patariczaforhelpfuldiscussions.

References

Bieniasz,L.K.,2000.Useofdynamicallyadaptivegridtechniquesforthesolutionofelec- trochemicalkineticequations:part5.Afinite-difference,adaptivespace/timegrid strategybasedonapatch-typelocaluniformspatialgridrefinement,forkineticmod- elsinone-dimensionalspacegeometry.J.Electroanal.Chem. 481(2),115–133.

Bieniasz,L.K.,2004a.Useofdynamicallyadaptivegridtechniquesforthesolutionof electrochemicalkineticequations.Part14:extensionofthepatch-adaptivestrategy totime-dependentmodelsinvolvingmigration-diﬀusiontransportinone-dimensional spacegeometry,andit.J.Electroanal.Chem. 565(2),251–271.

Bieniasz,L.K.,2004b.Useofdynamicallyadaptivegridtechniquesforthesolutionof electrochemicalkineticequations. Part15:patch-adaptivesimulationofexample transientexperimentsdescribedbyNernst–Planck–electroneutralityequationsinone- dimensionalspacegeometry.J.Electroanal.Chem. 565(2),273–285.

(11)

Cao,W.,Huang,W.,Russell,R.D.,2003.Approachesforgeneratingmovingadaptive meshes:locationversusvelocity.Appl.Numer.Math. 47(2),121–138.

Chun,H.,Chung,T.D.,2015.Iontronics.Annu.Rev.Anal.Chem. 8,441–462.

Cockburn,B.,1999.TheGraduateStudent’sGuidetoNumericalAnalysis’98.Springer, Berlin.

Coimbra,M.D.C.,Sereno,C.,Rodrigues,A.,2004.Movingﬁniteelementmethod:appli- cationstoscienceandengineeringproblems.Comput.Chem.Eng. 28(5),597–603.

Farmaga,I.,Shmigelskyi,P.,Spiewak,P.,Ciupinski,L.,2011.Evaluationofcomputational complexityofﬁniteelementanalysis.In:201111thInternationalConferencetheEx- perienceofDesigningandApplicationofCADSystemsinMicroelectronics(CADSM), pp. 213–214.

Heged˝us,L., Kirschner,N., Wittmann,M.,Simon, P.,Noszticzius,Z., Amemiya,T., Ohmori,T.,Yamaguchi,T.,1999.Nonlineareﬀectsofelectrolytediodesandtran- sistorsinapolymergelmedium.Chaos 9(2),283–297.http://www.ncbi.nlm.nih. gov/pubmed/12779826.

Heged˝us,L.,Wittmann,M.,Kirschner,N.,Noszticzius,Z.,1996.Reaction,diﬀusion,elec- tricconductionanddeterminationofﬁxedionsinahydrogel.Prog.Colloid&Polym.

Sci. 102,101–109.

Huang,W.,Ren,Y.,Russell,R.D.,1994.Movingmeshpartialdiﬀerentialequations(MM- PDES)basedontheequidistributionprinciple.J.Numer.Anal. 31(3),709–730.

Iván,K.,Kirschner,N.,Wittmann,M.,Simon,P.L.,Jakab,V.,Noszticzius,Z.,Merkin, J.H.,Scott,S.K.,2002.Directevidenceforﬁxedionicgroupsinthehydrogelofan electrolytediode.Phys.Chem.Chem.Phys. 4(8),1339–1347.

Iván,K.,Wittmann,M.,Simon,P.L.,Noszticzius,Z.,Vollmer,J.,2004.Electrolytediodes andhydrogels:determinationofconcentrationandpKvalueofﬁxedacidicgroupsin aweaklychargedhydrogel.Phys.Rev.E 70(61),061–402.

Kirschner,N.,Heged˝us,L.,Wittmann,M.,Noszticzius,Z.,1998.Nonlinear“salt-eﬀect”in anelectrolytediode.Theoryandexperiments.ACH- ModelsChem. 135,279–286.

Koncz,V.,Noszticzius,Z.,Roszol,L.,Kály-Kullai,K.,2017.Time-dependentmodelingof salt-contaminatedacid-basediodes.J.Electrochem.Soc. 164(4),H257–H265.http://

jes.ecsdl.org/lookup/doi/10.1149/2.1701704jes.

Lang,J.,1998.AdaptiveFEMforreaction—diﬀusionequations.Appl.Numer.Math. 26 (1–2),105–116.

LeVeque,R.J.,1992.NumericalMethodsforConservationLaws.Birkhäuser,Basel.

Lindner,J.,Šnita,D.,Marek,M.,2002.Modellingofionicsystemswithanarrowacid–

baseboundary.Phys.Chem.Chem.Phys. 4(8),1348–1354.

Marcus,Y.,1997.IonProperties.MarcelDekker,NewYork.

Merkin,J.H.,Simon,P.L.,Noszticzius,Z.,2000.Analysisoftheelectrolytediode.Electro- diﬀusionandchemicalreactionwithinahydrogelreactor.J.Math.Chem. 28,43–58.

Noszticzius,Z.,Schubert,A.,1973.“Electrolytediode”analysisofisothermaltransport processesintheinterfacesofaqueoussolutions.Period.Polytech.,Chem.Eng. 17, 165–177.

Pˇribyl,M.,Šnita,D.,Kubíˇcek,M.,2006.Adaptivemeshsimulationsofionicsystemsin microcapillariesbasedontheestimationoftransporttimes.Comput.Chem.Eng. 30 (4),674–685.

Roszol,L.,Várnai,A.,Lorántfy,B.,Noszticzius,Z.,Wittmann,M.,2010.Negativesalt eﬀectinanacid-basediode:simulationsandexperiments.J.Chem.Phys. 132(6), 064902.

Schubert,A.,Noszticzius,Z.,1977.Electrolytediode- anexperimentalstudypolarization phenomenaatthejunctionoftheaqueoussolutionsofanacidandabase.PartII.

Period.Polytech.,Mech.Eng. 21,279–283.

Schwab,C.,1999.p- andhp- FiniteElementMethods.ClarendonPress,Gloucestershire.

Trompert,R.A.,Verwer, J.G.,1991.A static-regriddingmethodfortwo-dimensional parabolicpartialdiﬀerentialequations.Appl.Numer.Math. 8(1),65–90.

Wang,R.,Keast,P.,Muir,P.,2004.Ahigh-orderglobalspatiallyadaptivecollocation methodfor1-DparabolicPDEs.Appl.Numer.Math. 50(2),239–260.

Zegeling,P.A.,Blom,J.G.,1992.Anevaluationofthegradient-weightedmoving-ﬁnite- elementmethodinonespacedimension.J.Comput.Phys. 103(2),422–441.

Zegeling,P.A.,deBoer,W.D.,Tang,H.Z.,2005.Robustandeﬃcientadaptivemoving meshsolutionofthe2-DEulerequations.Contemp.Math. 383,375–386.

Zegeling,P.A.,Keppens,R.,2001.AdaptiveMethodofLinesforMagnetohydrodynamic PDEModels.AdaptiveMethodofLines,pp. 117–137.

Zegeling,P.A., Kok,H.P.,2004. Adaptive movingmesh computationsfor reaction- diﬀusionsystems.J.Comput.Appl.Math. 168(1–2),519–528.

Zimmerman,W.B.J.,2006.MultiphysicsModelingwithFiniteElementMethods.World ScientiﬁcPublishingCompany,Singapore.