An immersed boundary method based on domain decomposition

Delft University of Technology

An immersed boundary method based on domain decomposition

Krishnan, Navaneetha; Viré, Axelle; Schmehl, Roland; van Bussel, Gerard

DOI

10.1016/j.compfluid.2020.104500

Publication date

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Krishnan, N., Viré, A., Schmehl, R., & van Bussel, G. (2020). An immersed boundary method based on

domain decomposition. Computers and Fluids, 202, [104500].

https://doi.org/10.1016/j.compfluid.2020.104500

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Computers and Fluids 202 (2020) 104500

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

An

immersed

boundary

method

based

on

domain

decomposition

Navaneetha

Krishnan,

Axelle

Viré

∗

,

Roland

Schmehl,

Gerard

van

Bussel

Wind Energy Section, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS, Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 21 February 2019 Revised 6 January 2020 Accepted 6 March 2020 Available online 12 March 2020 Keywords:

Immersed boundary method Finite element method Fluid-structure interactions

a

b

s

t

r

a

c

t

A novelimmersed boundarymethod basedona domaindecompositionapproach isproposed in the contextofafiniteelementdiscretisationmethod.Itisapplicabletoincompressibleflowspastrigid, de-forming,ormovingbodies.Inthismethod,unlikemostimmersedboundarymethods,strongboundary conditionsareimposedintheregionsofthecomputationaldomainthatareoccupiedbythestructure.In ordertoachievethis,theproposedformulationdecomposesthecomputationaldomainbysplittingthe finiteelementtestfunctionsintosolidandfluidparts.InthecontinuousGalerkinformulation,this pro-ducesasmearedrepresentationofthefluid-structureinterface.Theabsenceofanimmersedboundary forcingtermimpliesthatthemethoditselfhasnoinfluenceontheCFLstabilitycriterion.Furthermore, thestiffnessmatrixinthemomentumequationissparserthancomparedwithotherforcingimmersed boundarymethods, and symmetryand positive-definitenessof theLaplacian operator inthe pressure equationispreserved.Asshowninthispaper,stabilityandaccurateimpositionofboundaryconditions makethemethodpromisingforhighReynoldsnumberflows.Themethodisappliedtothesimulations oftwo-dimensionallaminar flowoverstationary and moving cylinders,as well as amoderately high Reynoldsnumberflowpastanaerofoil.Goodresultsareobtainedwhencomparedwiththosefrom pre-viousexperimentalandnumericalstudies.

© 2020ElsevierLtd.Allrightsreserved.

1. Introduction

Fluid-structureinteractions (FSI)are encountered inmany en-gineering applications, such as the flow past large wind turbine blades, the deployment of spacecraft parachutes, or blood flows througharteries.Immersedboundarymethodsareattractiveto nu-merically simulatetheseproblems,especially whenthegeometry ofthestructureiscomplexandarbitrarily deforming[1–3].These methods immerse the structure ina fluid domain andmimic its effectonthefluiddynamics.Thiscontrastswithbody-conforming methods, in which the fluid domain excludes the structure and the governing equations of the flow field are solved only in re-gions surroundingthestructure. Theoriginal immersedboundary method was introduced by Peskin [4] to simulate cardiac flows. Sincethen,severalvariationsofthisapproachhavebeenproposed, asdetailedhereafter.

Based on how the fluid-structure interface is represented on the discrete mesh, immersed boundary methods can be classified as sharp-interface and diffused-interface methods. In diffused-interfacemethods,thediscreterepresentationofthe fluid-structure interfacespreadsover oneormoregridcells.Therefore,

∗ _{Corresponding author.}

E-mail address: a.c.vire@tudelft.nl (A. Viré).

theinterfaceissmearedanddiffersfromtheactualsurface bound-aryof the structure. Methods mimicking the effectof the struc-turethroughan additionalforcing terminthefluid’sequationsof motionoftenfall underthiscategory. Thesimplest oftheforcing methods are the penalty methods [5–7]. In that case, the force termthat is addedto the right-hand-sideof the fluid’s momen-tumequation isproportional to thedifference betweenfluid and structurevelocitieswithacertainpenaltyfactor.Thelatter,which isdifficulttointerpretphysically,canchangethestrengthatwhich the penalty condition is enforced. Such penalised forcing meth-odsarealsopopularinthecontextoffiniteelementmethods [8– 11].The penaltyforce isalways maskedbyasteppingfunction to limit theregion in whichit acts. The stepping function can be a Diracfunction across the interface element [6,8], weighted func-tions spreadacross multipleelements toavoid numerical oscilla-tions[7]oraconservativelyprojectedsolidconcentrationfieldthat marksthepresenceofthebody[11].Aslightlymodifiedversionof the penalty method is the feedback control method [12]. In this method, in addition to the damping term, the forcing term also has a stiffness term. From a physical point of view, the bound-aryconditionisthenimposedasadampedoscillator[1].A differ-encebetweenthe penalty forcing methodand thefeedback con-trolmethod lies in the waythe boundary conditionis satisfied: it isalways lagging behindin the penalty forcing method, while itisoscillatory inthefeedbackcontrolmethod.Anothercommon

https://doi.org/10.1016/j.compfluid.2020.104500

forcing method is the direct forcing method,in which the addi-tionalforcingtermbalancesouttheconvective,viscous,and pres-suregradienttermsintheNavier–Stokesequations[13,14].Inthis classof immersed boundarymethods, the physicalinterpretation oftheforcingtermisclearerthaninthepenaltymethod.A draw-back,however, isthat the implicitevaluationofthe forcing term leadsto adenser stiffnessmatrixinthediscrete systemof equa-tions.Inordertosolvethis,variationsoftheoriginaldirectforcing methodhavebeenproposedandevaluatetheforcing termat dif-ferenttimelevels.Forexample,asemi-implicitstaggeredapproach wasproposed andcomputedan intermediate velocity field with-out the forcing term, then calculated the forcing term usingthe intermediatevelocityfield,andeventuallymadeavelocity correc-tionusingtheforcingterm[15,16].Implicitapproachesareusually preferredinterms ofaccuracy,since theboundaryconditionsare imposedwhile solving themomentum equations.However, inall their variants, theno-slip conditionat the interface could be vi-olated inthe staggered stepto projectvelocity onto a solenoidal field.Also, solutionto thePoisson’s systemthat hastobe solved forthe projection is affected by the quality ofthe grid, both in-sideandoutsidetheimmersedboundary.Thiscanbe particularly challengingforproblemsinvolvingthinmembranes.

Oneoftheundesirablefeaturesofdiffused-interfacemethodsis thatthediscretefluid-structureinterfaceissmearedacrossseveral meshnodes.Cut-cellmethods areatypeofsharp-interface meth-ods that eliminate this issue by locally reconstructing the mesh toalign it withthe boundary [17–20].However, this approach is mostlyapplied totwo-dimensional problems asthe mesh recon-structioninthree-dimensionalcasesisinherentlydifficult.Evenin two-dimensionalproblems,itwasshownthatamajorchallengeis toproperly define andmanage thetopology of cut cells[3]. An-other sharp-interface approach is the hybrid Cartesian immersed boundarymethod,whichuses aCartesiangridandapplies strong boundary conditions at the nodes located in the structure [21]. Thismethodcanbemoreaccurate inapplyingthefluid-structure interfaceconditionsthan theprevious immersedboundary meth-ods.However, the useof aCartesian grid isinefficient toresolve boundarylayers nearcomplexgeometries, hencelimitingthe ap-plicability of the method to flows past simple geometries. The curvilinearimmersedboundarymethod[22]isanextensionofthis methodandusescurvilineargrids.Thelatterarebetterthan Carte-siangrids to resolve boundarylayers for relativelysimple curved geometries.

Finally,acompletelydifferentapproach,relevantforbothsharp anddiffusedinterfaces,isbasedontheideathat theforcingterm can be viewed as a Lagrange multiplier that imposes a velocity boundary condition in regions occupied by the immersed body

[23–25]. The role of the multiplier is then similar to therole of thepressurefieldinanincompressibleflow,wherepressureis ef-fectivelyusedtosatisfythedivergence-freecondition.Thismethod presentsseveraladvantagesasitisabletoimposeboundary con-ditions accurately, while preserving the symmetry of the Pois-son’ssystem.Also,theoperatorscanbetweakedtoimpose sharp-interfaceordiffused-interfaceconditions.

In this paper, a novel method is proposed to accurately im-posetheboundaryconditionatthefluid-structureinterface, with-outlimitationstoCartesiangridsandwithoutrequiringadditional equationsto be solved. The method is developed in the context ofthefinite-elementdiscretisation andapplies theno-slip condi-tionattheinterface asa strongboundaryconditioninthefluid’s momentumequations.Thisisachievedbydecomposingthe com-putationaldomainintotwo distinctregions,one forthefluidand anotherone for thestructure. The domain decompositionis per-formedbysplittingthefiniteelementtestfunctionsintotwoparts. The momentum equation is solved foracross the whole compu-tationaldomain witha strongimpositionofthe boundary

condi-Fig. 1. Schematic representation of the computational domain ( ), and solid ( s ) and fluid ( f ) sub-domains. _{∂is the boundary of the computational domain. ∂}f and ∂s are the boundaries of the fluid and solid sub-domains. The blue curve is the immersed boundary that is shared between the solid and fluid sub-domains. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tioninside thesolid phase.ThePoisson’sequation forpressureis solvedonthecomputationalsub-domainthatbelongstothefluid phase.Thepaperisorganisedasfollows.Section2outlinesthe ba-sicprincipleofthe methodandpresentsthegoverningequations in the solid and fluid sub-domains. Section 3 converts the con-tinuousequationsintodiscreteforms.Section 4demonstratesthe performanceofthemethodonseveraltestcases.First,three lam-inarcasesareconsidered:Taylor-Couetteflow,flowpasta station-ary two-dimensional cylinder, and flow past an oscillating cylin-der. Second, flow past a NACA0012 aerofoil at a Reynolds num-ber of Re=1000 issimulated to evaluate the stabilityand accu-racyoftheproposedmethodatamoderatelyhighReynolds num-ber.Results arecomparedagainstanalytical solutions, experimen-taldataandothernumericalsimulations.Inbothflowregimes, re-sultsshow good convergencerates andexcellent agreementwith existingstudies.

2. Continuousformulation

2.1. Strongformofthegoverningequations

As stated above, the computational domain in immersed boundaryproblemsencompassesbothfluidandsolidregions.This is illustrated by Fig. 1, which shows how the computational do-main



isdecomposedintotwosub-domains:afluiddomain(



f₎ andasoliddomain (



s_),with

_

₌

_

f_∪

_

s_and

_

f_∩

_

s_=0.The boundary of the computational domain is denoted by

∂



, while boundariesof the fluid andsolid sub-domainsare

∂



f _and

_∂

_

s_, respectively.Fig.1alsoshowsthattheimmersedboundaryofthe solid domain (

∂



s _{shown by the blue curve)is a sub-set of the} boundaryofthefluiddomain,i.e.

∂



f₌

∂

_

_∪

∂

_

s_.Theboundary

∂



canbefurtherdecomposedinto

∂



D_,

_∂

_

N_,and

_∂

_

R _toapply Dirichlet,Neumann,andRobinboundaryconditions,respectively.

Inthefluiddomain,theadvection-diffusionequation that gov-ernsthetimeevolutionofascalarfieldcisgiveninthe conserva-tiveformby,

 f



∂

c

∂

t +

∇

·

(

uc

)

−

∇

·

(

κ∇

¯¯ c

)

− F



dV=0, (1)

inwhichuisthetransportvelocityvector,

κ

¯¯ isthediffusivity ten-sor,andFrepresentsanyforcingterms.Thescalarfieldcisdefined

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 3 onthewholecomputationaldomainandcanbe,forexample,one

ofthevelocitycomponents,aspeciesconcentration,the tempera-ture, theturbulentkinetic energyk,theturbulentfrequency

ω

,or theturbulent dissipation



.Thedensity

ρ

ofthefluid isassumed tobeconstant,sothatthecontinuityequationforvelocityisgiven by,



f

(

∇

· u

)

dV =0. (2)

In the solid domain



s_{, the transport variable c is set to the} structuralvaluecs_{,whichgenerallyspeakingiseitherprescribedor} obtainedfromtheresolutionofthestructuraldynamicsequations, suchthat  s

∂

c

∂

tdV=  s

∂

cs

∂

t dV. (3)

In thispaper,results willbe shownfor fixedstructures or struc-turesundergoingaprescribedmotion.

Eqs.(1) –(3)arethe governingequationsofthe problemwith thefollowingboundaryconditionsinspaceandtime.

c

(

·,0

)

=c0

(

·

)

in



f, (4a) c

(

·,0

)

=cs 0

(

·

)

in



s, (4b) c=c˜ in

∂



D_{, (4c)}

κ∇

c· n=q in

∂



N_{, (4d)} −

κ∇

c· n=h

(

c− ca

)

in

∂



R, (4e)

where,c˜isthefixedconditionappliedon

∂



D_{,qisthefluxacross}

∂



N_{,histhesurfacetransfercoefficienton}

∂

_

R_,andc

aisan am-bientfieldvalue.

2.2. Weakformofthegoverningequations

The weak form of Eq. (1)is obtainedby multiplying its left-hand-sidewitha testfunction

ψ

f _{andusingGreen’s firstidentity} tointegratebypartsoverthedomain



f_,suchthat,

 f



ψ

f

∂

c

∂

t −

∇ψ

f_{· uc+}

_∇ψ

f_{· ¯¯}

_κ

_·

_∇

_c−

_ψ

f_F



dV +  ∂f



ψ

f

₍

_{uc− ¯¯}

κ∇

_c

₎

_{· n}



_{dS=0, (5)}

wherenisthesurfacenormalvector.Notethat

∂



f_includesboth the boundariesof thecomputational domain(



) andthe bound-arybetweenthefluiddomain(



f_{)andthesoliddomain(}

_

s_). Ad-ditionally,theintegraloverthepartsof

∂



f_{thatcoincidewiththe} boundaries ofthe computational domain (



) is often subdivided into various components depending on the boundary conditions thatareapplied.Forsimplicityitisassumedthatthepartsof

∂



f coincidingwiththeboundaries

∂



areclosed(withu· n=0)and that

∇

c· n=0 (Neumann condition). Furthermore,the test func-tion(

ψ

f_{)isconstructedsuchthatitsvalueiszeroatthefluid-solid} interface(bluecurveinFig.1).Undertheseconditions,the bound-aryintegralinEq.(5)cancelsout,whichyields

 f



ψ

f

∂

c

∂

t −

∇ψ

f_{· uc+}

_∇ψ

f_·

_κ

_·

_∇

_c−

_ψ

f_F



dV=0. (6)

Eq. (6) holds in the region



f _{covered by the fluid. In order to} write the equations of motionof the transport variable c in the wholecomputationaldomain



=



f_∪

_

s_{,theweakformsofthe} fluid- and solid-governing equations need to be added together.

TheweakformofEq.(3)isobtainedbytestingitwithatest func-tion

ψ

s_{.AddingthistoEq.(6)yields}

 



(

ψ

f₊

_ψ

s

₎

∂

c

∂

t −

∇ψ

f_{· uc+}

_∇ψ

f_{· ¯¯}

_κ

_·

_∇

_c −

ψ

f_F



dV=  



ψ

s

∂

cs

∂

t



dV. (7)

TheweakformofEq.(2)isobtainedbymultiplyingthe equa-tionwithanothertestfunction

ξ

f_,as,

 f

ξ

f

₍

_∇

_{· u}

₎

_{dV =0. (8)}

NotethatEq.(8)isnotintegratedbypartsandcontinuityisonly evaluatedin



f_{.Thischoicewillbereasonedinthefollowing} sec-tion.Eqs. (8)and(7) formtheweak formulation associatedwith

Eqs.(1)–(3).

3. Discreteformulation

3.1. Generalformulation

Inthispaper,continuouspiecewise-linearfiniteelementshape functions(

φ

i) are usedto approximatethe continuousflow vari-ablesonthediscretiseddomain,sothat

c≈ N  i=1 ci

φ

i, u≈ N  i=1 ui

φ

i, p≈ N  i=1 pi

φ

i, (9)

inwhichN isthenumberofnodesinthediscretiseddomainand pisthefluidpressure.

Several key points have to be noted before defining the test functions forthe problem. First, although



f_∩

_

s_{=0 in a} con-tinuoussense, inthepresentdiscrete formulationasindicated in

Fig.2, in which the shaded interface element is part of both



f and



s_,

_

f_∩

_

s_{= 0.Becauseofthisoverlap,thetestfunctions}

_ψ

f and

ψ

s_{usedinthetransportequationhavetobecomplementsof} eachothertoavoiddoublingthecommonterm,

∂

c/

∂

t.Second,the pressuregradient operatorwillbe testedusing

ψ

f_,andthus lim-ited to



f _{only. Thus,to makethe Laplacianoperator symmetric,} evaluationofcontinuity,andthereby

ξ

f_{, isalsolimitedto}

_

f_. Fi-nally, unlike

ψ

f_{,the function}

_ξ

f _{doesnot havea complementary} solid part.Thus, to fullyevaluate continuity in



f_,

_ξ

f _{should not} beforcedtozeroatthesolidendoftheinterface.Forthisreason,

Eq.(8) is not integratedby partssince thesurface integral term wouldnotvanishbasedonthisdefinitionof

ξ

f_.

The test functions (

ψ

i,

ξ

i) are chosen from a function space thatisappropriatefortheproblem[26].ForacontinuousGalerkin approach this coincides with the function space from which the shapefunctionsarederived.Wedefine

ψ

_if,

ψ

s

i,and

ξ

f i as,

ψ

f i =

φ

i

(

1−

α

in+1

)

,

ψ

s i =

φ

i

α

in+1, (10a)

∃

x∈

α

n+1 e :x<1⇒

ξ

f k =

φ

k. (10b)

Here,

α

(blue line inFig.2) is a solidconcentration field that representsthepositionofthestructureonthefluidmesh,andthe superscript n+1 denotes that the solid position fromthe forth-comingtime stepisconsidered. Itisobtainedby aconsistent in-terpolationoraboundedconservativeprojectionofaunitaryfield fromthesolidmeshontothefluidmesh[27,28].Furthermore,the indices eand k in Eq.(10b) iterate over the total numberof el-ements andthe local node number in the element, respectively. In thismethod, the test functions

ψ

_if are the nodal shape func-tionsofthe fluidnodes(•,Fig.2).Theyhaveanon-zerovalue in the shaded region inFig. 3b. Similarly, thetest functions

ψ

s

i are thenodalshape functionsofthesolid nodes(,Fig. 2)andhave

Fig. 2. Schematic representation of the solid ( ψs ), and fluid ( _ψf , _ξf ) linear test functions in a one-dimensional case. The shaded element represents the interface region between the fluid and solid domains. The blue line on top represents the variation of αacross the one-dimensional computational domain. ψs and _ψf are used in the transport equations in the solid and fluid domains, respectively. ξf is used to enforce the continuity constraint in the fluid domain. Note that _ξf is discontinuous at the solid end of the interface ( i + 3 ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

anon-zero value in the shaded region in Fig.3c.Unlike

ψ

_if and

ψ

s

i whicharedefinedontheglobalnodebasis,thefunctions

ξ

f i in

Eq.(10b)aredefinedonthelocalnodebasis. Itisessentially

ψ

_if, withouttherequirementthatitsvalueiszeroatthefluid-solid in-terface.Asexplainedbefore,thefunction

ξ

f_{isdiscontinuousatthe} solidendofthedomain(Fig.2,nodei+3).

Eqs. (8) and (7) can be written in discrete acceleration form as,

(

M+

θ

t

(

A− D

))

cn+1

− c_t n =−

(

A− D

)

cn+f+fs, (11a)

¯CTun+1=0. (11b)

In Eq.(11), M is the mass matrix, A is the advection matrix, andDisthediffusionmatrix.Additionally,¯CTisthediscrete diver-genceoperator(superscript Tdenotesthetranspose).The variable

tdenotesthetimestepthatcanvaryinthecourseofthe simula-tiontosatisfyacertainCFLvalue,and

θ

isatimedifferencing pa-rametertochoosebetweendifferenttemporalintegrationschemes (

θ

=0foranexplicitscheme,

θ

=0.5forthesemi-implicit

Crank-Nicolsonscheme,and

θ

=1fora fullyimplicitdiscretisation). Fi-nally,fs_{istheimmersedbodyforcing termandfrepresentsother} forcingterms(e.g.gravityforce).Morespecifically,thediscrete op-eratorsandthe forcing terms inthe discrete equationsare given by, Mi j=  

(

ψ

f i +

ψ

s i

)

φ

jdV, (12a) Ai j=−  

∇ψ

f i · u n

_φ

jdV, (12b) Di j=  

∇ψ

f i · ¯¯

κ

·

∇φj

dV, (12c) ¯CTi j=  

ξ

f i

∇φj

dV, (12d) fi=  

ψ

f i

φ

jFjndV, (12e) fs i =  

ψ

s i

φ

j cs_j(n+1) − cn j

t dV. (12f)

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 5

Fig. 3. Schematic representation of the solid (a) and fluid (b, c, d) domains. Figure (a) shows the solid grid ( s ). Figures (b, c, d) show the fluid grid ( _f ), with circular markers ( ◦) representing fluid nodes ( αn+1 = 0 ) and square markers ( ) representing solid nodes ( _αn+1 > 0 ), where strong velocity boundary conditions are imposed. Shaded

elements represent: (b) the domain in which the fluid operators act, (c) the domain in which the solid operator acts, and (d) the interface region between the solid and fluid domains.

Fig. 3b and c show the decomposed domains in which these operators act. The operator M acts on the whole domain, whilst A,D,CT,andfactsintheshadedfluid domaininFig.3b.The im-mersed body forcing term, fs_{,actsin theshaded solid domainin}

Fig.3c.

3.2. ApplicationtotheNavier–Stokesequations

The discrete Navier–Stokesmomentum equationsareobtained from Eq. (11a) by treating velocity as a transport variable and addingCpn+θ _{totheright-hand-side,as}

(

M+

θ

t

(

A− D

)

)

u∗

− u_t n =−

(

A− D

)

un_+Cpn+θ_+f+fs_{, (13)} where, Ci j=−  

ψ

f i

∇ξ

f jdV, C T i j=  

ξ

f i

∇ψ

f jdV. (14)

Boldletters denote that thematrices (M,A,D,f,and fs) are now inblockmatrix formandu∗ isthepredictedvelocity.The matrix –C is the discrete pressure gradient operator. Note that we used

ξ

f _{aspressurebasisfunctions.Thisisbecause,justlikecontinuity,} pressureisalsodefinedonlyin



f_.

WhensolvingtheNavier–Stokesequations,itiscommonto de-couplethesolutionprocedures forvelocity andpressure. Chorin’s projection method [29] isused here forfractional time stepping. TheprocedurestartsbysolvingthepressurePoisson’sequation to computean initialpressurefield. ThePoisson’sequation for

pres-sureisobtainedbytakingthedivergenceofEq.(13)andapplying theincompressibility condition(Eq.(11b)) onthepredicted veloc-ity,i.e. ¯Lpθ=− ¯CT



M−1RHS+

u0 t



+ ¯C T u∗

t . (15)

Here, RHS contains all the terms on the right hand side of

Eq.(13)exceptthepressuregradientterm(Cpn+θ_),and¯Lisa dis-creteLaplacianexpressed as ¯L= ¯CTM−1C.However, ¯L isnot sym-metric.AsymmetricPoisson’ssystemofequationscanbeobtained byreformulatingthecontinuityequationforvelocityusingCT,as

CTu∗=−CT(s) _u∗_{≡ −C}T(s) _us(n+1) _{, (16a)} Lpθ=−CT



M−1RHS+

u0 t



−C T(s) _us(1)

t . (16b) Here, CT_{i j}(s) =  

ξ

f i

∇ψ

s jdV ≡  

ξ

f i

∇

(

φ

j−

ψ

f j

)

dV ≡ ¯CTi j− C T i j, (17a) L=CTM−1C. (17b)

An intermediatevelocity field (u∗) isthen obtainedbysolving themomentum equation (Eq. (13)). Ingeneral, the predicted ve-locity field is not divergence-free and will not satisfy Eq. (16a). Theprojectionmethodtacklesthisissuebydefininganirrotational pressure correction, determined by solving the following Poisson equation, L

pn+θ₌−

(

CT+CT(s)

)(

θ

u∗+

(

1−

θ

)

un

)

t ≡− ¯C T

(

θ

u∗+

(

1−

θ

)

un

₎

t , (18)

andcorrectingthevelocityfieldas,

un+1_=u∗₊

_tM−1_C

_pn+θ_{. (19)}

It can be seen from Eqs. (16b) and (18) that, at the interface, a smearedzerogradient (∂_∂np=0) boundaryconditionis appliedon thepressurefield.Theboundaryconditionissmearedinthesense thatthegradientgradually decreasesacross theinterfaceand be-comeszero attheend ofthe fluiddomain. Finally,the Laplacian operatorisdefinedonlyin



f_{.Therowsandcolumnsinthe} dis-creteLaplacian,corresponding tothenodes in



s_ࢨ

_∂

_

s_,wouldbe filled withzeros. Hence, while solving Eqs. (16b) and (18)those rowsandcolumnshavetobesetinactive.Thiscanbedoneintwo ways.Theefficient methodisto discardtheserows andcolumns fromthesolutionprocedureandsolveasubsysteminstead.Inthat case,theconditionnumberofthematrixremainsunchangedand thesizeofthe linearsystemtobe solvedis reduced.Thesecond optionistoadd unityto thediagonalentries ofthosenodesand solvethefullsystem.However,dependingonthenumberofnodes insidethe solid domain theconditionnumber ofthe matrixwill beworsened.

Inthispaper,Eqs.(11),(13),(16),and(18)aresolvedusingthe fluiddynamics modelFluidity, an open-source finiteelement nu-merical tool that solves the Navier–Stokes equations on unstruc-turedgrids[30–32].

3.3.Treatmentofmovinginterfaces

In a deforming/moving body problem, elements and nodes moveacrossthefluid-structure interface.It isknownthat,for di-rectforcing immersedboundarymethods,thiscancause numeri-cal oscillations [33]. Onesource of theseoscillations comes from

Fig. 4. Decomposed solid and fluid domains of a moving body on a finite ele- ment grid at time level n + 1 . Circular markers ( ◦) represent fluid nodes ( αn+1 = 0 ).

Square markers ( ) represent solid nodes ( αn+1 = 1 ), where strong velocity bound-

ary conditions are imposed. Shaded elements represent the domain in which regu- lar fluid operators act.

thedifferenceintheorderoftheerrorbetweenthefluidandsolid domains.Inthepresentmethod,boththesolidandfluiddomains arespatiallydiscretisedusingthesameshapefunctions.Hence,the numericalerrorsassociatedwithbothdomainsareofthesame or-derofmagnitude.However,thetemporaldiscontinuityinthe pro-jectionstep (Eq. (18)) is amplifiedby a factorof1/

t.This is il-lustratedbyFig.4.Thecontinuouscirclehighlightsthepositionof the body attime leveln+1,while the dashed circleoutlines its positionattimeleveln.Atthistimelevel,thenodethatismarked withadarksquare()isinternaltothesoliddomain.Pressurehas nodefinedvalueinsidethesolidandissettozero.However,when thebodymoves,thisnodebecomespartofthefluid-solidinterface andstartsinteractingwiththefluidequationsviathepressure gra-dientterm.Sincevelocityhasdefinedvaluesinsidethesolid, oper-atorsthatact onthevelocityfield areunaffected.However,errors areintroducedintheevaluationofthepressuregradientandwill propagateontothevelocityfield.Toremedythis,thepressure gra-dient termhas to be moved entirely to the projection step. The modifiediterationschemeisthen

(

M+

θ

t

(

A− D

)

)

u∗

− u_t n =−

(

A− D

)

un_+f+fs_{, (20)} Lpn+θ₌− ¯C T

(

θ

u∗+

(

1−

θ

)

un

₎

t , (21) un+1_=u∗₊

_tM−1_Cpn+θ_{. (22)}

Inthat case,thePoisson’ssystemsolves thefullpressure field (p)andnotacorrectionforpressure(

p).

4. Results

The method described in the previous section is applied to flowspastfixedandmovingstructures.Allthecasesconsideredin thissectionaresolvedontwo-dimensionalunstructuredgridswith triangular elementsanda CFL number equal to0.95. The Crank-Nicolsonmethodisusedfortemporalintegration.

4.1. Taylor-Couetteflow

Taylor-Couetteflowistheflowbetweentwoconcentrichollow cylinders[34].TheReynoldsnumberforthisflowisdefinedbased

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 7

Fig. 5. Taylor-Couette flow results. (a) Projected solid concentration field on the fluid mesh, the domain shaded in red being the inner solid cylinder whilst the region of zero solid concentration (shaded in blue) being filled with fluid. (b) Computed velocity on the fluid mesh, where inside the solid domain, the velocity is equal to the projected solid velocity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

on the relative angular velocity of the cylinders and the diame-ter of one of the cylinders. At low Reynolds numbers, the flow issteadyandpurely azimuthal.The criticalvalue oftheReynolds numbercorresponds tothe pointatwhichthe flowbecomes un-stableandthree-dimensional.Belowthecriticalpoint,ananalytical solutionexistsandisgivenby,incylindricalcoordinates,

u_θ=Ar+B r, ur=uz=0, A=



in

μ

−

η

2 1−

η

2, B=



inR 2 in 1−

μ

1−

η

2

μ

=



out



in,

η

= Rout Rin, (23) whereRistheradius,



istheangularvelocity, andsubscriptsin andoutdenotethecylindertowhich theparameterbelongs.The existence of an analytical solution makes Taylor-Couette flow an ideal case for grid convergence study. For the current test case, flow parameters are chosen as follows:



in=4,



out =0,Rin= 0.5,Rout=1,andthekinematicviscosity

ν

=1.TheReynolds num-berbasedonthediameteroftheinnercylinderisRe_in=2_,which is wellbelow theknowncritical Reynoldsnumberforthissetup, i.e.Recritical

in =68.4.

To modelthis problemusing the currentimmersed boundary method,thehollowcylindersarereplacedbysolid discs.The mo-tion ofthe inner cylinderis prescribed usingan immersedsolid, whichrotatesatanangularvelocityof



in=4.Thedomainshaded in red inFig. 5a representsthe solid inner discthat isprojected onto the fluid mesh. Inside this solid domain, velocity is pre-scribedandvarieslinearlyalongtheradialdirection,asshownby

Fig.5b.Thefixedoutercylinderismodelledusingaregular body-conformingwall.Thefluidfieldisinitialisedasatrestandthe im-mersedbodyisimpulsivelysetintomotion.

Resultsareplottedatsteady-state.Fig.6aandbshowthe varia-tionintangentialvelocityintheradialdirection.InFig.6a,thered lineshowsthe exactsolution andtheblack linesshow the com-putedprofilesatvaryingspatialresolutions.Duringthegrid refine-mentstudy,theelementedgelengthisvariedonlyattheinterface located at r=0.5, that is atthe boundary of the inner cylinder.

Bycontrast,neartheouter-wall,theedgelengthiskeptata con-stantvalue of

x=2× 10−2_{.The profile shows good agreement} withtheanalyticalsolutionandconvergenceisobservedwhenthe meshdensityisincreased.Byfocusingtheincrementinspatial res-olutionattheinterface,itcanbeinferredfromFig.6athatthe er-rorisamaximumattheinterface(r=0.5).Thisisaconsequence oftheconservativeprojectionontoanon-matchinggrid.Inthe in-setofFig.6a, itcan be seenthat thepoint fromwhichthe solid concentrationfield (bluelines)startstodropvariesdependingon themeshresolution.Becauseoftheconservativeprojection,the ra-dius of the projected cylinder is slightly smaller than the actual solidbody.Since,inthisparticularcase,thevelocityboundary con-ditionthat hastobe imposedatthe fluid-structure interfaceis a function ofthe radius of the cylinder(u_θ(r)), the imposed angu-larvelocityattheinterfaceisreduced.Asshowninthefigure,this effectdisappearsforsufficientlyfinegridsatthefluid-structure in-terface.

Fig. 6b shows the L_∞ (◦) and L2 () spatial velocity error (



ue



) norms plotted against the element edge length at the fluid-structure interface.The spatialerrornorms are evaluated at steady-statebycomparingtheresultswiththeanalyticalsolution. Convergencestudyshowsan orderofabout1.31and2.06forthe L_∞ andL2errornorms,respectively.

4.2.Flowpastastationarycylinder

Inthissection, flow past a stationarycircularcylinderis con-sideredasa test problem. Asshown inFig. 7,the computational domain hasthe dimensions25d × 20d(dbeing the cylinder di-ameter)andthecylinderis centredatadistance5dfromthe in-let and10dfromthe sides. Onthe left boundary,that is marked with arrows pointing inwards, a uniform velocity profile is pre-scribedasinletboundarycondition.Openboundaryconditionsare applied on all the other domain boundaries. Fig. 7 also shows thecharacteristic dimensionsofthewake structure asdefinedin CoutanceauandBouard[36].Dimensionsl,a,andbrepresentthe length ofthe recirculationzone, distance to thevortex core, and gap betweenthe vortex cores, respectively. The separation angle

Fig. 6. (a) Variation of tangential velocity in the radial direction in a Taylor-Couette flow. Red line shows the analytical solution [34] . Black lines show results using the current method at different grid resolutions close to the fluid-structure interface. Blue lines in the inset show the variation of the solid concentration field ( α). (b) The L ∞ :

 and L 2 :  spatial velocity error (  u e ) norms plotted against the element edge length at the interface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Dimensions of the computational domain in terms of the cylinder diameter d and characteristic dimensions of the wake structure, as defined in [24] .

isdenoted by

θ

. TwoReynolds numbers are considered:Re=20 andRe₌40. The steady-stateresults fromthe presentimmersed boundarymethodarecomparedwithseveralresultsfromthe

lit-Fig. 9. Non-dimensionalised inline force acting on the cylinder in one cycle at Re = 100 and KC = 5. Black lines are results obtained using the immersed boundary method presented in this paper at varying grid resolutions near the interface. Black cross markers are from experimental data [40] .

Fig. 8. Time averaged coefficients of pressure and friction from flow past a cylinder at a Reynolds number of 200 with varying edge lengths near the interface and compared to DNS data from existing literature [39] . θis measured clockwise from the stagnation point. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 9

Fig. 10. Comparison of non-dimensionalised velocity components at four cross-sections between present computation (lines) and experimental data [40] (markers) at three different phases within one cycle of oscillation: (a) 180 ◦_{, (b) 210} ◦, and (c) 330 ◦_.

erature[24,35–38].In particular, the characteristiclengths ofthe flow (non-dimensionalised by the cylinder diameter d) and the dragcoefficientarepresentedinTable1.Itisshownthattheflow profiles andthe force coefficient fromthe present methodagree wellwiththosereportedintheliterature.

Furthermore,anunsteadytestcaseisconsideredataReynolds number of 200. Fig. 8 shows time averaged pressure and fric-tion coefficients along the periphery of the cylinder. Coefficients are averaged using a simple mean and

θ

is measured from the stagnationpoint. Results are plotted at three different grid edge lengths near the boundary of the cylinder: le=0.01d (dash-dot line), le=0.0075d (dash-dash line), and le=0.005d (solid black line).Resultsshowgridcovergenceandagreesverywellwithdata availablefrom abody-conforming DNSsimulation(redline) [39]. Aswiththesteadycases,itcanbenoticedfromthecoefficientof frictioncurvesthattheaveragepointofseparationispredicted ac-curately.

Fig. 11. Velocity field and streamlines at three different phases within one cycle of oscillation: (a) 180 ◦_{, (b) 210} ◦, and (c) 330 ◦_.

Table 1

Steady state non-dimensionalised wake dimensions and coefficient of drag from flow past a cylinder at Reynolds numbers of 20 and 40.

l / d a / d b / d θ CD Re = 20 Dennis and Chang [35] 0.94 – – 43.7 ◦ _2.05

Coutanceau and Bouard [36] 0.93 0.33 0.46 45.0 ◦ _–

Tritton [37] – – – – 2.09

Linnick and Fasel [38] 0.93 0.36 0.43 43.5 ◦ _2.06

Taira and Colonius [24] 0.94 0.37 0.43 43.3 ◦ _2.06

Present study 0.94 0.34 0.43 43.8 ◦ _2.05

Re = 40 Dennis and Chang [35] 2.35 – – 53.8 ◦ _1.52

Coutanceau and Bouard [36] 2.13 0.76 0.59 53.8 ◦ _–

Tritton [37] – – – – 1.59

Linnick and Fasel [38] 2.28 0.72 0.60 53.6 ◦ _1.54

Taira and Colonius [24] 2.30 0.75 0.60 53.7 ◦ _1.54

Present study 2.30 0.71 0.61 53.3 ◦ _1.55

4.3. Flowpastanoscillatingcylinder

In this section, the vortices generated by an oscillating cylin-derplacedinstillfluidaresimulatedandcomparedagainst exper-imental data[40].The cylinderisrigidandirrotational.The hori-zontaloscillationofthecylinderisdefinedby,

x

(

t

)

=−Asin

(

2

π

ft

)

, u

(

t

)

=−2

π

f Acos

(

2

π

ft

)

, (24) where,x(t)isthepositionofthecentreofthecylinderandu(t)is thevelocityofthecylinder.TheparametersAandfarethe ampli-tudeandfrequencyoftheprescribedoscillation.Reynoldsnumber andKeulegan-Carpenternumber(KC)arefurtherdefinedas,

Re=umaxd

ν

, KC=

umax

f d . (25)

For the simulation, the cylinder withdiameter d is placed at thecentreofacomputationaldomainofsize55d_{× 35}d.Reynolds numberissetat100 andKeulegan-Carpenternumberissetat5. Time(t)andhorizontalforce(F)arenon-dimensionalisedast˜=t f andF˜₌F_/

ρ

dumax,respectively.Fig.9showsthevariationof non-dimensionalisedinlineforce(F˜)against non-dimensionalisedtime (t˜), for one period of oscillation, and at varying grid resolutions near the interface: le=0.05d (dash-dot line), le=0.025d (dash-dash line), and le=0.0125d (solid black line). The results have goodagreementwiththereferencedata[40]andshowsgrid con-vergence.

Fig.10compares velocitycomponentsoftheflow fieldatfour cross-sectionsalongthelineofoscillationwithexperimentaldata

[40].Forthemostpart,theprofiles agreeverywell. Whereit

dif-Fig. 12. Coefficient of pressure on the surface of the cylinder at a phase angle of 270 ◦_{with varying edge lengths near the interface and compared to DNS data from}

existing literature [42] . θis measured clockwise from the stagnation point. (For in- terpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 11

Fig. 13. Time evolution of the coefficients of lift (top) and drag (bottom) on a NACA0012 aerofoil at Re = 10 0 0 and an angle of attack 10 ◦_{. Black lines are results obtained}

using the immersed boundary method presented in this paper at varying grid resolutions near the interface. Red line is a body-conforming simulation result using the solver Fluidity, and black cross markers are DNS results from existing literature [43] . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

fers, the mismatch is only on one side of the symmetric plane (y=0) andis dueto theasymmetry intheexperimentaldata.In

Fig.11,streamlinesinthevicinityofthecylinderandthevelocity fieldareplottedatdifferentphasesinthecycle.Streamlinecurves areagoodmatch,andcanbecomparedagainstplotsfromthe ref-erence[40]atthesamephaseinthecycle.

As an extension, we now add a highenough background ve-locitytothe previouslyquiescentfluid suchthat thecylindercan naturallyshed vortices even whenit is stationary.From hereon, we switchback totheflowvelocity baseddefinitionforReynolds number(insteadofumaxoftheoscillation),anduseStrouhal num-ber(insteadofKeulegan-Carpenternumber)todefineourtestcase. For our test case, we set the Reynolds number at 185, and this gives us a Strouhal number of 0.193 [41]. Additionally, the am-plitude ofoscillation (A) isset as0.2d, and the ratio of the fre-quencies of excitationand naturalvortexshedding(f/fn) is setto 0.9. Fig.12 plots thecoefficient ofpressureon thesurfaceofthe cylinderataphaseangleof270◦,(centreofthecylinder,x

(

t

)

=A) for three different grid resolutions near the immersed interface: le=0.01d (dash-dot line), le=0.0075d (dash-dashline), andle= 0.005d (solidblackline).Theresultsmatchwiththereferencedata (redline)[42]andshowsgridconvergence.

4.4. FlowpastaNACA0012aerofoil

Flow pasttheNACA0012aerofoilis consideredtodemonstrate the capabilities of the method with moderately turbulent flows. The chord-based Reynolds numberis 1000.The size ofthe com-putationaldomainisthesameasinFig.7,exceptthatdimensions are here expressed interms of the chord length (c). The bound-ary conditions are also the same. Fig. 13 shows the time evolu-tion oftheliftanddragcoefficientsatanangleof attack10◦ us-ing the present immersed boundary method at three grid edge lengths near the aerofoil boundary: le=0.0025c (dash-dot line), le=0.00175c (dash-dash line), and le=0.001c (solid black line). The results are further compared with a body-conforming DNS

available inthe literature [43] (black crossmarkers) anda body-conformingsimulationusingFluidityandle=0.0025c (redline).It canbeseenfromtheplotsthat thereisadifferenceinamplitude oftheoscillations betweenthe bodyconforming simulationfrom Fluidity(redline)andtheDNSdataobtainedfromliterature(black crossmarkers). The relativedifference inmagnitude betweenthe body-conforming results is maximum4% for the dragcoefficient and1% fortheliftcoefficient.Thiscould beexplainedby theuse oflinear shape functionsfor both velocity andpressure fields in thepresentsimulationsandtheuseoffullyunstructuredmeshes. Theagreement betweenthe presentimmersedboundary method andthebody-conformingsimulationusingFluidityishoweververy good.Inparticular,theliftcoefficientsfromalltheFluidity simula-tionsalmostmatchperfectlywithoneanother.Thedragcoefficient ismoresensitivetothemeshresolutionusedneartheaerofoil.The dragcoefficientgivenbytheimmersedboundarymethoddeviates by ~ 5% withthe coarsermesh (dash-dot line)compared to the body-conformingcase(redline),andthispercentagedecreasesas themeshisrefined,to ~ 3% for the intermediatemeshand ~ 1% forthefinestmesh.

5. Conclusions

Thispaperpresentsanewformulationoftheimmersed bound-arymethod that isspecific to finite-element discretisation meth-ods. The proposed methodmodifies the test functionsin such a waythatthe spacesoccupiedby thefluid andsolid, respectively, are decomposedinto two sub-domains.Unlikemanyof the forc-ing immersed boundary methods available in the literature, the presentmethodsatisfiesboththedivergence-freeandno-slip con-ditionswithoutintroducinganadditionalforcing terminthe mo-mentumequationsofthefluid.Theabsenceofaforcing term im-plies that the overall CFL limit is not affected by the immersed boundary implementation. Furthermore, the modified test func-tionspreservesymmetryandpositivedefinitenessoftheLaplacian

operator,sothatthePoisson’ssystemcanbeefficientlysolved us-ingexistingnumericalmethods.

It is shown that the method provides accurate results on a series of test cases with both fixed and moving bodies, as well as laminar and moderately turbulent flows. In particular, excel-lentagreement was found between the presentmethod and ex-perimental data for low-Reynolds numbers flows. For flow past an aerofoil at a higher Reynolds number, an accurate prediction oftheliftcoefficient wasobtainedatthe samegridsize require-mentofabody-conformingsimulation.Thepresentmethod how-everrequiredaslightlyfinermeshresolutioncomparedtoa body-conformingmethod,inordertoprovideadragcoefficientofequal accuracy.Futureworkwillfocusonextendingthepresentmethod tohighlyturbulentflowsofinterestforengineeringapplications.

DeclarationofCompetingInterest

Theauthorsdeclarethattheyhavenoknowncompeting finan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto influencetheworkreportedinthispaper.

CRediTauthorshipcontributionstatement

Navaneetha Krishnan: Conceptualization, Methodology, Writ-ing-originaldraft.AxelleViré:Conceptualization, Writing- orig-inaldraft, Supervision,Fundingacquisition. RolandSchmehl: Su-pervision,Writing-review &editing,Fundingacquisition.Gerard vanBussel:Resources,Supervision.

Acknowledgements

A. Viré and N. Krishnan are financially supported by the

European Union SeventhFramework Programme (FP7/2007-2013) under a Marie Curie Career Integration Grant (Grant Agree-mentPCIG13-GA-2013-618159).A.Viré alsoacknowledgessupport fromtheRijksdienstvoorOndernemendNederland(RVO)through the TSE Hernieuwbare Energie funding scheme (ABIBA project - TEHE116332). R. Schmehl is financially supported by the Eu-ropean Commission Horizon 2020 framework programme under AWESCO(H2020-ITN-642682) andREACH(H2020-FTIPilot-691173) projects.ThisworkwaspartlycarriedoutontheDutchnational e-infrastructurewiththesupportofSURFCooperative(project num-ber17215).

References

[1] Iaccarino G, Verzicco R. Immersed boundary technique for turbulent flow sim- ulations. Appl Mech Rev 2003;56(3):331–47. doi: 10.1115/1.1563627 .

[2] Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech 2005;37(1):239–61. doi: 10.1146/annurev.fluid.37.061903.175743 .

[3] Sotiropoulos F, Yang X. Immersed boundary methods for simulating fluid- structure interaction. Prog Aerosp Sci 2014;65:1–21. doi: 10.1016/j.paerosci. 2013.09.003 .

[4] Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972;10(2):252–71. doi: 10.1016/0 021-9991(72)90 065-4 .

[5] Arquis E , Caltagirone JP . Sur les conditions hydrodynamiques au voisinage d’une interface milieu fluide-milieux poreux: application à la convection na- turelle. In: Comptes rendus de l’academie des sciences, 299; 1984. p. 1–4 .

[6] Kolomenskiy D, Schneider K. A fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles. J Comput Phys 2009;228(16):5687–709. doi: 10.1016/j.jcp.2009.04.026 .

[7] Gazzola M, Chatelain P, van Rees WM, Koumoutsakos P. Simulations of sin- gle and multiple swimmers with non-divergence free deforming geometries. J Comput Phys 2011;230(19):7093–114. doi: 10.1016/j.jcp.2011.04.025 .

[8] Boffi D, Gastaldi L. A finite element approach for the immersed bound- ary method. Comput Struct 2003;81(8):491–501. doi: 10.1016/S0045-7949(02) 00404-2 .

[9] Zhang L, Gerstenberger A, Wang X, Liu WK. Immersed finite element method. Comput Methods Appl Mech Eng 2004;193(21):2051–67. doi: 10.1016/j.cma. 2003.12.044 .

[10] Liu WK, Liu Y, Farrell D, Zhang L, Wang XS, Fukui Y, et al. Immersed finite element method and its applications to biological systems. Comput Methods Appl Mech Eng 2006;195(13):1722–49. doi: 10.1016/j.cma.2005.05.049 .

[11] Viré A, Xiang J, Pain CC. An immersed-shell method for modelling fluid- structure interactions. Philos Trans R Soc Lond A 2015;373(2035):1–14. doi: 10. 1098/rsta.2014.0085 .

[12] Goldstein D, Handler R, Sirovich L. Modeling a no-slip flow boundary with an external force field. J Comput Phys 1993;105(2):354–66. doi: 10.1006/jcph.1993. 1081 .

[13] Mohd-Yusof J . Combined immersed-boundary/b-spline methods for simula- tions of flow in complex geometries. In: Center for turbulence research annual research briefs. Stanford University; 1997. p. 317–27 .

[14] Fadlun E, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 20 0 0;161(1):35–60. doi: 10.10 06/jcph.20 0 0.6484 .

[15] Uhlmann M. An immersed boundary method with direct forcing for the sim- ulation of particulate flows. J Comput Phys 2005;209(2):448–76. doi: 10.1016/j. jcp.2005.03.017 .

[16] Su S-W, Lai M-C, Lin C-A. An immersed boundary technique for simulating complex flows with rigid boundary. Comput Fluids 2007;36(2):313–24. doi: 10. 1016/j.compfluid.20 05.09.0 04 .

[17] Ye T, Mittal R, Udaykumar H, Shyy W. An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries. J Comput Phys 1999;156(2):209–40. doi: 10.1006/jcph.1999.6356 .

[18] Kirkpatrick M, Armfield S, Kent J. A representation of curved boundaries for the solution of the Navier-Stokes equations on a staggered three-dimensional Cartesian grid. J Comput Phys 2003;184(1):1–36. doi: 10.1016/s0021-9991(02) 0 0 013-x .

[19] Cheny Y, Botella O. The ls-stag method: a new immersed boundary/level- set method for the computation of incompressible viscous flows in com- plex moving geometries with good conservation properties. J Comput Phys 2010;229(4):1043–76. doi: 10.1016/j.jcp.2009.10.007 .

[20] Meyer M, Devesa A, Hickel S, Hu X, Adams N. A conservative immersed inter- face method for large-eddy simulation of incompressible flows. J Comput Phys 2010;229(18):6300–17. doi: 10.1016/j.jcp.2010.04.040 .

[21] Kim J, Kim D, Choi H. An immersed-boundary finite-volume method for sim- ulations of flow in complex geometries. J Comput Phys 2001;171(1):132–50. doi: 10.10 06/jcph.20 01.6778 .

[22] Ge L, Sotiropoulos F. A numerical method for solving the 3D unsteady in- compressible Navier–Stokes equations in curvilinear domains with complex immersed boundaries. J Comput Phys 2007;225(2):1782–809. doi: 10.1016/j.jcp. 2007.02.017 .

[23] Glowinski R, Pan T, Priaux J. Distributed Lagrange multiplier methods for in- compressible viscous flow around moving rigid bodies. Comput Methods Appl Mech Eng 1998;151(1):181–94. doi: 10.1016/s0 045-7825(97)0 0116-3 .

[24] Taira K, Colonius T. The immersed boundary method: a projection approach. J Comput Phys 2007;225:2118–37. doi: 10.1016/j.jcp.2007.03.005 .

[25] Boffi D, Cavallini N, Gastaldi L. The finite element immersed boundary method with distributed lagrange multiplier. SIAM J Numer Anal 2015;53(6):2584–604. doi: 10.1137/140978399 .

[26] Elman HC, Silvester DJ, Wathen AJ. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Numerical Mathematics and Scientific Computation. Oxford University Press; 2005. doi: 10.1093/acprof: oso/9780199678792.001.0001 .

[27] Farrell PE, Piggott MD, Pain CC, Gorman G, Wilson C. Conservative interpola- tion between unstructured meshes via supermesh construction. Comput Meth- ods Appl Mech Eng 2009;198(33–36):2632–42. doi: 10.1016/j.cma.2009.03.004 .

[28] Viré A, Xiang J, Milthaler F, Farrell PE, Piggott MD, Latham J-P, et al. Modelling of fluid-solid interactions using an adaptive mesh fluid model coupled with a combined finite discrete element model. Ocean Dyn 2012;62(10–12):1487–501. doi: 10.1007/s10236- 012- 0575- z .

[29] Chorin AJ. Numerical solution of the Navier-Stokes equations. Math Comput 1968;22:745–62. doi: 10.1090/S0025- 5718- 1968- 0242392- 2 .

[30] Pain CC, Umpleby AP, de Oliveira C, Goddard AJH. Tetrahedral mesh op- timisation and adaptivity for steady-state and transient finite element cal- culations. Comput Methods Appl Mech Eng 2001;190:3771–96. doi: 10.1016/ S0 045-7825(0 0)0 0294-2 .

[31] Pain C, Piggott M, Goddard A, Fang F, Gorman G, Marshall D, et al. Three- dimensional unstructured mesh ocean modelling. Ocean Modell 2005;10(1– 2):5–33. doi: 10.1016/j.ocemod.20 04.07.0 05 .

[32] Piggott MD, Gorman GJ, Pain CC, Allison PA, Candy AS, Martin BT, et al. A new computational framework for multi-scale ocean modelling based on adapting unstructured meshes. Int J Numer Methods Fluids 20 08;56:10 03–15. doi: 10. 1002/fld.1663 .

[33] Luo H, Dai H, de Sousa PJF, Yin B. On the numerical oscillation of the direct- forcing immersed-boundary method for moving boundaries. Comput Fluids 2012;56:61–76. doi: 10.1016/j.compfluid.2011.11.015 .

[34] Taylor GI. Stability of a viscous liquid contained between two rotating cylin- ders. Philos Trans R Soc A 1923;223(605–615):289–343. doi: 10.1098/rsta.1923. 0 0 08 .

[35] Dennis SCR, Chang G-Z. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J Fluid Mech 1970;42(3):471–89. doi: 10.1017/S0 022112070 0 01428 .

[36] Coutanceau M, Bouard R. Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform trans- lation. Part 1. Steady flow. J Fluid Mech 1977;79(2):231–56. doi: 10.1017/ S0 0221120770 0 0135 .

[37] Tritton DJ. Experiments on the flow past a circular cylinder at low Reynolds numbers. J Fluid Mech 1959;6(4):547–67. doi: 10.1017/S0 0221120590 0 0829 .

N. Krishnan, A. Viré and R. Schmehl et al. / Computers and Fluids 202 (2020) 104500 13 [38] Linnick M, Fasel H. A high-order immersed interface method for simulat-

ing unsteady incompressible flows on irregular domains. J Comput Phys 2005;204(1):157–92. doi: 10.1016/j.jcp.2004.09.017 .

[39] Rajani B, Kandasamy A, Majumdar S. Numerical simulation of laminar flow past a circular cylinder. Appl Math Model 2009;33(3):1228–47. doi: 10.1016/j. apm.2008.01.017 .

[40] Dütsch H, Durst F, Becker S, Lienhart H. Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers. J Fluid Mech 1998;360:249–71. doi: 10.1017/s0 0221120980 0860x .

[41] Williamson CHK. Defining a universal and continuous Strouhal-Reynolds num- ber relationship for the laminar vortex shedding of a circular cylinder. Phys Fluids 1988;31(10):2742–4. doi: 10.1063/1.866978 .

[42] Guilmineau E, Queutey P. A numerical simulation of vortex shedding from an oscillating circular cylinder. J Fluids Struct 2002;16(6):773–94. doi: 10.1006/jfls. 2002.0449 .

[43] Mittal S, Tezduyar TE. Massively parallel finite element computation of incom- pressible flows involving fluid-body interactions. Comput Methods Appl Mech Eng 1994;112(1):253–82. doi: 10.1016/0 045-7825(94)90 029-9 .