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

Analysis of thoracic aorta hemodynamics using 3D particle tracking velocimetry and computational fluid dynamics.

Parallel to the massive use of image-based computational hemodynamics to study the complex flow establishing in the human aorta, the need for suitable...
741KB Sizes 0 Downloads 4 Views