Animmersedinterfacemethodforsimulatingtheinteractionofa - PDF document

Animmersedinterfacemethodforsimulatingtheinteractionofa”uidwithmovingboundariesShengXu,Z.JaneWangDepartmentofTheoreticalandAppliedMechanics,CornellUniversity,KL214KimballHall,Ithaca,NY14853,USAReceived21March2005;receivedinrevisedform1December2005;accepted16December2005Intheimmersedinterfacemethod,boundariesarerepresentedassingularforceintheNavier…Stokesequations,whichentersanumericalschemeasjumpconditions.Recently,wesystematicallyderivedallthenecessaryspatialandtemporaljumpconditionsforsimulatingincompressibleviscous”owssubjecttomovingboundariesin3Dwithsecond-orderspatialandtemporalaccuracyneartheboundaries[ShengXu,Z.JaneWang,Systematicderivationofjumpconditionsfortheimmersedinterfacemethodinthree-dimensional”owsimulation,SIAMJ.Sci.Comput.,2006,inpress].Inthispaperweimplementtheimmersedinterfacemethodtoincorporatethesejumpconditionsina2Dnumericalscheme.Westudytheaccuracy,e
ciencyandrobustnessofourmethodbysimulatingTaylor…Couette”ow,”owinducedbyarelaxingballoon,”owpastsingleandmultiplecylinders,and”owarounda”appingwing.Ourresultsshowthat:(1)ourcodehassecond-orderaccuracyinthein“nitynormforboththevelocityandthepressure;(2)theadditionofanobjectintroducesrelativelyinsigni“cantcomputationalcost;(3)themethodisequallyeectiveincomputing”owsubjecttoboundarieswithpre-scribedforceorboundarieswithprescribedmotion.2006ElsevierInc.Allrightsreserved.Immersedinterfacemethod;Immersedboundarymethod;Cartesiangridmethod;Movingdeformableboundaries;Complexgeometries;Flowaroundmultipleobjects;Singularforce1.IntroductionItisimportanttoaccuratelyande
cientlyresolvemovingboundariesandtheireectson”uidsinnumer-icalsimulationsof”uid…structureinteractions.Body-“ttedgridmethods,whichemploystructuredorunstruc-turedbody-“ttedgrids,aredesignedtoresolveboundariesandtheireectswithhigh-orderaccuracy,butitiscomputationallycostlytoupdategridsinmovingboundaryproblems.VariousCartesiangridmethodshavebeendevelopedtoavoidthegridregeneration.Theyallowforfast”owsolversandhavetheadvantagesofsimplicityande
ciency.OneclassofCartesiangridmethodsaresuitableforsolving”owswithmovingboundariesundergoingprescribedmotion.Examplesincludethevirtual/immersedboundarymethod0021-9991/$-seefrontmatter2006ElsevierInc.Allrightsreserved.doi:10.1016/j.jcp.2005.12.016 Correspondingauthor.Tel.:+16072548593.E-mailaddresses:sx12@cornell.edu(S.Xu),jane.wang@cornell.edu(Z.J.Wang). JournalofComputationalPhysicsxxx(2006)xxx…xxx www.elsevier.com/locate/jcp ARTICLEINPRESS [9,28,6,14,29],thesharpinterfacemethod[34,39],theghost”uidmethod[7,32,10],PHYSALIS[25,31],Cal-hounsmethodmethodandRussellandWangsmethodmethod.AnotherclassofCartesiangridmethodsaremoresuitabletosolve”owswithboundariesmovedanddeformedbydrivingforce,mostnotably,Peskinsimmersedboundarymethod[21,22]andthemorerecentimmersedinterfacemethod[17…19,16].Ourgoalofthispaperistodeveloptheimmersedinterfacemethodforsimulating”owwithbothtypesofmovingboundariesinanaccurateande
cientway.Inhissimulationofblood”owintheheart,Peskinintroducedtheimmersedboundarymethod[21,22].Inthismethod,theboundaryofanimmersedobjectistreatedasasetofLagrangian”uidparticles.Thecon“g-urationoftheseLagrangian”uidparticlesdeterminesaforcedistributioninthe”uidtorepresenttheeectoftheobject.Theforcedistributionisdeterminedbyaprescribedconstitutivelawofmaterial,anditappearsintheformoftheDiracdeltafunction.Inthissense,Peskinsimmersedboundarymethodisamathematicalformulation,inwhichsingularforceisaddedtotheNavier…Stokesequationsformodelingimmersedbound-aries.Thisformulationnaturallycouplesa”uidwithmultiplemovingobjects,andcanbecomputedinane
-cientway.Arigorousderivationoftheformulationfromtheprincipleofleastactioncanbefoundinin.Thenumericalimplementationoftheimmersedboundarymethodemploysa“xedCartesiangridfora”uidandaLagrangiangridforanimmersedboundary.Thecommunicationbetweenthe”uidandtheimmersedboundaryisachievedthroughthespreadingofthesingularforcefromtheLagrangiangridtotheCartesiangridandtheinterpolationofthevelocityfromtheCartesiangridtotheLagrangiangrid.AdiscreteDiracdeltafunctionisusedtospreadthesingularforceandinterpolatethevelocity.TherearemultiplechoicesofthediscreteDiracdeltafunction,whichisconstructedtopreservevariousmoments.TheuseofthediscreteDiracdeltafunctionremovesthesingularityinthegoverningequationsandthereforethediscontinuitiesofthe”ow“eld.Thus,standarddiscretizationschemescanbeadoptedwithoutmodi“cations.Tocon“nethethick-nessoftheinterfacebetweena”uidandanobject,thediscreteDiracdeltafunctionhasanarrowsupport,andisdependentonthegridsize.Sincetheimmersedboundarymethodwasintroducedin19721972,ithasbeenwidelyusedinthesimulationof”uid…structureinteractions,especiallyinbiological”uiddynamics.Examplescanbefoundinin.Despiteitse
ciencyandrobustness,theinitialimplementationsofthemethodhadthefollowingshortcomings.First,itwasonly“rst-orderaccurateinspace;Second,therewasasystematictendencyforaclosedimmersedboundarytoslowlylosevolumeasthoughthe”uidwereleakingout;Third,asolutionwhichispiecewisesmoothacrossanimmersedboundarywassmearedout.Inthepast30years,therehavebeenvariousresearcheortstoanalyzeandimprovetheimmersedbound-arymethod.BeyerandLeVequeLeVequeexaminedtheaccuracyofthemethodfortheone-dimensionaldiusionequation,andfoundthatadditionaltermsforthediscreteapproximationoftheDiracdeltafunctionaresome-timesnecessaryinordertoachievesecond-orderaccuracy,butitisunclearhowtomaintainthesecond-orderaccuracyfor”owsimulationinhigherdimensions.Aformallysecond-orderimmersedboundarymethodwasproposedbyLaiandPeskinin,butthesecond-orderaccuracyisachievedonlyiftheDiracdeltafunctionisapproximatedbyagrid-independent“xedsmoothfunction.Inpractice,itisstill“rst-orderaccurate.RealizingthattheprojectionofadiscreteDiracdeltafunctionontoadivergence-freespacemaybecomputedanalyt-ically,CortezandMinionMiniondevisedtheblobprojectionimmersedboundarymethod,whichdisplayedfor-mallyfourth-orderconvergenceratesoftheirbackground”owsolver.However,theanalyticalformoftheprojectiondependsonthevelocityboundaryconditionsimposedonthecomputationaldomain.Sincethepressurewasnotreported,itisalsounclearhowaccuratelythepressurecanberecovered.Analyzingtheleakageproblemintheimmersedboundarymethod,PeskinandPrintzPrintzfoundthatthevolumelossenclosedbyanimmersedboundaryiscausedbytheviolationofthedivergence-freeconditionatLagrangian”uidparticles.Theyruledouttheobviousexplanationthat”uidescapesbetweendiscreteLagrangian”uidparticleswhichde“netheimmersedboundary,andintroducedarecipefortheconstructionofa“nite-dierencedivergenceoperatortoachievebettervolumeconservation.TheblobprojectionimmersedboundarymethodbyCortezandMinionMinionalsoachievedgoodvolumeconservation.OtheranalysisandimprovementoftheimmersedboundarymethodincludethestabilityanalysisbyTuandandandStockieandWettonn,theadaptiveversionbyRomaetal.al.,andtheinclusionofbound-arymassbyZhuandPeskinPeskin.However,despiterecentimprovements,themethodstillsuersfromsomeshortcomings,inparticular“rst-orderaccuracyingeneral.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Motivatedbythegoaltoeventuallyobtainsecond-orderaccuracyinPeskinsimmersedboundarymethod,LeVequeandLiidevelopedtheimmersedinterfacemethod.Thekeyideaoftheimmersedinterfacemethod,whichisalsothemaindierencefromtheimmersedboundarymethod,istoincorporatethejumpconditionscausedbytheDiracdeltafunctioninto“nitedierenceschemes.TheapproximationoftheDiracdeltafunctionbyasmoothfunctionisthereforeavoided.TheimmersedinterfacemethodwasoriginallyproposedforellipticequationsuationsandStokesequationsuations.Later,itwasextendedtoone-dimensionalnon-linearparabolicequationsuations,PoissonequationswithNeumannboundaryconditionsitions,andthetwo-dimensionalincompressibleNavier…Stokesequations[19,16]For”uiddynamicsproblems,theimmersedinterfacemethodisbasedonthesamemathematicalformula-tionasPeskinsimmersedboundarymethod,butthecouplingbetweena”uidandanobjectisnowhandledbyincorporatingjumpconditions.Ifnecessaryjumpconditionsareallknown,theimmersedinterfacemethodcanachievesecond-orderorevenhigher-orderaccuracy.Thesharpnessofaninterfacecomputedbythemethoddoesnotdependongridresolutions.Furthermore,themethodshowsverygoodconservationofthemassenclosedbyano-penetrationboundary.Thus,theimmersedinterfacemethodsharestheadvantagesoftheimmersedboundarymethodwiththesamemathematicalformulation,whileovercomingsomeofitsshortcomingsbyeliminatingtheuseofdiscreteDiracdeltafunctions.Theapplicabilityoftheimmersedinterfacemethoddependsonwhetherthenecessaryjumpconditionsareknown.Recently,wesystematicallyderivedthejumpconditionsofall“rst-,second-andthird-orderspatialderivativesofthevelocityandthepressureaswellas“rst-andsecond-ordertemporalderivativesofthevelocityforthe3DincompressibleNavier…Stokesequationssubjecttosingularforceforce.Usingthesejumpconditions,theimmersedinterfacemethodcanbeappliedtothesimulationof3Dincompress-ibleviscous”owswithlocal“rst-orsecond-orderspatialandtemporaldiscretizationaccuracy.Inthispaper,weobtainthejumpconditionsin2Dfromour3Dresultstsbytakingonedirectionin3Dasuni-form,andimplementtheimmersedinterfacemethodtosimulatetheinteractionofa”uidwithmovingboundaries.Themethodwedevelopinthispapercansimulatetwoclassesof”owproblems,onewithboundariesmovedanddeformedbydrivingforceandtheotherwithboundariesprescribedwithknownmotion.LiandLaiLaiandLeeandLeVequeuehaveimplementedtheimmersedinterfacemethodforsometwo-dimensional”owsandachievedsecond-orderspatialaccuracyintheirsimulations.Thecurrentpaperdiersfromtheirworkinthefollowingaspects:(1)wederivethejumpconditionsinthispaperfromourthree-dimen-sionaltheoreticalresultsts,sothenumericaltestsinthecurrentpaperserveinparttoverifyourprevioustheoreticalderivation;(2)inadditiontomorespatialjumpconditions,wealsoderivetemporaljumpcondi-tionsandexaminetheireectontemporalintegration;(3)wesimulate”owproblemswithboundariesmovedanddeformedbydrivingforceandwithboundariesinprescribedmotion,andinvestigatethespatialandtem-poralconvergenceandaccuracyofourmethodagainstknownanalytical”owsolutionsandcanonical”owexamples;(4)weinvestigatethee
ciencyandrobustnessofthemethodtosimulate”owwithmultiplemovingboundaries.Thispaperiswrittenwithsu
cientdetailssothataninterestedreadercanprogramandtestthemethod.InSection,wepresentthemathematicalformulationofgoverningequationsintheimmersedinterfacemethod.InSection,wedescribetheimmersedinterfacemethodandgive“nitedierenceschemeswithjumpcondi-tionsincorporated.InSection,wepresentthenecessaryjumpconditionsforachieving“rst-orderlocalspa-tialandtemporaldiscretizationaccuracyin2Dsimulation.Thosejumpconditionsarefunctionsofthesingularforce.InSection,wediscusstheconstructionofthesingularforce.InSection,weimplementtheMACschemeschemeasabasic”owsolverandalsogiveinterpolationschemesand”uidforcecalculations.InSection,weprovidenumericalteststoexaminetheaccuracy,e
ciencyandrobustnessoftheimmersedinterfacemethod.Sectionisconclusions.2.GoverningequationssubjecttosingularforceBoththeimmersedinterfacemethodandtheimmersedboundarymethodmodelobjectsassingularforceintheNavier…Stokesequations.Thenondimensional2DNavier…StokesequationssubjecttosingularforceS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Þ¼r istime,isthevelocity,isthepressure,istheReynoldsnumber,isthenumberofobjects,andisthesingularforcefromobject.Takingthedivergenceofthemomentumequationtheequationforpressure Dot2vD ReDDþ2 ouox ovoy ouoy ¼r
isthedivergenceofthevelocityandistheCartesiancoordinates.WekeeptermswithdivergenceinEq.toenforcethedivergence-freeconditionintheMACscheme,whichwewillusetosolveEqs.(1)and(3)Sinceournumericaltreatmentofeachobjectisthesame,weomitsubscriptassociatedwithobjectinlaterpresentation.Weconsiderthesituationinwhichthesingularforceisappliedtotheclosedsmoothboundaryof.ReferringtoFig.1,theclosedboundaryisaclosedcurvein2D,andwedenoteitbyanditscoor-dinatesby.WeparametrizecurveatareferencetimebynondimensionalizedLagrangianparam-.Thesingularforcecanbeexpressedbyisthenondimensionalsingularforcedensity,and)isthenondimensionalDiracdeltafunction.Weassumearesmoothfunctionsof3.Finitedi
erenceswithjumpcontributionsBecauseofthesingularforceappliedoncurve,a”ow“eldgovernedbyEqs.(1)and(2)isgenerallynotsmoothacrossthecurve.Theimmersedinterfacemethoddiersfromtheimmersedboundarymethodinitshandlingoftheforcesingularity.InsteadofapproximatingtheDiracdeltafunctionbysmoothfunctions,theimmersedinterfacemethodmodi“esthestandarddiscretizationschemestoincludethediscontinuitiesofthe”ow“eld.Themodi“cationisbasedonthegeneralizedTaylorexpansion.Lemma1(GeneralizedTaylorexpansion).Assumefunctiong(z)hasdiscontinuitypointsoftheÞrstkindatin(z),z,andandz1;z2Þ[


[ð.g(z)canbeeithercontinuousordiscontinuousatzandz.LetÞ¼...;.Then gðnÞðzþ0Þn!ðzmþ1z0ÞnþXml¼1X1n¼0 ½gðnÞðzlÞn!ðzmþ1zlÞn.ð5Þ x
(X, Y)Š+Lagrangian poin t Fig.1.Geometricdescriptionofanobjectsurface.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Theproofoftheabovelemmawaspresentedinin.Weapplytheabovelemmatoconstructcentral“nitedierenceschemeswherethediscontinuitymayoccuratmostoncebetweentwoadjacentgridpoints.Lemma2(Generalizedcentral“nitedierences).Letx0andx.Supposeu(x)isinÞnitelysmoothexceptatdiscontinuitypointsoftheÞrstkind,.u(x)canbeeithercontinuousordiscontinuousatxandx.Then duðxiÞdx¼ uðxiþ1ðxþi1Þ2hþ 12hX2n¼0 ½uðnÞðnÞn!ðxi1nÞnX2n¼0 ½uðnÞðnÞn!ðxiþ1gÞnOðh2Þ;ð6Þ d2uðxiÞdx2¼ uðxiþ1uðxiðxþi1Þh2 1h2X3n¼0 ½uðnÞðnÞn!ðxi1nÞnþX3n¼0 Other“nitedierenceschemeswithdierentorderscanalsobeconstructedbasedonLemma1,buttheavail-abilityofjumpconditionssetsanupper-limitontheorderofaccuracy,asstatedinthefollowingproposition.Proposition1.ThehighestorderofaÞnitedifferenceschemeforu(x)withastencilcontainingadiscontinuityismn+1,wheremtakesamaximumnumbersuchthatjumpconditions[u)](l=0,1,2,,m)areallknown.Thehighestorderofspatialderivativesinmomentumequationandpressureequationis2.Accord-ingtoProposition1,todiscretizethesederivativeswith“rst-orderaccuracyatgridpointsnearanobjectboundary,thejumpconditionsofthevelocity,thepressureandtheir“rst-orderandsecond-orderspatialderivativesareneeded.Toachievesecond-orderaccuracy,thejumpconditionsoftheirthird-orderspatialderivativesarealsoneeded.Ifanobjectismoving,thetemporalderivativesofthevelocityatagridpointmayhavejumpswhenevertheobjectboundarycrossesthatgridpoint.Supposeaboundarypassesagridpointattimebetweentime.Thentherelationforvelocityatthegridpointbetweentimeisgiven on~vðt0Þotn ðtmþ1t0Þnn!þXml¼1X1n¼0 on~vðtlÞot where[[]]denotesatemporaljumpattime¼ð
Þð
Þ.Eq.followsdirectlyfromLemma1.Thehighestorderoftemporalderivativesinmomentumequationis1.AccordingtoProposition1,todiscretizethesederivativeswith“rst-orderaccuracyatgridpointscrossedbyanobjectboundary,thejumpconditionsofthe“rst-ordertemporalderivativesofthevelocityareneeded.Toachievesecond-orderaccuracy,thejumpconditionsofsecond-ordertemporalderivativesofthevelocityareneededaswell.4.SpatialandtemporaljumpconditionsInthissection,wegivethenecessaryjumpconditionsforachieving“rst-orderlocalspatialandtemporaldiscretizationaccuracyinthe2DNavier…Stokesequations.Theyarederivedfromthe3Dresultsininbytakingonedirectionin3Dasuniform.ForcurveFig.1,unittangentialvectorandunitnormalvectorarede“nedas J oXoa; where oXoaÞ2 .Asbefore,weuse[]todenoteaspatialjump,i.e.¼ð
Þð
Þ,whereisatthesideofinthedirectionofnormal,andisattheothersideof.IntheCartesiancoordinatesystem,tangentialforcedensityandnormalforcedensityarede“nedasS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS fs¼ fxsxþfysyJ;fn¼ 4.1.SpatialjumpconditionsThespatialjumpconditionsusedfor2Dsimulationare o~voxsy~sfs; o~voyResx~sfs; o2~vox2~ru1ðs2xs2yru2ð2sxsyru3ðs2yÞ; o2~voy2~ru1ðs2ys2xru2ð2sxsyru3ðs2xÞ; ½¼ opox sxJ ofnoaþ syJ ofsoa; opoy syJ ofnoa sxJ ofsoa; o2pox2rp1ðs2xs2yp2ð2sxsyp3ðs2yÞ; o2poy2rp1ðs2ys2xp2ð2sxsyp3ðs2xÞ; o2poxoyrp1ð2sxsyp2ðs2ys2xp3ðsxsyÞ;where~ru1,~ru2,~ru3,rp1,rp2andrp3are~ru1 J2 oJsxoa o~vox oJsyoa o~voy~ru2 JRe ofs~soaþ onxoa o~vox onyoa o~voy~ru3¼Re½rp;rp1¼ 1J2 o2fnoa2 oJsxoa opox oJsyoa opoyrp2¼ 1J ooa 1J ofsoa onxoa opox onyoa opoyrp3¼2 ouox ovoy2 ouoy Termswiththeformof[]appearhereinandalsolater.Notethat[[a]Æ[b].Thecorrectcalcu-lationof[]is
½
½canbeobtainedthroughinterpolation.WegivetheinterpolationschemeinSection4.2.TemporaljumpconditionsWhencurvepassesa“xedpointinspaceattime,usingtodenotethepointonwhichcoincideswiththepoint,wehavethefollowingrelationbetweenÞ¼ðÞ¼ðS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS ½¼0,wecanapproximatetemporalderivativesatbythoseatwith[[]]=0.Accord-ingtoEq.,weneedspatialjumpconditions inour2Dsimulation.Thelatterisgivenby
½r5.ConstructionofsingularforceThejumpconditionsgiveninSectionarefunctionsofthesingularforcedensity.Theconstructionofthesingularforcedependsonwhetherthemotionofanimmersedboundaryisprescribedordrivenbyaforcelawbasedonitsdeformation.WeconsiderthegeneralcaseinwhichtheimmersedboundaryhasLagrangianmassdensity).Wecanwritethedensity)ofthewholesystemaswhere)isthe”uiddensity.NotethattheDiracdeltafunctionisnondimensional.Takinganin“ni-tesimalsegmentofcurveasshowninFig.2,weapplyonittheNewtonssecondlawnondimensionalizedbythesamescalesasthoseusedtonondimensionalizeEqs.(1)and(2)toobtain qsqf whereisthevelocityofthesegment,istheresultant”uidforceactingonthesegment,andistheforcegeneratedbytheobjectitself,whichcanbemuscleforceorcontrolforce.Wesimplycallnon”uidforce.The”uidforceisrelatedtothesingularforcedensitybyIfanobjectisdeformableandiftheforcelawrelatesthemechanicforceandthedeformationoftheobjectisknown,thesingularforcecanbereadilyobtained.Therelaxationofastretchedmembraneina”uidisatyp-icalexample.Ifthemembraneiselastic,wehave wheretensionisgivenby inwhichisaconstantandfortheunstretchedmembrane. Fig.2.SurfacesegmentofvelocitysubjecttoforcesS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Ifthemotionofanobjectisprescribed,aninverseproblemneedstobesolvedtoobtainasingularforcedensitydistributionsuchthattheresultingmotionoftheobjectmatchestheprescribedone.Here,wecon-structthesingularforcedensitydistributionbasedonafeedbackcontrol.Welet,givenby qmqf areconstants,arethesimulatedandprescribedvelocityoftheboundaryseg-ment,andanditssimulatedandprescribedcoordinates.Wemayalsocombineasolidmodelwithfeed-backcontrol,i.e.Whenisfromonly,afterpluggingEq.(13)and(16)intoEq.,weobtain dD~VdtþKdD~VþKsZD~Vdt¼~fþ~fs;ð17ÞwhereKm¼ qsþqmqf,D~V¼~Ve~Vand~fs¼ qsqf .ThuswehaveafeedbackcontrolsystemasgovernedbyandsketchedinFig.3.ThisfeedbackcontrolsystemhasalinearPID(ProportionalIntegralDeriv-ative)controllerandanonlinearplantoperatedthroughtheNavier…Stokesequations.Ifcurveismassless=0),andwelet=0,thelinearcontrollerbecomesaPI(ProportionalIntegral)controllerwithinEq.ThenonlinearplantmakestheanalyticaldesignofthePIDorPIcontrollerverydi
cult.TotunethePIDorPIcontrollerforreducingtheerror,,andmaintainingnumericalstability,wecurrentlyresorttoatrialanderrorapproach.Ourexperienceindicatesthathavetotakeverysmallpositivevaluesorzerotoensurenumericalstability.ThemainparametertobetunedinthePIDorPIcontrolleristhus.Sincethepressurejumpacrossaboundaryissubjecttoanarbitraryconstant,wetunebasedontheorderofmag-nitudeanalysisfortheviscousforce.TheviscoustermsinEq.haveorderof1intheboundarylayerofthickness alongtheboundary.FromSection,weestimate 1Re o~vo~n LetthetolerancesforinEq.,wheresubscriptdenoteatangentialcomponentofavector.Wehave 6.ImplementationoftheimmersedinterfacemethodWehavedescribedthethreemaincomponentsoftheimmersedinterfacemethod,whicharethesingularforcedensity,thejumpconditionsintermsofthesingularforcedensity,and“nitedierenceschemesincor-poratedwiththejumpconditions.Nextweimplementthemethodina”owsolverbasedontheMACscheme Fig.3.Blockdiagramofthefeedbackcontrolintheforceconstructionforanobjectinprescribedmotion.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 6.1.MACgridandboundaryrepresentationAMACgridisastaggeredCartesiangrid.WeuseauniformMACgridassketchedinFig.4(a)forthe“nite-dierenceapproximationtogoverningequations(1)…(3).Thecenterofacellforpressureis(),where{1,2,}and{1,2,}.Thecenterofacellforvelocitycomponentis()andthecenterofacellforvelocitycomponentis(),where()correspondsto 12;jþ {0,1,}and{0,1,}.ThedimensionsofacellareWecalculategeometricquantitiesandthesingularforcedensityattheLagrangianpointsnumberedbyoncurve,asmarkedbycirclesinFig.4(b),where{0,1,}.WeuseperiodiccubicsplinestointerpolatethevaluesofthegeometricquantitiesandthesingularforcedensityattheirregularpointsmarkedasXinFig.4(b),whicharetheintersectionsofcurveandgridlines.Thecostofcubic…splineinterpolationisoforder.Wecanthenidentify“nitedierencestencilsthatcontaintheirregularpoints,andcomputethejumpcontributionsforthe“nitedierenceschemesonthosestencils.TomoveLagrangianpoints(),we“rstinterpolatethevelocityofalltheirregularpointswithaninterpolationschemegiveninSection,andthenuseperiodiccubicsplinestointerpolatevelocity(attheLagrangianpoints.Whentwoirregularpointsareveryclosetoeachother,weuseoneoftheminordertoavoidfailureofthecubicsplines.6.2.SpatialdiscretizationWiththecellde“nitionsforandshowninFig.5,wecanwritethesecond-ordercentral“nitedierenceapproximationstoEqs.(1)and(3)inthefollowingform: u(I,j)v(i,J)p(i,j)Š1)u(IŠ1,j)i,I j ,Jxy Fig.4.Discreterepresentationofthecomputationaldomainandtheobjectsurface.(a)TheMACgridwithstaggeredarrangementof”owvariables;(b)Lagrangianpoints(opencircle)andirregularpoints(x-mark)ontheobjectsurface. Fig.5.Cellde“nitionsfor(a)velocitycomponent,(b)velocitycomponent,and(c)pressureandvelocitydivergenceS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS ouot;j¼CuI;j u2iþ1;ju2i;jDxþ uI;JvI;JuI;J1vI;J1Dy piþ1;jpi;jDxþ 1Re uIþ1;j2uI;jþuI1;jDx2þ uI;jþ12uI;jþuI;j1Dy2ð18Þ ovot;J¼Cvi;J uI;JvI;JuI1;JvI1;JDxþ v2i;jþ1v2i;jDy pi;jþ1pi;jDyþ 1Re viþ1;J2vi;Jþvi1;JDx2þ vi;Jþ12vi;Jþvi;J1Dy2ð19Þ oDot;j¼Cpi;jþ piþ1;j2pi;jþpi1;jDx2þ pi;jþ12pi;jþpi;j1Dy22 uI;jDI;juI1;jDI1;jDxþ vi;JDi;Jvi;J1Di;J1Dy 1Re Diþ1;j2Di;jþDi1;jDx2þ Di;jþ12Di;jþDi;j1Dy22 uI;juI1;jDx vi;Jvi;J1Dy ui;Jui;J1Dy arethetermsduetojumpcontributionsin“nitedierences,andiscomputedas uI;juI1;jDxþ withthejumpcontributiontermdenotedby.SomeofthevelocityanddivergencevaluesinEqs.arenotcenteredintheircellsshowninthecelldiagraminFig.5.Weinterpolatethemfromadjacentvalues.Again,theinterpolationschemesarepresentedlaterinSection6.4isnonzeroonlyifthestencilofa“nitedierenceinthecorrespondingequationcontainsirregularpoints.Togiveanexample,wecalculateforthecaseshowninFig.6(a).Wedenotethecoordi-natesofthetwoirregularpointsas()and(),respectively.Ajumpconditionatpointnowde“nedas½
¼ð
Þð
Þ,andajumpconditionatpoint½
¼ð
Þð
Þisthesumofthejumpcontributionfromeach“nitedierenceinEq..Inthiscase,thejumpcontributionassociatedwith u2iþ1;ju2i;jDxis 1Dx ou2oxxiþ1X 2 thejumpcontributionassociatedwith uI;JvI;JuI;J1vI;J1Dyis Fig.6.Schematicsofexamplesfor(a)thecalculationofjumpcontributiontermand(b)thetreatmentofboundaryconditionsonafar-“eldrigidwall.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 1Dy ouvoyyJY 2 thejumpcontributionassociatedwith piþ1;jpi;jDxis 1Dx½p poxxiþ1X 2 thejumpcontributionassociatedwith 1 Iþ1;j2uI;jþuI1;jDx2is 1ReDx2 ouoxxIþ1X 2 thejumpcontributionassociatedwith 1 I;jþ12uI;jþuI;j1Dy2is 1ReDy2 ouoyyjþ1Y 2 (22)and(23)containtermsoftheformof[],forexample ou2ox½u ,whichcanbecalculatedaccordingtoEq.InournumericaltestsinSection,wesometimesplotthestreamfunctionofavelocity“eld,whichisobtainedbysolving o2wox2þ whereisthestreamfunctionandisthevorticity.Byde“nition,wehave woy;v¼ owox;ð28Þx¼ ovox Thuswecanderivethefollowingjumpconditions: owox0; owoy0; o2wox2 ovox o2woy2 Thecentral“nitedierenceapproximationstoEq.andde“nition wiþ1;j2wi;jþwi1;jDx2þ wi;jþ12wi;jþwi;j1Dy2þCwi;j¼xi;j;xi;j¼ vI;jvI1;jDx wherearefromthejumpcontributionsin“nitedierences.Theyarecomputedinthesimilarwayasthecalculationofgiveninthepreviousexample.Wenowdiscussourimplementationoftheno-slip,no-penetrationandthepressureboundaryconditionsonafar-“eldrigidwall.ReferringtoFig.6(b),weenforcetheno-slipandno-penetrationboundaryconditionsasfollows: uðI;0ðI;2Þ2¼0; Withthetreatmentof aspresentedinSection,Eq.forthepressureisadiscretePoissonequa-tion.WecurrentlysolvethediscretepressurePoissonequationandstreamfunctionPoissonequationusingFFT,whichhascostoforder,whereisthenumberofnodesinaCartesiangrid.ReferringFig.6(b),wederivefromEq.thefollowingNeumannboundaryconditionforthepressureatarigidS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS opoy;j¼1¼ 1Re iscomputedbasedonthedivergencefreeconditionatgridnode(=0)asfollowing: uI;J¼0uI1;J¼0Dxþ 6.3.TemporalintegrationWeemployanexplicitfourth-orderRunge…Kuttaschemeinthetemporalintegration.High-orderexplicitschemesareappropriateforthemoderateReynoldsnumbersthatweconsiderhere,asshownbyJohnstonandandandEandLiuLiu.Inthe”owregimeofmoderatetohighReynoldsnumbers,viscoustimestepconstraintismuchlessrestrictivethantheconvectiveone.Ahigh-ordertemporalintegrationschemewhosestabilityregionincludesaportionoftheimaginaryaxiscanensurestability,andanimplicittreatmentofvis-coustermsisnotnecessary.Inthissection,wefocusonhowtoincorporatetemporaljumpconditionsintothefourth-orderRunge…Kuttascheme.Duetolimitedtemporaljumpconditions,thetemporalaccuracyisonly“rst-orderforagridpointattheinstantwhenthegridpointiscrossedbyanimmersedboundary,buttheoveralltemporalaccuracyisnotaectedmuchsincethenumberofsuchgridpointsismuchlessthanthetotalnumberofgridpoints.Wede“nevectorsandaretheright-handsidesofEqs.(18)and(19),respectively.Thenwehavetheordinarydierentialequationasfollowing: supplementedwithEq.,whichisrewrittenbelowwithsubscriptsomitted. Withsuperscriptdenotingatimelevelandatimestep,thesequenceoftemporalintegrationofEq.usingtheforth-orderRunge…Kuttaschemeisasfollows:(1)solvepressure andthencompute 121: 0Dn Dt2Dpnþ 121þRpðqnÞ;qnþ 121¼qnþ Dt2Rðqn;pnþ (2)solvepressure andthencompute 122: 0Dn Dt2Dpnþ 122þRpðqnþ 121Þ;qnþ 122¼qnþ Dt2Rðqnþ 121;pnþ (3)solvepressureandthencompute 0DnDt¼Dpnþ13þRpðqnþ 122Þ;qnþ13¼qnþDtRðqnþ S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS (4)solvepressureandthencompute 0DnDt¼Dpnþ14þRpðqn3Þ;qnþ1¼qnþDt 16Rðqn;pnþ 121Rðqnþ 121;pnþ 122Rðqnþ wherearejumpcontributionsforthetemporaldiscretizationof inRunge…Kuttasub-steps.Theyarenonzeroonlyatthosegridnodeswhicharepassedbyaboundary.Theyarecomputedasfollows:treatRunge…Kuttasub-step(1)asaforward“nitedierenceininterval ,andcalculateforagridnodeifthecurvepassesthegridnodeattimebetweentimeand 12:cn1¼ þ oqotnþ treatRunge…Kuttasub-step(2)asabackward“nitedierenceininterval ,andcalculateforagridnodeifthecurvepassesthegridnodeattimebetweentimeand 12:cn2¼ þ treatRunge…Kuttasub-step(3)asacentral“nitedierenceininterval,andcalculateforagridnodeifthecurvepassesthegridnodeattimebetweentimeand 12:cn3¼ þ orifthecurvepassesthegridnodeattimebetweentime 12andtncn3¼ þ treatRunge…Kuttasub-step(4)asacombinationofoneforward,onebackward,andtwocentral“nitedif-ferencesininterval,andcalculateforagridnodeifthecurvepassesthegridnodeattimebetween 12:cn4¼ þ oqottnþ1tÞ þ oqottnþ1tÞ þ orifthecurvepassesthegridnodeattimebetweentime 12andtncn4¼ þ oqottnþ1tÞ þ oqottntÞ þ andajumpcondition[[]]attime,whichisde“neas½½
¼ð
Þð
Þ,areobtainedbyinterpolationusingknowninformationattimelevel +1.WepresenttheinterpolationschemesinSection6.4.FilteringandinterpolationForamovingobjectproblem,Lagrangianpointsontheobjectboundaryaremovedusingthevelocityinterpolatedfromthesurrounding”uidvelocityonCartesiangridpoints.Theboundaryshapethuscontainstheirregulartruncationerrorsoftheinterpolation.Sincethejumpconditionsinvolvethedierentiationofgeometricquantitiesandthesingularforcedensityalongaboundary,itisimportanttomaintaintheshapesmoothnessoftheobject.WeemployFourier“lteringtosmooththeshape.Itisnecessarytointerpolatequantitiessuchas,[[)andthevelocityatirregularpoints.Wecategorizetheinterpolationsintothreescenarios,asshowninFig.7S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Todeterminetimewhenaninterfacecrossesagridpoint,wemonitorthepositionofairregularpointalongagridline,i.e.intheinsetofFig.7(a).Sincethetimeisasmoothfunctionofthechangingcoordinateoftheirregularpoint,theinterpolateschemeisgivenby istimeandisthechangingcoordinateinFig.7(a)forthecurrentcase.Itcanalsobeshownthatjumpcondition[[)isasmoothfunctionoftimealongthegridline,andthereforetheaboveinterpolationschemeappliesto[[Thespatialinterpolationfor)andthevelocityatirregularpointsisdescribedinFig.7Theinterpolationschemeis hgCþhþgAh h½gBh hþhh2 ogBozOðh2Þ;ð38Þg¼ hgCþhþgAhþ hþ½gBh hþhh where[]=0andforthevelocityinterpolationsincethevelocityiscontinuous.TheinterpolationforisdescribedinFig.7(c).WhenpointDfallsbetweenpointsAandBasinFig.7(c),wehave gAþgC2þ 12½gD 2 ogDozþ 14 WhenpointDfallsbetweenpointsBandC,wehave gAþgC2 12½gD 2 ogDozþ 14 6.5.ForcecalculationWehereuserespectivelytodenotethecomponentsofthenondimensionalforceappliedbythe”uidtotheobject.TheyarenondimensionalizedbyhalfofthepressurescaleinEq..ReferringtoFig.1,theycanbecalculatedasfollows: h Fig.7.Interpolationscenariosassociatedwith(a)asmoothfunction,(b)adiscontinuousfunctionatitsdiscontinuitypoint,and(c)apiecewisesmoothfunction.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Cx¼2ICpþnxþ 1Re ouonCð42ÞICfxdaþ2ICpnxþ 1Re ouonC;ð43ÞCy¼2ICpþnyþ 1Re ovonCð44ÞICfydaþ2ICpnyþ 1Re Foracentrallysymmetricboundary,decomposingitsprescribedmotiontothesuperpositionofatranslationandarotationaroundthegeometriccenter,we“nd duTdt;ICðpnyÞdC¼S dvTdt;IC ouondC¼0;IC where()arethetranslationalvelocityandistheareaoftheobject.SoinEqs.(43)and(45)canbesimpli“edto duTdt;CyICfydaþ2S 7.ResultsInthissection,wetestthemethodinseveraldistinct”uid…structureinteractions,including”owsofknowanalyticalsolutions,”owinducedbyarelaxingballoon,”owpastacylinder,a”apper,andmultiplecylinders.Inparticular,weinvestigatethespatialandtemporalconvergenceratesofthemethod,anddemonstratetherobustnessande
ciencyofthemethodinhandlingsingleormultipleboundariesprescribedwithknownmotionordrivenbyaforcelaw.TheReynoldsnumber,,ofthe”owsrangesfrom1to200.Inthesetests,objectboundariesaremassless,i.e.=0inEqs.(11)and(12),andwehave.Welet=0inEq..So=0and=0inEq.whenthefeedbackcontrolisusedtoconstructforce7.1.FlowwithoutimmersedboundariesAsa“rsttest,wesimulatea”owwithoutanyimmersedboundariestocheckthe”owsolver.Weconsiderthefollowing”owwhichsatis“estheNavier…Stokesequationsexactly:          17sin    15cos    
Wesimulatedthe”owat=1inthedomainde“nedby    .Periodicboundaryconditionswereusedinthistest.Werunthesimulationuptotime=10withtimestep=0.01anddierentspatialresolutionstocheckthespatialconvergencerate.Thein“nitynormofnumericalerrorbasedontheanalyticalsolutionisgiveninTable1,wheretheorderofaccuracyisde“nedasS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS order k

kisthein“nitynormofa“eldwitharesolutiontwicetheresolutionofa“eldassociatedwithin“nity.Clearly,second-orderspatialconvergencerateisachieved,asexpectedfromthecentral“nitedif-ferenceschemesinspace.We“x=10and=3232,butusedierenttimestepstocheckthetemporalconvergencerate.Wecomputeareferencesolutionusingasmalltimestep,=0.001,andsubtractthereferencesolutionfromtheothernumericalsolutionstocancelouttheerrorduetospatialdiscretization.Thein“nitynormofnumer-icalerrorbasedonthereferencesolutionisgiveninTable2.Only“rst-ordertemporalconvergencerateisachieved,whichislowerthantheexpectedfourth-orderconvergencerateoftheRunge…Kuttascheme.Thisisduetothenumericaltreatmentofthetemporalderivativeofthedivergence, ,inpressureequation,wherewesetdivergenceattimelevels and+1tobezero,whichdoesnotfollowtheerrorcan-cellationmechanismintheRunge…Kuttascheme.Ifwediscardterm inEq.,wecanachievenearlyfourth-ordertemporalconvergencerate,asseenTable3.Thesimulationisrunto=25.6,andthereferencesolutioniscomputedwith=0.0005=3232.However,keeping inpressureequationimprovesthedivergence-freecondition,asindicatedbythecomparisonsinTable4.InTable4,thesimulationwasrunto=10with=0.01and Table1Spatialconvergenceanalysisforthe”owsolverwithoutimmersedboundaries161.40323.472.018.662.011.03648.652.002.162.002.571282.162.005.402.006.43Thein“nitynormisbasedontheanalyticalsolution. Table2Temporalconvergenceanalysisforthe”owsolverwithoutimmersedboundaries0.013.280.026.921.083.591.086.170.041.431.057.391.041.270.082.901.021.511.032.59Temporaldivergenceterm inpressureequationisincluded.Thein“nitynormisbasedonareferencesolutioncomputedwith=0.001. Table4Eectoftemporaldivergenceterm inpressureequationonnumericalaccuracy oDot3.46·1048.66·1051.03·1051.25·108Discard oDot4.71·1042.22·1047.46·1044.13·104 Table3Temporalconvergenceanalysisforthe”owsolverwithoutimmersedboundaries0.045.840.089.714.061.193.301.520.161.754.171.914.002.550.324.084.543.084.014.27Temporaldivergenceterm inpressureequationisnotincluded.Thein“nitynormisbasedonareferencesolution.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS =3232.FromTable4,wenotethattheaccuracybasedontheanalyticalsolutionisaboutthesamewithorwithout ,whichindicatesthattheerrorisdominatedbythespatialerrorinsteadofthetemporaloneevenatthislarge=0.01.Inourlatersimulations,wehaveaboutthesamespatialresolutionasthecurrentcasebutsmaller.Sowekeep tobetterenforcethedivergence-freeconditionwithoutlosingtheaccuracy.7.2.Taylor…Couette”owWenextconsiderTaylor…Couette”owbetweentworotatingandtranslatingconcentriccylinders.Thegeometryofthe”owisshowninFig.8,with=0.5,=1,=1and1.TheReynoldsnumberbasedontheradiusandtheangularvelocityoftheoutercylinderis 10.Totestmovingbound-ariescrossingtheCartesiangrid,weallowthecenterofthecylinderstooscillateaccordingtowhereisaconstant.Theanalyticalsolutiontothe”owbetweenthetwocylindersisgivenby Br2yycÞ;v¼d
cosðtþ Br2xxcÞ;p¼d
sinðtþy 2r22 where X2r22X1r21r22r21,B¼ .Theanalyticalsolutiontothe”owinsidetheinnercylinderis Weuseperiodicboundaryconditionsatthefar-“eldboundaries.Thenon”uidforcemodelforbothcylindersisgivenby=1000. Fig.8.Geometryof”owbetweentworotatingconcentriccylinders.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 7.2.1.d=0Inthesteadystate,thenormalforcedensity,,iszeroandthetangentialforcedensity,,isaconstantattheinnercylinder.ReferringtothejumpconditionsinSection,weknowthejumpconditionsacrosstheinnercylinderforallthe“rst-orderandsecond-ordervelocityderivativesandallthesecond-orderpressurederivativesarenonzero.Thereforethesimulationofthiscaseisagoodtesttocheckthespatialconvergencerateoftheimmersedinterfacemethod.Inparticular,welookatthein“nityerrornormtocheckthesimulationaccuracylocally.Weuseaverysmalltimestep,=0.0001,toensurethetemporaldiscretizationerrorisneg-ligiblecomparedwiththespatialone.TheresultsaregiveninTable5,whichindicatesecond-orderspatialconvergencerateforboththevelocityandthepressure.AsindicatedbyLemma2,thediscretizationoftheLaplaceoperatornearthecylindersinEqs.(1)and(3)only“rst-orderaccuratebecauseonlylimitedjumpconditionsareusedandthejumpconditionsforvelocityandpressurederivativesofhigherordersarenonzero.Interestingly,westillachievesecond-orderaccuracyofthe”ow“eldevennearthecylinders.ThesamephenomenonwasnotedbyLiandLaiintheirsimulationlation.Wecalculatejumpconditionsacrosstheinnercylindersurfacebasedontheanalyticalsolutions,andcom-parethemwiththosecomputedfromtheformulagiveninSectionFig.9showsthecomparisonfor ¼ð
Þð
Þ,and ¼ð
Þð
Þ,indicatingverygoodagreement.Forthesecomparisons,thespatialresolutionofthesimulationis=647.2.2.d=0.5Inthiscase,thevelocityacrossthetwocylindersarenotsmooth,whichcausesjumpsoftemporalvelocityderivativesandthereforejumpcontributionsintemporalvelocitydiscretizationsforagridpointwhenitiscrossedbythecylinders.Thusthiscaseprovidesatestoftheeectofthetemporaljumpcontributionsonthetemporalaccuracy.We“rstrunsimulationswiththetemporaljumpcontributionsupto=10witha“xedspatialresolution,=128256.Wecomputeareferencesolutionusingaverysmalltimestep,andsubtractthereferencesolutionfromtheothernumericalsolutionstocanceloutspatialdiscretizationerror Table5SpatialconvergenceanalysisforsteadyTaylor…Couette”ow32,642.0964,1286.031.795.331.639.18128,2561.561.951.351.982.34256,5124.121.924.231.674.68Thein“nitynormisbasedontheanalytical”owsolution. 60 120 180 240 300 360 Š20 Š10 0 10 20 30 analytical []]2ux2versus 0 60 120 180 240 300 360 Š15 Š10 Š5 0 5 10 15 analytical versus Fig.9.Comparisonsbetweenanalyticalandnumericalresultsforjumpconditions.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS andobtaintemporalerror.Thein“nitynormofthetemporalerrorisshowninTable6.Wethenrunanothersetofsimulationswithoutthetemporaljumpcontributions.TheresultsareprovidedinTable7Tables6and7indicatethatthetemporalconvergencerateswithandwithoutthetemporalcontributionsareaboutthesame.Botharenearly“rst-orderinsteadoffourth-orderbecauseofthenumericaltreatmenttothetemporaldivergenceterminthepressureequation(seeSection).Inaddition,theinclusionofthetem-poraljumpcontributionshasnegligibleeectonthelocalandoverallaccuracyofthenumericalsolutionsbasedontheanalyticalsolution,asindicatedbyTables8and97.3.FlowinducedbyarelaxingballoonInthistest,wecomputethemovingboundaryproblemconsideredbyLiandLaiLai,wherea2Ddistortedpressurizedballoonimmersedinanincompressible”uidrelaxestoitscircularequilibriumshape.Theinitialvelocityandthepressurearesetzero,andtheonlydrivingforceistheballoontension.The”uid”owandtheballoonmotionarefullycoupled.Atequilibrium,thevelocityiszeroandthepressureispiecewiseconstantinsideandoutsidetheballoonwithajumpacrosstheballoon.Theinitialshapeofthedistortedballoonexpressedinthecylindricalcoordinates()iswhereandareconstants,andisaninteger.Thenormalandthetangentialforcedensitiesofthedrivingforcearewhereisaconstantandthecurvature.Atequilibrium,theradiusoftheballoon,,isandthepressurejumpacrosstheballoonis.TheareaoftheballoonisconservedanditisequaltoWe“rstrunthetestwiththesameparametervaluesasusedbyLiandLai(aftercorrectingsometypoerrorsintheirpaperpaper1).Theyare=0.5,=0.4,=5,=0.05and=1.HereReynoldsnumbercorrespondstoviscosity=1.Theinitial(=0)andtheequilibrium()balloonshapesinthecompu-tationaldomainaregiveninFig.10.Theboundariesofthecomputationaldomainarerigidwalls.Wesimulatethecaseupto=98with=64128and Theinitialinterfaceshouldbe)with=0.2.Thevalueofshouldbe1insteadof0.1 Table6TemporalconvergenceanalysisforoscillatingTaylor…Couette”ow1.065.941.158.700.672.171.454.762.147.671.501.00Temporaljumpcontributionsareincluded.Thein“nitynormisbasedonareferencesolution. Table7TemporalconvergenceanalysisforoscillatingTaylor…Couette”ow0.906.091.128.690.841.581.944.762.006.001.401.01Temporaljumpcontributionsarenotincluded.Thein“nitynormisbasedonareferencesolution.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Fig.10plotstheballoonshapes,thevorticitycontoursandthevelocityvectorsattime=7,21and35.TheballoonshapesatthesetimeinstantsareveryclosetothoseobtainedbyLiandLaiLai.Weregardtheballoonisnearequilibriumattime=98.Thepressureattime=98isplottedinFig.11.TheconservationoftheareaoftheballoonischeckedinFig.12=1and=100,whichindicatesverylittleleakage,about0.1%.Thesimulationfor=100isrunwith=0.0001andthesamespatialresolutionas=1.Fig.13plotsthetemporalvariationoftheballoonradiusat/10at=1and=100.Wecanobserveat=100thefastrelaxationoftheballoonthroughdampedoscillations,whichareabsentat=1. Table8NumericalaccuracywithandwithouttemporaljumpcontributionsWithtemporaljumpcontributions2.04Withouttemporaljumpcontributions2.04Thein“nitynormisbasedontheanalyticalsolution. Table9NumericalaccuracywithandwithouttemporaljumpcontributionsWithtemporaljumpcontributions6.94Withouttemporaljumpcontributions6.94The2-normisbasedontheanalyticalsolution. 1 Š0.5 0 0.5 1 Š1 Š0.5 0 0.5 1 t=0,t= 1 Š0.5 0 0.5 1 Š1 Š0.5 0 0.5 1 t=7 Š0.5 0 0.5 1 Š1 Š0.5 0 0.5 1 1 Š0.5 0 0.5 1 Š1 Š0.5 0 0.5 1 (a)(b)Fig.10.Numericalevolutionoftheballoonshapeandthevorticityandthevelocity“eldsfor”owinducedbytherelaxationofadistortedS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Weperformspatialconvergenceanalysisforthe”ow“eldof=10attime=2.Inthisanalysis,thetimestepis=0.0001.Wecalculatetheerrorbetweenasimulated”ow“eldandtheoneobtainedbythe“negrid=256Table10givethein“nitynormoftheerroragainstthespatialresolution.Nearsecond-orderconvergencerateforbothvelocitycomponentsisobserved.Whenthespatialconvergence 0 1 –1 0 1 –0.05 0 0.05 0.1 p(x,y) Š1 Š 0.1 0.15 analytical (x=0, Fig.11.Equilibriumpressure:(a)pressure“eld,(b)comparisonofpressuredistributionalongaxisat=0betweentheoreticalandnumericalresults. 20 40 60 80 100 0.847 0.8475 0.848 0.8485 0.849 timearea numerical 5 0.847 0.8475 0.848 0.8485 0.849 numerical =100Fig.12.Thetemporalevolutionoftheareaenclosedbytheballoon. 20 40 60 80 100 0.5 0.55 0.6 0.65 0.7 radius 5 10 15 0.45 0.5 0.55 0.6 0.65 0.7 radius=100Fig.13.ThetemporalevolutionoftheradiusofaLagrangianpointontheballoon.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS rateforthepressureiscalculated,specialcareisneededneartheinterface.Atdierentspatialresolutions,agridpointcanbeinsidetheballoonatonegridlevelandoutsideatanother,suchasthemiddlegridpointsketchedinFig.14.Wethusdiscardthistypeofpointstoexcludethepressurejumpinthecalculationofthein“nitynorm.Thesametreatmentisnotnecessaryforthevelocitysincethevelocityiscontinuousacrosstheinterface.Wedenotethemodi“edin“nitynormsforinsteadof,respectively.Theresultsoftheconvergenceanalysisbasedonthemodi“edin“nitynormaregiveninTable11,whichindicatetheconvergencerateforthepressureisaboutthesameasthevelocity.7.4.FlowpassingamovingcylinderInsimulationsusingtheimmersedinterfacemethodandtheimmersedboundarymethod,”exibleobjectsareoftenconstructedasanetworkofsprings.Springsarealsousedtotethertheobjects.Thespringconstantisgivenbythematerialpropertyanditde“nesacharacteristictimescaleassociatedwiththevibrationmodeofanobject.Toinvestigatetheeectofthistimescaleonatime-dependent”ow,wesimulate”owpassingacylinderwhichisacceleratedfromresttoauniformvelocity.ThegeometryisgiveninFig.15,withrigidwallsasfar-“eldboundaries.Thevelocityofthecylinder,,isgivenby 11þtanhð2Þtanh 4ttc2tanhð2Þ Table10Spatialconvergenceanalysisforthe”owinducedbyarelaxingballoon16,326.1432,641.492.042.9164,1283.562.077.25128,2567.332.287.80Thein“nitynormisbasedonareferencesolution. Table11Spatialconvergenceanalysisforthe”owinducedbyarelaxingballoon16,322.3332,649.661.272.221.162.6264,1282.561.916.211.841.10128,2565.682.176.063.362.15Thein“nitynormisbasedonareferencesolutionandismodi“edtoexcludethegridpointsneartheinterface. Fig.14.Schematictoshowthatagridpointisatonesideofapressurediscontinuitypointatonegridresolutionlevelandattheothersideatanotherlevel.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS whereisthecharacteristicaccelerationtime.Fig.16.Wede“neReynoldsnumbersbasedontheuniformvelocity.ThespringforcemodelisgivenbyEq.,andsketchedontheleftinFig.17.Thetem-poralresolutionofthesimulationissetbytheCFLnumberswithCFL Dt 1Dx2þ 1andCFL uxþ Wevarytobe0.1,0.2,0.4,0.8,1.6,3.2and6.4with“xed=20and=1000.AsseeninFig.18,theLagrangianpointscanfollowtheprescribedmotionwhen1.6.When=0.1,largeoscillationsareseenindraghistory,asshowninFig.19(a).ItspowerspectruminFig.19(b)hasadominantfrequency,=7.=1.6,oscillationsinthedraghistoryhasthesamefrequencybutsmalleramplitude,asseeninFig.20.Thesamefrequencyisobservedforallvaluesof.Theamplitudesoftheoscillationsdecreaseasincreases.Thiscylinderandspringsystemcanbeeectivelyapproximatedasaspring-mass…dampersystem,asdescribedontherightofFig.17.Thefrequencyofthedampedoscillationsisgivenby (0,0)xy1WNSERe Fig.15.Geometryfor”owpassingastationarycylinder. Fig.17.Schematicshowingthephysicalinterpretationofthenon”uidforcemodel. velocity Fig.16.Schematicshowingcylindervelocityversustime.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS fs¼  comesfromtheupperlimitofnondimensionalLagrangianparameter)istheeectivespringstiness, isthe”uidmass,and isthedampingratio,whereisthedampingcoe
cientofthedamper.AsshowninFig.21(a),theoscillationfrequencyisproportionalto,asexpectedforalin-earspring.De“ne.FromEq.,wehave 2pMc1f2¼ ,whichisonlyafunctionofthedamping.WenowlookattherelationbetweenandReynoldsnumber,asplottedinFig.21(b),where 0.5 1 1.5 0 0.5 1 1.5 timevelocitytc= 0.1 0 1 2 3 0 Fig.18.VelocityofLagrangianpoints…solidlines:simulated(coveredbyopencirclesandinvisiblein(b)),dashedlineswithopencircles:prescribed. 1 2 3 Š150 Š100 Š50 0 50 100 timedrag 0 20 30 40 50 x 104 Fig.19.Dragassociatedwith=0.1:(a)timesignaland(b)frequencycontent. 2 4 6 Š8 Š6 Š4 Š2 0 2 timedrag 0 10 20 30 40 0 5 10 15 20 25 Fig.20.Dragassociatedwith=1.6:(a)timesignaland(b)frequencycontent.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS =20,40,50,100,150,200with“xed=0.2andthecorresponding=1000,500,400,200,150,100.WeobservelittledependenceofontheReynoldsnumber,,intherange.Insummary,we“ndthataslongasthetypicaltimescaleofthe”owisabout10timeslargerthanthechar-acteristictimescaleofthespringsystem,theLagrangianpointscanfollowtheprescribedmotion.Wewillusethisasourruleofthumbforsimulating”owpastacylinderanda”appingwinginthefollowingsections.7.5.FlowpassingastationarycylinderFlowpassingastationarycylinderisacanonicalexampletotestnumericalmethods.WeusethegeometryshowninFig.15forthetestofourimmersedinterfacemethod.Initially,the”owissettobeuniformwith=1andthecylindermoveswiththesamevelocityasthe”ow.Thefar-“eldboundaryconditionsusedforthetestare=1,=0and opox¼ 1 8(boundaryW); ouox¼0, 0and pox¼ 1 =24(boundaryE); =0and opoy¼ 1 8(boundaryS); =0and opoy¼ 1 =8(boundaryN).ThediscreteformsoftheaboveboundaryconditionsaresimilartothosedescribedinEqs.(30)and(31).Inparticular,thedivergence-freeconditionisincorporatedinthediscretepressureboundaryconditions.Wetakethefar-upstreamvelocityandthecylinderdiameterasthevelocityandthelengthscales.Wecomputethe”owat“veReynoldsnumbers,=20,40,50,100and200.Thespatialresolutionforthesecomputationsis=640128.ThetemporalresolutionissetbyCFL=CFL=0.2.Non-”uidforceforthecylinderisacombinationofafeedbackcontrolandaspring-supportedmem-brane,asshownschematicallyinFig.22.Thus,wehavesolid,withandgivenby ooa JJe1Es1 whereandareconstants,istheradius,andisthedesiredradius,i.e.=0.5.Table12liststhevaluesofforceparametersfortheconsideredReynoldsnumbers.Thedominantparameterinthesetestsis.FortheReynoldsnumbersconsideredhere,wecandeterminefromthefollowingempiricalformula: 20000Re. 500 1000 1500 2000 0 20 40 60 80 100 Ksfs2 100 150 200 250 20 30 40 50 (a)(b)Fig.21.(a)Thefrequency,,versusthespringstiness,;(b)theconstant(equivalenttothedampingratio),,versustheReynoldsS.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS Weremarkthatthecombinationofthesolidmodelandthefeedbackcontrolinthesetestsistodemon-stratethe”exibilityoftheforceconstruction.Incasesof=20,50and200,weonlyusethefeedback Lagrangian point(X, V)Fig.22.Schematicoftheforceconstructionforastationarycylinderin”ow. Table12Parametervaluesintheforceconstructionforastationarycylinderin”owatdierentReynoldsnumbers=20=40=50=100=200040010004002000.10.10001000400500160100 2 4 6 8 10 –100 –50 0 50 dragversus 40 60 80 100 2.22 2.24 2.26 2.28 dragversus 5 10 15 20 0.05 0.1 versus 20 40 60 80 100 –5 –4.8 –4.6 –4.4 –4.2x 10 versus(a)(b)(c)(d)Fig.23.Dragandliftcoe
cientsversustimefor”owpassingastationarycylinderat=20:(a)draghistoryintransient;(b)draghistoryfromtransienttosteadystate;(c)lifthistoryintransient;(d)lifthistoryfromtransienttosteadystate.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 7.5.1.Re=20,40and50Fig.23showsthedragandthelifthistoryfor=20.Thelargeoscillationsatthestartareduetotheimpulsivestopofthecylinder,whichareanalyzedintheprevioussection.The”owcanbeconsideredsteady=100,longafterthetransient.Thenonzeroliftoforder10isascribedtonumericalerror.Between=40and=50,the”owbecomesunstableanddevelopsavonKarmanwake.Figs.24andshowthedragandthelifthistoryfor=40and=50.At=40theoscillationsinliftdampout,and=50theyamplify.Figs.26…28showthe”owdetailsaroundthecylinderfor=20at=100and=40at=120.Table13compareslengthofthetrailingbubble,angleofseparation,anddragcoe
cientpreviousexperimentalandcomputationalresultssummarizedinin.Geometricquantitiescom-parefavorablywithothers,thoughthevalueofdragcoe
cientisslightlyhigherthanpreviousstudies.RussellandWangWangsuggestedthattheirsimpli“cationofthefar-“eldboundaryconditionsisasourceofincreaseddrag.Webelievethereasonforthedragincreaseinourcaseisthesame.Totestthis,wehavechan-gedthedomainsizefor=20from3216to4824withthecorrespondingchangeof320to960480.ThenewdragandlifthistoryisshowninFig.29.Atthesteadystate,thedragcoef-“cientsettlesdownto2.14,whichisintherangeofpreviouslyreportedvalues.The”owinsidethecylinderisstaticinourcase.Thusthevorticityandthepressuredistributionsoverthecylindersurface,denotedasrespectively,canbecalculatedasfollows: ovox ouoy¼ ovox ouoy
fs;psn; 100 120 1.654 1.656 1.658 dragversus 40 60 80 100 120 –2 –1.5 –1 –0.5x 10 versus(a)(b)Fig.24.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=40. 60 80 100 120 140 1.518 1.52 1.522 1.524 1.526 1.528 dragversus 40 60 80 100 120 140 –5 0 5x 10 versus(a)(b)Fig.25.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=50.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS –1 0 1 2 3 4 –1 –0.5 0 0.5 1 0 1 2 3 4 –1 –0.5 0 0.5 1 Fig.26.Streamlinesof”owpassingastationarycylinderatthesteadystate. 0 5 10 15 20 –5 0 5 Re = 20 0 5 10 15 20 –5 0 5 Fig.27.Vorticitycontoursof”owpassingastationarycylinderatthesteadystate.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 0 5 15 20 5 0 5 5 Fig.28.Pressurecontoursof”owpassingastationarycylinderatthesteadystate. Table13Summaryof”owcharacteristicsfor”owpassingastationarycylinderat=20and=40=20=40Previous0.73:42.3:2.00:1.89:52.8:1.48:1.48:0.9445.32.222.3554.2Present0.9244.22.232.2153.51.66 100 150 200 2.12 2.13 2.14 2.15 2.16 2.17 2.18 dragversus 50 100 150 200 0x 10 versus(a)(b)Fig.29.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=20inadomainofsize48S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS ¼ð
Þð
ÞFigs.30and31comparethevorticityandthepressuredistributionsoverthecylindersurfacewiththepreviouscomputationalresultstsforRe=20and=40,indicatingagoodagreement.Wealsomonitorthesurfacepositionandthesurfacevelocitydistributionforanobjecttoseewhethertheforcemodelenforcesthedesiredmotionwhilemaintainingtheshape.Fig.32(a)showsthetemporalvariationoftherelativeerrorforthesurfaceposition,whereisde“nedas 90 180 270 360 0 5 10 vorticityersus 0 90 180 270 360 0 0.5 1 1.5 versus (a)(b)Fig.30.(a)Vorticityand(b)pressuredistributionsonthecylindersurfacein”owat=20…lines:currentresults,opencircles:numericalresultsfrom 40 60 80 100 0.025 0.03 0.035 0.04 0.045 0.05 0.055 timeerror 0 90 180 270 360 0.005 0 0.005 0.01 alphavelocity v Fig.32.Flowpastastationarycylinderat=20:(a)relativeshapeerrorversustime;(b)thevelocitydistributiononthecylindersurfaceattime=100. 90 180 270 360 0 5 10 15 vorticityversus 0 90 180 270 360 0 0.5 1 1.5 versus(a)(b)Fig.31.Vorticity(a)andpressure(b)distributionsoncylindersurfacein”owat=40…lines:currentresults,opencircles:numericalresultsfromfrom.30S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS er¼100 Therelativeerrorofthesurfacepositionatsteadystateislessthan0.04%.Fig.32(b)showsthevelocitydistributionaroundthecylindersurfaceattime=100.7.5.2.Re=100and200Figs.33and34arethedragandthelifthistoryfor=100and=200.Figs.35and36showthe”owsaroundthecylinderintermsofvorticityandpressure.TheexpectedtrailingVonKarmanvortexstreetdevelopsinthewake.Table14providesasummaryofresultsfor=100and=200,whereistheStrouhalnumber,i.e.thenondimensionalvortexsheddingfrequency.Ourresultsarewithintherangesofthepreviouslyreported 100 150 200 1.1 1.2 1.3 1.4 1.5 dragversust 100 150 200 0 0.2 0.4 versustFig.33.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=100. 50 60 70 80 90 100 110 120 130 1.2 1.3 1.4 1.5 dragversus 0 20 40 60 80 100 120 0 0.5 1 versusFig.34.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=200.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 0 5 15 20 5 =100 5 15 20 5 =200Fig.35.Vorticitycontoursof”owpassingastationarycylinder. 0 5 10 15 20 0 5 0 5 10 15 20 0 5 Fig.36.Pressurecontoursof”owpassingastationarycylinder. Table14Summaryof”owcharacteristicsfor”owpassingastationarycylinderat=100and=200=100=200Previous1.33±0.014:±0.25:0.164:1.17±0.058:±0.47:0.192:0.192:1.43±0.009±0.3390.1751.45±0.036±0.750.202Present1.423±0.013±0.340.1711.42±0.04±0.660.202S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS values.Wehavealsochangedthedomainsizefor=100from3216to4824withthecorrespondingchangeoffrom640320to960480.ThenewdragandlifthistoryisshowninFig.37,whichshows=1.379±0.009,=±0.31,and=0.167.7.6.Flowarounda”appingwingInthisexample,wesimulatethe”owaroundahoveringwinginsidearigidbox,asdescribedinFig.38=512512and=0.001.Thenon”uidforceisgivenbyEqs.(49)and(50)=2,=2,=0.1and=160.Thewingisanellipsewithchordlengthandaspectratio.Wetakeasthelengthscale, thevelocityscale,and thetimescale,whereistheamplitudeofthewingcentertranslationandisthewing”appingperiod.Thenondimensional”appingperiodis .Thenondimen-sionalizedwingmotionisgovernedbycos 2ta01hðt01sin 2ta0þ/ 170 180 190 200 1.365 1.37 1.375 1.38 1.385 1.39 dragversus 170 180 190 200 0.4 versus(a)(b)Fig.37.(a)Dragand(b)liftcoe
cientsversustimefor”owpassingastationarycylinderat=100inadomainofsize48 1yx(0,0) (6,6) a00.25 Fig.38.Geometryfor”owaroundahovering”apper.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS )isthedisplacementofthewingwithanamplitude)istheangleofattack(weuseinsteadof,sinceisusedastheLagrangianparameterofasurface)withamplitude2,andthephasedierence.TheReynoldsnumberisde“nedas .Werunthesimulationfor=4,=2.5, =0and=157,whicharethesameasusedbyWangang.Fig.39showsfoursnapshotsofthecomputedvorticity“eldsnearthewingduringone”appingperiod.TheyareverysimilartothoseobtainedbyWangWang,wherethephysicalinterpretationwasgiven.WeplotinstantaneousdragandliftinFig.40andcomparethemwiththeresultsofWangWang.Consideringthedierenceinthefar-“eldboundaryconditionsandtheinevitabledierenceofthewingmotionduetotheoscillationsofLagrangianpoints,theagreementisquitegood.Wede“netheshapedistortiontobe Fig.39.Vorticity“eldsaroundahoveringwingof=157atfourdierentinstantsina”appingperiod.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS where()istheprescribedcoordinatesofLagrangianpoint.Theshapedistortionofthewingisverysmall,asseeninFig.41Fig.42(a)and(b)showthevelocityevolutionofLagrangianpointsattheleadingedgeandinthemiddleofthewing,whichindicatesthatitismoredi
culttocontroltheLagrangianpointsattheleadingandtrailingedgestofollowtheprescribedvelocity.7.7.FlowpassingmultiplecylindersWeprovidetwoexamplesbelowtodemonstratethecapabilityande
ciencyofourmethodtosimulate”owswithmultiplemovingobjects. 70 75 80 85 90 95 100 105 110 115 1 versus 70 75 80 85 90 95 100 105 110 115 1 2 3 versusFig.40.(a)Dragand(b)liftcoe
cientsversustimefor”owaroundahovering”apperof=157„solidlines:currentresults,dashedlines:previousresults 125 130 135 140 velocity 125 130 135 140 velocity(a)(b)Fig.42.Velocityonthewingsurface:(a)thevelocityversustimeofaLagrangianpointat=0;(b)thevelocityversustimeofaLagrangianpointat p2. 125 130 135 140 0 1 2 3x 10 Fig.41.Theshapeerrorofahoveringwingat=157.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 7.7.1.TwocylindersmovingwithrespecttoeachotherThisexamplewaspreviouslystudiedbyRussellandWangg.Theinitialgeometry,showninFig.43,isthesameastheirs.Inthissimulation,=640128,and=0.0005.Thefar-“eldbound-ariesarerigidwalls.Thenon”uidforceforbothcylindersisgivenbyEq.=800.Toavoidtheimpulsivestartofthecylinders,weleteachcylinderoscillateaboutitsinitialpositionfortwoperiodsandthenmovetowardtheotherat=40.Themotionofthelowercylinderisgivenby 4psin andthemotionoftheuppercylinderisgivenby 0 5 10 15 20 0 5 0 5 10 15 20 0 5 Fig.44.Flow“eldsaroundtwocylindersmovingwithrespecttoeachotherat=40whentwocylindersareclosest.Contoursof(a)vorticityand(b)pressure. (0,0)xy1WNSE1(16,1.5)u=1 (24,8) Fig.43.Geometryfor”owaroundtwocylindersmovingwithrespecttoeachother.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS 0 5 5 0 5 5 Fig.45.Flow“eldsaroundtwocylindersmovingwithrespecttoeachotherat=40whentwocylindersareseparatedbyadistanceof16.Contoursof(a)vorticityand(b)pressure. 1 1.5 2 2.5 dragversus 20 22 24 26 28 30 1820222426283032 0 0.5 1 versusFig.46.(a)Dragand(b)liftcoe
cientsversustimefortheuppercylinderin”owaroundtwocylindersmovingwithrespecttoeachother=40.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS x16 4psin Fig.44showsthevorticityandpressure“eldsattime=24,whenthetwocylindersareclosesttoeachother,andFig.45showsthevorticityandpressure“eldsattime=32.Fig.46showsthetemporalevolutionofthedragandliftcoe
cientsfortheuppercylinder.WeseethesamedragbehaviorasobtainedbyRussellandWangWang,thatisanincreaseindragasthetwocylindersapproacheachotherandadecreaseindragastheypassincloseproximity.Thetwocylindersarerepulsivewhenapproachingeachotherandbecomeattractiveafterpassedeachother.7.7.2.CylinderstranslatingalongacircleInthislastexample,wesimulatecylinderstranslatingwiththesamespeedalongacircletoexaminethee
ciencyofourmethodforsimulatingmultiplemovingobjects. 0 2 4 0 2 4 0 2 4 0 1 2 3 4 Fig.48.(a)Vorticityand(b)pressure“eldsaround“vecylinderstranslatingaboutacommoncenterwith=20. Table15ComputationalcostintermsofthenumberofcylindersNumber12345678Hours2.582.833.263.453.453.883.894.19Intel(R)Xeon(TM)3.20GHzCPU,512KBcache,1GRAM;CompiledbyIntelFortranCompilerwithcommand:ifort-ftz-fpe0-O3-ip-parallel;=256=0.0002,40,000timeintegrationsteps. (–4,–4)(4,–4)Fig.47.Geometryfor”owaroundcylinderstranslatingaroundacenter.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS ThegeometryofthesimulationisgiveninFig.47.TheReynoldsnumberbasedonthediameterandthetranslationalspeedofacylinderis=20.Thenon”uidforceforeachcylinderismodeledagainbyEqs.and(50)with=20,=20,=0and=800.Simulationsarecarriedoutwith=256256,and=0.0002.Table15liststhecomputationaltimeversusthenumberofcylindersin40,000timeintegrationsteps.Sincethecomputationaltimeassociatedwithacylinderisoforder,thereisnosigni“cantincreaseofthecom-putationaltimewhenoneortwomorecylindersareadded.MostcomputationaltimeisspentontheFFTPoissonsolverforthepressure,whichhascostoforderFig.48showsthevorticityandthepressurecontoursforthe”owcontaining“vecylinders.8.ConclusionsTheimmersedinterfacemethodsharesthesamemathematicalformulationandthereforetheadvantageofPeskinsimmersedboundarymethod.Withtheappropriateinclusionofjumpconditionsin“nitedierenceschemes,theimmersedinterfacemethodasdescribedherecanachievehigherorderaccuracy,maintainasharpinterfaceandpreservevolumes.Speci“cally,ournumericaltestsshowthatsecond-orderspatialaccuracyintermsofthein“nitynormforboththevelocityandthepressurecanbeachieved.Becauseoftheuseofa“xedCartesiangrid,themethodissimpleande
cientfor”owswithmovingobjects.Thecosttotreatanobjectisoforder,whereisthenumberofLagrangianpointsintheobjectrepresentation.Themethodcanbeappliedequallywelltoobjectswithprescribedforceorobjectswithprescribedmotion.Inthecurrentimple-mentationofthemethod,wechoosetoemployuniformgridsandFFTPoissonsolvers,butthemethodisamenabletononuniformgridsandotherPoissonsolvers.Wearenowimplementingthemethodin3Dwiththejumpconditionsthathavebeenpresentedinin.AcknowledgementsWethankProfessorRandallJ.LeVequeandProfessorZhilinLiforhelpfuldiscussionsduringthecourseofthiswork.TheworkissupportedbyAFOSRandPackardFoundation.References[1]R.P.Beyer,R.J.Leveque,Analysisofaone-dimensionalmodelfortheimmersedboundarymethod,SIAMJ.Numer.Anal.29(2)(1992)332…364.[2]M.Braza,P.Chassaing,H.HaMinh,Numericalstudyandphysicalanalysisofthepressureandvelocity“eldsinthenearwakeofacircularcylinder,J.FluidMech.165(1986)79…130.[3]DonnaCalhoun,ACartesiangridmethodforsolvingthetwo-dimensionalstreamfunction-vorticityequationsinirregularregions,J.Comput.Phys.176(2002)231…275.[4]R.Cortez,M.Minion,TheBlobprojectionmethodforimmersedboundaryproblems,J.Comput.Phys.161(2000)428…453.[5]WeinanE,Jian-GuoLiu,Vorticityboundaryconditionandrelatedissuesfor“nitedierenceschemes,J.Comput.Phys.124(1996)[6]E.A.Fadlun,R.Verzicco,P.Orlandi,J.Mohd-Yusof,Combinedimmersed-boundary“nite-dierencemethodsforthree-dimensionalcomplex”owsimulations,J.Comput.Phys.161(2000)35…60.[7]RonaldP.Fedkiw,TariqAslam,BarryMerriman,StanleyOsher,Anon-oscillatoryEulerianapproachtointerfacesinmultimaterial”ows(theGhostFluidMethod),J.Comput.Phys.152(1999)457…492.[8]AaronL.Fogelson,JamesP.Keener,ImmersedinterfacemethodforNeumannandrelatedproblemsintwoandthreedimensions,SIAMJ.Sci.Comput.22(5)(2000)1630…1654.[9]D.Goldstein,R.Handler,L.Sirovich,Modelingano-slip”owboundarywithanexternalforce“eld,J.Comput.Phys.105(1993)[10]YueHao,AndreaProsperetti,Anumericalmethodforthree-dimensionalgas…liquid”owcomputations,J.Comput.Phys.196(2004)[11]FrancisH.Harlow,J.EddieWelch,Numericalcalculationoftime-dependentviscousincompressible”owof”uidwithfreesurface,Phy.Fluids8(12)(1965)2182…2189.[12]HansJohnston,Jian-GuoLiu,Finitedierenceschemesforincompressible”owbasedonlocalpressureboundaryconditions,J.Comput.Phys.180(2002)120…154.[13]HansJohnston,Jian-GuoLiu,Accurate,stableande
cientNavier…Stokessolversbasedonexplicittreatmentofthepressureterm,J.Comput.Phys.199(2004)221…259.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS [14]JungwooKim,DongjooKim,HaecheonChoi,Animmersed-boundary“nite-volumemethodforsimulationsof”owincomplexgeometries,J.Comput.Phys.171(2001)132…150.[15]Ming-ChihLai,CharlesS.Peskin,Animmersedboundarymethodwithformalsecond-orderaccuracyandreducednumericalviscosity,J.Comput.Phys.160(2000)705…719.[16]LongLee,RandallJ.LeVeque,AnimmersedinterfacemethodforincompressibleNavier…Stokesequations,SIAMJ.Sci.Comput.25(3)(2003)832…856.[17]RandallJ.LeVeque,ZhilinLi,Theimmersedinterfacemethodforellipticequationswithdiscontinuouscoe
cientsandsingularsources,SIAMJ.Numer.Anal.31(4)(1994)1019…1044.[18]RandallJ.LeVeque,ZhilinLi,Immersedinterfacemethodsforstokes”owwithelasticboundariesorsurfacetension,SIAMJ.Sic.Comput.18(3)(1997)709…735.[19]ZhilinLi,Ming-ChihLai,TheimmersedinterfacemethodfortheNavier…Stokesequationswithsingularforces,J.Comput.Phys.171(2001)822…842.[20]ZhilinLi,PrivateCommunication.[21]CharlesS.Peskin,Flowpatternsaroundheartvalves:anumericalmethod,J.Comput.Phys.10(1972)252…271.[22]CharlesS.Peskin,Numericalanalysisofblood”owintheheart,J.Comput.Phys.25(1977)220…252.[23]CharlesS.Peskin,BethFellerPrintz,Improvedvolumeconservationinthecomputationof”owswithimmersedelasticboundaries,J.Comput.Phys.105(1993)33…46.[24]CharlesS.Peskin,Theimmersedboundarymethod,ActaNumer.(2002)1…39.[25]A.Prosperetti,H.N.Oguz,PHYSALIS:anewo()methodforthenumericalsimulationofdispersesystems:potential”owofspheres,J.Comput.Phys.167(2001)196…216.[26]A.M.Roma,C.S.Peskin,M.J.Berger,Anadaptiveversionoftheimmersedboundarymethod,J.Comput.Phys.153(1999)509…534.[27]DavidRussell,Z.JaneWang,ACartesiangridmethodformodelingmultiplemovingobjectsin2Dincompressibleviscous”ow,J.Comput.Phys.191(2003)177…205.[28]E.M.Saiki,S.Biringen,Numericalsimulationofacylinderinuniform”ow:applicationofavirtualboundarymethod,J.Comput.Phys.123(1996)450…465.[29]A.L.F.Lima,E.Silva,A.Silveira-Neto,J.J.R.Damasceno,Numericalsimulationoftwo-dimensional”owsoveracircularcylinderusingtheimmersedboundarymethod,J.Comput.Phys.189(2003)351…370.[30]JohnM.Stockie,BrianR.Wetton,Analysisofstinessintheimmersedboundarymethodandimplicationsfortime-steppingschemes,J.Comput.Phys.154(1999)41…64.[31]S.Takagi,H.N.Oguz,Z.Zhang,A.Prosperetti,PHYSALIS:anewmethodforparticlesimulationPartII:two-dimensionalNavier…Stokes”owaroundcylinders,J.Comput.Phys.187(2003)371…390.[32]Yu-HengTseng,JoelH.Ferziger,Aghost-cellimmersedboundarymethodfor”owincomplexgeometry,J.Comput.Phys.192(2003)593…623.[33]C.Tu,C.S.Peskin,Stabilityandinstabilityinthecomputationof”owswithmovingimmersedboundaries,SIAMJ.Statist.Comput.13(1992)70…83.[34]H.S.Udaykumar,R.Mittal,P.Rampunggoon,A.Khanna,AsharpinterfaceCartesiangridmethodforsimulating”owswithcomplexmovingboundaries,J.Comput.Phys.174(2001)345…380.[35]Z.JaneWang,Twodimensionalmechanismforinsecthovering,Phys.Rev.Lett.85(10)(2000).[36]AndreasWiegmann,KennethP.Bube,Theimmersedinterfacemethodfornonlineardierentialequationswithdiscontinuouscoe
cientsandsingularsources,SIAMJ.Numer.Anal.35(1)(1998)177…200.[37]AndreasWiegmann,KennethP.Bube,Theexplicit-jumpimmersedinterfacemethod:“nitedierencemethodsforPDEswithpiecewisesmoothsolutions,SIAMJ.Numer.Anal.37(3)(2000)827…862.[38]ShengXu,Z.JaneWang,Systematicderivationofjumpconditionsfortheimmersedinterfacemethodinthree-dimensional”owsimulation,SIAMJ.Sci.Comput.,2006,inpress.[39]T.Ye,R.Mittal,H.S.Udaykumar,W.Shyy,AnaccurateCartesiangridmethodforviscousincompressible”owswithcompleximmersedboundaries,J.Comput.Phys.156(1999)209…240.[40]LuodingZhu,CharlesS.Peskin,Simulationofa”apping”exible“lamentina”owingsoap“lmbytheimmersedboundarymethod,J.Comput.Phys.179(2002)452…468.S.Xu,Z.J.Wang/JournalofComputationalPhysicsxxx(2006)xxx…xxx ARTICLEINPRESS