• No results found

Pulsatile aortic blood flow – A critical assessment of boundary conditions

N/A
N/A
Protected

Academic year: 2022

Share "Pulsatile aortic blood flow – A critical assessment of boundary conditions"

Copied!
55
0
0

Loading.... (view fulltext now)

Full text

(1)

Postprint

This is the accepted version of a paper published in . This paper has been peer-reviewed but does not include the final publisher proof-corrections or journal pagination.

Citation for the original published paper (version of record):

Fuchs, A., Berg, N., Prahl Wittberg, L. [Year unknown!]

Pulsatile aortic blood flow – A critical assessment of boundary conditions

ASME Journal of Engineering and Science in Medical Diagnostics and Therapy (JESMDT)

Access to the published version may require subscription.

N.B. When citing this work, cite the original published paper.

Permanent link to this version:

http://urn.kb.se/resolve?urn=urn:nbn:se:kth:diva-278873

(2)

1 JESMDT-20-1007, Fuchs

Pulsatile aortic blood flow – A critical assessment of boundary

conditions

Alexander Fuchs1

Department of Radiology, Karolinska University Hospital, Stockholm.

Department of Medical and Health Sciences, Linköping University, Linköping, Sweden.

FLOW & BioMEx, Department of Engineering Mechanics, Royal Institute of technology, KTH,

Osquars Backe 18, SE-100 44 Stockholm, Sweden.

e-mail: alex@mech.kth.se Niclas Berg

FLOW & BioMEx, Department of Engineering Mechanics, Royal Institute of technology, KTH,

Osquars Backe 18, SE-100 44 Stockholm, Sweden.

e-mail: niber@mech.kth.se Lisa Prahl Wittberg

FLOW & BioMEx, Department of Engineering Mechanics, Royal Institute of technology, KTH,

Osquars Backe 18, SE-100 44 Stockholm, Sweden.

e-mail: prahl@mech.kth.se

1 Corresponding author

(3)

2 JESMDT-20-1007, Fuchs ABSTRACT

Patient specific (PS) blood flow studies have become popular in recent years but have thus far had limited clinical impact. This is possibly due to uncertainties and errors in the underlying models and simulations set-up. This study focuses on the sensitivity of simulation results due to in- and outflow boundary conditions (BC:s). Nine different inlet- and seven different outlet BC:s were applied to two variants of a healthy subject’s thoracic aorta. Temporal development of the flow is essential for the formation and development of helical/spiralling flow where the commonly observed clockwise helical motion may change direction during the heart-cycle. The sensitivity to temporal and spatial variations in the inlet conditions is significant both when expressed in terms of mean and maximal wall shear stress (WSS) and its different averaged variables, e.g. Time-Averaged WSS (TAWSS), Oscillating Shear Index (OSI) and Relative Residence Time (RRT). The simulation results are highly sensitive to BC. For example, the maximal WSS may vary over 3 orders of magnitude (1 to 1000 Pa) depending on particular combinations of BC:s. Moreover, certain formulations of outlet boundary conditions may be inconsistent with the computed flow field if the underlying assumptions of the space-time dependence are violated. The results of this study show that CFD simulations can reveal flow details that can enhance understanding of blood flows. However, the results also demonstrate the potential difficulties in mimicking blood flow in clinical situations.

INTRODUCTION

Patient-specific (PS) flow simulations using computational fluid dynamics (CFD) has evolved successively over recent years, aiming at improving patient care. A large number of scientific publications describe the application of CFD to specific clinical cases. PS simulations has by some

(4)

3 JESMDT-20-1007, Fuchs been defined as simulations using measured values of the arteries under consideration along with measured data at each boundary. Another approach, adopted in this work, is to use the term PS only to characterize the geometry of the vessel but the specified boundary conditions are determined from other, non-patient specific data and/or assumptions. The underlying argument to strive for patient-specific studies is the potential of gaining more efficient diagnostic tools and treatments that in turn would reduce patient anxiety and societal costs. The advantage of CFD is capturing details of blood flow having a temporal and spatial resolution that cannot be matched by diagnostic tools.

Studies using vessel geometries derived from medical images to study blood flow associated with certain groups of pathologies (e.g. aortic aneurysms and arterial stenosis) are commonly found in literature. In particular, wall shear stress (WSS) and related expressions have been used to assess the effects of blood flow on the vessel walls and the potential risk associated with arterial wall

pathologies.

Although its advantages, there are several uncertainties and potential sources of errors using CFD, one being setting appropriate boundary conditions (BC). In blood flow simulations, a chain of assumptions has to be made in order to make the problem tractable. Heart- and flow rates may vary by a factor of more than three between rest and exercise/stress [1]. The temporal contraction rate of the left ventricle may also vary among the heart beats themselves, leading to natural cycle to cycle variation that depends on for example breathing or psychological factors. Other, although smaller, natural variations are associated with blood compositions and in the mechanical properties of the arterial walls. In addition to water, the blood contains components of different sizes ranging from

(5)

4 JESMDT-20-1007, Fuchs nm (ions and small molecules) to tens of m (cells). The red blood cells (RBC) have high concentration in the plasma where the mean distance among individual cells is of the order of the cells themselves.

Often, a continuum model is used for the plasma and the RBC phases with constant viscosity for the former and variable for the latter. Other blood constituents of different sizes are of much lower concentration and the interaction between similar substances is weak. Thus, the other blood components are often treated as individual particles [2]. The main drawback using a continuum approach is that the concentration and mechanical properties of the RBC vary from one person to another. Furthermore, these quantities also vary for each individual over time and depend on individual conditions at the given instant of time. Another common approach is assuming that the blood is a perfect mixture where the apparent viscosity of the mixture satisfies some empirical correlation based on the local shear rate and RBC volume fraction (hematocrit). Moreover, assuming that blood is an incompressible fluid and that the local shear is larger than some threshold, the RBC clustering (rouleaux) can be disabled allowing the mixture to be assumed to have a constant viscosity, approximately 3-4 times that of water. Besides the fluid model used, another major point of debate in blood flow simulations is the handling of the arterial walls. Similarly to the above, the mechanical properties of the arterial wall vary along the arterial tree and over time (e.g. aging) where

characteristics of arterial wall compliance is an active field of research [3-5]. However, since aging and pathologies of the arteries (atherosclerosis) reduce the distensibility of the arteries, it is common to consider arterial walls as rigid. Another source of uncertainty coupled to the vessel walls is

encountered when deriving the patient specific vessel geometry. The different imaging techniques applied, such as computed tomography (CT) or magnetic resonance imaging (MRI), all have different

(6)

5 JESMDT-20-1007, Fuchs limits in spatial and temporal resolution. Consequently, the segmentation procedure needed for generating a computerized shape of the vessel geometry may introduce further approximations and errors.

The flow in larger arteries has commonly been assumed to be laminar or transitional, unless the patient has severe stenosis or obstruction leading to the formation of a high speed jet [6]. Such jets form a shear layer with fast growth rate of disturbances, enabling the formation of turbulence that may persist over a considerable distance along the aorta, an effect that can be heard using a common stethoscope. Without the presence of a significant stenosis, the flow in a normal aorta is laminar with possibly intermittently isolated transitional episodes [6-7]. However, recent 4D-MRI based aortic flow measurements in healthy individuals have reported significant levels of turbulence [8 - 12]. If applying adequate temporal and spatial resolution, it is possible to capture transitional and intermittently turbulent flows numerically without the need for additional (turbulence) modeling.

As mentioned above, boundary as well as initial conditions are a crucial part of performing CFD simulations. In principle, blood flow in large arteries can be measured by MRI (e.g. phase contrast, PC-MRI, as a common method) and ultrasound (US) with Doppler, providing flow details that in turn can be used to derive BC. Recent MRI measurements of the time-dependent flow near the aortic root have been performed [13 - 17]. However, besides being restricted to reflecting the blood flow conditions at the time of measurement, measured data is associated with uncertainties and measurement errors. Assessment of the result’s sensitivity and in particular the parameters used

(7)

6 JESMDT-20-1007, Fuchs for evaluating results to uncertainties and errors is crucial if simulation data is to be used for quantitative evaluation.

Different assessments of the effects of boundary data on certain parameters are found in literature [14, 18 - 21]. The issue of outflow BC has attracted considerably more attention as compared to inflow BC [18, 21 -24]. Inlet conditions to the aorta can be measured by the abovementioned techniques and is therefore commonly assumed to be “known”, or at least

measurable. On the other hand, it is practically impossible to measure the spatial and temporal inlet velocity with an accuracy in par with the resolution needed for high fidelity blood flow simulation.

Sankaran et al. [25] proposed a data-driven framework for modeling uncertainties and studied the impact on blood flow in stenoted pipe flow, using incompressible Navier–Stokes equations with an adaptive stochastic collocation method. The results showed that uncertainties associated with computed results are not additive but are only slightly greater than the highest single parameter uncertainty when put together. Bozzi et al. [26] studied propagation of inflow boundary condition uncertainties using the steady-state incompressible Navier-Stokes equations. A main conclusion was that uncertainty in inlet boundary condition implies larger uncertainty when computing the pressure and the WSS. Several studies focused on investigating the effects of using PS inlet velocity profiles for CFD simulated blood flows in different arterial segments/branches [27] and carotid arteries [28].

Morbiducci et al. [29] considered the effect of outlet BC and in later work [30], the effect of inlet BC:s on the thoracic aorta of a healthy subject, reporting that that idealized BC:s may lead to misleading results.

(8)

7 JESMDT-20-1007, Fuchs Several studies have carried out aortic flow simulations using PS inlet flow data [14 - 23].

Benim et al. [19] investigated the flow in an aorta including the main aortic branches assuming solid arterial walls and constant blood viscosity, applying the SST turbulence model. The time dependent flow-rate was specified at the inlet (aortic root) and a certain fraction of the inlet flow was set on each aortic branch. As an alternative BC, Benim et al. [19] also applied a constant pressure, showing that there are significant differences in the flow distribution among the different branches when compared to clinically observed data. Madhavan and Kemmerling [21] used similar types of outlet conditions, although characterized by somewhat different flow-rate among distributions among the major aortic branches. As an alternative, the effects of using two- and three-element Windkessel BC were investigated. At the inlet, the time dependent flow-rate was specified, but with different spatial distributions applied over the inlet plane. Madhavan and Kemmerling [21] found that the solutions due to different inlet spatial distribution, were qualitatively similar beyond 1.75 inlet diameters. Inlet secondary flow had most impact at less than one dimeter downstream and in particular during diastole. Expressed in terms of WSS, the differences between various conditions were about 26%.

The parameters of the Windkessel model were tuned to match the outlet flow-rates obtained from measured data. In one of the PS simulations, the model parameters were obtained from the

literature. The effect of the outlet conditions was estimated to be approximately 18% in terms of Time Averaged WSS (TAWSS). The outlet conditions showed significant effect on the flow in the vicinity of the outlets. The overall conclusion of Madhavan and Kemmerling’s [21] study was that it is important to use PS conditions, especially if details of the flow in the proximity of the inlet and outlet boundaries are of interest.

(9)

8 JESMDT-20-1007, Fuchs Frydrychowicz et al. [13] assessed in vivo WSS spatial distribution and Oscillatory Shear Index (OSI) in healthy human aortas by 4D-flow MRI and found that the mean absolute TAWSS was rather small (0.25 and 0.33 Pa) but contained an approximately 10% circumferential WSS component in all subjects investigated. The lowest absolute WSS and highest OSI were found to differ significantly from local mean values. The distribution of low WSS and high OSI resembled typical locations of atherosclerotic lesions at the lesser curvature of the aortic arch and its adjacent 3 large branches.

Youssefi et al. [18] considered the flow in the thoracic aorta using PS inflow velocity profile as compared to an idealized Womersley-type profile. Simulations having the same flow-rate, but different inlet profiles showed considerably different peak and radial velocities as well as flow symmetry. However, some similarity between patient derived and parabolic inflow profiles was observed. Youssefi et al. [18] concluded that fully patient-specific inflow boundary conditions must be used in order to produce meaningful results.Wei et al. [31] discussed the issue of time-averaged boundary conditions for surgical planning. The study comprehensively investigated the sensitivity of Fontan hemodynamic (with low cardiac output) metrics when using time-averaged BC:s, quantified by the difference between simulations using the time-averaged and pulsatile BC:s. Although time- averaged BC:s can save computational time, it is only possible to use such conditions for clinical purposes when the estimated sensitivity of the results to time-averaged BC has been established.

Muñoz et al. [16] provided results stressing the need to assess the sensitivity of the simulation results to BC:s. In that study, Color-Doppler and MRI were used to investigate the intracardiac flow in a group of healthy individuals, characterizing the flow in the proximal part of the ascending aorta. The results provide an illustrative qualitative information about the vortical flow in the left ventricle and

(10)

9 JESMDT-20-1007, Fuchs the ascending aorta. The volume flow-rate as function of time were found to show larger levels of variation among the individuals during the early ejection phase as compared to the ventricle relaxation phase, along with the observation that the flow through the aortic valve is non-uniform.

Apart from setting boundary conditions, given the fact that CFD simulations of blood flow in arteries require additional assumptions and models, the task of carrying out a true patient-specific flow “reconstruction” is very challenging. The emphasis of this study is to assess the uncertainties associated with varying time-dependency of (fixed) inlet flow- and heart-rates and different types of outlet BC:s. The flow in the thoracic aorta is characterized in terms of commonly used wall shear stress (WSS) related parameters. The impact of geometrical simplifications at the inlet, the need and length of extensions at the outlets, effects of flow-rate variations at the inlet velocity profile are addressed.

METHODOLOGY

Numerical methodology

In this study the blood flow in the human aorta is considered. The blood is modelled as an incompressible fluid characterized by a non-Newtonian viscosity depending on the RBC

concentration. The walls of the arteries are assumed to be rigid, thus, pressure wave propagation effects are negligible. Moreover, neglecting the density variation of the flowing blood, the Navier- Stokes equations is as follows;

(11)

10 JESMDT-20-1007, Fuchs i 0

i

u x

 

 (1a)

( i j) 1

i i

j i j j

u u u p u

t x x xx

       

     (1b)

where ui is the velocity component in the i:th direction, p,  and are the pressure, density and kinematic viscosity of the blood, respectively. Here, blood viscosity is assumed to obey the Quemada

model ( / , )

i j RBC

u x

where RBC is the concentration of the Red Blood Cells (RBC). In the results presented herein, RBC is assumed to be constant, having a value of 0.45. van Wyk et al 2015 [32] shows the details of the Quemada model.

Boundary condition requirement

The system of equations (1) is characterized as “incompletely parabolic”. In order to solve the system, equal number of conditions as the dimensions (D) of the problem need to be specified on all boundaries along with specification of initial conditions. These boundary conditions could be of Dirichlet, Neumann, mixed or of more complex character. Additionally, a necessary condition for the existence of a solution for the differential problem is that global mass conservation is satisfied (i.e.

the total mass flow into the domain equals the total mass flow out of the domain under

consideration). Similar BC:s for the differential problem are needed to solve Eqs. 1 numerically.

However, as almost all numerical schemes use collocated grids, an additional (numerical) BC need to be introduced on all boundary points. This additional condition is most often given in terms of the

(12)

11 JESMDT-20-1007, Fuchs pressure or its gradients. The numerical solution is obtained by solving the discrete momentum equations in time and updating the pressure in each time step by solving a Poisson equation (found by taking the divergence of equation (1b) and utilizing equation (1a)):

2 ( i j)

i i i j

p u u

x x x x

  

    (1c)

The Poisson equation is also elliptic, requiring one condition on all boundary points. Equation (1c) also implies that

i.) A change in flow variables on any boundary point propagates into the entire domain ii.) The maximum principle implies that the change in pressure BC will be largest at the boundary and propagate into the domain having smaller amplitude, depending on the used BC for (1c).

Focusing on boundary conditions, the following discussion is restricted to flow in the thoracic aorta and geometries that are topologically similar and where the inlet orifice of the thoracic aorta is considered to be at the inlet plane of the aortic valve or the sinotubular junction. The thoracic aorta has three main branches leading blood to the head, neck and upper limbs. The flow-rate in these vessels may vary substantially. In the literature, time-resolved volumetric flow-rates at the different branches are available [18-19,21]. Such BC:s assures total mass conservation, given that the

numerical mass fluxes on the boundaries add up to zero. Failing to reach such a level of accuracy, the iterative process will converge only to a level related to the discretization error in contrast to the former situation where convergence is possible to machine round-off.

(13)

12 JESMDT-20-1007, Fuchs Due to the lack of adequate spatial resolution, several studies examined different inlet spatial distributions with given time-dependent influx [18, 21, 30]. Commonly used spatial distributions include a flat profile (“plug flow”), parabolic flow profile (valid for slow stationary flow in pipes) or Womersley type profiles (valid for laminar time-dependent flow in a circular pipe). Youssefi [18] and Morbiducci [30] respectively found significant differences in hemodynamics when different inlet conditions were used. However, it was also found that single-component (plane average) inflow profiles performed similar when compared to the results obtained using three-component inflow profiles. The importance of temporal variation of the inlet profile was not examined.

Considering outlet BC, the following three types conditions are commonly used in literature:

i. Flow-rate based [19, 21],

ii. Pressure or pressure gradient based,

iii. So-called Windkessel type models based on analogy of the distal circulatory system with electrical circuit [32]. The Appendix describes of how different assumptions may lead to the different types of conditions.

Literature provides measurement data of the time-dependent flow-rate at different segments of the aorta and its branches providing data that can be applied as BC by assuming the spatial profile (e.g. plug flow or a parabolic profile). However, it should be noted that such data is valid for the given patient and for the specific time when the measurements were performed. Benim et al. and

Madhavan and Kemmerling [19, 21] both presented reference values for flow-rate distribution in the

(14)

13 JESMDT-20-1007, Fuchs aortic branches, reproduced in Table 1. For a given flow-rate, the outflow profile is commonly assumed to vary weakly over the cross section.

Often the inlet and outlet sections are extended in order to allow for the use of general (generic) velocity profiles and eliminating the risk for retrograde flow at an outlet boundary. Under such conditions it may be justified to impose simpler type of BC:s, such as assumed spatial

distribution of the velocity vector with a given flow-rate. Simplified BC are justifiable when the length scale (L) in the streamwise direction is much larger than in the cross-plane (l). When l/L << 1, different simplifications to the system of equations (1) are allowed. See Appendix for an extended discussion.

An approach for outlet BC that have gained support in recent years is based in patient-specific adaption of 2 or 3 element electric circuit analogues [33]. The 3-element electrical-circuit analog uses an outlet resistor (representing the peripheral arterial resistance) in serial with parallel

resistor/capacitor element (corresponding to the vessel compliance and peripheral resistance). The pressure and flow-rate are related through the electrical-circuit analog [21]:

1 dp ( 1 2) 1 2 dQ

R C p R R Q CR R

dt     dt (2) The three model coefficients (R1, R2 and C) are determined by fitting the model parameters to available data, see Table A1 (Appendix) in which the time-scale associated with the RC-circuit is provided. The RC-circuit has to account for vessel compliance distal to the artery. Hence, the response time is of the order of the diastole period allowing exponential decay during diastole. As may be noted, not all the Windkessel parameters in Table A1 satisfy this criterion and consistency of

(15)

14 JESMDT-20-1007, Fuchs the results with the underlying assumptions needed in order to make the computed results admissible.

Case setup: geometry and boundary conditions

The aortic geometry use in this study was derived from a CTA data of a healthy subject.

Segmentation of the CTA Dicom data was carried out using a combination of centerline based and threshold methods with manual adjustments afterwards to segment as much of the vessel lumen as possible [33]. The walls of geometry were smoothed.

In the following, the thoracic aorta is considered, maintaining only the main arterial branches;

brachiocephalic, left common carotid and left subclavian. The two geometries had the same shape and branches but differed in the inlet segment. The first geometry, G1, is characterized by an extension segment at the inlet, where the shape of the sinotubular junction was extended as a straight tube of five cross sectional diameters. The second thoracic aorta, G2, included the aortic sinus without extension. Both geometries are shown in Figure 1. All arterial branches were further extended (at least five mean diameters) by straight tubes having the shape as derived from the CTA derived cross sectional shape of each branch. In order to study potential difficulties associated with different boundary conditions at the outlet, the thoracic aorta outlet plane positioned at level of the diaphragm was left without extension.

Table 2 shows the references providing the inlet and outlet BC:s investigated. The temporal variation of the different inlet flow-rates is depicted in Fig 2. It needs to be stresses that not all profiles were measured aortic flows [19]. However, the underlying reasoning for including these

(16)

15 JESMDT-20-1007, Fuchs profiles was to having them serve as indicators of the sensitivity to and importance of the temporal variation of the inlet flow-rate profile. The instantaneous flow was distributed uniformly (so called plug flow) over the inlet plane aligned in a direction normal to the inlet plane. All flow-rate profiles were scanned manually and scaled after smoothing to produce a flow-rate of 5 L/min at a rate of 60 BPM.

Regarding the outlets, flow-rate and Windkessel BC:s were investigated (Table 2). The 3- element circuit Windkessel model parameters used are provided in Table A1 (Appendix) and according to the respective study found in literature.

Considering Figure 2, it is observed that some of the flow-rate profiles have retrograde flow in late systole. Furthermore, the inlet data also differ in the distance of the measuring planes from the aortic valve plane and the measuring technique. The purpose of including these inlet data in the study was that, in spite of the differences, it is instructive to assess the inflow conditions against the same outlet BC:s to quantify BC sensitivity and to understand the importance of the temporal variation of the flow rate. Similarly, using a given inlet flow-rate profile with different outlet BC:s allows for assessing the sensitivity of the different outlet BC:s. Moreover, it can be noted that the flow-rate as presented by [21] is a smoothed version of the profile of [35] with a minimum of zero flow-rate during diastole. Furthermore, as pointed out by [19], the inlet profile marked as BEN is in fact flow-rate profile into the Carotid artery (after the thesis of David N. Ku [36]), which is clearly noticed as it has a weaker systolic phase and a considerable higher flow during diastole as compared to the other profiles.

(17)

16 JESMDT-20-1007, Fuchs Altogether, 9 different inlet profiles (Figure 2) and 7 different outlet BC:s (Table 2) were assessed, leading to 63 cases for each of the two geometries G1 and G2. Simulations were also carried out using refined grids and different heart-rates. The governing equations were discretized using formally a second order discretization in space and time. The discrete system of equations was integrated in time using an implicit method with time steps varying between 0.05 ms to 0.1 ms, depending on the grid size, as described below. The numerical integration was carried out using OpenFOAM 5.0 and post-processing was done using Paraview along with the needed python and MATLAB scripts.

The research was approved by the Swedish regional ethical vetting board in Linköping (Project identification code: DNR 2017/258-31, Prof. Anders Persson).

RESULTS Grid sensitivity

In order to assess the accuracy of the numerical simulations, each of the different geometries were studied using a set of three sequentially refined grids. The grids for the thoracic aorta consisted of 2.5, 5 and 10 Million cells, respectively. The accuracy studied showed that after three heart cycles, the computed velocities differed approximately 0.2 % between the “coarse” and the medium grid and about 0.1 % between the finest and medium grid (by L2-norm averaged from a sample of 500 points in space).

(18)

17 JESMDT-20-1007, Fuchs The length scales near the inlet and outlet sections were estimated using the computed solution for each case. The axial (velocity based) length scale was assessed by Lu=|UCL|/|dUCL/dCL| where UCL and CL represent the axial velocity and the distance along the center-line, respectively. As the flow develops into an ideal flow in a closed vessel (under steady flow conditions and for a very long vessel), the term |dUCL/dCL| decreases, implying a co-current increase of the length scale. The large scales in the cross-planes are bounded from above by the size of the vessel. For axial length- scales larger by about a factor of 20-50 relative the vessel diameter, approximate BC:s relying on neglecting axial variations may be used. Different levels of approximations and assumptions are discussed in the appendix. In this study, vanishing pressure gradient was considered reasonable when the length scale based of the first to second derivative of the pressure in the axial direction relative to the vessel size was O(104) or larger.

General flow characteristics and Helicity - sensitivity of results to imposed boundary conditions

A general impression of the time-dependent flow in geometries G1 and G2 are provided as animations of the streamwise velocity profiles (supplementary material). Figure 3 depicts the flow in the aorta using the inlet conditions of HOPE and outlet condition of Ben. Fig 3a shows four time instances of the flow for G1 at mid systole, early diastole, mid diastole and start of systole. Fig 3b, depicts the corresponding cases, i.e. using the same inlet and outlet conditions, for G2. The effect of the two inlet and two outlet conditions can be elucidated in terms of defects (non-uniformities) in the streamwise velocity. Moreover, intermittent flow separation in the inner/lesser curvature of the

(19)

18 JESMDT-20-1007, Fuchs aortic arch near peak systole, retrograde flow in the early diastole is also observed. The formation of helical motion is observed as the region of large axial velocity changes with time. The corresponding animations (supplementary material) clearly show this spiraling motion.

Combinations of different inlet conditions and the shape of the ascending aorta and the aortic arch may lead to helical/spiraling flow in both clockwise and counter-clockwise directions. Most commonly, the flow is spiraling in the clockwise direction. Figure 4 presents stream tracers colored by the helicity (defined as the scalar product of the local velocity vector and the local vorticity vector). It is shown that the spiraling motion is generated during the later accelerating part of systole and the strength of vorticity and helicity start to decay during the deceleration part of systole and early diastole. Depending on the time-variation of the inlet flow-rate, the spiraling motion may change direction in the ascending aorta/aortic arch. The strength of both helicity and vorticity decay substantially during diastole. For the given flow rate, the level of small-scale vorticity decreases at end-diastole and may be up to two orders of magnitude smaller at the descending part of the aorta as compared to the peak value at peak systole. Moreover, depending on the temporal variation in the flow-rate, the helical motion may change direction between late systole and early diastole, e.g. Fig 4b and 4c, second and third columns, (BERT-Ben and REGN-Ben). The change of direction starts at the aortic arch and propagates down to the descending aorta. Later in diastole, the weak spiraling motion may take the original direction at the descending aorta (fig 4, second to fourth columns).

(20)

19 JESMDT-20-1007, Fuchs Wall Shear Stress – sensitivity of results to boundary condition

To assess the effect of change in inlet geometry and sensitivity to BC:s, the following

commonly used WSS related parameters were used: Time Averaged Wall Shear Stress (TAWSS), the Oscillatory Shear Index (OSI), and the Relative Reference Time (RRT) [37-38].

The effect of the aortic sinus with the BEN inlet and Ben outlet flow-rates, were 8%, 26% and 63% for OSI, RRT and TAWSS, respectively. The corresponding numbers using the MAD inlet with Ben outlet conditions were, 71%, 59% and 165% for OSI, RRT and TAWSS, respectively. Table 3 displays in more detail the (spatial) mean values of these parameters, where Table 2 contains combinations of inlet- and outlet-BC:s for G1. It should be noted that although the two outlet BC:s used are similar, the different WSS parameters show a substantial deviation among the cases investigated, even in terms of mean values. This deviation is even more significant if expressed in terms of peak values.

Moreover, the spatial fluctuations are also significant, being of the same order as the mean values.

This effect can be observed qualitatively in the Figs 5 and 6 as well as the corresponding animations in the supplementary material.

Typical distribution of OSI for a selection of different cases is depicted in Figures 5 and 6. Figs 5a and 6a relate to different inlet flow-rate profiles for a given set of outlet BC:s for geometry G1 and G2, respectively. The corresponding cases for fixed inlet condition but having different outlet BC:s are shown in Figs 5b and 6b. As expected, the main effect of the BC is larger closest to the location of the corresponding boundary. However, the effects are not limited to those regions and the extent of the BC:s influence depends strongly on the deviation from the ideal conditions assumed by each BC.

(21)

20 JESMDT-20-1007, Fuchs The strength of the WSS in terms of peak and root mean square values depends rather strongly on the BC:s, see Table 4 for three of the inlet flow-rates combined with six possible outlet BC:s. As observed, the peak values of the WSS may be larger than the mean by as much as three orders of magnitude. The peak values depend more strongly on the BC:s as compared to the

corresponding mean values. The mean values of WSS of the different combinations of BC, are of the same order as measured by [13]. As the spatial resolution in their MRI data was of the order of about 2 mm (more than one order of magnitude larger than the coarse cells in our simulations), the

difference between the measured and computed WSS and OSI can be attributed to the smoothing effect of the MRI resolution and the filtering of the data.

Windkessel model sensitivity

The Windkessel models were applied together with the different inlet conditions for both thoracic aorta configurations (G1 and G2). For some inlet profiles, the different Windkessel based models converges. For others, the simulations enter a limit-cycle or diverge. The latter cases occur when the underlying assumptions of the 0D models (see Appendix) are not satisfied, for example when the length-scale ratio in the streamwise and spanwise directions is of order one or when there is flow reversal in some parts of the vessel not too far from the outlet. In order to assess the

sensitivity of the Windkessel parameters found in literature, measured flow and pressure in a canine aorta were used [36]. Given the temporal flow-rate and using the 3-element circuit along with equation (2), it is possible to compute the temporal variation of the pressure as predicted by the different Windkesssel model parameters. In these calculations, the same initial condition (t = 0) was applied on the pressure as given by the measured pressure at this time. Figure 7 depicts the temporal

(22)

21 JESMDT-20-1007, Fuchs variations of pressure predicted by equation (2) using the different model parameters (Table A1) as well as the measured pressure and flow-rate versus time. The computed pressure was rescaled by a factor 0.1 for the parameters provided by [22-23,33].

DISCUSSION AND CONCLUSIONS

This study shows the importance of the time-variation of the cardiac flow-rate even with given mean flow- and heart-rates. By rescaling nine different arterial velocity profiles to a cardiac output of 5 LPM (at 60 BPM), the effect of cardiac contraction and relaxation rates in aortic flow, were assessed. The helical motion, commonly clockwise, in the thoracic aorta and well documented by MRI techniques was found to depend on the temporal variation of the inlet flow rate. This in turn lead to that the helical motion disappeared or changed direction in the ascending aorta. The strength of the helical motion was found to reach its peak close to the peak flow rate in systole, characterized by a quickly decreasing amplitude during diastole. The location of the change in helical motion was closely related to the location where helicity changed sign (Figs 4). The time-dependent

helical/spiraling motion also induced a time-dependent WSS strength and direction (parallel to the wall). Regarding helical motion, one approach to study the underlying physics in more detail could be to consider the different terms in the vorticity transport equation (not discussed in the paper).

Similarly, the potential relation between the helical motion and its impact on the local concentration of blood components (cells and macromolecules) and the transport process into and beyond the endothelium of the artery were not considered in this study.

(23)

22 JESMDT-20-1007, Fuchs The sensitivity of the WSS related parameters to the seven different outlet BC:s was also assessed. The results show that the CFD based simulations must be checked explicitly for possible errors and potential inconsistencies originating from the underlying assumptions of the applied outlet BC:s. The appendix of the paper details these assumptions for several common outlet BC:s. The impact of geometrical simplifications at the inlet, the need and length of extensions at the outlets, effects of flow-rates and the effect of strength of shear in the inlet velocity profile were presented shortly. The results presented were aimed at stressing the need for caution when using CFD tools for predicting blood flow in clinical situations. The main reason for such caution is due to the large natural variability of cardiac flow also over relatively short periods of time. The variability may lead to strong impact on the hemodynamics (spiraling motion, retrograde flow, WSS variation in space and time) and therefore a single or a small number of simulations are inadequate for capturing potential pathologies. Hence, simulations should be carried out for the “flight envelope” of the heart and a smaller range of possible flow distribution among the main arteries branching from the aorta.

It should be noted that the assumptions made in this study imply inherent limitations to the conclusions of the paper. Two of the assumptions that were made were related to blood rheology and assuming that the aorta had rigid walls. The effects of RBC concentration on blood viscosity is known to be rather small for larger shear-rates. With increasing age, the stiffness of the arteries increases and thereby the physiological importance of vessel compliance decreases. Outlet BC:s may be challenging from computational point of view when flow reversal occurs close to an outlet BC:s.

Such difficulties were encountered for certain combinations of inlet flow-rate profiles and certain 3-

(24)

23 JESMDT-20-1007, Fuchs element Windkessel model parameter set-ups. The reason for this originated from the underlying assumptions for the BC models being violated and therefore, the numerical solutions were not admissible due to inconsistency.

The effect of inlet BC:s, for a given set of outflow BC:s and the effects of different outlet BC:s for each given inlet were assessed using different parameters. OSI, RRT and TAWSS are commonly believed to be related to the formation of atherosclerosis, thus relevant as rough estimative indicators for judging the sensitivity of the numerical results. The results presented are mainly in terms of OSI. However, in the supplementary material additional data regarding the temporal development of the flow itself is found. Such animations expose the complexity of the time-

dependent flow field, formation of time-varying larger structures (separated regions, retrograde and helical flows). As mentioned above, considering the vorticity equation, BC:s can be related to the formation of such structures. However, this aspect was not explored here as the main goal of the study was to assess the sensitivity of certain types of BC:s found in the literature. The challenge still remains on how to interpret and to quantitatively show the relationship between flow structures and vessel pathologies. The simulated results indicate that these parameters are rather sensitive to the temporal distribution of the inlet flow-rate. The changes in these parameters varied in the range of up to 260% in terms of mean values and up to a factor 6 for peak values. Commonly used

parameters, such as OSI, RRT, TAWSS and other averaged WSS related expression can only indicate potential pathology but not possible progress, nor its underlying mechanisms.

(25)

24 JESMDT-20-1007, Fuchs Due to the large range of and quick changes in cardiac output, simulations for clinical

purposes has to cover the entire range of BC:s, in contrast to current, single shot (set-up of geometry, BC:s, rheological and wall models) in order to give truly patient-specific simulations.

ACKNOWLEDGMENT

We thank Associate Professor Chunliang Wang and Professor Örjan Smedby in Royal Institute of Technology (Dept. of Biomedical Engineering and Health systems) for allowing us to use the segmentation software Mialab. The computations were carried out on computer resources at NSC at Linköping University and HPC2N at Umeå University through a SNIC grant. Fuchs, A., acknowledges the support from the department of radiology at Linköping University Hospital.

FUNDING

This research used faculty funding at the Royal Institute of Technology, KTH (Prahl Wittberg and Berg) and partial support by the regional county of Östergötland (Fuchs).

SUPPLEMENTARY MATERIAL

All the additional files show an animation of how the velocity component, tangential to the aortic centreline, varies over a heart cycle at certain planes. Each animation is from simulations on the thoracic aorta with aortic sinus (G2) using different inlet- and outlet boundary conditions (denoted as in main body of the text).

(26)

25 JESMDT-20-1007, Fuchs Additional_file_1:

Tangential velocity over a heart cycle. Inlet condition: Benim. Outlet conditions: Benim.

Additional_file_2:

Tangential velocity over a heart cycle. Inlet condition: Benim. Outlet conditions: Madhavan.

Additional_file_3:

Tangential velocity over a heart cycle. Inlet condition: Hope. Outlet conditions: Benim.

Additional_file_4:

Tangential velocity over a heart cycle. Inlet condition: Rengier. Outlet conditions: Benim.

(27)

26 JESMDT-20-1007, Fuchs NOMENCLATURE

ui Velocity

p Pressure

ρ Density

Ν Kinematic viscosity αRBC

4D-MRI

Concentration of Red Blood Cells (hematocrit) Four Dimensional Magnetic Resonance Imaging BC Boundary Condition

BPM CFD

Beats Per Minute

Computational Fluid Dynamics CTA

LPM

Computed Tomography Angiography Liters Per Minute

MRI Magnetic Resonance Imaging

OSI PC-MRI

Oscillatory Shear Index

Phase Contrast Magnetic Resonance Imaging PS Patient Specific

RBC Red Blood Cell

(28)

27 JESMDT-20-1007, Fuchs RRT Relative Residence Time

TAWSS Time Averaged Wall Shear Stress

US Ultrasound

WSS Wall Shear Stress

REFERENCES

[1] Rodeheffer, R. J., Gerstenblith, G., Becker, L. C., Fleg, J. L., Weisfeldt, M. L. and Lakatta, E. G., 1984, “Exercise cardiac output is maintained with advancing age in healthy human subjects:

cardiac dilatation and increased stroke volume compensate for a diminished heart rate,”

Circulation., 69(2), pp. 203-213. DOI: 10.1161/01.CIR.69.2.203

[2] Massoudi, M., Kim, J. and Antaki J. F., 2011, “Modeling and numerical simulation of blood flow using the Theory of Interacting Continua,” Int. J. Non. Linear. Mech., 47(5), pp. 506–520. DOI:

10.1016/j.ijnonlinmec.2011.09.025.

[3] Segers, P., Stergiopulos, N., and Verdonck, P., 1997, “A non-invasive pulse pressure method for the estimation of total arterial compliance,” Computers in Cardiology 1997. IEEE., pp. 171-174.

DOI: 10.1109/CIC.1997.647858

[4] Stergiopulos, N., Meister, J. J., and Westerhof, N., 1995, “Evaluation of methods for estimation of total arterial compliance,” Am. J. Physiol., 268(4): H. 1540-1548. DOI:

10.1152/ajpheart.1995.268.4.H1540

[5] Joseph, J., Venkataraman, J., Kumar, V.J., and Suresh, S., 2011, “Non-invasive estimation of arterial compliance,” IEEE International Instrumentation and Measurement Technology Conference, 5, pp. 1-5. DOI: 10.1109/IMTC.2011.5944282

[6] Stein, P. D., and Sabbah, H. N., 1976, “Turbulent blood flow in the ascending aorta of humans with normal and diseased aortic valves,” Circ. Res., 39(1), pp. 58-65. DOI: 10.1161/01.res.39.1.58 [7] Prahl Wittberg, L., van Wyk, S., Fuchs, L., Gutmark, E., Backeljauw, P., and Gutmark-Little I., 2016,

“Effects of aortic irregularities on blood flow,” Biomech. Model. Mechanobiol., 15(2), pp. 345-360.

DOI:10.1007/s10237-015-0692-y.

(29)

28 JESMDT-20-1007, Fuchs [8] Ha, H., Ziegler, M., Welander, M., Bjarnegård, N., Carlhäll, C.-J., Lindenberger, M., Länne, T.,

Ebbers, T., and Dyverfeldt, P., 2018, “Age-Related Vascular Changes Affect Turbulence in Aortic Blood Flow,” Front. Physiol., 9(36) DOI: 10.3389/fphys.2018.00036

[9] Ziegler, M., Lantz, J., Ebbers, T., and Dyverfeldt, P., 2017, “Assessment of turbulent flow effects on the vessel wall using four-dimensional flow MRI,” Magn. Reson. Med., 77(6), pp. 2310-2319. DOI:

10.1002/mrm.26308

[10] Zajac, J. Eriksson, J., Dyverfeldt, P., Bolger, A. F., Ebbers T., and Carlhall C.-J., 2015, “Turbulent kinetic energy in normal and myopathic left ventricles,” J. Magn. Reson. Imaging., 41(4), pp.

1021–1029, DOI: 10.1002/jmri.24633

[11] Andersson, M., Lantz, J., Ebbers, T., Karlsson, M., 2015, “Quantitative assessment of turbulence and flow eccentricity in an aortic coarctation: impact of virtual interventions,” Cardiovasc. Eng.

Technol., 6(3), pp. 281–293. DOI: 10.1007/s13239-015-0243-9

[12] Fredriksson, A., Trzebiatowska-Krzynska, A., Dyverfeldt, P., Engvall, J., Ebbers T., and Carlhall C.-J., 2018, “Turbulent kinetic energy in the right ventricle: Potential MRmarker for risk stratification of adults with repaired tetralogy of fallot,” J. Magn. Reson. Imaging., 47(4), pp. 1043–1053. DOI:

10.1002/jmri.25830

[13] Frydrychowicz, A., Stalder, A. F., Russe, M. F., Bock J., Bauer, S., Harloff, A., Berger, A., Langer, M., Hennig, J., and Markl, M., 2009, “Three‐dimensional analysis of segmental wall shear stress in the aorta by flow‐sensitive four‐dimensional‐MRI,” J. Magn. Reson. Imaging., 30(1), pp. 77–84. DOI:

10.1002/jmri.21790

[14] Rengier, F., Delles, M., Unterhinninghofen, R., Ley, S., Muller-Eschner, M., Partovi, S., Geisbusch, P., Dillmann, R., Kauczor, H. U., and von Tengg-Kobligk, H., 2012, “In vivo and in vitro validation of aortic flow quantification by time-resolved three-dimensional velocity-encoded MRI,” Int. J.

Cardiovasc. Imaging., 28(8), pp. 1999–2008. DOI:10.1007/s10554-012-0027-3.

[15] Bertelsen, L., Svendsen, J. H., Køber, L., Haugan, K., Højberg, S., Thomsen, C., and Vejlstrup, N., 2016, “Flow measurement at the aortic root - impact of location of through-plane phase contrast velocity mapping,” J. Cardiovasc. Magn. Reson., 18(1):55. DOI:10.1186/s12968-016-0277-7.

[16] Muñoz, R. D., Markl, M., Mur, J. L. M., Barker, A., Fernandez-Golfın, C., Lancellotti, P., and Gomez J. L. Z., 2013, “Intracardiac flow visualization: current status and future directions,” Eur. Heart. J.

Cardiovasc. Imaging., 14(11), pp. 1029–1038. DOI:10.1093/ehjci/jet086

[17] Hope, T. A., Markl, M., Wigström, L., Alley, M. T., Craig Miller, D., and Herfkens R. J., 2007,

”Comparison of Flow Patterns in Ascending Aortic Aneurysms and Volunteers Using Four- Dimensional Magnetic Resonance Velocity Mapping,” J. Magn. Reson. Imaging, 26(6), pp. 1471- 1479. DOI:10.1002/jmri.21082.

[18] Youssefi, P., Gomez, A., Arthurs, C., Sharma, R., Jahangiri, M., and Figueroa C. A., 2018, “Impact of Patient-Specific Inflow Velocity Profile on Hemodynamics of the Thoracic Aorta,” J. Biomech. Eng.

140(1). DOI: 10.1115/1.4037857

(30)

29 JESMDT-20-1007, Fuchs [19] Benim, A. C., Nahavandi, A., Assmann, A., Schubert, D., Feindt, P., and Suh, S.H., 2011, “Simulation

of blood flow in human aorta with emphasis on outlet boundary conditions,” Applied Mathematical Modelling, 35(7), pp. 3175–3188. DOI: 10.1016/j.apm.2010.12.022

[20] Liu, X., Fan, Y., Deng, X., and Zhan, F, 2011, “Effect of non-Newtonian and pulsatile blood flow on mass transport in the human aorta,” J. Biomech., 44(6), pp. 1123-1131.

DOI:10.1016/j.jbiomech.2011.01.024

[21] Madhavan, S., and Kemmerling, E. M. C., 2018, “The effect of inlet and outlet boundary conditions in image-based CFD modeling of aortic flow,” Biomed. Eng. Online., 17(1):66. DOI:

10.1186/s12938-018-0497-1

[22] Pirola, S., Cheng, Z., Jarral, O. A., O’Regan, D.P., Pepper, J.R., Athanasiou T., and Xu X.Y., 2017, “On the choice of outlet boundary conditions for patient-specific analysis of aortic flow using

computational fluid dynamics”. J. Biomech., 60, pp. 15–21. DOI: 10.1016/j.jbiomech.2017.06.005 [23] Zhu, Y., Chen, R., Juan, Y.-H., Li, H., Wang, J., Yu, Z., and Liu, H., 2018, ”Clinical validation

and assessment of aortic hemodynamics using computational fluid dynamics simulations

from computed tomography angiography”. BioMed. Eng. OnLine., 17(1):53. DOI: 10.1186/s12938- 018-0485-5

[24] Guan, D., Liang, F., and Gremaud, P. A., 2016, “Comparison of the Windkessel model and structured-tree model applied to prescribe outflow boundary conditions for a one-dimensional arterial tree model”. J. Biomech., 49(9), pp. 1583–1592. DOI: 10.1016/j.jbiomech.2016.03.037 [25] Sankaran, S., Kim, H. J., Choi, G. and Taylor, C. A., 2016, “Uncertainty quantification in coronary

blood flow simulations: Impact of geometry, boundary conditions and blood viscosity,” J.

Biomech., 49(12), pp. 2540–2547. DOI: 10.1016/j.jbiomech.2016.01.002

[26] Bozzi, S., Morbiducci, U., Gallo, D., Ponzini, R., Rizzo, G., Bignardi, C., and Passoni. G.,2017,

”Uncertainty propagation of phase contrast MRI derived inlet boundary conditions in

computational hemodynamics models of thoracic aorta,” Computer Methods in Biomechanics and Biomedical Engineering, 20(10), pp. 1104-1112. DOI: 10.1080/10255842.2017.1334770 [27] Chandra, S., Raut, S.S., Jana, A., Biederman, R.W., Doyle, M., Muluk, S.C., and Finol, E.A., 2013,

“Fluid-Structure Interaction Modeling of Abdominal Aortic Aneurysms: The Impact of Patient- Specific Inflow Conditions and Fluid/Solid Coupling,” ASME J. Biomech. Eng., 135(8). DOI:

10.1115/1.4024275

[28] Campbell, I.C., Ries, J., Dhawan, S.S., Quyyumi. A.A., Taylor, W.R., and Oshinski, J.N., 2012, “Effect of Inlet Velocity Profiles on Patient-Specific Computational Fluid Dynamics Simulations of the Carotid Bifurcation,” ASME J. Biomech. Eng., 134(5). DOI: 134:5:051001.

[29] Morbiducci, U., Gallo, D., Massai, D., Consolo, F., Ponzini, R., Antiga, L., Bignardi, C., Deriu, M., and Redaelli, A., 2010, “Outflow Conditions for Image-Based Hemodynamic Models of the Carotid Bifurcation: Implications for Indicators of Abnormal Flow,” J. Biomech. Eng., 132(9). DOI:

10.1115/1.4001886.

(31)

30 JESMDT-20-1007, Fuchs [30] Morbiducci, U., Ponzini, R., Gallo D., Bignardi, C., and Rizzo, G., 2013, “Inflow Boundary Conditions

for Image-Based Computational Hemodynamics: Impact of Idealized Versus Measured Velocity Profiles in the Human Aorta,” J. Biomech., 46(1), pp. 102–109. DOI:

10.1016/j.jbiomech.2012.10.012

[31] Wei, Z., Trusty, P.M., Tree, M., Haggerty C. M., Tang E., Fogel M., and Yoganathan A.P., 2017, “Can time-averaged flow boundary conditions be used to meet the clinical timeline for Fontan surgical planning?,” Journal of Biomechanics., 50, pp. 172–179. DOI: 10.1016/j.jbiomech.2016.11.025 [32] van Wyk, S., Prahl Wittberg, L., Bulusu, K. V., Fuchs, L., and Plesniak, M. W., 2015, “Non-

Newtonian perspectives on pulsatile blood-analog flows in a 180 curved artery model,” Physics of Fluids., 27(7), 071901. DOI: 10.1063/1.4923311

[33] Vignon-Clementel, I .E., Figueroa, C. A., Jansen, K. E., and Taylor, C .A., 2010, “Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries,” Comp. Meth. Biomech. and Biomedical Eng., 3(5), pp. 625-640. DOI:

10.1080/10255840903413565

[34] Schaap, M., Metz, C. T., van Walsum, T., van der Giessen, A. G., Weustink, A. C., Mollet, N. R., Bauer, C., Bogunović, H., Castro, C., Deng, X., Dikici, E., O’Donnell, T., Frenay, M., Friman O., Hernández Hoyos, M., Kitslaar, P. H., Krissian, K., Kühnel, C., Luengo-Oroz, M. A., Orkisz, M., Smedby, Ö., Styner, M., Szymczak, A., Tek, H., Wang, C., Warfield, S. K., Zambal, S., Zhang, Y., Krestin, G. P., and Niessen, W. J., 2009, “Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms,” Medical Image Analysis, 13(5), pp. 701-714. ISSN 1361-8415. DOI: 10.1016/j.media.2009.06.003.

[35] Fuster, V., Walsh, R. A., and Harrington, R. A., 2011, “Pathophysiology of Heart Failure”. Hurst's The Heart 13th Ed. McGraw Hill., 721.

[36] Ku, D. N., 1997, “Blood flow in arteries,” Annu. Rev. Fluid Mech., 29, pp. 399–434. DOI:

10.1146/annurev.fluid.29.1.399

[37] He, X., and Ku, D. N., 1996, “Pulsatile flow in the human left coronary artery bifurcation: average conditions”. J Biomech. Eng., 118, pp. 74–82. DOI: 10.1115/1.2795948

[38] Fuchs, A., Berg, N., and Prahl Wittberg, L., 2019, “Stenosis Indicators Applied to Patient-Specific Renal Arteries without and with Stenosis,” Fluids, 4(1):26. DOI:10.3390/fluids4010026.

[39] Taylor, C. A., Hughes, T. J. R., and Zarins, C. K., 1998, “Finite element modeling of three‐

dimensional pulsatile flow in the abdominal aorta: relevance to atherosclerosis,” Ann. Biomed.

Eng., 26, pp. 975–987. DOI: 10.1114/1.140

(32)

31 JESMDT-20-1007, Fuchs APPENDIX: OUTFLOW BC

We assume that the section of the aorta in proximity to inlet boundary have no branches, the vessel curvature is mild and the flow in the streamwise direction has a much larger length scale relative to that in the spanwise directions. As noted before, equation (2b), the momentum equations, require BC:s on all boundary points. It can be shown that if the flow direction is laminar and “uni- directed” the effects of the downstream data is limited to a distance proportional to the square root of the Reynolds number. This effect occurs when the length scale in the streamwise direction is considerably larger than in the spanwise directions. By using such a scale relation, the size of the different terms in equations (2b) can be estimated.

Two main types of approximations may be derived. One is based on neglecting the smallest terms and the second is based on maintaining the largest terms in the governing equations. Further assumptions on two set of approximations leads ultimately to a simpler 1D and 0D approximations. In the following, a short description of the different approximations that can applied to set appropriate outlet BC:s on arteries are provided.

2.5D model (BC)

More systematic BC:s may be derived by making assumptions on the flow in the vicinity of the outlet boundaries. Commonly, the inlet and outlet sections are extended so that generic velocity profiles are applicable and at the same time allowing for the flow at the region of interest to develop

(33)

32 JESMDT-20-1007, Fuchs with less interference from the error due to the boundary conditions. This behavior is based on the fact or assumption that the length scales are directionally dependent and the ratio of scales is not close to unity. Suppose that the coordinates () are used, where  is aligned with the main flow (streamwise) direction and (are coordinates in the plane normal to the vector in the direction

If assuming that the length scale in  is L and in the other direction is l, such that l/L<<1, then equation system (1) may be simplified in different ways, depending on further assumptions.

One possibility is to neglect solely the smallest term of the system, namely the diffusive term in the direction. This assumption implies that the diffusion of momentum in the streamwise direction is neglected. Thus, the viscous term would contain only derivatives in the spanwise plane, implying that no BC have to be specified in the streamwise direction if the pressure gradient in that direction is known. This pressure gradient can be found by assuming that the pressure gradient in the streamwise direction is the same throughout the spanwise plane (due to the weak spanwise flow that is supported by the spanwise components of the pressure gradient). With this assumption, the

pressure p can be split, depending on the size of the pressure gradients; p=P()+p’(). Fuchs and Zhao [1-2] showed that under such conditions, the streamwise momentum equations can be

integrated, depending only on the upstream conditions and adjusting P() to enforce the mass flow- rate through the cross-section. To simplify the discussion, the approach is applied to the flow in Cartesian coordinates with the x-direction be aligned with the streamwise direction () and the y-z- plane to be the spanwise plane.

The reduced momentum equations become:

(34)

33 JESMDT-20-1007, Fuchs

u ( u u u) dP u u

u v w

t x y z dx y y z z

   

              

        (A1)

v ( v v v) p v v

u v w

t x y z y y y z z

   

              

         (A2)

( )

w w w w p w w

u v w

t x y z z y y z z

   

              

         (A3)

where all

xx

were neglected and the pressure was split into P(x)+p(y,z). The pressure P may be computed by imposing global continuity; i.e. the flow through each spanwise plane is a given constant, C; i.e. the mass flow-rate through the vessel is provided by;

( , , )

A

u x y z dydz C

 



(A4)

with u being the velocity in the x-direction. Note also that the continuity equation must be satisfied.

u v w 0

t x y z

   

   

   

    (A5)

The simplified equations still needs BC:s on all boundaries except on the outlet boundary (x=Xout).

The new system of equations has the benefit that it may be integrated more easily. With a given inlet velocity profile, the inlet mass flow-rate and hence the constant C can be obtained.

Dividing the x-direction into planes with small distance in-between (denoted by x, not

(35)

34 JESMDT-20-1007, Fuchs necessarily a constant), u and dP/dx (equations (A1) and (A4) can be computed at plane x=xi

given that u,v,w and p are known at a the spanwise plane (x=xi-1) considered. Once u and dP/dx are computed at x=xi, the 2D (Navier-Stokes) equation for v.w and p at the same spanwise plane (x=xi) can be solved. Repeating this procedure enables integration of the 3D Navier- Stokes equation without the need for downstream BC. However, it is important that the solution is examined for consistency with the assumptions. This implies the following has to be verified:

a. The length scale in the x-direction is (much) larger than in the y- and z-directions and b. The pressure p is a weak function of y and z and the gradients are small as compared to the size of dP/dx.

Finally, the size of the neglected terms can be estimated:

; ;

u v v

xx xx xx

     

     

relative to the advection term and/or the time-derivatives. When the consistency criteria are not met, the computational domain (extending it in the x-direction) may need to be prolonged until the estimated error of the approximation reaches an acceptable level.

When reached an acceptable relative error in the assumptions, the algorithm described above, enables determining the outflow BC without the need for additional assumptions.

(36)

35 JESMDT-20-1007, Fuchs 1D model

A more restrictive assumption can be made on the system of equations. Instead of neglecting only the smallest terms, only the largest terms are kept. This mean that the flow in

the spanwise plane is negligible as compared to the streamwise flow. Altogether, such assumptions lead to the so called 1D model for flow in compliant vessels [3-4]. Hence, from equations (A1) and (A5), the following is obtained:

w

u u dP

t u x dx x

  

 

and u 0

t x

 

  

  (A6a & A6b)

Integrating equations (A1) and (A5b) in the plane (denoting the cross-sectional area as A and the plane averaged streamwise flow-rate by Q) results in the following (as done by [5]):

( 2/ )

w

Q Q A A dP

t x dxdl

   

 

(A7a)

A Q 0

t x

(A7b)

where the last term in equation (A6a) is integral of the wall shear-stress (w) over the boundary () of the spanwise (cross-sectional) plane. In some studies, the last term on the right hand side is simplified by assuming a power law streamwise velocity profile in a pipe having a radius R:

w

dl KQ

A

(A7c)

(37)

36 JESMDT-20-1007, Fuchs Here, K < 0. This may be easily observed for the limiting case for which the x-derivative vanishes and the dissipation due to viscosity will lead to reduction in the kinetic energy over time if K is negative.

In the equations above, possible compliance effects of the vessel wall were omitted.

However, such effects have been incorporated in several 1D simulations [5,6], although for the boundary model above, the compatibility of the assumptions with the underlying assumptions need to be assessed.

0.5D model

Furthermore, the 1D equation can be used to further simplify the model. For a vessel with constant or slowly varying cross-section and cross-sectional flow-rate in the streamwise direction, equations (A7a) and (A7c) simplify to

dQ Q dP

K A

dt A dx (A8)

Relation (A8) indicates that the pressure gradient has to be negative in order to maintain the flow. It can also be observed that the spatial variation of the pressure gradient may be assumed to be small whereby the pressure gradient (taken from the upstream) can be specified and through equation (A8) detail the flow-rate on the outlet boundary.

Assuming dP/dx=G*Q(t-t0), where t0 is the lag in time between the flow and pressure. The linear problem can be solved for each Fourier mode to reconstruct the pressure gradient if t0 is given.

(38)

37 JESMDT-20-1007, Fuchs

0D models (Windkessel)

Further assumptions may lead to the so called Windkessel model, originally developed in analogy to passive electrical circuits. Zero-dimensional (0D) models imply that the total pressure

( ) ( )

2

( ) /

P t

t

P tQ tA

is independent of x. It is assumed that, downstream the vessel, there exists a system to which a force is exerted on the boundary. Thus, the right-hand side of equations (A6a and 7a) would be non-zero;

( 2/ )

( BC)

Q Q A A dP Q

K x x F

t x dx A

     

  (A8)

where (x-xBC) is the Dirac delta function, and F is the force acting on the boundary point xBC. It could be argued that F should contain a pressure term and its temporal derivative, and no spatial derivative (to be 0D). This assumption may be clarified by considering the 1D equation. If Q(t) is expressed as Q(t) = qeit then the time-derivative of Q and Q itself (in A8) are out of phase relative to each other and can be balanced only at steady state. To balance the equation, F should be proportional to (P- p0)+dP/dt where p0 and  are model parameters. This assumption allows for satisfying the simplified 1D equation for time-dependent flows with Q=Q(t).

Altogether, such 0D model would take the form of:

( 0)

dQ dP

Q P p

dt dt (A9)

References

Related documents

A thermal study of an inlet textured contact has been per- formed for a flow of lubricant with temperature dependent viscosity, density and specific heat capacity submitted

In addition to the AAIM and corresponding atlases, we have developed and described a method for intracranial 4D flow MRI vessel segmentation and flow quantification

The New- tonian case has an average friction Reynolds number Re τ equal to 183, while the polymer one provides a value of 155 which corresponds to a drag reduction DR of

Ingold 2007b; Latour, 1988, 2005; Orlikowski, 2007; Scott &amp; Orlikowski, 2014), I argue in this thesis that traditional dualistic assumptions, and the consequential tendency

Transient pressure measurements in the high and low-pressure regions of the turbine at different loads point out the unsteadiness of the flow in both regions, especially in the

Left ventricular blood flow in acquired heart disease (Paper III) Previous 4D flow CMR studies have suggested marked abnormalities in 4D flow volume and KE in severely dilated

Quantification of Wall Shear Stress using Large Eddy SimulationT.

Den första besvarar han med att det är de sex globaliseringsfaserna som fört oss hit och ur dessa kommer den inriktning som vi kallar modernismen, som spridit sig över plane- ten