A coupled immersed boundary and immersed interface method for interfacial flows with soluble surfactant

ContentslistsavailableatScienceDirect

Computers and Fluids

journalhomepage:www.elsevier.com/locate/compfluid

A coupled immersed boundary and immersed interface method for interfacial flows with soluble surfactant

Wei-Fan Hu

^a,∗

, Ming-Chih Lai

^b

, Chaouqi Misbah

^c

a Department of Applied Mathematics, National Chung Hsing University, 145, Xingda Road, Taichung City 402, Taiwan

b Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

c Université Grenoble Alpes, LIPHY, Grenoble F-380 0 0, France

a rt i c l e i n f o

Article history:

Received 30 October 2017 Revised 25 January 2018 Accepted 6 April 2018 Available online 9 April 2018 Keywords:

Immersed interface method Moving irregular domain Soluble surfactant Immersed boundary method Navier–Stokes equations

a b s t r a c t

Inthispaper, wedevelop acoupledimmersedboundary(IB)and immersedinterfacemethod(IIM)to simulatetheinterfacialflowproblemswithsolublesurfactant.Thatis,acoupledsystemofsurface-bulk convection-diffusionequations mustbesolvednot onlyonamoving interfacebutalsoinanevolving irregulardomain.Basedontheimmersedinterfaceframework,wefirstbeginwiththenumericaldevel- opmentforthediffusionequationinafixedirregulardomainandthenextendtheschemetosolvethe convection-diffusionequationinanevolvingdomain.Thefluidmotiongovernedbytheincompressible Navier–Stokesequationsissolvedbyusingtraditionalimmersedboundarymethodwhilethebulksur- factantissolvedbytheproposedimmersedinterfacemethod.Aseriesofnumericaltestsforthepresent schemehavebeenconductedtoillustratetheaccuracyandapplicabilityofthemethod.Wefirstperform theaccuracyandefficiencytestsforthepresentIIMsolver.Wethenchecktheconvergenceofthesurface andbulksurfactantandthefluidvariablesintheinterfacialflowproblems.Wefurtherrun aseriesof numericalsimulationsforasuspendeddropletunder shearflowwithpresenceofsolublesurfactantto studytheeffectsofthedimensionlessBiotnumberandthebulkPecletnumberonthedropletdeforma- tionindetails.

© 2018ElsevierLtd.Allrightsreserved.

1. Introduction

Manyproblemsinbiological,physicalandmaterialsciencesin- volve solving partial differentialequations (PDEs)in complexdo- mains or deformable interfaces.It isknown that solving PDEs in complex domains or deformable interfaces numerically is quite challengingespeciallywhentheinterface(ortheinteriorboundary ofdomains)is moving.Evenin thecaseofonlysurfacematerial, developing numericalmethods for convection-diffusionequations onanevolvinginterfaceisstillofmajorinterestinscientificcom- putingcommunity.

The immersed interface method is a sharp interface method andhasbeensuccessfullyappliedtoproblemsarisingfromawide range of area. For instance, LeVeque and Li proposed the first IIM to solve elliptic interface problems with discontinuouscoef- ficients[15];further applicationssuchasfluid-structureproblems [11,28],ACdielectrophoresis[8],dropletelectrohydrodynamics[9], andvesicleelectrohydrodynamics[10,13],justtonameafew.Vari-

∗ Corresponding author.

E-mail addresses: wfhu@nchu.edu.tw (W.-F. Hu), mclai@math.nctu.edu.tw (M.-C.

Lai), chaouqi.misbah@univ-grenoble-alpes.fr (C. Misbah).

ousversionsofIIMhavealsobeenimplementedsuchasimmersed finite element method [5], particle-mesh methods [18], coupling immersedinterface andlevelsetmethod[29],etc.Seethereview inthebook[16].

In thispaper, we develop a coupled immersed boundary and immersedinterface methodtostudytheinterfacialflowproblems withsolublesurfactant. Thatis,thefluid systemgovernedby the incompressibleNavier–Stokes equationsis solved by using tradi- tionalIB method while theconvection-diffusion equation for the bulk surfactantis solved by a sharp immersed interface method.

Surfactant molecules typically consist of a hydrophilic head and a hydrophobic tail may adsorb and desorb between bulk fluids andthe interface so that the interfacial tensioncan be reduced.

Meanwhile, thisnon-uniformdistribution of surfactantmolecules along the interface produces an extra force (Marangoni force) in thetangentialdirectionto affectthe dynamics.The effectdue to insoluble surfactant has been extensively investigated since past two decades. For the insoluble surfactant case we only need to solve the surface concentration equation on an evolvingsurface, whileforthe solublecasewe haveto considera coupledsystem ofsurface-bulkequationsincomplexdomainsordeformableinter- faceswiththeinteractionthroughadsorptionanddesorptionpro- https://doi.org/10.1016/j.compfluid.2018.04.013

0045-7930/© 2018 Elsevier Ltd. All rights reserved.

Fig. 1. The computational domain consists of the physical domain ⁺ and the extended domain ⁻ with the embedded interface . The unit outward normal vector on is defined by n .

cesses.Therefore,tacklingsuchproblemvianumericalmethodsis indeedessentialandcanbeaquitechallengingissue.Thenumeri- calsimulationbyusingboundaryintegralmethodwasproposedby EggletonandStebe[4],they assumedasimplifiedmodelthat the bulkdiffusivity dominatesthe bulk surfactanttransportation (i.e., thebulk Pecletnumberisassumedto be zero)so thatthe effect oftheconstantbulksurfactantonlyappearsinthesourcetermof thesurfacesurfactant equation. BootyandSiegel [2] investigated the interfacial fluid flow at large values of Peclet number via a numericalmethodthatincorporatesasingularperturbationanaly- sisintoa full numericalsolution,whereas themethods proposed in[2,4]are limitedto theflow inStokes regime.Forthe Navier–

Stokesflow,dependingontheinterfacerepresentation,Muradoglu andTryggvason introduced a front-tracking method to studythe interfacialflows withsolublesurfactantin axisymmetric[20] and 3D multiphase flows [21]. Teigen et al. [26] proposed a diffuse- interfacemethodforsolvingtwo-phase flowswithsolublesurfac- tantinboth2Dand3Dspaces.ChenandLai[3]developedacon- servativeschemeforthecoupledsystemofsurface-bulkconcentra- tionequation. Theirextension work onthe dropletbouncing and coalescencewithsolublesurfactantcanbefoundin[23].Recently, Xu et al. [30] proposed a level set method for two-phase flows withsolublesurfactant.However,mostoftheaforementionedref- erencesemploy the numerical methods using a smoothing tech- niquesothatthebulkconcentrationisregularizednearthevicin- ity of interfaces which fail to capture the solution discontinuity across interfaces.Until recently, Khatri and Tornberg[12] used a two-dimensionalsegment projection methodto representthe in- terfacewhichiscapableofcapturingthejumpdiscontinuityofthe bulkconcentrationsolution.

In this paper, we take advantage of IIM in which a uniform Cartesiangridislayoutinaregular computationaldomainsothat thespatialdiscretizationsusingfinitedifferenceforderivative op- eratorscanbe performedin astraightforward manner exceptfor thegridpointsnearaninterface.Firstly,weproposeasimpleand efficientmethodtosolvethediffusionequationinanirregulardo- main(seethe setupinFig.1). Byextendingtheirregular domain toaregularone,theboundaryconditionofthesolutiononthein- terfaceisalternativelyregardedasthejumpconditionwhichisin- corporatedintheimmersedinterface discretization.Althoughthe idea of using IIM to solve the elliptic type of PDEs in irregular domains is not new as already described in [16], the treatment of jump conditions in the present scheme is different from the one in [16] where the augmented approach is adopted. We also provideboth implicitbackward EulerandCrank-Nicolson scheme andthen check theconvergenceofthemethod.Secondly,we ex- tendthe scheme forsolving the convection-diffusionequation in anevolvingdomain.Asimplesplittingmethodisadoptedinwhich we first update the intermediate solution by a convection step andthensolvethediffusionequation.Thenumericaldiscretization nearthe interface of the convection partcan be trickyand here we propose a simple explicit method which follows the idea of methodofcharacteristicsto trackthediscrete solutionaccurately viatheusage ofthejumpinformation.Againaconvergencestudy

forthe solveris performed.Lastly, we apply theimmersedinter- face solverto solvethe coupledsurface-bulk concentration equa- tionwithNavier–Stokesflow. Thefluid systemissolvedby tradi- tionalimmersedboundarymethod.Thiscoupledsolversharesthe samespiritwithourpreviousworksforthesimulationsofdroplet electrohydrodynamics[9]andvesicleelectrohydrodynamics[10].

The paper is organized as follows. An efficient immersed in- terface methodfor solving the diffusion equation in an irregular domainisdevelopedinSection2.Theextensionofthenumerical schemeforsolvingtheconvection-diffusionequation inamoving irregulardomainispresentedinSection3.Thenewnumericalal- gorithmthat couplesimmersedboundaryandimmersedinterface methodforsolvingtheinterfacialproblemswithpresenceofsolu- blesurfactantisoutlinedinSection4.Aseriesofnumericalsimu- lationsto investigatethestudyofhow theBiotnumberandbulk Pecletnumberaffectthe dropletdynamicsinshear floware pre- sented. Some concluding remarks and future works are given in Section5.

2. Asimpleversionofimmersedinterfacemethodforthe diffusionequationinanirregulardomain

Inthissection,weproposeanumericalmethodforsolvingthe diffusionequation inanirregulardomainbasedontheframework oftheimmersedinterfacemethod.ConsideradomainⁱⁿR²and asimpleclosedinterface^immersedin^{whichseparates}^into twodomains;namely,^{+(exteriortotheinterface)and}^−(inte- riorto theinterface)sothat =⁺∪⁻,seetheillustrationof thesedomainsinFig.1.Throughoutthispaper,wefocusonsolv- ingequationsintheexteriordomain^{+althoughourmethodcan} beappliedtosolvetheequationsintheinteriordomain^{− with-} out anydifficulty.Theinterface  ^isrepresented bya Lagrangian parametricformofX(^s,t)=(^X(^s,t),Y(^s,t)),0≤ s≤ 2

π

^,wheresis

theparameteroftheinitialconfigurationoftheinterface.Here,we assumethat theinterface ^{isnotevolving,andthusthedomain}

^+(or^{−)isfixed.}

Letusconsiderthefollowingdiffusionequationin^+as

∂φ ∂

^{t =}

1

Pe

 φ

ⁱⁿ



⁺,

φ|

∂=

φ

b, (1)

wherePeisthedimensionlessPecletnumber.Ofcourse,theabove equation shouldbe accompaniedwithaninitial condition

φ

^(x,0)

andoneinnerboundaryconditionon^{;namelytheDirichlet(}

φ

=

φ

D), or Neumann (^∂φ_∂n=

φ

N), orRobin type (^∂φ_∂n=

αφ

+g,

α

>0), wheren isthe unit outward normal vector pointingfrom^{− to}

^{+alongtheinterface}^{.Itisimportanttomentionthat,thetype} of boundarycondition on  ^{justresults ina slightdifference of} thepresentnumerical schemeaswe cansee fromthenumerical implementationlaterinnextsubsection.

The present immersedinterface method forsolving Eq.(1) is tosetthesolutionbeingidenticalzerointheinteriordomain⁻, andreplacetheinnerboundaryconditionon ^{by thejumpcon-} dition so that the equation in irregular domain ^{+ becomes an} equationinthewholeregulardomain^as

∂φ ∂

^{t =}

1

Pe

 φ

ⁱⁿ



,

φ|

∂=

φ

b, (2)



φ

=

φ

⁺

(

s,t

)

, 

φ

ⁿ=

φ

n⁺

(

s,t

)

^on



, (3) where

φ

n= ^∂φ_∂n standsfor the normal derivative andthe bracket

· denotes the jump discontinuity across the interface  ^{of the} quantityapproachingfromthe^+sidetothe^{−side.Oneshould} notice that,either one ofthe above jump conditionsis givenfor imposingDirichletorNeumannboundaryconditionon^;whereas forRobinboundarycondition,bothjumpsareunknownbutrelated

Fig. 2. The five-point Laplacian of the irregular point x i,j .

to eachother.Thesetwo jumpconditionswill beincorporatedin thenumerical discretizationforspatialderivatives.The remaining issue is how to solve the above diffusion Eq. (2) withthe jump discontinuities(3)in^{whichweaddressasfollows.}

2.1. Animmersedinterfaceframework

To proceed,let us firstlayout a uniform Cartesian grid inthe computational domain^{withmeshwidthh}=^x=^yforsim- plicity. The grid point x_{i, j}=(^xi,y_j) ^{is defined at the grid center} where the discrete solutions

φ

i, j=

φ

(^xi, j) ^{are located (here we} temporarilyomitthetimevariable).Theinterfaceisembeddedand cutsthroughsome gridcellssothesolutionisnot smoothacross theinterfaceaswecanseefromthejumpconditionsinEq.(3).We then classify thegrid point aseithera regular orirregularpoint.

Foraregularpoint,itmeansthatthestandardfive-pointLaplacian discretization(denotedbyh) doesnot cutthrough theinterface sothesecond-orderlocaltruncationerrorisachieved.Ontheother hand, atan irregular point, thefive-point Laplacian cutsthrough theinterface sothegridpointsusedinvolvebothinsideandout- sidetheinterface.Sincethesolutionanditsderivativeshavejumps acrosstheinterface,acorrectionfortheLaplaciandiscretizationis thus needed atthe irregularpoint to maintainthe desiredaccu- racy. Thus,thediscretizationofLaplaceoperatoratthegridpoint x_i,jcanbegenerallywrittenintheformof



h

φ

ⁱ, j+C_{i, j}

h², (4)

whereC_i,jisthecorrectiontermwhichisnonzeroonlyifthegrid pointisirregular.

Let us describe what the correction term is at the particular irregularpointasdepictedinFig.2.Whenweapplythefive-point Laplacianto x_i,j,thegridpoint x_{i−1, j} fallsindifferentsideofthe interface sothe correction of discretization comes only fromthe point x_{i−1, j}. Toderive the correctionterm, one needs to findthe orthogonal projection of x_{i−1, j} atthe interface (say X∗ in Fig.2), andthenapply theTaylor’sexpansion alongthe normaldirection atX∗.Thecorrectiontermthusbecomes

C_{i, j}=

φ

+d_

φ

ⁿ+d²

2

(



 φ

−

κ



φ

ⁿ−



^s

φ



) 



X_∗

, (5)

where thevalue disthe signed distancebetween thegrid point x_{i−1, j}andtheorthogonalprojectionX∗,

κ

^{isthelocalcurvatureof}

the interface, ands isthe surfaceLaplace (or Laplace-Beltrami)

operatordefinedby



^s=

|

^X1s

| ∂

∂

^s



1

|

^Xs

| ∂

∂

^s



. (6)

Here,the subscriptin Xs denotesthe partialderivative ofX with respecttos.Itcanbeseenthatthefive-pointLaplaciandiscretiza- tion in Eq.(4) gives the second-order accuracy at regular points butonlyfirst-orderaccuracyatirregularones. Weleavethecom- putationof theinterfacialderivativessuch aslocal curvatureand thesurfaceLaplaceoperatortothenextsubsection.Themorede- tailedderivationofthecorrectionterminEq.(5)canbe foundin [7,24].

2.2.Implementationdetails

Since we assume that the interface is closed, we can use the Fourierspectral discretizationto representtheconfiguration X(s).

We first choose a collectionof discrete Lagrangian markers X_k= X(^sk)^withequallydistributedinterfacialcoordinatess_k=k^s,k= 0,1,· · · ,M withs=2

π

/M.Theinterface thus canbe represented intruncatedFourierseries(invectorform)as

X

(

^s

)

=

M/2−1

=−M/2

Xe^is, (7)

whereXaretheFouriercoefficientsforX(s),andcanbecomputed veryefficientlyusingthe FastFourierTransform(FFT). Underthis Fourierrepresentation,thederivativeswithrespecttoscanalsobe computedquiteeasilybyusingthepseudospectralmethod[25].In thisfashion,allgeometricalquantitiessuchascurvature

κ

^andthe

surfaceLaplaceoperatorsinEq.(6)canbecomputedwithspec- tralaccuracy(inaddition,theunittangentvectorcanbecomputed by

τ

=(

τ

1,

τ

2)=Xs/

|

^Xs

|

^{and the} normal vectoris n=(

τ

2,−

τ

1)^).

Sinceall theinterfacialquantities aredefinedinLagrangianman- ner so the jumps of 

φ

, 

φ

ⁿ, 

φ

, and s

φ

 are all defined atX_k.Thanks to theFourierrepresentation,those jumps andthe curvatureattheorthogonal projectionpointsX∗ inthecorrection terminEq.(5)canbeeasilyinterpolatedthroughthevaluesatthe LagrangianmarkersX_k.

Nowwearereadytodiscretizethedifferenceequationwiththe jumpconditionsinEqs. (2)and(3)basedonthe aforementioned immersedinterface framework. Here we usethe Backward Euler (BE)andCrank–Nicolson(CN)methodforthetimeintegration.De- noting^, 0, 1,and 2 asthesolutionvectorsformedby

φ

i,j,



φ

X_k,

φ

X_k, and

φ

nX_k,respectively, the BackwardEuler dis- cretizationforEq.(2)isgivenby



ⁿ⁺¹−



ⁿ



^{t =} 1 Pe

 

h



ⁿ⁺¹+C0

0ⁿ⁺¹+C1

1ⁿ⁺¹+C2

2ⁿ⁺¹



, (8)

where the superscript n stands for the time step tn=n^{t with} thetime step size^t;C0, C₁ andC₂ are theformal matrix oper- atorsinvolvingtheFourierinterpolationofthosejumpsinthecor- rectiontermasdescribedpreviously.Letusdescribehowtosolve theabovematrixequationinthefollowing.SupposethattheNeu- mannboundaryconditionisimposedon^{,ityieldsthat} 2ⁿ⁺¹ is givenbut ₁^{n+1 istreatedasanunknowntobe}determined.Inad- dition,unlikeour previous worksonsolving Laplaceequation for theelectricpotential (

φ

=0) forthedroplet andvesicleelec- trohydrodynamics in [9,10], the term 0ⁿ⁺¹(=

φ

ⁿ⁺¹X_k) ^{is also} anunknown here.We canapproximatethistermusingtheorigi- naldifferentialequationbyapplyingonthebothsidesofinterface with Backward Euler discretization to obtain ( 1ⁿ⁺¹− 1ⁿ)/^t= ₀ⁿ⁺¹/Pe.It isinteresting toseethat, althoughtheapproximation for ₀^{n+1 givesonlyfirst-orderaccuracyintime, thelocaltrunca-} tionerrorforthefive-pointLaplacianoperatorremains first-order accuracyinspaceatirregularpointsaslongaswetake^t=O(^h)

Fig. 3. A diagram showing least squares approximation for one-sided interpolation or normal derivative at X k .

(onecanseefromEqs.(4)to(5)).Proceedingtorearrangethedis- cretizationinEq.(8)givesthematrixequation



h− Pe



^tI



ⁿ⁺¹+

_Pe



^tC0+C1

1ⁿ⁺¹

=−Pe



^t



ⁿ+ Pe



^tC0

1ⁿ− C2

2ⁿ⁺¹. (9) Ontheotherhand,fortheCrank-NicolsondiscretizationofEq.(2), wehave ( 1ⁿ⁺¹− 1ⁿ)/^t=( 0ⁿ+ 0ⁿ⁺¹)/2Pe.Thus, theresultant matrixequationbecomes



h−2Pe



^tI



ⁿ⁺¹⁺

_2Pe



^tC0+C1

1ⁿ⁺¹

=−2Pe



^t



ⁿ−



h



ⁿ+C0

_2Pe



^t

1ⁿ+

0ⁿ

− C2

2ⁿ⁺¹. (10)

One should notice that, the computation of h^{n at irregular} pointsonthe righthand side ofEq. (10)indeed involvesthe cor- rectionterm at previous step but we use the same notation for simplicity. At this stage, it becomes clear that we have to add onemorematrixequationtoapproximatetheaugmentedvariable ₁^{n+1. However, aswe can seefrom Eq.(3),a one-sidedapprox-} imationfor

φ

⁺(^Xk,t_n+1) ^{isneeded. To}approximate that, we first constructaleastsquarescubic polynomialby usingthe valuesof

φ

iⁿ, j⁺¹ at the gridpoints satisfying



^xi, j− Xk



≤ 4hin the domain

^{+ as depictedin Fig. 3. In general, the number of chosen grid} pointsis larger than the number of unknown coefficients inthe leastsquaresapproximationsothecubicpolynomialP(x,y)canbe obtained.FortheinterfacepointX_kinahighcurvatureregion,we haveto increase thespatialresolution tomake surethat enough grid points fall into the same region so the least squares cubic polynomial is achievable. Once we have the polynomial P(x, y), theone-sidedinterpolation

φ

⁺(^Xk,t_n+1)^{canbecomputeddirectly.}

This least squares cubic polynomial approach roughly has third- order accuracy to the approximation of the function itself. With thisapproximation,we canwritethejumpcondition

φ

=

φ

^{+ in}

Eq.(3)inamatrixformby

B



ⁿ⁺¹⁼

1ⁿ⁺¹ (11)

whereB denotes the formal matrix arising fromthe above least squarescubicpolynomialapproximation.

Now,we coupleeitherthematrixEq.(9)(BE)orEq.(10)(CN) withEq.(11)intoonelinearsystemas



H C

B −I

 

_n+1

₁ⁿ⁺¹



=



F 0



. (12)

Inpractice,wedonotformthematricesCandBexplicitlyaswe can see fromthe iterativeprocedure of ourmatrix solver below.

To solve the above linear system, we first eliminate ^{n+1 from} Eq. (12) by Schur complement techniqueto obtain a new linear systemfor 1ⁿ⁺¹as



BH⁻¹C+I



1ⁿ⁺¹=BH⁻¹F.

Thisisa M× Msystemfor 1ⁿ⁺¹ whichismuch smallerthanthe originaloneinEq.(12).WethenusetheGMRESiterativemethod tosolvetheabovelinearsystem.SincetheGMRESmethodonlyre- quiresthematrix-vectormultiplication,itisnotnecessarytocon- structthematricesH⁻¹,C,andBexplicitly. Notethat,thematrix Hcomesfromthe five-pointLaplaciandiscretizationso theinver- sionofHcanbeperformedefficientlybyapplyingthefastPoisson solverprovidedbyFishpackpublicsoftwarepackage[1].

In summary, the detailednumerical algorithm for solving the linearsystemEq.(12)tofind^n+1and 1ⁿ⁺¹canbesplitintothe followingthreesteps.

Step1. ApplyonefastPoissonsolvertoobtain^∗in H



^∗=F.

Step2. UseGMRESiterationtosolve 1ⁿ⁺¹in



BH⁻¹C+I



1ⁿ⁺¹=B



^∗.

Step3. ApplyonefastPoissonsolvertosolveⁿ⁺¹ⁱⁿ H



ⁿ⁺¹=F− C

1ⁿ⁺¹.

DuringtheGMRESiterationinStep2,wesetthestoppingcri- terion as h². The overall computational cost for Step 1-3 in our presentscheme canbe evaluated in termsof thenumber offast Poissonsolverbeingapplied.Inthenextsubsection,weshallpro- videa moresystematicnumericalcheckontheaccuracyandeffi- ciencyforourpresentnumericalalgorithm.

We should remark that, for the case ofimposing Dirichletor Robinboundaryconditionontheinterface^{,thereisaslightdif-} ferenceoftheBEandCNdiscretizations fromEqs.(9)to(10).For the Robin boundary condition case, one still needs to apply the one-sidedinterpolationfor

φ

⁺(^Xk),whereasfortheDirichletcase the approximation for

φ

n⁺(^Xk) ^{is required and it can be simply} done by takingnormalderivative oftheleastsquares polynomial

∇

^P(Xk)· n(X_k)(thusthederivativeapproximatelyhassecond-order accuracy).Theresultant linearsystemfromothertypesofbound- aryconditionson^canberepresentedandsolvedsimilarlyasin Eq.(12).

2.3. Numericalaccuracyandefficiencystudy

Inthissubsection,we demonstratetheaccuracyandefficiency test of the present IIM scheme for solving the diffusion equa- tion developed in previous subsection. We construct an analyt- ical solution for Eq. (1) with a three-leaved interface as X(^s)= (^r(^s)^coss,r(^s)^sins)^{withthe radiusfunction r}(^s)=0.1(⁴+cos3s) as

φ (

x,y,t

)

=1+0.5exp



−2

π

^2t

Pe



cos

( π

^x

)

^cos

( π

^y

)

, x∈



⁺.

(13) Here, we test three cases by imposing different inner boundary conditions asDirichlet,Neumann, andRobintype boundary con- dition on ^{. The} computational domain is set by =[−1,1]× [−1,1]and thePeclet numberis chosen asPe=10. We setN to bethegridsize andthusthemeshsizeish=2/N. Thetimestep fortheBEandCNdiscretizationofEq.(9)or(10)ischosenasthe same of mesh width as^t=h and the simulation is run up to

Table 1

Mesh refinement results for accuracy and efficiency for the BE (upper half-panel) and CN (lower half-panel) discretization. The numerical solution is denoted by φh and exact solution by φe .

N Dirichlet Neumann Robin

φh −φe∞ Rate Iter. φh −φe∞ Rate Iter. φh −φe∞ Rate Iter.

128 7.627E − 04 – 3.25 9.051E − 04 – 2.10 7.935E − 04 – 2.04 256 3.789E − 04 1.01 3.42 4.489E − 04 1.01 2.10 3.902E − 04 1.02 2.02 512 1.889E − 04 1.00 3.50 2.233E − 04 1.01 2.09 1.935E − 04 1.01 2.01 1024 9.423E − 05 1.00 3.75 1.113E − 04 1.00 2.23 9.633E − 05 1.01 2.19 128 7.121E − 05 – 3.10 3.825E − 05 – 2.17 2.369E − 05 – 2.04 256 1.720E − 05 2.05 3.47 5.571E − 06 2.78 2.17 5.958E − 06 1.99 2.02 512 4.648E − 06 1.89 3.16 1.395E − 06 2.00 2.18 1.490E − 06 2.00 2.01 1024 1.145E − 06 2.02 3.43 3.460E − 07 2.01 2.24 3.722E − 07 2.00 2.06

Fig. 4. (a) The snapshots of the numerical solution at T = 0, 0.25, 0.5, 0.75 and 1. The numerical solution is subject to Dirichlet boundary condition on the outer boundary

∂and Robin boundary condition on the inner boundary . (b) The snapshots of the numerical error at T = 0.25, 0.5, 0.75 and 1.

the terminal time T=1. We presentthe rateof convergence for the maximum error between the numerical solution

φ

h and the exact solution

φ

^{e in Table1. As expected, the first-and second-}

order convergence are achieved for BE and CN scheme, respec- tively.As forthe computational complexity,one can seethat the average number ofGMRES iterationsin the time interval [0,1] is justlessthan4andbecomessteadyevenasthegridnumberdou- bles. Moreover, we plotthe snapshots ofnumericalsolution (ob- tained by CNdiscretization withRobinboundaryconditionon  with N=512) and snapshots of their errors

φ

^e−

φ

h at different timesinFig.4.Onecanseethatthepresentmethodcapturesthe jump discontinuityacrossthe interfaceaccurately andthelargest numericalerrormainlyresultsfromthecomputationaldomain⁺ rather than the extended domain ^{−. In fact, for each case we} tested here, the numerical errors at grid points in the extended domain^{−areconsistentwiththespatial}discretizationofmagni- tudeO(h²).

3. Animmersedinterfacemethodfortheconvection-diffusion equationinamovingirregulardomain

WithsuccessfuldevelopmentofIIMfordiffusionequationina fixed irregular domain,in thissection, we extendour methodto the following convection-diffusionequation ina moving irregular

domain⁺(^t)^as

∂φ ∂

^{t +u·}

∇φ

⁼_P¹_e

 φ

ⁱⁿ



^+,

φ|

∂=

φ

b, (14)

where u is a given velocity or is obtained fromsolving Navier–

Stokes equations (see in next section). As before, besides of the boundarycondition on

∂

^{,one of innerboundary conditions on}

^{(t)(Dirichlet(}

φ

=

φ

D),Neumann(^∂φ_∂n=

φ

N),orRobintype(^∂φ_∂n =

αφ

+g,

α

>0))must be given. Again, assuming thezero solution in the interior domain ⁻(^t) ^{and regarding the inner boundary} condition as jump conditions across the interface ^{(t), we can} rewrite the convection-diffusionEq.(14) in the fixed regular do- main=⁺(^t)∪⁻(^t)^{(referringthesetup ofthosedomainsin} Fig.1)as

∂φ ∂

^{t +u·}

∇φ

⁼_P¹_e

 φ

ⁱⁿ



,

φ|

∂=

φ

b, (15)



φ

=

φ

⁺

(

^s,t

)

^, 

φ

ⁿ=

φ

n⁺

(

^s,t

)

^on

(

^t

)

^{. (16)}

Certainly,Eq.(15) issubjectto an initialconditionin ⁺(⁰)^{. We} remindthereadersthattheabovejumpconditionsarerequiredfor thecalculationofthecorrectiontermswhereaseitheroneorboth ofthemcanbeunknown(dependingonthetypeoftheboundary conditionon^).

Withpresenceoftheappliedvelocityfiledu,theinterface ^is assumedtomovewiththelocalvelocityby

∂

^X

∂

^t

⁽

^s,t

⁾

^=u

⁽

^X

⁽

^s,t

⁾

^,t

⁾

^{. (17)}

Unlikeinprevioussectionwheretheirregulardomain^+isfixed, heredue to the fact that the interface ^{(t) is moving according} to Eq.(17) (thus ⁺(^t) ^and ⁻(^t) ^{are both evolving, but}  ^re- mainsasa fixed domain),we haveto carefullytrackthediscrete solutionsin theneighborhood of theinterface. Therefore,instead ofdiscretizingEqs.(15)and(16)directly,here,we proposeasim- pleoperatorsplittingscheme.Thatis,wefirstsolvetheintermedi- atesolutionfortheconvectionstep,andthen updatethesolution bysolving the diffusion part.The fullnumerical algorithm isde- scribedinthefollowingsubsections.

3.1.Anexplicitimmersedinterfacemethodfortheconvection equation

Inthissubsection,wefocusondevelopinganumericalscheme basedonIIMforsolving theintermediate solutionoftheconvec- tionequation

∂φ ∂

^{t +u·}

∇φ

=0 in



. (18)

Since the interface ^{(t) moves with the velocity field u, captur-} ingthejumpdiscontinuitybehaviorofthediscretesolutionsinthe vicinityof^{(t)becomesthemajortask.Eq.(18)is}discretizedby thesecond-orderLax–Wendroff scheme(9-pointstencil)[14]as

φ

iⁿ, j⁺¹=

φ

ⁿi, j

+ 1

2

α

ⁿi, j^+1/2

(

¹+

α

ⁿi−1/2, j

) φ

iⁿ−1, j−1

2

α

ⁿi, j^+1/2

( α

iⁿ−1/2, j+

α

iⁿ+1/2, j

) φ

ⁿ_{i, j}

+ 1

2

α

ⁿ_{i, j}^+1/2

(

−1+

α

ⁿi+1/2, j

) φ

iⁿ+1, j+1

2

β

_{i, j}^n+1/2

(

¹+

β

iⁿ, j−1/2

) φ

iⁿ, j−1

− 1

2

β

iⁿ, j^+1/2

( β

iⁿ, j−1/2+

β

iⁿ, j+1/2

) φ

iⁿ, j+1

2

β

iⁿ, j^+1/2

(

−1+

β

iⁿ, j+1/2

) φ

iⁿ, j+1

+ 1 8

α

ⁿi, j^+1/2

( β

iⁿ+1/2, j−

β

iⁿ−1/2, j

) φ

iⁿ, j+1−

( β

iⁿ+1/2, j−

β

iⁿ−1/2, j

) φ

ⁿi, j−1



+ 1 8

β

iⁿ, j^+1/2

( α

_{i, j+1}ⁿ _/2−

α

_{i, j−1}ⁿ _/2

) φ

ⁿi+1, j−

( α

_{i, j+1}ⁿ _/2−

α

_{i, j−1}ⁿ _/2

) φ

iⁿ−1, j



+ 1 8

 α

ⁿ_{i, j}^+1/2

β

iⁿ+1/2, j+

α

iⁿ, j+1/2

β

_{i, j}^n+1/2

 φ

iⁿ+1, j+1

− 1 8

 α

ⁿi, j^+1/2

β

iⁿ+1/2, j+

α

iⁿ, j−1/2

β

iⁿ, j^+1/2

 φ

iⁿ+1, j−1

− 1 8

 α

ⁿi, j^+1/2

β

iⁿ−1/2, j+

α

iⁿ, j+1/2

β

iⁿ, j^+1/2

 φ

iⁿ−1, j+1

+ 1 8

 α

ⁿ_{i, j}^+1/2

β

iⁿ−1/2, j+

α

_{i, j−1}ⁿ _/2

β

_{i, j}^n+1/2



φ

iⁿ−1, j−1, (19)

where

α

ⁿi, j=uⁿ_{i, j}^t/^{x and}

β

iⁿ, j=

v

ⁿ_{i, j}^t/^{y. We assume that all} thecoefficients onthe righthand side of theabove discretization are given fora prescribed velocity field. Notice that the velocity componentu andv appearing in those coefficientsare not only evaluated at the cell center x_i,j but also the cell-normal edges x_{i± 1/2,j± 1/2}.

As before,thegridpoint isidentifiedaseitheraregular orir- regularpoint.Foraregularpoint,wemeanthatallthegridpoints usedontherighthand side ofEq.(19)fallsintothesameside of the interface;whereas an irregular point involves using the grid pointsfrombothinside andoutsidethe interface. Itis clearthat weneedto modifyEq.(19)atirregularpointstoachieve thede- siredaccuracy.

Weillustratehowtofindthecorrectiontermviathetwocases showninFig.5.As afirstcasedepictedinFig.5(a),theinterface movesform ^(tn) to (^tn+1)^{. For a certain small time step} ^t,

Eq.(19)canberegardedasaninterpolationforsomegridvaluely- ingintheregion[x_i−1,x_i+1]× [y_j−1,y_j+1]attheprevioustimestep tn,andthisgrid pointflows along somecharacteristic linetox_i,j attime stept_n+1. Thisis exactly theconcept of methodof char- acteristics.Since thegridpoint x_i,jstays insidetheinterface, one needs tocorrect thegrid valueswhichcome fromtheoutsideof theinterface(⁺(^tn)^{).Forinstance,}

φ

iⁿ+1, j−1correct

===

φ

iⁿ+1, j−1+C_iⁿ_{+1, j−1},

in which the above correction term follows the same computa- tion asin Eq. (5). Of course, one also has to correct other grid values

φ

iⁿ+1, j,

φ

iⁿ, j+1, and

φ

iⁿ+1, j+1 as depicted in Fig. 5(a). For the second case shown in Fig. 5(b), as the interface moves, the grid point x_i,j fallsinto thedifferent sidesofthe interface (from

⁺(^tn)^to⁻(^tn+1)^{).Inthiscase,thegridvalue}

φ

_{i, j}^{n+1 comesfrom}

some grid value lyinginside the interface attime steptn.There- fore,one shouldcorrectthose gridvalueswhichoriginatingfrom

⁺(^tn)^side(say

φ

ⁿ_{i+1, j−1},

φ

_{i, j}ⁿ,

φ

_iⁿ_{+1, j},

φ

_iⁿ_{−1, j+1},

φ

_{i, j+1}ⁿ ,and

φ

_iⁿ_{+1, j+1}

inFig.5(b)).

From thesetwocases,we canconcludea ruleforthecalcula- tionofcorrectionterms:ifthegridpointx_i,jremainsonthesame side of theinterface when the interface is moving, one needs to correct thegrid value(s) which ison the different side of x_i,j at theprevioustimesteptn.Ontheotherhand,ifthegridpointx_i,j crosses the interface at differenttime steps, one have to correct thegridvalue(s)whichliesinthedifferentsideofx_i,jatt_n+1.As aconsequence,theLax–Wendroff schemeinEq.(19)canbegener- allymodifiedbyamatrixformof



^n+1=A



^n+C0

0ⁿ+C₁

1ⁿ+C₂

2ⁿ,

where A denotesthe coefficient matrix constructed by theright- handside ofEq.(19)andC₀ 0ⁿ+C₁ 1ⁿ+C₂ 2ⁿ arethe correction terms.OneshouldnoticethatthecorrectionmatrixC₀,C₁,andC₂ aredifferentfromthosematricesC₀,C₁,andC₂inprevioussection since theycontain correction termsmultiplying thosecoefficients inEq.(19).Asone cansee fromthe aboveequation, thesolution isexplicitlyupdated bythe solutionatprevioustime step,which canbedoneveryefficiently.

3.2. Fullnumericalalgorithm

Once we know how to update the intermediate solution for the convection part, now we are ready to solve Eq. (15) with Eq.(16)viathesplittingtechniqueasfollows.Giventhe interface position,thesolutionandjumpsatthetimestepn^t,Xn,n, 0ⁿ, ₁ⁿ,and ₂ⁿ,weupdate^{n+1byusingthefollowingtwosteps:}

Step1. AdvancetheinterfacetothepositionXⁿ⁺¹andthenapply themodifiedLax–Wendroff schemetoobtaintheinterme- diatesolution^{∗duetotheconvectionpart}



^∗=A



ⁿ+C0

0ⁿ+C1

1ⁿ+C2

2ⁿ.

Step2. Updatethesolution^{n+1bysolvingthediffusionpartus-} ingthebackwardEulerscheme



ⁿ⁺¹−



^∗



^{t =} 1 Pe

 

h



ⁿ⁺¹+C0

0ⁿ⁺¹+C1

1ⁿ⁺¹+C2

2ⁿ⁺¹



andupdatethediscretejumpconditions ₀ⁿ⁺¹, ₁ⁿ⁺¹,and 2ⁿ⁺¹.

InStep1,wenotonlyobtaintheintermediate solution^∗but alsothejumpconditionfor^∗(say ₁^{∗)ontheinterface} (^tn+1) as 1^∗= ₁^{n.Thisisdueto thefactthat solutiononthe interface} alsoflows along some characteristic lineso that 

φ

^∗=

φ

⁺(^s,tn)^. InStep2, ₀^n+1isthusapproximatedbyusing( ⁿ₁⁺¹− ₁^∗)/^t=

Fig. 5. Two cases for the interface moves from ( t n ) (left figure) to (t n+1) (right figure).

0ⁿ⁺¹/Pe.Thus, solvingfor^{n+1 isnothingbutfollowingthe dif-} fusion equation solver developed in previous section. The above splittingprocedureisexpectedtohavefirst-orderconvergencedue to the backward Eulermethod in Step2. Ofcourse, one can en- hancetheoverallaccuracyofthepresentmethodby employinga higherorderschemesuchasStrangsplitting[12],however,inthis paperweaimtoimplementthenumericalsolverinasimpleman- ner.Next,we shallcarry outanumericaltestfortheconvergence andefficiencyofourproposednumericalmethod.

3.3. Numericalaccuracyandefficiencystudy

In this subsection, we perform a numerical test for the con- vergence andefficiencyofthepresentmethodforEq.(14)inthe exterior of a moving interface inthe computational domain = [−2,2]× [−2,2].Here,wechoosetheinitialinterfaceasaunitcir- cle centered atthe originand apply the shear flow u=(^y,0) ^as our prescribed velocity so theevolving interface ^{(t) can be ex-} plicitly described as X(^s,t)=(^cos(^s)+tsin(^s),sin(^s)) ^{by solving} Eq.(17).Inthismoving irregulardomain,wecanconstructanan- alyticalsolutionforEq.(14)as

φ (

^x,y,t

)

^=1+0.5exp



−

(

^t+t³/3

)

Pe



sin

(

^{x− yt}

)

^{, x∈}



^+. (20)

Again,threecasesaretestedbyimposingdifferentinnerboundary conditionsasDirichlet,Neumann,andRobintypeboundarycondi- tionon ^{.In thesecases,the solutionis subjectto theDirichlet} boundaryconditionon

∂

^.

We choosethe Pecletnumber Pe=1 withdifferentgrid sizes N=128,256,512, 1024.andthetimestepsize^t=h/2.Table2 showsthe maximumerrorsforthenumericalsolutionattheter- minaltime T=1. One can seethat thepresent numericalresult performsclean first-orderconvergenceforthenumerical solution

φ

^{in each case. Theaverage numberof GMRESiterationsin Step}

2isjustwithintwoevenasthegridnumberdoubles.Wefurther showthesnapshotsforthenumericalsolution(obtainedbyusing innerRobinboundaryconditionon ^{withthegridsizeN}=512) atdifferenttimesinFig.6.

4. Interfacialflowswithsolublesurfactant

Inthissection, weapply theIIMsolver developedinprevious section tosimulateinterfacialflow problemswithsolublesurfac- tant.Consider an incompressibleNavier–Stokesflow consistingof two-phase fluids in a fixed rectangular domain  ^{with an im-} mersed droplet interface  represented by the parametric form X(s,t) with the parameter s∈[0, 2

π

^{] (see Fig. 1). It is assumed}

that the surfactantexists on theinterface asa monolayerand is adsorbed fromor desorbed to the bulk fluid in ^{+; that is, the} surfactant is soluble in the exterior bulk but not in the interior droplet. The interface is contaminated by the surfactant so that