Fluid Structural Analysis of Human Cerebral Aneurysm Using Their Own Wall Mechanical Properties

19 Pages·2015·3.36 MB·English

HindawiPublishingCorporation
ComputationalandMathematicalMethodsinMedicine
Volume2013,ArticleID293128,18pages
http://dx.doi.org/10.1155/2013/293128
Research Article
Fluid Structural Analysis of Human Cerebral Aneurysm Using
Their Own Wall Mechanical Properties
AlvaroValencia,1PatricioBurdiles,1MiguelIgnat,1JorgeMura,2
EduardoBravo,2RodrigoRivera,2andJuanSordo2
1DepartmentofMechanicalEngineering,UniversidaddeChile,8370448Santiago,Chile
2InstituteofNeurosurgeryDr.Asenjo,7500691Santiago,Chile
CorrespondenceshouldbeaddressedtoAlvaroValencia;[email protected]
Received29May2013;Revised31July2013;Accepted1August2013
AcademicEditor:NestorV.Torres
Copyright©2013AlvaroValenciaetal.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttributionLicense,
whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
ComputationalStructuralDynamics(CSD)simulations,ComputationalFluidDynamics(CFD)simulation,andFluidStructure
Interaction(FSI)simulationswerecarriedoutinananatomicallyrealisticmodelofasaccularcerebralaneurysmwiththeobjective
ofquantifyingtheeffectsoftypeofsimulationonprincipalfluidandsolidmechanicsresults.EightCSDsimulations,oneCFD
simulation,andfourFSIsimulationsweremade.Theresultsallowedthestudyoftheinfluenceofthetypeofmaterialelementsin
thesolid,theaneurism’swallthickness,andthetypeofsimulationonthemodelingofahumancerebralaneurysm.Thesimulations
usetheirownwallmechanicalpropertiesoftheaneurysm.ThemorecomplexsimulationwastheFSIsimulationcompletelycoupled
withhyperelasticMooney-Rivlinmaterial,normalinternalpressure,andnormalvariablethickness.TheFSIsimulationcoupledin
onedirectionusinghyperelasticMooney-Rivlinmaterial,normalinternalpressure,andnormalvariablethicknessistheonethat
presentsthemostsimilarresultswithrespecttothemorecomplexFSIsimulation,requiringone-fourthofthecalculationtime.
1.Introduction Theaneurysmwallisgenerallycomposedofonlyintimaand
adventitiaoflayeredcollagen.Wallstrengthisrelatedtoboth
Ananeurysmisalocalizeddilationofthewallofanartery; collagen fiber strength and orientation. The average break
it appears most frequently in the abdominal aorta or in strength of aneurysm wall ranged from 0.7MPa to 1.9MPa
the brain vasculature. Intracranial cerebral aneurysms are [2]. Seshaiyer et al. [3] reported the mechanical properties
formed preferentially in abrupt curvatures or bifurcations of cerebral aneurysms. It was found that the lesions were
of arteries belonging to the circle of Willis. In general, its stifferthanpreviouslyreportedexperimentaldataofScottet
geometry resembles a projecting dome on the wall of the al. [4] and To´th et al. [5] and also stiffer as the mechanical
artery. The formation of aneurysms represents the loss of properties of cerebral arteries reported by Monson et al.
thestructuralintegrityofthewall,butthereasonsfortheir [6].
formation and growth are still not clear. A subarachnoid Wall shear stress (WSS) modulates endothelial cell
hemorrhageduetotheruptureofanintracranialaneurysm remodelling via realignment and elongation. Consequently,
isadevastatingeventassociatedwithlargeratesofmorbidity fluid dynamics play important roles in the growth and
andmortality.Approximately12%ofthepatientsdiebefore ruptureofcerebralaneurysms.Aneurysmruptureisrelated
receivingmedicalattention,40%ofhospitalizedpatientsdie toalowlevelofWSSandthereforeisassociatedwithlowflow
within one month of the hemorrhage, and more than one- conditions.Theaneurysmregionwithlowflowconditionsis
third of the patients that survive are left with an important normally the fundus, [7]. High WSS is regarded as a major
neurologicaldeficit[1]. factorinthedevelopmentandgrowthofcerebralaneurysms
A cerebral aneurysm wall has a thin tunica media, and [8].ItisassumedthataWSSofapproximately2Paissuitable
the internal elastic lamina is normally severely fragmented. formaintainingthestructureoftheaneurysmwall,whereas
2 ComputationalandMathematicalMethodsinMedicine
z
𝜇)m 0 0760.454 m)
( 320 1200 𝜇(
900
640
960 600 m)
𝜇
(
(𝜇 1280 300
m)
1600
x
(a) (b)
Figure1:(a)Microtractiontestingmachinemountedunderamicroscope.(b)Measurementofthetissuethicknessofthewallofaneurysm
1[13].
alowerWSSresultsinthedegenerationofendothelialcells Table1:Fittingcoefficientsthroughfive-parameterMooney-Rivlin
viatheapoptoticcellcycle[9]. modelandlinearelasticmodel.
This problem can be investigated using fluid structure
Model Coefficients(MPa)
interaction (FSI) simulations. However, few FSI investiga-
tions have been performed on patient specific aneurysms MooneyRivlin, 𝐶10 =0.3848,𝐶01=−0.0891,
models[10,11].Toriietal.[10]investigatedtheFSIincerebral five-parameters 𝐶11=0.5118,𝐶20=0.5109,𝐶02=0.4912
Linearelastic 𝐸=1.7742
aneurysmsundernormalandhypertensivebloodpressures.
Toriietal.[11]havecalculatedthatthewalldisplacementina
modelwithvariablewallthicknesswas60%largerthanthat
data of the stress versus stretch ratio. The strain energy
oftheuniformwallmodel.
function𝑤isgiveninthemodelby
OurgrouphasdoneFSIsimulationworkontheinterac-
tion between the wall of the aneurysm and the blood. The
𝑤(𝜆)= 𝐶 (𝐼 −3)+𝐶 (𝐼 −3)
simulations have been made using computational geome- 10 1 01 2
trreiseusltsobotfaitnheedsifmroumlataionngsioignrdamicasteofa rreealaltiaonneubreytswmees.nThthee +𝐶11(𝐼1−3)(𝐼2−3)+𝐶20(𝐼1−3)2 (1)
shear stress on the wall caused by the blood flow and the +𝐶 (𝐼 −3)2,
02 2
zones of the aneurysm’s wall that are subjected to greater
stress [12]. All the simulations used mechanical behavior where𝐼1 and𝐼2 arethefirstandsecondstraininvariantsof
models of the aneurysm’s wall developed by other authors the Cauchy-Green deformation tensor 𝐶𝑖𝑗, 𝐶𝑖𝑗 = 2𝜀𝑖𝑗 + 𝛿𝑖𝑗,
[3,4]. where𝛿𝑖𝑗 istheKroneckerdelta;𝐶10,𝐶01,𝐶11,𝐶20,and𝐶02
As part of previous work done by our team [13], the
arematerialconstants,ADINAmanual[14].
mechanicalbehaviorundertractionofacerebralaneurysm
A linear elastic model was also used in which Young’s
extractedfromapatientwascharacterized,andamodelofthe
modulus was obtained from the first two parameters of the
mechanical behavior of the aneurysm’s wall was developed.
Mooney-Rivlinmodel,followingtherecommendationsofthe
Figure1 shows the traction testing machine used and the
ADINAmanual[14]:
measurement of the thickness of one of the samples. The
present work consists in an advanced FSI modeling of the 𝑇=𝐸∗𝑒, (2)
mechanicalbehaviorofcerebralaneurysms,consideringthe
backgroundinformationonthistypeofbehaviorobtainedin 𝐸=6∗(𝐶10+𝐶01), (3)
[13] in terms of experimental characterization and models.
Only a single geometry will be used in the simulations, where 𝑇 is the stress, 𝑒 is the engineering strain, and 𝐸 is
correspondingtothedigitalizationoftheaneurysmusedin Young’smodulusofthematerial.
[13]tocreateitsmodel. Table1 shows the five parameters of the hyperelastic
The simulations were made using the ADINA 8.8.0 Mooney-Rivlin model and Young’s modulus obtained from
commercialsoftware.Afive-parameterhyperelasticMooney- the experimental data from [13]. Figure2(a) shows curves
Rivlinmodelwasusedtocharacterizethehyperelasticbehav- obtainedfromtheexperimentaltestsin[13],withthegreen
ioroftheaneurysmsamples;thismodelwithfiveparameter curvecorrespondingtothefitmadewiththefive-parameter
wasfoundin[13]asthemodelthatbestfittheexperimental Mooney-Rivlin model of the experimental data from [13]
ComputationalandMathematicalMethodsinMedicine 3
1.4 5
1.2
4
1
3
a) 0.8 Pa)
P M
Stress (M 0.6 Stress ( 12
0.4
0.2 0
0 −1
1 1.05 1.1 1.15 1.2 1.25 1.3 0.9 0.95 1 1.05 1.1 1.15
Stretch ratio (𝜆) Stretch ratio (𝜆)
Linear elastic Toth Costalat
Mooney-Rivlin Seshaiyer 1 Valencia
Seshaiyer 2
(a) (b)
Figure2:(a)ModelofthewallmaterialoftheaneurysmfittedthroughMooney-RivlinandthelinearelasticmodelbyContente[13].(b)
HyperelasticmodelsbySeshaiyeretal.[3],To¨thetal,[5],CostalatandSanchez[15]andContente[13].
and the blue curve to that made with the linear elastic oftheaneurysmsthatcanbeusedintheADINAsimulation
model approximation. Figure2(b) shows a comparison of software. For that purpose a case delivered by the Instituto
thehyperelasticcurvemodeledthroughMooney-Rivlinwith de Neurocirug´ıa Asenjo was available. The examinations
othermodelsfoundintheliterature,Seshaiyer[3],To´thetal. were obtained with a Philips Integris Allura 3D Rotational
[5],andCostalatandSanchez[15]. Angiograph, and they deliver dimensional computational
CostalatandSanchez[15]havecharacterizedthemechan- filesinVRMLformat,asseeninFigure3(a).SincetheVRML
icalpropertiesofsixteenintracranialaneurysms,forruptured files cannot be used directly with the ADINA simulation
and unruptured samples. They have found a significant software, a meticulous reconstruction of each case must be
modification in biomechanical properties, the unruptured madetoobtainanadequatefiletomakethesimulation.The
aneurysmsweresoft,theymeasuredonlytoastrainof10%, methoddevelopedin[17]allowsthistasktobeperformedin
and therefore they do not report rupture strain and stress. a couple of hours. The reconstructed CAD model is shown
Valencia’smodelfallswithintherangesdescribedbyCostalat inFigure3(b).Thedifferencebetweentheoriginalgeometry
andSanchez[15]andotherauthors[3,5]. andthereconstructedoneislessthan5%.Theaneurysmsin
Inthisinvestigation,wepresentdetailednumericalsim- thereconstructedgeometryassignedthethicknessmeasured
ulationsoffluidandsolidmechanicsinananatomicallyreal- in [13], while the artery assigned a theoretical thickness
istic cerebral aneurysm model using their own mechanical that corresponds to 10% of its diameter. The aneurysms are
properties.Thepatientspecificgeometrywasreconstructed joined with the artery through a special section that varies
from3Drotationalangiographyimagedata.Thepredictions from the thickness of the aneurysm to the thickness of
using CFD, CSD, and FSI in fluid and solid variables on the artery. Figure4 shows a detail of the variable thickness
the aneurysm 1 are compared. The effects of hypertensive section.
pressureload,aneurysmwallthickness,andwallmodelare The complete geometry is a single volume and is used
alsoreported.Theprincipalgoalofthisworkistoquantify fortheCFDandFSIsimulations.Ahollowgeometryisused
the differences in prediction of fluid and solid variables to carry out the CSD simulations; it has artery thickness
between a time expensive complete FSI simulation with forthearteries,aneurysmthicknessfortheaneurysms,and
CSD simulations. In addition, preprocedural planning for a variable thickness for the junction between artery and
cerebral aneurysm treatments will benefit from an accurate aneurysm. Table2 summarizes the data of interest of the
assessmentofflowpatterns,effectivewallstress,andstrainin completegeometry,andTable3summarizestherelevantdata
theaneurysm,aspresentedherebymeansofcomputational ofbothaneurysms.
simulations. Table4showsthethicknessesoftheaneurysm,theartery,
and the section that joins both of them (aneurysm-artery
2.Methodology junction or A-A junction). The thickness of the artery is
considered to be 10% of its diameter. Some simulations use
2.1.Reconstruction. Tocarryoutthefluiddynamicssimula- halftheaneurysmthicknessshowninTable4,inordertosee
tionsofrealcases,itisnecessarytogenerate3Dgeometries theresultsinathinneraneurysm.
4 ComputationalandMathematicalMethodsinMedicine
Aneurysm 1 Outflow
Aneurysm 2
Inflow
(a) (b)
Figure3:(a)Imageobtainedbyangiography.(b)CADobtainedfromthereconstructionoftheoriginalgeometry.
Table3:Relevantaneurysmdimensions.
Neckdiameter
Aneurysm Width(mm) Length(mm)
(mm)
1(larger) 9.2 16.1 20.8
2(smaller) 3.42 3.33 4.39
Table 4:Thickness oftheaneurysm, theartery,and thejunction
betweenthem.
Section Thickness(mm)
Aneurysm 0.35
A-Ajunction 0.38
Artery 0.4
whereVisthevelocity,𝜌isthefluiddensity,𝜏arethestress,
and𝑓arethebodyforces.
The fluid is considered incompressible, and the flow is
Figure4:Detailofthethicknessoftheaneurysm,theartery,and consideredlaminar.Carreau’smodel(6)isusedtomodelthe
thesectionbywhichtheyarejoined. viscosity𝜇:
𝜇(𝛾)̇ =𝜇 +(𝜇 −𝜇 )(1+𝐴𝛾2̇ )𝑛. (6)
∞ 0 ∞
Table2:Relevantdimensionsforthereconstructedgeometry.
TheCarreaubloodmodelpredictsdecreasingviscosityat
Artery high strain, where 𝜇0 and 𝜇∞ are low and high shear rate
Entrancediameter(mm) 4.98 asymptotic values and the parameters 𝐴 and 𝑛 control the
Entrancearea(mm2) 19.46 transitionregionsize.Wehavetakenthevaluesusedin[16]
TotalCADvolume(mm3) 3862 as 𝜇∞ = 0.00345Ns/m2, 𝜇0 = 0.056Ns/m2, 𝑘 = 10.976,
and𝑚 = −0.3216.Thedensityofbloodwasassumedtobe
constant𝜌=1050kg/m3.
The wall of the aneurysm is considered to be a linear
2.2. Constitutive Equations. The Navier-Stokes equations
elasticmaterial(7)orahyperelasticmaterial(8):
allowmodelingfluids.Theseequationsconsidertheconser-
vationofmass(4)andtheconservationofmomentum(5)in 𝑆=𝐶𝜀, (7)
thefluid:
1 𝜕𝑊
𝑆= ( ), (8)
∇⋅V=0, (4) 2 𝜕𝜀
𝜕V where𝑆isthesecondPiola-Kirchhoffstress,𝜀istheGreen-
𝜌𝜕𝑡 +∇⋅(𝜌VV−𝜏)=𝑓𝐵, (5) Lagrangestrain,and𝑊isthestrainenergyofthematerial.
ComputationalandMathematicalMethodsinMedicine 5
0.7 18000
16000
0.6
14000
0.5
12000
d (m/s) 0.4 ure (Pa) 10000
Spee 0.3 Press 68000000
0.2
4000
0.1
2000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Time (s) Time (s)
Speed pulse Pressure pulse
(a) (b)
Figure5:(a)Bloodvelocitypulseappliedtotheentrancesectionoftheartery.(b)Bloodpressurepulseappliedtotheoutletsectionofthe
artery.
2.3. Boundary Conditions. For the simulations an averaged path.Toreplicatethiseffect,anoscillatingpressureresistance
bloodpulseisusedinthiscase.Inthestudy,Valenciaelal.[12] between80and120(mmHg)inphasewiththeheartbeatsis
usedonepulseforeachpatient.Althoughthisallowsastudy appliedattheoutletsofthesection.Thisallowsthereplication
closertoeachpatient’sreality,itaddsanimportantvariableto oftheclosedcircuitfollowedbythebloodflow.Figure5(b)
beconsideredintheanalysisoftheresults:theheartrate.In showstwopressurepulses.Itisimportanttospecifythatthis
viewofthelargedifferenceinheartbeatfromonepatientto pressureisrelative,notabsolute.
another,itisnotdifficulttounderstandthatthiswillinfluence The boundary conditions for CSD simulations imposed
theresults,andwhencomparingonecasewithanother,itwill on the models was a time-dependent pressure on the inner
notbeknownexactlywhetherthedifferencesfoundaredue wallrepresentativefornormalhumanpressurevariationwith
todifferencesintheheartrateorintheshapeoftheaneurysm. a heart rate of 70beats/min; see Figure5(b). The effects
Thatis,whyitwasdecidedtouseanaveragedpulseinthis of hypertension are reported using the temporal pressure
case.Thisaveragedpulsecannotbejustanyone,becausethe variation between 180mmHg and 100mmHg. The outside
pulse of a healthy person is different from that of a sick pressureduethecerebrospinalfluidwasconsideredconstant
person,soanaveragepulsewasobtainedfromcolorduplex as3(mmHg)(400(Pa))wasusedbyValenciaetal.[12].The
Dopplerimagesof36patientswithcerebralaneurysms. modelwasfixedontheinflowandoutflow.
Tomakeacorrectsimulationofbloodflowintheartery,it CSDsimulationsweremadewithapressurepulse200%
isnecessarytousethecorrectvelocityprofileattheentrance. greaterthanthatshowninFigure5(b),attemptingtoelimi-
For pulsating flows in the arteries the classical parabolic natetheincreasedinternalpressureasacausefortherupture
profile is not sufficient to describe the flow of the velocity ofaneurysms.Thecerebralarteriesarenotfoundinanempty
profile. Womersley’s solution for a pulsating flow in a rigid environment,buttheyareratherimmersedincerebralfluid,
tubehasbeenusedpreviously[16]withgoodresults.Tobe whichproducesaconstantexternalpressurethatcompresses
abletoimplementWomersley’sprofileinthesimulations,use thearteriesradially.
will be made of the method developed in [16], which use On the FSI interface states that (i) displacements of the
theMATLABcalculationsoftwaretodeveloptheprofileand fluid and solid must be compatible, (ii) tractions at this
forlaterexportingtoADINA.Figure5(a)showstheaverage boundarymustbeatequilibrium,and(iii)fluidobeystheno-
velocity profile in two pulses. Numericat tests performed slipcondition.Theseconditionsaregivenasfollows:
by Valencia et al. in [12] have showed that after the second
𝛿 =𝛿 ,
cardiac cycle, the results of velocity and pressure do not 𝑠 𝑓 (9)
changeinsimilargeometriesofcerebralaneurysms,andthe 𝜎 ⋅𝑛̂ =𝜎 ⋅𝑛̂ ,
𝑠 𝑠 𝑓 𝑓 (10)
resultsofthefirstcyclearedifferentbecauseallvariablesstart
thesimulationsatzero. 𝑢=𝑢 ,
𝑔 (11)
To simulate correctlythepath ofthe blood throughthe
cerebral arteries, an outlet condition must be used for the where 𝛿, 𝜎, and 𝑛̂ are displacement, stress tensor, and
bloodflow.Ifanoutletconditionisnotimposed,wewould boundary normal with the subscripts 𝑓 and 𝑠 indicating a
haveaflowthatissetfreegoingoutoftheartery,therefore property of the fluid and solid, respectively. The condition
not representing reality, because the blood follows a closed of (10)doesnotrequireidenticalmatchingmeshesbetween
6 ComputationalandMathematicalMethodsinMedicine
Table5:Numberingandcharacterizationofthesimulationsmade.
Simulation Type Elements Pressure Thicknessaneurysm Material Coupling
1 CFD Tetrahedral Normal — — —
2 CSD Shell Normal Normal L.E —
3 CSD Shell H.T Normal L.E —
4 CSD Shell Normal 1/2 L.E —
5 CSD Shell Normal Normal M.R —
6 CSD Shell H.T Normal M.R —
7 CSD Shell Normal 1/2 M.R —
8 CSD Tetrahedral Normal Normal L.E —
9 CSD Tetrahedral Normal Normal M.R —
10 FSI Tetrahedral Normal Normal L.E Complete
11 FSI Tetrahedral Normal Normal L.E Onedirection
12 FSI Tetrahedral Normal Normal M.R Complete
13 FSI Tetrahedral Normal Normal M.R Onedirection
H.Tishypertensionbloodpressurerangedbetween13200Paand23800Pa.NormalpressureisthebloodpressurepulseshowninFigure5(b),L.E.isthelinear
elasticwallmodel,andM.R.isthehyperelasticMooney-Rivlinwallmodel.
the two domains and instead supports the use of solution Plane 4
mappingtoestablishtheequilibrium.
Plane 3
2.4.GeneralConsiderationsoftheSimulations. Thirteensim-
ulationsweremadeinthisresearch,andtheircharacteristics
are summarized in Table5. One CFD simulation (rigid
walls),eightCSDsimulations,andfourFSIsimulationswere
Plane 2
performed.
CSD simulations 2, 3, and 4 use linear elastic material,
CSDsimulations5,6,and7usehyperelasticMooney-Rivlin
material with the approximation of shell elements; in these
Plane 1
simulationweinvestigatethewallthicknessandhypertensive
pressureeffects.IntheCSDsimulations8and9weinvestigate
theeffectof3Dtetrahedralelements.
TheFSIsimulations10and11uselinearelasticmaterial;
we investigate the effect of the method coupling, and in
the FSI simulations 12 and 13 we investigate the effect of
the method coupling using a hyperelastic Mooney-Rivlin
Figure 6: Control planes. Plane 1 or entrance plane, plane 2 or
material.
middleplane,plane3orupperplane,andplane4ortransverseplane
The coupling indicates how the solid and fluid models
ofaneurysm1.
interactwithoneanother.Itispossibletohaveacompletely
coupledFSIsimulationinwhichanequilibriumofstrength
anddisplacementsbetweensolidandfluidissoughtforevery
time step, requiring several iterations between them and Itisimportanttospecifythatthecontrolelementcontains
thereforelongercalculationtime.ItisalsopossibletouseFSI all the control points, so that the control element contains
simulations coupled in one direction. The results obtained thedataforthefournodesoftheelement.Thedatumforthe
in theCFD simulationareappliedto thesolidat each time controlpointisthenodewiththegreatervalues,becausein
step,butinthiscasethesoliddoesnotaffectthefluid.This thecaseoftheCFDthenodesthatareonthewall(1,2,and
decreases calculation time and the precision of the results. 3 in Figure7) always have values equal to zero. In the case
MeshtestsweremadefortheCSDandCFDsimulations. of CSD simulations with shell-type elements the element is
Figure6showsthecontrolplanesforaneurysm1.Plane triangular,soitonlyhasnodes1,2,and3thatareshownin
4 is a transverse plane of aneurysm 1 perpendicular to its Figure7.
entranceareaandparalleltothedirectionofflowintheartery. The rupture point in saccular cerebral aneurysm is
Planes1,2,and3aretheentrance,middle,andupperplanes locatedinthezonewithmaximumstressanddeformationor
of aneurysm 1, and they are parallel to the entrance area of intheaneurysmfundus;forthisreasonwedefineandstudy
aneurysm1.Thecontrolpointsforeachplanecorrespondto indetailstheresultsonthecontrolpointshowninFigure7
theelementsubjectedtothemaximumvalueofthevariable related with the wall shear stress. Averaged values on the
inquestion. aneurysmsurfacedonotprovideusefulinformationtostudy
ComputationalandMathematicalMethodsinMedicine 7
Control point on
aneurysm 1 surface
3
2
1
4
Figure7:Controlpointinthefundusofaneurysm1.
Table6:Meshdensityandcharacteristicsoftheelements.
Element Geometry Nodes/element Elementsize(mm) Meshdensity(elements/mm3)
3Dsolid Tetrahedral 4 0.2 1400
3Dfluid Tetrahedral 4 0.33 110
Shell Triangular 3 0.33 25
thegrowthandruptureofsaccularaneurysms.Wereported 3.Results
thelocalmaximumvaluesofseveralvariablesonaneurysm1
andthetemporalvariations. Thefollowingresultsanddiscussionsarebasedonaneurysm
1 due to this aneurysm that is the principal pathology in
this patient. The aneurysm 2 is a part of the computational
2.5.NumericalMethod. TheCFD,FSI,andCSDmodelswere domain, ant it is located after the principal aneurysm, so
solvedbyacommercialfinite-elementpackageADINA8.8.0. that the secondary aneurysm does not affect the aneurysm
The Finite-Element Method (FEM) is used to solve the 1 or principal pathology. The results of the aneurysm 2 are
governingequations.TheFEMdiscretizesthecomputational not discussed due to this an incipient geometrical change
domain into finite elements that are interconnected by ele- of the artery, and in this stage of development it cannot be
mentnodalpoints.ThefluiddomainemploysspecialFlow- classifiedasfullsaccularaneurysm.Thepatientisclinically
Condition-Based-Interpolation(FCBI)tetrahedralelements. treatedduetothepresenceoftheaneurysm1orpathology;
Wehaveusedtheformulationwithlargedisplacementsand thesecondaneurysmisinthiscasenotpartofthetreatment.
small strains in the FSI calculation available in ADINA. Theaneurysm2isnotanalyzed.
The unstructured grids were composed of tetrahedral with The differences between the FSI simulations and the
four-nodeelementsinthefluidandfour-nodeisoparametric CFDandCSDsimulationswerestudiedwiththepurposeof
elementsfortheshellstructureofthesolid. observingthedifferencesthataregeneratedbynotconsider-
We have performed a sensitivity analysis of the mesh ingtheinteractionbetweenthesolidandthefluid.
size in CFD and CSD simulations. Using a grid size with Theresultsareshownintwoways.Thefirstcorresponds
3
110elements/mm ,thedifferenceonpredictionsofmaximum to distribution graphics of simulation 12 only, because it
wall shear stress at peak systole respect to a grid with was considered as the most complete. The second way of
130elements/mm3wereonly1%.FortheCSDsimulationswe showing the results is by means of temporal graphics of
haveused1400elements/mm3,andthedifferencewithagrid significantvariablesonspecificcontrolpointsonaneurysm
3 1consideringseveralsimulations.Thesignificantvariablesin
of1760elements/mm onmaximumdisplacementwasonly
thefluidarepressure,wallshearstress,andvelocity,whilein
1%. The mesh density used in each case is summerized in
thesolidtheyarethedisplacementofthewalloftheaneurysm
Table6.Fortheintegrationofthistime-dependentproblem,
1, Von Mises effective stress, and first principal stress. The
we used the full implicit Euler method with a time step of
Δ𝑡 = 0.01s.Thetoleranceforalldegreesoffreedomwasset distribution graphics are shown at the times at which the
significantvariablesaremaximumorminimum.
to0.001.
Figure8showsthepressuredistributionoverthegeom-
Theworkstationusedtoperformthesimulationsisbased
etryfor0.92(s)and1.2(s)forsimulation12.Inbothcasesit
onanIntelXeondualcore64bitsprocessorof3.0Ghzclock
canbeseenthatthepressuredropsinthedirectionoftheflow.
speed,8.0GbRAMmemory.ThesimulationtimefortheFSI
Thepressuredropis∼10(kPa)duringsystole.Thepressureon
case based on 2 pulsatile flow cycles was around 29CPU
aneurysmisconstant.
hours.
8 ComputationalandMathematicalMethodsinMedicine
Pressure Pressure
RST calc RST calc
Time2 Time2.28
12900 24667
12600 23333
12300 22000
12000 20667
11700 19333
11400 18000
11100 16667
Maximum Maximum
12972 24992
EG1, EL145793, IPT111 EG1, EL119478, IPT111
Minimum Minimum
10935 15743
EG1, EL36742, IPT11 EG1, EL137734, IPT 111
(a) (b)
Figure8:Simulation12.Pressuredistributioninthecompletegeometryduringdiastole(0.89(s))ontheleftandduringsystole(1.16(s))on
theright.
×104 1isoneorderofmagnitudelowerthanthemaximumvelocity
Max an. 1
3 attheentrancetothisaneurysm.
Thetemporalevolutionofthevelocityintheaneurysm1
2.5 inthreeplanesandinthecontrolpointisshowninFigure12;
theeffectsofthemodelareveryimportantinthemagnitude
2
a) ofthevelocity.ThedifferencesbetweenCFDsimulation1and
P
Pressure ( 1.51 dFiSaIsFstoiimgleuuralenat1di3osnsyhssot1o0wl,es1i1tn,h1et2hw,eaancllodsm1h3peaalerrteesrtgreeelesovsmadneittsr.tryi.bIuttiisosneednurtihnagt
bothaneurysmspresentlowerwallshearstresscomparedto
0.5 the rest of the geometry. Figure14 shows in detail the wall
shearstressinaneurysm1withothercolormapranges.The
0 greatest wall shear stress is concentrated on the wall where
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
thebloodflowenterstheaneurysm.Thetemporalevolution
Time (s)
ofwallshearstressatthecontrolpointofaneurysm1isshown
Lin O-W M.R. inFigure15;theeffectsofthemodelareveryimportantinthe
Linear FL wallshearstress.TheFSIsimulations12and13usingMooney-
M.R. O-W Rivlin show similar values of wall shear stress, except at
systole,andtheCFDsimulation1showslowwallsheartstress.
Figure 9: Temporal evolution of the maximum pressure in
Figure16 shows the displacement distribution for the
aneurysm1forsimulation1(FL),simulation10(linear),simulation
wall of the geometry during systole. It can be seen that the
11(LinO-W),simulation12(M.R),andsimulation13(M.RO-W).
maximumdisplacementoccursontheleftsideofaneurysm1
withrespecttotheentrance-outletdirectionofthebloodflow.
Themaximumdisplacementis4.1(mm),whichisimportant
The temporal variation of pressure shown in Figure9 consideringthesizeoftheaneurysm1.
shows that the differences on pressure between the CFD The temporal evolution of maximum displacement in
simulation 1 and FSI simulations 10, 11, 12, and 13 are low. aneurysm1forCSDsimulations2,5,8,and9andFSIsimula-
Figure10 shows the velocity distribution for the transverse tions10,11,12,and13isshowninFigure17;theCSDpredicts
andentrance(plane1)controlplanes,whileFigure11shows lowerwalldisplacementscomparedwiththeFSIresults.
themiddle(plane2)andupper(plane3)planesofaneurysm Figure18 shows the distribution of effective Von Mises
1duringsystole.Inthetransverseplane(Figure10)itisseen stress on the wall of the artery and of aneurysm 1. The
how the blood enters aneurysm 1 and goes to the left side high stress zone occurs near the dome of aneurysm 1. In
of the wall due to the direction of the flow. In the entrance aneurysm 1 the maximum stress occurs in the same zone
(Figure10), middle and upper (Figure11) planes, it is seen as the maximum displacement. The overall maximum is
that the flow loses velocity as it starts recirculating around 100(kPa) higher than the maximum stress in aneurysm 1,
aneurysm1.Themaximumvelocityattheupperofaneurysm whichis713(kPa).
ComputationalandMathematicalMethodsinMedicine 9
Velocity Velocity
magnitude magnitude
Time2.28 Time2.28
0.39 0.65
0.33 0.55
0.27 0.45
0.21 0.35
0.15 0.25
0.09 0.15
0.03 0.05
Maximum Maximum
0.4322 0.6946
Node44854(0.4346) Node40391(0.6743)
Minimum Minimum
0.000928 0.0002415
Node278(0.0009467) Node14056(0.0002373)
(a) (b)
Figure10:Simulation12.Distributionofvelocityatthetransversecontrol(left)andentranceplanes(right)foraneurysm1atpeaksystole
(1.16(s)).
Velocity Velocity
magnitude magnitude
Time2.28 Time2.28
0.2925 0.08667
0.2475 0.07333
0.2025 0.06
0.1575 0.04667
0.1125 0.03333
0.0675 0.02
0.0225 0.00667
Maximum Maximum
0.3061 0.09558
Node50380(0.3033) Node37486(0.09408)
Minimum Minimum
0.0002972 0.001252
Node18849(0.0003648) Node17228(0.001271)
(a) (b)
Figure11:Simulation12.Distributionofthevelocityinthemiddle(left)andupper(right)controlplanesforaneurysm1atpeaksystole
(1.16(s)).
The temporal evolution of the effective maximum Von principalmaximumstressisconcentratedcoincideswiththe
Misesstressinaneurysm1forCSDsimulations2,5,8,and areaonwhichthemaximumeffectiveVonMisesstressand
9andFSIsimulations10,11,12,and13isshowninFigure19; themaximumdisplacementareconcentrated.
theCSDpredictslowerVonMisesstresscomparedwiththe The temporal evolution of the first principal maximum
resultsofFSImodelswithMoney-Rivlin. stressinaneurysm1forCSDsimulations2,5,8,and9and
Figure20 shows the distribution of the first principal FSI simulations 10, 11, 12, and 13 is shown in Figure21; the
stress on the wall of the artery and of aneurysm 1 during CSDpredictssimilarstresscomparedwiththeresultsofFSI
systole.Again,themaximumstressforthewholegeometry modelsexceptforthesimulation9.
occurs in the neck of aneurysm 1 and is 150(kPa) higher Figure22 shows the first principal stretching for the
thanthefirstprincipalmaximumstressofaneurysm1,which geometry and aneurysm 1 under systole. For the geometry
is 725(kPa). For aneurysm 1 the area on which the first the maximum stretching is concentrated in the neck of
10 ComputationalandMathematicalMethodsinMedicine
0.8 0.35
0.7 0.3
0.6
m/s) 0.5 m/s) 00.2.25
d ( 0.4 d (
pee 0.3 pee 0.15
S S
0.1
0.2
0.1 0.05
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Time (s) Time (s)
Lin O-W M.R. Lin O-W M.R.
Lin FL Lin FL
M.R. O-W M.R. O-W
(a) Plane1 (b) Plane2
0.16 0.045
0.14 0.04
0.12 0.035
m/s) 0.1 m/s) 00.0.0235
d ( 0.08 d ( 0.02
e e
e 0.06 e
p p 0.015
S S
0.04 0.01
0.02 0.005
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Time (s) Time (s)
Lin O-W M.R. Lin O-W M.R.
Lin FL Lin FL
M.R. O-W M.R. O-W
(c) Plane3 (d) ControlPoint
Figure12:Temporalevolutionofvelocityinaneurysm1intheentranceplane(a),middleplane(b),upperplane(c),andcontrolpoint(d)
forsimulation1(FL),simulation10(Lin),simulation11(LinO-W),simulation12(M.R),andsimulation13(M.RO-W).
Smoothed max Smoothed max
shear stress shear stress
RST calc RST calc
Time2 Time2.28
14 32.5
12 27.5
10 22.5
8 17.5
6 12.5
4 7.5
2 2.5
Maximum Maximum
15.47 35.78
Node7591 Node12682
Minimum Minimum
0.0009514 0.0008122
Node137 Node137
(a) (b)
Figure13:Simulation12.Distributionofwallshearstressduringdiastoleontheleftandsystoleontheright.