SPECIAL ISSUE PAPER - NUMERICAL METHODS AND APPLICATIONS OF MULTI-PHYSICS IN BIOMECHANICAL MODELING A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models

SUMMARY We present a systematic comparison of computational hemodynamics in arteries between a one-dimensional (1-D) and a three-dimensional (3-D) formulation with deformable vessel walls. The simulations were performed using a series of idealized compliant arterial models representing the common carotid artery, thoracic aorta, aortic bifurcation, and full aorta from the arch to the iliac bifurcation. The formulations share identical inﬂow and outﬂow boundary conditions and have compatible material laws. We also present an iterative algorithm to select the parameters for the outﬂow boundary conditions by using the 1-D theory to achieve a desired systolic and diastolic pressure at a particular vessel. This 1-D/3-D framework can be used to efﬁ- ciently determine material and boundary condition parameters for 3-D subject-speciﬁc arterial models with deformable vessel walls. Finally, we explore the impact of different anatomical features and hemodynamic conditions on the numerical predictions. The results show good agreement between the two formulations, especially during the diastolic phase of the cycle. © 2013 The Authors. International Journal for Numerical Methods in Biomedical Engineering published by John Wiley & Sons, Ltd.


INTRODUCTION
One-dimensional (1-D) and three-dimensional (3-D) formulations have been used extensively to simulate arterial hemodynamics. Landmark contributions in 1-D modeling include the works of Hughes and Lubliner [1], Stergiopulos et al. [2,3], Olufsen et al. [4], Formaggia et al. [5], Sherwin et al. [6], Bessems et al. [7], Mynard and Nithiarasu [8], Alastruey et al. [9], and Müller and Toro [10]. Key contributions to 3-D blood flow modeling in deformable vessels include the works of Perktold et al. [11], Taylor et al. [12], Quarteroni et al. [13], Steinman et al. [14], Cebral et al. [15], Gerbeau et al. [16], and Figueroa et al. [17]. 1-D methods have been used to improve our theoretical understanding of hemodynamics, in particular to study the mechanisms underlying pulse wave propagation [18,19] and also clinically in applications such as wave intensity analysis [20]. Accurate predictions can be made when the flow is predominantly unidirectional and there are no sudden changes in the cross-sectional area. However, 1-D models require the introduction of additional empirical laws to account for recirculation and pressure losses in the presence of vessel curvatures, stenoses, aneurysms, etc. [2,21,22]. These geometric complexities are intrinsically captured with 3-D models, which can provide localized hemodynamic quantities such as wall shear stress, particle residence time, etc. Additionally, 3-D modeling enables the use of structurally motivated biaxial constitutive laws [23] and circumferentially varying mechanical properties for the arterial wall [24] and the simulation of complex processes such as the interactions between the arterial wall and medical devices [25,26]. Nevertheless, 1-D models typically contain far fewer DOFs in comparison with 3-D models (on the order of 10 3 vs. 10 6 ), and the simulations can be executed in a matter of minutes on a personal laptop computer.
Previous works have compared 1-D models and 3-D models with rigid walls in the context of cerebral arterial flow [27,28] and with deformable vessels in the aorta under steady flow [29]. However, a systematic comparison between 1-D and 3-D models in a variety of deformable arterial configurations was still missing in the literature. In general, a good agreement between the two modeling techniques is expected when the flow is mostly unidirectional and there are no sudden changes in cross-sectional area, whereas larger differences are expected in more complex configurations, such as in curved vessels and areas with highly dynamic flows. For a proper comparison to be made, a framework containing the 1-D and 3-D formulations that share equivalent boundary conditions and material laws must be developed. In general, this framework can be used to quickly estimate boundary condition parameters and distributions of material coefficients for an extensive multi-branched 3-D model (e.g., a full-body scale arterial network model [30]), a task that is more time-consuming in the isolated 3-D setting. The purpose of this paper is twofold: 1. To perform a systematic comparison between cross-sectionally averaged hemodynamics (i.e., average flow, average pressure, and radial deformation) obtained with two specific implementations of the 1-D and 3-D formulations in a series of four increasingly complex idealized geometries representing the common carotid artery (CCA), the thoracic aorta, the aorto-illiac bifurcation, and the full aorta up to the first generation of branches. This simplified set of geometries, as opposed to patient-specific geometries, allows for a more direct comparison between the 1-D and 3-D schemes. We investigated the differences in the hemodynamic predictions introduced by increased Reynolds number, the presence of tapering, curvature, and bifurcation angles. To ensure a consistent comparison between formulations, we used equivalent constitutive laws in the 1-D and 3-D formulations, imposed identical flow waveforms and velocity profiles at the inlets, and coupled the same threeelement Windkessel models at the outlets. In these comparisons, the 1-D geometric parameters were obtained directly using dimensions taken from the idealized 3-D geometries ( Figure 1). 2. To develop a computational framework where the 1-D formulation is used to quickly determine boundary condition and material parameters to match patient data such as flow and pressure waveforms and localized measurements of distensibility. Once evaluated, the parameters are fed directly to the 3-D model to perform localized studies. Our aim is to exploit the advantages of both schemes: a computationally efficient 1-D model combined with a full 3-D model sharing identical boundary condition and constitutive laws. This overlapping 3-D/1-D approach can potentially accelerate the solution turn-around time of complex 3-D models, therefore improving their clinical applicability, and differs from previous efforts [13,[31][32][33], where 3-D and 1-D models were coupled to represent spatially distinct parts of the arterial tree.
The structure of this paper is as follows: we first present the details of the 1-D and 3-D formulations utilized in this work, paying special attention to the description of boundary conditions and material laws that are consistent between the two approaches. We then compare the numerical predictions of both methods in a series of idealized geometries. Finally, we describe and discuss the similarities and differences in the results between the two formulations and finish with conclusions and future work.

One-dimensional formulation
Conservation of mass and momentum applied to a 1-D impermeable and deformable tubular control volume of an incompressible Newtonian fluid yields [6,19,34] : where x is the axial coordinate along the vessel, t is the time, A.x, t/ is the cross-sectional area of the lumen, U.x, t/ is the axial blood flow velocity averaged over the cross-section, P .x, t/ is the blood pressure averaged over the cross-section, f is the density of blood assumed to be constant, and f .x, t/ is the frictional force per unit length. The momentum correction factor in the convection acceleration term of Equation (1) was assumed to be equal to one, following the work of Stergiopulos et al. [2]. Equation (1) can also be derived by integrating the incompressible Navier-Stokes equations over a generic cross section of a cylindrical domain [1,[35][36][37]. In this work, the velocity profile is assumed to be constant in shape and axisymmetric. The axial velocity (u) is assumed to have the shape where r.x, t/ is the lumen radius, is the radial coordinate, and is the polynomial order. Following [35], D 9 provides a good compromise fit to experimental data obtained at different points in the cardiac cycle.
Integration of the Navier-Stokes equations of an incompressible Newtonian fluid for axisymmetric vessels yields f .x, t/ D 2 r @u @ j Dr [35]. For the velocity profile given by Equation (2), we have f D 2 . An explicit algebraic relationship between P and A (or tube law) is also required to close Equation (1) and account for the fluid-structure interaction of the problem. In general, we have P D A.A.x, t/I x, t/, where the function A depends on the model used to describe the dynamics of the arterial wall. We solved the system of equations in (1) with the elastic tube law described in Section 2.3 using a DG scheme with a spectral/hp spatial discretization and a second-order Adams-Bashforth time discretization [9,19]. The initial conditions are .A.x, 0/, U.x, 0/, P .x, 0// D .A 0 .x/, 0, 0/, where A 0 is the initial area that yields A d at P D P d , with P d as the diastolic pressure. In all our 1-D simulations, we implemented boundary conditions and solved matching conditions at bifurcations by taking into account the correct propagation of the characteristic information and neglecting energy losses [19,38].

Three-dimensional formulation
The coupled-momentum method [17] is used to formulate the 3-D fluid-structure interaction problem. Here, the DOFs for the vessel wall (displacements u) are described as a function of the fluid velocities at the fluid-wall interface, †, using an 'enhanced' membrane formulation. Under the assumption of a fixed fluid domain, i.e. linearized kinematics of the vessel wall, the strong form is defined in an Eulerian configuration: where f denotes the fluid domain, v is the fluid velocity, p is the pressure, b is a body force (here, assumed to be zero), and f D .r v C .r v/ T / is the viscous stress tensor for a Newtonian fluid. The arterial wall is treated as a thin membrane with thickness h and density s where v s is the solid velocity, s is the membrane stress tensor, and n s is the outward normal at the fluid-wall interface. The value of s was 0.001 g/mm 3 for all models.
To account for the mechanical forces exerted by the external tissues on the arterial walls, an additional term, f support D .k s uCc s v/, was included in Equation (4), which approximates the mechanical behavior of the external tissue [39], where the parameters k s and c s are the stiffness and damping coefficients. This additional force can be used to eliminate spurious and nonphysiological oscillations in specific cases where the geometry is elongated, and the vessel wall experiences asymmetric loads. For the simulations presented in this paper, unless explicitly stated, the values of k s and c s are set to zero.
With the assumption that the vessel wall thickness (h) is small and with appropriate choices for the solution and test function spaces, S, P, and W, the strong form gives rise to the variational equation: where w 2 W and q 2 P are test functions. in is a Dirichlet boundary where the test functions w vanish and where the fluid velocities are prescribed, that is, at the inlet. The traction term h.v, p/ on the outlet boundary out is given as a function of v and p depending on the reduced-order model (i.e., three-element Windkessel) chosen to represent the downstream vasculature, as described by Vignon-Clementel et al. [40,41]. The augmented Lagrangian scheme described by Kim et al. [42] was used to stabilize the outflow velocities in the presence of flow reversal at the outflow boundaries.
To discretize and solve Equation (5), we employed a stabilized semi-discrete FEM on the basis of the work of Brooks and Hughes [43], Franca and Frey [44], Taylor et al. [12], and Whiting and Jansen [45], using equal-order interpolation (P1/P1 elements) for the velocity and pressure fields. The generalized˛-method [46,47] was used to integrate the system of ordinary differential equations resulting from the finite element discretization.

Material laws
The arterial wall was modeled as a thin, incompressible, homogeneous, isotropic, linear elastic membrane characterized by an elastic modulus E, a Poisson's ratio D 0.5, and a thickness h.
In the 1-D model, the arterial wall is assumed to deform axisymmetrically, each cross-section independently of the others, so that the relationship between circumferential hoop stress, T Â , and radial displacement is where r d is the radius at diastolic pressure (P d ). Applying Laplace's law, T Â D .P P d /r= h, and assuming that 1=r can be approximated by 1=r d , we obtain the tube law [9] where A d .x/ is the luminal area at diastolic pressure. In Equation (7), both E and h are functions of the position x. The pulse wave speed c.x, t/ is related to A through The initial area A 0 is calculated by replacing P D 0 and A D A 0 in Equation (7), which leads to In the 3-D formulation, no assumptions regarding axisymmetry were made. The enhanced membrane stress tensor s is given as a function of a tensor Q K of material parameters, a tensor Q P describing the prestress of the wall, and the infinitesimal strain tensor Q ; s D Q K Q , with where k is a transverse shear factor, and u is the displacement vector [48]. In general, the prestress tensor can be specified using a variety of methods [39,48,49]. In the case of linear elasticity, the 'prestress' can be incorporated by subtracting a reference displacement (i.e., at diastole) from the displacement field.

Boundary conditions
The inlet and outlet boundary conditions were chosen to be consistent between the 1-D and 3-D schemes. At the inlet, we prescribed a known flow rate with an axisymmetric velocity profile (Equation (2)). At the outlet of each terminal vessel, we coupled a three-element Windkessel model ( Figure 1). This zero-dimensional (0-D) electrical circuit analog of the downstream vasculature consists of a resistance (R 1 ) connected in series with a parallel combination of a second resistance, (R 2 ), and a compliance, (C ); P out is the pressure at which flow to the microcirculation ceases and is assumed to be zero in all models. The pressure and the flow at an outlet of the 1-D or 3-D domain is related by The numerical implementation of this 0-D model is detailed in [19,38] for the 1-D formulation and [41] for the 3-D formulation.

Iterative procedure for the determination of outflow Windkessel parameters
The parameters of the three-element Windkessel outflow models were calculated as described in the succeeding text. Given a target diastolic (P d ) and systolic (P s ) pressure and flow rate at the inlet (Q in .t /), the initial estimate for the net peripheral resistance (R T ) was calculated as [50] where Q in is the mean flow rate, and P m is the mean blood pressure, assumed uniform throughout the arterial network. We then calculated the resistance R 1 C R 2 at the outlet of each terminal vessel that yields the desired flow distribution and satisfies where M is the number of terminal branches, and j D 1 corresponds to the aortic root. For each individual outlet, the proximal resistance (R 1 ) is assumed to be equal to the characteristic impedance of the upstream 1-D domain, that is, where c d and A d are, respectively, the wave speed and area at diastolic pressure (P d ). This choice of R 1 minimizes the magnitude of the waves reflected at the outlet of the 1-D domain [38]. The total compliance (C T ) was calculated from either: (i) the time constant D 1.79 s of the exponential fall-off of pressure during diastole given in [51] or (ii) using an approximation to C T D dV dP , where V .t/ is the total blood volume contained in the systemic arteries. According to [50], which can be calculated once R T is determined using Equation (13). Alternatively, C T D dV dP can be approximated by [50] where Q max and Q min are the maximum and minimum flow rates at the inlet, and t is the difference between the time at Q max and the time at Q min . We use both Equations (16) and (17) depending on the available input data.
According to [52], we have where C c is the total arterial conduit compliance, C p is the total arterial peripheral compliance, N is the total number of vessels in the 1-D domain, M 1 < N is the number of terminal branches (j D 1 denotes the inlet and is not included in the sum), R 1 , R 2 , and C are the parameters of the three-element Windkessel model (Figure 1), and C 0D is the compliance of each vessel, which is calculated as where L is the length of the vessel. We calculated C p D C T C c and distributed it following the methodology described by Alastruey et al. [38]. More specifically, we have where Q C j is the terminal compliance of each branch distributed in proportion to flow as described by Stergiopulos et al. [2]. We then introduced a correction factor to arrive at the final value of C j : This expression follows from a linear analysis of the 1-D equations in a given arterial network in which each terminal branch is coupled to a three-element Windkessel model [52]. For all of the simulations, the Windkessel compliances and resistances .C j , j D 2, : : : , M / R j 1 andR j 2 , j D 2, : : : , M Á were iteratively calculated to achieve physiologically realistic pressure ranges. To reach a desired pulse pressure (P pulse D P s P d ) and diastolic pressure (P d ) at a particular vessel, we calculated R 0 T and C 0 T given by Equations (13) and (16) or (17) using the iterative formulae where the superscript n is the iteration number of the Windkessel parameter estimation process performed using the 1-D formulation, and P n d and P n pulse are the diastolic and pulse pressure, respectively, at a specific target location in the 1-D model, typically the inlet, at each iteration. Equations (22) and (23) follow from a first-order Taylor expansion of Equations (13) and (17) around the current mean and pulse pressures P n m and P n pulse , respectively, with P n m approximated using the change in diastolic pressure. The total compliance was adjusted by altering the total peripheral compliance C p , because the total conduit compliance C c is a function of the vessel geometry and wall stiffness. This process was iterated using the 1-D model until P n d and P n pulse were smaller than 1% of the target P d and P pulse , respectively. Figure 2 shows the evolution of the systolic, mean and diastolic pressure, net peripheral resistance, and total compliance calculated using the 1-D formulation to match the target systolic and diastolic pressures for the baseline aorta model. The final values of the Windkessel compliances and resistances were used in the 3-D counterparts of the 1-D models.
Other methods have been proposed in the literature to estimate the parameters of the outflow boundary conditions. A root-finding method is described by Spilker and Taylor [53] in the context of 3-D models with compliant arterial walls. Devault et al. proposed a Kalman-filter-based methodology in a 1-D model of the circle of Willis [54].

Error calculations
The numerical solutions of pressure P and volumetric flow Q were compared between the 1-D and 3-D models by using the following relative error metrics [9]: where N t is the number of time points where the comparison is made (typically around several thousand points over a single cardiac cycle depending on the time step size), P 1D i and Q 1D i are the pressure and flow results at each time point i from the 1-D simulation at a single spatial location, and P 3D i and Q 3D i are the cross-sectional averaged pressure and flow at each time point i from the 3-D model at a single cross section perpendicular to the vessel centerline. " P ,avg and " Q,avg are the average relative errors for pressure and flow, respectively, over one cardiac cycle, and " P ,max and " Q,max are the maximum relative error in pressure and flow. " P ,sys and " Q,sys are the errors for systolic pressure and flow, and " P ,dias and " Q,dias are the errors for diastolic pressure and flow. In order to avoid division by small values of flow, we normalized the flow error metrics by the maximum value of flow over the cardiac cycle, max i Q 3D i . The error metrics were calculated over a single cardiac cycle once both numerical solutions achieved periodic behavior.

RESULTS
We investigated the differences in the numerical prediction of flow and pressure between the 1-D and 3-D formulations (where flow and pressure refer to cross-sectional averages in planes perpendicular to the vessel centerline) in a series of test cases by using idealized geometries. The physical dimensions, vessel wall properties, and inflow and outflow boundary conditions were made to be consistent between the 1-D models and their 3-D counterparts. Mesh independence studies were undertaken to ensure that the results in the final meshes were not affected by inadequate mesh resolution. The geometry, inflow rate, and mechanical properties of the blood were taken from [17], the Young's modulus E from [65], and the pressures from [55, p. 343]. The parameters of the RCR Windkessel model were calculated as described in the text. The resulting wave speed at mean pressure is c m D 6.74 m s 1 .

Baseline common carotid
We considered a straight cylindrical vessel with representative dimensions of the common carotid artery (CCA). The initial total peripheral resistance (R T ) and compliance (C T ) were calculated from Equations (13) and (17), respectively, using a reference value of mean blood pressure (P m ) measured in a 23-year-old human [55, p. 343] and a reference inflow waveform [17]. The final values of R T and C T were then computed as described in Section 2.5, requiring a single iteration to achieve convergence to the target pressures at the outlet. The parameters of the CCA model are given in Table I. The 1-D simulation was run using six elements with a quadrature and polynomial order of 10 and a time step of 0.1 ms. The velocity profile order was D 2. The initial area A 0 D 0.22 cm 2 that yields the reference diastolic area (A d ) at P D P d was calculated using Equation (9). The 3-D geometry was constructed with the same dimensions as in the 1-D segment. The velocity on the inlet boundary was prescribed using Equation (2) with D 2 and with the same time-averaged flow rate as at the inlet of the 1-D model. The 3-D simulation was run with a mesh containing 792,559 linear tetrahedral elements (145,394 nodes) and a time-step of 0.2 ms. Figure 3 shows the numerical predictions of flow rate, pressure, changes in luminal radius, and velocity profiles at several sites in the 1-D and 3-D models. The flow waveforms exhibit the typical attenuation observed in vivo between the inlet and outlet. We also include the pressure gradient (given as the difference between inlet and outlet pressures) predicted by the two numerical schemes. The results show excellent agreement between the two models, with average relative errors smaller than 1%. However, even though the flow waveforms in 1-D and 3-D models are virtually identical, differences can be observed in the velocity profiles. These differences result from the linearized kinematics assumption (i.e., fixed computational domain) of the 3-D formulation, in contrast to the 1-D model, which accounts for changes in the cross-sectional area. In all the examples considered in this manuscript, the 1-D velocity profiles were plotted in the reference configuration given by the nominal (diastolic) diameter of the vessel to provide a qualitative comparison between the profile shapes.

Baseline aorta
We considered a straight cylinder with diameter representative of the average diameter of the thoracic aorta. Here, the initial peripheral resistance was calculated using systolic and diastolic pressures from the study by Simon et al. [51], and the total compliance C T was calculated from the time constant D 1.79 s [51] of the exponential fall-off of pressure during diastole using Equation (16). The final values of R T and C T were obtained after seven iterations to reach convergence to the target pressures at the outlet. The parameters of the baseline aortic model are given in Table II.  The inflow rate was taken from [50] and the pressures from [51]. The parameters of the RCR windkessel model were calculated as described in the text. The resulting wave speed at mean pressure is c m D 5.17 m s 1 .
For the 1-D model, the simulation was run using 12 elements with a quadrature and polynomial order of 5 and a time step of 0.1 ms. The initial area was calculated using Equation (9) to be A 0 D 3.06 cm 2 . The polynomial order of the velocity profile was chosen to be D 9 on the basis of [35]. For the 3-D model, the velocity boundary condition at the inlet was prescribed with a profile of order D 9. The 3-D simulation was run using a mesh containing 1,480,048 linear tetrahedral elements (261,912 nodes) and a time step of 0.2 ms. Figure 4 shows the numerical predictions of flow rate, pressure, changes in luminal radius, and velocity profiles at several sites in the 1-D and 3-D models. Figure 4 also contains a plot of the pressure gradient between the inlet and outlet. The overall agreement between the two formulations is still reasonably good, particularly during diastole, with average relative errors smaller than 2%. The largest differences in pressure are seen at the inlet during acceleration, and the largest differences in flow are observed during peak systolic flow at the midpoint and outlet locations. The velocity contours at deceleration (t D 0.33 s) are markedly different: although the 1-D solution preserves the fixed ninth-order velocity profile, the 3-D results show near-wall flow reversal, resulting in a Womersley-like profile. Furthermore, the axisymmetry of the 3-D profile breaks down in the second half of the vessel during peak systole (t D 0.15 s). The pressure gradient plot shows larger absolute differences than in the carotid geometry and a slight phase-lag between the 3-D and 1-D predictions.

Low flow and larger diameter aorta
To study the impact of Reynolds number and flow inertia, we considered two more cases in addition to the baseline straight aorta model. In the first case, the flow prescribed at the inlet was reduced  by a factor of 9.54 to match the peak Reynolds number to that of the CCA model (Re D 748).
In the second case, the diameter of the cylindrical domain was increased to 1.5 cm from 1.2 cm, whereas the flow rate and vessel length remained unchanged (Re D 5713 at peak systole). Figures 5  and 6 show the flow rate, pressure, change in luminal radius, and pressure gradient at several sites in the 1-D and 3-D models. Although reducing the Reynolds number via decreasing the total flow did not decrease the error in the inlet pressures during acceleration ( Figure 5), increasing the vessel diameter without altering the total flow ( Figure 6) did bring the 3-D and 1-D predictions closer.

Tapered carotid
We considered a linearly tapered cylinder with the same length as the previous CCA geometry. On the basis of reference values for the degree of tapering of the left CCA in [56], the diameter was set to 8 mm at the inlet and 4 mm at the outlet. Although R T remained unchanged from the non-tapered CCA model, the resistance R 1 (Equation (15)) was recalculated with the values of c d and A d at the location of the outlet, giving a value of R 1 D 6.8548 10 8 Pa s m 3 and R 2 D 1.4330 10 9 Pa s m 3 . The remaining parameters were unchanged from those given in Table I. For the 1-D model, the velocity profile was specified with D 2, and the simulation was carried out using six elements with a quadrature and polynomial order of 10 and a time step of 0.1 ms. For the 3-D model, the velocity boundary condition at the inlet was prescribed using a profile of order D 2, and the simulation was run using a mesh containing 840,899 linear tetrahedral elements (153,539 nodes) and a time step of 0.2 ms. Figure 7 shows the flow rate, pressure, pressure gradient, velocity profiles, and change in luminal radius at several sites in the 1-D and 3-D models. Slightly larger differences in 3-D/1-D predicted pressure are observed at the inlet and in the center of the vessel. Compared with the non-tapered CCA geometry, there is also a significant increase in pressure gradient discrepancy.

Tapered aorta
We considered a linearly tapered cylinder with length identical to that of the baseline non-tapered aorta model. On the basis of typical reported aortic dimensions [56], the diameter was set to 3 cm at the inlet and 2 cm at the outlet. R T remained unchanged from the baseline non-tapered aorta model, but the resistance R 1 was recalculated with the values of c d and A d at the location of the outlet. We obtained R 1 D 1.8503 10 7 Pa s m 3 and R 2 D 1.0492 10 8 Pa s m 3 . The remaining parameters were unchanged from those in Table II.
For the 1-D model, the polynomial order of the velocity profile was chosen to be D 9, and the 1-D simulation was run using 12 elements with a quadrature and polynomial order of 5 and a time step of 0.1 ms. For the 3-D model, the velocity boundary condition at the inlet boundary was prescribed using a profile of order D 9, and the simulation was run using a mesh containing 1,663,772 linear tetrahedral elements (293,667 nodes) and a time step of 0.2 ms. Figure 8 shows the flow rate, pressure, changes in luminal radius, velocity profile, and pressure gradient at several sites in the 1-D and 3-D models. Results show larger differences in pressure and flow waveforms at the midpoint and inlet locations, respectively, compared with those found in the tapered carotid case. The differences between the velocity contours are also more noticeable, with the 3-D profile displaying near-wall flow-reversal patterns.

Aorta with curvatures
We considered two additional aortic models using the same inflow and outflow boundary conditions as the baseline aorta model with the purpose of studying the impact of vessel curvature on cross-sectional average pressure and velocity. In the first case, curvature was introduced in a single plane (sagittal) to represent an idealized aortic arch. In the second case, curvatures were introduced in three planes (sagittal, coronal, and transverse) as described by Yearwood and Chandran [57]. Lengths and diameters were kept unchanged relative to the baseline aorta case. Because the 1-D model is independent of curvature, the 1-D solution is identical to that obtained in Section 3.2. The boundary conditions and material parameters were identical to those in the baseline aorta model (Table II).
Because of the centrifugal loading experienced by the vessel wall in the curved portion of the 3-D geometry, a small value for the external damping coefficient (c s D 300 Pa s m 1 ) was chosen to reduce nonphysiological oscillatory modes. The external damping smooths the pressure and flow waveforms predominantly where high frequencies are dominant, with distal waveforms being smoother than proximal ones. These results are in agreement with the linear analysis of wall viscoelasticity described in [50]. The simulation of the 3-D model containing the single plane of curvature was run using a mesh containing 1,691,525 linear tetrahedral elements (297,472 nodes) and a time step of 0.2 ms. The model containing three planes of curvature was simulated with a mesh containing 1,688,467 linear tetrahedral elements (297,083 nodes) and a time step of 0.2 ms. Figures 9 and 10 show the flow rate, pressure, changes in luminal radius, and pressure gradient for the single-curvature and three-curvature cases, respectively, and compares them with the corresponding 1-D predictions, which are the same as those shown in Figure 4. In-plane and out-of-plane velocity profiles at three cross sections in the 3-D models illustrate the skewing of the velocity profile toward the outer wall during systole and the presence of backflow along the inner wall during diastole. Nevertheless, the relative errors in flow, pressure, and pressure gradient are similar to those observed in the baseline aorta case.

Aortic bifurcation
Here, we considered an idealized model of the aortic bifurcation, containing a single parent segment, representing the abdominal aorta, and two branch segments representing the iliac arteries. The lengths and diameters are listed in Table III. The initial net peripheral resistance and total compliance were calculated from Equations (13) and (17) by using the systolic and diastolic pressures from [51]. The final values of R T and C T were obtained after three iterations to reach convergence to the target pressures at the inlet. The peripheral resistance R 1 CR 2 at the outlet of the iliac arteries was calculated assuming equal flow distribution. For the 1-D model, the simulation was run using 12 elements with a quadrature and polynomial order of 5 and a time step of 0.1 ms. The initial area was calculated using Equation (9) to be A 0 D 1.8062 cm 2 in the abdominal aorta and A 0 D 0.9479 cm 2 in the iliac branches. The velocity profile was chosen to have a polynomial order D 9. For the 3-D model, the geometry was constructed from three cylinders with the same diameters and lengths as in the 1-D network. The angle between the two iliac branches was set to 47.9 degrees [58]. The velocity boundary condition at the inlet was prescribed with a profile of order D 9. The simulation was run using a mesh containing 1,799,117 linear tetrahedral elements (321,651 nodes) and a time step of 0.2 ms. Figure 11 shows the flow rate, pressure, velocity, and changes in luminal radius at several sites in the 1-D and 3-D models. The results show excellent agreement between the two models (average relative errors smaller than 2%), with similar relative errors to the baseline common carotid artery case.

Full aorta
We considered an idealized geometry representing the aorta and the first generation of main branches from just above the sinuses to the aortic bifurcation, neglecting the coronary and intercostal vessels. The curvature and angulation from the sagittal plane in the aortic arch were based on published measurements from human cadavers [57,59]. The dimensions and spacing between branch vessels were based on those reported in [3]. The parameters are summarized in Table IV, and the geometry is illustrated in Figure 12. The spatially varying elastic moduli (E) were calculated from the pulse wave velocity, c, given by the empirical relationship [56] c D a 2 =.2r d / b 2 , where c is given in m/s, r d is the radius at diastolic pressure expressed in mm, a 2 D 13.3 and b 2 D 0.3. For each vessel segment, once c was determined, E was calculated using Equation (8) with A D A d . The spatially varying wall thickness (h) was chosen to be 10% of r d [55] at the inlet of each segment. The total peripheral resistance R T D 1.1583 10 8 Pa s m 3 was obtained from Equation (13) with Q in D 6.17 l min 1 , P s D 16.8 kPa, and P d D 9.5 kPa. At the outlet of each terminal branch, the resistance R 1 C R 2 was calculated to yield the flow distribution given in [2] (Table V). Each resistance R 1 follows from Equation (15) with c d and r d calculated at the outlet of the terminal branch. The total compliance C T was calculated from the time constant D 1.79 s [51] of the exponential decline in pressure during diastole using Equation (16); C T D =R T D 1.5453 10 8 m 3 Pa 1 . The compliance C of each terminal RCR Windkessel model was then calculated as described in Section 2.5. The final values of R T and C T were obtained after eight iterations to reach convergence to the target pressures at the outlet of the abdominal aorta.
For the 1-D model, the simulation was run using 46 elements with a quadrature and polynomial order of 5 in each vessel and a time step of 0.05 ms. The initial areas A 0 that yield the diastolic areas A d at P D P d were calculated using Equation (9) (A 0 D 4.7203 cm 2 at the aortic root). The The dimensions are based on data given in [3]. The flow rate is taken from [66] and the pressures from [51]. The parameters of the two RCR models were calculated as described in the text. The resulting wave speed at mean pressure is c m D 6.26 m s 1 in the abdominal aorta and c m D 7.4 m s 1 in both iliac arteries.
velocity profile was chosen with polynomial order D 9. For the 3-D model, the external damping coefficient was assigned a small value (c s D 300 Pa s m 1 ) to reduce nonphysiological oscillations. The simulation was run using a mesh containing 2,554,521 linear tetrahedral elements (475,000 nodes) and a time-step of 0.2 ms. The velocity boundary condition at the inlet was prescribed with a profile of order D 9. Figure 13 shows flow rate and pressure at a number of representative locations in the 1-D and 3-D models. The differences in the pressure and flow predictions occur predominantly in peak systole, and these discrepancies are greater in distal sites than in the proximal aortic arch branches. Additionally, the results from the 3-D model exhibit fewer high-frequency features in the flow and pressure waveforms than in the 1-D model as a consequence of external damping.

Carotid versus aortic geometry: combined impact of Reynolds number, wall strain and velocity profile
In the baseline carotid case, the assumptions of the 1-D formulation, that is, predominantly unidirectional flow, hold well. The relative errors in flow and pressure were small between the 1-D and 3-D carotid models, with " P ,avg and " Q,avg generally less than 0.3%, as shown in Figure 3. In the baseline aorta case, however, the differences between the 1-D and 3-D models were greater (Figure 4), particularly at the inlet, where the 3-D model shows a more pronounced systolic 'shoulder' (at t 0.1 s) in the pressure waveform compared with the 1-D model. Furthermore, the pressure gradient between the inlet and outlet is greater in the 3-D model during peak systole and also exhibits a phase lag between the 1-D and 3-D models, that is, the time at which the pressure gradient from the 3-D model reverses direction during systole is delayed compared with the 1-D model.
Inertial forces play a larger role in the aorta model compared with the carotid model, with the peak Reynolds number in the aorta model being nearly an order of magnitude greater than in the carotid model (Re D 7, 140 vs. Re D 748). To investigate the impact of the Reynolds number in the pressure waveforms, we scaled down, by a factor of 9.54, the aorta model inflow to match the peak Reynolds number of the baseline carotid model. We observed that the average relative error in the pressure waveforms at the midpoint and outlet locations decreased, but the relative errors at the shoulder of the inlet pressure waveform and the phase lag in the pressure gradient were not reduced (compare Figures 4 and 5). To investigate the origin of these discrepancies further, we must reflect on the different treatments that our specific 1-D and 3-D implementations make regarding geometric nonlinearities. Although both formulations consider equivalent linear material laws, the 3-D formulation assumes a linearized kinematics behavior for the vessel wall (i.e., a fixed computational grid is maintained throughout the simulation), whereas the 1-D formulation accounts for changes in cross-sectional area (Equation (7)). Therefore, the fixed-grid assumption of the 3-D formulation introduces errors in the velocity profile that are proportional in magnitude to the vessel wall strain. The aorta experiences both greater flow velocities and larger wall strains ( D r max =r d ) compared with the carotid r in ! r out : diastolic cross-sectional radii at the inlet and outlet of the arterial segment. c in ! c out : wave speed at diastolic pressure at the inlet and outlet of the arterial segment.
setting (i.e., aortic D 12.5% versus carotid D 6%). These two phenomena are responsible for the larger discrepancies between the 1-D and 3-D results during systole in the baseline aorta compared with the carotid example. This idea is supported by the larger diameter aorta model. Despite being a more compliant model than the baseline aorta (because of a larger cross-sectional area), it experiences smaller wall strains. The centerline velocities are also reduced compared with the baseline aorta model. We therefore have a situation where both the Reynolds number and the wall strains are reduced. This ultimately results in a smaller relative error in the shoulder of the inlet pressure waveform ( Figure 6) compared with the baseline aorta model (Figure 4).

Impact of tapering
The presence of tapering introduces errors between the 1-D and 3-D models in the systolic pressure (but not during the diastolic decay) in both carotid and aortic models (Figures 7 and 8). The 3-D model more accurately captures the spatial changes in the velocity profile and therefore the viscous friction losses in the tapered tube. With tapering, we observed greater errors between the 1-D and 3-D model in the pressure gradient between the inlet and outlet, as the 3-D model predicts that a greater pressure gradient is needed to drive flow during peak systole. In the tapered aorta case, inertial and frictional effects are naturally better captured in the 3-D model, resulting in greater errors in the pressure waveform than in the tapered carotid case.

Impact of curvature
The velocity field in the curved 3-D domain is complex and not axisymmetric (Figures 9 and 10). During systole, the region of higher velocity is located near the outer wall of the arch, and during deceleration, significant backflow occurs near the inner wall. As such, the 3-D model is inherently more accurate in capturing the inertial and frictional forces in the curved domain, which leads to differences in the integrated flow and pressure quantities compared with the 1-D model. Hence, we observed that the addition of a single plane of curvature slightly increased the relative errors in the pressure and flow waveforms ( Figure 9) compared with the uncurved baseline aorta case. In particular, the error in the systolic shoulder of the inlet pressure waveform was notably increased (compare Figures 4 and 9). External damping was necessary in the 3-D model to produce flow and pressure waves without spurious oscillations. This damping was not included in the 1-D formulation and contributes to some of the discrepancies observed. The added dissipation attenuated the high-frequency features of the flow and pressure waveforms in the 3-D model, whereas these features are more prominent in the 1-D model, that is, the dicrotic notch and the early-diastolic flow reversal phase. We further observed that the addition of the second and third planes of curvature to the 3-D model did not further introduce significant differences between the 3-D and 1-D predictions ( Figure 10) compared with the single curvature case (Figure 9). In Section 4.5, we investigate the relationship between external damping and curvature in the 3-D setting.

Bifurcation
The 1-D and 3-D models are in good agreement in the aortic bifurcation case. We observed small relative errors (average relative errors smaller than 2%) between the 1-D and 3-D models (Figure 11), which was expected considering the lower flow rate in comparison with the aorta case. Furthermore, changing the bifurcation angle from 47.9 to 27.5 degrees resulted in a negligible change in the relative errors. A similar good agreement was obtained with a geometric model of the abdominal aorta featuring two 90-degree branches representing the renal arteries. These results suggest that the effect of energy losses at the aortic bifurcation on the pressure waveform is minor, as was previously observed using in-vitro data [22].

Full aorta comparison
The full aorta model is essentially a combination of all the previously considered geometric features, including tapering, curvature, branching, and a variety of vessel diameters. As such, we expected to observe larger differences in the prediction of flow and pressure waveforms between the two theories. In general, the largest differences occur in systole, whereas the diastolic predictions are much closer. We can explain the closer agreement between models during diastole using a 0-D Windkessel model for pressure. Starting with the 1-D formulation, a space-independent Windkessel pressure, p w .t /, can be derived by neglecting nonlinearities, flow inertia and viscous dissipation, and assuming that wall compliance and fluid peripheral resistance are the dominant effects [52]. Under these premises, p w .t / is given by where Q in .t / is the flow waveform at the inlet of the ascending aorta, p w .T 0 / is the pressure p w at t D T 0 , T 0 is the time corresponding to the beginning of systole (T 0 D 0 in Figure 14), M 1 < N is the number of terminal branches, and q j out .t / is the outflow in the terminal Segment j . The net peripheral resistance (R T ) is given by Equation (14), and the total arterial compliance (C T ) is given by Equation (18). The parameters of the resistor-capacitor-resistor Windkessel models are R 1 , C , R 2 , and P out (Figure 1). The validity of Equation (28) is also supported by in vivo studies, which show that the aortic pressure waveform is approximately uniform in space during approximately the last two thirds of diastole [60,61]. Figure 14 (left) illustrates that the agreement between 3-D, 1-D, and analytical pressures improves as diastole progresses, indicating that the physics of flow are fundamentally linear and inertiafree during this phase of the cardiac cycle. Conversely, systolic flow is fundamentally nonlinear and advection/inertia-dominated, and therefore, larger differences between the two methods were observed.
The discrepancies observed in systole are partially due to the different treatment of the arterial wall mechanics in the two formulations: although both theories consider equivalent purely elastic material laws, the 3-D case requires a viscous dissipation term provided by the external tissue support model. Without this dissipation, the 3-D flow and pressure waveforms contain spurious, large amplitude high-frequency oscillations resulting from unconstrained rigid-body motion modes that are generated in non-axisymmetric flows. This viscous damping, although necessary in the 3-D setting, attenuated the sharper features of the flow and pressure waveforms, most noticeably in the dicrotic notch and early-diastole flow reversal phases, compared with the predictions seen in the 1-D setting ( Figure 13). It is worth noting that external tissue support can also be included in 1-D modeling, as recently demonstrated by Formaggia et al. [33]. That work focuses on the importance of including external tissue support in the 1-D formulation to reduce reflections at the 1-D/3-D interfaces. Furthermore, although the external tissue support will never be engaged in a 1-D setting to limit spurious wall motion as in the 3-D model, its presence will affect the deformation of the arterial wall because of the additional stiffness. Figure 14 (right) illustrates the impact of external tissue damping by introducing a new model consisting of a modified 'straightened' full aorta geometry. We compared infrarenal flow and pressure waveforms in four different cases: (i) 3-D formulation in original full aorta geometry considering tissue damping; (ii) 3-D formulation in straightened full aorta geometry considering tissue damping; (iii) 3-D formulation in straightened full aorta geometry without tissue damping; and (iv) 1-D formulation. The results in the straightened 3-D aorta models (dashed blue and green lines) demonstrated sharper features than those observed in the 3-D curved aorta model (solid red lines) and, notwithstanding the larger pulse pressures in the 3-D model, resembled more closely the results from the 1-D model (dashed black lines). Furthermore, the exclusion of external tissue damping (green dashed lines) resulted in almost no difference in the pressure and flow waveforms in the straightened model. This behavior suggests that eliminating the aortic curvature reduced the rigid-body modes in the thoracic region and, even though the external viscous damping is still present, it is engaged to a lesser degree than in the original curved aorta model.

Limitations
It is important to note that the comparisons presented in this paper were between two particular implementations of 1-D and 3-D theories with equivalent boundary conditions and material laws. However, a one-to-one comparison of hemodynamics between the two theories is not always possible.
A few examples of fundamental modeling differences that are not easy to overcome are the inability to account for secondary flows and the assumption of a fixed velocity profile in the 1-D formulation; the linearized kinematics approach and the need for outflow boundary condition stabilization and external tissue support in the 3-D formulations.
To facilitate an improved comparison, the 1-D formulation could incorporate space-varying and time-varying velocity profiles [7,56,[62][63][64]. Furthermore, external damping modeled using a viscoelastic tube law [9] and wall inertia [5] could be included in the 1-D model as well. On the other hand, the 3-D formulation would need to incorporate moving domains/meshes [11,16] to enable a more consistent description of cross-sectional area changes between the two methods.

Implications of fundamental modeling differences between 3-D and 1-D theories
An inconsistency between our specific 1-D and 3-D implementations was introduced by assigning zero velocity boundary conditions at the inlet and outlet rings in the 3-D model, effectively clamping the vessels at those locations. Consequently, the conduit compliance (C c ) of the vessels in the 3-D models is slightly reduced relative to the 1-D models; this becomes more significant as the length of the vessel domain decreases. The clamping has a negligible local effect on pressure and flow in the 3-D domain, even though the wall displacements near the inlets and outlets are reduced. A better treatment of the arterial wall boundary conditions in the 3-D formulation can be performed using external tissue support formulations [39] and time-resolved medical-image data to inform both the radial and longitudinal components of the vessel wall motion. Lastly, an important inconsistency between our two modeling approaches is introduced by the need for external damping in the 3-D model in complex and curved geometries, such as the curved aorta and full aorta.

CONCLUSIONS
In this article, we have compared 1-D and 3-D hemodynamics in a series of idealized compliant arterial models of the carotid artery, thoracic aorta, aortic bifurcation, and full aorta. The 1-D and 3-D formulations share common inflow and outflow boundary conditions and have equivalent material laws. We also presented an iterative algorithm to determine the parameters of each outflow boundary condition module to achieve desired systolic and diastolic pressures at a particular vessel. We have demonstrated good agreement between 1-D and 3-D predictions, especially during diastole. The larger differences in systole can be explained by the following: (i) the inability to account for secondary flow features, vessel curvature, and a prescribed velocity profile shape in the 1-D model and (ii) the linearized treatment of the kinematics of the vessel wall and the external tissue support (which introduces viscous damping) in the 3-D model.
Our 1-D/3-D framework demonstrates the use of a 1-D model to determine boundary condition and material parameters that are subsequently fed to a corresponding 3-D model. Indeed, the relatively good agreement between the numerical predictions in the full aorta setting suggests that the 1-D model is a reasonable representation of the 3-D system in terms of the global behavior of the spatially averaged pressure and flow waveforms. It follows that the 1-D model can be used to perform tasks such as: (i) estimation of the parameters of outflow boundary conditions to reproduce clinical measurements and (ii) sensitivity studies under different hemodynamic conditions to gain an understanding of the behavior of the arterial system. This overlapping 1-D/3-D approach will accelerate the selection of model parameters for 3-D subject-specific models. and the NIHR grant above. The 3-D simulations were supported in part by an NSF Extreme Science and Engineering Digital Environment (XSEDE) startup allocation. Lastly, the authors acknowledge Emeritus Professor Kim H. Parker for his critical review of the manuscript, Desmond Dillon-Murphy for his assistance with the generation of the CAD geometry of the full aorta, and Simmetrix, Inc. (http://www.simmetrix.com/) for their MeshSim mesh generation library.