Boundary conditions for flow simulations of abdominal aortic aneurysms

International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Boundary conditions for flow simulations of abdominal aortic aneurysms Tamás István Józsa ⇑, György Paál } egyetem rkp. 3., 1111 Budapest, Hungary Department of Hydrodynamic Systems, Faculty of Mechanical Engineering, Budapest University of Technology and Economics, Mu

a r t i c l e

i n f o

Article history: Received 22 August 2013 Received in revised form 10 June 2014 Accepted 11 September 2014 Available online xxxx Keywords: Aortic aneurysm Fluid–structure interaction Boundary conditions Hemodynamic simulations

a b s t r a c t Boundary conditions for abdominal aortic aneurysm simulations are problematic both on the fluid and the solid side. In this paper improvements are suggested on existing methodology in both respects. First, a derivation of a hyperelastic wall model is given, taking into account the wall stresses at the diastolic instant. It is demonstrated that this model can be approximated with a simplified linear wall model in the physiologically interesting range. Then a new method for inlet and outlet boundary condition generation is introduced on the fluid side, based on a one-dimensional transient simulator. Finally, the effect of spine support on the intra-aneurysmal flow is studied. Good agreement was found between rigid wall flow simulations on the ‘‘systolic’’ geometry and the fluid–structure interaction simulations. Other authors found much larger differences because earlier the ‘‘diastolic’’ geometry had been used for comparisons and the stresses in the diastolic state were neglected. It was concluded that the spine support does not have a great impact on the flow field. Significant differences were found between the flow behaviour of artificially generated and real aneurysm geometries. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Cardiovascular diseases (including aneurysms) are the leading cause of death before the age 65 in Europe (Allender et al., 2008). An aneurysm is a localised blood-filled bulge on the wall of a blood vessel. The abdominal aorta is one of the most frequent location of aneurysms. The rupture of abdominal aortic aneurysms (AAAs) leads to death in almost 90% of the cases (Egelhoff et al., 1999). According to the traditional view, when the aneurysm sac reaches a certain size, there is an increasing risk of rupture. A long term goal of aneurysm research is to facilitate the estimation of the risk of rupture for medical doctors more objectively, based on statistical data and simulations, rather than on maximum transverse diameter (Gerouslakos and Nicolaides, 1992; Long et al., 2012). Some authors approach the problem from the wall mechanics point of view because rupture occurs when the stress in the wall reaches the failure strength (Raghavan et al., 2000). Others try to find a correlation between the fluid field variables and the degeneration of the vessel wall (Salsac et al., 2006; Sheidaei et al., 2011). The distensibility can also be a useful parameter representing the biological state of the aortic wall (Blaha et al., 2009).

⇑ Corresponding author. Tel.: +36 1 463 2991. E-mail addresses: [email protected] (T.I. Józsa), [email protected] (G. Paál). URL: http://www.hds.bme.hu (G. Paál).

Computational fluid dynamics (CFD) simulations (Sheidaei et al., 2011), and finite element (FE) analyses (Venkatasubramaniam et al., 2011) based on biomedical imaging techniques are a promising way to determine the fluid mechanical forces and mechanical stresses in aneurysms. A correct technique with rigid wall assumption was worked out for the case of intracranial aneurysms (Szikora et al., 2008; Kulcsár et al., 2012). Probably the rigid wall assumption is legitimate for intracranial aneurysms but it is questionable whether it is appropriate for AAAs. Gaillard et al. (2007) and Vavourakis et al. (2011) state, based on experimental and computational evidence, respectively, that the deformation there is too large to be ignored. Owing to the development of information technology and computational algorithms in recent decades it has become feasible to carry out fluid–structure interaction (FSI) simulations on the interaction between blood flow and vessel (Di Martino et al. (2001)). This is a very resource-intensive method but offers a chance to assess the effect of wall motion on the blood flow more accurately. The aim of this paper is fourfold: first to investigate the effect of the wall deformation on the blood flow. The questions asked are: (i) is the flow field significantly different for FSI simulations than for simulations with rigid walls? (ii) If yes, is the change of the geometry or the wall dynamics the major cause for the differences? To this end FSI simulations were compared with rigid wall simulations of the minimum deformed (‘‘diastolic’’) and the maxi-

http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004 0142-727X/Ó 2014 Elsevier Inc. All rights reserved.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

2

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

mum deformed (‘‘systolic’’) geometries. The influence of geometrical features were examined in the healthy aorta (Tse et al., 2013) and in intracranial aneurysms too (Szikora et al., 2008). Second, the effect of spine support on the flow field in aortic aneurysms was investigated. A supported and a free wall FSI simulation were compared in each case. Third, it was examined analytically and numerically whether the widely used Mooney–Rivlin material model can be approximated linearly in the physiological pressure range. This material model plays an important role in the wall stress prediction (Polzer et al., 2013). The easiest way to define the wall material model is the linear approximation. A better approach is the hyperelastic material model, for example Neo-Hookean (1 parameter) (de Putter et al., 2007) or Mooney–Rivlin (multiple parameter) (Raghavan and Vorp, 2000). In a recent study the patient-specific distribution of the material properties was also taken into account (Reeps et al., 2013). The most realistic hyperelastic material models of the aorta can even handle anisotropy (Geest et al., 2006). The parameters of these models were derived from biaxial measurements and some of them take the bidirectional fibre reinforcement (due to collagen fibres) into account. The right implementation of these anisotropic models is complicated in real geometries; furthermore, the necessary angle of the fibres cannot be determined in vivo. The method to take the axial and tangential residual stresses into consideration was worked out only for simplified geometries (Ogden et al., 2003,). We note that the uniaxial and biaxial measurements are implemented on undeformed sample tissues of the wall and the parameters of the fitted curves are used for simulations on deformed aneurysm geometries (Leung et al., 2006; Khanafer et al., 2009; Kelly and O’Rourke, 2012). The CT images visualise a geometry that is between a maximum and minimum deformed shape. The minimum deformed shape belongs to the diastolic trough, and the maximum deformed geometry manifests itself at the systolic peak. A better way to compute wall deformations and stresses would be to track back the undeformed, zero-pressure geometry based on inverse elastostatics (Vavourakis et al., 2011; de Putter et al., 2007; Raghavan et al., 2006; Polzer et al., 2013) and starting from that to use a full material model (e. g. the Mooney–Rivlin model). Fourth, consistent boundary conditions (BCs) are presented for three-dimensional (3D) flow simulations based on a one-dimensional (1D) transient simulator. A problematic part of hemodynamic simulations is the prescription of the BCs. The generally used method is to prescribe the inlet volumetric flow rate and outlet pressure based on experimental data (Xenos et al., 2010). We propose another solution: a 1D in-house code was used to obtain feasible inlet and outlet BC data in the required section of the arterial system. This way it is not necessary to make invasive measurements to obtain BCs for 3D CFD simulations of arteries. Pereira et al. (2013) used a similar method in case of cerebral aneurysms with rigid walls to define the inlet volumetric flow rate based on the work of Reymond et al. (2009). In the current study the 1D simulator solved the governing equations based on the method of characteristics using a viscoelastic wall model, while Reymond et al. (2009) applied the finite difference discretization technique.

2. Methods One artificial and two real aneurysm geometries were used for the 3D simulations. The artificial geometry was the same as the non-axisymmetric model of Egelhoff et al. (1999). In Fig. 1 the sections of the three investigated geometries are shown. The flow direction is from left to right. To define the material model of the vessel, we look at the possibility of approximating the Mooney– Rivlin hyperelastic wall model in the physiological pressure range

(80–120 mmHg = 10.7–16 kPa) with a linear model. A pseudo Young’s modulus was determined from the hyperelastic material model and was used for later calculations. The outer BCs may also play a great role. Fusiform AAAs are often supported by the spine, so that a rigid constraint should be prescribed on that side. The aorta and the spine in a cross section of the real geometry B is displayed in Fig. 2(a). The aorta is often elongated. In these cases the whole outer surface can move freely. The role of the spine is studied in artificial and real aneurysm geometries using FSI simulations. In case of the real geometry the supported elements on the solid mesh surface were chosen and fixed. The arrows in Fig. 2(b) represent the support in a cross-section. In case of the artificial geometry every element on the flat side was fixed. In an incompressible flow the absolute level of the pressure is irrelevant and only the pressure difference matters. Therefore, on the fluid side the zero level of the relative pressure was associated with the diastolic trough (Di Martino et al., 2001). An in-house 1D code, (developed in MATLAB, based on method of characteristics (Bárdossy and Halász, 2011)) was used to get the transient volumetric flow rate history for the inlet and the pressure history for the outlet boundary conditions. Its validation procedure with in vivo measured data is in progress. We propose the 1D simulator to generate BCs on the fluid side not only for AAAs but in an arbitrary part of the arterial system. For the wall deformation, the zero level of the relative pressure has to be the atmospheric pressure. The calculation of the correct deformation curve involves the determination of the unloaded state or the estimation of the residual stresses in the loaded state. Since these methods are computationally expensive, the authors suggest another approach to compute the deformation. 2.1. Solid domain The original vessel wall model was a hyperelastic Mooney–Rivlin model, based on the particular form of the generalised polynomial hyperelastic strain energy function

W ¼ aðI1  3Þ þ bðI1  3Þ2 ;

ð1Þ

and with a wall density of qw ¼ 2000 kg=m3 (Leung et al., 2006; Scotti et al., 2005). In Eq. (1) a ¼ 17:4 N=cm2 , and b ¼ 188:1 N=cm2 . These parameters were derived by Raghavan and Vorp (2000) from uniaxial measurements, and they were used by several authors (Leung et al., 2006; Khanafer et al., 2009) who performed FSI simulations. I1 denotes the first invariant of the left Cauchy-Green deformation tensor. In the hyperelastic material model incompressibility was assumed because of the high Poisson’s ratio (m ¼ 0:49) of the blood vessel wall (Oshima et al., 2006). To seek simplifications in the generally used Mooney–Rivlin model, a deformed simplified model of the aorta was defined as an L ¼ 170 mm long cylindrical tube with Rdia in ¼ 10 mm inner radius (measured by Allison et al. (2008)) belonging to the diastolic pressure (80 mmHg). The zero-pressure wall thickness was s ¼ 1:5 mm (belonging to Rinitial Þ). The two ends of the tube were in fixed, and the tube was loaded with an internal pressure. There exists a steady state analytical solution for tubes with fixed ends, loaded with internal pressure (Ogden, 1997). The thick-walled solution had to be used here (Zendehbudi and Kazemi, 2007), since the wall thickness-diameter ratio is around 6% (Pearson et al., 1994). The analytical calculations were carried out in Maple 13. The pressure as a function of the tangential stretch in steady state was

p¼

Z

kt;out

kt;in

1 @W dkt ; k2t  1 @kt

ð2Þ

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

3

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

(a) Artificial geometry

(b) Real geometry A

(c) Real geometry B

Fig. 1. The meshes on the solid–fluid interfaces and on the outlet surfaces of the used geometries.

(a) CT image about the spine and the aorta

(b) Spine support on the solid mesh

Fig. 2. The spine support in a cross section (real geometry B).

4

x 10 Mooney−Rivlin model Linear model 2 Finite Element Method

p intraluminal [Pa]

2.5

← systole

1.5

(120 mmHg)

diastole →

1

(80 mmHg)

0.5 0 19

19.5

20

20.5

21

D inner [mm] (a) Sketch about the notations

(b) Pressure as the function of the inner diameter

Fig. 3. Drawings for the analytical solution and the linearisation.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

4

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

where kt is the principal tangential stretch. The principal tangential stretch of the inner surface of the tube was written as where Rin is the instantaneous inner radius and kt;in ¼ Rin =Rinitial in Rinitial is the zero-pressure radius. A similar expression was used in for the principal tangential stretch of the outer surface too (kt;out ¼ Rout =Rinitial out ). The Rout outer radius could be expressed in terms of the inner radius because of the incompressible assumption was the only unknown variable. for the wall material. Thus Rinitial in Fig. 3(a) helps in understanding the notations and imagine the sketched deformation. The first invariant of the left Cauchy–Green tensor as a function of the principal stretches is

I1 ¼ k21 þ k22 þ k23 ;

ð3Þ

where k1 ¼ kt and k2 ¼ 1=kt . The cylindrical tube was inflated without any other load and the two ends of the tube were fixed, the axial stretch k3 was equal to one (plane strain case). The partial derivative of the strain energy function from Eqs. (3) and (1) was

@W 12; 000ð627k8t  1225k6t þ 1225k2t  627Þ ¼ : @kt k5t

ð4Þ

After performing the integration in Eq. (2) an expression of the primitive function as a function of the principal tangential stretch kt is obtained (Eq. (5)).

p ¼ 3:762  106 k2t þ

3:588  106 k2t



1:881  106 k4t

 7:176  106 lnðkt Þ: ð5Þ

Using the Newton–Leibniz rule and the primitive function (5) one can get an expression from Eq. (2) for the pressure as a function of the inner surface tangential stretch. Since the diastolic taninitial gential stretch kdia ¼ Rdia corresponds to the diastolic in =Rin t

pressure

pdia ¼ 80 mmHg

and

Rdia in ¼ 10 mm

is

known,

the

 9:5 mm could be unknown zero-pressure inner radius Rinitial in determined iteratively. Substituting Rinitial ¼ 9:5 mm into Eqs. (2) and (5) yields Eq. (6), in thus a relationship between the pressure and the inner radius has been derived. The systolic inner radius Rsys in ¼ 10:2 mm corresponding to psys ¼ 120 mmHg was calculated simply from Eq. (6) where the SI units were used for all quantities.

p ¼ 4:1968  1010  R2in þ   þ

0:0151 R4in

This material model yields a 4% pulsatile cross-sectional area change. Zhang et al. (2007) measured a cross-sectional area change between 2.97% and 5.46% on the aorta, using an ECG-gated CT, and our value falls within this range. It can be seen in Fig. 3(b) that between the usual diastolic and systolic pressure values, and with the above parameters of the Mooney–Rivlin material model the vessel wall behaviour can be well approximated linearly. The solution was linearised so that the systolic-diastolic pressure difference Dp ¼ 40 mmHg caused exactly the deformation of the inner surface between the systolic and diastolic state. The deformation from the systolic and diastolic radius difference was dia DRin ¼ Rsys in  Rin  0:2 mm. The generalised plane strain solution with a linear elastic material could be written as

321:6 R2in

 7:176  106  lnð105:6Rin Þ

 3:762  106 ð8313R2in þ 0:2548Þ

3:588  106 8313R2in þ 0:2548

þ 3:588  106  lnð8313R2in þ 0:2548Þ

1:881  106 ð8313R2in þ 0:2548Þ

2

ð6Þ

: 0

If the calculations were started from Rinitial ¼ 10 mm as the in undeformed inner diameter with the Mooney–Rivlin material 0 model, we would get Rsys in  10:8 mm for the systolic value (120 mmHg pressure load) meaning  300% error in the deformation. This is the case when CT or MR images are treated as initial geometries and pressure load is defined between 80 and 00 120 mmHg. If the load were only 40 mmHg, and Rdia ¼ 10 mm in were taken as deformed diameter but the existing stresses were 00 ignored, the Mooney–Rivlin model would yield Rsys  10:35 mm in would be 75% error in the deformation. This is the case when the stresses, coming from the diastolic pressure, are neglected.

DRin ¼

1þm n

3  2  2 2 Rdia Rdia in out R Dp 7 6Dp  2  2 4 þ ð1  2mÞ   2 5: R dia dia Rout  Rin Rdia out

ð7Þ

The derivation of Eq. (7) can be found in Bowe (2011). In the above expression n is the parameter of the linearised solution, and R is the radial coordinate along the vessel wall. n can be expressed easily from Eq. (7). Substituting the results from the hyperelastic analytical solution, and the inner diameter R ¼ Rdia in yields the linearised solution (shown in Fig. 3(b)). Based on the above linearisation the hyperelastic material model could be changed to a linear material model which behaves similarly to the hyperelastic one from the deformation point of view in the pressure range 80–120 mmHg. For FSI simulations a pseudo Young’s modulus E ¼ n can be defined in the pressure range 0–40 mmHg which means that only the additional deformations and stresses are modelled, in addition to the ‘‘diastolic’’ mechanical state of the blood vessel wall. The deformations are reproduced very well this way, the wall stresses less so but they do not concern us in this paper. Similar linearisations can be performed for more complex hyperelastic material models as well. 2.2. Fluid domain The inlet volumetric flow rate was calculated with an in-house 1D transient simulator (Bárdossy and Halász, 2013), whose detailed description can be found in Bárdossy and Halász (2011). The section of interest of the 1D human artery model is shown in Fig. 4(a), where the ‘‘A04’’ branch corresponds to the abdominal aorta, and ‘‘A37I’’ and ‘‘A37II’’ to the arteria iliaca communis. The outflow through resistance ‘‘R’’ models the smaller arteries, e. g. the arteria mesenterica inferior and the arteria sacralis mediana. The important branch for our 3D simulations is ‘‘A04’’. The inlet volumetric flow rate in the 3D simulations corresponds to the volumetric flow rate in node ‘‘23a’’. The outlet BC in the 3D simulations was ‘‘opening’’ with prescribed time-dependent pressure obtained from the 1D simulations. BCs too close to the region of interest could influence the flow field, so that it is recommended to use an attached tube section after the examined geometry part too. Accordingly, we should not use the pressure curve from node ‘‘20a’’ but rather that in node ‘‘21ae’’ as an outlet BC. For this we need an equivalent outlet tube length that we can get from the pressure loss Dp between ‘‘20a’’ and ‘‘21ae’’. Based on the symmetry of the human body, the left and right arteria iliaca communis were treated as two parallelly connected tubes, and they were replaced by one hydraulically equivalent tube. The outlet tube section was extruded in every case from the end of the examined part, so that the diameter of the incremental tube was given, and the length of the tube section was unknown. The calculation, to determine the equivalent tube length, can be carried out assuming rigid wall, using the unsteady Bernoulli-equation.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

5

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

account with outflows through resistances. It is possible to take into consideration the drained blood in the 3D model by prescribing a negative volumetric mass source in the outlet tube section. The quantity of drained blood per cycle was the same as that in the 1D transient simulator. (If the mentioned small arteries are clogged by an intraluminal thrombus, this source term is not necessary.) The solver for the 3D blood flow was ANSYS CFX 14.0. The inlet BC on the fluid side was in every case volumetric flow rate with a pulsatile, spatially parabolic velocity distribution. A ten diameter long tube section was attached before the examined geometry part, in order to allow the formation of an accurate Womersleyprofile (Ugron, 2011; Ugron and Paál, 2014). For the segmentation itk-SNAP, ParaView, and VMTK were used. Fig. 1 displays the meshes. The block-structured grid was generated with ICEM CFD 14.0. The parameters of the applied grids can be found in Table 1. The equivalent diameter was calculated from the artery cross-section area as

Deqv ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Asurf  4 :

ð8Þ

p

Based on the ECG-gated CT images, the authors chose the least deformed instant for the segmentation. Therefore it was assumed that the segmentation gave the ‘‘diastolic’’ geometry. The ‘‘systolic’’ geometry was obtained from the FSI simulations as the lumen corresponding to the maximum deformation. 2.3. The coupled model (a) The examined part of the 1D systemic circuit model (Ba´rdossy and Hala´sz, 2013)

(b) The volumetric flow rate in “23a”, and the pressure in “21ae” Fig. 4. The 1D simulator model, and the prescribed boundary conditions.

There are some smaller arteries (arteria mesenterica inferior, arteria sacralis mediana), which branch off the abdominal aorta near of the examined part. The 1D simulator takes them into

Based on the computations in Section 2.1 the vessel wall was treated during the coupled simulations as a linear material with a pseudo Young’s modulus E ¼ 1:9 MPa, and with a Poisson’s ratio m ¼ 0:49. The solver for the solid blood vessel deformation was ANSYS Mechanical 14.0. The solid mesh was extruded from the solid–fluid interface with a uniform 1:5 mm thickness (Jeltsch et al., 2009; Molony et al., 2009), so that it was a shell around the examined fluid section. There were two elements in the radial direction. A mesh convergence study was carried out with steady state simulations and it was found that this spatial resolution is satisfactory from the deformation point of view. The intraluminal thrombus was neglected. On the solid side fixed boundaries were prescribed at the two ends of the artery section (Chandra et al., 2013) and a solid–fluid interface was prescribed on the lumen side. The flow was assumed to be laminar, the blood a Newtonian fluid with a dynamic viscosity l ¼ 0:003 Pa s and with a density qb ¼ 1050 kg=m3 . (It was shown by Fournier (2011) that in large vessels the non-Newtonian behaviour of the blood does not play a large role.) A similar computational model was validated by Ugron et al. (2012) by LDA and PIV measurements for intracranial aneurysms. Fig. 4(b) shows the applied BCs. The pressure difference between the systolic peak and diastolic trough was 6510 Pa. The mean value of the Reynolds number ðRe ¼ Din qb v =lÞ in the rigid wall case (Din ¼ 20 mm) was Remean  520, and the maximum value was Remax  2260. The Womersley number pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (Wo ¼ 0:5Din 2pqb =T=l) is 14.8, where T ¼ 1 s is the period time of a heart cycle. These characteristics are in good agreement with

Table 1 Parameters of the applied grids (art. stands for artificial). Geom.

Nodes in the circumferential direction

Elements in the fluid zone

Elements in the solid zone

Max. element thickness at the wall (mm)

In/outlet equivalent diameter (mm)

Art. A B

80 114 100

471,744 738,208 551,680

26,928 19,928 14,904

0.05 0.05 0.05

20/20 21.7/18.6 17.2/16

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

6

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

residual target on the fluid side was 2:5  105 . The convergence target on the solid side was the default value 102 . The CFX solver uses the fourth-order Rhie–Chow method for the pressure–velocity coupling as default. The temporal scheme was the second order backward Euler differencing scheme, and the available high resolution scheme was chosen to discretize the advection term. The imbalances were checked during the run, and the convergence of the velocity, the pressure and the mesh displacement in selected points were also monitored. A mesh convergence study was performed with steady state on the fluid side of the artificial geometry. The applied spatial resolution was found satisfactory. 3. Results and discussion

Fig. 5. The calculated inlet pressures from the 1D and 3D simulations.

literature data (Mills et al., 1970), and with the Womersley number 12–16.1 measured by Stalder et al. (2011) so that our BCs can be considered realistic. Each simulation contained three heart cycles to let the transient effects attenuate. Only the last cycle was evaluated. The time step was 0.02 s in all cases both on the fluid and the solid side. The

(a) WSS history at the location of the negative peak

(c) WSS history at the location of the maximum deformation

The validity of the linear approximation for the wall model has been demonstrated in Section 2.1. It was checked whether the output data of the transient 1D simulator are compatible with those of the 3D simulations. From the 1D inlet volumetric flow rate and pressure outlet history as BCs a 3D FSI simulation was performed on a single cylindrical tube. The calculated pressure histories at the inlets of the 1D and 3D models were compared, and a very good agreement was found (Fig. 5). Then, the effect of the wall motion was investigated. To this end FSI simulations were carried out and they were compared to the

(b) WSS history at the location of the positive peak

(d) Area-averaged WSS history on the fluidsolid interface

Fig. 6. WSS histories of the artificial geometry.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

simulations with the rigid wall ‘‘diastolic’’ and ‘‘systolic’’ geometry. The ‘‘diastolic’’ geometry was the baseline of the FSI simulations, and the ‘‘systolic’’ geometry was obtained from the FSI simulations as the geometry with the maximum deformation. The pressure distributions on the solid–fluid interfaces are almost the same in the rigid wall and the FSI simulations. For the comparison of the three simulations, three typical points were selected in each geometry. These are the locations of the positive and negative WSS peaks in the flow direction (maximum value in space and time), and the location of the maximum deformation. These nodal values were identified on the solid–fluid interface of the FSI simulations, and the same nodes on the rigid wall boundaries were used for the comparison. During most of the cycle the characteristic points were represented by the same nodes. The time histories of all these values are displayed in Figs. 6– 8. In addition, the time history of the area-averaged WSS on the interface is presented. In Fig. 4(b) the inlet volumetric flow rate is  0 m=s after 0.7 s, and the WSS tends to zero so that the diagrams below show only the WSS histories in the first two-thirds of the cycle. Fig. 6 shows the results for the artificial geometry. On the flat side the deformation is unrealistically large, compared with the measurement results of Arko et al. (2007) (see also Table 2). It is clearly visible in the further diagrams of Fig. 6 that the rigid wall simulations are closer to each other than to the FSI simulation. Typical deformations in the real geometries are much smaller than

7

in the artificial case. There it was found that the positive peak of the WSS in the flow direction is on the less degenerated section of the distal part of the wall, while a negative peak, caused by the backflipping vortices, can be found typically somewhere on the bulge surface. Around the point of the maximum deformation larger differences are detected between the rigid wall and the FSI simulations (Figs. 7(c), 8(c)) but there the WSS values are relatively low and these points can be considered as the most extreme ones. On the one hand, Figs. 7(a) and (b) and 8(a) and (b) show that the rigid wall ‘‘systolic’’ geometry results in a maximum relative error of  10% at the characteristic points compared with the FSI. On the other hand, there is more than 20% relative error in the case of the ‘‘diastolic’’ geometry as found by Vavourakis et al. (2011) as well. According to Figs. 7 and 8 there are no noteworthy differences between the FSI and the rigid wall ‘‘systolic’’ simulations in the characteristic points and this is valid for the whole geometry. Based on Figs. 7 and 8 the rigid wall simulations with the ‘‘systolic’’ geometries provide good estimates for the WSS. The question arises, why the FSI simulations do not agree with the ‘‘systolic’’ rigid wall simulations in the case of the artificial geometry, while they agree well in the real geometries. The differences may appear because of the simplified, too regular shape. The lumen surface of a real aneurysm includes more irregularities. Another cause could be the unnatural flatness of the artificial geometry. Last but not least we have to bear in mind that rigid wall simulations neglect the so called Windkessel effect. Since the

(a) WSS history at the location of the negative peak

(b) WSS history at the location of the positive peak

(c) WSS history at the location of the maximum deformation

(d) Area-averaged WSS history on the fluidsolid interface

Fig. 7. WSS histories of the real geometry A.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

8

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

(a) WSS history at the location of the negative peak

(c) WSS history at the location of the maximum deformation

(b) WSS history at the location of the positive peak

(d) Area-averaged WSS history on the fluidsolid interface

Fig. 8. WSS histories of the real geometry B.

Table 2 Maximum deformations and deformation velocities on the solid–fluid interfaces (art. stands for artificial). Geom.

Max. def. (mm)

Max. def. on the spine side (mm)

Area-averaged def. (mm)

Max. def. vel. (m/s)

Art. A B

6.2 3.8 2.9

6.2  0:7  0:8

0.9 0.7 0.5

0.070 0.032 0.021

periodic deformation of the vessel wall modifies the volumetric flow rate, it has clearly an effect on the velocity field, thus on the WSS distribution too. The artificial geometry went through a slightly different deformation history which may lead to a different Windkessel mechanism. The authors could not identify unequivocally the main cause of the difference. Next, the effect of the spine support was investigated. The spine-supported flow field was compared with the free geometry with FSI simulations to look at whether that the spine support has a significant influence on the flow field. Table 2 contains some important data about the deformation of the free geometries. These data correlate well with the FE simulation of Gee et al. (2009) who used the same nonlinear material model with different implementation and the maximum displacement was between 1.35 and 2.5 mm. In that paper the intraluminal thrombus was also modelled, which can be responsible for the smaller deformations.

The in vivo ultrasonography measurements carried out by Arko et al. (2007) result in an average deformation 1:7  0:6 mm on the aneurysm surface, which is also in good agreement with our values. To compare the supported and free simulations several characteristic points were investigated but they display the same trend. Here only the area-averaged histories of the WSS on the solid–fluid interface are presented (Fig. 9). The diagrams in Fig. 9 show that the effect of the support is not negligible in the case of the artificial geometry. In this case there is an unrealistically large deformation on the supported (flat) side, and it is responsible for the differences. In contrast to the artificial geometry, the supported simulations of the real aneurysms are in great agreement with those without the support (Fig. 9(b) and (c)). The current study overlooks biological and physical factors which proved to be important in earlier research. The authors assume that these factors do not play a major role from the actual comparisons point of view but find important to mention some of them. It was attempted to find aneurysms with relatively little thrombus (see Fig. 2(a)), hence modelling of the intraluminal thrombus was omitted such as in earlier studies (Vavourakis et al., 2011; Piccinelli et al., 2013). It had been demonstrated that the varying wall thickness is important in risk of rupture prediction (Reeps et al., 2013; Scotti et al., 2005). Furthermore, it is a fact that the diseased vessel wall has a different mechanical behaviour than

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

(a) Area-averaged WSS history of the artificial

(b) Area-averaged WSS history of the real

geometry

geometry A

9

(c) Area-averaged WSS history of the real geometry B Fig. 9. Area-averaged WSS histories on the solid–fluid interface.

the healthy one (Oshima et al., 2006). All these factors have been neglected in this paper. Moreover, although the authors are aware of the importance of validation experiments to substantiate these findings, they are not aware of any in vivo experiments of such kind. The continuous improvement of time resolved medical imaging techniques could provide at least an indirect validation in the near future. Until then these results can be compared only with those from other computational studies. 4. Conclusions An analytical solution was presented for the simplified model of the aorta with the commonly used Mooney–Rivlin model. It was found that the model behaves nearly linearly in the physiologically relevant pressure range. This way a linear model can be used for investigation of the flow field with FSI simulations, thereby significantly simplifying them. A recent 1D numerical model yields the BCs for our 3D calculations. This method provides reliable inlet and outlet BCs for CFD simulations in an arbitrary section of the arterial system. As far as the necessity of FSI simulations vs. rigid wall simulations was concerned, we found remarkable differences between the behaviour of the artificial and those of the real geometries. In the former case the deformations and the deformation speeds were much higher, although the geometrical dimensions and the material model were the same as for the real geometries. Here also

significant qualitative and quantitative differences have been found between FSI and rigid wall simulations. In the real cases, however, the deformations were lower and the agreement was satisfactory both qualitatively and quantitatively. The better agreement was provided always by the ‘‘systolic’’ geometry, i. e. the geometry frozen at the largest displacement. The results highlight that the geometry change is important rather than the wall dynamics in the differences between rigid wall and FSI simulations. If one wants to perform only ‘‘systolic’’ rigid wall simulations, it is not necessary to carry out an FSI simulation to obtain the systolic geometry; it is much easier to extract it directly from ECG-gated CT images. Other authors found that FSI is necessary because the differences between FSI and rigid wall results were impermissibly large. This can be attributed to two factors: first they did not take into account that the medical images do not correspond to the zeropressure state, so that they may overestimate the deformation by a factor of three. Second, even when the wall model was treated correctly (back tracking), the basis for the rigid wall simulations was the ‘‘diastolic’’ geometry, which, as our computations demonstrate, yields a much worse approximation of the FSI simulation results than the ‘‘systolic’’ geometry. There are some extreme singular points where this agreement is worse but still better than in the case of the artificial geometry. This leads to the conclusion that we should exercise extreme care when extrapolating results from such artificial models of aortic aneurysms to real aneurysms.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004

10

T.I. Józsa, G. Paál / International Journal of Heat and Fluid Flow xxx (2014) xxx–xxx

Similar conclusion is reached for the spine support. We found that the spine support has practically no influence on the flow field for the real geometries while it has a noticeable influence for the artificial geometry. Examining the velocity fields we found similar trends. Note, however, that on the solid side supports of the environmental tissues have great influence on the wall stresses. Acknowledgements We thank Dr. Attila Kossa and Dr. Gergely Bárdossy for useful discussions. We thank Dr. Péter Sótonyi for providing us the CT images of AAAs. This work was financially supported by the New Széchenyi Development Plan (TÁMOP-4.2.1/B-09/1/KMR-20100002). References Allender, S., Scaborough, P., Peto, V., Rayner, M., Leal, J., Luengo-Fernandez, R., Gray, A., 2008. European Cardiovascular Disease Statistic. European Heart Network, Brussels. Allison, M.A., Kwan, K., DiTomasso, D., Wright, C.M., Criqui, M.H., 2008. The epidemiology of abdominal aortic diameter. J. Vasc. Surg. 48, 121–127. Arko, F.R., Murphy, E.H., Davis, C.M., Johnson, E.D., Smith, S.T., Zarins, C.K., 2007. Dynamic geometry and wall thickness of the aortic neck of abdominal aortic aneurysms with intravascular ultrasonography. In: Presented at the Thirty-First Annual Meeting of the Southern Association for Vascular Surgery. Puerto Rico, January 18–21. Bárdossy, G., Halász, G., 2011. Modeling blood flow in the arterial system. Period. Polytech. 55, 49–55. Bárdossy, G., Halász, G., 2013. A ‘‘backward’’ calculation method for the estimation of central aortic pressure wave in a 1D arterial model network. Comp. Fluids 73, 134–144. Blaha, M.J., Budoff, M.J., Rivera, J.J., Katz, R., O’Leary, D.H., Polak, J.F., Takasu, J., Blumenthal, R.S., Nasir, K., 2009. Relationship of carotid distensibility and thoracic aorta calcification: multi-ethnic study of atherosclerosis. Hypertension 54, 1408–1415. Bowe, A.F., 2011. Applied Mechanics of Solids. CRC press, Chapter 4. Chandra, S., Raut, S.S., Jana, A., Biederman, R.W., Doyle, M., Muluk, S.C., Finol, E.A., 2013. Fluid–structure interaction modeling of abdominal aortic aneurysms: the impact of patient-specific inflow conditions and fluid/solid coupling. J. Biomech. Eng., 135. Di Martino, E.S., Guadagni, G., Fumero, A., Ballerini, G., Spirito, R., Biglioli, P., Redaelli, A., 2001. Fluid–structure interaction within realistic threedimensional models of the aneurysmatic aorta as a guidance to assess the risk of rupture of the aneurysm. Ann. Biomed. Eng. 23, 647–655. Egelhoff, C.J., Budwig, R.S., Elger, D.F., Khraishi, T.A., Johansen, K.H., 1999. Model studies of the flow in abdominal aortic aneurysms during resting and exercise conditions. J. Biomech. 32, 1319–1329. Fournier, R.L., 2011. Basic Transport Phenomena in Biomedical Engineering. CRC Press, United States of America. Gaillard, E., Bergeron, P., Deplano, V., 2007. Influence of wall compliance on hemodynamics in models of abdominal aortic aneurysm. J. Endovasc. Therapy 14, 593–599. Gee, J.M.W., Reeps, C., Eckstein, H.H., Wall, W.A., 2009. Prestressing in finite deformation abdominal aortic aneurysm simulation. J. Biomech. 42, 1732–1739. Geest, J.P.V., Sacks, M.S., Vorp, D.A., 2006. The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta. J. Biomech. 39, 1324–1334. Gerouslakos, G., Nicolaides, A., 1992. Infrarenal abdominal aortic aneurysms less than five centimetres in diameter: the surgeon’s dilemma. Euro. J. Vasc. Surg. 6, 616–622. Jeltsch, M., Klass, O., Klein, S., Feuerlein, S., Aschoff, A.J., Brambs, H.J., Hoffmann, M.H., 2009. Aortic wall thickness assessed by multidetector computed tomography as a predictor of coronary atherosclerosis. Int. J. Cardiovasc. Imag. 25, 209–217. Kelly, S., O’Rourke, M., 2012. Fluid, solid, and fluid–structure interaction simulations on patient-based abdominal aortic aneurysm models. J. Eng. Med. 226, 288–304. Khanafer, K.M., Bull, J.L., Berguer, R., 2009. Fluid–structure interaction of turbulent pulsatile flow within a flexible wall axisymmetric aortic aneurysm model. Euro. J. Mech. B/Fluid 28, 88–102. Kulcsár, Z., Ugron, Á., Marosf’i, M., Berentei, Z., Paál, G., Szikora, I., 2012. Hemodynamics of cerebral aneurysm initiation: the role of wall shear stress and spatial wall shear stress gradient. Am. J. Neuroradiol. 32 (3), 587–594. Leung, J.H., Wright, A.R., Cheshire, N., Crane, J., Thom, S.A., Hughes, A.D., Yun, Xu, 2006. Fluid structure interaction of patient specific abdominal aortic aneurysms: a comparison with solid stress models. BioMed. Eng. 5, OnLine. Long, A., Rouet, L., Lindholt, J.S., Allaire, E., 2012. Measuring the maximum diameter of native abdominal aortic aneurysms: review and critical analysis. Eur. J. Vasc. Endovasc. Surg. 43, 515–524. Mills, C., Gabe, I., Gault, J., Mason, D., Ross, J., Braumwald, E., Shillingford, J., 1970. Pressure–flow relationships and vascular impedance in man. Cardiovasc. Res. 4, 4:405–417.

Molony, D.S., Callanan, A., Kavanagh, E.G., Walsh, M.T., McGloughlin, T.M., 2009. Fluid–structure interaction of a patient-specific abdominal aortic aneurysm treated with an endovascular stent-graft. BioMed. Eng. 8, OnLine. Ogden, R.W., 1997. Non-Linear Elastic Deformations. Dover Publications, United Kingdom. Ogden, R.W., . In: Holzapfel, G.A., Ogden, R.W. (Eds.), Biomechanics of Soft Tissue in Cardiovasular Systems, CISM Courses and Lectures Series, vol. 441. Springer, Austria, pp. 65–108. Pearson, A.C., Guo, R., Orsinelli, D.A., Binkley, P.F., Pasierski, T.J., 1994. Transesophageal echocardiographic assessment of the effects of age, gender, and hypertension on thoracic aortic wall size, thickness, and stiffness. Am. Heart J. 128, 344–351. Pereira, V.M., Brina, O., Gonzales, A.M., Narata, A.P., Bijlenga, P., Schaller, K., Lovblad, K.O., Ouared, R., 2013. Evaluation of the influence of inlet boundary conditions on computational fluid dynamics for intracranial aneurysms: a virtual experiment. J. Biomech. 46, 1531–1539. de Putter, S., Wolters, B.J., Rutten, M.C., Breeuwer, M., Gerritsen, F.A., van de Vosse, F.N., 2007. Patient-specific initial wall stress in abdominal aortic aneurysms with a backward incremental method. J. Biomech. 40, 1081–1090. Raghavan, M.L., Vorp, D.A., Federle, M.P., Makaroun, M.S., Webster, M.W., 2000. Wall stress distribution on three dimensionally reconstructed models of human abdominal aortic aneurysm. J. Vasc. Surg. 31, 760–769. Raghavan, M.L., Vorp, D.A., 2000. Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysms: identification of a finite strain constitutive model and evaluation of its application. J. Biomech. 33, 475–482. Raghavan, M.L., Ma, B., Fillinger, M.F., 2006. Non-Invasive determination of zeropressure geometry of arterial aneurysms. Ann. Biomed. Eng. 34, 1414–1419. Reeps, C., Maier, A., Pelisk, J., Hrtl, F., Grabher-Meier, V., Wall, W.A., Essler, M., Eckstein, H.H., Gee, M.W., 2013. Measuring and modeling patient-specific distributions of material properties in abdominal aortic aneurysm wall. Biomech. Model. Mechanobiol. 12, 717–733. Reymond, P., Merenda, F., Perren, F., Rfenacht, D., Stergiopulos, N., 2009. Validation of a one-dimensional model of the systemic arterial tree. Am. J. Physiol. – Heart Circul. Physiol. 297, 208–222. Salsac, A.V., Sparks, S.R., Chomaz, J.M., Lasheras, J.C., 2006. Evolution of the wall shear stresses during the progressive enlargement of symmetric abdominal aortic aneurysms. J. Fluid Mech. 560, 19–51. Scotti, C.M., Shkolnik, A.D., Muluk, S.C., Finol, E.A., 2005. Fluid–structure interaction in abdominal aortic aneurysms: effects of asymmetry and wall thickness. BioMed. Eng. 4, OnLine. Sheidaei, A., Hunley, S.C., Zeinali-Davarani, S., Raguin, L.G., Baek, S., 2011. Simulation of abdominal aortic aneurysm growth with updating hemodynamic loads using a realistic geometry. Med. Eng. Phys. 30, 80–88. Stalder, A.F., Frydrychowicz, A., Russe, M.F., Korvink, J.G., Hennig, J., Li, K., Markl, M., 2011. Assessment of flow instabilities in the healthy aorta using flow-sensitive MRI. J. Magn. Reson. Imaging 33, 839–846. Polzer, S., Gasser, T.C., Bursa, J., Staffa, R., Vlachovsky, R., Man, V., Skacel, P., 2013. Importance of material model in wall stress prediction in abdominal aortic aneurysms. Med. Eng. Phys. 35, 1282–1289. }i, M., Berentei, Zs., Kulcsár, Szikora, I., Paál, G., Ugron, Á., Nasztanovics, F., Marosfo Zs., Lee, W., Bojtár, I., Nyáry, I., 2008. Impact of aneurysmal geometry on intraaneurysmal flow: a computerized flow simulation study. Neuroradiology 50, 411–421. Oshima, M., Torii, R., Takagi, T., 2006. In: Holzapfel, G.A., Ogden, R.W. (Eds.), Mechanics of Biological Tissue. Springer, Germany, pp. 323–335. Piccinelli, M., Vergara, C., Antiga, L., Forzenigo, L., Biondetti, P., Domanin, M., 2013. Impact of hemodynamics on lumen boundary displacements in abdominal aortic aneurysms by means of dynamic computed tomography and computational fluid dynamics. Biomech. Model. Mechanobiol. 12 (6). Tse, K.M., Chang, R., Lee, H.P., Lim, S.P., Venkatesh, S.K., Ho, P., 2013. A computational fluid dynamics study on geometrical influence of the aorta on haemodynamics. Eur. J. Cardiothorac. Surg. 43, 829–838. Ugron, Á., 2011. Flow Simulation in Intracranial Aneurysms, Ph.D. thesis. Budapest University of Technology and Economics, Budapest. Ugron, Á., Farinas, M.I., Kiss, L., Paál, G., 2012. Unsteady velocity measurements in a realistic intracranial aneurysm model. Exp. Fluids 52 (1), 37–52. Ugron, Á., Paál, G., 2014. On the boundary conditions of cerebral aneurysm simulations. Period. Polytech. Mech. Eng. 58 (1), 37–45. Vavourakis, V., Papaharilaou, Y., Ekaterinaris, J.A., 2011. Coupled fluid–structure interaction hemodynamics in a zero-pressure state corrected arterial geometry. J. Biomech. 44, 2453–2460. Venkatasubramaniam, A.K., Fagan, M.J., Mehta, T., Mylankal, K.J., Ray, B., Kuhan, G., Chetter, I.C., McCollum, P.T., 2011. A comparative study of aortic wall stress using finite element analysis for ruptured and non-ruptured abdominal aortic aneurysms. Euro. J. Vasc. Endovasc. Surg. 28, 168–176. Xenos, M., Rambhia, S.H., Alemu, Y., Einav, S., Labropoulos, N., Tassiopoulos, A., Ricotta, J.J., Bluestein, D., 2010. Patient-based abdominal aortic aneurysm rupture risk prediction with fluid structure interaction modeling. Med. Eng. Phys. 23, 647–655. Zendehbudi, G.R., Kazemi, A., 2007. The accuracy of thin-shell theory in estimation of aneurysm rupture. J. Biomech. 40, 3230–3235. Zhang, J., Fletcher, J.G., Vrtiska, T.J., Manduca, A., Thompson, J.L., Raghavan, M.L., Wentz, R.J., McCollough, C.H., 2007. Large-vessel distensibility measurement with electrocardiographically gated multidetector CT: phantom study and initial experience. Radiology 245, 258–266.

Please cite this article in press as: Józsa, T.I., Paál, G. Boundary conditions for flow simulations of abdominal aortic aneurysms. Int. J. Heat Fluid Flow (2014), http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.09.004