Effects of Reynolds and Womersley Numbers on the Hemodynamics of Intracranial Aneurysms

17 Pages·2016·6.48 MB·English

Hindawi Publishing Corporation
Computational and Mathematical Methods in Medicine
Volume 2016, Article ID 7412926, 16 pages
http://dx.doi.org/10.1155/2016/7412926
Research Article
Effects of Reynolds and Womersley Numbers on the
Hemodynamics of Intracranial Aneurysms
HafezAsgharzadehandImanBorazjani
DepartmentofMechanicalandAerospaceEngineering,UniversityatBuffalo,TheStateUniversityofNewYork,Buffalo,NY,USA
CorrespondenceshouldbeaddressedtoImanBorazjani;[email protected]
Received3May2016;Revised29July2016;Accepted10August2016
AcademicEditor:GiancarloFerrigno
Copyright©2016H.AsgharzadehandI.Borazjani. This is an open access article distributed under the Creative Commons
AttributionLicense,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkis
properlycited.
TheeffectsofReynoldsandWomersleynumbersonthehemodynamicsoftwosimplifiedintracranialaneurysms(IAs),thatis,
sidewallandbifurcationIAs,andapatient-specificIAareinvestigatedusingcomputationalfluiddynamics.Forthispurpose,we
carriedoutthreenumericalexperimentsforeachIAwithvariousReynolds(Re = 145.45to378.79)andWomersley(Wo = 7.4
to9.96)numbers.Althoughthedominantflowfeature,whichisthevortexringformation,issimilarforalltestcaseshere,the
propagationofthevortexringiscontrolledbybothReandWoinbothsimplifiedIAs(bifurcationandsidewall)andthepatient-
specificIA.ThelocationofthevortexringinalltestedIAsisshowntobeproportionaltoRe/Wo2 whichisinagreementwith
empiricalformulationsforthelocationofavortexringinatank.InsidewallIAs,theoscillatoryshearindexisshowntoincrease
withWoand1/Rebecausethevortexreachedthedistalwalllaterinthecycle(higherresidenttime).However,thistrendwas
notobservedinthebifurcationIAbecausethestressesweredominatedbyparticletrappingstructures,whichwereabsentatlow
Re=151.51incontrasttohigherRe=378.79.
1.Introduction Therefore, the investigation of their effects on the hemody-
namics of IAs is strongly required. However, contradictory
The rupture of intracranial aneurysms (IAs) is highly asso- conclusionshavebeenmadeintheliteratureonthistopic.In
ciatedwithmortalityandmorbidity[1].Hemodynamicshas fact,JiangandStrother[16]concludedthatincreaseofWom-
a significant role in the growth and rupture of IAs [2– ersley number can significantly increase the complexity of
4]. Among the hemodynamic factors, vortical structures theflowpatternandvortexstructuresintwopatient-specific
determinethecomplexityandstabilityoftheflowpatternin intracranial aneurysms based on their CFD simulations. In
an IA dome, which plays an important role in the rupture contrast, Le et al. [17] stated that Womersley number does
ofIAs[3,5–7].Computationalfluiddynamics(CFD)holds notaffecttheflowfeatureandstructuresbasedontheirCFD
animportantpositionintheinvestigationofhemodynamic simulationsofanIAfromarabbit.Furthermore,Gopalakr-
factorsinaneurysmsbecauseofitshigherresolutionnearthe ishnanetal.[18]statedthatwhileWomersleynumberdoes
wallsrelativetoexperimentalmethodssuchaslaserDoppler not change the vortex mode, but high Womersley number
velocimetry,particleimagevelocimetry,andmagneticreso- is associated with weak vortex rings in their simulation on
nanceimaging,whichisrequiredtocomputehemodynamic abdominal aortic aneurysms. Bouillot et al. [10] concluded
factors such as shear stress correctly [8, 9]. Many investi- that Re has a negligible effect on the flow structures of an
gations have been carried out on the hemodynamics of IAs idealizedsidewallintracranialaneurysmforasteadyinflow
usingexperimentalmethods[10–12]andCFD[7,13,14]. accordingtotheirPIVmeasurements.Asimilarconclusion
Reynolds and Womersley numbers (explained in Sec- was made by Le et al. [17] and Cebral et al. [19] based on
tion3)aretheonlytwonondimensionalparametersrequired theirCFDsimulationsofcerebralaneurysmsinrabbitsand
for full dynamic similarity in pulsatile internal flows [15]. humans,respectively.Incontrast,Gopalakrishnanetal.[18]
2 ComputationalandMathematicalMethodsinMedicine
statedthat,byincreasingRe,thestrengthofthemainvortex of the overset-CURVIB method can be found in [20]. The
structureincreases. methodhasbeenvalidatedagainstexperimentalandbench-
Previous works [7, 19] have shown the significant effect marksolutions[20,24]andhasbeenappliedtoavarietyof
of aneurysm geometry on the aneurysm hemodynamics, problems such as cardiovascular flows [17, 25–28], aquatic
whichisnotinvestigatedfurtherhere.Ouraimistocompare swimming[23,29],andrheology[30].Furthermore,wehave
the effects of the Reynolds and Womersley numbers on validatedourmethodforflowsinsideanimmersedbodyby
thehemodynamicsofbothsidewallandbifurcationsIAsby comparingourresultswiththemeasurementsofthepulsatile
flowthrougha90∘bendinAppendix.AsshowninAppendix,
keepingotherparameters,forexample,inletflowwaveform
the computational results are in excellent agreement with
and geometry, constant. Furthermore, their effect on the
ourprevioussimulationsusingbody-fittedgrids[20]andthe
hemodynamicsisformulatedusingdimensionalanalysisand
comparedwiththeliterature. experimentalresults[31].
Wehaveassumedrigidwallssimilartoprevioussimula-
tions [2, 14, 17, 19] because the displacement of aneurysm’s
2.GoverningEquationsand
wallistypicallysmallandtheflowpatternsofsmalldisten-
theNumericalMethod sibleandrigidmodelsinthecarotidregionareverysimilar
[32]. In addition, we have assumed Newtonian fluid in our
In this section, Einstein’s tensor notation, where repeated
simulationsbecausethenon-Newtonianeffectsarenegligible
indicesimplysummation,isusedunlessotherwiseindicated in larger (>500𝜇m) arteries [33] and previous simulations
(𝑘, 𝑑, and 𝑞 = 1,2,3). The governing equations are the
3𝐷, unsteady incompressible Navier-Stokes equations for a of both Newtonian and non-Newtonian fluids have shown
similarflowpatterns[19].
Newtonianfluidincurvilinearcoordinatesasfollows[20]:
𝜕 𝑈𝑞
( )=0, 3.DescriptionofSimulatedTestCases
𝜕𝜉𝑞 𝐽
Numerical simulations have been carried out on two sim-
1 𝜕 𝜕 1 𝜕 1
(𝑈𝑚)+(𝜉𝑚){ (𝑈𝑘𝑢 )+ (𝑃𝜉𝑘 ) plified geometries (sidewall and bifurcation), which can
𝐽𝜕𝑡 𝑥𝑞 𝜕𝜉𝑘 𝐽 𝑞 𝜕𝜉𝑘 𝐽 𝑥𝑞 (1) be considered as simplified models of IA, and a patient-
specificIAgeometry.Figure1showstheconfigurationofthe
− 1 𝜕2 (1𝑔𝑘𝑑𝑢 )}=0, simplifiedmodelsfor(a)sidewalland(b)bifurcationIAs.The
Re𝜕𝜉𝑘𝜕𝜉𝑑 𝐽 𝑞 inletandoutletofthegeometriesareconstructedfromapipe
with the diameter 𝐷. The aneurysm dome is modeled by a
where 𝑈𝑞, 𝜉𝑚, and 𝑥𝑞 are the contravariant velocity, curvi- cone-likeshapewithellipticalbaseandlocus,whoseradiuses
linear coordinate, and Cartesian coordinates components, are𝑎=𝐷/2,𝑏=𝐷,and𝑐=0.6𝐷.
respectively. 𝑃 and 𝑡 are the nondimensional pressure and Figure2(a)showstheimmersedbodyandoversetgrids
time,respectively.𝐽istheJacobianofthegeometrictransfor- layout for sidewall type IAs and Figure 2(b) shows the
mation,𝐽 = 𝜕(𝜉1,𝜉2,𝜉3)/𝜕(𝑥1,𝑥2,𝑥3),and𝜉𝑥𝑘 = 𝜕𝜉𝑘/𝜕𝑥𝑞 are immersedbodyandaCartesiangridlayoutforthepatient-
𝑞 specific IA geometry. The inlet and outlet(s) in simplified
themetricsofthetransformation.𝑢𝑞 isthenondimensional models are meshed by body-fitted curvilinear grids; the
Cartesian velocity, 𝑔𝑘𝑑 is the contravariant metric tensor, circularbaseismeshedby41×41gridnodesandtheaxial
𝑔𝑘𝑑 =𝜉𝑥𝑘 𝜉𝑥𝑑 ,andReistheReynoldsnumberoftheflowbased (flow) direction is meshed by 81 grid nodes. The geometry
𝑞 𝑞
oncharacteristiclengthandvelocityscales. of the aneurysm is placed as an immersed boundary onto
We use curvilinear/immersedboundary (CURVIB) and theuniformgrid,andallsubdomains(Figure2(a))aresolved
overset grid methods, which are extensively described simultaneously[20].Thesizeofthedomainthatcontainsthe
and validated [20]. A sharp-interface immersed boundary immersed body is 𝐿𝑥 = 1.4𝐷, 𝐿𝑦 = 2𝐷, and 𝐿𝑧 = 3𝐷
methodisusedtohandlethe3𝐷,arbitrarycomplexbound- and 𝐿𝑥 = 1.4𝐷, 𝐿𝑦 = 3𝐷, and 𝐿𝑧 = 3𝐷 for sidewall and
aries(IAgeometryinthisstudy)insidetheflowdomain[21]. bifurcation simplified IAs, respectively. The uniform grid is
Thenodesthatareoutsidetheflowdomainareblankedout meshedby73×109×153and73×153×153nodesforsidewall
and do not affect the solution. These nodes are identified andbifurcationsimplifiedIAs, respectively.Suchresolution
using an efficient ray-tracing algorithm [22]. The boundary wasfoundtobeenoughtoobtaingridindependentsolutions
conditions are reconstructed on the fluid nodes in the basedonourextensivegridrefinementstudies.Furthermore,
immediate vicinity of the immersed boundary along the the results obtained using a similar resolution for the pul-
normal direction to the boundary [21]. The method has satile flow through a 90∘ bend show good agreement with
been shown to be second-order accurate for a variety of experimentalresults(Appendix).Thesizeofthedomainthat
flows [21, 23]. The overset grid approach is implemented to contains the immersed body of a patient-specific geometry
reduce wasted nodes in a domain, which are blanked out is 𝐿𝑥 = 4.64𝐷, 𝐿𝑦 = 4.04𝐷, and 𝐿𝑧 = 5.8𝐷, where 𝐷
by the immersed boundary method [20]. In this approach, (= 0.0048m) is the hydraulic diameter of the cross section
a complex flow domain is divided into several arbitrary oftheinletvesselandCartesiangridlayout.Theuniformgrid
subgrids with overlaps. To solve the governing equations forthepatient-specificgeometryismeshedby233×201×289
at each subgrid, boundary conditions at the interfaces are nodes,whichhasasimilargridresolutionwiththesimplified
constructed by interpolation from host subgrid. The details IAs.
ComputationalandMathematicalMethodsinMedicine 3
x
a
z 2 x 2 z
( ( +( ( =1 c z 2 y 2
a c ( ( +( ( =1
a b
D
b
5D y 5D
z
x
a z Outlet D Outlet
(az(2+(cx(2=1 c (az(2+(by(2=1 Δs
D
b D
5D y 5D
z
Inlet Δs Outlet
D
Inlet
(a) (b)
Figure1:SchematicillustratingofIAsimplifiedmodels(a)sidewalland(b)bifurcationIA,whereellipseradiusesare𝑎=𝐷/2,𝑏=𝐷,and
𝑐=0.6𝐷).Δ𝑠denotesthelocationofthevortexringfromtheproximalwallandinlet/outletcrosssectioninthesidewallandbifurcationIA,
respectively.
Lz
Lx
Outlet Plane
L Immersed body
z
Body-fitted
curvilinear mesh Ly Immersed body
L
y
Outlet
y L Body-fitted
x
z curvilinear mesh
Cartesian mesh
x Cartesian mesh Inlet
(a) (b)
Figure2:Immersedbodyandoversetgridslayoutfor(a)sidewallIAsimplifiedmodeland(b)patient-specificIAgeometry.Theimmersed
bodyismeshedwithtriangularelements.
TheNeumannboundaryconditionisappliedattheoutlet the blood density and dynamic viscosity, resp. [7]). The
boundaries, and the volumetric flow rate of each outlet is Womersley number is defined as Wo = 𝐷√2𝜋/𝑇𝜗, where
calculatedbasedontheprincipleofoptimalwork[36].The 𝑇 is the period of the waveform calculated based on the
inlet boundary condition is a uniform (plug) flow in space heart rate. The simulations on the sidewall and bifurcation
but pulsatilein timewith a waveform as shown in Figure3 IAs are carried out with two Womersley and two Reynolds
withpulsatilityindexPI = Δ𝑈/𝑈 = 1.43(whereΔ𝑈isthe numbers. The Womersley number is modified by changing
differencebetweenthepeaksystolicandminimumdiastolic theheartrateandReynoldsnumberismodifiedbyaltering
velocities).Thepulsatilevelocitywaveformisfromapatient the inlet bulk velocity in the physiologically relevant range
with a cerebral aneurysm [34]. The numerical simulations [37]. The nondimensional time-step (Δ𝑡) is calculated by
were carried out for three cycles to obtain quasi-steady dividingthenondimensional𝑇into3200equaltimeinstants
results,whicharenotaffectedbytheinitialcondition. forallsimulatedcases.Table1showsthesimulatedcases,their
The Reynolds number is defined as Re = 𝑈𝐷/𝜗, where specifications,anddenotednames.Itshouldbenotedthatthe
𝑈 is the average bulk inlet velocity and 𝐷 = 0.005m and calculatedReandWoinTable1arewithinthephysiological
0.0048maretheinletdiameterforsimplifiedIAsandthepa- rangeinIAs(173<Re<914and5<Wo<30[7,14,38]).
tient-specific IA, respectively, and 𝜗 = 3.3 × 10−6m2s−1 An (aneurysm number) in Table 1 is the ratio of the
is the blood kinematic viscosity (calculated based on 𝜌 = transport to vortex formation time scales in IAs [17]. The
1056kgm−3 and 𝜇 = 0.0035kg(ms)−1, where 𝜌 and 𝜇 are transporttimescaleinsidewallIAsisthetimeittakesforthe
4 ComputationalandMathematicalMethodsinMedicine
Table1:Specificationsofsimulations.H,high;L,low;𝐷,diameter;𝑇,heartbeatperiod;Re,Reynoldsnumber;Wo,Womersleynumber;𝑈,
theaveragebulkinletvelocity;An,aneurysmnumber;Δ𝑡,thenondimensionaltime-step.
Simulation IAtype 𝐷(m) Heartbeatsper 𝑇(s) Wo 𝑈(m/s) Re An Δ𝑡
minute
ReLWoL Simplifiedsidewall 0.005 75 0.8 7.71 0.1 151.51 1.43 0.002
ReHWoL Simplifiedsidewall 0.005 75 0.8 7.71 0.25 378.79 1.43 0.005
ReHWoH Simplifiedsidewall 0.005 125 0.48 9.96 0.25 378.79 1.43 0.003
ReLWoL Simplifiedbifurcation 0.005 75 0.8 7.71 0.1 151.51 2.86 0.002
ReHWoL Simplifiedbifurcation 0.005 75 0.8 7.71 0.25 378.79 2.86 0.005
ReHWoH Simplifiedbifurcation 0.005 125 0.48 9.96 0.25 378.79 2.86 0.003
ReLWoL Patient-specific 0.0048 75 0.8 7.4 0.1 145.45 2.4 0.002
ReHWoL Patient-specific 0.0048 75 0.8 7.4 0.25 363.64 2.4 0.005
ReHWoH Patient-specific 0.0048 125 0.48 9.56 0.25 363.64 2.4 0.003
2 investigated. Because of the substantial difference between
vortical structures in sidewall and bifurcation IAs, they are
investigatedindifferentsections,thatis,Sections4.1and4.2
forsidewallandbifurcationIAs,respectively.InSection4.3,
1.5
whetherobservedeffectsofReynoldsandWomersleynum-
bers (on flow patterns and vorticalstructures) in simplified
IAscanbeseeninapatient-specific(complex)IAisexplored.
/U 1 In Section 4.4, the location of the vortex ring in simplified
u
IAs and the patient-specific IA is investigated using the
nondimensionalizationmethod.Finally,theeffectofReand
Woontheoscillatoryshearindexandwallshearstressofall
0.5 testcasesisinvestigatedinSection4.5.Thethree-dimensional
vortical structures are visualized by the isosurfaces of 𝑞-
criteria[35].
0
0 0.2 0.4 0.6 0.8 1 4.1.SidewallIASimplifiedModels. Inthissection,theeffectof
t/T WomersleyandReynoldsnumbersontheflowpatternsand
vorticalstructures of simplified sidewall IAs is investigated.
Figure3:Theinletvelocitywaveformduringoneheartbeat.Adapt-
Figure4showsthetimeevolutionofthenondimensionalout-
edfrom[34].
of-planevorticityandin-planevelocityvectorsforsimplified
sidewall IAs with different Re and Wo. In Figure 4, the
parentarteryflowtotransportafluidparticleacrosstheIA dominantflowpatternsaresimilarfordifferentReandWo.In
neckandthetransporttimescaleinbifurcationIAsisthetime allcasesavortexring,denotedasR1,startstoformfromthe
ittakesfortheinletarteryflowtotransportafluidparticleto proximal wall in the early acceleration phase (Figure 4(a)).
theoutlets.Ancanbedefinedasfollows: Beforetheinletvelocityreachesthepeaksystole,thevortex
𝑊 ring is convected across the IA neck (Figures 4(b)–4(d)).
An=𝛼 𝐷modPI, (2) Eventually it hits the distal wall and diffuses just after the
peak systole and in the deceleration phase (Figures 4(e)
where 𝑊mod is the aneurysm neck and the outlet diameter and 4(f)). A similar behavior can be observed in Figure 5,
in sidewall and bifurcation IAs, respectively, defined based which shows the time evolution of the three-dimensional
on the transport time scale. 𝛼 = 1 and 2 is determined vortical structures visualized by the isosurface of 𝑞-criteria
basedonthetransporttimescaledefinitionforsidewalland for simplified sidewall IAs with different Re and Wo. In
bifurcation IAs, respectively. If vortex formation time scale all cases the vortex ring forms at early acceleration phase
is smaller, that is, An > 1, the vortex is formed before it is (Figure 5(a)) and is convected across the IA neck (Figures
advectedoutofthedomearea.Incontrast,ifAn<1theflow 5(b)–5(d)) and hits the distal wall before the peak systole
is dominated by a stationary shear layer from the proximal phase.Subsequently,itrollsuptoformarecirculationvortical
walltothedistalwallthroughtheIAneck. structure, denoted as VS1, which stays until the end of the
cycle. The size of vortex ring for all cases is approximately
4.ResultsandDiscussion similar(Figure5).TheAnforallthesimplifiedsidewallcases
isequalto1.43,whichisconsistentwiththepreviouswork
Inthissection,theeffectsofReynoldsandWomersleynum- [11,17,26]statingthatwhenAn>1theflowisdominatedby
bers on the flow patterns and vortical structures of IAs are thevortexformation.
ComputationalandMathematicalMethodsinMedicine 5
Nondimensional out-of-plane vorticity
−8 −4 0 4 8
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
u/ 0 u/ 0 u/ 0 u/ 0 u/ 0 u/ 0
0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
t/T R1 t/T R1 t/T R1 t/T R1 t/T t/T
ReHWoH
R1 R1 R1 R1
ReHWoL
R1 R1 R1 R1
ReLWoL
y
z
(a) (b) (c) (d) (e) (f)
Figure4:Evolutionofthenondimensionalout-of-planevorticityandin-planevelocityvectorsforsimplifiedsidewallIAswithdifferentRe
andWo;(upperrow)Re=378.79andWo=9.96;(middlerow)Re=378.79andWo=7.71;and(lowerrow)Re=151.51andWo=7.71at
varioustimeinstantsduringthecycle:(a)𝑡/𝑇 = 0.21,(b)𝑡/𝑇 = 0.24,(c)𝑡/𝑇 = 0.27,(d)𝑡/𝑇 = 0.29,(e)𝑡/𝑇 = 0.44,and(f)𝑡/𝑇 = 0.75.R1
depictsthevortexringinsimplifiedsidewallIAandthegreencrossshowsthelocationofthemaximum𝑞-criteria[35]onthemidplaneof
IAs.
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
u/ 0 u/ 0 u/ 0 u/ 0 u/ 0 u/ 0
0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
VS1
t/T VS1 t/T R1 t/T R1 t/T R1 t/T VS1 t/T
H
o
W
H
e
R
R1
VS1
R1 R1 R1 VS1 VS1
L
o
W
H
e
R
R1
oL VS1 R1 R1 R1 VS1 VS1
W
L
e
R
y
R1
x z (a) (b) (c) (d) (e) (f)
Figure5:Evolutionandtopologyofthethree-dimensionalvorticalstructuresforsidewallIAswithdifferentReandWo;(upperrow)Re=
378.79andWo=9.96;(middlerow)Re=378.79andWo=7.71;and(lowerrow)Re=151.51andWo=7.71atvarioustimeinstantsduring
thecycle:(a)𝑡/𝑇=0.21,(b)𝑡/𝑇=0.24,(c)𝑡/𝑇=0.27,(d)𝑡/𝑇=0.29,(e)𝑡/𝑇=0.44,and(f)𝑡/𝑇=0.75.
6 ComputationalandMathematicalMethodsinMedicine
Nondimensional out-of-plane vorticity
−10 −5 0 5 10
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
u/ 0 u/ 0 u/ 0 u/ 0 u/ 0 u/ 0
0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
t/T t/T t/T t/T t/T t/T
ReHWoH
R2 R2 R2 R2
ReHWoL
R2 R2 R2 R2
ReLWoL
y R2 R2 R2 R2
z (a) (b) (c) (d) (e) (f)
Figure6:Evolutionofthenondimensionalout-of-planevorticityandin-planevelocityvectorsforsimplifiedbifurcationIAswithdifferent
ReandWo;(upperrow)Re=378.79andWo=9.96;(middlerow)Re=378.79andWo=7.71;and(lowerrow)Re=151.51andWo=7.71
atvarioustimeinstantsduringthecycle:(a)𝑡/𝑇 = 0.19,(b)𝑡/𝑇 = 0.22,(c)𝑡/𝑇 = 0.24,(d)𝑡/𝑇 = 0.25,(e)𝑡/𝑇 = 0.44,and(f)𝑡/𝑇 = 0.75.
R2depictsthevortexringinsimplifiedbifurcationIAsandthegreencrossshowsthelocationofthemaximum𝑞-criteriaonthemidplaneof
IAs.
The only difference in the vortex ring formation and anditseffectonhemodynamicstressesisdiscussedindetail
evolutionisthelocationofthevortexringfordifferent Wo inSections4.4and4.5,respectively.
and Re (Figures 4 and 5). To further investigate the effect
of Re and Wo, the location of the center of the vortex ring 4.2. Bifurcation IA Simplified Models. In this section, the
at different time instants in the cycle is identified by the effect of Womersley and Reynolds numbers on the flow
maximum 𝑞-criteria on the midplane of simplified sidewall patternsandvorticalstructuresofsimplifiedbifurcationIAs
IAsandisshownbythegreencrossesinFigure4.Asitcan is investigated. Figure 6 shows the time evolution of the
be observed by comparing different rows of Figures 4(a)– nondimensionalout-of-planevorticityandin-planevelocity
4(d), the location of vortex ring is different at specific time vectors for simplified bifurcation IAs with different Re and
instants. The effect of Wo on the convection of the vortex Wo.FromFigure6,thedominantflowpatterns,thatis,vortex
ringcanbeobservedbycomparingmiddleandupperrows formation,aresimilarinallsimulationsofsimplifiedbifurca-
of Figures 4(a)–4(d) at specific time instants; that is, the tionIAs.Thisisconsistentwithourprevioussimulationswith
convection of the vortex ring corresponding to Wo = 7.71 An>1intermsofvortexformation.Thevortexring,denoted
(middle row) is higher than Wo = 9.96 (upper row) with asR2,startstoformatinlet/outletcrosssectionsintheearly
similarRe=378.79becausethedistanceofvortexringfrom acceleration phase (Figure 6(a)). Subsequently, the vortex
proximalwallishigherforWo=7.71(middlerow)atsimilar ringgrowsinsize,whileitmovestotheoutlets(Figures6(b)–
𝑡/𝑇 in the cycle. The effect of Re on the convection of the 6(d)). Finally in the peak systole, the complicated vortical
vortexringcanbeobservedbycomparingmiddleandlower structuresareobserved(Figure6(e)).
rowsofFigures4(a)–4(d)atspecifictimeinstants;thatis,the At peak systole more segments of vortical structures in
convectionofthevortexringcorrespondingtoRe = 378.79 the dome can be observed for the higher Re (middle and
(middle row) is higher than Re = 151.51 (lower row) with upper rows of Figure 6(e)) in comparison to the lower Re
similar Wo = 7.71. This behavior can be observed by the (thelowerrowofFigure6(e)).Thisistheconsequenceofthe
locationofthenondimensionalout-of-planevorticityaswell. vortexringinteractionswithothervorticalstructures,which
TheconvectionofthevortexringinsimplifiedsidewallIAs cannotbeseenintwo-dimensionalvorticityonthemidplane.
ComputationalandMathematicalMethodsinMedicine 7
Nondimensional vorticity-x
−8 −4 0 4 8
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
u/ 0 u/ 0 u/ 0 u/ 0 u/ 0 u/ 0
0 t0/.T5 1 VS20 t0/.T5 1 VS2 0 t0/.T5 1 VS20 t0/.T5 1 VS2 0 t0/.T5 1 0 t0/.T5 1
ReHWoH
R2 R2 R2 R2
VS2 VS2 VS2 VS2
ReHWoL
R2 R2 R2 R2
ReLWoL
y
R2 R2 R2 R2
z
x (a) (b) (c) (d) (e) (f)
Figure7:Evolutionandtopologyofthethree-dimensionalvorticalstructurescoloredbyvorticity-𝑥magnitudeforsimplifiedbifurcationIAs
withdifferentReandWo;(upperrow)Re=378.79andWo=9.96;(middlerow)Re=378.79andWo=7.71;and(lowerrow)Re=151.51
andWo=7.71atvarioustimeinstantsduringthecycle:(a)𝑡/𝑇=0.19,(b)𝑡/𝑇=0.22,(c)𝑡/𝑇=0.24,(d)𝑡/𝑇=0.25,(e)𝑡/𝑇=0.44,and(f)
𝑡/𝑇=0.75.
In order to show the three-dimensional vortical structures, forRe < 200.Asaconsequence,morecomplexandhighly
Figure 7 plots the time evolution of the three-dimensional asymmetric flow patterns are observed for Re = 387.79
vortical structures visualized by the isosurfaces of the 𝑞- (upper and middle row of Figure 7(e)) in contrast to Re =
criteria for simplified bifurcation IAs with different Re and 151.51becauseofthevortexringinteractionwiththeparticle
Wo.Otherthanthevortexring,anothervorticalstructurecan trappingvorticalstructure.
beseenattheinlet/outletcrosssectionclosetotheIAneck, The green crosses in Figure 6 show the location of the
denotedasVS2,forRe = 387.79(middleandupperrowsin center of the vortex ring by the maximum 𝑞-criteria on
Figures5(a)–5(d)).VS2isobservedbyVigoloetal.[39]for the midplane of simplified bifurcation IAs at different time
thefirsttimeandlaterbyChenetal.[40]in𝑇-junctionpipes instantsinthecycle.TheeffectofWoontheconvectionofthe
and named “particle trapping vortical structures.” Particle vortexringcanbeobservedbycomparingmiddleandupper
trapping vortical structures are distinguished from other rowsofFigures6(b)–6(d)atspecifictimeinstants;thatis,the
structuresbythevorticitymagnitudeinnormaldirectionto convection of the vortex ring corresponding to Wo = 7.71
the midplane of models (𝑥-direction) in Figure 7, because (middlerow)ishigherthanWo=9.96(upperrow)atsimilar
particletrappingstructureshavealmostzerovorticityinthe Re=378.79and𝑡/𝑇.TheeffectofReontheconvectionofthe
normal direction to the midplane of the models. Because vortexringcanbeobservedbycomparingmiddleandlower
of the existence of an IA in the current study, the particle rows of Figures 6(b)–6(d) at specific time instants; that is,
trappingisdistortedincomparisontothatinthe𝑇-junction theconvectionofvortexringcorrespondingtoRe = 378.79
pipe[39]. (middle row) is higher than Re = 151.51 (lower row) with
One effect of Re in simplified bifurcation IAs is that similarWo = 7.71atthesame𝑡/𝑇.Thisbehaviorissimilar
particle trapping vortical structures do not form at the towhatwasobservedinsidewallIAsintheprevioussection
accelerationphaseofRe=151.51(thelowerrowofFigure7). (Figure 4). The convection of the vortex ring in simplified
This behavior is also observed by Vigolo et al. [39] who bifurcation IAs and its effect on hemodynamic stresses is
reportedthatparticletrappingvorticalstructuresdonotform discussedindetailinSections4.4and4.5,respectively.
8 ComputationalandMathematicalMethodsinMedicine
Nondimensional out-of-plane vorticity
−8 −4 0 4 8
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
u/ 0 u/ 0 u/ 0 u/ 0 u/ 0 u/ 0
0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
t/T t/T t/T t/T t/T t/T
R1 R1
ReHWoH
Δs
R1
R1
R1 R1
ReHWoL
R1
R1 R1
ReLWoL
R1
R1
(a) (b) (c) (d) (e) (f)
Figure8:Evolutionofthenondimensionalout-of-planevorticityandin-planevelocityvectorsforpatient-specificIAswithdifferentReand
Wo;(upperrow)Re=363.64andWo=9.56;(middlerow)Re=363.64andWo=7.4;and(lowerrow)Re=145.45andWo=7.4atvarious
timeinstantsduringthecycle:(a)𝑡/𝑇=0.19,(b)𝑡/𝑇=0.22,(c)𝑡/𝑇=0.24,(d)𝑡/𝑇=0.25,(e)𝑡/𝑇=0.44,and(f)𝑡/𝑇=0.75.R1depictsthe
vortexringinpatient-specificIAsandthegreencrossshowsthelocationofthemaximum𝑞-criteriaonthemidplaneofIAs.
4.3.Patient-SpecificIA. Inthissection,theeffectofWomer- by the vortex formation. The only effect of Wo and Re on
sleyandReynoldsnumbersontheflowpatternsandvortical thevortexringformationandevolutionisthelocationofthe
structuresofapatient-specificIA,whichisofsidewalltype, vortexringatspecifictime(Figures8and9).
is investigated. Figure 8 shows the time evolution of the Bycomparingtheflowstructuresinpatient-specificIAs
nondimensionalout-of-planevorticityandin-planevelocity (Figure 9) with simplified sidewall IAs (Figure 5), it can be
vectorsforpatient-specificIAswithdifferentReandWo.The observedthatmorecomplexandbrokenvorticalstructures
planeisselectedthewaythatvorticalstructuresinthedome are formed in the patient-specific IA in comparison to the
canbeobservedclearly(theplaneisshowninFigure2(b)). simplifiedsidewallIA,whilethedominantflowpattern(vor-
From Figure 8, the dominant flow patterns are similar to texringformationandevolution)issimilar.Thevortexring
sidewallIAs(Figure4).InallcaseswithdifferentReandWo inpatient-specificIAisdeformedbecauseofthecomplexity
avortexring,denotedasR1,startstoformfromtheproximal ofitsgeometry.
wallintheearlyaccelerationphase(Figure8(a)).Thevortex TofurtherinvestigatetheeffectofReandWo,thelocation
ring evolves as it is convected across the IA neck (Figures ofthecenterofthevortexringatdifferenttimeinstantsinthe
8(b)–8(d)). Finally it hits the distal wall and diffuses after cycle is identified by the maximum 𝑞-criteria on the plane
thepeaksystoleandinthedecelerationphase(Figures8(e) (shown in Figure2(b)) ofpatient-specificIAs and is shown
and8(f)).AsimilarbehaviorcanbeobservedfromFigure9 by the green crosses in Figure 8. The similar behavior with
fordifferentReandWo,whichshowsthetimeevolutionof simplified models (Figures 4 and 6) can be observed in the
the three-dimensional vortical structures visualized by the patient-specific geometry (Figure 8); that is, the convection
isosurfaceof𝑞-criteriaforpatient-specificIAs.Inallcasesthe ofthevortexringcorrespondingtoWo=7.4(middlerow)is
vortexringformsatearlyaccelerationphase(Figure9(a))and higherthanWo=9.56(upperrow)withsimilarRe=363.64
isconvectedacrosstheIAneck(Figures9(b)–9(d))andhits because the distance of vortex ring from proximal wall is
thedistalwallbeforethepeaksystolephase.Anforallcases higherforWo=7.4(middlerow)atsimilar𝑡/𝑇inthecycle.
is equal to 2.4, which is consistent with the previous work The convection of the vortex ring corresponding to Re =
[11,17,26]statingthatwhenAn > 1theflowisdominated 363.64(middlerow)ishigherthanRe = 145.45(lowerrow)
ComputationalandMathematicalMethodsinMedicine 9
2 2 2 2 2 2
U 1 U 1 U 1 U 1 U 1 U 1
/ / / / / /
u 0 u 0 u 0 u 0 u 0 u 0
0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
t/T t/T t/T t/T t/T t/T
R1 R1 R1 R1
H
o
W
H
e
R
R1
R1 R1
L
o
W
H
e
R
R1 R1 R1
R1
L
o
W
L
e
R
(a) (b) (c) (d) (e) (f)
Figure9:Evolutionandtopologyofthethree-dimensionalvorticalstructuresforpatient-specificIAswithdifferentReandWo;(upperrow)
Re=363.64andWo=9.56;(middlerow)Re=363.64andWo=7.4;and(lowerrow)Re=145.45andWo=7.4atvarioustimeinstants
duringthecycle:(a)𝑡/𝑇=0.19,(b)𝑡/𝑇=0.22,(c)𝑡/𝑇=0.24,(d)𝑡/𝑇=0.25,(e)𝑡/𝑇=0.44,and(f)𝑡/𝑇=0.75.
withsimilarWo = 7.4.Theconvectionofthevortexringin wallsooner.Thisbehaviorcanbeexplainedbydimensional
thepatient-specificIAsanditseffectofhemodynamicstresses analysis.Thedimensionlessdistanceofthevortexringfrom
arediscussedindetailinSections4.4and4.5,respectively. theproximalwallcanbedefinedas
Δ𝑠 𝑢𝑡
4.4.TheEffectofReandWoontheLocationoftheVortexRing. = , (3)
𝐷 𝐷
Inthissection,theeffectofWoandReontheconvectionof
the vortex ring in simplified IAs (sidewall and bifurcation) where 𝑢 and 𝑡 are the vortex convection velocity and time,
and a patient-specific IA is investigated, and a formula for respectively.Δ𝑠isthedistanceofthevortexringlocationfrom
thelocationofthecenterofthevortexringisderivedusing the proximal wall in simplified sidewall IAs (Δ𝑧), also the
dimensional analysis. Figure 10 plots the distance of the proximalwallinpatient-specificIAs(Δ𝑠),andthedistanceof
maximum𝑞-criteria(locationofthecenterofthevortexring) thevortexringlocationfromtheinlet/outletcrosssectionin
onthemidplaneofsimplifiedIAs(sidewallandbifurcation) bifurcationIAs(Δ𝑦)asdenotedinFigure1.Here,wesimply
and the plane of patient-specific IAs at different Re and approximatethevortexringconvectionvelocitybyassuming
Wo. These distances (Δ𝑠) are calculated from the proximal thatitisequaltothebulkflowvelocityintheparentartery.
wall for sidewall (in 𝑧-direction shown in Figure 1(a)) and Applyingthedimensionalanalysisin(3)leadsto
also from the proximal wall for patient-specific IAs (in 𝑠-
Δ𝑠 𝑈𝑇
direction shown in upper row of Figure 8(f)) and from the =𝑢𝑡 , (4)
inlet/outletcross section forbifurcationIAs (in𝑦-direction 𝐷 𝐷
shown in Figure 1(b)). Δ𝑠 is normalized by the diameter of
theinletparentartery(𝐷)inalltestcases.Itcanbeobserved where𝑢=𝑢/𝑈and𝑡=𝑡/𝑇.Re=𝑈𝐷/𝜗andWo2 =𝐷22𝜋/𝑇𝜗
that the vortex ring formation starts in 𝑡/𝑇 ≈ 0.19 and definitionscanbeappliedin(4),whichresultsin
0.17forsimplifiedsidewallandbifurcationIAs,respectively, Δ𝑠 Re
regardless of Re and Wo. The location of the vortex ring 𝐷 =2𝜋𝑢𝑡Wo2. (5)
is different in various Re and Wo at similar instants in
the cycle. Based on Figure 10 and previous discussion in Based on (5), decreasing Re and increasing Wo result in
Figures4,6,and8,highReandlowWoareassociatedwith decreasingdistanceofthevortexringfromtheproximalwall
faster convection of the vortex ring since the vortex ring atspecificinstants(𝑡),whichisinagreementwithFigure10.It
corresponding to high Re and low Wo reaches the distal isnotedthat𝑢(=𝑢/𝑈)isthedimensionlessvelocity,whichis
10 ComputationalandMathematicalMethodsinMedicine
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
D D
/ 0.5 / 0.5
z y
Δ Δ
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0.15 0.2 0.25 0.3 0.35 0.15 0.2 0.25 0.3 0.35 0.4
t/T t/T
ReLWoL ReHWoH ReLWoL ReHWoH
ReHWoL ReHWoL
(a) (b)
1.8
1.6
1.4
1.2
1
D
/
s
Δ 0.8
0.6
0.4
0.2
0
0.15 0.2 0.25 0.3 0.35
t/T
ReLWoL ReHWoH
ReHWoL
(c)
Figure10:Thedistanceofthemaximum𝑞-criteria[35]onthemidplanefordifferentReandWofor(a)sidewallIAs,wherethevortexring
distanceismeasuredfromtheproximalwallinthe𝑧-direction;(b)bifurcationIAs,wherethevortexringdistanceismeasuredfromthe
inlet/outletcrosssectionin𝑦-direction;(c)patient-specificIAs,wherethevortexringdistanceismeasuredfromtheproximalwallinthe
𝑠-directionatvarioustimeinstantsduringapartofcyclethatvortexringformsandbreaksdown.ThisdistanceisnormalizedbytheIAwidth
(𝐷),distancebetweenproximalanddistalwall(𝑊),andoutletdiameter(𝐷)insimplifiedsidewall,patient-specific,andsimplifiedbifurcation
IAs,respectively.ValuesforReHWoH,ReHWoL,andReLWoLcanbefoundinTable1fordifferentgeometries.
similarforallsimulationsatthespecificdimensionlesstime Thetrendobservedin(6)agreeswellwiththatofFigure10,
(𝑡/𝑇)regardlessofReandWo(ascanbeseenfromFigure3). meaning that Δ𝑠(ReHWoL) > Δ𝑠(ReHWoH) and Δ𝑠(ReHWoH) >
TheeffectivenessofRe andWo onthenormalizeddistance Δ𝑠(ReLWoL)atspecificinstants(𝑡/𝑇).Thetrendofthelocation
ofthevortex ringfromtheproximalwallcan be compared of the vortex ring for various Re and Wo agrees well with
using(5)asfollows: the formulation developed from experiments on thin core
Δ𝑠
(Re=378.79,Wo=7.71) =1.67, ringsgeneratedbyapistonguninwater[41].Bysubstituting
Δ𝑠 simulationspecificationsofcurrentstudyontheformulations
(Re=378.79,Wo=9.96)
(6) expressed in [41] Δ𝑠(Re=378.79,Wo=7.71)/Δ𝑠(Re=378.79,Wo=9.96) ≈
Δ𝑠(Re=378.79,Wo=9.96) =1.5. 1.62 and Δ𝑠(Re=378.79,Wo=9.96)/Δ𝑠(Re=151.51,Wo=7.71) ≈ 1.73 can
Δ𝑠(Re=151.51,Wo=7.71) be reached which shows promising agreement with the