Numerical simulation of fluid-structure interaction with application to aneurysm hemodynamics

15 Pages·2009·1.88 MB·English

Fluid-StructureInteraction.
Theory,NumericsandApplications
pp.283–294
HerrschingamAmmersee,29.9.-1.10.2008
Numerical simulation of fluid-structure interaction with application to
aneurysm hemodynamics
M.Razzaq,S.Turek,J.Hron,J.F.Acker,F.Weichert,M.Wagner,I.Q.Grunwald,C.Roth,andB.F.Romeike
Asanexampleforfluid-structureinteractioninbiomedicalproblems,theinfluenceofendovascularstentimplan-
tationontocerebralaneurysm hemodynamicsis numericallyinvestigated. The aim is to study the interactionof
theelasticwallsoftheaneurysmwiththegeometricalshapeoftheimplantedstentstructureforprototypical2D
configurations. This study can be seen as a basic step towards the understandingof the resulting complex flow
phenomenaso thatin future aneurysm ruptureshallbe suppressed by an optimalsetting forthe implantedstent
geometry.Fromthemathematicalside,numericaltechniquesforsolvingtheproblemoffluid-structureinteraction
withanelasticmaterialinalaminarincompressibleviscousflowaredescribed.AnArbitraryLagrangian-Eulerian
(ALE)formulationisemployedinafullycoupledmonolithicway,consideringtheproblemasonecontinuum.The
mathematicaldescription andthe numericalschemes are designedin such a way thatmore complicatedconsti-
tutiverelations(andmorerealisticforbiomechanicsapplications)forthefluidaswellasthestructuralpartcan
beeasilyincorporated. Weutilizethewell-knownQ P finiteelementpairfordiscretizationinspacetogainhigh
2 1
accuracyandperformastime-steppingthe2ndorderCrank-Nicholson,resp.,Fractional-Step-q -schemeforboth
solidandfluidparts. TheresultingnonlineardiscretizedalgebraicsystemissolvedbyaNewtonmethodwhichap-
proximatestheJacobianmatricesbythedivideddifferencesapproach,andtheresultinglinearsystemsaresolved
by iterative solvers, preferably of Krylov-multigrid type. Preliminary results for the stent-assisted occlusion of
cerebralaneurysmarepresented. Sincetheseresultsarecurrentlyrestrictedto2Dconfigurations,theaimisnot
to predict quantitatively the complex interaction mechanisms between stents and elastic walls of the aneurysm,
buttoanalysequalitativelythebehaviourofthe elasticityofthewallsvs. thegeometricaldetailsofthestentfor
prototypicalflowsituations.
1 Introduction
Inthiscontribution,weconsiderthegeneralproblemofviscousflowinteractingwithanelasticbodywhichisbeing
deformedby the fluid action. Such a problem is of great importancein many real life applications, and typical
examplesofthistypeofproblemaretheareasofbiomedicalfluidswhichincludetheinfluenceofhemodynamic
factorsinbloodvessels, cerebralaneurysmhemodynamics,jointlubricationanddeformablecartilageandblood
flow interactionwith elastic veins(Appanaboyinaet al., 2008), (Valenciaet al., 2008), (Fernandezet al., 2008),
(Tezduyaretal.,2007),(Tezduyaretal.,2008).Thetheoreticalinvestigationoffluid-structureinteractionproblems
iscomplicatedbytheneedofa mixeddescriptionforbothparts: While forthesolidpartthenaturalviewisthe
material(Lagrangian)description,forthefluiditisusuallythespatial(Eulerian)description. Inthecaseoftheir
combinationsomekindofmixeddescription(usuallyreferredtoastheArbitraryLagrangian-Euleriandescription
or ALE) has to be used which brings additionalnonlinearity into the resulting equations (see (Hron and Turek,
2006b)).
Thenumericalsolutionoftheresultingequationsofthefluid-structureinteractionproblemposesgreatchallenges
sinceitincludesthefeaturesofstructuralmechanics,fluiddynamicsandtheircoupling.Themoststraightforward
solutionstrategy, mostly used in the available software packages(see for instance (Hron et al., 2002)), is to de-
coupletheproblemintothefluidpartandsolidpart,foreachofthosepartsusingsomewellestablishedmethod
ofsolution;thentheinteractionprocessisintroducedasexternalboundaryconditionsineachofthesubproblems.
Thishastheadvantagethattherearemanywelltestednumericalmethodsforbothseparateproblemsoffluidflow
andelasticdeformation,whileontheotherhandthetreatmentoftheinterfaceandtheinteractionisproblematic
dueto highstiffnessand sensitivity. In contrast, the monolithicapproachdiscussed here treats the problemas a
283
singlecontinuumwiththecouplingautomaticallytakencareofasinternalinterface.
Besideashortdescriptionoftheunderlyingnumericalaspectsregardingdiscretizationandsolutionprocedurefor
thismonolithicapproach(see(Razzaqetal.,2008),(HronandTurek,2006a)),weconcentrateonprototypicalnu-
mericalstudiesfor2Daneurysmconfigurations. Thecorrespondingparametrizationwasbasedonabstractionsof
biomedicaldata(i.e.,cutplanesof3DspecimensfromNewZealandwhiterabbitsaswellascomputertomographic
andmagneticresonanceimagingdataofhumanneurocrania).Inourstudies,weallowthewallsoftheaneurysmto
beelasticandhencedeformingwiththeflowfieldinthevessel. Moreover,weexamineseveralconfigurationsfor
stentgeometrieswhichclearlyinfluencetheflowbehaviorinsideoftheaneurysmsuchthataverydifferentelastic
displacementofthewallsisobservedtoo. Wedemonstratethateithertheelasticmodelingoftheaneurysmwalls
as well as the properdescription of the geometricaldetails of the shape of the aneurysmand particularlyof the
stentsisofgreatimportanceifthecomplexinteractionbetweenstructureandfluidshallbequantitativelyanalyzed
in future, especially in view of more realistic blood flow models and anisotropicconstitutive laws of the elastic
walls.
2 Fluid-structureinteractionproblemformulation
Thegeneralfluid-structureinteractionproblemconsistsofthedescriptionofthefluidandsolidfields,appropriate
interfaceconditionsat the interface andconditionsfor the remainingboundaries,respectively. In this paper, we
considerthe flow of an incompressibleNewtonianfluid interactingwith an elastic solid. We denote the domain
occupiedbythefluidbyW bandthesolidbyW satthetimet∈[0,T]. LetG 0=W¯b∩W¯sbethepartoftheboundary
t t t t t
where the elastic solid interacts with the fluid. In the following, the description for both fields fields and the
interfaceconditionsareintroduced. Furthermore,discretizationaspectsandsolutionproceduresarepresentedin
thenextsection.
2.1 Constitutiverelationsforthefluid
ThefluidisconsideredtobeNewtonian,incompressibleanditsstateisdescribedbythe velocityand pressure
fieldsvb, pbrespectively.Theconstantdensityofthefluidisr bandthekinematicviscosityisdenotedbyn b. The
balanceequationsare:
Dvb
r b =divs b, divvb=0 in W b (1)
Dt t
Inordertosolvethebalanceequationsweneedtospecifytheconstitutiverelationsforthestresstensors. Forthe
fluidweusetheincompressibleNewtonianrelation
s b=−pbI+m ((cid:209) vb+((cid:209) vb)T), (2)
where m representsthedynamicviscosityofthefluidand pb istheLagrangemultipliercorrespondingto thein-
compressibilityconstraintin(1). Thematerialtimederivativedependsonthechoiceofthereferencesystem.There
arebasically3alternativereferencesystems: theEulerian,theLagrangian,andtheArbitraryLagrangian-Eulerian
formulation. Themostcommonlyuseddescriptionforthefluid-structureinteractionistheALEdescription. For
the ALE formulationpresentedin thispaper, correspondingdiscretizationtechniquesare discussed in section 3.
Letusremarkthatalsononnewtonianflowmodelscanbeusedformodelingbloodflow,forinstanceofPowerLaw
typeorevenincludingviscoelasticeffects(see(Damaniketal.,2008))whichisplannedforfutureextensions.
2.2 Constitutiverelationsforthestructure
The structure is assumed to be elastic and compressible. Its configurationis described by the displacementus,
withvelocityfieldvs= ¶ us. Thebalanceequationsare:
¶ t
¶ vs
r s +r s((cid:209) vs)vs=divs s+r sg, in W s. (3)
¶ t t
WritteninthemorecommonLagrangiandescription,i.e.withrespecttosomefixedreference(initial)stateW s,we
have
¶ 2us
r s =div(Js sF−T)+r sg, in W s. (4)
¶ t2
284
Theconstitutiverelationsforthestresstensorsforthecompressiblestructurearepresented,however,alsoincom-
pressiblestructurescanbehandledinthesameway(see(HronandTurek,2006b)).Thedensityofthestructurein
theundeformedconfigurationisr s. Thematerialelasticityischaracterizedbyasetoftwoparameters,thePoisson
ration sandtheYoungmodulusE. Alternatively,thecharacterizationisdescribedbytheLame´coefficientsl sand
theshearmodulusm s. Theseparameterssatisfythefollowingrelations
l s m s(3l s+2m 2)
n s= E = (5)
2(l s+m s) (l s+m s)
E n sE
m s= l s= , (6)
2(1+n s) (1+n s)(1−2n s)
where n s =1/2 for a incompressible and n s <1/2 for a compressible structure. In the large deformationcase
it is common to describe the constitutive equation using a stress-strain relation based on the Green Lagrangian
straintensorE andthe2.Piola-KirchhoffstresstensorS(E)asafunctionofE.The2.Piola-Kirchhoffstresscanbe
obtainedfromtheCauchystresss s as
Ss=JF−1s sF−T, (7)
andtheGreen-LagrangetensorE as
1
E = (FTF−I). (8)
2
Inthispaper,thematerialisspecifiedbygivingtheCauchystresstensors s bythefollowingconstitutivelawfor
theSt.Venant-Kirchhoffmaterialforsimplicity
1
s s= F(l s(trE)I+2m sE)FT Ss=l s(trE)I+2m sE. (9)
J
J denotesthedeterminantofthedeformationgradienttensorF,definedasF=I+(cid:209) us. Similarasinthecaseof
morecomplexbloodflowmodels,alsomorerealisticconstitutiverelationsfortheanisotropicbehaviorofthewalls
ofaneurysmscanbeincludedwhichhoweverisbeyondthescopeofthiscontribution.
2.3 Interactionconditions
Theboundaryconditionsonthefluid-solidinterfaceareassumedtobe
s bn=s sn, vb=vs, on G 0, (10)
t
where n is a unit normalvector to the interface G 0. This implies the no-slip condition for the flow and that the
t
forcesontheinterfaceareinbalance.
3 Discretizationandsolutiontechniques
Inthisstudy,werestrictatthemomenttotwodimensionswhichallowssystematictestsoftheproposedmethods
forbiomedicalapplicationsinaveryefficientwaysuchthatthequalitatitivebehaviourcanbecarefullyanalyzed.
Thecorrespondingfullyimplicit,monolithictreatmentofthefluid-structureinteractionproblemsuggeststhatan
A-stablesecondordertimesteppingschemeandthatthesamefiniteelementsforboththesolidpartandthefluid
regionshould be utilized. Moreover, to circumventthe fluid incompressibility constraints, we have to choose a
stablefiniteelementpair. Forthatreason,theconformingbiquadratic,discontinuouslinearQ P pair,seeFigure
2 1
1forthelocationofthedegreesoffreedom,ischosenwhichwillbeexplainedinthenextsection.
3.1 Spacediscretization
LetusdefinetheusualfinitedimensionalspacesUfordisplacement,V forvelocity,Pforpressureapproximation
asfollows
U ={u∈L¥ (I,[W1,2(W )]2),u=0on¶ W },
V ={v∈L2(I,[W1,2(W )]2)∩L¥ (I,[L2(W )]2),v=0on¶ W },
t t
P={p∈L2(I,L2(W ))},
285
y
v ,u
h h
x
p ,¶ ph,¶ ph
h ¶ x ¶ y
Figure1:LocationofthedegreesoffreedomfortheQ P element.
2 1
thenthevariationalformulationofthefluid-structureinteractionproblemistofind(u,v,p)∈U×V×Psuchthat
theequationsaresatisfiedforall(z ,x ,g )∈U×V×Pincludingappropriateinitialconditions.ThespacesU,V,P
onaninterval[tn,tn+1]wouldbeapproximatedinthecaseoftheQ ,P pairas
2 1
U ={u ∈[C(W )]2,u | ∈[Q (T)]2 ∀T ∈T ,u =0on¶ W },
h h h h T 2 h h
V ={v ∈[C(W )]2,v | ∈[Q (T)]2 ∀T ∈T ,v =0on¶ W },
h h h h T 2 h h
P ={p ∈L2(W ),p | ∈P(T) ∀T ∈T }.
h h h h T 1 h
Letusdenotebyun theapproximationofu(tn),vn theapproximationofv(tn)and pn theapproximationof p(tn).
h h h
ConsiderforeachT ∈T thebilineartransformationy :Tˆ →T totheunitsquareT. Then,Q (T)isdefinedas
h T 2
Q (T)= q◦y −1:q∈span<1,x,y,xy,x2,y2,x2y,y2x,x2y2> (11)
2 T
(cid:8) (cid:9)
withninelocaldegreesoffreedomlocatedatthevertices,midpointsoftheedgesandinthecenterofthequadri-
lateral. ThespaceP(T)consistsoflinearfunctionsdefinedby
1
P (T)= q◦y −1:q∈span<1,x,y> (12)
1 T
(cid:8) (cid:9)
with the function value and both partial derivatives located in the center of the quadrilateral, as its three local
degrees of freedom, which leads to a discontinuous pressure. The inf-sup condition is satisfied (see (Boffi and
Gastaldi,2002));however,thecombinationofthebilineartransformationy withalinearfunctiononthereference
squareP (T)wouldimplythatthebasisonthereferencesquaredidnotcontainthefullbasis. So,themethodcan
1
atmostbefirstorderaccurateongeneralmeshes(see(Arnoldetal.,2002),(BoffiandGastaldi,2002))
b
kp−p k=O(h). (13)
h
The standard remedy is to consider a local coordinate system (x ,h ) obtained by joining the midpoints of the
opposingfacesofT (see(Arnoldetal.,2002),(RannacherandTurek,1992),(Turek,1999)).Then,wesetoneach
elementT
P(T):=span<1,x ,h >. (14)
1
Forthiscase,theinf-supconditionisalsosatisfiedandthesecondorderapproximationisrecoveredforthepressure
aswellasforthevelocitygradient(see(BoffiandGastaldi,2002),(Gresho,1990))
kp−p k=O(h2) and k(cid:209) (u−u )k =O(h2). (15)
h h 0
Forasmoothsolution,theapproximationerrorforthevelocityintheL -normisoforderO(h3)whichcaneasily
2
bedemonstratedforprescribedpolynomialsorforsmoothdataonappropriatedomains.
3.2 Timediscretization
Inviewofamorecompactpresentation,theappliedtimediscretizationapproachisdescribedonlyforthefluidpart
(see(Razzaq,2009)formoredetails). Inthefollowing,werestricttothe(standard)incompressibleNavier-Stokes
equations
v −n D v+v·(cid:209) v+(cid:209) p=f, divv=0, in W ×(0,T], (16)
t
286
forgivenforcefandviscosityn ,withprescribedboundaryvaluesontheboundary¶ W andaninitialconditionat
t=0. Then,theusualq -schemefortimediscretizationreads:
Basicq -scheme: GivenvnandK=t −t ,thensolveforv=vn+1and p=pn+1
n+1 n
v−vn
+q [−n D v+v·(cid:209) v]+(cid:209) p=gn+1, divv=0, in W (17)
K
withrighthandsidegn+1:=q fn+1+(1−q )fn−(1−q )[−n D vn+vn·(cid:209) vn].
The parameter q has to be chosen dependingon the time-stepping scheme, e.g., q =1 for the Backward Euler
(BE), or q =1/2 for the Crank-Nicholson-scheme(CN) whichwe prefer. The pressure term (cid:209) p=(cid:209) pn+1 may
bereplacedbyq (cid:209) pn+1+(1−q )(cid:209) pn,butwithappropriatepostprocessing,bothstrategiesleadtosolutionsofthe
sameaccuracy.Inallcases,weendupwiththetaskofsolving,ateachtimestep,anonlinearsaddlepointproblem
ofgiventypewhichhasthentobediscretizedinspaceasdescribedabove.
Thesetwomethods,CN andBE,belongtothegroupofOne-Step-q -schemes. TheCNschemecanoccasionally
sufferfromnumericalinstabilities becauseof its onlyweak dampingproperty(notstronglyA-stable), while the
BE-schemeisoffirstorderaccuracyonly(however: itisagoodcandidateforsteady-statesimulations). Another
methodwhichhasproventohavethepotentialtoexcelinthiscompetitionistheFractional-Step-q -scheme(FS).
It uses three different values for q and for the time step K at each time level. In (Razzaq et al., 2008), (Turek
etal.,2006)weadditionallydescribedamodifiedFractional-Step-q -schemewhichparticularlyforfluid-structure
interactionproblemsseemstobeadvantageous.Adetaileddescriptionwillappearinthethesis(Razzaq,2009).
3.3 Solutionalgorithms
Thesystemofnonlinearalgebraicequationsarisingfromthegoverningequationsdescribedabovereads
Suu Suv 0 u fu
 Svu Svv kB  v = fv  (18)
c BT c BT 0 p f
 u s v f    p 
whichisatypicalsaddlepointproblem,whereSdescribesthediffusiveandconvectivetermsfromthegoverning
equations. Theabovesystemofnonlinearalgebraicequations(18)issolvedusingNewtonmethodasbasicitera-
tion. ThebasicideaoftheNewtoniterationistofindarootofafunction,R(X)=0,usingtheavailableknown
functionvalueanditsfirstderivative,whereX=(u ,v ,p )∈U ×V ×P. OnestepoftheNewtoniterationcan
h h h h h h
bewrittenas
¶ R −1
Xn+1=Xn− (Xn) R(Xn). (19)
(cid:20)¶ X (cid:21)
1. LetXnbesomestartingguess.
2. SettheresiduumvectorRn=R(Xn)andthetangentmatrixA= ¶ R(Xn).
¶ X
3. Solveforthecorrectiond X
Ad X=Rn.
4. Findoptimalsteplengthw .
5. UpdatethesolutionXn+1=Xn−wd X.
Figure2: OnestepoftheNewtonmethodwithlinesearch.
This basic iteration can exhibit quadratic convergenceprovided that the initial guess is sufficiently close to the
solution. To ensure the convergenceglobally, some improvementsof this basic iteration are used. The damped
Newton method with line search improves the chance of convergenceby adaptively changing the length of the
correctionvector.ThesolutionupdatestepintheNewtonmethod(19)isreplacedby
Xn+1=Xn−wd X, (20)
287
where the parameterw is determinedsuch thata certain error measure decreases(see (Turek,1999), (Hronand
Turek, 2006a) for more details). The Jacobian matrix ¶ R(Xn) can be computed by finite differences from the
¶ X
residualvectorR(X)
¶ R [R](Xn+a e )−[R](Xn−a e )
(Xn)≈ i j j i j j , (21)
(cid:20)¶ X(cid:21) 2a
ij j
wheree aretheunitbasisvectorsinRnandthecoefficientsa areadaptivelytakenaccordingtothechangeinthe
j j
solutionin the previoustime step. Sincewe knowthe sparsity patternof theJacobianmatrixin advance,which
is given by the used finite element method, this computation can be done in an efficient way so that the linear
solverremainsthedominantpartintermsoftheCPUtime(see(Turek,1999),(TurekandSchmachtel,2002)for
more details). A good candidate, at least in 2D, seems to be a direct solver for sparse systems like UMFPACK
(see (Davis and Duff, 1999)); while this choice provides very robust linear solvers, its memory and CPU time
requirements are too high for larger systems (i.e. more than 20.000 unknowns). Large linear problems can be
solvedbyKrylov-spacemethods(BiCGStab,GMRes,see(Barrettetal.,PA1994))withsuitablepreconditioners.
OnepossibilityistheILUpreconditionerwithspecialtreatmentofthesaddlepointcharacterofoursystem,where
weallowcertainfill-inforthezerodiagonalblocks,see(BramleyandWang,1997).
Asanalternative,wealsoutilizeastandardgeometricmultigridapproachbasedonahierarchyofgridsobtained
by successive regular refinement of a given coarse mesh. The complete multigrid iteration is performed in the
standarddefect-correctionsetupwiththeVorF-typecycle. Whileadirectsparsesolver(DavisandDuff,1999)
isusedforthecoarsegridsolution,on finerlevelsa fixednumber(2or 4)ofiterationsbylocalMPSC schemes
(Vanka-likesmoother)(Turek,1999),(Vanka,1985),(HronandTurek,2006a)isperformed. Suchiterationscan
bewrittenas
uvll++11=uvll−w (cid:229) SSuvuu||WW ii SSuvvv||WW ii kB0|W i−1ddeefflulv.
pl+1 pl elementW icuBTs|W i cvBTf|W i 0  deflp
The inverse of the local systems (39×39) can be done by hardware optimized direct solvers. The full nodal
interpolationisusedastheprolongationoperatorPwithitstransposedoperatorusedastherestrictionR=PT (see
(Hronetal.,2002),(Turek,1999)formoredetails).
4 Problemdescription
Inthefollowing,weconsiderthenumericalsimulationofspecialproblemsencounteredintheareaofcardiovas-
cularhemodynamics,namelyflowinteractionwiththick-walleddeformablematerial,whichcanbecomeauseful
toolfordeeperunderstandingoftheonsetofdiseasesofthehumancirculatorysystem,asforexamplebloodcell
andintimaldamagesinstenosis,aneurysmrupture,evaluationofthenewsurgerytechniquesofheart,arteriesand
veins(see (Appanaboyinaetal., 2008),(Lo¨hneretal., 2008)(Valenciaet al., 2008)andthereincited literature).
Inthiscontribution,prototypicalstudiesareperformedforbrainaneurysm.Theword‘aneurysm’comesfromthe
latinwordaneurysmawhichmeansdilatation.Aneurysmisalocaldilatationinthewallofabloodvessel,usually
anartery,duetoadefect,diseaseorinjury. Typically,astheaneurysmenlarges,thearterialwallbecomesthinner
andeventuallyleaksorruptures,causingsubarachnoidhemorrhage(SAH)(bleedingintobrainfluid)orformation
ofabloodclotwithinthebrain.Inthecaseofavesselrupture,thereisahemorrhage,andwhenanarteryruptures,
thenthehemorrhageismorerapidandmoreintense.Inarteriesthewallthicknesscanbeupto30%ofthediameter
anditslocalthickeningcanleadtothecreationofananeurysmsothattheaimofnumericalsimulationsistorelate
the aneurysmstate (unruptureor rupture)with wall pressure, wall deformationand effectivewall stress. Sucha
relationshipwouldprovideinformationforthediagnosisandtreatmentofunruptureandruptureofananeurysm
byelucidatingtheriskofbleedingorrebleeding,respectively.
Inordertousetheproposednumericalmethodsforaneurysmhemodynamics,simplifiedtwo-dimensionalexam-
ples,whichhoweverincludetheinteractionoftheflowwiththedeformablematerial,areconsidered.Flowthrough
adeformableveinwithelasticwallsofabrainaneurysmissimulatedtoanalysequalitativelythedescribedmeth-
ods; here, the flow is drivenby prescribingthe flow velocityatthe inflow partof the boundarywhile the elastic
partoftheboundaryiseitherfixedorstress-free. Bothendsofthewallsarefixedattheinflowandoutflow,and
theflowisdrivenbyaperiodicalchangeoftheinflowattheleftend.
288
4.1 Geometryoftheproblem
Forconvenience,thegeometryofthefluiddomainunderconsiderationiscurrentlybasedon2Dmodels(seeFig.
3) which allows us to concentrate on the detailed qualitative evaluation of our approachbased on the described
monolithicALEformulation. Theunderlyingconstructionofthe(2D)shapeoftheaneurysmcanbeexplainedas
follows:
• Thebentbloodvesselisapproximatedbyquartercirclesaroundtheorigin.
• Theinnermostcirclehastheradius6mm,thenexthas8mm,andthelastonehas8.25mm.
• Thisresultsinonerigidinnerwallandanelasticwallbetween8mmand8.25mmofthickness0.25mm.
Figure 3: Left: Schematic drawing of the measurement section. Middle: Mesh without stents (776 elements).
Right:Meshwithstents(1431elements)whicharepartofthesimulations.
The aneurysmshape is approximatedby two arcs and lines intersecting the arcs tangentially. The midpointsof
thearcsarethesame(-6.75;6),theyhavetheradius1.125mmand1.25mm. Theyareintersectedtangentiallyby
linesatangularvalue1.3radians. Thisresultsinawallthicknessof0.125mmfortheelasticaneurysmwalls(see
Fig.3).Theexaminedstentsareofcircularshape,placedontheneckoftheaneurysm,andweusethree,resp.,five
stents (simplified ‘circles’ in 2D as cutplanesfrom 3D configurations)of differentsize and position. The stents
alsoconsistofagrid, immersedin thebloodflow, whichislocatedatthe inletofthe aneurysmso thatin future
elasticdeformationsofthestentscanbeincluded,too,sinceinreallife,thestentisamedicaldevicewhichconsists
ofawiremetalwiretube. Stentsaretypicallyusedtokeeparteriesopenandarelocatedonthevesselwallwhile
thisstentisimmersedinthebloodflow(Fig.3). Thepurposeofthisdeviceistoreducethefluxintoandwithin
theaneurysminordertooccludeitbyaclotorrupture.Theaneurysmisthenintersectedwiththebloodvesseland
allmissingangularvaluesandintersectionpointscanbedetermined.
4.2 Boundaryandinitialconditions
The(steady)velocityprofile,toflowfromtherighttotheleftpartofthechannel,isdefinedasparabolicinflow,
namely
vb(0,y)=U¯(y−6)(y−8). (22)
Correspondingly,thepulsatileinflowprofileforthenonsteadytestsforwhichpeaksystoleanddiastoleoccurfor
D t=0.25sandD t=0.75srespectively,isprescribedas
vb(t,0,y)=vb(0,y)(1+0.75sin(2p t)). (23)
The natural outflow condition at the lower left part effectively prescribes some reference value for the pressure
variable p, here p =0. While this value could be arbitrarily set in the incompressible case, in the case of a
compressible structure this might have influence onto the stress and consequently the deformation of the solid.
Theno-slipconditionisprescribedforthefluidontheotherboundaryparts,i.e.topandbottomwall, stentsand
fluid-structureinterface.
289
5 Numericalresults
The newtonian fluid used in the tests has a density r b =1.035×10−6kg/mm3 and a kinematic viscosity n b =
3.38mm2/swhichissimilartothepropertiesofblood.IfweprescribetheinflowspeedU¯ =−50mm/s,thisresults
ina ReynoldsnumberRe≈120basedon the prescribedpeaksystole inflowvelocityandthewidth ofthe veins
which is 2mm such that the resulting flow is within the laminar region. Parameter values for the elastic vein in
thedescribedmodelareasfollows: Thedensityoftheupperelasticwallisr s =1.12×10−6kg/mm3,solidshear
modulus is m s =42.85kg/mms2, Poisson ratio is n p =0.4, Young modulus is E =120kN/mm2. As described
before,theconstitutiverelationsusedforthematerialsaretheincompressibleNewtonianmodel(2)forthefluid
and a hyperelastic neo-Hookean material for the solid. This choice includes most of the typical difficulties the
numericalmethodhastodealwith,namelytheincompressibilityandsignificantdeformations.
Fromamedicalpointofview,theuseofstentsprovidesanefficienttreatmentformanagingthedifficultentityof
intracranialaneurysms. Here, thethicknessoftheaneurysmwallisattenuatedandtheaneurysmhemodynamics
changessignificantly.Sincethepurposeofthisdeviceistocontrolthefluxwithintheaneurysminordertoocclude
itbyaclotorrupture,theresultingflowbehaviorintoandwithintheaneurysmisthemainobjective,particularly
inviewofthedifferentstentgeometries. Therefore,wedecidedforthe2Dstudiestolocatethe(2Dpartsofthe)
stentsonlyindirectconnectiontotheaneurysm.
ComparingourstudieswiththeCFDliterature(see(Fernandezetal.,2008),(Appanaboyinaetal.,2008),(Valencia
etal.,2008),(Torrietal.,2007a),(Torrietal.,2007b)), severalresearchgroupsfocusonCFD simulationswith
realistic3Dgeometries,buttypicallyassumingrigidwalls. Incontrast,weconcentrateonthecomplexinteraction
between elastic deformations and flow perturbations induced by the stents. At the moment, we are only able
to perform these simulations in 2D, however, with these studies we should be able to analyse qualitatively the
influence of geometrical details onto the elastic material behavior, particularly in view of more complex blood
modelsandconstitutiveequationsforthestructure.Therefore,theaimsofourstudiescanbedescribedasfollows:
1. Whatistheinfluenceoftheelasticityofthewallsontotheflowbehaviorinsideoftheaneurysm,particularly
w.r.t.theresultingshapeoftheaneurysm?
2. Whatistheinfluenceofthegeometricaldetailsofthe(2D)stents,thatmeansshape,size,position,ontothe
flowbehaviorintoandinsideoftheaneurysm?
3. Do both aspects, small-scale geometrical details as well as elastic fluid-structure interaction, have to be
consideredsimultaneouslyorisoneofthemnegligibleinfirstorderapproximation?
4. AremodernnumericalmethodsandcorrespondingCFDsimulationstoolsabletosimulatequalitativelythe
multiphysicsbehaviorofsuchbiomedicalconfigurations?
Inthefollowing,weshowsomecorrespondingresultsforthedescribedprototypicalaneurysmgeometry,firstfor
thesteadystateinflowprofile,followedbynonsteadytestsforthepulsatileinflow,bothwithrigidandelasticwalls,
respectively.
5.1 Steadyconfigurations
Duetothegiveninflowprofile,whichisnottime-dependent,andduetothelowRenumbers,theflowbehaviour
leadstoasteadystatewhichonlydependsontheelasticityandtheshapeofthestents. Moreover,forthefollowing
simulations, we only treat the aneurysm wall as elastic structure. Then, the aneurysm undergoes some slight
deformationswhichcanhardlybeseeninthefollowingfigures. Howevertheyresultinadifferentvolumeofthe
flowdomain(seeFig.6)andleadtoasignificantlydifferentlocalflowbehavioursincethespacingbetweenstents
andelasticwallsmaychange(seethesubsequentcolorpictures).
290
Figure 4: Deformed mesh for steady configuration without stents, with elastic wall (left). Mesh for rigid wall
(right).
Figure5:Deformedmeshforsteadyconfigurationwithstents: 3stents(left)and5stents(right).
26.65
no stents, elastic fundus
no stents, rigid walls
3 stents, elastic fundus
5 stents, elastic fundus
26.6
26.55
volume 26.5
26.45
26.4
26.35
1 2 3 4 5 6 7 8
time step
Figure6:Resultingvolumeofthefluiddomainfordifferentconfigurations.
Inthefollowingpictures,wevisualizethedifferentflowbehaviourbycoloringduethevelocitymagnitudeandby
showingcorrespondingvectorplotsinsideoftheaneurysm.Particularlytheinfluenceofthenumberofstentsonto
thecompletefluidflowthroughthechannelincludingtheaneurysmcanbeclearlyseen.
291
Figure7:Rigidwallwithoutstents.
Figure8:Elasticaneurysmwallwithoutstents.
Figure9:Elasticaneurysmwallwith3stents.
Figure10:Elasticaneurysmwallwith5stents.
292