INTERNATIONALJOURNALFORNUMERICALMETHODSINENGINEERINGInt.J.Numer.Meth.E - PDF document

INTERNATIONALJOURNALFORNUMERICALMETHODSINENGINEERINGInt.J.Numer.Meth.Engng0000;00:1–26PublishedonlineinWileyInterScience(www.interscience.wiley.com).DOI:10.1002/nmeHybridnitedifference/niteelementversionoftheimmersedboundarymethodBoyceE.Grifth*andXiaoyuLuoLeonH.CharneyDivisionofCardiology,DepartmentofMedicine,NewYorkUniversitySchoolofMedicineandSchoolofMathematicsandStatistics,UniversityofGlasgowSUMMARYTheimmersedboundary(IB)methodisaframeworkformodelingsystemsinwhichanelasticstructureisimmersedinaviscousincompressibleuid.TheIBformulationofsuchproblemsdescribestheelasticityofthestructureinLagrangianformanddescribesthemomentum,viscosity,andincompressibilityoftheuid-structuresysteminEulerianform.InteractionsbetweenLagrangianandEulerianvariablesaremediatedbyintegraltransformswithdeltafunctionkernels.Whendiscretized,theLagrangianequationsareapproximatedonacurvilinearmesh,theEulerianequationsareapproximatedonaCartesiangrid,andaregularizedversionofthedeltafunctionisusedinapproximationstotheLagrangian-Eulerianinteractionequations.Here,weemployaversionoftheIBmethodthatallowsustodiscretizethestructureviastandardLagrangianniteelement(FE)methods.UnlikemostotherextensionsoftheIBmethodthatuseFEstructuraldiscretizations,however,ourapproachretainsanitedifferencediscretizationoftheincompressibleNavier-Stokesequations.AkeyfeatureofournumericalschemeisthatitenablestheuseofLagrangianmesheswithmeshspacingsthatareindependentofthegridspacingofthebackgroundEuleriangrid.Resultsfromcomputationalexperimentsareincludedthatdemonstratetheaccuracyandefciencyofourmethodology.Copyrightc\r0000JohnWiley&Sons,Ltd.Received...KEYWORDS:immersedboundarymethod;uid-structureinteraction;niteelementmethod;nitedifferencemethod;incompressibleelasticity;incompressibleow1.INTRODUCTIONSinceitsintroduction[1,2],theimmersedboundary(IB)methodhasbeenwidelyusedtosimulatebiologicaluiddynamics[3]andotherproblemsinwhichastructureisimmersedinauidow[3–5].Inthispaper,weconsidertheIBmethodforuid-structureinteractionproblemsinvolvinganincompressibleelasticbodythatisimmersedinaviscousincompressibleuid[3].TheIBformulationofsuchproblemsusesaLagrangiandescriptionoftheimmersedstructureandanEuleriandescriptionofthemomentum,viscosity,andincompressibilityoftheuid-structuresystem.LagrangianandEulerianvariablesarecoupledbyintegraltransformswithdeltafunctionkernels.Whenthecontinuousequationsarediscretized,theLagrangianequationsareapproximatedonacurvilinearmesh,theEulerianequationsareapproximatedonaCartesiangrid,andtheLagrangian-Eulerianinteractionequationsareapproximatedbyreplacingthesingulardeltafunctionwitharegularizeddeltafunction.OneadvantageoftheIBformulationisthatitenablestheuseofefcientCartesiangridsolvers,suchasthosebasedonfastFouriertransform(FFT)ormultigrid Correspondenceto:BoyceE.Grifth,LeonH.CharneyDivisionofCardiology,DepartmentofMedicine,NewYorkUniversitySchoolofMedicine,550FirstAvenue,NewYork,NY10016USA,boyce.griffith@nyumc.orgCopyrightc\r0000JohnWiley&Sons,Ltd.Preparedusingnmeauth.cls[Version:2010/05/13v3.00] 2B.E.GRIFFITHANDX.Y.LUOmethods.Anotherstrengthofthisapproachisthatitpermitsnonconformingdiscretizationsoftheuidandstructure.Specically,theIBmethoddoesnotrequiredynamicallygeneratedbody-ttedmeshes,apropertythatisespeciallyusefulforproblemsinvolvinglargestructuraldeformationsorcontactbetweenstructures.InmanyapplicationsoftheIBmethod,theelasticityoftheimmersedstructureisdescribedbysystemsofbersthatresistextension,compression,orbending.Suchdescriptionscanbewellsuitedforthehighlyanisotropicmaterialsencounteredinbiologicalapplications.Fibermodelsarealsoconvenienttouseinpracticebecausetheypermitanespeciallysimpledenitionanddiscretization.Theber-basedapproachtoelasticitymodelingalsopresentscertainchallenges.Forinstance,itcanbedifculttoincorporaterealisticshearpropertiesorexperimentallybasedconstitutivelawsintobermodels.Nonetheless,theber-basedmodelingapproachtypicallyemployedwiththeIBmethodhasfacilitatedsignicantworkinbiouiddynamics[6–11],includingthree-dimensionalsimulationsofcardiacuiddynamics[12–21].Overthepastdecade,severaldistinctresearcheffortshavesoughttousemoregeneralmaterialmodelswiththeIBmethod.Forinstance,recentworkhasledtothedevelopmentofanenergyfunctional-basedversionoftheconventionalIBmethodthatobtainsanodalapproximationtotheelasticforcesgeneratedbyanimmersedhyperelasticmaterialviaaniteelement(FE)typeapproximationtothedeformationofthematerial[22,23].Otherworkhasalsoledtothedevelopmentoftheimmersedstructuralpotentialmethod,whichusesameshlessmethodtodescribethemechanicsofanincompressiblehyperelasticstructureimmersedinuid[24].Theformulationusedbyeachofthesenumericalschemesreliesontheavailabilityofanelasticenergyfunctional,i.e.,thesemethodsarespecializedtohyperelasticmaterialmodels.Aseparatelineofresearchledtothedevelopmentoftheimmersedniteelement(IFE)method[25–27],whichcanbeviewedasageneralizationoftheIBmethodinwhichFEmethodsareusedforboththeuidandthestructure.LiketheIBmethod,theIFEmethodcouplesLagrangianandEulerianvariablesbydiscretizedintegraltransformswithregularizeddeltafunctionkernels,althoughdifferentfamiliesofsmoothedkernelfunctionsaregenerallyusedbythetwomethods.OtherworkhasledtothedevelopmentofafullyvariationalIBmethodthatavoidsregularizeddeltafunctionsaltogether[28,29].Inthispaper,weintroduceanalternativeapproachtousingFEmechanicsmodelswiththeIBmethodthatcombinesaCartesiangridnitedifferencemethodforincompressibleuiddynamicswithanodalFEmethodfornonlinearelasticity.Todevelopournumericalscheme,weusethecontinuousIBformulationintroducedbyBofetal.[28].ThemathematicalformulationandnumericalmethodrequireonlythespecicationofaLagrangianstresstensortodescribethematerialresponseoftheimmersedstructure.Ournumericalmethodalsoreadilyaccommodatesarbitrarychoicesofnodalniteelements,anditcanbeeasilyextendedtotreatnoninterpolatoryFEbases(e.g.,NURBS).WeconsidertwoweakformulationsoftheequationsofmotionsuitableforusewithstandardC0FEmethodsforstructuralmechanics.Oneoftheseformulations,referredtohereinastheuniedweakform,issimilartothoseusedbyearlierIB-likemethods[25–29].Theother,whichdoesnotappeartohavebeenusedpreviouslytoconstructanumericalscheme,partitionstheelasticforcedensityintoaninternalforcethatissupportedthroughouttheimmersedelasticstructureandatransmissionforcethatissupportedonlyonthesurfaceoftheimmersedstructure.Weprovidenumericalexamplesthatdemonstratethatimprovementsinaccuracycanbeobtainedbyemployingthepartitionedformulationinsteadoftheuniedformulation,especiallyforcasesinwhichtheLagrangianmeshisrelativelycoarsecomparedtotheEuleriangrid.Conventionally,regularizeddeltafunctionsareusedbytheIBmethodbothtospreadtheforcesgeneratedbythestructuredirectlyfromthenodesoftheLagrangianmeshtotheEuleriangridandalsotointerpolatevelocitiesfromtheEuleriangriddirectlytothenodesoftheLagrangianmesh.ThisdiscreteLagrangian-EuleriancouplingapproachisalsoadoptedbyIB-likemethodssuchastheimmersedstructuralpotentialmethod[24]andtheimmersedniteelementmethod[25–27].AsignicantlimitationofthisapproachisthatifthephysicalspacingoftheLagrangiannodesistoolargeincomparisontothebackgroundEuleriangrid,severe“leaks”willdevelopatuid-structureinterfaces.AnempiricalrulethatgenerallypreventssuchleaksistorequiretheLagrangianmeshCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD3 U(U;t):(U;t)7!\n\nFigure1.LagrangianandEuleriancoordinatesystems.TheLagrangianmaterialcoordinatedomainisU,andtheEulerianphysicalcoordinatedomainis\n.Thephysicalpositionofmaterialpointsattimetis(s;t),thephysicalregionoccupiedbythestructureis(U;t),andthephysicalregionoccupiedbytheuidis\nn(U;t).tobeapproximatelytwiceasneastheEuleriangrid[3].BecausehighEulerianresolutionisrequiredtocapturethethinboundarylayerscharacteristicofmoderate-to-highReynoldsnumberows,however,followingthisempiricalrulecanrequireusingverydenseLagrangianmeshesthatmaygeneratesignicantnumericalstiffnessinthediscretizedequations.Moreover,forproblemsinvolvinglargedeformations,astructuralmeshthatisinitially“watertight”maydeforminamannerthatultimatelyyieldsleaksasthesimulationprogresses.AkeycontributionofthispaperisthatitintroducesanewapproachtodiscretizingtheequationsofLagrangian-EulerianinteractionthatovercomesthislongstandinglimitationoftheIBmethod.Specically,ratherthanspreadingforcesfromandinterpolatingvelocitiestothenodesoftheLagrangianmesh,weinsteadspreadforcesfromandinterpolatevelocitiestoquadraturepointswithintheinteriorsoftheLagrangianniteelements.ThisapproachtakesadvantageoftheadditionalgeometricalinformationprovidedbytheFEdescriptionofthestructure,whichyieldsapproximationsnotonlytothecurrentpositionsofthenodesoftheLagrangianmesh,butalsotothepositionsofanymaterialpointswithinthestructure.Suchgeometricalinformationisnotreadilyavailableinconventionalber-basedIBmethods,andthisinformationdoesnotappeartohavebeenfullyexploitedinmostearlierextensionsoftheIBmethodthatincorporateFE-typestructuraldiscretizations.Animportantfeatureofourdiscretizationapproachisthatobtainingawatertightstructuresimplyrequiresaquadratureschemewithsufcientlymanyquadraturepointstopreventleaks.Moreover,becauseitisstraightforwardtolocallyadaptthequadratureschemetoaccountforstructuraldeformations,thisapproachcanensurethatthestructureremainswatertighteveninthepresenceofverylargestructuraldeformations.WepresentnumericalteststhatdemonstratethatourschemecanyieldaccurateresultsevenforLagrangianmeshesthataresignicantlycoarserthanthebackgroundEuleriangrid,andwedemonstratethepotentialgainsinthecomputationalefciencyofferedbythisapproach.2.CONTINUOUSFORMULATIONSIntheIBformulationofproblemsofuid-structureinteraction,themomentum,velocity,andincompressibilityofthecoupleduid-structuresystemaredescribedinEulerianform,whereastheelasticityoftheimmersedstructureisdescribedinLagrangianform.Letx=(x1;x2;:::)2\nRdd=2or3,denoteCartesianphysicalcoordinates,with\ndenotingthephysicalregionthatisoccupiedbythecoupleduid-structuresystem;lets=(s1;s2;:::)2URddenoteLagrangianmaterialcoordinatesthatareattachedtothestructure,withUdenotingtheLagrangiancoordinatedomain;andlet(s;t)2\ndenotethephysicalpositionofmaterialpointsattimet.Thephysicalregionoccupiedbythestructureattimetis(U;t)\n,andthephysicalregionoccupiedbytheuidattimetis\nn(U;t).SeeFig.1.WedonotassumethattheLagrangiancoordinatesaretheinitialcoordinatesoftheelasticstructure(i.e.,wedonotassumethat(s;0)=s),nor,moregenerally,dowerequirethatU\nTouseanEuleriandescriptionoftheuidandaLagrangiandescriptionoftheelasticityoftheimmersedstructure,itisnecessarytodescribethestressoftheuid-structuresysteminbothCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 4B.E.GRIFFITHANDX.Y.LUOEulerianandLagrangianforms.Let=(x;t)denotetheCauchystresstensorofthecoupleduid-structuresystem.Then(x;t)=(f(x;t)+e(x;t)forx2(U;t)f(x;t)otherwise,(1)inwhichf(x;t)isthestresstensorofaviscousincompressibleuid,ande(x;t)isthestresstensorthatdescribestheelasticityoftheimmersedstructure.Theuidstresstensoristheusualoneforaviscousincompressibleuid,f=pI+hru+(ru)Ti;(2)inwhichp=p(x;t)isthepressure,isthedynamicviscosityoftheuid,andu=u(x;t)istheEulerianvelocityeld.TodescribetheelasticityofthestructurewithrespecttotheLagrangianmaterialcoordinatesystem,itisconvenienttousetherstPiola-KirchhoffelasticstresstensorPe(s;t),whichisdenedsothatZ@VPe(s;t)NdA(s)=Z@(V;t)e(x;t)nda(x)(3)foranysmoothregionVU,inwhichN=N(s)istheoutwardunitnormalalong@V,andn=n(x;t)istheoutwardunitnormalalong@(V;t).Forsimplicity,weshallrestrictourattentiontohyperelasticconstitutivemodels,whichmaybecharacterizedbyastrain-energyfunctionale=e(F),inwhichF=F(s;t)=rs(s;t)=@ @s(s;t)isthedeformationgradientassociatedwiththemapping:(U;t)\n.Forsuchconstitutivelaws,Pe(s;t)=@We @F(s;t).Weremark,however,thatthepresentformulationcanaccommodateincompressiblematerialsthataredescribedonlyintermsofaLagrangianstresstensor,andthatthisformulationisthereforenotrestrictedtohyperelasticmaterialmodels.2.1.StrongformulationThestrongformoftheequationsofmotionfortheuid-structuresystemis:@u @t(x;t)+u(x;t)
ru(x;t)=r
(x;t)=rp(x;t)+r2u(x;t)+fe(x;t);(4)r
u(x;t)=0;(5)fe(x;t)=ZUrs
Pe(s;t)(x(s;t))dsZ@UPe(s;t)N(s)(x(s;t))dA(s);(6)@ @t(s;t)=Z\nu(x;t)(x(s;t))dx;(7)inwhichisthemassdensityofthecoupleduid-structuresystem,fe(x;t)istheEulerianelasticforcedensity,and(x)=Qd=1(x)isthed-dimensionaldeltafunction.Weomitthederivationofthisformulationoftheequationsofmotion,whichwasintroducedbyBofetal.[28].Intheequationsofmotion,Eqs.(4)and(5)arethestandardEulerianincompressibleNavier-Stokesequations,exceptthattheright-handsideofthemomentumequation,Eq.(4),isaugmentedbyanEulerianelasticforcedensity,fe(x;t)=r
e(x;t),thatisdeterminedfromtheLagrangiancongurationoftheimmersedstructureviaEq.(6).TheconversionfromtheLagrangiantotheEuleriandescriptionoftheforcesgeneratedbytheelasticityoftheimmersedstructureismediatedbytwointegraltransformswithdeltafunctionkernels.TheseintegraltransformsappearinEq.(6)Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD5andconvertrs
Pe,theLagrangianinternalforcedensity,andPeN,theLagrangiantransmissionforcedensity,intoequivalentEuleriandensities.NoticethatthesetwoLagrangianforcedensitieshavetotallydifferentcharacters:TheinternalforcedensityissupportedthroughoutUandhasunitsofforceperunitvolume,whereasthetransmissionforcedensityissupportedonlyalong@Uandhasunitsofforceperunitarea.AssumingthatPe(s;t)issufcientlysmooth,theinternalforceactsasaregular(i.e.,nonsingular)bodyforceontheuidandmaybetreatedwithhigher-orderaccuracybytheIBmethod[28,30].Thetransmissionforceactsasasingularforcelayerontheuid,andalthoughthisforcewillinducejumpsinthepressureandshearstressalong@(U;t)suchforcelayersmayalsobereadilytreatedbytheIBmethod.InmanypresentationsoftheIBmethodforuid-structureinteraction(see,e.g.,Peskin[3]),thetransmissionforceisnotincludedintheequationsofmotion,althoughitseffectmustbeeitherimplicitlyorexplicitlyincludedinthenumericalscheme.TheexplicitinclusionofthetransmissionforceisthekeydifferencebetweenthisformulationandtheconventionalIBformulationandfacilitatesthedevelopmentofweakformsoftheequations.AnintegraltransformisalsousedinEq.(7)todeterminethevelocityoftheimmersedelasticstructurefromthematerialvelocityeldu(x;t).Thedeningpropertyof(x)impliesthatEq.(7)isequivalentto@ @t(s;t)=u((s;t);t),whichmaybeinterpretedastheno-slipandno-penetrationconditionsofaviscousincompressibleuid.Notice,however,thattheno-slipandno-penetrationconditionsdonotappearasconstraintsontheuidmotion.Instead,theseconditionsdeterminethemotionoftheimmersedstructure.2.2.WeakformulationsToobtainversionsoftheequationsofmotionthatallowustousestandardC0FEmethodsfornonlinearelasticity,weconsidertwodifferentformulationsthateachemployaweakformoftheLagrangianequations.Thesetwoformulationsareequivalentinthecontinuoussetting;however,whendiscretized,theyyieldnumericalschemesthataregenerallydifferent.BecauseweuseanitedifferenceschemetoapproximatetheincompressibleNavier-Stokesequations,wedonotemployaweakformulationfortheEulerianequations.Werefertoourrstweakformoftheproblemastheuniedweakformulationbecauseitincludesonlyasingle,uniedbodyforcingtermthataccountsforboththeregularinternalelasticforcedensityandthesingulartransmissionelasticforcedensityofthestrongformoftheequations.Thisformulationis:@u @t+u
ru=rp+r2u+f;(8)r
u=0;(9)f(x;t)=ZUF(s;t)(x(s;t))ds;(10)ZUF(s;t)
(s)ds=ZUPe(s;t):rs(s)ds;8(s);(11)@ @t(s;t)=Z\nu(x;t)(x(s;t))dx;(12)inwhichf(x;t)andF(s;t)aretheEulerianandLagrangiantotalelasticforcedensities,and(s)isanarbitraryLagrangiantestfunctionthatisnotassumedtovanishon@U.ThisweakformoftheequationsofmotionissimilartotheformulationemployedintheIFEmethod[25–27]andthefullyvariationalIBmethod[28,29].NoticethatF(s;t)accountsforboththeinternalandtransmissionforcedensitiesofthestrongformoftheequationsinthesensethatZUF(s;t)
(s)ds=ZUPe(s;t):rs(s)ds(13)=ZU(r
Pe(s;t))
(s)dsZ@U(Pe(s;t)N(s))
(s)dA(s);(14)Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 6B.E.GRIFFITHANDX.Y.LUOforany(s)Oursecondweakformoftheproblem,whichwerefertoasthepartitionedweakformulationissimilartotheuniedformulation,exceptthatittreatstheinternalandtransmissionelasticforcedensitiesseparately,asfollows:@u @t+u
ru=rp+r2u+g+t;(15)r
u=0;(16)g(x;t)=ZUG(s;t)(x(s;t))ds;(17)ZUG(s;t)
(s)ds=ZUPe(s;t):rs(s)ds+Z@UPe(s;t)N(s)
(s)dA(s);8(s);(18)t(x;t)=Z@UT(s;t)(x(s;t))dA(s);(19)T=PeN;(20)@ @t(s;t)=Z\nu(x;t)(x(s;t))dx;(21)inwhichg(x;t)andG(s;t)aretheEulerianandLagrangianinternalelasticforcedensities,t(x;t)andT(s;t)aretheEulerianandLagrangiantransmissionelasticforcedensities,and(s)isanarbitraryLagrangiantestfunction.Noticethathere,T(s;t)isequaltothetransmissionforcedensityofthestrongformulation,andthatG(s;t)isweaklyequivalenttotheinternalforcedensityofthestrongformoftheequations,ascanbeshownbyintegratingtheright-hand-sideofEq.(18)byparts.3.SPATIALDISCRETIZATIONWeconsideronlythetwo-dimensionalcasefortheremainderofthepaper.Theextensionofthenumericalschemetothecased=3isstraightforward,andimplementationsford=2and3areprovidedbytheopen-sourceIBAMRsoftware[31].3.1.EuleriandiscretizationTodiscretizetheincompressibleNavier-Stokesequationsinspace,weemployastaggered-gridnitedifferencescheme,anapproachthatyieldssuperioraccuracywhenusedwiththeIBmethodascomparedtocollocateddiscretizations(i.e.,purelycell-ornode-centeredschemes)[32].Tosimplifytheexposition,weassumethat\nistheunitsquareandisdiscretizedonaregularNNCartesiangridwithgridspacings
x1=
x2=h=1 N.Let(i;j)labeltheindividualCartesiangridcellsforintegervaluesofiandj0i;jN.ThecomponentsoftheEulerianvelocityeldu=(u1;u2)areapproximatedatthecentersofthex1-andx2-edgesoftheCartesiangridcells,i.e.,atpositionsx1 2;j=ih;j+1 2
h
andxi;j1 2=i+1 2
h;jh
,respectively.AstaggeredschemeisalsousedfortheEulerianbodyforcef=(f1;f2).Weusethenotation(u1)1 2;j(u2)i;j1 2(f1)1 2;jand(f2)i;j1 2todenotethediscretevaluesofuandf.ThepressurepisapproximatedatthecentersoftheCartesiangridcells,anditsdiscretevaluesaredenotedpi;jLetrh
ur
urhprp,andr2hur2urespectivelydenotestandardsecond-orderaccuratenitedifferenceapproximationstothedivergence,gradient,andLaplaceoperators[33].Inthisapproach,rh
uisdenedatcellcenters,whereasbothrhpandr2huaredenedatcelledges.Weuseastaggered-gridversion[32,33]ofthexsPPM7variant[34]ofthepiecewiseparabolicmethod(PPM)[35]todiscretizethenonlinearadvectionterms.Whereneeded,physicalboundaryconditionsarediscretizedandimposedalong@\nasdescribedbyGrifth[33].Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD7Ifuandvarediscretestaggered-gridvectorelds,wedenoteby[u]and[v]thecorrespondingvectorsofgridvalues.If\nhasperiodicboundaries,wedenethediscreteL2innerproducton\nby(u;v)x=[u]T[v]h2:(22)Minoradjustmentstothisdenitionarerequiredwhennonperiodicphysicalboundaryconditionsareused[33].3.2.LagrangiandiscretizationLetTh=[UbeatriangulationofUcomposedofelementsU.WedenotebyfslgMl=1thenodesofthemesh,andbyfl(s)gMl=1interpolatoryLagrangianbasisfunctions.Ourformulationisindependentofthechoiceofthenodalbasisfunctionsandmaybeusedwithanynodalniteelements.Inparticular,althoughournumericalexamplesconsideronlyQ1(bilinear)elements,itisstraightforwardtoemployhigher-orderelements.Ageneralizationtononinterpolatorybasisfunctions(e.g.,NURBS)wouldbestraightforwardbuthasnotyetbeenattempted.Wedenotethetime-dependentphysicalpositionsofthenodesoftheLagrangianmeshbyfl(t)gMl=1.UsingtheLagrangianbasisfunctions,wedeneanapproximationto(s;t)byh(s;t)=MXl=1l(t)l(s):(23)BecauseweuseinterpolatoryFEbasisfunctions,h(sl;t)=l(t).AnapproximationtothedeformationgradientisgivenbyFh(s;t)=@ @sh(s;t)=MXl=1l(t)@ @sl(s):(24)UsingFh(s;t),wecomputedirectlyPeh(s;t)andTh(s;t),approximationstotherstPiola-KirchhoffstresstensorandtheLagrangiantransmissionforcedensity,respectively.ForC0Lagrangianbasisfunctions,h(s;t)isacontinuousfunctionofs,butFh(s;t)isgenerallydiscontinuousatinternalelementboundaries.Hence,PehandTharealsogenerallyonlypiecewisecontinuous.WeapproximatetheLagrangianforcedensitiesF(s;t)andG(s;t)byFh(s;t)=MXl=1Fl(t)l(s);and(25)Gh(s;t)=MXl=1Gl(t)l(s):(26)ThenodalvaluesfFlgMl=1andfGlgMl=1mustbedeterminedfromPeh(s;t).ByrestrictingthetestfunctionstobelinearcombinationsoftheLagrangianbasisfunctions,wehavethatEq.(11)becomes,afterrearrangingterms,MXl=1ZUl(s)m(s)dsFl(t)=ZUPeh(s;t)rsm(s)ds;(27)foreachm=1;:::;M.Similarly,Eq.(18)becomesMXl=1ZUl(s)m(s)dsGl(t)=ZUPeh(s;t)rsm(s)ds+Z@UPeh(s;t)N(s)m(s)dA(s);(28)Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 8B.E.GRIFFITHANDX.Y.LUOforeachm=1;:::;MLetting[F]denotethevectorofnodalcoefcientsofFh,wewriteEq.(27)as[M][F]=[B];(29)inwhich[M]isthemassmatrixthathasentriesoftheformRUl(s)m(s)ds.Eq.(28)mayberewrittensimilarly.Themassmatrix[M]canalsobeusedtoevaluatetheL2innerproductofLagrangianfunctionsonU.Inparticular,foranyUh(s;t)=PlUl(t)l(s)andh(s;t)=Pll(t)l(s)(Uh;h)s=[U]T[M][]:(30)Differentchoicesofmassmatrices(e.g.,lumpedmassmatrices)inducedifferentdiscreteinnerproductsonUTosimplifynotation,intheremainderofthispaperwedropthesubscript“h”fromournumericalapproximationstotheLagrangianvariables.3.3.Lagrangian-EulerianinteractionAsintheconventionalIBmethod,weapproximatethesingulardeltafunctionkernelappearingintheLagrangian-Eulerianinteractionequationsbyasmoothedd-dimensionalDiracdeltafunctionh(x)thatisofthetensor-productformh(x)=Qd=1h(x).Inthiswork,wetaketheone-dimensionalsmootheddeltafunctionh(x)tobethefour-pointdeltafunctionofPeskin[3].Tocomputeanapproximationtof=(f1;f2)ontheCartesiangrid,weconstructforeachelementU2ThaGaussianquadraturerulewithNquadraturepointssQ2Uandweights!QQ=1;:::;N.Wethencomputef1andf2ontheedgesoftheCartesiangridcellsvia(f1)1 2;j=XUe2ThNeXQ=1F1(sQ;t)h(x1 2;j(sQ;t))!Q;and(31)(f2)i;j1 2=XUe2ThNeXQ=1F2(sQ;t)h(xi;j1 2(sQ;t))!Q;(32)withF(s;t)=(F1(s;t);F2(s;t)).Weusetheshorthandf=SF;(33)inwhichS=S()istheforce-prolongationoperatorimplicitlydenedbyEqs.(31)and(32).Acorrespondingvelocity-restrictionoperatorR=R()isusedtodeterminethemotionofthenodesoftheLagrangianmeshfromtheCartesiangridvelocityeldviad dt=Ru:(34)TherearemanypossiblewaystoconstructR;however,wehavefoundthataneffectiveapproachistorequired d=Rutosatisfythediscretepoweridentity,F;d dts=(f;u)x;(35)whichimpliesthatthesemi-discreteuniedformulationconservesenergyduringLagrangian-Eulerianinteraction[3].Thispoweridentitycanberewrittenas(F;Ru)s=(SF;u)x;(36)i.e.,R=SCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD9A.B.C. Figure2.ProlongingtheelasticforcedensityfromtheLagrangianmeshontotheEuleriangrid.StartingwithanapproximationtotheforcedensityatthenodesoftheLagrangianmesh(A),weusetheinterpolatoryFEbasisfunctionstoapproximatetheforcedensityatquadraturepointsintheelementinterior(B),andthenspreadtheinterpolatedforcesfromthequadraturepointstothebackgroundEuleriangridusingthesmootheddeltafunctionh(x)(C).ThisapproachpermitsLagrangianmeshesthataresignicantlycoarserthantheEuleriangridsolongasthe“net”or“submesh”ofquadraturepointsissufcientlydense.Densernetsofquadraturepointscanbeobtainedbyincreasingtheorderofthenumericalquadraturescheme.ToconstructRexplicitly,itisconvenienttousematrixnotation.Identifying[S]and[R]withthematrixrepresentationsofSandR,wehavethat[f]=[S][F];and(37)d[] dt=[R][u]:(38)Eq.(36)thenbecomes[F]T[M][R][u]=([S][F])T[u]h2:(39)IfEq.(39)istoholdforany[F]and[u],then[R]mustbedenedvia[R]=[M]1[S]Th2:(40)Inourtime-steppingscheme,whichisstatedbelowinSec.4,noticethatweneedonlytoapply[R]todiscretevelocityeldsdenedontheCartesiangrid.Specically,wedonotneedtocompute[R]explicitly.ItisstraightforwardtoshowthatthisconstructionofRimpliesthatd d(s;t)isanapproximationtotheL2projectionoftheLagrangianvectoreldUIB(s;t)=(UIB1(s;t);UIB2(s;t)),withUIB1(s;t)=Xi;j(u1)1 2;jh(x1 2;j(s;t))h2;and(41)UIB2(s;t)=Xi;j(u2)i;j1 2h(xi;j1 2(s;t))h2:(42)BecausethecomponentsofUIB(s;t)arenotgenerallylinearcombinationsoftheLagrangianbasisfunctions,andthereforegenerallyd d=UIBForthesemi-discretizedpartitionedweakformulation,giscomputedontheCartesiangridviag=SG.TheEuleriantransmissionforcedensitytiscomputedinasimilarmanner,butinthiscase,thenumericalintegrationoccursonlyoverthoseelementboundariesthatcoincidewith@U.Weusetheshorthandt=S@UTtodenotethisoperation.Weusethesameregularizeddeltafunctionh(x)toconstructbothSandS@U.Forsimplicity,weusethesamevelocity-restrictionoperatorforbothformulations.ThischoiceensuresthatthetwoformulationscoincidewheneverT0TheLagrangian-EulerianinteractionoperatorsintroducedinthisworkaredifferentfromanalogousoperatorsgenerallyusedinstandardIBmethods.StandardIBmethodsandschemessuchastheIFEmethoduseregularizeddeltafunctionstoapplynodalforcesdirectlytotheCartesiangridandtointerpolateCartesiangridvelocitiesdirectlytotheLagrangiannodes[3].Insuchschemes,Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 10B.E.GRIFFITHANDX.Y.LUOf(x;t)wouldbeapproximatedontheEuleriangridbyexpressionssimilarto(fIB1)1 2;j=MXl=1(F1)l(t)h(x1 2;jl(t))!IBl;and(43)(fIB2)i;j1 2=MXl=1(F2)l(t)h(xi;j1 2l(t))!IBl;(44)inwhichfIBdenotestheEulerianforcedeterminedbythestandardIBforce-spreadingoperatorand!IBldenotesthevolumeassociatedwithLagrangiannodel.Inthisapproach,eachnodalforceFlisappliedonlytoCartesiangridcellswithinthesupportofh(xl),andtheLagrangianmeshmustthereforebenerthantheCartesiangridtoavoidleaks[3].Thecorrespondingapproachtovelocityrestrictionusedbysuchmethodswouldbetosetdl d=UIB(sl;t)Ourforce-prolongationoperatorcanbeseenasthecompositionoftwooperations:rst,thenodalvaluesofFareinterpolatedfromthemeshnodestoa“net”or“submesh”quadraturepointswithineachelement;then,h(x)isusedtospreadthoseinterpolatedforcedensitiesfromthequadraturepointstotheCartesiangrid.SeeFig.2.Velocityrestrictionissimilar.First,theCartesianvelocityeldisinterpolatedtothequadraturepointsusingh(x);thesesamplesofUIBarethenusedtocomputetheright-hand-sideofanL2-projectionequation,andthesolutionofthisprojectionequationdeterminesd d.Inourapproach,theLagrangianstructureiswatertightsolongasthenetofquadraturepointsissufcientlydense.Densernetsofquadraturepointscanbeobtainedbyincreasingtheorderofthequadraturerule,andthismaybedoneinanadaptivemannerasthesimulationprogresses.Inournumericaltests,weuseadaptiveGaussianquadraturerulestoconstructSandRthatprovide,onaverage,atleasta33netofquadraturepointsperCartesiangridcell.4.TEMPORALDISCRETIZATION4.1.Basictime-steppingschemeLetnun,andpn1 2denotetheapproximationstothevaluesofanduattimetn,andtothevalueofpattimetn1 2,respectively.Weadvancethesolutionvaluesforwardintimebythetimeinterval
tasfollows.First,wedetermineapreliminaryapproximationtothedeformedstructurecongurationattimetn+1via~n+1n
t=Rnun;(45)withRn=R(n),andwedeneanapproximationtoattimetn+1 2byn+1 2=~n+1+n 2:(46)Then,wesolveun+1un
t+n+1 2=rhpn+1 2+r2hun+1+un 2+fn+1 2;(47)rh
un+1=0;(48)fn+1 2=Sn+1 2Fn+1 2;(49)n+1n
t=Rn+1 2un+1+un 2;(50)forn+1un+1,andpn+1 2,inwhichn+1 2=3 2un
rhun1 2un1
rhun1iscomputedviaaPPM-typeapproximationtothenonlinearadvectionterm[32,33].Thetime-steppingschemeforthepartitionedweakformulationisanalogous.NoticethatsolvingEqs.(47)–(50)forn+1Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD11un+1,andpn+1 2requiresthesolutionofaCrank-Nicolson-typediscretizationofthetime-dependentincompressibleStokesequations.WesolvethissystemofequationsviatheexibleGMRES(FGMRES)algorithm[36],usingunandpn1 2asinitialapproximationstoun+1andpn+1 2,andusingapressure-freeprojectionmethodwithmultigridsubdomainsolversasapreconditioner[33].4.2.InitialtimestepBecausetimestep-laggedvaluesofuandpareusedbythetime-steppingschemeofSec.4.1,wecannotusethatschemefortheinitialtimestep.Instead,weuseatwo-steppredictor-correctormethod:First,wesolve~un+1un
t+n=rh~pn+1 2+r2h~un+1+un 2+fn;(51)rh
~un+1=0;(52)fn=S(n)F(n);(53)~n+1n
t=R(n)un;(54)for~n+1~un+1,and~pn+1 2,withn=un
rhun.Becausewedonothaveaninitialvalueforthepressure,weusep=0asaninitialguessforpn+1 2.Second,wesetn+1 2=~n+1+n 2andsolveEqs.(47)–(50)forn+1un+1,andpn+1 2,exceptthatweusen+1 2=un+1 2
rhun+1 2withun+1 2=1 2~un+1+un
5.IMPLEMENTATIONThisversionoftheIBmethodisimplementedintheopen-sourceIBAMRsoftware[31],whichisaC++frameworkfordevelopinguid-structureinteractionmodelsthatusetheIBmethod.IBAMRprovidessupportfordistributed-memoryparallelismandadaptivemeshrenement(AMR).IBAMRreliesupontheSAMRAI[37–39],PETSc[40–42],hypre[43,44],andlibMesh[45,46]librariesformuchofitsfunctionality.6.NUMERICALRESULTS6.1.ThickellipticalshellThissetofnumericaltestsusesathickellipticalshelllikethatemployedinthecontextoftheconventionalIBmethod[18,30]andthefullyvariationalIBmethod[28]todemonstratethattheIBmethodcanobtainhigher-orderconvergenceratesforcertainproblems.Inthesecomputations,thephysicaldomainis\n=[0;1][0;1]withperiodicboundaryconditions,and,followingBofetal.[28],theLagrangiancoordinatedomainisU=[0;2R][0;w],withR=0:25andw=0:0625andwithperiodicboundaryconditionsinthes1direction.Thecoordinatemapping:(U;t)\nisgivenattimet=0by(s;0)=(cos(s1=R)(R+s2)+0:5;sin(s1=R)(R+\r+s2)+0:5):(55)For\r=0,theinitialcongurationoftheshellisacircularannuluswithinnerradiusRandthicknessw.Weuse\r=0forstaticproblemsand\r=0:1fordynamicproblems.Ineithercase,wediscretize\nusinganN-by-NCartesiangrid,andwediscretizeUusinga28M-by-Mmeshofbilinear(Q1)elements,withM=1 MfacN 16.TheLagrangiandiscretizationisconstructedsothatthenodesoftheLagrangianmesharephysicallyseparatedbyadistanceofapproximatelyMfac
x.RepresentativenumericalresultsareshowninFig.3.Althoughthisisarelativelysimplebenchmarkproblem,thestaticversion(\r=0)isoneoftheonlytestproblemsavailablefortheIBmethodthatweknowofthatpermitsasimpleanalyticCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 12B.E.GRIFFITHANDX.Y.LUOA. B. C. Figure3.Representativeresultsfromthedynamic(\r=0:1)versionoftheorthotropicshellproblemofSec.6.1.2forN=128andthepartitioned(split)weakformulationoverthetimeinterval0t1:25.ThecomputedpressureandstructuredeformationforMfac=4areshowninAandB.ThecomputeddeformationsobtainedwithMfac=1andMfac=4arecomparedinC.NoticethehighlevelofagreementbetweenthecomputeddeformationsinC.solution.Moreover,becausecertainchoicesofPeyieldproblemsforwhichtheIBmethodisabletoattainhigher-orderconvergencerates,thisproblemallowsustoverifythatourimplementationattainsitsdesignedorderofaccuracy.6.1.1.AnisotropicshellWerstconsideranidealizedanisotropicshellthatisdenedintermsofthestrain-energyfunctional[28],e(F)=e 2w\r\r\r\r@ @s1(s;t)\r\r\r\r2=e 2wF1F1;(56)withrepeatedindicesimplyingsummation.Inthiscase,Pe(s;t)=@We @F(s;t)=e w @1 @s10@2 @s10!=e wF110F210:(57)Thisenergyfunctionalandstresstensorcorrespondtoanidealizedanisotropicelasticmaterialcomposedofacontinuousfamilyofextension-resistantbersthatwrapthethickshell.BecauseUisperiodicinthes1direction,PeN0along@U.Thisreectsthefactthatnoneofthebersterminatealongtheboundaryofthestructure.Becausethetransmissionforcevanishesinthiscase,theuniedandpartitionedweakformulationsareidentical.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD13 L1erroruL2errorL1errorNpN 641282565126412825651264128256512641282565126412825651264128256512102101100103102101104103102106105104106105104106105104Figure4.ErrorsinuandpindiscreteL1,L2,andL1normsforthestatic(\r=0)versionoftheidealizedanisotropicshellproblemofSec.6.1.1.ErrorsforMfac=1areplottedinblue,anderrorsforMfac=4areplottedinred.Referencelineswithslope-2areprovidedfortheuerrordata.Forp,referencelineswithslopes-2,-1.5,and-1areprovidedfortheL1,L2,andL1normdata,respectively.Whenweset\r=0,theinitialcongurationoftheelasticshellisanequilibriumcongurationofthecontinuousproblem,sothatu(x;t)0.RequiringR\np(x;t)dx=0,itcanbeshown[28]thatp(x;t)=8�&#x-0.4;掔က:p0+e RrR;p0+e 1 R(R+wr)RrR+w;p0R+wr;(58)withr=kx(0:5;0:5)kandp0=e 3R2(R+)3 R.Weset=1=1,ande=1,andweconsiderthetimeinterval0t3.Fig.4summarizestheerrordataattimet=3forN=64128256,and512andMfac=1and4,with
t=0:25
x.Second-orderconvergenceratesareobservedinthediscreteL1L2,andL1normsforthevelocityeld.Second-orderconvergenceratesarealsoobservedforthepressureinthediscreteL1norm;however,becausethepressureeldisC0butnotC1,onlyrst-orderconvergenceratesareobservedforthepressureinthediscreteL1norm,andintermediateconvergenceratesareobservedintheL2norm.Wealsoconsiderthecaseinwhich\r=0:1,sothattheinitialcongurationoftheshellisnotinequilibrium.Weset=1=0:01,ande=1,yieldingaReynoldsnumberofapproximately50.Weconsiderthetimeinterval0t0:75,whichcorrespondstoapproximatelyonedampedoscillationoftheshell.Fig.5summarizestheerrordataattimet=0:75forN=64128256and512andMfac=1and4,with
t=0:25
x.Essentiallysecond-orderconvergenceratesareobservedinthediscreteL1L2,andL1normsforthevelocityeld.Essentiallysecond-orderconvergenceratesarealsoobservedforthepressureinthediscreteL1norm,whereasonlyrst-orderconvergenceisobservedintheL1norm,andintermediateconvergenceratesareobservedintheL2norm.Convergenceratesforthedeformationaresomewhatlessregular,withnearlysecond-orderconvergenceratesbeingobservedinthediscreteL1normandbetweenrst-andsecond-orderconvergenceratesobservedinthediscreteL2andL1norms.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 14B.E.GRIFFITHANDX.Y.LUO L1erroruL2errorL1errorNpNN 641282565126412825651264128256512641282565126412825651264128256512641282565126412825651264128256512104103102105104103105104103102101103102101103102103102101104103102104103102Figure5.Errorsinu,p,andindiscreteL1,L2,andL1normsforthedynamic(\r=0:1)versionoftheidealizedanisotropicshellproblemofSec.6.1.1.ErrorsforMfac=1areplottedinblue,anderrorsforMfac=4areplottedinred.Referencelineswithslope-2areprovidedfortheuerrordata.Forp,referencelineswithslopes-2,-1.5,and-1areprovidedfortheL1,L2,andL1normdata,respectively.For,referencelineswithslope-2areprovidedfortheL1normdata,andlineswithslopes-1.5areprovidedfortheL2andL1normdata.Noticethatinallcases,accuracyislargelysimilarforbothchoicesofMfac,suggestingthattheschemeyieldsresultsthatarelargelyindependentoftherelativecoarsenessoftheLagrangianmesh.Inparticular,becausetheconvergenceratesforuaresimilarforbothchoicesofMfac,theseresultssuggestthattheschemedoesnotallowleaksatuid-structureinterfaces,evenforLagrangianmeshesthatarerathercoarsecomparedtotheEuleriangrid.6.1.2.OrthotropicshellThesecondelasticitymodelthatweuseforthisproblemisanorthotropicmaterialmodeldenedintermsofthestrain-energyfunctional[28],e(F)=e 2w \r\r\r\r@ @s1(s;t)\r\r\r\r2+\r\r\r\r@ @s2(s;t)\r\r\r\r2!=e 2wF:F;(59)forwhichPe(s;t)=@We @F(s;t)=e wF(s;t):(60)Thisenergyfunctionalandstresstensorcorrespondtoanelasticmaterialcomposedoftwocontinuousfamiliesofbers.Therstfamilyofberswrapstheellipticalshellcircumferentially,andthesecondfamilyiscomposedofradialbersthatareorthogonaltothecircumferentialbers.Becauseonefamilyofbersterminatesalongtheuid-structureinterfaces,therearesingularforcelayersalong@(U;t)thatmustbebalancedbydiscontinuitiesinthepressureandviscousstressTherefore,inthiscasethediscretizeduniedandpartitionedformulationsyielddifferentresults.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD15 L1erroruL2errorL1errorNpN 64128256512641282565126412825651264128256512641282565126412825651281012100101100104103104104Figure6.ErrorsinuandpindiscreteL1,L2,andL1normsforthestatic(\r=0)versionoftheorthotropicshellproblemofSec.6.1.2.ErrorsforMfac=1areplottedinblue,anderrorsforMfac=4areplottedinred;errorsforthepartitionedformulationappearassolidlines,anderrorsfortheuniedformulationappearasdashedlines.Referencelineswithslope-1areprovidedfortheuerrordata.Forp,referencelineswithslopes-1and-0.5areprovidedfortheL1andL2normdata,respectively.(Becausethisproblempossessesdiscontinuitiesinthepressureatuid-structureinterfaces,thepresentmethoddoesnotyieldconvergenceinpinthediscreteL1norm.)Werstset\r=0,whichisanequilibriumcongurationofthecontinuousproblem,sothatu(x;t)0.RequiringR\np(x;t)dx=0,itcanbeshown[28]thatp(x;t)=8�&#x-0.4;掔က:p0+e1 R1 R+
rR;p0+e 1 R(R+wr)+R R+
RrR+w;p0R+wr;(61)withr=kx(0:5;0:5)kandp0=e 33wR+R2(R+)3 R.Weset=1=1,ande=1andweconsiderthetimeinterval0t3.Fig.6summarizestheerrordataattimet=3forN=64128256,and512andMfac=1and4,with
t=0:25
x.First-orderconvergenceratesareobservedforuinallnorms.First-orderconvergenceratesarealsoobservedforpintheL1norm.Becauseppossessesdiscontinuitiesatuid-structureinterfacesforthisproblem,however,thepresentmethodyieldsconvergenceratesof0.5intheL2normanddoesnotconvergeintheL1norm.Wealsoconsiderthecaseinwhich\r=0:1,sothattheinitialcongurationoftheshellisnotinequilibrium.Becausethereisnoexactanalyticsolutionavailableforthisproblem,weestimatetheconvergenceratesinastandardwaybycomputingthenormsofthedifferencesofquantitiescomputedonanN-by-NCartesiangridandcorrespondingLagrangianmesh,andthosecomputedona2N-by-2NCartesiangridandcorrespondingLagrangianmesh.Weset=1=0:01,ande=1,yieldingaReynoldsnumberofapproximately50.Weconsiderthetimeinterval0t1:25,whichcorrespondstoapproximatelyonedampedoscillationoftheshell.Fig.7summarizetheerrordataattimet=0:75forN=64128256,and512andMfac=1and4Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 16B.E.GRIFFITHANDX.Y.LUO L1erroruL2errorL1errorNpNN 641282565126412825651264128256512641282565126412825651264128256512641282565126412825651264128256512104103102104103105104103357101100101100102101103102101103102101Figure7.Errorsinu,p,andindiscreteL1,L2,andL1normsforthedynamic(\r=0:1)versionoftheorthotropicshellproblemofSec.6.1.2.ErrorsforMfac=1areplottedinblue,anderrorsforMfac=4areplottedinred;errorsforthepartitionedformulationappearassolidlines,anderrorsfortheuniedformulationappearasdashedlines.Referencelineswithslope-1areprovidedfortheuerrordata.Forp,referencelineswithslopes-1and-0.5areprovidedfortheL1andL2normdata,respectively.(Becausethisproblempossessesdiscontinuitiesinthepressureatuid-structureinterfaces,thepresentmethoddoesnotyieldconvergenceinpinthediscreteL1norm.)For,referencelineswithslope-1areprovided.with
t=0:25
x.Essentiallyrst-orderconvergenceratesareobservedforuandinallnorms,whereaspexhibitsrst-orderconvergenceinonlytheL1norm.Forthisproblem,wendthattheuniedandpartitionedformulationsyieldsimilaraccuracyinmostcasesforu,andthattheuniedformulationcanoffermodestlybetteraccuracyforBycontrast,thepartitionedformulationofferssignicantlybetteraccuracyforthepressureforrelativelycoarseLagrangianmeshes.Thispropertyappearsalsotoresultinimprovementsinvolumeconservation;seeSec.6.2.6.2.SoftelasticdiscThesetestsdemonstratethevolume-conservationpropertiesofourmethodbyconsideringtheproblemofasoftelasticdiscimmersedinalid-drivencavityow.Thephysicaldomainis\n=[0;1]2,andthephysicalboundaryconditionsareu0alongtheleft(x1=0),right(x1=1),andbottom(x2=0)boundariesof\n,andu(1;0)alongthetop(x2=1)boundaryof\n.Theimmersedstructureisadiscofradius0:2thatisinitiallycenteredabouttheposition(0:6;0:5)Theowinducedbythedrivenlidbringsthestructurenearlyintocontactwiththemovingupperboundaryofthedomain.ThisnearcontactishandledautomaticallybytheIBformulationoftheproblem.Closeto@\n,weuseamodiedregularizeddeltafunctionformulation[19]toensurethatforceandtorqueareconservedandthatvelocityisinterpolatedaccuratelyduringLagrangian-Eulerianinteraction.Noadditionalspecializedmethodsarerequiredbythepresentschemetohandlethiscase.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD17 Figure8.ResultsfromtheelasticdiscproblemofSec.6.2shownatequallyspacedtimesduringtheinterval3t8. Figure9.SimilartoFig.8,buthereshowinganenlargedviewofthestructureattimest=4,5,and6.Weuseanisotropicneo-Hookeanmaterialmodele(F)=e 2F:F;(62)forwhichPe(s;t)=eF:(63)BecausegenerallyPeN60,thesolutionwillgenerallypossessdiscontinuitiesinthepressureandviscousstressatuid-structureinterfaces,andweexpecttheIBmethodtoyieldnobetterthanrst-orderconvergenceratesforthisexample.Thisowalsopossesseswell-knowncornersingularitiesthatalsoacttoreducetheconvergencerateoftheincompressibleowsolver.Althoughitispossibletodevisenumericalschemesthataccuratelytreatthecornersingularitiespresentintheclassicallid-drivencavityow[47],wedonotemploysuchamethodinthiswork.Asinpreviousstudiesofthisproblem[48,49],weset=0:01=1,ande=0:1.Theinitialvelocityisu0,andtheinitialdeformationis(s;0)s.ThephysicaldomainisdiscretizedonanNNCartesiangrid,andtheLagrangiancoordinatedomainisdiscretizedusingasemistructuredmeshofbilinear(Q1)elementsthatisconstructedsothatthecoarsestelementisapproximatelyafactorofMfaccoarserthanthebackgroundEuleriangrid.Weconsiderthetimeinterval0t10duringwhichthediscissubjectedtoslightlymorethanonerotationwithinthecavity.Thestructurebecomesentrainedintheshearingowalongthecavitylidfromt4untilt6,andduringthistimeissubjectedtoextremelylargedeformations.SampleresultsareshowninFigs.8and9.Fig.10showstherelativechangeindiscareafordifferentvaluesofNandMfac.Overthistimeinterval,Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 18B.E.GRIFFITHANDX.Y.LUO percentchange percentchange split,Mfac=1split,Mfac=4unied,Mfac=1unied,Mfac=4tpercentchange 24681024681024681000.020.0400.050.100.10.2Figure10.Percentchangeinareaasafunctionoftimefortheelasticdiscproblem.Resultsareshownfortheuniedandpartitioned(split)formulationsforN=128,256,and512forMfac=1and4.themaximumareachangeyieldedbytheuniedformulationislessthan0.5%forN=64andMfac=4,andisapproximately0:2%forN=64andMfac=1.Forthepartitionedformulation,themaximumareachangeisapproximately0:2%forN=64andbothMfac=1andMfac=4.Witheitherformulation,themaximumvolumechangeconvergestozeroatarst-orderrate.Ingeneral,thepartitionedformulationappearstoyieldsuperiorvolumeconservationforthisproblem,anditsvolume-conservationpropertiesseemtoberelativelyinsensitivetotherelativecoarsenessoftheLagrangianmesh.Usingeitherformulation,volumeerrorsconvergetozeroatessentiallyarst-orderrate.TheseresultscomparefavorablytothoseobtainedbytheIFEmethod,which,evenforrelativelyneLagrangianmeshes,yieldsvolumelossesofupto20%whenappliedtothesametestproblemwithoutusingavolume-conservationalgorithm,andwhichexhibitsvolumelossesofapproximately2.5%whenusingavolume-conservationalgorithm[49].6.3.CollapsiblechannelowThesetestsconsidertheproblemofsteadypressure-drivenowinatwo-dimensionalcollapsiblechannel[50].Theuidpropertiesare=1gm=cm3and=0:01ba
s.ThechannelisoflengthL=40cmandwidthD=1cmandiscenteredintheregion\n=[0;L][0;H]withH=3DInthereferenceconguration,thechannelwallshavethicknessesw=0:01Dandarestraight.AportionoftheupperwalloflengthLf=5cmisexible,andanupstreamportionoftheupperwalloflengthLu=5cmandadownstreamportionoflengthLd=30cmarerigid.Theentirelengthofthelowerwallisrigid.Thetangentialvelocityisclampedtozeroalong@\n.Atsteady-state,theupstreampressureinthechannelispu=0:45Pa,thedownstreampressureisp=0Pa,andtheexternalpressureispe=1:755Pa.SeeFig.11.Near@\n,weagainuseamodiedregularizeddeltafunctionformulation[19]toensurethatforceandtorqueareconservedandthatvelocityisinterpolatedaccuratelyduringLagrangian-Eulerianinteraction.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD19 p=pup=pdp=peLDLuLfLdFigure11.Apartially-collapsiblechanneloflengthLandwidthD.Thechanneliscenteredintheregion\n=[0;L][0;H]withH=3D.AportionoftheupperwalloflengthLfisexible,andanupstreamportionoftheupperwalloflengthLuandadownstreamportionoflengthLdarebothrigid.Theentirelengthofthelowerwallisrigid.Wesettheupstreampressureinthechanneltopu,thedownstreampressuretopd,andtheexternalpressuretope. t(s)position(cm) IB,N=48IB,N=96ALE 02040601.51.61.71.81.92Figure12.Trajectoryofamaterialpointontheexiblesectionofthecollapsiblechannelfortwodifferentgridresolutions.Theypositionofmaterialpoints=(7:5;2:0),whichisinitiallylocatedinthecenteroftheexibleportionoftheupperwall,isplottedasafunctionoftimeforN=48and96.Althoughtherearedifferencesinthedynamics,noticethatthesteady-statecongurationsarethesameforbothgridspacings. Figure13.ComparisonbetweenIBandALEresultsforthecollapsiblechannelbenchmark.Thesteady-statepositionoftheupperwallinthevicinityoftheexibleregionisshownascomputedbytheIB(grey)andALE(black)methods.Despitedifferencesbetweenthetwoschemes,thedeformationspredictedbybothmethodsareessentiallyidentical.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 20B.E.GRIFFITHANDX.Y.LUOIntheexibleportionofthewall,weuseaSt.Venant-Kirchhoffelasticitymodel,forwhiche=e 2tr(E)2+etr(E2);(64)withE=1 2(CI)andC=FTF.Theconstitutiveparametersaregivenbye=1 2Ee=(1+e)ande=(Eee)=((1+e)(12e)),inwhichEe=71:8kPaistheYoung'smodulusofthematerialande=0:45isthePoissonratio.Rigidityconstraintsareenforcedbyafeedback-forcingapproach.Specically,anadditionalLagrangianconstraintforceF(s;t)isaddedtotheelasticforcedensity.TheadditionalconstraintforceisdeterminedviaF(s;t)=e(s)(s(s;t));(65)inwhiche(s)�0isapenaltyparameterthatissettozerowithintheexiblepartofthechannelwall.Noticethatase!1(s;t)sWecomparetheresultsproducedbytheuniedversionofthepresentIBmethodtoresultsthatwereobtainedbyanarbitraryLagrangian-Eulerian(ALE)method[50,51].AnimportantdifferencebetweenthephysicalproblemsconsideredbytheALEmethodandthepresentIBmethodisthatthephysicalsystemtreatedbyALEmethodneglectstheinertiaassociatedwithanymaterialexteriortothechannelandinsteadtreatstheexteriorasaconstant-pressurereservoir.Bycontrast,theIBmodelconsidersthechanneltobetotallyimmersedinuid.(Inprinciple,neglectingtheinertiaoftheexteriorregionispossiblewithintheframeworkoftheIBmethod,butdoingsorequiresavariable-densityincompressibleNavier-Stokessolver.Ourpresentimplementationonlyprovidesaconstant-densitysolver.)Consequently,thedynamicsproducedbythetwophysicalmodelswillgenerallybedifferent.Despitethesedifferences,however,bothphysicalmodelsadmitidenticalsteady-stateequilibriumcongurations.Wetherebycompareonlythestaticcongurationsproducedbythetwonumericalmethods.Wediscretize\nusinganL HN-by-NCartesiangridwithN=48and96,andwediscretizeeachchannelwallusingasinglelayerofbilinearelements,eachoflength
xandwidthwinthereferenceconguration.NoticethatforeithervalueofNw
x.Weuse
t=4:05e-3
x,andintherigidportionsofthechannelwalls,wesete=1:0e3
x=
t2,whichisapproximatelythelargestvalueofethatisstableforthephysicalandnumericalparametersconsidered.Weincreasethedrivingandloadingpressuresfromzerooveratimeintervalof1s.Thesystemundergoeslarge,oscillatorydeformationsduringtheinitialpartofthesimulation,butbyt=100sthesystemhasreachedequilibrium.Atsteadystate,Re80.TheIBmethodyieldsconvergedsteady-statedeformationsforN=48;seeFig.12.ThenalcongurationsgeneratedbytheIBandALEmethodsareshowninFig.13.Itcanbeseenthatthedeformationsgeneratedbythetwomethodsarevirtuallyidentical.Themaximumydisplacementsdeterminedbythetwomethodsalsoagreetowithin0:5%:intheIBmodel,themaximumydisplacementofthecenterlineoftheupperwallis0:3654cm,whereasintheALEmodel,itis0:3641cm.Recallthatthethicknessofthebeamisw=0:01cm.Giventhedifferencesbetweenthephysicalmodelsandnumericalschemes,theagreementbetweentheIBandALEmethodsisexcellent.6.4.Turek-Hronuid-structureinteractionbenchmarkThissectionpresentspreliminaryresultsobtainedbyourIBmethodfortheTurek-Hronuid-structureinteractionbenchmarkproblem[52],whichconsiderschannelowpastastructurecomposedofarigidcylinderandatrailingelasticbar;seeFig.14.WeusetheversionofthisproblematRe=200,forwhichtheuidpropertiesare=1000kg=2and=1Pa
s,andforwhichthemassdensityoftheelasticbaristhesameasthatoftheuid.Thechannelis\n=[0;L][0;H]withL=2:5mandH=0:41m.No-slipandno-penetrationboundaryconditions(u0)areusedalongthebottom(x2=0)andtop(x2=H)boundariesof\n.Aparabolicinowproleoftheformu1(x2;t)=1:5U(t)y(Hy) (H=2)2,withU(t)=(2:0=s)(1 2(1cos(t=2))ift2;1otherwise;(66)Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD21 Figure14.SampleresultsfromtheIBversionoftheTurek-Hronuid-structureinteractionbenchmarkproblem[52]atRe=200.Theuidvorticityrangesfrom30m=s2(blue)to+30m=s2(red).Thestructureisshowninblack. xdisplacement(m)ydisplacement(m)t(s)lift(N)t(s)drag(N) 19.619.719.819.919.619.719.819.919.619.719.819.919.619.719.819.9420440460480-1000100-0.0200.02-4e-3-3e-3-2e-3-1e-30 Figure15.xandydisplacementsofthecontrolpointAandliftanddragforcesfortheuid-structureinteractionbenchmarkproblemforN=328.isimposedalongthex1=0boundary,andzero-tangential-slipandzero-normal-tractionboundaryconditionsareimposedatx1=LCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 22B.E.GRIFFITHANDX.Y.LUO Mfac conventional 0.5124 IB N64 0.3680.1560.1000.089 0.312128 1.6300.6490.3830.318 1.525256 6.4412.9531.7621.416 6.065TableI.Averagewall-clocktime(inseconds)pertimestepforthedynamicversionoftheproblemofSec.6.1.2.WealsoprovidetimingsforaversionofthismethodsimilartothestandardIBmethod,inwhichweuseadiagonally-lumpedmassmatrix,setMfac=0:5,andusethetrapezoidalruletoconstructSandR.ThedischascenterC=(0:2;0:2)andradiusR=0:05mandisheldrigidviaapenaltyapproachlikethatusedinSec.6.3.Theelasticbarisattachedtothebackofthediscandhaslengthl=0:35mandwidthh=0:02m.ASt.Venant-Kirchhoffmodelisusedtodescribetheelasticityofthebarwithe=2:0e6MPaande=0:4.AcontrolpointAwithinitialposition(0:6;0:2)isattachedtothecenterofthetrailingedgeofthebar.ThereisadifferencebetweenthebenchmarkproblemdenedbyTurekandHron[52]andthepresentIBmodel:Intheoriginalproblemspecication,theelasticbariscompressible.ThepresentIBformulation,however,requiresthatthecoupleduid-structuresystembeincompressible.Consequently,inthepresentwork,theexiblebaristreatedasanincompressiblematerial.(TreatingacompressiblestructurewithintheframeworkoftheIBmethodappearstorequireemployingacompressibleNavier-Stokessolverandadoptingacompressibleformulationofboththeuidandthestructure.Suchanextensionofthepresentformulationisbeyondthescopeofthiswork.)Theextenttowhichthisproblemissensitivetothisdifferencehasnotyetbeenfullycharacterized;however,ourresultssuggestthatthedifferencesbetweenthecompressibleandincompressibleversionsofthisbenchmarkproblemmayberelativelyminor(seebelow).Wediscretize\nusinganL HN-by-NCartesiangrid.Thestructureisthesuperpositionofacirculardiscandanelasticbar,whicharebothdiscretizedusingbilinearelementswithedgelengthsthatarenolargerthan
x.Weset
t=7:8125e-3
xandusee=1:0e5
x=
t2,whichisapproximatelythelargestvalueofethatisstableforthesephysicalandnumericalparameters.SampleresultsobtainedusingtheuniedformulationwithN=328areshowninFigs.14and15.Overthetimeinterval[19:5s;20s],thexdisplacementofcontrolpointAhasmean-2:18e-3mandamplitude2:19e-3m;theydisplacementhasmean1:34e-3mandamplitude3:16e-2m;theliftforcehasmean2:80Nandamplitude1:84e+2N;andthedragforcehasmean4:47e+2Nandamplitude40:9N.ThefrequencyofthexdisplacementsofAisapproximately11:0=s,andthefrequencyoftheydisplacementsisapproximately5:4=s.Althoughtheredoesnotyetappeartobeareferencesolutionforthisbenchmarkproblem,thesevaluesareallwithintherangeofcomparableresultsobtainedbyvariousothernumericalmethodsthathavebeenappliedtothisproblem[52,53].Noticethattheseearlierresultsallconsideredtheelasticbartobecompressible.6.5.TimingresultsWeconcludebypresentingpreliminarytimingdataobtainedforthedynamicversionoftheproblemofSec.6.1.2.Weremark,however,thatwehavenotyetfullyoptimizedourimplementationofthismethod.Forthesetests,weperformtherst64timestepsforvariouschoicesofNandMfacusingthepartitionedversionofthemethod.WealsosimulatetheconventionalIBmethodwithinourimplementationbyusingtheuniedformulation,usingadiagonally-lumpedmassmatrix,settingMfac=0:5,andusingthetrapezoidalruletoconstructSandR.Werecordthewall-clocktimerequiredusingtheCfunctionclock gettime()withtheCLOCK REALTIMEtimer.Timingswereperformedusingasinglecoreofacomputeserverequippedwithfourquad-core2.3GHzAMDBarcelonaCPUsand32GBRAM.ResultsaresummarizedinTableI.NoticethatlargervaluesofMfac,whichcorrespondtocoarserLagrangianmeshes,alwaysyieldsmallerruntimes.AlsonoticethattheLagrangian-EulerianinteractionoperatorsintroducedinthisworkgenerallyyieldbetterperformancethanthoseoftheconventionalIBapproachexceptforthenestLagrangianmeshes.ThesepreliminaryresultssuggestthatthereductionsinthenumberofLagrangiandegreesCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD23offreedompermittedbytheLagrangian-Eulerianinteractionoperatorsintroducedhereincanyieldsignicantimprovementsinperformance.Weexpectthatthesedifferenceswouldbeevenmorepronouncedinthree-dimensionalmodels.7.CONCLUSIONSInthispaper,wehavedescribedaversionoftheIBmethodforproblemsofuid-structureinteractionthatusesgeneralincompressibleelasticitymodelswithunstructuredFEdiscretizationswhileretainingaCartesiangridnitedifferenceschemefortheincompressibleNavier-Stokesequations.AfeatureofthepresentmethodisthatitusesstandardmethodsforboththeLagrangianandEulerianequations.Inpractice,itshouldbefeasibletousethisapproachtocoupleexistingstructuralanalysisandincompressibleowcodesbypassingforcesfromthestructuralcodetotheuidsolver,andbypassingvelocitiesfromtheuidsolverbacktothestructuralcode.Althoughweonlyconsidercasesinwhichtheelasticityofthestructureisdescribedbyahyperelasticconstitutivelawdenedintermsofastrain-energyfunctional,ournumericalalgorithmdoesnotrelyontheavailabilityofsuchanenergyfunctionalandrequiresonlyaroutinetoevaluatethematerialstresstensor.OurmethodisbasedonaformulationofthecontinuousIBequationsintroducedbyBofetal.[28].WeconsidertwodifferentstatementsoftheseequationsthateachuseaweakformulationoftheLagrangianequationsofnonlinearelasticity.Oneoftheseformulationstreatstheinternalandboundarystresseswithintheimmersedelasticstructureusingasingle,unied,volumetricelasticforcedensity.Theotherpartitionsthesestressesintoainternalelasticforcedensitysupportedontheinteriorofthestructure,andatransmissionelasticforcedensitysupportedonuid-structureinterfaces.Bothweakformulationsaredemonstratedtoyieldsimilarconvergencerates,butthepartitionedformulationisseentoyieldhigheraccuracyforLagrangianmeshesthatarerelativelycoarseincomparisontothebackgroundEuleriangrid.Thepartitionedformulationisalsoseentoyieldsuperiorvolumeconservationincomparisontotheuniedformulation.AlimitationofthepartitionedformulationisthatitdoesnotsatisfyadiscretepoweridentitythatimpliesthatenergyisconservedduringLagrangian-Eulerianinteraction.Suchapoweridentityissatisedbytheuniedformulation,andmaybenecessarytoobtainanunconditionallystableimplicittime-steppingscheme[54].AkeycontributionofthispaperisthatitintroducesanewapproachtoLagrangian-EulerianinteractionthatovercomesalongstandinglimitationoftheIBmethod,namelytherequirementthattheLagrangianmeshbeapproximatelytwiceasneasthebackgroundEuleriangridtoavoidleaksatuid-structureinterfaces.NumericalexamplesdemonstratethatourschemepermitstheuseofLagrangianmeshesthatareatleastfourtimesascoarseasthebackgroundEuleriangridwithoutleak.WespeculatethattheLagrangianmeshmaybearbitrarilycoarsewithrespecttothebackgroundEuleriangridwithoutleakaslongasthenumericalquadratureschemeusedtodiscretizetheLagrangian-Eulerianinteractionequationshasasufcientlydensenetofquadraturepoints.Intwospatialdimensions,ourschemehasbeendemonstratedtoreducethenumberofLagrangiandegreesoffreedombyafactorofatleast82=64incomparisontotheconventionalIBmethodforselectedproblems.Inthreespatialdimensions,weexpectthatwewouldbeabletoseesavingsontheorderof83=512ormoreinsomecases.ThesereductionsinthenumberofLagrangiandegreesoffreedomarealsoshowntotranslateintoimprovedwall-clocktimings.Inpractice,whetheritispossibletouseacoarsestructuralmeshisclearlyproblemdependent.TheadvantageofthepresentapproachisthatitenablestheanalysttoemployspatialdiscretizationsthataretailoredtotherequirementsofthestructuralmodelratherthandictatedbythebackgroundEuleriandiscretization.Finally,wespeculatethatthepartitionedformulationcouldbeusefulinconstructinghigher-orderversionsoftheIBmethod.Inparticular,thepartitionedformulationiswellsuitedfordevelopingahybridapproachinwhichtheIBmethodisusedtotreatthevolumetricinternalforcedensity,andanothermethodisusedtotreatthesingulartransmissionforcedensity.Forinstance,becausethetransmissionforcedensityofthepartitionedformulationisdenedonaclosedsurface,itmaybepossibletotreatitwithhigher-orderaccuracybyaversionoftheimmersedinterfacemethodCopyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 24B.E.GRIFFITHANDX.Y.LUO[55–60].Suchanextensionofthismethodcouldyieldafullysecond-orderaccurategeneralizationoftheIBmethodforproblems,likethoseconsideredinthiswork,inwhichtheimmersedstructureisofcodimensionzerowithrespecttotheuid.ACKNOWLEDGEMENTSThisworkwassponsoredinpartbyanawardfromtheRoyalSocietyofEdinburgh.B.E.G.acknowledgesresearchsupportfromtheAmericanHeartAssociation(AHAaward10SDG4320049)andtheNationalScienceFoundation(USNSFawardsDMS-1016554andOCI-1047734).HealsogratefullyacknowledgesdiscussionsofthisworkwithCharlesPeskinandDavidMcQueen.X.Y.L.acknowledgesresearchsupportfromtheEngineeringandPhysicalSciencesResearchCouncil(UKEPSRCawardsEP/G015651andEP/I029990).ComputationswereperformedatNewYorkUniversityusingcomputerfacilitiesfundedinlargepartbyagenerousdonationbySt.JudeMedical,Inc.REFERENCES1.C.S.Peskin.Flowpatternsaroundheartvalves:Adigitalcomputermethodforsolvingtheequationsofmotion.PhDthesis,AlbertEinsteinCollegeofMedicine,1972.2.C.S.Peskin.Numericalanalysisofbloodowintheheart.JComputPhys,25(3):220–252,1977.3.C.S.Peskin.Theimmersedboundarymethod.ActaNumer,11:479–517,2002.4.G.IaccarinoandR.Verzicco.Immersedboundarytechniqueforturbulentowsimulations.ApplMechRev,56(3):331–347,2003.5.R.MittalandG.Iaccarino.Immersedboundarymethods.AnnuRevFluidMech,37:239–261,2005.6.L.A.MillerandC.S.Peskin.Whenvorticesstick:anaerodynamictransitionintinyinsectight.JExpBiol,207(17):3073–3088,2004.7.L.A.MillerandC.S.Peskin.Acomputationaluiddynamicsof`clapanding'inthesmallestinsects.JExpBiol,208(2):195–212,2005.8.A.L.FogelsonandR.D.Guy.Immersed-boundary-typemodelsofintravascularplateletaggregation.ComputMethApplMechEng,197(25–28):2087–2104,2008.9.X.Yang,R.H.Dillon,andL.J.Fauci.Anintegrativecomputationalmodelofmulticiliarybeating.BullMathBiol,70(4):1192–1215,2008.10.C.-Y.HsuandR.Dillon.A3Dmotilerod-shapedmonotrichousbacterialmodel.BullMathBiol,71(5):1228–1263,2009.11.L.A.MillerandC.S.Peskin.Flexibleclapandingintinyinsectight.JExpBiol,212(19):3076–3090,2009.12.C.S.PeskinandD.M.McQueen.Fluiddynamicsoftheheartanditsvalves.InH.G.Othmer,F.R.Adler,M.A.Lewis,andJ.C.Dallon,editors,CaseStudiesinMathematicalModeling:Ecology,Physiology,andCellBiology,pages309–337.Prentice-Hall,EnglewoodCliffs,NJ,USA,1996.13.D.M.McQueenandC.S.Peskin.Shared-memoryparallelvectorimplementationoftheimmersedboundarymethodforthecomputationofbloodowinthebeatingmammalianheart.JSupercomput,11(3):213–236,1997.14.J.D.LemmonandA.P.Yoganathan.Three-dimensionalcomputationalmodelofleftheartdiastolicfunctionwithuid-structureinteraction.JBiomechEng,122(2):109–117,2000.15.J.D.LemmonandA.P.Yoganathan.Computationalmodelingofleftheartdiastolicfunction:Examinationofventriculardysfunction.JBiomechEng,122(4):297–303,2000.16.D.M.McQueenandC.S.Peskin.Athree-dimensionalcomputermodelofthehumanheartforstudyingcardiacuiddynamics.ComputGraph,34(1):56–60,2000.17.D.M.McQueenandC.S.Peskin.Heartsimulationbyanimmersedboundarymethodwithformalsecond-orderaccuracyandreducednumericalviscosity.InH.ArefandJ.W.Phillips,editors,MechanicsforaNewMillennium,Proceedingsofthe20thInternationalConferenceonTheoreticalandAppliedMechanics(ICTAM2000).KluwerAcademicPublishers,2001.18.B.E.Grifth,R.D.Hornung,D.M.McQueen,andC.S.Peskin.Anadaptive,formallysecondorderaccurateversionoftheimmersedboundarymethod.JComputPhys,223(1):10–49,2007.19.B.E.Grifth,X.Luo,D.M.McQueen,andC.S.Peskin.Simulatingtheuiddynamicsofnaturalandprostheticheartvalvesusingtheimmersedboundarymethod.IntJApplMech,1(1):137–177,2009.20.B.E.Grifth.Immersedboundarymodelofaorticheartvalvedynamicswithphysiologicaldrivingandloadingconditions.IntJNumerMethBiomedEng,28(3):317–345,2012.21.X.Y.Luo,B.E.Grifth,X.S.Ma,M.Yin,T.J.Wang,C.L.Liang,P.N.Watton,andG.M.Bernacca.Effectofbendingrigidityinadynamicmodelofapolyurethaneprostheticmitralvalve.BiomechModelMechanobiology,11(6):815–827,2012.22.C.S.Peskin.TheImmersedBoundaryMethod.III.EnergyFunctionsfortheRepresentationofImmersedElasticBoundariesandMaterials,2007.Handwrittenlecturenotesavailablefromhttp://math.nyu.edu/faculty/peskin/ib\_lecture\_notes/index.html.23.D.DevendranandC.S.Peksin.Anenergy-basedimmersedboundarymethodforincompressibleviscoelasticity.JComputPhys,231(14):4613–4642,2012.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme FINITEDIFFERENCE/FINITEELEMENTIMMERSEDBOUNDARYMETHOD2524.A.J.Gil,A.ArranzCarre˜no,J.Bonet,andO.Hassan.TheImmersedStructuralPotentialMethodforhaemodynamicapplications.JComputPhys,229(22):8613–8641,2010.25.L.Zhang,A.Gerstenberger,X.Wang,andW.K.Liu.Immersedniteelementmethod.ComputMethApplMechEng,193(21–22):2051–2067,2004.26.W.K.Liu,Y.Liu,D.Farrell,L.Zhang,X.S.Wang,Y.Fukui,N.Patankar,Y.Zhang,C.Bajaj,J.Lee,J.Hong,X.Chen,andH.Hsu.Immersedniteelementmethodanditsapplicationstobiologicalsystems.ComputMethApplMechEng,195(13–16):1722–1749,2006.27.L.T.ZhangandM.Gay.Immersedniteelementmethodforuid-structureinteractions.JFluidStruct,23(6):839–857,2007.28.D.Bof,L.Gastaldi,L.Heltai,andC.S.Peskin.Onthehyper-elasticformulationoftheimmersedboundarymethod.ComputMethApplMechEng,197(25–28):2210–2231,2008.29.L.HeltaiandF.Costanzo.Variationalimplementationofimmersedniteelementmethods.ComputMethApplMechEng.Toappear.30.B.E.GrifthandC.S.Peskin.Ontheorderofaccuracyoftheimmersedboundarymethod:Higherorderconvergenceratesforsufcientlysmoothproblems.JComputPhys,208(1):75–105,2005.31.IBAMR:Anadaptiveanddistributed-memoryparallelimplementationoftheimmersedboundarymethod.http://ibamr.googlecode.com.32.B.E.Grifth.Onthevolumeconservationoftheimmersedboundarymethod.CommunComputPhys,12:401–432,2012.33.B.E.Grifth.AnaccurateandefcientmethodfortheincompressibleNavier-Stokesequationsusingtheprojectionmethodasapreconditioner.JComputPhys,228(20):7565–7595,2009.34.W.J.Rider,J.A.Greenough,andJ.R.Kamm.Accuratemonotonicity-andextrema-preservingmethodsthroughadaptivenonlinearhybridizations.JComputPhys,225(2):1827–1848,2007.35.P.ColellaandP.R.Woodward.Thepiecewiseparabolicmethod(PPM)forgas-dynamicalsimulations.JComputPhys,54(1):174–201,1984.36.Y.Saad.Aexibleinner-outerpreconditionedGMRESalgorithm.SIAMJSciComput,14(2):461–469,1993.37.SAMRAI:StructuredAdaptiveMeshRenementApplicationInfrastructure.http://www.llnl.gov/CASC/SAMRAI.38.R.D.HornungandS.R.Kohn.ManagingapplicationcomplexityintheSAMRAIobject-orientedframework.ConcurrencyComputPractEx,14(5):347–368,2002.39.R.D.Hornung,A.M.Wissink,andS.R.Kohn.ManagingcomplexdataandgeometryinparallelstructuredAMRapplications.EngComput,22(3–4):181–195,2006.40.S.Balay,K.Buschelman,W.D.Gropp,D.Kaushik,M.G.Knepley,L.C.McInnes,B.F.Smith,andH.Zhang.PETScWebpage,2009.http://www.mcs.anl.gov/petsc.41.S.Balay,K.Buschelman,V.Eijkhout,W.D.Gropp,D.Kaushik,M.G.Knepley,L.C.McInnes,B.F.Smith,andH.Zhang.PETScusersmanual.TechnicalReportANL-95/11-Revision3.0.0,ArgonneNationalLaboratory,2008.42.S.Balay,V.Eijkhout,W.D.Gropp,L.C.McInnes,andB.F.Smith.Efcientmanagementofparallelisminobjectorientednumericalsoftwarelibraries.InE.Arge,A.M.Bruaset,andH.P.Langtangen,editors,ModernSoftwareToolsinScienticComputing,pages163–202.Birkh¨auserPress,1997.43.hypre:Highperformancepreconditioners.http://www.llnl.gov/CASC/hypre.44.R.D.FalgoutandU.M.Yang.hypre:alibraryofhighperformancepreconditioners.InP.M.A.Sloot,C.J.K.Tan,J.J.Dongarra,andA.G.Hoekstra,editors,ComputationalScience-ICCS2002PartIII,volume2331ofLectureNotesinComputerScience,pages632–641.Springer-Verlag,2002.AlsoavailableasLLNLTechnicalReportUCRL-JC-146175.45.libMesh:C++FiniteElementLibrary.http://libmesh.sourceforge.net.46.B.Kirk,J.W.Peterson,R.H.Stogner,andG.F.Carey.libMesh:AC++libraryforparalleladaptivemeshrenement/coarseningsimulations.EngComput,22(3–4):237–254,2006.47.O.BotellaandR.Peyret.Benchmarkspectralresultsonthelid-drivencavityow.ComputFluid,27(4):421–433,1998.48.H.Zhao,J.B.Freund,andR.D.Moser.Axed-meshmethodforincompressibleow-structuresystemswithnitesoliddeformations.JComputPhys,227(6):3114–3140,2008.49.X.WangandL.T.Zhang.Interpolationfunctionsintheimmersedboundaryandniteelementmethods.ComputMech,45(4):321–334,2010.50.H.F.Liu,X.Y.Luo,andZ.X.Cai.Stabilityandenergybudgetofpressure-drivencollapsiblechannelows.JFluidMech.(inpress).51.Z.X.CaiandX.Y.Luo.Auidbeammodelforowincollapsiblechannels.JFluidStruct,17(1):123–144,2003.52.S.TurekandJ.Hron.Proposalfornumericalbenchmarkingofuid-structureinteractionbetweenanelasticobjectandlaminarincompressibleow.InH.J.BungartzandM.Sch¨afer,editors,Fluid-StructureInteraction:Modelling,Simulation,Optimization,number53inLectureNotesinComputationalScienceandEngineering,pages371–385.Springer,2006.53.S.Turek,J.Hron,M.Razzaq,H.Wobker,andM.Sch¨afer.Numericalbenchmarkingofuid-structureinteraction:Acomparisonofdifferentdiscretizationandsolutionapproaches.Technicalreport,Fakult¨atf¨urMathematik,TUDortmund,2010.ErgebnisberichtedesInstitutsf¨urAngewandteMathematik,Nummer405.54.E.P.Newren,A.L.Fogelson,R.D.Guy,andR.M.Kirby.Unconditionallystablediscretizationsoftheimmersedboundaryequations.JComputPhys,222(2):702–719,2007.55.R.J.LeVequeandZ.Li.ImmersedinterfacemethodsforStokesowwithelasticboundariesorsurfacetension.SIAMJSciComput,18(3):709–735,1997.56.Z.-L.LiandM.-C.Lai.TheimmersedinterfacemethodfortheNavier-Stokesequationswithsingularforces.JComputPhys,171(2):822–842,2001.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme 26B.E.GRIFFITHANDX.Y.LUO57.L.LeeandR.J.LeVeque.AnimmersedinterfacemethodforincompressibleNavier-Stokesequations.SIAMJSciComput,25(3):832–856,2003.58.S.XuandZ.J.Wang.Animmersedinterfacemethodforsimulatingtheinteractionofauidwithmovingboundaries.JComputPhys,216(2):454–493,2006.59.S.XuandZ.J.Wang.Systematicderivationofjumpconditionsfortheimmersedinterfacemethodinthree-dimensionalowsimulation.SIAMJSciComput,27(6):1948–1980,2006.60.S.XuandZ.J.Wang.A3Dimmersedinterfacemethodforuid-solidinteraction.ComputMethApplMechEng,197:2068–2086,2008.Copyrightc\r0000JohnWiley&Sons,Ltd.Int.J.Numer.Meth.Engng(0000)Preparedusingnmeauth.clsDOI:10.1002/nme