A reduced-order model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition

Journal of Biomechanics 54 (2017) 33–43

Contents lists available at ScienceDirect

Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com

A reduced-order model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition Gary Han Chang a, Clemens M. Schirmer b, Yahya Modarres-Sadeghi a,⇑ a b

Department of Mechanical and Industrial Engineering, University of Massachusetts, Amherst, MA 01003, United States Department of Neurosurgery and Neuroscience Institute, Geisinger Health System, Wilkes-Barre, PA 18711, United States

a r t i c l e

i n f o

Article history: Accepted 21 January 2017

Keywords: Abdominal aortic aneurysms Computational fluid dynamics Reduced-order model Wall shear stress

a b s t r a c t In this paper, we introduce a method to construct a Reduced-Order Model (ROM) to study the physiological flow and the Wall Shear Stress (WSS) conditions in Abdominal Aortic Aneurysms (AAA). We start the process by running a training case using Computational Fluid Dynamics (CFD) simulations with timevarying flow parameters, such that these parameters cover the range of parameters that we would like to consider in our ROM. We use the inflow angle as the variable parameter in the current study. Then we use the snapshot Proper Orthogonal Decomposition (POD) to construct the reduced-order bases, which are subsequently enhanced using a QR-factorization technique to satisfy the relevant fluid boundary conditions. The resulting ROM enables us to study the flow pattern and the WSS distribution over a range of system parameters computationally very efficiently. We have used this method to show how the WSS varies significantly for an AAA with a simplified geometry, over a range of inflow angles usually considered mild in clinical terms. We have validated the ROM results with CFD results. This approach enables comprehensive analysis of the model system across a range of inflow angles and frequencies without the need to re-compute the simulation for small changes. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Computational Fluid Dynamics (CFD) is a useful tool in analyzing patient-specific Abdominal Aortic Aneurysm (AAA) models. It allows to study the hemodynamic conditions which cannot be directly measured and has been utilized in surgical treatment planning (Pekkan et al., 2008; Taylor and Steinman, 2010; Tse et al., 2011). Simulation and modeling of physiological flows remain underutilized in routine clinical practice in part due to computational difficulties and the limitations of the image-based models. In order to save computational cost, aneurysms are often treated as an isolated model and the computational domain may include the aneurysm body and several segments of the connected arteries, while in reality an aneurysm is part of the global and complicated circulatory system and far-field influences can originate both from upstream and downstream (Quarteroni et al., 2014). CFD simulations aiming at realistic results require appropriate initial assumptions and boundary conditions for fluid properties and pressure. In vivo, the exact boundary conditions are difficult if not impossible to obtain right at the site of the computational ⇑ Corresponding author. E-mail address: [email protected] (Y. Modarres-Sadeghi). http://dx.doi.org/10.1016/j.jbiomech.2017.01.035 0021-9290/Ó 2017 Elsevier Ltd. All rights reserved.

domain (Les et al., 2010). Often the analytical solutions or the measurements from another site of the patient are used (Olufsen et al., 2000; Quarteroni et al., 2014). For example, three-dimensional ultrasound imaging is a time-saving and cost-effective alternative to Computed Tomography (CT) imaging and, as opposed to CT imaging, allows to measure blood flow velocities in real-time. The limited field-of-view, however, usually does not capture the full geometry of the aneurysm and its connected arteries, which increases the uncertainties of the assigned boundary conditions (van Disseldorp et al., 2016). The computational cost of the CFD simulation can increase even more when Fluid-Structure Interactions (FSI) and the aneurysm’s material properties are considered (Li and Kleinstreuer, 2007; Scotti et al., 2008), which makes the CFD runs even less practical for clinical usage. This draws our interest to deriving a Reduced-Order Model (ROM) in order to reduce the computational cost while obtaining reliable results. ROMs have successfully been used for systems with high nonlinearity and geometrical non-uniformity to replace the full, nonlinear model with a low degree-of-freedom ROM, thus reducing the overall computational cost (Dowell and Hall, 2001). In our specific problem, this method provides a means to investigate a wide range of different physiological flow conditions without conducting CFD simulations repetitively. We hypothesize

34

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

that this approach will reduce challenges that impede the use of modeling and CFD results in routine clinical practice and constitutes a component of a putative comprehensive patient-based risk prediction model that may aid clinicians in their decisionmaking process. The balloon-shaped geometries of AAAs combined with complicated vortex formations induced by pulsatile flow result in complex flow patterns. Finol and Amon (2001) found that the presence of pulsatile flow increases the magnitude of Wall Shear Stress Gradients (WSSG) significantly when compared with steady flow at the same time-averaged Reynolds number. Also, as asymmetry is introduced to an AAA’s geometry, the maximum magnitude of WSS is increased, while the diastolic flow is dominated by secondary flows (Finol et al., 2003; Li and Kleinstreuer, 2007). Patient-specific FSI models (Rissland et al., 2009; Xenos et al., 2010) have also been used to estimate the rupture potential index (ratio between the transmural stress and the material strength) and the potential location of rupture. Grinberg et al. (2009) studied the transient turbulent flow in carotid artery with stenosis by applying the POD method to direct numerical simulation (DNS) results. They analyzed the time evolution of POD modes to provide indications for the onset of transition from laminar to turbulent flow during the systolic phase. Arzani and Shadden (2012) studied the flow topology of AAAs using imaged-based patient specific models. By identifying the Lagrangian Coherent Structures (LCS), they found that the vortex formations can change dramatically from one patient to another, and it can also be influenced by exercise (Arzani et al., 2013). Furthermore, Arzani et al. (2016) considered the mass transfer processes of AAAs by studying the complex near-wall flow field and WSS patterns. They showed that the LCS of WSS can provide a template of near-wall transport, which can initiate biological events at the aneurysmal wall. In the present paper, we construct an ROM by the Proper Orthogonal Decomposition (POD) method to estimate the flowinduced Wall Shear Stress (WSS) of a simplified AAA with varying inflow angles. The variation of the inflow angle as well as aortic aneurysm can be found in a number of segments of the aortic artery on its course between the aortic root and the descending aortic artery. In this paper we focus on descending aortic artery segment as the angulation of the proximal neck has been known to be a predictive factor for the surgical outcome of Endovascular Abdominal Aortic Aneurysm (EAAA) repair (Sternbergh et al., 2002; Isselbacher, 2005). Vessel asymmetry with large inflow angles is also known to be a significant indicator of an aneurysm’s potential to rupture (Doyle et al., 2009; Li and Kleinstreuer, 2007; Scotti et al., 2005). The regular method to measure the angulation of aneurysm neck involves observer variations and often leads to over- or under-estimation (de Vries, 2012; van Keulen et al., 2010). Limited field-of-view of ultrasonic images can also increase the difficulty in estimating the angulation (van Disseldorp et al., 2016). Therefore, a numerical model which is robust enough to examine the physiological flow conditions over a range of neck angulation is advantageous. We also study the effectiveness of the ROM when modeling varying frequencies in the hemodynamic waveforms.

2. Reduced order modeling by proper orthogonal decomposition 2.1. Snapshot Proper Orthogonal Decomposition (POD) The snapshot POD method is used here to construct the ROM (Berkooz et al., 1993). Consider the incompressible Navier-Stokes equations in the flow domain of an abdominal aneurysm model:

* *

!

r  u ðx ; tÞ ¼ 0 for x 2 Xf ; *

*

*

qf @ u =@t þ u r u ¼ rp þ

ð1Þ 1 2* r u Re

!

for x 2 Xf ;

ð2Þ

where the boundaries are: * *

*

*

*

u ðx ; tÞ ¼ u in ðx ; tÞ for x 2 Cin ; * *

!

ð3Þ

*

r u ðx ; tÞ  n ¼ 0;

*

pðx ; tÞ ¼ pout ðtÞ for x 2 Cout ; *

ð4Þ *

in which t is time, x is a vector containing the spatial coordinates, u is the flow velocity, p is the flow pressure, qf is the fluid density, Xf is the flow domain, Cin is the inlet, Cout is the outlet, Re is the Rey*

nolds number, u in is the time-dependent inlet boundary condition, and pout is the time-dependent outlet boundary condition. The inflow angle, u, determines the ratio between the lateral and axial =uaxial inlet velocity components, namely, u ¼ tan1 ðulateral in in Þ. Assume a collection of discrete snapshots are sampled at a fixed *

*

time interval, u k ðx ; tk Þ; k ¼ 1    N, and the collection of snapshots contains enough information to describe the behavior of the full* *

scale simulation, u ðx ; tÞ. By calculating the time-correlation coefficients between every two snapshots and solving the discretized eigenvalue problem, the eigenvalues of the POD modes, r2j , and the transformation matrix, akj , from the POD modes to the snapshots can be found as (Rapún and Vega, 2010): N X * * ðui ; uk Þakj ¼ r2j aij ;

ð5Þ

k¼1

where i, j and k are ranged from 1 to N, the number of POD modes used in an ROM, and ð; Þ denotes the usual L2 -inner product. The *

original snapshots, uj , can be expressed in terms of the POD modes, *

U j , as: *

uk ¼

N X

*

*

rj akj U j or U j ¼

j¼1

N 1X

rj k¼1

!

akj u k ;

ð6Þ

and the reduced-order solution is found to be:

^¼ u

n X ai ðtÞ  U i ;

ð7Þ

i¼1 !

!

where ai ¼ ðU i ; u Þ, and n  N is the order of the reduced-order solution. Using the snapshot method results in avoiding the requirement for calculating the eigenmodes from the CFD results, which is rather complicated in geometrically-complex problems. Appropriate sampling frequency, location and number of snapshots depend on the system’s overall behavior, and should be chosen prudently to ensure a converged and physically meaningful result. The magnitude of the eigenvalues of the unused POD modes can be used a priori as an estimator of the Root Mean Square (RMS) error in the reconstructed velocities (Rapún and Vega, 2010): N X j¼1

*

^  uj j ¼ ju

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XN N

j¼nþ1

r2j :

ð8Þ

This relation is critical for deciding on the accuracy of the resulting POD model and how the snapshots are chosen from the simulation results, as we will discuss later in the text corresponding to Fig. 3.

35

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

2.2. POD ROMs with multiple time dependent boundary conditions In our method, we should be able to adapt multiple timedependent boundary conditions of the physiological flow at the inlets and outlets of the computational flow domain. We do this by using the QR-decomposition algorithm (Gunzburger et al., 2007). *

* *

P þnH nP nX X *P * *H * ai ðtÞU i ðx Þ þ ai ðtÞU i ðx Þ;

i¼1

ð9Þ

i¼nP þ1

where nH and nP are the number of homogenous and particular modes, respectively, and n ¼ nP þ nH is the total degrees of freedom. The particular modes have non-zero values at the boundaries, and are used to satisfy the prescribed inflow condition. The discretized Navier-Stokes equations then will have the following form (Gunzburger et al., 2007): N X N N X dai X aj ak Bijk þ aj C ij þ Di ¼ dt j¼1 k¼1 j¼1

for i ¼ nP þ 1    n;

ð10Þ

where

* *  * Bijk ¼  U i ; U j  rU k ; *  * 1 C ij ¼ Re U i ; r2 U j ;

ð11Þ

*  Di ¼ U i ; rpj ; in which hi is the temporal average, Bijk is the convection force tensor, C ij is the viscous force matrix, and Di is the discretized pressure force vector. Eq. (10) is found by applying the Galerkin projection onto the Navier-Stokes equation (Eq. (2)). For example, the i-th discretized Navier-Stokes equation is found by calculating the spatial integration:

Z

N 1X

rj k¼1

akj pk ;

ð13Þ

*

where pk ¼ pk ðx ; tÞ  pout ðtÞ. The resulting pressure POD modes are used to model the fluctuating pressure distribution without additional pressure-velocity equation (Lassila et al., 2014).

*

The eigenmodes consist of particular, U pi , and homogenous. U Hi , modes, and the reconstructed flow velocity is written as:

U ðx ; tÞ ¼

pj ¼

0 * 1  Z *  * * @U j * 1 @ Ui  U i  rp þ r2 U j dV; þ U j  rU k AdV ¼ Re @t V V *

ð12Þ *

where U i is the i-th POD mode. Because the flow is strongly driven by time-dependent inlet pressure, pressure POD modes are directly calculated from the fluctuating pressure values by subtracting the time-dependent outlet pressure from each snapshot,

3. A simplified abdominal aortic aneurysm model with various inflow angles Fig. 1 shows the simplified aneurysm model considered here. The blood flows from proximal to distal in the aorta which then splits into the two iliac arteries. The diameters at the proximal neck, the distal neck, and the iliac arteries are 16 mm, 13.56 mm, and 8 mm, respectively. The maximum diameter of the aneurysm body is 42 mm. The fluid is assumed to exhibit Newtonian behavior and be laminar, with a viscosity of l ¼ 4  103 Pa  s and density of q = 1050 kg/m3. Fig. 2 shows the waveforms of the inlet flow velocity and the outlet pressure boundary conditions with a maximum Reynolds number of Remax = 1700 and a maximum pressure of Pmax = 120 mmHg (16 kPa). The simulation is conducted in OpenFOAM (www.openfoam.org) using the PISO incompressible transient flow solver with a 24-core workstation. We use the Euler scheme in the CFD simulation in OpenFOAM with a time step of 104 sec. 4. The local-POD method To construct an ROM that is capable of predicting the system behavior over a range of different parameters, we first obtain a set of snapshots from a flow simulation with multiple timedependent parameters. We call these snapshots the training set. The training set should be simulated when the parameters of interest are changing, to create the richest dynamical behavior for the snapshots and the resulting POD modes. In the present study, other than the time-dependent inlet flow velocity and the outlet pressure boundary conditions, we varied the incoming flow angle, u, from 30° to +30° over 60 sec in the original simulation, as shown in Fig. 2(c). We chose a range of ±30° to capture the situation encountered in most patients with mild clinical angulation. Cases classified as severely angulated in the literature also often come with significantly asymmetric aneurysm sac geometries and may not be amenable to reasonable generalizations. In this study, we create a total number of 3000 snapshots taken from 60 sec of training CFD simulation using a fixed sampling time interval of 0.02 sec. The eigenvalue problem is solved in the entire region by calculating the time-correlation coefficients between every two snapshots with the spatial average using the cell volume of each node.

Iliac bifurcation

Proximal neck

Inflow angle

Distal neck

Fig. 1. A simplified model for abdominal aortic aneurysm.

36

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

0.5

Inlet flow velocity

(a)

0

0

T

0.5

0.4 m/s

1

1.5

t (s)

2

2.5

3

150

Outlet pressure

(b) 120 mmHg

0

0

0.5

1

1.5

2

2.5

3

t (s) 30

(c) Inflow angle

20 10 0 -10 -20 local-POD 1

local-POD 2

local-POD 3

-30 0

15

30

45

60

t (s) Fig. 2. (a) The inlet flow velocity and (b) the outlet pressure boundary conditions, where T is the period of the blood waveform (in sec), and (c) the Inflow angle over time in the training CFD simulation.

Fig. 3. The error in L2-norm over the mode numbers of the local and complete POD ROM models.

The local POD method together with the Galerkin projection method were first introduced by Rapún and Vega (2010). Projecting the Navier-Stokes equations onto the POD manifolds includes calculating the cubic convection term, and the computational effort increases rapidly with the number of POD modes. Therefore, it is computationally efficient to use local sets of POD modes in multiple intervals in the parameter space, as opposed

to using a whole set of POD modes covering the entire parameter space. It is shown in Fig. 3 that when attempting to construct a POD model using the snapshots from the entire training set, the RMS error in velocities reduces rapidly in the first 100 modes, but the speed of convergence decreases by further increasing the mode number, which results in increasing computational cost for higher

37

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

0

1

(a)

Mode 1

2

3

4

5

6

Mode 1

2

3

4

5

6

Mode 1

2

3

4

5

6

(b)

(c)

Fig. 4. The POD modes of the (a) first, (b) second, and (c) third local-POD models.

modes, and therefore an inefficient ROM. Instead, three local-POD models in time intervals of t e [0, 16]s, [16, 44]s, and [44, 60]s result in even smaller errors, and are computationally far more efficient than the original POD model. Fig. 4(a) shows the first set of local-POD modes that are generated from snapshots within the first interval (t e [0, 16]s). The first set of local-POD modes represents the dominating inflow with a positive inflow angle. Fig. 4(b and c) shows the second (for t e [16, 44]s) and the third (for t e [44, 60]s) sets of the local-POD modes, respectively. The dominating inflow with zero and negative inflow angles can be observed in their first modes, respectively.

5. Flow reconstruction results 5.1. Reconstruction of the flow velocity In this section, the ROMs derived based on local-PODs are used to reconstruct the flow behavior with various inflow angles, and the reconstructed results are compared with the exact CFD results (referred to as the reference results). Fig. 5 shows the comparison between the reconstructed and reference flow velocity distributions for a = 30° and T = 1.0 s, where T is the period of the blood waveform, using 200 POD modes at 4 time instances. It can be observed that the flow velocities are well reconstructed

38

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

(a) t = 0.2

Reconstructed

(b) t = 0.3

Reference

(c) t = 0.4

Reconstructed

0

Reconstructed

0.4 m/s

Reference

(d) t = 0.6

Reference

Reconstructed

Reference

Fig. 5. Reconstructed and reference flow velocities for a = 30° and T = 1 s, at (a) t/T = 0.2, (b) 0.3, (c) 0.4, and (d) 0.6.

for all stages of the systolic cycle in terms of the amplitude and direction of the main flow. This also includes the large vortex that is generated near the aneurysm’s sidewall, between t/ T = 0.4 and 0.6, which is captured in the reconstructed flow field as well. In order to conduct a more quantitative comparison, we extract

structed successfully in a cycle. Fig. 6(b) shows the reconstructed total kinetic energy using 50, 100 and 200 modes, as sample cases. It can be observed that in the present case the reconstructed kinetic energy converges to the value found based on the reference results with relatively low number of modes (100 modes).

the amplitudes of POD modes, aðtÞref j , directly from the CFD (reference) results as:

5.2. Reconstruction of the wall shear stress

 aðtÞref j

¼ *

 u ref ðx ; tÞ; U j ðx Þ ;

*

*

*

*

*

ð14Þ

where u ref ðx ; tÞ is the reference result. Then we compare these amplitudes with the amplitudes of POD modes found in the reconstructed results. Fig. 6(a) shows the comparison between the reconstructed and reference total kinetic energy. At the systolic peak, the kinetic energy is reconstructed well with less than 3% of relative error. The relative error increases in the diastolic phase (up to 22%), but since the flow velocity is much smaller in this phase than in the systolic phase, the majority of the kinetic energy is recon-

We use the ROMs to calculate the WSS distribution. The distribution of WSS is closely related to the near-wall flow pattern and the first N POD modes, which are sorted by the flow’s kinetic energy, are not guaranteed to contain the most dominant information for the WSS reconstruction. Therefore, we replace the last few POD modes by the originally unused POD modes with the highest !

!

time-averaged contribution of WSS (hjai ðtÞ  s i ji, where s i are the WSS distributions in each POD mode). In particular, we replace the last 20 modes of a 200-mode model by the modes with higher WSS contribution, which were originally truncated. The WSS modes are calculated directly using the WSS distributions from the training CFD simulation by using the same transfor-

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

39

Fig. 6. (a) Reconstructed and reference total kinetic energy and their relative error. (b) Reconstructed total kinetic energy with various number of modes in the ROM for a = 30° and T = 1 s.

mation matrix, akj , which transforms the flow velocity distributions to the flow velocity POD modes. The snapshots of WSS are collected at the same fixed time intervals,

*

*

s k ðx ; tk Þ; k ¼ 1    N, and

! ! P ! the WSS modes can be found as Cj ¼ r1j Nk¼1 akj s k , where Cj are the WSS modes. Fig. 7 shows the reconstructed and reference WSS distributions for a = 30° and T = 1.0 s. While the largest magnitude of WSS occurs around the aneurysm’s distal neck and iliac artery bifurcation right at the systolic peak, the large-angle inflow induces a large vortex near the sidewall (t/T = 0.3–0.4) and generates large WSS on the sidewall. The magnitudes of these values are comparable with those previously discussed in the literature. For example, Les et al. (2010) have reported values in the range of 0–1 Pa for time-averaged WSS. The magnitudes obtained in our results are within 0–3 Pa. To quantify the evolution of WSS distribution over the inflow angles, Fig. 8 shows the reconstructed and reference time evolution of the total WSS (WSS integrated over the corresponding area) on the proximal and distal neck (green) and on the sidewall (blue1) with various inflow angles. The reconstructed total WSS values obtained using the original set of POD modes (before updating them following the method discussed in the first paragraph of this section) are shown in black. The total WSS distribution on both the proximal and distal neck area is mainly generated by the main inflow, and is reconstructed accurately by the ROM. The overall time evolution of the total WSS applied

1 For interpretation of color in Fig. 8, the reader is referred to the web version of this article.

on the aneurysm’s sidewall is also well predicted using the ROM, although the maximum magnitude of the WSS is slightly overestimated. From these results, it is clear that while using the updated modes improves the predicted total WSS, it does not have a significant influence on the overall predicted response, which implies that even without updating the modes, only by using several POD modes, one can obtain reliable reconstruction of the total WSS. 5.3. The influence of the inflow angle The fact that the present ROM method is computationally very efficient makes it possible for us to cover a wide range of inflow angles and calculate the WSS as the inflow angle varies. This is of practical significance, because the exact prediction of the inflow angle in a real case is very difficult and therefore it is desired to understand how the WSS varies over a range of inflow angles so that the worst-case scenario can be estimated. Fig. 9 summarizes the influence of inflow angle on the average magnitude of the total WSS on the aneurysmal wall found by ROM. We have run CFD simulations for sample inflow angles to validate the ROM-based results. The continuous line found based on the ROMs shows that the magnitude of the total WSS on the aneurysm’s sac increases with the inflow angle. The present ROM model makes it possible to observe the gradual increase of the total WSS on the aneurysmal sac as the magnitude of the inflow angle changes from 5° to 30°, as well as the sudden increase of the total WSS when the inflow angle deviates from 0°. The rather sudden increase from 0° to 5° is interesting to note from a clinical perspective as almost no patient has a completely straight anatomic configuration of the aorta along the length of the body. This could

40

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

0

(a) t = 0.2

Reconstructed

(b) t = 0.3

Reference

(c) t = 0.4

Reconstructed

3Pa

Reconstructed

Reference

(d) t = 0.6

Reference

Reconstructed

Reference

Fig. 7. The reconstructed and reference WSS for a = 30° and T = 1 s, at (a) t/T = 0.2, (b) 0.3, (c) 0.4, and (d) 0.6.

explain the preponderance of AAAs toward the distal end of the aorta where it is less likely to be subject to an ideal flow regimen. In the results discussed here, the training CFD takes 200 CPU hours. (The time step for the CFD simulation is 0.0001 sec and the maximum cell size is 1.1 mm.) The reconstruction takes only around 1 min for 50 modes, around 5 min for 100 modes, and around 35 min for 200 modes, for a 10 sec reconstruction. The CFD run for each of the reconstructed cases (to make comparison with the POD results) takes 75 CPU hours (with 0.5 million cells). This means that to conduct a parametric study based on say 100 different sets of parameters, we will need 100  75 = 7500 CPU hours for pure CFD results compared with 200 + 0.5  100 = 250 CPU hours for POD results. 5.4. Other waveform frequencies Considering the fact that the frequency of the cardiovascular waveform used in the original training set was set to

1 Hz, the question arises whether or not we can use the ROMs derived based on these simulations and construct the WSS for cases with slightly different frequency. The frequency of 1 Hz, used in the original training simulation is close to the minimum frequency of the cardiovascular waveform, which can be increased to 1.5 Hz (Bowker et al., 2010). Fig. 10 shows the comparison between the reconstructed and reference WSS distributions with three different frequencies of the hemodynamics waveform. In all three cases, the reconstructed results overlay with the reference results. This means that it is possible for the ROMs derived here to predict the WSS for systems with frequencies varying from 1 Hz to 1.5 Hz. This is of significance from an applied point of view, because according to these results, one needs the CFD simulations for one frequency only, and once the POD modes are derived from this set of training simulations, the resulting ROMs can be used to estimate the WSS for systems over a range of frequencies.

41

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

7

x 10

−6

Sidewall WSS, reference Sidewall WSS, reconstructed Sidewall WSS, reconstructed without updated modes

(a) Total WSS (Pa·m2)

6 5

Neck WSS, reference Neck WSS, reconstructed Neck WSS, reconstructed without updated modes

4 3 2 1 0 0 −6 x 10 7

0.5

1

t (s) Time

1.5

2

2.5

3

0.5

1

t (s) Time

1.5

2

2.5

3

0.5

1

t (s)

2

2.5

3

(b) Total WSS (Pa·m2)

6 5 4 3 2 1 0 0

7

x 10

−6

(c) Total WSS (Pa·m2)

6 5 4 3 2 1 0 0

Fig. 8. The reconstructed and reference total WSS on the aneurysm’s sidewall and the neck area for (a) a = 30°, (b) a = 6°, (c) a = 0°, and T = 1 s.

0.95

x10-6

0.9

Total WSS (Pa·m2)

0.85 0.8 0.75 0.7 0.65 0.6

Reference Reconstructed

0.55 0.5 −30

−20

−10

0

10

20

30

Inflow angle Fig. 9. Reconstructed and reference mean WSS on the aneurysmal sac over different inflow angles.

42

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

7

x 10

−6

Sidewall WSS, reference Sidewall WSS, reconstructed Sidewall WSS, reconstructed without updated modes

(a) Total WSS (Pa·m2)

6 5

Neck WSS, reference Neck WSS, reconstructed Neck WSS, reconstructed without updated modes

4 3 2 1 0 7

0 −6 x 10

0.5

1

t (s)

1.5

2

2.5

3

0.5

1

t (s)

1.5

2

2.5

3

0.5

1

t 1.5 (s)

2

2.5

3

(b) Total WSS (Pa·m2)

6 5 4 3 2 1 0

7

0

x 10

−6

(c) Total WSS (Pa·m2)

6 5 4 3 2 1 0

0

Fig. 10. Reconstructed and reference total WSS on the aneurysm’s sidewall and neck area for (a) 1/T = 1.0 (Hz), (b) 1.3 (Hz), (c) 1.5 (Hz), and a = 30°.

6. Conclusions A reduced order model is constructed to study the physiological flow in Abdominal Aortic Aneurysms. The method of snapshot POD is used to construct the reduced-order bases based on a highly robust training CFD simulation, which makes it possible to capture a rich dynamical behavior in the POD modes. Reconstructed distributions of the flow velocity and the wall shear stress are calculated for different physiological flow parameters and are in agreement with the results from the CFD simulations, while maintaining the ROM’s computational efficiency. The main flow pattern, the kinetic energy, and the wall shear stress distribution can be reconstructed from the ROMs very efficiently. The reconstruction of the wall shear stress distribution requires higher POD modes. The ROM efficiently predicts the location and magnitude for the maximum WSS to appear on the aneurysm’s sidewall with various inflow angles and frequencies, corresponding to the anatomic configuration of mild angulation along with a resting heart rate regimen found in the majority of patients. We demonstrate that even in a population where no heightened suspicion of unfavorable flow pattern exists,

a significant deviation from the idealized flow regimen appears to exist which may contribute to significant clinical changes over the long term. The ability for the ROMs to predict the flow velocity and WSS distributions with various parameters is promising for clinical applications. A robust ROM, trained by a well-designed CFD simulation, can be used to investigate the flow conditions in patientspecific image-based aneurysm models with various uncertainties in flow conditions, and is able to capture results across a whole range of values for a given parameter. Even if high-resolution meshes are used for patient-specific aneurysm CFD simulations, the computational cost for a comprehensive parametric study using ROMs will remain low compared with a parametric study conducted purely based on the CFD results. It ought to be mentioned here that the major goal of this paper is to discuss a method for efficient reconstruction of flow pattern inside an AAA, and to this end we have chosen the geometry of the computational domain to be a simple one, so that the focus could be on the POD method and its applicability, rather than CFD methods needed to correctly model the fluid. Also the symme-

G.H. Chang et al. / Journal of Biomechanics 54 (2017) 33–43

try of the system suggests that the training CFD run could have been conducted only for positive (or negative) angles. However, since the long-term goal is to apply this technique to asymmetric systems with more complicated boundary conditions, we decided to focus on a generalized methodology, which will be applicable to the future cases as well. Clearly in patient-specific AAAs, the inflow velocity is far from uniform, and the geometry far from symmetric. The method we have introduced in this paper, however, can be used in the future to model cases with non-uniform incoming flow velocities and asymmetric geometries, where incoming flow profile could be treated as the parameter of interest and POD modes could be created using the method discussed here over a range of flow profiles. Conflict of interest There is no conflict of interest to report. References Arzani, A., Gambaruto, A.M., Chen, G., Shadden, S.C., 2016. Lagrangian wall shear stress structures and near-wall transport in high-Schmidt-number aneurysmal flows. J. Fluid Mech. 790, 158–172. http://dx.doi.org/10.1017/jfm.2016.6. Arzani, A., Les, A.S., Dalman, R.L., Shadden, S.C., 2013. Effect of exercise on patient specific abdominal aortic aneurysm flow topology and mixing. Int. J. Numer. Meth. Biomed. Eng. 30, 280–295. http://dx.doi.org/10.1002/cnm.2601. Arzani, A., Shadden, S.C., 2012. Characterization of the transport topology in patient-specific abdominal aortic aneurysm models. Phys. Fluids 24, 081901. http://dx.doi.org/10.1063/1.4744984. Berkooz, G., Holmes, P., Lumley, J.L., 1993. The proper orthogonal decomposition in the analysis of turbulent flows. Annu. Rev. Fluid Mech. 25, 539–575. http://dx. doi.org/10.1146/annurev.fl.25.010193.002543. Bowker, T.J., Watton, P.N., Summers, P.E., Byrne, J.V., Ventikos, Y., 2010. Rest versus exercise hemodynamics for middle cerebral artery aneurysms: a computational study. Am. J. Neuroradiol. 31, 317–323. http://dx.doi.org/10.3174/ajnr.A1797. de Vries, J.P.P.M., 2012. The proximal neck: the remaining barrier to a complete EVAR world. YSVAS 25, 182–186. http://dx.doi.org/10.1053/j.semvascsurg. 2012.09.001. Dowell, E.H., Hall, K.C., 2001. Modeling of fluid-structure interaction. Annu. Rev. Fluid Mech. 33, 445–490. http://dx.doi.org/10.1146/annurev.fluid.33.1.445. Doyle, B.J., Callanan, A., Burke, P.E., Grace, P.A., Walsh, M.T., Vorp, D.A., McGloughlin, T.M., 2009. Vessel asymmetry as an additional diagnostic tool in the assessment of abdominal aortic aneurysms. YMVA 49, 443–454. http://dx.doi.org/10.1016/j. jvs.2008.08.064. Finol, E.A., Amon, C.H., 2001. Blood flow in abdominal aortic aneurysms: pulsatile flow hemodynamics. J. Biomech. Eng. 123, 474–484. http://dx.doi.org/10.1115/ 1.1395573. Finol, E.A., Keyhani, K., Amon, C.H., 2003. The effect of asymmetry in abdominal aortic aneurysms under physiologically realistic pulsatile flow conditions. J. Biomech. Eng. 125, 207–217. http://dx.doi.org/10.1115/1.1543991. Grinberg, L., Yakhot, A., Karniadakis, G.E., 2009. Analyzing transient turbulence in a stenosed carotid artery by proper orthogonal decomposition. Ann. Biomed. Eng. 2037, 2200. http://dx.doi.org/10.1007/s10439-009-9769-z. Gunzburger, M.D., Peterson, J.S., Shadid, J.N., 2007. Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data. Comput. Methods Appl. Mech. Eng. 196, 1030–1047. http://dx.doi.org/10.1016/j. cma.2006.08.004.

43

Isselbacher, E.M., 2005. Thoracic and abdominal aortic aneurysms. Circulation 111, 816–828. Lassila, T., Manzoni, A., Quarteroni, A., Rozza, G., 2014. Model order reduction in fluid dynamics: challenges and perspectives. In: Reduced Order Methods for Modeling and Computational Reduction. Springer International Publishing, Cham, pp. 235–273. http://dx.doi.org/10.1007/978-3-319-02090-7_9. Les, A.S., Shadden, S.C., Figueroa, C.A., Park, J.M., Tedesco, M.M., Herfkens, R.J., Dalman, R.L., Taylor, C.A., 2010. Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics. Ann. Biomed. Eng. 38, 1288–1313. http://dx. doi.org/10.1007/s10439-010-9949-x. Li, Z., Kleinstreuer, C., 2007. A comparison between different asymmetric abdominal aortic aneurysm morphologies employing computational fluid–structure interaction analysis. Eur. J. Mech. B/Fluids 26, 615–631. http://dx.doi.org/ 10.1016/j.euromechflu.2007.03.003. Olufsen, M.S., Peskin, C.S., Kim, W.Y., Pedersen, E.M., Nadim, A., Larsen, J., 2000. Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Ann. Biomed. Eng. 28, 1281–1299. http://dx.doi.org/10.1114/1.1326031. Pekkan, K., Whited, B., Kanter, K., Sharma, S., de Zelicourt, D., Sundareswaran, K., Frakes, D., Rossignac, J., Yoganathan, A.P., 2008. Patient-specific surgical planning and hemodynamic computational fluid dynamics optimization through free-form haptic anatomy editing tool (SURGEM). Med. Biol. Eng. Comput. 46, 1139–1152. http://dx.doi.org/10.1007/s11517-008-0377-0. Quarteroni, A., Tuveri, M., Veneziani, A., 2014. Computational vascular fluid dynamics: problems, models and methods. Comput. Visual Sci. 2, 163–197. http://dx.doi.org/10.1007/s007910050039. Rapún, M.-L., Vega, J.M., 2010. Reduced order models based on local POD plus Galerkin projection. J. Comput. Phys. 229, 3046–3063. http://dx.doi.org/ 10.1016/j.jcp.2009.12.029. Rissland, P., Alemu, Y., Einav, S., Ricotta, J., Bluestein, D., 2009. Abdominal aortic aneurysm risk of rupture: patient-specific FSI simulations using anisotropic model. J. Biomech. Eng. 131, 031001. http://dx.doi.org/10.1115/1.3005200. Scotti, C.M., Jimenez, J., Muluk, S.C., Finol, E.A., 2008. Wall stress and flow dynamics in abdominal aortic aneurysms: finite element analysis vs. fluid–structure interaction. Comput. Methods Biomech. Biomed. Eng. 11, 301–322. http://dx. doi.org/10.1080/10255840701827412. 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. Online 4. http://dx.doi.org/10.1186/1475-925X-4-64. 64–22. Sternbergh III, W.C., Carter, G., York, J.W., Yoselevitz, M., Money, S.R., 2002. Aortic neck angulation predicts adverse outcome with endovascular abdominal aortic aneurysm repair. J. Vasc. Surg. 35, 482–486. http://dx.doi.org/10.1067/ mva.2002.119506. Taylor, C.A., Steinman, D.A., 2010. Image-based modeling of blood flow and vessel wall dynamics: applications, methods and future directions. Ann. Biomed. Eng. 38, 1188–1203. http://dx.doi.org/10.1007/s10439-010-9901-0. Tse, K.M., Chiu, P., Lee, H.P., Ho, P., 2011. Investigation of hemodynamics in the development of dissecting aneurysm within patient-specific dissecting aneurismal aortas using computational fluid dynamics (CFD) simulations. J. Biomech. 44, 827–836. http://dx.doi.org/10.1016/j.jbiomech.2010.12.014. van Disseldorp, E.M.J., Hobelman, K.H., Petterson, N.J., van de Vosse, F.N., van Sambeek, M.R.H.M., Lopata, R.G.P., 2016. Influence of limited field-of-view on wall stress analysis in abdominal aortic aneurysms. J. Biomech. 1–8. http://dx. doi.org/10.1016/j.jbiomech.2016.01.020. van Keulen, J.W., Moll, F.L., Tolenaar, J.L., Verhagen, H.J.M., van Herwaarden, J.A., 2010. Validation of a new standardized method to measure proximal aneurysm neck angulation. YMVA 51, 821–828. http://dx.doi.org/10.1016/j.jvs.2009.10. 114. 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. Ann. Biomed. Eng. 38, 3323–3337. http://dx.doi.org/10.1007/s10439-010-0094-3.