Pulsatile Aortic Blood Flow—A Critical Assessment of Boundary Conditions

Patient-specific (PS) flow simulations using computational fluid dynamics (CFD) have 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 have by some been defined as simulations using measured values of the arteries under consideration along with the 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 (BCs) are determined from other, nonpatient 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 the 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 there are advantages, there are several uncertainties and potential sources of errors using CFD, one being setting appropriate BCs. 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 the mechanical properties of the arterial walls. In addition to water, the blood contains components of different sizes ranging from nm (ions and small molecules) to tens of μm (cells). The red blood cells (RBCs) 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 characteristic 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 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 and 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 are associated with uncertainties and measurement errors. Assessment of the result's sensitivity and in particular the parameters used for evaluating results to uncertainties and errors is crucial if simulation data are to be used for quantitative evaluation.

Different assessments of the effects of boundary data on certain parameters are found in the 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 BCs on the thoracic aorta of a healthy subject reported that idealized BCs may lead to misleading results.

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 diameter 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 the 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.

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 the 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 BCs, quantified by the difference between simulations using the time-averaged and pulsatile BCs. Although time-averaged BCs 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 BCs. 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 the ascending aorta. The volume flow rate as a function of time was 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 nonuniform.

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 BCs. The flow in the thoracic aorta is characterized in terms of commonly used 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.

In this study, the blood flow in the human aorta is considered. The blood is modeled 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 are as follows:

where u_i is the velocity component in the ith direction, p, ρ, and ν are the pressure, density, and kinematic viscosity of the blood, respectively. Here, blood viscosity is assumed to obey the Quemada model $ν=ν(∂ui/∂xj,αRBC)$⁠, where α_RBC is the concentration of the RBCs. In the results presented herein, α_RBC is assumed to be constant, having a value of 0.45. van Wyk et al. [32] shows the details of the Quemada model.

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

• A change in flow variables on any boundary point propagates into the entire domain

• The maximum principle implies that the change in pressure BC will be the largest at the boundary and propagates into the domain having smaller amplitude, depending on the used BC for Eq. (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 BCs 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.

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 et al. [18] and Morbiducci et al. [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 of conditions are commonly used in the literature:

• flow rate based [19,21],

• pressure or pressure gradient based,

• the 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 are 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 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 BCs, such as assumed spatial distribution of the velocity vector with a given flow rate. Simplified BCs 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 Eq. (1) are allowed. See the Appendix for an extended discussion.

The three model coefficients (R₁, R₂, and C) are determined by fitting the model parameters to available data, see Table 6 (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 6 satisfy this criterion and consistency of the results with the underlying assumptions needed in order to make the computed results admissible.

The aortic geometry used in this study was derived from a computed tomography angiography (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 [34]. 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 Fig. 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 the level of the diaphragm was left without extension.

Table 2 shows the references providing the inlet and outlet BCs investigated. The temporal variation of the different inlet flow rates is depicted in Fig. 2. It needs to be stressed that not all profiles were measured aortic flows [19]. However, the underlying reasoning for including these 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 BCs were investigated (Table 2). The three-element circuit Windkessel model parameters used are provided in Table (Appendix) and according to the respective study found in the literature.

Considering Fig. 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 BCs 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 BCs allows for assessing the sensitivity of the different outlet BCs. Moreover, it can be noted that the flow rate as presented by Madhavan and Kemmerling [21] is a smoothed version of the profile of Fuster et al. [35] with a minimum of zero flow rate during diastole. Furthermore, as pointed out by Benim et al. [19], the inlet profile marked as BEN is in fact flow rate profile into the Carotid artery (after the thesis of Ku [37]), which is clearly noticed as it has a weaker systolic phase and a considerable higher flow during diastole as compared to the other profiles.

Altogether, nine different inlet profiles (Fig. 2) and seven different outlet BCs (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 postprocessing 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).

In order to assess the accuracy of the numerical simulations, each of the different geometries was studied using a set of three sequentially refined grids. The grids for the thoracic aorta consisted of 2.5, 5, and 10 × 10⁶ 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).

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 L_u = |U_CL|/|dU_CL/dξ_CL|, where U_CL and ξ_CL represent the axial velocity and the distance along the centerline, 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 |dU_CL/dξ_CL| decreases, implying a cocurrent 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 BCs 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(10⁴) or larger.

A general impression of the time-dependent flow in geometries G1 and G2 is provided as animations of the streamwise velocity profiles (supplementary material available in the Supplemental Materials on the ASME Digital Collection). Figure 3 depicts the flow in the aorta using the inlet conditions of HOPE and outlet condition of Ben. Figure 3(a) shows four time instances of the flow for G1 at midsystole, early diastole, mid-diastole, and start of systole. Figure 3(b) 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 (nonuniformities) in the streamwise velocity. Moreover, intermittent flow separation in the inner/lesser curvature of the aortic arch near peak systole and retrograde flow in the early diastole are also observed. The formation of helical motion is observed as the region of large axial velocity changes with time. The corresponding animations (supplementary material available in the Supplemental Materials on the ASME Digital Collection) 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 counterclockwise 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 the direction between late systole and early diastole, e.g., Figs. 4(b) and 4(c), 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).

To assess the effect of change in inlet geometry and sensitivity to BCs, the following commonly used WSS related parameters were used: TAWSS, the OSI, and the relative residence time (RRT) [38,39].

The effect of the aortic sinus with the BEN inlet and Ben outlet flow rates was 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. Tables 3 and 4 display in more detail the (spatial) mean values of these parameters, where Table 2 contains combinations of inlet- and outlet-BCs for G1. It should be noted that although the two outlet BCs 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 Figs. 5 and 6 and the corresponding animations in the supplementary material available in the Supplemental Materials on the ASME Digital Collection.

Typical distribution of OSI for a selection of different cases is depicted in Figs. 5 and 6. Figures 5(a) and 6(a) relate to different inlet flow rate profiles for a given set of outlet BCs for geometry G1 and G2, respectively. The corresponding cases for fixed inlet condition but having different outlet BCs are shown in Figs. 5(b) and 6(b). 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 BCs influence depends strongly on the deviation from the ideal conditions assumed by each BC.

The strength of the WSS in terms of peak and root-mean-square values depends rather strongly on the BCs; see Table 5 for three of the inlet flow rates combined with six possible outlet BCs. 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 BCs 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 Frydrychowicz et al. [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.

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 converged. For others, the simulations entered a limit-cycle or diverge. The latter cases occur when the underlying assumptions of the 0D models (see the  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 [37]. Given the temporal flow rate and using the three-element circuit along with Eq. (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 variations of pressure predicted by Eq. (2) using the different model parameters (Table 6) 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 Refs. [22,23], and [33].

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 (Fig. 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.

The sensitivity of the WSS-related parameters to the seven different outlet BCs 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 BCs. The appendix of the paper details these assumptions for several common outlet BCs. 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 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 are 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 BCs may be challenging from computational point of view when flow reversal occurs close to an outlet BCs. Such difficulties were encountered for certain combinations of inlet flow rate profiles and certain three-element Windkessel model parameter setups. 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 BCs, for a given set of outflow BCs, and the effects of different outlet BCs 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, available in the Supplemental Materials on the ASME Digital Collection, 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, BCs 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 BCs 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.

Due to the large range of and quick changes in cardiac output, simulations for clinical purposes have to cover the entire range of BCs, in contrast to current, single shot (setup of geometry, BCs, rheological and wall models) in order to give truly patient-specific simulations.

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.

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).

u_i =

velocity

p =

pressure

ρ =

density

ν =

kinematic viscosity

α_RBC =

concentration of red blood cells (hematocrit)

4D-MRI =

four dimensional magnetic resonance imaging

BC =

boundary condition

BPM =

beats per minute

CFD =

computational fluid dynamics

CTA =

computed tomography angiography

LPM =

liters per minute

MRI =

magnetic resonance imaging

OSI =

oscillatory shear index

PC-MRI =

phase contrast magnetic resonance imaging

PS =

patient specific

RBC =

red blood cell

RRT =

relative residence time

TAWSS =

time averaged wall shear stress

US =

ultrasound

WSS =

wall shear stress

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, Eq. (2b), the momentum equations, require BCs on all boundary points. It can be shown that if the flow direction is laminar and “unidirected” 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 Eq. (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 one-dimensional (1D) and 0D approximations. In the following, a short description of the different approximations that can applied to set appropriate outlet BCs on arteries are provided.

More systematic BCs 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 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 [40,41] 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.

where all $∂∂xμ∂∂x$ 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;

The simplified equations still needs BCs on all boundaries except on the outlet boundary (x = X_out).

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 necessarily a constant), u and dP/dx (Eqs. (A1) and (A4) can be computed at plane x = x_i given that u, v, w, and p are known at a the spanwise plane (x = x_i−1) considered. Once u and dP/dx are computed at x = x_i, the two-dimensional (Navier–Stokes) equation for v, w, and p at the same spanwise plane (x = x_i) can be solved. Repeating this procedure enables integration of the three-dimensional 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:

• The length scale in the x-direction is (much) larger than in the y- and z-directions and

• The pressure p is a weak function of y and z and the gradients are small as compared to the size of dP/dx.

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

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 [42,43]. Hence, from Eqs. (A1) and (A5), the following is obtained:

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 [44,45], although for the boundary model above, the compatibility of the assumptions with the underlying assumptions need to be assessed.

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 Eq. (A8) detail the flow rate on the outlet boundary.

Assuming dP/dx = G*Q(t − t₀), where t₀ 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 t₀ is given.

where δ(x − x_BC) is the Dirac delta function, and F is the force acting on the boundary point x_BC. 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) = qe^iωt then the time-derivative of Q and Q itself (in Eq. (A9)) 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 − p₀) + β dP/dt where p₀ and β are model parameters. This assumption allows for satisfying the simplified 1D equation for time-dependent flows with Q = Q(t).

where R_p, C, and R_c are the distal resistance, capacitance and outlet resistance, respectively. Altogether, the role of the Windkessel model is to determine the phase angle between the pressure and the flow rate. This function, the phase angle, is known by measurement and can be implemented directly in the 1D, 0.5D, or the 0D models.