Author's Accepted Manuscript
Analysis of thoracic aorta hemodynamics using 3D particle tracking velocimetry and computational fluid dynamics Diego Gallo, Utku Gülan, Antonietta Di Stefano, Raffaele Ponzini, Beat Lüthi, Markus Holzner, Umberto Morbiducci
www.elsevier.com/locate/jbiomech
PII: DOI: Reference:
S0021-9290(14)00356-X http://dx.doi.org/10.1016/j.jbiomech.2014.06.017 BM6699
To appear in:
Journal of Biomechanics
Accepted date: 12 June 2014 Cite this article as: Diego Gallo, Utku Gülan, Antonietta Di Stefano, Raffaele Ponzini, Beat Lüthi, Markus Holzner, Umberto Morbiducci, Analysis of thoracic aorta hemodynamics using 3D particle tracking velocimetry and computational fluid dynamics, Journal of Biomechanics, http://dx.doi.org/10.1016/j.jbiomech.2014.06.017 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
Analysis of Thoracic Aorta Hemodynamics Using 3D Particle Tracking
2
Velocimetry and Computational Fluid Dynamics
3 4 5
Diego Gallo1, Utku Gülan2, Antonietta Di Stefano1, Raffaele Ponzini3, Beat Lüthi2, Markus Holzner2,
6
Umberto Morbiducci1
7 8
1 Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Turin, Italy
9
2 Institute of Environmental Engineering, ETH Zurich, Zurich, Switzerland
10
3 HPC and Innovation Unit, CINECA, Milan, Italy
11 12 13
Address for correspondence:
14
Umberto Morbiducci, Ph.D
15
Department of Mechanical and Aerospace Engineering, Politecnico di Torino
16
Corso Duca degli Abruzzi, 24
17
10129 Turin, Italy
18
Tel.: +39 011 0903394
19
Fax: +39 011 5646999
20
E-mail:
[email protected] 21 22 23
24
ABSTRACT
25 26
Parallel to the massive use of image-based computational hemodynamics to study the complex flow
27
establishing in the human aorta, the need for suitable experimental techniques and ad hoc cases for the
28
validation and benchmarking of numerical codes has grown more and more.
29
Here we present a study where the 3D pulsatile flow in an anatomically realistic phantom of human
30
ascending aorta is investigated both experimentally and computationally. The experimental study uses
31
3D particle tracking velocimetry (PTV) to characterize the flow field in vitro, while finite volume
32
method is applied to numerically solve the governing equations of motion in the same domain, under
33
the same conditions. Our findings show that there is an excellent agreement between computational
34
and measured flow fields during the forward flow phase, while the agreement is poorer during the
35
reverse flow phase.
36
In conclusion, here we demonstrate that 3D PTV is very suitable for a detailed study of complex
37
unsteady flows as in aorta and for validating computational models of aortic hemodynamics.
38
In a future step, it will be possible to take advantage from the ability of 3D PTV to evaluate velocity
39
fluctuations and, for this reason, to gain further knowledge on the process of transition to turbulence
40
occurring in the thoracic aorta.
41 42
Keywords: Aortic Flow; Computational Fluid Dynamics; Particle Tracking Velocimetry;
43
Hemodynamics; Ascending Aorta.
44 45
46
INTRODUCTION
47 48
The complex hemodynamics observed in the human aorta play important roles in the regulation of
49
vascular homeostasis and the localization of atherosclerosis [Debakey et al. 1985]. In recent years, the
50
coupling of imaging techniques and computational fluid dynamics (CFD) has been applied to
51
reconstruct highly resolved blood flow patterns in more and more realistic and fully personalized
52
simulations of thoracic aorta hemodynamics [Liu et al., 2011; Lantz et al., 2012; Brown et al., 2012;
53
Morbiducci et al. 2013; Caballero & Laín, 2013], thus allowing the study of disease mechanisms and
54
aiding the design and evaluation of medical devices. However, computational methods require
55
assumptions and operator-dependent specifications (like temporal and spatial discretization) that might
56
influence the predicted hemodynamic scenario [Morbiducci et al. 2010, 2011a; Gallo et al. 2012;
57
Valen-Sendstad & Steinman 2013]. The accuracy of CFD simulations despite these issues can be
58
evaluated by the cross-validation of numerical results with experimental results.
59
Among the techniques that can be used to such purpose, one possibility is given by 4D phase contrast
60
MRI (PC MRI), recently applied in vivo to access the velocity field in the aorta [Kilner et al. 1993;
61
Hope et al. 2008; Frydrychowicz et al. 2011; Morbiducci et al. 2011b; Hope et al. 2013] and to
62
validate numerical results [Marshall et al. 2004; Wen et al. 2010; Kung et al. 2011; Edelhoff et al.
63
2013; Lantz et al. 2013]. However, PC MRI usually requires long acquisition times and suffers from
64
low spatial and temporal resolution and a low signal-to-noise ratio. Anemometric techniques have
65
been also applied for in vitro characterization of the fluid dynamics in aortic phantoms [Browne et al
66
2000; Balducci et al. 2004; Dasi et al. 2007; Doyle et al. 2008; Chen et al. 2014]. Among them, 2D
67
and stereoscopic Particle Image Velocimetry achieve high spatial resolution but are inherently limited
68
to bi-dimensional observation planes. Laser Doppler velocimetry can achieve high spatial and
69
temporal resolution but it is a point measurement technique and it is time consuming to fully resolve
70
the 3D velocity field in a large domain such as the human aortic phantom. An alternative is
71
represented by 3D Particle Tracking Velocimetry (PTV), an optical technique based on imaging of
72
flow tracers. It has been successfully used to obtain Lagrangian velocity fields in a wide range of
73
complex and turbulent flows [Malik et al. 1993; Balducci et al. 2004; Cenedese et al. 2005; Lüthi et al.
74
2005; Holzner et al. 2008; Boutsianis et al. 2009], and very recently applied to characterize fluid
75
structures in the ascending aorta [Gülan et al. 2012; Knobloch et al. 2013].
76
In this study, the 3D unsteady flow in an anatomically realistic rigid phantom of human ascending
77
aorta (AAo) is investigated both experimentally and computationally. The feasibility of unsteady
78
velocity measurements by 3D PTV, when applied to aortic flows, is demonstrated, and the in vitro
79
measured data are used as boundary conditions in high-resolution CFD simulations to investigate the
80
hemodynamics numerically. The numerical solution is then cross-validated by PTV experiments.
81
82 83
METHODS
84
Anatomical Model Construction
85
A high-resolution MRI scan of a healthy subject was performed to construct a 3D computer
86
anatomical model of an aorta. The anatomical model was used (1) to physically manufacture a silicon-
87
made rigid phantom (Elastrat, Geneva, Switzerland) for the in vitro experiment [Gülan et al. 2012] and
88
(2) as geometry input for the CFD study.
89 90
In Vitro Setup
91
The experimental setup resembles the one recently proposed by Gülan and coworkers (2012). Briefly,
92
it comprises the aortic phantom, a ventricular assist device (VAD, MEDOS, Stolberg, Germany), a
93
pump system to drive the VAD and the optical part. A scheme of the test bench is presented in Figure
94
1. Technically, the VAD, used to mimic the function of the heart by pumping the working fluid into
95
the aortic phantom, was driven by a custom built PWG-01-ST pressure and vacuum pump (Innovative
96
Technology and Engineering Solutions, Haifa, Israel). To avoid optical disturbances affecting the
97
measurements, matching the index of refraction of the materials is crucial. In this study, a mixture of
98
37 % glycerol, 48 % water and 15 % NaCl (kinematic viscosity equal to 4.85 x 10-6 m2/s) was used to
99
match the refractive index of silicon (n=1.41) [Gülan et al. 2012].
100
The hydraulic system was tuned to mimic a waveform with a large forward flow peak and a lower
101
reverse flow phase, as in Gülan et al. (2012).
102
The non-intrusive, image based measurement technique 3D PTV was applied to measure the flow
103
field. 3D PTV allows to characterize the 3D flow field from both a Eulerian and Lagrangian point of
104
view and it has been used for the analysis of different flow fields [Balducci et al. 2004; Lüthi et al.
105
2005; Holzner et al. 2008; Gülan et al. 2012]. The optical system used for 3D PTV comprises: (1) a 25
106
W continuous laser with a wavelength of 514 nm (BeamLok 2080, Spectra Physics), used to
107
illuminate the investigation domain; (2) a high speed camera (Photron SA5, Tokyo, Japan),
108
synchronized with the pump system to trigger the recordings at the beginning of the pulse wave and
109
allowing to record 1.56 s at full resolution of 1024 x 1024 pixels and 7000 fps; (3) a system of lenses
110
allowing to generate a 3 cm thick laser volume; (4) an orange bandpass filter with a wavelength of 514
111
nm, used to block the optical noise; (5) mirrors to guide the laser path; (6) a four-way image splitter
112
consisting of four fixed primary slanted mirrors assembled on a regular pyramid and four secondary
113
mirrors [Hoyer et al. 2005], with mirrors aligned to obtain a proper stereoscopic vision (Fig. 1, top).
114
An exhaustive description of the setting of the optical system and of its calibration can be found in
115
Gülan et al. (2012).
116
Rhodamine particles (180-200 m diameter range; excitation/emission spectra of 558 - 582 nm) were
117
used as fluorescent tracers [Gülan et al. 2012].
118
The dimension of the investigated domain are 50 x 25 x 25 mm3. The flow velocities in the ascending
119
aorta and aortic arch were measured separately, thereafter the data were combined to produce a single
120
investigation domain covering both regions. A total of 45 cycles (8000 image frames per cycle) were
121
recorded. As we are interested in the large scale aortic fluid structures, only phase averaged results are
122
presented in this study. In detail, the Lagrangian-Eulerian transformation of the PTV measurements
123
was performed by mapping the Lagrangian information into a Eulerian grid of 34 x 20 x 25 voxels of
124
2 x 2 x 2 mm3. The average number of particles in a single voxel after the phase averaging process was
125
226 for the ascending aorta and 201 for the aortic arch, corresponding to a particle density of 28 and
126
25 particles/mm3, respectively.
127
The accuracy of the 3D PTV measurements was estimated elsewhere [Gülan et al. 2012]. Briefly, the
128
position accuracy of the particles was in the range of 0.15, 0.15 and 0.3 mm in x, y and z directions,
129
respectively. The velocity uncertainty of the raw PTV data, as obtained from the calibration, was 0.033
130
m/s. An ad hoc filtering strategy [Gülan et al. 2012] reduced the velocity uncertainty to 0.0072 m/s. As
131
for the velocity assigned to a voxel, the final accuracy was estimated to be 4.09 x 10-4 m/s [Gülan et al.
132
2012].
133
The application of a strategy for the spatial registration of experimental and numerical flow fields was
134
dictated by (1) the need to apply measured 3D velocity profile as condition at the inflow boundary of
135
the computational model and (2) for comparison purposes. Since PTV does not allow to identify walls
136
precisely because of the poor near-wall resolution, the zero velocity isosurface extracted from the
137
measured data was considered here as representative of the vessel wall. In this way, a 3D anatomic
138
model of aorta was obtained from PTV data and used only for the purpose of spatial registration.
139
Then, centerlines of the anatomic models obtained from PTV measurements and from MRI were (1)
140
extracted using the open-source Vascular Modeling ToolKit (VMTK, www.vmtk.org [Antiga et al.
141
2008]) and (2) used to realign the two models, considering the bifurcation points of the centerlines of
142
the aorta and of the supra-aortic vessels as anatomic landmarks.
143 144
In Silico Setup
145
We performed the numerical simulation using a commercial code (ANSYS Fluent 13) based on finite
146
volume method to solve the unsteady-state, incompressible Navier–Stokes equations. The 3D
147
anatomical model was discretized into a tetrahedral mesh of cardinality 19.8×106 (average edge size of
148
0.4 mm) using a commercial mesh generation software (ANSYS ICEM-CFD). This discretization was
149
dictated by the peak Reynolds number at the inlet section. Walls were assumed to be rigid, while fluid
150
was assumed to be Newtonian with the same kinematic viscosity as the working fluid of the
151
experimental set up. Here the same computational setting as in recent studies was applied [Gallo et al.
152
2012; Morbiducci et al. 2013]. In particular, a second-order upwinding method was applied to obtain
153
the solution at each time step (backward Euler implicit time integration scheme, with a fixed time
154
increment equal to 4 ms). The measured cardiac rate of 52 bpm was simulated. As for conditions at
155
boundaries, the time dependent 3D velocity profile prescribed as inflow boundary condition was
156
extracted from 3D PTV measurements, according to recent findings [Morbiducci et al. 2013], as
157
shown in Fig. 2. The measured flow rate waveforms were applied as outflow condition at the supra-
158
aortic vessels and a stress free condition was prescribed at the descending aorta [Gallo et al. 2012].
159
Straight flow extensions were added to the outlet faces of the model. The simulation was carried out
160
for a total of four cardiac cycles to dampen initial transients.
161 162
163 164
RESULTS
165
From the 3D-PTV measurements at the ascending aorta, it was found that the flow regime under
166
investigation is characterized by a mean Reynolds number of 406 and by a Womersley number of 11.
167
To perform the comparison between computational and experimental data, we focused on the largest
168
scale of motion and the simulation data was coarse grained to match the spatial resolution of the
169
measurement. The comparison was performed at different locations along the ascending aorta, i.e.,
170
three cross-sections (A, B and C in Figs. 3 to 5) and one longitudinal section (L in Fig. 6).
171
Considering the cross-sections A-C at the ascending aorta, Figs. 3 and 4 show that, in terms of velocity
172
magnitude, there is a very satisfactory agreement between numerical and experimental results all along
173
the whole forward flow phase (time points from T1 to T4). In particular, in the first part of the forward
174
flow phase ( T1 and T2, Fig. 3), the velocity profile is slightly skewed towards the outer wall and,
175
moving downstream, a slower velocity area grows along the inner wall, as clearly identified by both
176
approaches. At time point T2 (Fig. 3), a higher velocity region is more evident in PTV measurements
177
than in CFD. At mid forward flow phase, the high velocity region and the large slow velocity regions
178
developing toward the outer wall are well captured by both techniques (T3, Fig. 3). Also at late
179
forward flow phase (T4, Fig. 3) the velocity magnitude profiles, even if characterized by more
180
irregular patterns, as expected, still show satisfactory agreement. In the reverse flow phase of the cycle
181
poor agreement is observed between numerical and experimental velocity profiles (T5, Fig. 3).
182
To get a better overview of the flow structures in cross-sections A, B and C, secondary flows and
183
through-plane velocity components around peak forward flow (T2) are presented in Fig. 4. Numerical
184
and experimental data show an overall satisfactory agreement both in the structure of secondary flows
185
and in the profile of the through-plane velocity component, although some discrepancies can be
186
noticed in the in-plane vectors orientation.
187
The velocity magnitude profiles at cross-section A, averaged along the forward and reverse flow
188
phases (Fig. 5), confirm the observations reported above: the agreement between numerical and
189
experimental results is excellent during the forward flow phase of the cycle, while it is poorer during
190
the phase of the cycle characterized by reverse flow.
191
Snapshots of the velocity magnitude in the longitudinal section L are presented in Fig. 6. Also in this
192
case, the onset and development of a fast jet-like structure in the distal ascending aorta is successfully
193
replicated in the numerical model, as well as the slow moving region at the inner wall (T1-T4, Fig. 6).
194
Starting from late forward flow (T4), the computed map of the velocity magnitude in the downstream
195
region of section L shows a poor agreement with the experimentally measured counterpart, especially
196
along the reverse flow phase (T5, Fig. 6).
197
A more quantitative comparison between in silico and in vitro results is reported in Table 1. The
198
analysis performed for velocity magnitude values at cross-sections A, B, C at five time points
199
confirms (1) the good (and statistically significant, p