Table Of ContentHIGH ORDER FINITE ELEMENT SIMULATIONS FOR FLUID
DYNAMICS VALIDATED BY EXPERIMENTAL DATA FROM THE
FDA BENCHMARK NOZZLE MODEL
V.Chabannes,C.Prud’homme,M.Szopos,andR.Tarabay
Universite´ deStrasbourg,CNRS,IRMAUMR7501,F-67000Strasbourg,France,
{chabannes,prudhomme,szopos,tarabay}@math.unistra.fr
SUMMARY
7 Theobjectiveofthepresentworkistoconstructasoundmathematical,numericalandcomputational
1
framework relevant to blood flow simulations and to assess it through a careful validation against
0
2 experimental data. We perform simulations of a benchmark proposed by the FDA for fluid flow in
an idealized medical device, under different flow regimes. The results are evaluated using metrics
n
a proposedintheliteratureandthefindingsareinverygoodagreementwiththevalidationexperiment.
J
1
Keywords: CFD,validation,medicaldevice,opensourcefiniteelementsoftware
1
]
P
A
1 INTRODUCTION
.
h
t
a A challenging benchmark was proposed by the US Food
m
and Drug Administration (FDA) in [1] in order to assess
[ the stability, accuracy and robustness of computational
2 methods in different physiological regimes. The findings
v of 28 blinded investigations were reported in [2] and as
9
criticallyanalyzedin[3],practicallyallCFDsolversfailed
7
1 topredictresultsthatagreedinasatisfactorymannerwith
2 the experimental data. Several subsequent papers tack-
0
led this question, by employing different numerical ap-
.
1
proaches: forinstanceafinite-elementbaseddirectnumer-
0
7 ical simulation method in [4] or a large-eddy simulation
1 methodin[5].
:
v
Weaimatcontributingtotheeffortofimprovingtherelia-
i
X bilityandreproducibilityofcomputationalstudiesbyper-
r formingathoroughvalidationofthefluidsolverdeveloped
a
intheopensourcefiniteelementlibraryFeel++[6]. Inthe
Figure 1: Computational domain (top)
current investigation, we present results corresponding to
andflowregimespecifications(bottom);
three Reynolds numbers 500,2000 and 3500 obtained by
Re and Re : Reynolds number in the
i t
usingadirectnumericalsimulationmethodoftheNavier-
inlet section and throat section, respec-
Stokesequations. Inparticularweimplementandcompare
tively.
low order as well as high order approximations including
for the geometry and we discuss some issues not previ-
ouslyreportedintheliterature.
2 METHODOLOGY
Benchmark description. The FDA benchmark nozzle model provides a comprehensive dataset of
experimental measures using a well-defined geometry corresponding to an idealized medical device
(see Figure 1 for a schematic sketch of the domain and [1, Sec. 2.1] for the precise dimensions of
eachpart). Fivesetsofdataspanninglaminar,transitionalandturbulentregimesaremadeavailable;
wefocusinthecurrentworkontheflowregimespecificationsdescribedinFigure1.
Thecomparisonwithexperimentaldataismadeintermsof(i)wallpressuredifference(normalizedto
themeanthroatvelocity)versusaxialdistance; and(ii)axialcomponentofthevelocity(normalized
tothemeaninletvelocity)alongthecenterline:
p −p u 4Q 4Q
∆pnorm = z z=0 and unorm = z, where u2 = , u = , (1)
1ρ u2 z u t πD2 i πD2
2 f t i t i
and Q is the volumetric flow rate retrieved from Re (see Figure 1, right panel). Furthermore, we
t
presentresultsontwovalidationmetricsreportedin[2],alsoassessedin[4]: aconservationofmass
errormetricE (onapercentagebasis)andageneralvalidationmetricE comparingaverageexper-
Q z
imentalvelocitydatawithcomputedaxialvelocities.
Fluid equations and numerical approach. We now turn to the mathematical and the numerical
setting. We consider the homogeneous, incompressible, unsteady Navier-Stokes equations, which
readinconservativeform: find(u,p)suchthatρ(cid:0)∂u +(u·∇)u(cid:1)−µ∆u+∇p = 0, div(u) = 0,
∂t
in Ω × I. The set Ω ⊂ R3 represents the spatial domain described in Figure 1. I = (0,T) is the
time interval, u and p are the velocity and pressure of the fluid and ρ and µ are the density and the
dynamic viscosity of the fluid, respectively. We supplement the equations with initial and boundary
conditions. Att = 0s,thefluidisconsideredtobeatrest,u(x,t) = 0. APoiseuillevelocityprofile
isimposedonΓ ,homogeneousDirichletconditiononΓ andafreeoutflowonΓ .
inlet wall outlet
Wereferto[8,Sec. 2]regardingthevariationalformulation,thefiniteelementdiscretizationincluding
lowtohighordergeometryaswellasthetimediscretization. WechoosethegeneralizedTaylor-Hood
finite element for the velocity-pressure discretization; the notation P P G is used to specify
N+1 N kgeo
exactlythediscretizationspacesforthevelocity,pressureandgeometry,respectively.
The benchmark hereafter is developed in the framework of the Finite Element Embedded Library in
C++, Feel++[6], that allows to use a very wide range of Galerkin methods and advanced numeri-
cal techniques such as domain decomposition. The ingredients include a very expressive embedded
language, seamless interpolation, mesh adaption and seamless parallelization. Regarding the com-
putational domain, we used GMSH. The construction used the following steps: (i) start with a 2D
geometry embedding the benchmark metric locations and customizing characteristic mesh size de-
pending on the region and (ii) extrude by rotation to obtain the device geometry. Finally we use the
PETSC interfacedevelopedinFeel++andinparticulartheFieldSplitpreconditioningframeworkto
implementblockpreconditioningstrategiessuchasPCD[7]. NotethatPCDrequiresspecifictuning
withrespecttoboundaryconditions.
3 RESULTSANDCONCLUSIONS
We perform simulations for three
Reynolds numbers evaluated in the
h h h N
min max average elt
throat Re = 500,2000,3500, with
t
several mesh refinements and poly- M0 1.9·10−4 2.9·10−3 1.3·10−3 412575
nomial order approximations. The M1 1.6·10−4 1.8·10−3 7.6·10−4 830000
fluid’s prescribed density is ρ = M2 1.4·10−4 1.96·10−3 6·10−4 3400000
1056kg andviscosityµ = 0.0035Pa.s. M3 8.5·10−5 1.7·10−3 3.5·10−4 7000000
The mm3esh characteristics are de- M4 6.3·10−5 2.0·10−3 5.8·10−4 2879365
scribedinTable1. AtRe = 500,the M5 1.4·10−4 2.6·10−3 4.1·10−4 3200000
t
simulation is carried out until t = 3s,
Table 1: Characteristic lengths of the different meshes:
time reasonably close to the steady
h ,h ,h are respectively the minimum, maxi-
state, and we choose the time step min max average
equal to ∆t = 10−3. At Re = 2000 mum and average edge length in the meshes and Nelt is the
t
numberoftetrahedra.
(respRe = 3500), thenumericalex-
t
periments were carried out until t =
0.45s(respt = 0.4s),timewhentheturbulentregimewasfullydevelopedandweset∆t = 10−4.
Figure2showstheresultsinthethreeflowregimesforthenormalizedaxialvelocityandthenormal-
izedpressuredifferencealongthezaxis,respectively. Ineachcase,wecanseesatisfactoryagreement
withtheexperimentaldata. However,forRe = 2000,weobservethatthenumericaljetbreakdown
t
pointiscapturedfurtherdownstreamthantheexperimentallyobservedbreakdownpoint. Asrecently
highlighted in [5], the prediction of the axial location of the jet breakdown is extremely sensitive
to numerical parameters, therefore a possible explanation of this mismatch may be the accuracy of
the numerical integration. Finally, we illustrate in Figure 3 the computation of metrics E and E
z Q
for several mesh refinements at Re = 500. The metric E takes small values in each numerical
t z
experiment, identifying a good agreement betweeen computed and experimental data, and displays
only small variations with respect to mesh refinement. On the other hand, the metric E is more
Q
sensitivetothisfactor: errordoesn’texceedthe∼ 2%exceptforthecoarsemeshM0where,intwo
locations, the error increases up to ∼ 10%. Furthermore, we note that the P P G approximation
3 2 1
doesn’t improve the results for the coarse mesh, but that a satisfactory error below 2% is retrieved
whenusingaP P G approximation. Additionalteststocomplementthestudyoftheimpactofhigh
2 1 2
orderapproximationareongoing.
Conclusions and perspectives We validated our computation fluid dynamic framework against this
FDA benchmark for three different regimes and different discretization and solution strategies. Per-
spectivesincludeafullreportonourfindingsincludingintermsofiterationandtimingperformances
aswellextendingourresultstotheturbulentrange.
Acknowledgments The authors wish to thank Mourad Ismail for the fruitful discussions we had.
Moreover we would like to acknowledge the support of (i) Center of Modeling and Simulation of
Strasbourg (Cemosis), (ii) the ANR MONU-Vivabrain (iii) the LabEx IRMIA and (iv) PRACE for
awarding us access to resource Curie based in France at CCRT as well as GENCI for awarding us
accesstoresourceOccigenbasedinFranceatCines
REFERENCES
[1] P. Hariharan, M. Giarra,V. Reddy, S.W. Day, K.B. Manning, S. Deutsch, S.F.C. Stewart, M.
Y.Myers,M.R.Berman,G.W.Burgreen,E.GPaterson,andR.A.Malinauskas.Multilaboratory
particle image velocimetry analysis of the FDA benchmark nozzle model to support validation
ofcomputationalfluiddynamicssimulations.JournalofBiomechanicalEngineering,133(4):1–
14,2011.
[2] S.F.C. Stewart, P. Hariharan, E.G. Paterson, G.W. Burgreen, V. Reddy, S.W. Day, K.B. Man-
ning, S. Deutsch, M.R. Berman, M. Y. Myers and R.A. Malinauskas. Results of FDAs first
interlaboratorycomputationalstudyofanozzlewithasuddencontractionandconicaldiffuser.
CardiovascularEngineeringandTechnology,4(4):374–391,2013.
[3] F. Sotiropoulos. Computational fluid dynamics for medical device design and evaluation: are
wethereyet?.CardiovascularEngineeringandTechnology,3(2):137–138,2012.
[4] T.Passerini,A.Quaini,U.Villa,A.VenezianiandS.Canic.Validationofanopensourceframe-
workforthesimulationofbloodflowinrigidanddeformablevessels.Internationaljournalfor
numericalmethodsinbiomedicalengineering,29(11):1192–1213,2013.
[5] V. Zmijanovic, S. Mendez, V. Moureau and F. Nicoud. About the numerical robustness of
biomedical benchmark cases: Interlaboratory FDA’s idealized medical device. International
journalfornumericalmethodsinbiomedicalengineering,33(1),2017.
[6] C. Prud’homme, V. Chabannes, V. Doyeux, M. Ismail, A. Samake and G. Pena. Feel++: A
Computational Framework for Galerkin Methods and Advanced Numerical Methods. ESAIM:
Proceedings,EDPSciences,38:429–455,2012.
[7] H.C. Elman, D.J. Silvester and A.J. Wathen. Finite elements and fast iterative solvers: with
applicationsinincompressiblefluiddynamics.OxfordUniversityPress(UK),2014.
data1 data2 6
e
15 data3 data4 nc
e
data5 M0-P2P1G2 er
uz M3-P2P1G1 4 diff
d re
e 10 u
z 2 s
mali pres
nor data1 data2 0 zed
5 data3 data4 ali
m
data5 M0-P2P1G2
r
−2 o
M3-P2P1G1 n
−0.1 −5·10−2 0 5·10−2 0.1−0.1 −5·10−2 0 5·10−2 0.1
arclength arclength
15
data1 4
e
data2 nc
e
uz 10 ddaattaa34 3 differ
d M4-P2P1G1 re
e u
z 2 s
mali data1 pres
nor 5 ddaattaa23 1 zed
data4 ali
m
data5 0 nor
M4-P2P1G1
0
−0.1 −5·10−2 0 5·10−2 0.1 −0.1 −5·10−2 0 5·10−2 0.1
arclength arclength
data1
12 data2 2 nce
e
uz 10 ddaattaa34 1.5 differ
d 8 data5 re
e u
z M5-P2P1G1 1 s
mali 6 data1 pres
nor 4 ddaattaa23 0.5 zed
data4 ali
m
2 data5 0 or
n
M5-P2P1G1
−0.1 −5·10−2 0 5·10−2 0.1 −0.1 −5·10−2 0 5·10−2 0.1
arclength arclength
Figure 2: Comparison between experimental data and numerical results for the normalized axial
velocity along z (left) and the normalized pressure difference along z (right), for Re = 500 (top)
t
Re = 2000(middle)andRe = 3500(bottom).
t t
101
M0-P2P1G1
M0-P2P1G2
M0-P3P2G1 0
M1-P2P1G1
100 M2-P2P1G1
z M3-P2P1G1 q
E M0-P2P1G1 −5 E
M0-P2P1G2
M0-P3P2G1
10−1 M1-P2P1G1
M2-P2P1G1 −10
M3-P2P1G1
−0.1 −5·10−2 0 5·10−2 −0.1 −5·10−2 0 5·10−2
z z
Figure3: ValidationmetricsE (left)andE (right)forRe = 500.
z Q t
[8] C. Caldini Queiros, V. Chabannes, M. Ismail, G. Pena, C. Prud’Homme, M. Szopos and
R. Tarabay. Towards large-scale three-dimensional blood flow simulations in realistic geome-
tries.ESAIM:Proceedings,EDPSciences,43:195–212,2013.
