A least squares augmented immersed interface method for solving Navier-Stokes and Darcy coupling equations

ContentslistsavailableatScienceDirect

Computers and Fluids

journalhomepage:www.elsevier.com/locate/compfluid

A least squares augmented immersed interface method for solving Navier–Stokes and Darcy coupling equations

Zhilin Li

^a,∗

, Ming-Chih Lai

^b

, Xiaofei Peng

^c

, Zhiyue Zhang

^d

a Center for Research in Scientific Computation (CRSC) and Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA

b Department of Applied Mathematics, National Chiao Tung University, Taiwan

c School of Mathematical Sciences South China Normal University, Guangzhou 510631, China

d School of Mathematical Sciences, Nanjing Normal University, Nanjing 210023, China

a rt i c l e i nf o

Article history:

Received 6 December 2016 Revised 2 August 2017 Accepted 9 March 2018 Available online 13 March 2018 MSC:

65M06 65M85 76M20 Keywords:

Navier-Stokes and Darcy coupling Coupled fluid flow with porous media Least squares augmented IIM Fast Poisson/Helmholtz solver Analytic solution for Navier–Stokes and Darcy coupling

Interface relations

a b s t ra c t

Anew finitedifferencemethodbasedonCartesianmeshes andfastPoisson/Helmholtz solversispro- posedtosolvethecouplingofafluidflowmodeledbytheincompressibleNavier–Stokesequationsand aporousmediamodeledbytheDarcy’slaw.Thefinitedifferencediscretizationintimeisbasedonthe pressurePoissonequationformulation.Ateachtimestep,severalaugmentedvariablesalongtheinterface betweenthefluidflowandtheporousmediaareintroducedsothatthecoupledequationscanbede- coupledintoseveralPoisson/Helmholtzequationswiththoseaugmentedvariablesactingasjumpsofthe unknownsolutionand somedirectionalderivatives.Theaugmentedvariablesshouldbechosensothat theBeavers–Joseph–Saffman(BJS)orBeavers–Joseph(BJ)andotherinterfaceconditionsaresatisfied.It hasbeentestedthatadirectextensionoftheaugmentedideain[27]doesnotworkwellwhenthefluid flowismodeledbytheNavier–Stokesequations.Oneofthenewideasofthispaperistoenforcethedi- vergencefreeconditionattheinterfacefromthefluidside.Inthisway,theSchurcomplementmatrixfor theaugmentedvariablesisover-determinedandtheleastsquaressolutionisusedforthecouplingprob- lem.ThenewaugmentedapproachenablesustosolvetheNavier–StokesandDarcycouplingefficiently withsecondorderaccuratevelocityandpressureintheL^∞normfortestedproblems.Theproposednew ideainenforcingthe divergencefreecondition attheinterfacefromthefluidside hasalsobeenuti- lizedtosolve theStokesand Darcycouplingequations and showntooutperformthe originalmethod in[27].Inadditionaltothedetailedaccuracycheckforthepresentmethod,someinterestingnumerical simulationsforNavier–StokesandDarcycouplinghavebeenconductedinthispaperaswell.

© 2018ElsevierLtd.Allrightsreserved.

1. Introduction

Inthispaper,weintendtofindnumericalsolutionsforthecou- plingofafluidflowandporousmedia.Forself-containedpurpose, westatethegoverningequations,theboundaryandinterfacecon- ditions,andthephysical parametersin detailalthough they have appearedintheliterature.

The fluidflowismodeledbytheincompressibleNavier–Stokes equations,

ρ

f

 ∂

^uf

∂

^{t +uf·}

∇

^uf



+

∇

^pf=

μ

f



^uf+F,

∇

^{· u}f=0,

x∈



f, (1)

∗ Corresponding author.

E-mail addresses: zhilin@math.ncsu.edu , zhilin@ncsu.edu (Z. Li).

whereu_fisthefluidvelocity,p_fisthepressure,andF=(^F1,F₂)^is an externalbody force.The physicalparameters

μ

f>0and

ρ

f>0 are the viscosityanddensity ofthe fluid,respectively, which are assumed to be constant in f. The porous media is modeled by theDarcy’slawas

uD=− K

μ

^D

∇

^pD=−K1

∇

^pD,

∇

· uD=0,

x∈



D, (2)

where u_D and p_D are the corresponding fluid velocity and the pressure in the porous media, and K is the permeability which is assumedto be a constant in thispaper. Forsimplicity, we set K₁=K/

μ

D.Notethat,thenumericalmethodpresentedinthispa- percanbeextendedtoaddafluidsourcetermsuchas

∇

· uD=

φ

intheporousmediaregionwithoutmuchdifficulties.

Fig.1 illustrates thecoupling probleminwhicha porous me- dia is modeled by the Darcy’s law inside of a closed interface

 ^{while an} incompressible fluid flow is modeled by the Navier–

https://doi.org/10.1016/j.compfluid.2018.03.032 0045-7930/© 2018 Elsevier Ltd. All rights reserved.

Fig. 1. A diagram of a fluid flow and porous media. In the diagram, n and τare the unit normal and tangential directions, respectively.

Stokesequationsoutsideoftheinterface.Weassumethat theve- locity isprescribed along the outer boundary

∂

f (either no-slip or givenflow boundaryconditions). Our methodcan also handle the case where the fluid flow is inside while the Darcy’s lawis outside. In such a case, the velocity normal component u_D· n is usually prescribed along the outer boundary

∂

D.Notice that, it is equivalent to impose the normal derivative of the pressure as

∂^pD

∂ⁿ =−uD· n/K1.The well-posedness andregularity oftheabove couplingproblemhasbeeninvestigatedintheliteraturesuchasin [25]andthereferencestherein.

Therelationofvelocityandpressureacrosstheinterfaceofthe fluid flow and the porous media has been studied and is well- knownnowadays.The normalvelocity iscontinuous,which leads to

u_f· n=uD· n, (3)

whereagainweassume thatnistheunitnormaldirectionatthe interface pointing outward tothe fluid side. The stress actingon the interface shouldbe balanced which leadstwomore interface relations

[p]=2

μ

fn· Df· n, (4)

√

α

K

(

^uf− uD

)

·

τ

⁼²

τ

^{· D}f· n or _√

α

Ku_f·

τ

⁼²

τ

^{· D}f· n, (5) where D_f=¹₂



∇

^uf+

∇

^Tuf



is the deformation tensor,

τ

^{is the}

unit tangent vector, and ^√_α^K is the friction constant arising from the experimental data and dimensional analysis [25]. For conve- nience,wedenoteK₂=^√α_K.Thejumpinthepressureisdefinedas [p]=p_f− pD andthisdefinitionwillbeappliedtoothervariables throughoutthepaper.Thuswehave[u· n]=0and[

∇

· u]=0.The second equation in (5) is so called the Beavers–Joseph–Saffman (BJS) condition while the first one is called the Beavers–Joseph (BJ) condition. Most numericalmethods inthe literature use the BJS equation as a valid approximation while our augmented ap- proachcanhandlethemorecomplicatedtangentialslipcondition as shownin the first equation in (5)or non-homogeneouscases asshowninsome ofournumericalexamples.Intheremainingof thepaper,forthesakeofsimplicity,we omitthesubscriptsfand D in notationsofu_f, u_D,p_f andp_D ifthere isno confusionsince thepresentcoupledequationswillbesolvedasonewholesystem.

Notethat ifwe know[u·

τ

^{] and[u}· n] (which iszero accord- ingtothefirstinterface condition),wecangetthejumprelations ofthevelocity foreach component. At apoint X=(^X,Y)^{on the} interface,let

θ

^{be theanglebetweenthe normaldirectionofthe}

interface andthex-axis, see thediagram inFig. 1.Then theunit normaldirectionisn=(^cos

θ

,sin

θ

)^.Ifwedefine

[u·

τ

^]=q5, (6)

andcombinewith[u· n]=0,weget

0=[u]cos

θ

^+[

v

]sin

θ

^{, (7)}

q₅=−[u]sin

θ

^+[

v

]cos

θ

^{. (8)}

By solving the jumps of [u] and [v] from the above system, we obtain

[u]=−q5sin

θ

^{, [}

v

]=q₅cos

θ

^{. (9)}

There are many applications of the interactions between a fluid flow and a porous media, which has attracted a lot of atten- tionsintheliterature.Earlierworkshavebeenfocusedonsolving StokesandDarcycouplingequations,particularlyusingfinite ele- mentsordomaindecompositionmethods,seeforexample,[1–5,7–

10,12,14–17,19,20,23–25,32–36]andthereferencestherein.Forthe Navier–StokesandDarcycoupling,theoreticalanalysisandnumer- icalmethods,particularlyfiniteelementmethodshavebeencatch- ingup,seeforexample,[5,11,13,18]foran incompletelistandthe referencetherein.

In the above mentioned literature, almost all the numerical methodsare basedon finiteelementformulations ordomainde- compositionmethods.In[27],an alternativeapproachbasedon a finitedifferencediscretizationandaCartesianmeshwasproposed for solving the Stokes and Darcy coupling. The idea of the ap- proachistotransformtheoriginalcouplingproblemtothreePois- son equationswithasource anddipoledistributions correspond- ingtothejumps inthesolutionandthe normalderivative ofthe Poissonequations.Someaugmentedvariablesareintroducedalong theinterfacesothatthesolutionsofthevelocityandpressuresat- isfy originalinterface conditionsasshownin Eqs.(3)–(5). Oneof advantagesisthatafastPoissonsolvercanbeutilizedontherect- angulardomain.The methodproduces second orderaccurate ap- proximationforthevelocitybutfirstorderaccurateapproximation forthepressure.

Inthispaper,we propose aCartesian finitedifference method to solve Navier–Stokes and Darcy coupling equations based on some of ideas in [27]. We first rewrite the Navier–Stokes equa- tions in the fluid domain using the equivalent pressure Poisson equationformulationintroduced by JohnstonandLiu[21,22],and then alsorewrite the Darcy’s equations inthe porous media do- mainbytakingthedivergenceoperatortothevelocityfieldsothat the couplingequations can be written asone whole system, see the details in next section. We attempted to use the augmented approach in[27] combined witha Navier–Stokes solver to simu- latethe Navier–StokesandDarcycouplingproblems.The attempt workedonlyforafewtimestepsbeforethecomputedsolutionei- ther blows up or loses accuracy. After some careful analysisand numericalexperiments,wecame totheideatoenforcethediver- gencefreeconditionattheinterfacefromthefluidsideasanextra augmentedequation(constraint).ThustheresultantSchurcomple- mentmatrixfortheaugmentedvariablesindiscretizationisover- determinedandtheleastsquaressolutionisusedforthecoupling problem.Thatis whythenewmethodis calledthe leastsquares augmentedimmersedinterfacemethod.

Belowwebrieflysummarizethenoveltyofthispaper.

• WiththisnewleastsquaresaugmentedIIMapproach,wehave developed a finite difference scheme to efficiently solve the

Navier–StokesandDarcycouplingequationswithsecond order accuracybothinthevelocityandpressure.

• Withthenewmethod,wehaveobtainedsomeinterestingsim- ulation results for time dependent problems that cannot be solvedusingthemethodin[27]forStokesandDarcycoupling.

The method in this paper is different from that in [27] for StokesandDarcycouplingequationsnot onlyforsolvingdifferent couplings,butalsointhe methodologywiththenewidea inen- forcingthedivergenceconditionalong theinterface.Inthispaper, wealso presentsome comparison resultsofthe two approaches.

For stationary Stokes and Darcy coupling equations, both meth- odscanwork.However,asthemeshgetsfiner,thefive-equations methodin[27]maybelessaccuratecomparedwiththatobtained fromthenew(six-equations)approachastheaugmented variable may loose accuracy due to the larger condition number of the Schurcomplement.

The restofthe paperisorganizedasfollows.Inthe nextsec- tion,we explain the theoretical aspects ofthe leastsquares aug- mented method for the Navier–Stokes andDarcy coupling equa- tions. The mathematical equivalence of the original and trans- formed problems will also be discussed. In Section 3, we intro- ducethe time discretization forthetransformed equationsbased onavariantofthepressurePoissonequationformulationandde- rivethe resultant Schur complementmatrix. Numerical accuracy checkandthecomparisonswithpreviousmethodwillbegivenin Section2.2,whilesome interesting flow simulationswithvarious shapesandorientationsoftheinterfacewillbegiveninSection5. Conclusionsandacknowledgementsare giveninthelast twosec- tions.

2. TheleastsquaresaugmentedmethodforNavier–Stokesand darcycouplingequations

Inthissection,weexplainthetheoreticalbasisoftheproposed numericalmethod.Wetransformtheoriginalgoverningequations (1) and (2) to an equivalent form so that we get several Pois- son/Helmholtzequationsthatare easierto solverseparately.Note that,the theoreticaljustificationof thetransformforthe Navier–

Stokes part is based on the pressure Poisson equation formula- tiondiscussed in[21,22]. Byapplying the divergenceoperator to themomentum equation ofthe Navier–Stokesequations andthe Darcy’slaw,theentireproblemcanbere-writtenas



^p=

 ∇

· F−

ρ∇

·

(

^u·

∇

^u

)

, x∈



f,

0, x∈



D, (10)

[p]=q1, [pn]=q2, on



;



^u=

⎧ ⎨

⎩

1

μ



px− F1+

ρ

 ∂

^u

∂

^{t +u·}

∇

^u



, x∈



f,

−K1



^px, x∈



D,

(11)

[u]=−q5sin

θ

, [un]=q₃, on



;

 v

=

⎧ ⎨

⎩

1

μ



py− F2+

ρ

 ∂v

∂

^{t +u·}

∇v



, x∈



f,

−K1



^py, x∈



D,

(12)

[

v

]=q₅cos

θ

, [

v

n]=q₄, on



,

∇

· u=0, (13)

where pn=

∇

^p· n=^∂_∂n^p is the normal derivative of pand [pn]=

∂^pf

∂ⁿ −^∂_∂^p_n^D. Similar definitions are applied to [un] and [vn]. Here,

q₁,...,q₅areintroducedaugmentedvariablesalongtheinterface^. Theaugmentedvariablesarepartofunknownswithco-dimension onewhosesolutionssatisfytheinterfaceconditions(3)–(5).Inthe Darcy’sregion,wealsohavetheDarcy’slaw

uD=−K1

∂

^pD

∂

^{x ,}

v

D=−K1

∂

^pD

∂

^{y, (14)}

in addition to

∇

· uD=0. In this way, we decouple the original problemtothreeseparatedHelmholtz/Poissonequationsliterately and an explicit expression (14). Note that from above equations, we cannot conclude

∇

· uf=0 untilit istruealong the interface andtheboundaryfromthefluidside.Thatisthestartingpointof thenewmethod.

2.1. Augmentedequations

The above Navier–Stokes and Darcy governing equations are coupledbysixinterfaceconditionsalongtheinterface^,

[p]=2

μ

ⁿ· Df· n, (15)

K2



u_f− uD



·

τ

⁼²

τ

^{· D}f· n, or K2u_f·

τ

⁼²

τ

^{· D}f· n, (16)

∂

^pf

∂

^{n =}

μ 

^uf· n+



F−

ρ 

∂

^uf

∂

^{t +uf·}

∇

^uf



· n, (17)

∂

^pD

∂

^{n =−}

1

K1

(

^uD· n

)

, (18)

[u·

τ

^]=q₅, (19)

∇

^{· u}f =0. (20)

The first five equations are similar to those for the Stokes and Darcycouplingdescribedin[27].However, we imposeextraaug- mented equation (20) on the interface. Since there are only five augmented variables q₁, ..., q₅, in discretization, we would have more equationsthan unknowns. Thatis whythe newmethod is calledtheleastsquaresaugmentedmethod.

2.2. Equivalenceoftheoriginalandtransformedgoverningequations

The well-posedness of theStokes and Darcy;andthe Navier–

Stokes andDarcycouplinghasbeen well addressedinthelitera- ture,forexample,[25].Assumethattheoriginalmodelsare well- posedwithsuitableregularityineachsub-domain;thatis,theso- lutions existandareunique exceptforthepressure upto acon- stant. Then the solutions to the transformed problem Eqs. (10)– (14)withaugmentedEqs.(15)–(20)alsoexistsincetheyareallde- finedfromoriginalvariablesandsatisfythetransformedequations which are exactly derived fromthe original governing equations.

Inother words,thesolutions totheoriginal problemare alsoso- lutionstothetransformedproblemiftheregularityconditionsare satisfied in each sub-domain. More precisely, let u and v, andp bethesolutiontotheoriginalproblem.From thesequantities,we define q₁=[p], q₂=[pn],...,q₅=[u]·

τ

^{according to (10)–(13).}

Fromthe definitionsofuandv,p,andq₁,...,q₅,they satisfythe transformedEqs.(10)–(14)andaugmentedEqs.(15)–(20).Thusthe existenceofthesolutionstothetransformedsystemshasbeenes- tablished.

The uniqueness of the solution is more subtle. The issue has beenaddressedto someextentfortheStokesandDarcycoupling in[27].Thekeyiswhethertheincompressibility conditionissat- isfieduptotheboundary.Theargumentin[27]isthattheincom- pressibilityconditionalong theboundary/intefaceisredundantat leastnumericallyandtherefore canbe discarded. Theoretically,it

isstill anopen questionthatthevelocity ofthetransformedsys- temisdivergencefree.Inournewapproach,wesimplyenforcethe divergence freecondition along theinterface forthe transformed equations. Thus the solutionto the transformed system isalso a solution to the original problemsince all original equations now are satisfied. Thesensitivity ofnew formulationstill needs to be investigated.

3. ThealgorithmfortheNavier–StokesandDarcycoupling

In this section, we explain the method for the Navier–Stokes andDarcycoupling.TheNavier–Stokessolveristheonebasedon thepressurePoissonequationformulationdevelopedin[21,22,31]. Thealgorithmisdescribedindetailinthissection.

Weassumetheentirecomputationaldomain=f∪Disa rectangle[a,b]× [c,d].Weuseauniformgrid,

xi=a+ihx, i=0,1,...M, hx= b− a

M ; (21)

yj=c+jhy, j=0,1,...N, hy= d− c

N , (22)

where thevelocity andpressure are all definedatthe collocated gridpoints(x_i,y_j). Notethat,itispossibletousetheMAC(mark- and-cell) grid,whichusually providesbetter resultsanddoesnot require the boundary condition for the pressure. However, the presence of an arbitraryinterface wouldmake the programming muchmorecomplicatedwhenusingtheMACgrid.Inournumeri- cal experimentsshowninlatersections,wedonotencounterany stabilityissuewhenweusethecollocatedgridsincetheReynolds number considered here is modest. The augmented equationsat theinterfaceandtheboundaryalsomakenumericalschemesmore stable. The interface between the flow and porous media is dis- cretized bya setofcontrolpoints(^Xl,Y_l),l=1,2...N_b connected byacubicspline[28].Thejumpandinterface conditions,andthe augmentedvariables arealldefinedatthose controlpointsinthe discretization.Weuseatimemarchingschemetoobtainapproxi- matesolutionsatdifferenttimelevelt^kwithaknowninitialcon- dition at t⁰=0. Here the time step size is a constant for sim- plicity although we can use an adaptive one. From time step t^k tot^k+1=t^k+^t,we denotethediscretesolution (approximation to p, u,andv) attime step t^k as

{

^P_{i j}^k

}

,

{

^U_{i j}^k

}

,and

{

^V_{i j}^k

}

^{which to-}

gether form a vector u˜^k withdimension O(3MN), where M and N arethenumberofgridlinesinthex-andy-directions,respec- tively.We alsodenotethe discretesolution (approximationtoq₁, q₂,..., q₅) at time stept^k asQ₁^k_,l, Q₂^k_,l,..., Q₅^k_,l, l=1,2,...,N_b, whichtogethertoformavectorQ^k∈R^5Nb.

3.1. Timemarchingschemefromtimet^ktot^k+1.

Assumethat wehavealreadycomputedthevelocity (^Ui j^k,V_{i j}^k), thepressureP_{i j}^k,attimelevelk,weusetheCrank–Nicholsontype discretizationandtheaugmentedideatocomputethesolutionsat thenexttimelevelt^k+1.

Given an augmented vector Q^k+1, we solve the pressure in Eq.(10)usingtheIIMtoget

{

^Pi j^k+1

}

^{fromthefollowing,}



^pk+1=

 ∇

· F^k+1/2−

ρ∇

·

(

^u·

∇

^u

)

^{k+1/2, x∈}



f,

0 x∈



D, (23)

[p^k+1]=q^k₁⁺¹,

∂

^pk+1

∂

ⁿ



=q^k₂⁺¹, on



. (24)

Intheexpressionsabove,we useasecondorderAdams-Bashforth schemetoapproximatethenon-lineartermas

∇

^·

(

^u·

∇

^u

)

^k+1/2=3₂

∇

^·

(

^u·

∇

^u

)

^k−1₂

∇

^·

(

^u·

∇

^u

)

^{k−1. (25)}

Onceweobtainthepressure,wecansolvetheHelmholtzequa- tions(11)–(13)togetthepredictedvelocityfromthefollowing,



^u∗− 2

ρ μ 

^tu∗=

⎧ ⎪

⎪ ⎪

⎪ ⎪

⎪ ⎨

⎪ ⎪

⎪ ⎪

⎪ ⎪

⎩

2

μ

p^k_x⁺¹−

ρ



^tuk



−



^uk

+2

μ 

ρ (

^u·

∇

^u

)

^k+1/2− F₁^k+1/2



, x∈



f,

−K1



^pkx⁺¹+2K1

ρ

μ 

^{tpkx+1, x∈}



D, (26)

[u^∗]=−q^k₅⁺¹sin

θ

^,

∂

^u∗

∂

ⁿ



=q^k₃⁺¹, on



^;

 v

^∗− 2

ρ μ 

^t

v

^∗=

⎧ ⎪

⎪ ⎪

⎪ ⎪

⎪ ⎨

⎪ ⎪

⎪ ⎪

⎪ ⎪

⎩

2

μ

p^k_y⁺¹−

ρ



^t

v

^k



−

 v

^k

+2

μ

 ρ (

^u·

∇v )

^{k+1/2− F}2^k+1/2



, x∈



f,

−K1



^pky⁺¹+2K1

ρ

μ 

^{tpky+1, x∈}



D, (27)

[

v

^∗]=−q^k₅⁺¹cos

θ

^,

∂v

^∗

∂

ⁿ



=q^k₄⁺¹, on



;

Notethat inthe Darcy’sregion, we add theterm−_μ^2ρ_tu^∗ to the left so that the right hand side becomes ²_μ^K1ρ_t

∇

^{pk+1 (since u}=

−K1

∇

^{p). Thus, the} intermediate velocity can be solved by a fast Helmholtzsolverintheentirerectangulardomain¹

Finally,wecanproceed withthe projectionsteptoensure the divergencefreeconditionforthecomputedvelocity,

⎧ ⎪

⎪ ⎨

⎪ ⎪

⎩

 ψ

^k+1=

∇

^{· u∗}



^{t , x∈}



∂ψ

^k+1

∂

ⁿ

 

∂

=0,

 ψ

^k+1



=0,

∂ψ

^k+1

∂

ⁿ



=0, on



; (28)

u^k+1=u^∗−



^t

∇ ψ

^{k+1. (29)}

Althoughwehavetheupdatedvelocityinthewholedomain^, westillreplacethevelocitybytheDarcy’slawinD as

u^k_D⁺¹=−K1

∂

^pkD⁺¹

∂

^{x ,}

v

^k_D⁺¹=−K1

∂

^pkD⁺¹

∂

^{y , (30)}

sothattheDarcy’slawcanbe exactlysatisfied.Theprevious step isto make use ofa fast Helmholtzsolver forthe velocity in the fluidregion.Thediscretizationofthefirstandsecondorderderiva- tivesinvolvedabove isbased onthestandard centered five-point finitedifferenceformulas.

3.2.AugmentedvariablesandtheSchurcomplementmatrix

The solutions to the three Poisson/Helmholtz equations (23)– (27)withjumps inthe solutionandthenormal derivativein the matrix-vectorformcanbewrittenas

Au˜^k+1+BQ^k+1=F˜^k₁⁺¹ (31)

1 A direct application of IIM will lead to O (1) local truncation errors near the in- terface because of O (h/ t) = O (1) , and thus a first order accurate solution globally.

forsomevector F˜^k₁⁺¹ andsparsematricesAandB.

The next step is to evaluate the residual of the augmented equations(15)–(20).This step involveslocal interpolations andit isequivalenttodiscretizetheinterfaceconditions(15)–(20).

At each control point X_l, we interpolate the discrete solution

{

^Pi j^k+1

}

,

{

^Ui j^k+1

}

,and

{

^Vi j^k+1

}

^{togettheirvalues,normalandtangen-}

tialderivativesfromeach side ofthe interface,P_l^k+1,±, (^∂_∂^P_n)^k_l^+1,±, (^∂_∂τ^P)^k_l^{+1,± and so on. Note that, for a} Poisson/Helmholtz equa- tionwithknownjumpconditionsinthesolution andthe normal derivative, the computed solution and first order partial deriva- tives at grid points have been proved to be second order accu- rateusing theIIM [6]. From the second order accurate values at gridpoints, we alsogetsecond orderaccurate valuesatthe con- trolpointssincethecoefficientsoftheinterpolationsareO(1).

The discretization ofthe interface conditions(15)–(20)can be formallywrittenas

R

(

^Qk+1

)

=Cu˜^k+1+DQ^k+1− ˜F^k₂⁺¹ (32)

forsomesparsematricesCandD.Thesolution(u˜^k+1,Q^k+1)satis- fiesR(^Qk+1)=0 andEq.(31) sowe havethefollowingsystemof equations,

A B

C D



_u˜k+1

Q^k+1



=



F˜^k₁⁺¹ F˜^k₂⁺¹



. (33)

Therefore,theSchurcomplementforQ^k+1is

(

^D− CA⁻¹B

)

^Qk+1=˜F^k₂⁺¹− CA⁻¹F˜^k₁⁺¹=¯F^k+1, or

SQ^k+1=¯F^k+1. (34)

It has been shown in [29] and other related papers, that the matrix-vectormultiplicationSQ^k+1 givenQ^k+1 issimply

SQ^k+1=R

(

^Qk+1

)

+¯F^k+1=R

(

^Qk+1

)

− R

(

⁰

)

, (35)

whereR(^Qk+1)=SQ^k+1− ¯F^k+1.TherighthandsideofSQ^k+1=¯F^k+1 canbecomputedfrom−R(⁰)correspondingtotheresidualofthe interfaceconditionswithzerovalueofaugmentedvariableswhich resultsfromthethreeregularPoisson/Helmholtzequationsonthe rectangulardomain.Inournumericaltests,we formthematrixS, butnot A, B,C, andD, by using Q^k+1=e_i, i=1,2,...,5N_b. The coefficient matrix has a rectangular form with the dimension of 6N_b× 5N_b.Wethenusethesingularvaluedecomposition(SVD)to findtheleastsquares solutionofthe system. Weoutline thenu- mericalalgorithmfromtimestept^ktot^k+1asfollows.

Step1:FindtherighthandsideoftheSchurcomplement ¯F^k+1 bysettingQ^k+1=0,solving thetransformedproblem(10)–(14)to get u˜^k+1(⁰),andusing u˜^k+1(⁰)^andQk+1=0tointerpolateinter- face quantities needed in Eqs. (15)–(20). The residual of the in- terfaceconditionsistherighthandsideoftheSchurcomplement withanegativesign.

Step2:Fori=1,2,...,5N_b,wesetQ^k+1=e_i,thei-thunitvec- tor,solvethe transformedproblem(10)–(14)toget u˜^k+1(^Q),and theninterpolate u˜^k+1(^Q)^{toget theresidualoftheinterface con-} ditions(15)–(20).Thei-thcolumnoftheSchurcomplementisthe residualof theinterface conditions(15)–(19) plus theright hand side ¯F^k+1.Note that, fora fixed interface andconstant time step size,thisstepjustneedstobedoneonceinitially.

Step3:SolvetheSchurcomplementsystemtogetthesolution Q^k+1usingtheSVDdecomposition.

Step 4: Solve the transformed problem (10)–(14) with com- putedQ^k+1togetthepressureandvelocityu˜^k+1.

3.3. Someimplementationdetails

There are some major differences between the algorithm for solving the Navier–Stokes and Darcy coupling and that for the StokesandDarcycoupling.Weoutlineafewasbelow.

FortheStokesandDarcycoupling,weneedtosolvethreePois- sonequations.FortheNavier–StokesandDarcycoupling,weneed to solve a Poisson equation forthe pressure, but two Helmholtz equationsforthepredictedvelocityu^∗.Thecoefficient−_μ^2ρ_t isof 1/^t∼ 1/hintheNavier–Stokesequationsbutitshouldbezerofor theDarcy’sequation.ThisO(1/h)jumpinthecoefficient oftheu^∗ terminthe Helmholtzequations wouldlead toO(h) errorinthe solution,seeforexample,[26,29].However,thereisnosuchissue fortheStokesandDarcycoupling.Inpresentapproach,weaddthe term−_μ^2ρ_tu^∗ to both sidesof the Darcy’s equation. In thisway, wecan getridofthediscontinuityinthecoefficient,andstill are abletousethefastHelmholtzsolverintheentirerectangulardo- main.Thevalueofu^∗intheDarcy’sregionisjustaghostvalueto makethecomputationinthewholedomainmoreaccurateandef- ficient.Afteru^∗intheflowregionhasbeencomputed,wereplace thevalueinDarcy’sregionbyu^k+1=−K1

∇

^pk+1.

UsingthepressurePoissonequation formulation,howtocom- pute^u· nalongthe boundary

∂

^{andtheinterface} ^iscrucial to thestabilityofthe algorithm.A directinterpolation fromuto get^{uoftenleadsto}instability.Thereareafewdiscussionsinthe literature about how to compute ^u· n if the velocity is known along rectangularboundaries, see forexample, [21,30]. Itis quite challengingforcurvedboundaries/interfacesanddifferentbound- ary/interfaceconditions.Weusetherelation

μ 

^uf=

ρ

 ∂

^uf

∂

^{t +uf·}

∇

^uf



+

∇

^pf− Ff (36)

togetthevalueof^ufatgridpoints,theninterpolatethemtoget thevalues atthecontrol pointsneededfortheaugmented equa- tions.Notethattherighthandsideonlyinvolvesfirstorderderiva- tivesofthosefluidquantities.

Again, usingthe pressurePoissonequation formulation, ithas been discussed in the literature that it is not necessary to per- form the projection step. We have tested both approaches with andwithoutthe projection step.The methodwiththe projection givesbetterresultseventhoughnotsignificantly.Furthermore,the enforcement ofthe incompressibility condition wouldensure the equivalenceoftheoriginalandtransformedproblemsasdiscussed earlier.Notethat,forfixed interfacesandtimeindependentinter- faceconditions(15)–(20),theSchurcomplementmatrixSisafixed time-independentmatrixso theSVDdecompositionjustneeds to beperformedonce.

4. Numericalaccuracytestsandcomparisons

Inthissection,weshowsomenumericalaccuracytestsforthe presentalgorithm andthe comparisonwithprevious one in[27]. We use a uniformmesh withhx=hy=h=4/N for various N so thatwe canutilizeafastPoisson/Helmholtzsolver.The computa- tionaldomainis[−2, 2]× [−2, 2].Theinterfaceisexpressedusing the cubicspline package [28].Throughout thissection, theinter- face^{isaunitcirclex2}+y²=1andallphysicalparameterssuch as

ρ

^,

μ

^{, K}1,K₂ are all set to be one forsimplicity, unlessother- wise stated. The errors of thecomputed solutions are definedin L^∞normas

^Ep

^N_∞=max

i, j



p

(

^xi,yj,T

)

− P_{i j}^K



_, (37)

^Eu

^N_∞=max_{i, j}

 

u

(

^xi,y_j,T

)

− U_{i j}^K



2

+



v ₍

x_i,y_j,T

)

− V_{i j}^K



2

,(38)

Fig. 2. Velocity plot of a fluid flow and a porous media (inside the domain) modeled by the Navier–Stokes and Darcy coupling. (a) and (b), the interfaces are ellipses. (c) and (d), the interface is r = 0 . 5 + 0 . 1 sin 5 θin the polar coordinates. For smaller K ’s, the porous media acts as some sort of obstacles while for larger ones, the flow can go through easily. In all the simulations, the final time is T = 6 . 5 . The inflow condition is u = 1 , v = 0 . K = K 1 .

forameshsizeNbyN,whereKisthefinaltimestepcorrespond- ing tothe final timeT=K^{t.The orderofaccuracy isgivenap-} proximatelyby

order=log₂



^Ep

^N∞/

^Ep

²∞^N



(39)

for the pressure and similarly for the velocity. Unless otherwise stated,wechoosethenumberofcontrolpointsN_b=N.

4.1. Accuracycheckagainstanalyticsolutions

Tovalidateournumericalalgorithm,wefirstcomparethecom- putedsolutions againstanalyticsolutions.Thefirstexamplehasa

discontinuoustangentialslipforthevelocity fieldwhilethepres- sure is continuous along the interface ^{. The normal derivative} ofthe pressureis alsodiscontinuousacross the interface.On the other hand, the second example has a continuous tangentialve- locity but the pressure is discontinuous along the interface. The normalderivativesofthevelocitycomponentsarealsodiscontinu- ousacrosstheinterface. Thediscontinuitiesarevaryingalong the interfacein thesecond example.Noticethat, bothexamples have continuousnormal velocity along the interface since it is one of therequiredinterface conditionsasdiscussed earlier.Wedemon- strate that our present algorithm can handle both cases well so moregeneralexamplescanworkwithoutanydifficulty.

Fig. 3. Velocity plot of a porous media and a fluid flow (inside the domain) modeled by (a) the Darcy and Stokes coupling; and (b) the Darcy and Navier–Stokes coupling.

The parameters are μ= 0 . 2 , K 1 = 0 . 05 , K ² = 1 . 3 with the mesh N = 100 and N b = 80 . The interface is r = 0 . 5 + 0 . 1 sin 5 θin the polar coordinates; (c)–(d) simulation results for the Darcy and Navier–Stokes coupling with the interfaces, x ² / 0 . 4 ² + y ² / 0 . 6 ² = 1 and r = 0 . 5 + 0 . 2 sin 5 θ, respectively. The parameters are set as μ= 0 . 02 , K 1 = 0 . 5 , and K 2 = 1 . 3 . The initial data is obtained using the solution of the Darcy and Stokes coupling.

4.1.1. Example1

TheanalyticsolutionoftheNavier–Stokesequationsis u_f =g

(

^t

)

y

(

^{x2+y2− 1}

)

^+2y



,

v

f =g

(

^t

)

− x

(

^x2+y²− 1

)

− 2x



, p_f =g

(

^t

)

x²+y²



, (40)

definedoutsideoftheunitcirclex²+y²=1whiletheanalyticso- lutionfortheDarcy’ssystemisu_D=

v

D=0,andp_D=g(^t)^.Atthe interfacex²+y²=1,wehave

n=[x,y]^T,

τ

^{=[−y, x]T, p}f =g

(

^t

)

^{, p}D=g

(

^t

)

^,

u_f=2yg

(

^t

)

,

v

f=−2xg

(

^t

)

, uD=0,

v

D=0,

∂

^uf

∂

^{n =4yg}

⁽

^t

⁾

^,

∂v

f

∂

^{n =−4xg}

⁽

^t

⁾

^,

∂

^uf

∂τ

^=2xg

⁽

^t

⁾

^,

∂v

f

∂τ

^=2yg

⁽

^t

⁾

^,

n· Df· n=0

τ

^{· D}f· n=0. (41)

It iseasy tocheck that theabove solutions u_f andu_D satisfy the Navier–StokesandDarcy’sequations,respectively,withthesource term F in the Navier–Stokes equations being computed directly fromtheanalyticsolution.Theinterfacecondition(4),[p]=2

μ

ⁿ· D_f· nwith

μ

=1issatisfied.TheBJS condition(thesecond inter- face condition in Eq.(5)) is also satisfied. However, the solution

Fig. 4. Navier–Stokes and Darcy coupling with a sharp angled interface with different parameters at T = 3 . 5 . In (a)-(b), the fluid flow is outside and the Darcy’s law is inside of the interface; while in (c)-(d), the fluid flow is inside while the Darcy’s law is outside of the interface. The parameters are μ= 0 . 1 , K 1 = 0 . 02 , and K 2 = 10 ⁻⁴ in (a);

μ= 0 . 1 , K 1 = 0 . 02 , and K 2 = 10 ⁻³ in (b); μ= 0 . 02 , K 1 = 0 . 25 , and K 2 = 1 . 3 in (c); and μ= 0 . 02 , K 1 = 0 . 5 , and K 2 = 10 ⁻⁴ in (d). The corners do have significant influences on the flow patterns.

doesnotsatisfythefirstcondition(BJ)inEq.(5)since

τ

· Df· n=0 while −(^uf− uD)·

τ

=2g(^t)^. Nevertheless, the present algorithm candealwitheitherinterfaceconditionsinEq.(5).

In Table1,we show agrid refinementanalysisforg(^t)=sint atthefinal time T=5.Inthe table,thefirstcolumnisthemesh size inboth x-andy-directions;the second, fourthcolumnsare themaximumerrorsofthepressure andvelocity, whilethethird andfifthcolumnsarecorrespondingorderofaccuracy,respectively.

The sixthandseventh columnsare theconditionnumbersofthe Schur complement matrix definedas cond(^S)= ^max_min^iσi

iσi, where

σ

i

arenon-zerosingularvaluesoftheSchurcomplementmatrix,with

andwithouttheaugmentedequation

∇

· uf

|

=0alongtheinter- face,respectively.Onecanobservethattheconditionnumberwith theconstraintissignificantsmallerthanthatofwithoutit,which wouldalsoaffecttheaccuracyforthecomputedvelocity andthe pressure.InTable1(a),wehavethesamenumberofcontrolpoints ontheinterfaceasthemeshsize,thatis,N_b=N.Theaverageor- dersforthepressure andvelocity are 2.0231and3.2327,respec- tively. With the cubic spline representation of the interface, we canusefewercontrolpointswithoutaffectingtheaccuracyforthe pressure andthe velocity.For an interface with a modest curva- ture,ourcubicsplinepackage hasorderofaccuracyO(^h3s),where

Fig. 5. Darcy and Navier–Stokes coupling with an interface and various parameters and different locations of the interface. The Darcy’s law is defined outside of the interface.

(a), μ= 0 . 1 , K 1 = 0 . 5 , K 2 = 10 ⁻⁴ . (b), μ= 0 . 002 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . The interface is rotated counterclockwise by π/6. (c), μ= 0 . 002 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . The interface is rotated by π/3 counterclockwise. (d), μ= 0 . 002 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . The interface is rotated by 3 π/4 clockwise.

Table 1

A grid refinement analysis of Example 1 at time T = 5 . The parameters are μ= 1 , ρ= 1 , K 1 = 1 , K 2 = 1 . (a) the re- sults with N b = N. (b) the results with fewer control points N b < N . The CPU times are 0.5462, 1.5884, 5.9216, 34.272, 226.9478 s for (b), respectively.

(a) The average order of accuracy for the pressure and velocity are 2.0231 and 3.2327, respectively. ( N b = N)

N E p ^N∞ order E u ^N∞ order cond-6eq cond-5eq

16 1.8890 1 . 9839 e − 01 8 . 0806 e + 01 1 . 2104 e + 02

32 5 . 6152 e − 01 1.7502 2 . 1319 e − 02 3.2182 2 . 9443 e + 02 1 . 3904 e + 03 64 1 . 4805 e − 01 1.9232 1 . 7205 e − 03 3.6313 1 . 0435 e + 03 1 . 7704 e + 04 128 3 . 5373 e − 02 2.0654 1 . 4974 e − 04 3.5223 7 . 0632 e + 03 3 . 6746 e + 05 256 6 . 9209 e − 03 2.3536 2 . 5414 e − 05 2.5588 7 . 0651 e + 04 1 . 2520 e + 07 (b) The average order of accuracy for the pressure and velocity are 2.0221 and 3.2110, respectively. ( N b < N )

N N b E p ^N∞ order E u ^N∞ order cond-6eq cond-5eq

16 16 1 . 8890 e − 00 1 . 9839 e − 01 8 . 0806 e + 01 1 . 2104 e + 02 32 22 5 . 7256 e − 01 1.7221 2 . 2640 e − 02 3.1314 2 . 0234 e + 02 7 . 3228 e + 02 64 30 1 . 4790 e − 01 1.9528 1 . 7365 e − 03 3.7046 7 . 7313 e + 02 6 . 7254 e + 03 128 42 3 . 5359 e − 02 2.0645 1 . 5333 e − 04 3.5015 3 . 0480 e + 03 8 . 4892 e + 04 256 68 6 . 9394 e − 03 2.3492 2 . 6986 e − 05 2.5064 1 . 0844 e + 04 1 . 0888 e + 06

Fig. 6. Navier–Stokes and Darcy coupling with an interface and various parameters and different locations of the interface. The Darcy law is defined inside of the interface.

(a), μ= 0 . 1 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . (b), μ= 0 . 1 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . The interface is rotated by π/3 clockwise. (c), μ= 0 . 1 , K 1 = 0 . 02 , K 2 = 10 ⁻³ . The interface is rotated by π/6 counterclockwise. (d), μ= 0 . 05 , K 1 = 0 . 1 , K 2 = 10 ⁻³ . The interface is rotated by π/4 counterclockwise.

hs is the meshsize of the spline interpolation. Thus to maintain secondorderaccuracy,weneedtohaveroughlyO(^h3s)∼ O(^h2),or N_b∼ N^2/3. In our simulation, we use roughly N_b∼ 16+N^2/3. The number 16 is used to make sure that there are enough initial pointsontheinterface.InTable1(b),weperformthesameexper- iments but with fewer control pointsN_b which leads to smaller Schur complement matrix and obviously the condition number becomes slightly smaller. However, the average accuracy orders 2.0221 forthe pressureand 3.2110forthe velocity remain about thesame.

InTable.2,wealsoshowhowwelltheinterfaceconditionsare satisfied.Wecanseethatalltheresidualsarerelativelysmalland decrease quadraticallyforthisexample.Ifonlyfive equationsare

used,the residualsare then proportional tocond(S)



^,where



^is

themachineprecisionandcond(S)istheconditionnumberofthe Schurcomplement.

4.1.2. Example2

Inthepreviousexample,althoughthevelocityhasatangential slipalongthe interface butitis aconstant inthe Darcy’sregion.

Here,wepresentanexamplewithacontinuoustangentialvelocity butdiscontinuouspressure alongthe interface.More importantly, thevelocity and pressureare non-trivial inboth regions and the normalderivativesofthevelocitycomponentsarealsodiscontinu- ousacrosstheinterface.TheanalyticsolutionoftheNavier–Stokes