On the periodicity of cardiovascular fluid dynamics simulations
Martin R. Pfaller, Jonathan Pham, Nathan M. Wilson, David W. Parker, Alison L. Marsden
OOn the periodicity of cardiovascular fluid dynamics simulations
Martin R. Pfaller, Jonathan Pham, Nathan M. Wilson, David W. Parker, Alison L. MarsdenFebruary 2, 2021
Abstract
Three-dimensional cardiovascular fluid dynamics simulations typically require computation of sev-eral cardiac cycles before they reach a periodic solution, rendering them computationally expensive.Furthermore, there is currently no standardized method to determined whether a simulation has yetreached that periodic state. In this work, we propose use of the asymptotic error measure to quan-tify the difference between simulation results and their ideal periodic state using lumped-parametermodeling. We further show that initial conditions are crucial in reducing computational time anddevelop an automated framework to generate appropriate initial conditions from a one-dimensionalmodel of blood flow. We demonstrate the performance of our initialization method using six patient-specific models from the Vascular Model Repository. In our examples, our initialization protocolachieves periodic convergence within one or two cardiac cycles, leading to a significant reduction incomputational cost compared to standard methods. All computational tools used in this work areimplemented in the open-source software platform SimVascular. Automatically generated initial con-ditions have the potential to significantly reduce computation time in cardiovascular fluid dynamicssimulations.
Three-dimensional (3D) blood flow simulations are commonly coupled with zero-dimensional (0D) lumpedparameter models, representing the downstream vasculature at the model’s boundary [1, 2, 3, 4]. Theselumped parameter models are analogous to an electric circuit, with resistors and capacitors modeling theviscosity of the blood and the elasticity of the vessel wall, respectively. A popular choice is the three-element Windkessel model, also known as the RCR model [5]. The Windkessel consists of a proximalresistance in series with a parallel distal resistance and capacitance (Figure 1).The RCR boundary condition, like many other outflow boundary conditions, contains a capacitorthat stores blood volume. While these capacitors are “charging”, it typically takes several cardiaccycles with a periodic pulsatile inflow to reach a periodic state. A periodic state is here defined astwo consecutive cardiac cycles yielding results for pressure and flow rate that agree within a giventolerance. Mathematically, this corresponds to the limit cycle of the model. It is essential to extractsimulation results only when they have reached this periodic state, e.g., when comparing them to invivo measurements. However, there are currently no clear guidelines on how to determine whether thisperiodic state has been achieved.Unfortunately, running several cardiac cycles of a 3D simulation is computationally expensive, typi-cally requiring a high-performance computer. The computation time scales linearly with the number ofcardiac cycles since they cannot be run in parallel. A common practice is initializing a simulation withresults from a steady-state solution [5]. Here, a computationally inexpensive simulation with constantinflow is computed first. Its solution is then used as an initial condition for the simulation with a pe-riodic pulsatile inflow condition. However, as we will show in this work, a simulation initialized with asteady-state solution often still requires several cardiac cycles to reach a reasonably periodic state.Our goal is twofold. First, we introduce metrics and tools to quantify whether a simulation hasyet reached a periodic state within a given tolerance. Second, we shorten the computation time of thecomputationally expensive 3D simulation by leveraging reduced-order modeling. It was previously shownthat one-dimensional (1D) models of cardiovascular fluid dynamics could accurately approximate integralquantities of 3D solutions such as velocity and pressure at the outlets [6, 7, 8, 9, 10, 11, 12]. In thiswork, we propose a novel method to initialize a 3D simulation using the results from an inexpensive 1Dsimulation. This framework is fully automated and requires no user interaction.1 a r X i v : . [ c s . C E ] J a n Methods
We begin by revisiting the governing equations of 3D, 1D, and 0D fluid dynamics in Sections 2.1, 2.2,and 2.3. Following theoretical considerations of the convergence of lumped parameter models, we definesuitable error metrics to determine the difference between a 3D simulation and its periodic state inSection 2.4. In Section 2.5, we introduce a tool to determine whether a simulation has yet reachedits periodic state. Finally, we develop a method to initialize a 3D simulation from a 1D simulation inSection 2.6 in order to jumpstart initialization and reduce computational cost.
The dynamics of blood flow in the cardiovascular system is mathematically governed by the incompress-ible 3D Navier-Stokes equations, ρ ( ˙ v + v ⋅ ∇ v ) = −∇ P + µ ∇ v + f , in Ω , (1) ∇ ⋅ v = , in Ω . (2)The first equation in this system represents conservation of momentum for a Newtonain, incompressiblefluid, where v is the velocity of the blood flow, P is pressure, f is a body force, ρ is the density of theblood, and µ is the dynamic viscosity. The second equation in this system represents conservation ofmass. In the computational cardiovascular modeling and simulation context, we typically numericallysolve the 3D Navier-Stokes equations in patient-specific models of vascular anatomies to simulate hemo-dynamics. Simulation results are used to elucidate the relationship between cardiovascular diseases andfluid mechanics, for personalized treatment planning, and to aid the development of novel biomedicaltechnologies. The initial conditions for velocity and pressure are v ( x , t = ) = v ( x ) , P ( x , t = ) = P ( x ) , (3)Boundary conditions that model the portion of the cardiovascular system not captured by the anatomicalmodel must be provided as well. A flow rate, Q , is commonly prescribed at the inlet surfaces of the 3Dvascular model, where the flow rate is computed via integration of the normal velocity over each inletsurface, Q ( x , t ) = ∫ Γ v ⋅ n d Γ = Q in ( t ) , in Γ on . (4)On the other hand, lumped parameter models, also known as 0D models, are commonly used as boundaryconditions at the outlets of the model. These lumped parameter models usually relate the pressure tothe flow rate via parametric differential-algebraic equations, P ( x , t ) = f ( x , t, Q out ( t ) , ˙ Q out ( t ) , φ ) , on Γ out , (5)where φ represents the set of variables parametrizing the differential equation. A discussion of 0D modelsand some commonly used outlet boundary conditions is provided in section 2.3.We generate 3D patient-specific models using SimVascular, an open-source, comprehensive modeling,and simulation software for vascular anatomies ( simvascular.org ) [13]. The models are simulated inour open-source solver svSolver using the Finite Element Method (FEM) ( github.com/SimVascular/svSolver ). It uses liner P1-P1 elements with a streamline upwind Petrov-Galerkin and pressure-stabilizing Petrov-Galerkin formulation (SUPG/PSPG) [14]. The pressure and momentum stabilizationis detailed in [15, 16]. Furthermore, it uses a linear solver with specialized pre-conditioners tailored tohandle large vascular resistances coupled at outflow boundaries [17]. All 3D simulations in this workwere run on Stanford’s Sherlock supercomputing cluster using four 12-core Intel Xeon Gold 5118 CPUs. In contrast to 3D models, one-dimensional (1D) models have only a single spatial dimension, the axialdimension along the centerline of the vessel [18, 6, 19]. Due to this lack of 3D spatial information, 1Dmodels are capable of simulating only bulk flow rate, and cross-sectionally averaged pressure at eachcenterline node of the 1D finite element model. We integrate incompressible 3D Navier-Stokes equations21) over the cross-section while assuming Newtonian fluid properties for the blood and an axisymmetricparabolic flow profile to obtain the governing equations for the 1D model. This process yields˙ S + ∂Q∂z = , (6)˙ Q + ∂∂z Q S + Sρ ∂P∂z = Sf − πν QS + ν ∂ Q∂z , (7)with flow rate Q , cross-sectional area S , pressure P , density ρ , body force f , and kinematic viscosity ν .The coordinate z represents the axial dimension of 1D model. To solve these equations, we also requirea constitutive law to relate the pressure to the cross-sectional area. In this work, we use the constitutiverelationship proposed by Olufsen [20], P ( z, t ) = P ( z ) + Ehr ( z ) ⎛⎜⎝ − ¿``(cid:192) S ( z ) S ( z, t ) ⎞⎟⎠ , Ehr ( z ) = k e k r ( z )+ k , (8)where E is the Young’s modulus of the blood vessel, h is the wall thickness, P is the reference pressure, r is the reference radius, and k , k , and k are empirically derived constants.Furthermore, for models with multiple vascular branches, mass conservation is obeyed at the junctionregions, and the pressure is assumed to be constant between the inlet and outlets of the junctions [6].As with 3D models, we need initial conditions to initialize the simulation, Q ( z, t = ) = Q ( z ) , P ( z, t = ) = P ( z ) , S ( z, t = ) = S ( z ) , (9)as well as inlet and outlet boundary conditions, discussed in section 2.3, to represent the portion of thecardiovascular not reflected in our 1D model, Q ( z, t ) = Q in ( t ) , in Γ in , (10) P ( z, t ) = f ( z, t, Q out ( t ) , ˙ Q out ( t ) , φ ) , in Γ in , (11)We generate and simulate the 1D centerline and finite element models using SimVascular and VMTK[21]. We solve the set of differential equations using our open-source solver svOneDSolver ( github.com/SimVascular/svOneDSolver ). We employ a stabilized space-time finite element method based onthe discontinuous Galerkin method in time [6]. The spatial discretization employs continuous piecewiselinear polynomials whereas we use a piecewise constant temporal discretization. For more backgroundinformation, see [22, 23, 24]. The third model fidelity we consider in this work is the zero-dimensional (0D) model. Unlike 3D and1D models, 0D models lack spatial information. However, as with 1D models, 0D models are capable ofaccurately simulating bulk flow rate and bulk pressure quantities in the cardiovascular system [25], [26].These 0D models are composed of individual lumped-parameter elements that connect to form an entirecomplex lumped parameter network. There are many lumped parameter elements commonly used in thecontext of cardiovascular modeling and simulation. Some of these elements include resistors, capacitors,and inductors [26]. Resistors model the viscous effects of blood flow, capacitors represent the elasticnature of blood vessels, while inductors capture the inertia of the blood flow. In lumped parameternetworks that model the heart, diode elements are also employed. These diodes mimic the behavior ofheart valves, where they allow flow to pass only when the valves are open [27], [28]. The flow ratesand pressures in each of these elements are respectively governed by the following linear differential andalgebraic equations, ∆ P = RQ, Q = C ∆ ˙ P , ∆ P = L ˙ Q, Q = ∣ Q ∣ + Q , (12)where R is the resistance, C is the capacitance, L is the inductance, and ∆ P is the pressure drop acrossthe element. Lumped parameter 0D models are also analogous to electrical circuits, where the flow rateand pressure are representative of current and voltage, respectively. In this work, we will focus ourattention to just resistors and capacitors. 3umped parameter networks are commonly used in two contexts. First, 0D models can be used assurrogate models of entire vascular anatomies [25]. In this case, each blood vessel in the vascular systemis represented by one or more lumped-parameter elements in the 0D model, where the value of eachelement is determined by the geometric and material properties of the blood and blood vessel. Second,0D models can be used to represent boundary conditions in 3D and 1D models [27], [29], [17], [30], [31].In this context, each 0D element reflects a different downstream (or upstream) anatomical feature ofthe cardiovascular system. To employ the 0D models as boundary conditions, they must be numericallycoupled to the 3D or 1D finite element models [17]. The governing equations for these 0D modelsthen can be numerically solved using traditional time-stepping schemes, such as an explicit fourth-orderRunge-Kutta method.Some of the 0D models most commonly employed as boundary conditions are the resistance modeland the 3-element Windkessel model [5, 31]. The resistance model is composed of a single linear resistorelement that captures the downstream resistance of the vascular network not portrayed in the 3D or 1Dmodel. On the other hand, the 3-element Windkessel model, also known as the RCR model, as shown inFigure 1, models the proximal resistance, distal resistance, and compliance of the downstream vasculatureusing two resistor elements and a capacitor. The RCR boundary condition is discussed further in thenext section.To simulate our 0D surrogate models, which represent the vessels using linear resistors, we solve thegoverning system of equations using an in-house, modular numerical solver that employs the implicitgeneralized- α scheme for time advancement [32].Figure 1: Flowchart of method to check periodic state of 3D simulation results.4 .4 The RCR boundary condition Thus far, we have introduced three different model fidelities commonly used in computational vascularmodeling and simulation. Before any of the models and simulation results can be used in predictiveand scientific applications, the quantities of interest, primarily the flow rates and pressures, must besimulated until they converge to a periodic state. Here, a simulated quantity of interest is consideredto be periodic if its values between 2 adjacent periods are the same, within some defined tolerance. Aperiod in the cardiovascular context is typically defined as a single cardiac cycle. We expound on theconcept of a periodic state in this section, using the RCR model to illustrate the primary concepts.The three parameters of the RCR boundary condition are commonly tuned to clinical measurements,i.e., phase-contrast magnetic resonance imaging and pressure measurements by solving an optimizationproblem to match minimal and maximal pressure over a cardiac cycle, measured flows, and flow distribu-tion among different outlets [33]. A common strategy is to tune the total resistance and capacitance inthe model to produce a physiologic pressure waveform, then distributed proportional to the vessel outletareas [34].We begin by reviewing the response of a single RCR model to pulsatile inflow. The governingdifferential equation for the inlet pressure, P , of the 3-element Windkessel model, as a function of theinlet flow rate, Q , is ˙ P + Pτ = R p ˙ Q + τ ( R p + R d ) Q, τ = R d C, (13)where R p is the proximal resistance, C is the capacitance, and R d is the distal resistance. We obtain thesemi-analytical solution [29] for this ordinary differential equation as P ( t ) = [ P ( ) − R p Q ( )] e − t / τ ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ I + R p Q ( t )·„„„„„„„„„„„‚„„„„„„„„„„¶ II + ∫ t e −( t − ˜ t )/ τ C Q ( ˜ t ) d˜ t ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ III , τ = R d C > , (14)which depends on the inflow Q and the time constant τ . We can identify three different terms in thisequation: (I) exponential decay of the initial solution, (II) pressure drop at proximal resistance, (III)pressure drop at sub-circuit R d C . Assuming a constant inflow ¯ Q > t >
0, we obtain the pressurestep response for the RCR boundary condition as P ( t ) = P ∞ + e − t / τ [ P − P ∞ ] , P = P ( ) , lim t →∞ P ( t ) = P ∞ = ¯ Q ( R p + R d ) , (15)starting at the initial pressure P and exponentially approaching the asymptotic pressure P ∞ for thelimit t → ∞ . We now define the periodic inflow Q ( t ) as Q ( t + T ) = Q ( t ) , ¯ Q = T ∫ T Q ( t ) d t, (16)where the period T is the length of a cardiac cycle and ¯ Q the time-averaged mean flow. The meanpressure in the n -th cardiac cycle is denoted by¯ P n = T ∫ ( n + ) TnT P ( t ) d t, n ∈ N + . (17)With this notation, Equation (15) can be reformulated for a non-constant, periodic pulsatile inflow as¯ P n = ¯ P ∞ + e − nT / τ [ ¯ P − ¯ P ∞ ] , ¯ P ∞ = lim n →∞ ¯ P n , (18)starting at the initial mean pressure ¯ P and approaching the asymptotic mean pressure ¯ P ∞ . Note that theasymptotic mean pressure ¯ P ∞ can in general not be determined analytically but depends on the function Q ( t ) . Notably, ¯ P ∞ is different from the asymptotic pressure P ∞ that P approaches for a solution withsteady mean flow ¯ Q .In the remainder of this section, we will define an error metric to quantify the difference between thepressure ¯ P n in the n -th cardiac cycle and the asymptotic pressure ¯ P ∞ . For simplicity of notation, weshow the following derivations for ¯ P n < ¯ P ∞ , i.e. the pressure approaches the asymptotic pressure “from5elow”. However, the conclusions hold for any choice of initial pressure. We define the asymptotic error (cid:15) ∞ as 0 < (cid:15) ∞ ≤ ¯ P ∞ − ¯ P n ¯ P ∞ = e − nT / τ ⋅ ( − ¯ P ¯ P ∞ ) . (19)Thus, we can calculate the number of cardiac cycles n ∞ required for the RCR boundary condition toreach a periodic state within the tolerance (cid:15) ∞ as n ∞ ≤ − τT ⋅ ln (cid:15) ∞ − ¯ P / ¯ P ∞ , n ∞ ∈ N + . (20)Several observations can be made from Equation (20). First, the number of cycles scales linearly with τ / T , the ratio of the length of the time constant τ to the length of one cardiac cycle. A two-fold increasein the time constant doubles the number of cardiac cycles required to reach the same tolerance (cid:15) ∞ .Second, the closer the initial pressure ¯ P to the asymptotic pressure ¯ P ∞ the fewer cardiac cycles n ∞ arerequired to achieve periodic convergence with an error (cid:15) ∞ . The error (cid:15) ∞ and the number of cardiac cycles n ∞ can in general not be determined a priori since the constant ¯ P ∞ cannot be evaluated analytically.Only in the special case of starting from zero initial conditions, i.e. P =
0, the number of cardiac cyclesto periodic convergence can be directly given as n ∞ ≤ − τT ⋅ ln (cid:15) ∞ . (21)We can easily compute the difference ¯ P n − ¯ P n − between two cardiac cycles numerically without thesolution having reached a periodic state. Using this difference, we define the cyclic error (cid:15) n between twoconsecutive cardiac cycles as0 < (cid:15) n ≤ ¯ P n − ¯ P n − ¯ P ∞ = e − nT / ¯ τ ⋅ ( − ¯ P ¯ P ∞ ) ⋅ [ e T / τ − ] , n ≥ . (22)We then define the ratio α between asymptotic error (cid:15) ∞ and cyclic error (cid:15) n as α = (cid:15) ∞ (cid:15) n = ¯ P ∞ − ¯ P n ¯ P n − ¯ P n − = e T / τ − , (cid:32) (cid:15) ∞ > (cid:15) n , for τT > ≈ .
44 (23)Equation (23) shows that in general (cid:15) n ≠ (cid:15) ∞ . While error metrics like (cid:15) n are commonly used in practiceto determine whether a simulation as reached a periodic state within a given tolerance, it is not anadequate metric to do so. For time constants τ larger than 1.44 times the length of the cardiac cycle T ,the cyclic error (cid:15) n will underestimate the asymptotic error (cid:15) ∞ . In simulations with a large time constant τ ≫ T the error (cid:15) n will appear small despite the solution being far from a periodic state due to theirslow convergence. As a remedy, we propose a method in Section 2.5 to estimate the asymptotic error (cid:15) ∞ as the simulation is being computed. The previous section introduced analytical and theoretical methods to check the periodicity of a singleRCR model subjected to a prescribed inflow condition. This section extends this discussion to considercomputational vascular models and present a method to check the periodicity of patient-specific, multi-outlet 3D models, where each outlet is coupled to a different RCR boundary condition. The steps forthis method are summarized in the flowchart shown in Figure 1.As mentioned in the previous section, ¯ P ∞ generally cannot be predicted analytically. However, avalue of ¯ P ∞ is needed to check if the simulated flow rate and pressure at each outlet of the 3D modelhave converged to a periodic state. To remedy this issue, we estimate ¯ P ∞ for each outlet by simulatinga separate 0D model of the attached RCR boundary condition. The inflow to each 0D RCR model isthe simulated 3D flow rate, Q D , corresponding to that outlet. We then simulate these simple 0D RCRmodels for many cardiac cycles, to guarantee periodic convergence of the 0D RCR model, and use thefinal simulated pressure values, P D , as our estimates of ¯ P ∞ .Furthermore, as previously discussed, (cid:15) ∞ represents the asymptotic error of the pressure solution.We therefore specify a desired value for (cid:15) ∞ and use this as our criteria for periodic convergence. Our 3D6ressure solution, P D , is considered to be periodically converged, as per Equation (19), if it matches P D within this error threshold. Mathematically, this condition is expressed as, P D − P D ( Q D ) P D ( Q D ) ≤ (cid:15) ∞ . (24)Note that Q D should be reasonably close to a periodic state before we can use it as the inflow to ourRCR model. Typically, Q D converges much faster, i.e., within one or two cardiac cycles, than ¯ P D does.This note will be further discussed in Section 3.1.3.Lastly, each outlet of our 3D models has in general a unique time constant τ . Although, in practice,parameter tuning yields time constants that are similar. It can be shown that the convergence of eachoutlet is determined by a single model time constant ¯ τ that can be approximated by the mean of allindividual time constants. As such, we can use this average time constant to estimate the number ofcardiac cycles for which our 3D models must be simulated to achieve periodic convergence. As previously discussed, 3D models must be simulated until the flow rates and pressures converge to aperiodic state before they can be applied in scientific or clinical investigations, which currently requiresrunning the simulation for several cardiac cycles. This limits the utility of computational 3D modelsand simulations in real-world clinical applications. To alleviate this bottleneck, we introduce a novelmethod to jumpstart the initialization of 3D simulations. In particular, we generate initial conditions tominimize the number of cardiac cycles required to reach periodic convergence within a chosen tolerance.This process is fully automated in SimVascular [13] and requires no user input. Starting from a periodic1D solution (Section 2.6.1) we create a map from the centerline, the 1D representation of the model, tothe 3D volume mesh (Section 2.6.2), and generate an artificial initial velocity vector field (Section 2.6.3)for simulation initialization.
We automatically generate a 1D flow model of the high-fidelity 3D flow model using the SimVascular1D-plugin. We then run the model until it achieves periodic convergence and extract the solution of thelast cardiac cycle.
To map the 1D solution to the 3D Finite Element mesh, we create a map I from nodes P on the1D centerline to nodes P in the 3D volume mesh. This allows us to map quantities defined on thecenterline to the volume mesh, such as 1D flow and pressure, cross-sectional area, and normal vectors.The iterative process is outlined in Algorithm 1 and visualized in Figure 2. In Line 3, we first create aset of seed points P seed consisting of volume mesh nodes P that are closest to the centerline nodes P and store the corresponding indices in I . To do this, we find the shortest Euclidean distances betweencenterline coordinates x p and 3D mesh coordinates x j . We then employ a region growing algorithm(Line 14) to grow the 1D-3D map outwards, starting from the centerline seed points P seed . The algorithmin Line 17 selects nodes P new in the new layer from the previous layer P old using cell-connectivity. Finally,the map is expanded in Line 30 by assigning nodes in the new layer P new to the same 1D node as theclosest 3D node in the previous layer P old . This results in centerline nodes being roughly assigned to 3Dmesh nodes within the same cross-section. For 3D meshes with O( ) nodes Algorithm 1 only takes afew seconds to complete on a single CPU. We generate this map only once and use it repeatedly to mapvarious centerline quantities to the volume mesh, see Figure 3. We directly map the pressure from the centerline to the volume mesh using the map generated inSection 2.6.2, resulting in a pressure that’s approximately constant over the cross-section of the vessels.However, as the 1D solution only provides a scalar flow along the centerline, we must generate a velocityvector field from scratch. The ingredients for the velocity vector field are visualized in Figure 3. Fromthe mapped 1D flow rate and the cross-sectional area we calculate the velocity magnitude. Assuming aPoiseuille flow, we apply a parabolic flow profile to the velocity magnitude. With the help of a normalizedradial coordinate, we prescribe the flow profile to be maximal on the centerline and zero on the boundary.7 lgorithm 1:
Mapping centerline nodes to nodes in the 3D volume mesh. // initialize empty map I ← [ , . . . , ] // initialize 3D seed points P seed ← {} // loop all 1D points for p ∈ P do // find closest 3D point q = arg min j ∈ P ∣∣ x p − x j ∣∣ if q ∉ P seed then // add point to 3D seed points P seed ← P seed ∪ { q } // assign map I q = p end end C all ← {} P all ← P seed P new ← P seed i ← // in each iteration, grow seed points one layer outward while ∣ P new ∣ > do P old ← P new P new ← {} // loop all 3D points in previous layer for p ∈ P old do // loop all 3D cells connected to point for c ∈ PointCells ( p ) do // skip 3D cells in previous layers if c ∉ C all then C all ← C all ∪ { c } // loop all 3D points connected to cell for q ∈ CellPoints ( c ) do // skip 3D points in previous layers if q ∉ P all then // add 3D points to new layer P new ← P new ∪ { q } P all ← P all ∪ { q } end end end end end // loop all 3D points in new layer for p ∈ P new do // assign map according to map of closest 3D point in previous layer I p = I q with q = arg min j ∈ P old ∣∣ x p − x j ∣∣ end i ← i + end P all for iterations i ∈ { , , , } until Algorithm 1 converges.Finally, we multiply the scalar velocity magnitude with the centerline tangent vector to generate a vectorfield. We aim to preserve the amount of flow through a cross-section of the vessel as obtained from the1D solution while approximating a somewhat physiological velocity field. This velocity vector field canthen be used for simulation initialization.Figure 3: Ingredients to generate the velocity field from a 1D solution. From left to right: Flow mappedfrom 1D, cross-sectional area, radial coordinate, and normal vectors Throughout this section, we consider a simulation periodically converged if the asymptotic error is (cid:15) ∞ ≤ vascularmodel.org ) [35]. Six out of these models are shown in Table 1. In this section, we detail the convergence behavior of cardiovascular fluid dynamics simulations withthree-element Windkessel, or RCR, boundary conditions. The models used in this section have differentnumbers of outlets, ranging from 4 in the aortic models to 33 in the pulmonary artery models. Foreach model, we automatically generate a reduced-order 0D model as outlined in Section 2.3. Using acomputationally inexpensive 0D model allows us to run a large number of simulations for many cardiaccycles and analyze their convergence behavior in detail. These 0D results are directly applicable to 1Dand 3D simulation models and form the foundation for generating initial conditions in Section 3.2.2.9 D Type aorta-femoral pulmonary aorta aorta aorta aorta-femoral
State normal Glenn end-to-endanastomosis coarctation aneurysm aneurysmTable 1: Subset of models from the Cardiovascular Model Repository.
The time constant τ is the metric of a boundary condition that determines the rate of periodic conver-gence. It is thus essential for all numerical experiments in this work. Throughout the remainder of thiswork, we normalize τ by the length of the cardiac cycle T , with all values reported in Figure 4. Allnormalized time constants fall within the range [ , ] , with time constants in pulmonary models beingthe lowest. Furthermore, the time constants of different outlets in one model all fall within a narrowrange. In addition to the outlet time constants, we also show the model time constants ¯ τ / T as crosses.We extract the model time constants from the slope of the exponential curves of the asymptotic error (cid:15) ∞ (see Equation (20)) when running the models for several cardiac cycles. Each model has one uniquemodel time constant ¯ τ / T that is approximately the mean value of the outlet time constants.Figure 4: Time constants for models from the Vascular Model Repository. Bars indicate the range oftime constants τ / T for all outlets, crosses indicate the model time constant ¯ τ / T . The colors indicate themodel category: aorta (blue), aorta-femoral (purple), and pulmonary (yellow). We estimated the number of cardiac cycles n ∞ required to reach a periodic pressure in Equation (20). Inthe special case of zero initial conditions, this relationship simplified to the inequality in Equation (21),scaling linearly with the time constant τ / T and the logarithm of the asymptotic error (cid:15) ∞ . This relation-ship holds regardless of model fidelity, i.e. for 0D, 1D, and 3D. In Figure 5a, we record the number ofcardiac cycles required to reach an asymptotic error of (cid:15) ∞ ≤ τ / T . For pressure (left), the results confirm the linear relationshipbetween the single model time constant ¯ τ / T and the number of cardiac cycles n ∞ from Equation (21)10black line). For a range of model time constants ¯ τ / T ∈ [ . , . ] we find the range of cardiac cycles n ∞ ∈ [ , ] . The number of cardiac cycles to reach a periodic flow solution (right) is not correlated tothe model time constant and is, in general, much lower than the number of cycles n ∞ to reach a periodicpressure. Flow rate commonly converges within one cardiac cycle, with a maximum of eight cardiaccycles in our cohort of models.The same study is repeated in Figure 5b for starting from steady state initial conditions. Comparedto zero initial conditions, the number of cardiac cycles required to each an asymptotic pressure withan error of (cid:15) ∞ ≤ n ∞ ∈ [ , ] . Flow is converged within a maximum of 5 cardiaccycles. When starting from steady state initial conditions, the number of cardiac cycles n ∞ cannot begiven analytically, it depends on the prescribed inflow profile. However, the number of cardiac cyclescommonly still increases with the model time constant. (a) Starting from zero initial conditions(b) Starting from steady state initial conditions Figure 5: Number of cardiac cycles n ∞ required to reach an asymptotic error e ∞ = We visualize the convergence of pressure and flow for n =
30 cardiac cycles in Figure 6 for model 0107 0001(normal aorta). The pressure curve (top left) builds up slowly in each cardiac cycle, starting from zero,whereas the flow curve (bottom) is close to periodic starting from the first cycle. Taking the meanvalue over each cardiac cycle, both solutions exponentially approach their periodic state, as shown inEquation (18). The logarithmic plots in the two rightmost columns show the exponential decay of thecyclic error (cid:15) n and the asymptotic error (cid:15) ∞ , as defined in Equations (20) and (22), respectively. Notethat the flow (bottom) exhibits smaller errors than the pressure (top) and converges faster during thefirst few cardiac cycles. After that, both flow and pressure at all outlets converge with the model time11onstant ¯ τ / T ≈ .
4. This model time constant yields a factor between cyclic and asymptotic error of α = (cid:15) ∞ / (cid:15) n ≈ .
9. The threshold is indicated by horizontal lines in the error plots. Here, the solutionis converged after n ∞ =
21 cardiac cycles, reaching errors of (cid:15) ∞ = (cid:15) n = (cid:15) ∞ / α ≈ .
26 %. Thisexample demonstrates that for simulations with ¯ τ / T > / ln 2 ≈ .
44 the cyclic error (cid:15) n underestimatesthe asymptotic error (cid:15) ∞ .Figure 6: Convergence of 0D pressure (top) and flow (bottom) solution in geometry 0107 0001 (normalaorta) for n =
30 cardiac cycles. From left to right: Solution, mean solution over one cardiac cycle(normalized by asymptotic solution), cyclic error (cid:15) n , and asymptotic error (cid:15) ∞ . Each color refers to oneof the outlets. As explained in Section 2.4, the periodic error (cid:15) ∞ , comparing the current cardiac cycle to a perfectlyperiodic cycle, can in general not be computed analytically. It requires the periodic solution of themodel which is not known a priori . We thus outlined a method in Section 2.5 to use the 0D lumped-parameter boundary condition to estimate the periodic cycle numerically. Figure 7 shows the pressureat all outlets of model 0107 0001 from cycle one to cycle n ∞ =
21 (from blue to red), starting from zeroinitial conditions as in Figure 6. The top row shows the pressure at the outlets of the 3D model in eachcardiac cycle. The bottom row shows the prediction of the periodic state using the flow at the outlets ofthe 3D model in each cardiac cycle. Since flow converges much faster than pressure in this model, seeFigure 6, the cycle-to-cycle variation is minimal. The periodic state can be accurately predicted evenfrom early cardiac cycles, where the actual pressure of the model has not yet converged.
To demonstrate the performance of our 3D initialization method, we compare results for the 6 modelsshown in Table 1. The chosen models contain a wide range of anatomies and physiological conditions andare specified by an eight-digit ID. They include aorta and femoral arteries (0003 0001), pulmonary arteriesafter a Glenn procedure (0097 0001), aortic coarctation post end-to-end anastomosis (0107 0001), un-treated aortic coarctation (0111 0001), aortic aneurysm in a patient with Marfan Syndrome (0130 0000),and abdominal aortic aneurysm (0156 0001).
In this section, we demonstrate the performance of our pipeline to generate 3D initial conditions thatgreatly reduce the number of cardiac cycles, as proposed in Section 2.6. As a ground truth, i.e. the idealinitial condition, we use an asymptotic 3D solution extracted after reaching periodic convergence. Usingthis initial condition, the 3D simulation would reach periodic converge within one cardiac cycle. We12igure 7: Pressure for model 0107 0001 for n ∞ =
21 cardiac cycles (blue to red) starting from zeroinitial conditions until reaching periodic convergence. Pressure at the 3D outlets (top) and peridiocstate predicted from 0D (bottom).compare this solution to initial conditions we generated from a periodic 1D solution, using the mappingtechnique from Section 2.6.2. By using initial conditions mapped from 1D to 3D, we introduce two kindsof errors. First, the 1D solution is computed on a highly simplified geometrical representation of the 3Dgeometry assuming a constant flow profile, see Section 2.2. Second, the mapping process tries to butcannot guarantee local preservation of the 1D solution characteristics, such as incompressibility.Figure 8 shows these errors by comparing different initial conditions: 3D ground truth (blue), 1Dsolution (green), and 3D mapped from 1D (orange). We integrate the 3D ground truth and the mapped3D initial conditions over the cross-section of the 3D geometry continuously along the centerline. Thisallows us to plot pressure (top) and flow (bottom) continuously over the vessel path. For a perfect 1Dapproximation of the 3D simulation, the blue and green lines would overlap. For a perfect mappingfrom 1D to 3D, the blue and orange lines would overlap. In general, pressure is approximated well bythe 1D model and mapped well to the 3D domain. The differences between the models ( < ∼
75 mmHg). Flow is approximated well by the 1D modelwhereas the the mapping results in oscillations and overestimates the actual flow. Since all three modelsshould represent a model with rigid walls and incompressible flow, flow should be constant along thevessel path. However, the 3D mapped from 1D solution still roughly represents the correct flow splits tothe different vessel branches.Figure 8: Initial conditions in aortic aneurysm model 0130 0000 for different models: 3D ground truth(blue), 1D solution (green), 3D mapped from 1D (orange).13he differences between 3D ground truth and 3D mapped from 1D initial conditions are visualizedin Figure 9a and are evident, especially within the aneurysm. The mapping from 1D to 3D results in apressure that is approximately constant over the cross-sections of the vasculature. However, the overallvariation, 74-75 mmHg in the ground truth and 73-76 mmHg in the extrapolated solution, is negligiblysmall compared to the overall pressure level. Similarly, the velocity field in Figure 9b is visibly differentfor the ground truth (left) and mapped solution (right). (a) Pressure: ground truth (left),mapped from 1D (right) (b) Velocity (mag.): ground truth(left), mapped from 1D (right)
Figure 9: Initial conditions mapped from 1D compared to the 3D asymptotic solution in aortic aneurysmmodel 0130 0000.
Finally, we demonstrate that 3D initial conditions mapped from periodic 1D solutions can greatly reducethe 3D simulation time. We compare three kinds of initial conditions: Zero, Steady, and 1D. For zeroinitial conditions, we prescribe uniformly P ( x ) = and v = − ⋅ [ , , ] . To generate steady initialconditions, we first run a steady state simulation with constant mean inflow ¯ Q . We compared starting thepulsatile flow simulation at different time steps for model 0003 0001. Common choices are either the firsttime step t = t = (cid:15) ∞ for pressure (top) and flow (bottom) atall outlets (colors) for all three initial conditions: zero (left), steady (middle), 1D (right). The threshold (cid:15) ∞ ≤ n ∞ required to reach thatthreshold is indicated in each simulation by a vertical line. All simulations are shown for a total of n = − as other numerical errors in the simulations outweigh the asymptotic error. As previously observed,pressure converges faster than flow when using zero or steady initial conditions. However, pressure andflow are converged in one and two cardiac cycles, respectively, when using 1D initial conditions. Dueto the drastic reduction in the number of cardiac cycles required for the pressure solution, flow is nowslightly slower to converge.Convergence results for all selected 3D geometries and are reported in Table 2. While initializing witha steady state solution considerably reduces the number of cardiac cycles to reach convergence, most of themodels still require at least five cardiac cycles for convergence. In all models, the simulations initializedwith the 1D solution converge in one or two cardiac cycles. This leads to speedup factors between oneand nine compared to the steady state initialization, saving between 0 and 53 h of computation time on48 CPU cores each per 3D model. Only in model 0097 0001 which has a very low model time constant14igure 10: Convergence of the 3D pressure in geometry 0003 0001 (aorta-femoral) for different initialconditions: zero (left), steady (middle), 1D (right). It shows pressure (top) and flow (bottom) overmultiple cardiac cycles. Each color refers to one of the outlets. The horizontal black line indicatesthe asymptotic error of (cid:15) ∞ ≤ ID t τ / T Zero 12 / 6 / 2 / 1 / 1 / 1 / 3 Steady 7 / 3 / 2 / 1 / 2 / 3 / 2 / 1 / 1 / 1 1 / t saved
11 h 0 9.5 h 25 h 53 h 19 hTable 2: Simulation time t for one cardiac cycle of the 3D simulation, model time constant ¯ τ / T , andnumber of cardiac cycles (pressure / flow) for different initial conditions (zero, steady, 1D) with (cid:15) ∞ ≤ We gave a detailed review of properties of lumped-parameter boundary conditions in cardiovascular fluiddynamics simulations. The speed of convergence to a periodic state solution is determined by a singlemodel time constant, that can be approximated as the mean of all individual time constants of multipleoutlet boundary conditions. We found that for 53 models of the Vascular Model Repository, the modeltime constant ¯ τ / T spanned from 0.3 in pulmonary models to 9.6 in aorta models. That means that inthe most extreme case, the time constant of the model is almost ten times as large as the length of thecardiac cycle. Using zero initial conditions, the number of cardiac cycles required to reach a periodicstate scales linearly with the model time constant, reaching a median of eight cardiac cycles for oursubset of models. Even when using initial conditions other than zero, e.g., from a steady state solution,15he number of cardiac cycles required still scales with model time constant.Figure 11: Pressure convergence of an outlet in model 0076 1001 with time constant ¯ τ / T = .
6, startingfrom zero initial conditions. The red cardiac cycle was selected at (cid:15) n = (cid:15) ∞ = a priori , we proposed to use a 0D model of the boundarycondition. Taking the flow at a 3D outlet as input, we estimate the periodic pressure and compare it tothe pressure at the 3D outlet. This method is a quick and easy validation tool for any 3D vascular fluiddynamics simulation with lumped-parameter boundary conditions. In models with time constants largerthan 1 / ln 2 ≈ .
44, the cycle-to-cycle difference is lower than the asymptotic error. A low cycle-to-cycledifference can erroneously lead the user to believe that the simulation has already reached a periodicstate when comparing two consecutive cardiac cycles. This is visualized in Figure 11. Using the sameerror tolerance of 1 % but different error norms, cyclic error (red) and asymptotic error (blue), yields apressure difference of 10 mmHg. Thus, only the asymptotic error is a suitable metric to determine thedistance of a simulation to its periodic state.To reduce 3D computation times, we proposed a method to generate good initial conditions in anautomated and computationally inexpensive way. We used SimVascular to automatically generate a1D replication of the 3D model, which we ran until we achieved a periodic state. We than mappedthe 1D solution onto the 3D Finite Element mesh. Note that the mapping does in general not resultin a fluid field that is divergence free. However, we are not interested in generating a physical orphysiological meaningful solution. Instead, the initial conditions are “washed out” by the first iterationof our numerical solver of the Navier Stokes equations. Similarly, the pressure field does not capture localvariations. Here, it should be noted that it is much more important to match the overall pressure levelof the model, which is in general much higher than any local variations. In a convergence analysis withsix different vascular models, we demonstrated that models with our 1D initialization method convergein one or two cardiac cycles. This greatly reduces the computation time for the 3D model over thestandard method of initializing pulsatile simulations with a steady-state solution, typically dozens ofhours in simulation time and hundreds of hours in CPU time. The code for the 1D initialization is freelyavailable on GitHub ( github.com/SimVascular/SimVascular ).We close with a discussion of the the limitations and future perspectives of our work. We did onlyconsider RCR boundary conditions in this work. However, there are many more examples of lumpedparameter 0D networks that are coupled to 3D models, such as coronary boundary conditions [36, 37]and 0D closed loop models used, e.g., in simulations for single ventricle [38, 39] and coronary arterybypass graft (CABG) patients [40]. In the general case of lumped parameter boundary conditions, themodel time constant cannot be determined analytically. However, it can still be calculated numericallyfrom the time it takes the boundary condition to respond to a step in the inflow condition. We plan toinclude these boundary conditions in future work. It should be noted that those boundary conditionsnot only require initialization of the velocity and pressure field in the 3D model but also in unknownsthat are internal to the 0D model. Furthermore, we only considered rigid-wall simulation in this work.We plan to test our 1D initialization method for deformable wall simulations as well, using the coupledmomentum method [41] or an arbitrary Lagrangian–Eulerian formulation [42]. To further improve ourmethod and guarantee periodic convergence within one cardiac cycle, it will be necessary to improve the16apping of the velocity field to the 3D model. This could be accomplished for example by solving aStokes flow problem, using the 1D results as boundary conditions.
We thank Dr. Luca Pegolotti and Dr. Irene Vignon-Clementel for insightful discussions. This workwas supported by NIH grants R01LM013120 and R01EB029362. The authors gratefully acknowledgethe Stanford Research Computing Center for providing the computational resources necessary to thenumerical simulations presented in this work.
References [1] Ethan Kung, Alessia Baretta, Catriona Baker, Gregory Arbia, Giovanni Biglino, Chiara Corsini, SilviaSchievano, Irene E. Vignon-Clementel, Gabriele Dubini, Giancarlo Pennati, Andrew Taylor, Adam Dorfman,Anthony M. Hlavacek, Alison L. Marsden, Tain-Yen Hsia, and Francesco Migliavacca. Predictive modelingof the virtual hemi-fontan operation for second stage single ventricle palliation: Two patient-specific cases.
Journal of Biomechanics , 46(2):423–429, jan 2013.[2] Jongmin Seo, Casey Fleeter, Andrew M. Kahn, Alison L. Marsden, and Daniele E. Schiavazzi. Multi-fidelityestimators for coronary circulation models under clinically-informed data uncertainty. arXiv , 2019.[3] Casey M. Fleeter, Gianluca Geraci, Daniele E. Schiavazzi, Andrew M. Kahn, and Alison L. Marsden. Mul-tilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics.
Computer Methods inApplied Mechanics and Engineering , 365:113030, jun 2020.[4] Christopher J. Arthurs, Nan Xiao, Philippe Moireau, Tobias Schaeffter, and C. Alberto Figueroa. A flex-ible framework for sequential estimation of model parameters in computational hemodynamics.
AdvancedModeling and Simulation in Engineering Sciences , 7(1), dec 2020.[5] I.E. Vignon-Clementel, C.A. Figueroa, K.E. Jansen, and C.A. Taylor. Outflow boundary conditions for3d simulations of non-periodic blood flow and pressure fields in deformable arteries.
Computer Methods inBiomechanics and Biomedical Engineering , 13(5):625–640, oct 2010.[6] Jing Wan, Brooke Steele, Sean A. Spicer, Sven Strohband, Gonzalo R. Feijo ´ o, Thomas J.R. Hughes, andCharles A. Taylor. A one-dimensional finite element method for simulation-based medical planning forcardiovascular disease. Computer Methods in Biomechanics and Biomedical Engineering , 5(3):195–206, jan2002.[7] S. M. Moore, K. T. Moorhead, J. G. Chase, T. David, and J. Fink. One-dimensional and three-dimensionalmodels of cerebrovascular flow.
Journal of Biomechanical Engineering , 127(3):440–449, sep 2004.[8] L. Grinberg, E. Cheever, T. Anor, J. R. Madsen, and G. E. Karniadakis. Modeling blood flow circulation inintracranial arterial networks: A comparative 3d/1d simulation study.
Annals of Biomedical Engineering ,39(1):297–309, jul 2010.[9] Philippe Reymond, Fabienne Perren, Fran¸cois Lazeyras, and Nikos Stergiopulos. Patient-specific mean pres-sure drop in the systemic arterial tree, a comparison between 1-d and 3-d models.
Journal of Biomechanics ,45(15):2499–2505, oct 2012.[10] Nan Xiao, Jordi Alastruey, and C. Alberto Figueroa. A systematic comparison between 1-d and 3-d hemo-dynamics in compliant arterial models.
International Journal for Numerical Methods in Biomedical Engi-neering , 30(2):204–231, sep 2013.[11] S. Pant, B. Fabr`eges, J-F. Gerbeau, and I. E. Vignon-Clementel. A methodological paradigm for patient-specific multi-scale CFD simulations: from clinical measurements to parameter estimates for individualanalysis.
International Journal for Numerical Methods in Biomedical Engineering , 30(12):1614–1648, nov2014.[12] Giulia Bertaglia, Adri´an Navas-Montilla, Alessandro Valiani, Manuel Ignacio Monge Garc´ıa, Javier Murillo,and Valerio Caleffi. Computational hemodynamics in arteries with the one-dimensional augmented fluid-structure interaction system: viscoelastic parameters estimation and comparison with in-vivo data.
Journalof Biomechanics , 100:109595, feb 2020.[13] Adam Updegrove, Nathan M. Wilson, Jameson Merkow, Hongzhi Lan, Alison L. Marsden, and Shawn C.Shadden. SimVascular: An open source pipeline for cardiovascular simulation.
Annals of Biomedical Engi-neering , 45(3):525–541, dec 2016.[14] Leopoldo P. Franca and S´ergio L. Frey. Stabilized finite element methods: II. the incompressible navier-stokesequations.
Computer Methods in Applied Mechanics and Engineering , 99(2-3):209–233, sep 1992.
15] Charles A. Taylor, Thomas J.R. Hughes, and Christopher K. Zarins. Finite element modeling of blood flowin arteries.
Computer Methods in Applied Mechanics and Engineering , 158(1-2):155–196, may 1998.[16] Christian H. Whiting and Kenneth E. Jansen. A stabilized finite element method for the incompressiblenavier-stokes equations using a hierarchical basis.
International Journal for Numerical Methods in Fluids ,35(1):93–116, 2001.[17] Mahdi Esmaily-Moghadam, Yuri Bazilevs, and Alison L. Marsden. A new preconditioning technique forimplicitly coupled multidomain simulations with applications to hemodynamics.
Computational Mechanics ,52(5):1141–1152, may 2013.[18] Thomas J.R. Hughes and J. Lubliner. On the one-dimensional theory of blood flow in the larger vessels.
Mathematical Biosciences , 18(1-2):161–170, oct 1973.[19] B.N. Steele, Jing Wan, J.P. Ku, T.J.R. Hughes, and C.A. Taylor. In vivo validation of a one-dimensionalfinite-element method for predicting blood flow in cardiovascular bypass grafts.
IEEE Transactions onBiomedical Engineering , 50(6):649–656, jun 2003.[20] Mette S. Olufsen. Structured tree outflow condition for blood flow in larger systemic arteries.
AmericanJournal of Physiology-Heart and Circulatory Physiology , 276(1):H257–H268, jan 1999.[21] Luca Antiga, Marina Piccinelli, Lorenzo Botti, Bogdan Ene-Iordache, Andrea Remuzzi, and David A. Stein-man. An image-based modeling framework for patient-specific computational hemodynamics.
Medical &Biological Engineering & Computing , 46(11):1097–1112, nov 2008.[22] Alexander N. Brooks and Thomas J.R. Hughes. Streamline upwind/petrov-galerkin formulations for con-vection dominated flows with particular emphasis on the incompressible navier-stokes equations.
ComputerMethods in Applied Mechanics and Engineering , 32(1-3):199–259, sep 1982.[23] Thomas J.R. Hughes and Michel Mallet. A new finite element formulation for computational fluid dynamics:III. the generalized streamline operator for multidimensional advective-diffusive systems.
Computer Methodsin Applied Mechanics and Engineering , 58(3):305–328, nov 1986.[24] Thomas J.R. Hughes, Leopoldo P. Franca, and Gregory M. Hulbert. A new finite element formulation forcomputational fluid dynamics: VIII. the galerkin/least-squares method for advective-diffusive equations.
Computer Methods in Applied Mechanics and Engineering , 73(2):173–189, may 1989.[25] Mehran Mirramezani and Shawn C. Shadden. A distributed lumped parameter model of blood flow.
Annalsof Biomedical Engineering , jul 2020.[26] Luca Formaggia, Alfio Quarteroni, and Alessandro Veneziani, editors.
Cardiovascular Mathematics . SpringerMilan, 2009.[27] H. J. Kim, I. E. Vignon-Clementel, C. A. Figueroa, J. F. LaDisa, K. E. Jansen, J. A. Feinstein, and C. A.Taylor. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model.
Annals of Biomedical Engineering , 37(11):2153–2169, jul 2009.[28] Daniele E. Schiavazzi, Alessia Baretta, Giancarlo Pennati, Tain-Yen Hsia, and Alison L. Marsden. Patient-specific parameter estimation in single-ventricle lumped circulation models under uncertainty.
InternationalJournal for Numerical Methods in Biomedical Engineering , 33(3):e02799, jun 2016.[29] Irene E. Vignon and Charles A. Taylor. Outflow boundary conditions for one-dimensional finite elementmodeling of blood flow and pressure waves in arteries.
Wave Motion , 39(4):361–374, apr 2004.[30] Mehran Mirramezani, Scott L. Diamond, Harold I. Litt, and Shawn C. Shadden. Reduced order modelsfor transstenotic pressure drop in the coronary arteries.
Journal of Biomechanical Engineering , 141(3), jan2019.[31] Irene E. Vignon-Clementel, C. Alberto Figueroa, Kenneth E. Jansen, and Charles A. Taylor. Outflowboundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries.
Computer Methods in Applied Mechanics and Engineering , 195(29-32):3776–3796, jun 2006.[32] Kenneth E. Jansen, Christian H. Whiting, and Gregory M. Hulbert. A generalized-alpha method for inte-grating the filtered navier–stokes equations with a stabilized finite element method.
Computer Methods inApplied Mechanics and Engineering , 190(3-4):305–319, oct 2000.[33] Ryan L. Spilker and Charles A. Taylor. Tuning multidomain hemodynamic simulations to match physiolog-ical measurements.
Annals of Biomedical Engineering , 38(8):2635–2648, mar 2010.[34] Yifang Zhou, Ghassan S Kassab, and Sabee Molloi. On the design of the coronary arterial tree: a general-ization of murray ' s law. Physics in Medicine and Biology , 44(12):2929–2945, nov 1999.[35] Nathan M. Wilson, Ana K. Ortiz, and Allison B. Johnson. The vascular model repository: A public resourceof medical imaging data and blood flow simulation results.
Journal of Medical Devices , 7(4), dec 2013.[36] H.J. Kim, I.E. Vignon-Clementel, C.A. Figueroa, K.E. Jansen, and C.A. Taylor. Developing computationalmethods for three-dimensional finite element simulations of coronary blood flow.
Finite Elements in Analysisand Design , 46(6):514–525, jun 2010.
37] H. J. Kim, I. E. Vignon-Clementel, J. S. Coogan, C. A. Figueroa, K. E. Jansen, and C. A. Taylor. Patient-specific modeling of blood flow and pressure in human coronary arteries.
Annals of Biomedical Engineering ,38(10):3195–3209, jun 2010.[38] E. L. Bove, M. R. de Leval, F. Migliavacca, R. Balossino, and G. Dubini. Toward optimal hemodynamics:Computer modeling of the fontan circuit.
Pediatric Cardiology , 28(6):477–481, aug 2007.[39] Edward L. Bove, Francesco Migliavacca, Marc R. de Leval, Rossella Balossino, Giancarlo Pennati, Thomas R.Lloyd, Sachin Khambadkone, Tain-Yen Hsia, and Gabriele Dubini. Use of mathematic modeling to compareand predict hemodynamic effects of the modified blalock–taussig and right ventricle–pulmonary artery shuntsfor hypoplastic left heart syndrome.
The Journal of Thoracic and Cardiovascular Surgery , 136(2):312–320.e2,aug 2008.[40] Abhay B. Ramachandra, Andrew M. Kahn, and Alison L. Marsden. Patient-specific simulations revealsignificant differences in mechanical stimuli in venous and arterial coronary grafts.
Journal of CardiovascularTranslational Research , 9(4):279–290, jul 2016.[41] C. Alberto Figueroa, Irene E. Vignon-Clementel, Kenneth E. Jansen, Thomas J.R. Hughes, and Charles A.Taylor. A coupled momentum method for modeling blood flow in three-dimensional deformable arteries.
Computer Methods in Applied Mechanics and Engineering , 195(41-43):5685–5706, aug 2006.[42] Kathrin B¨aumler, Vijay Vedula, Anna M. Sailer, Jongmin Seo, Peter Chiu, Gabriel Mistelbauer, Frandics P.Chan, Michael P. Fischbein, Alison L. Marsden, and Dominik Fleischmann. Fluid–structure interactionsimulations of patient-specific aortic dissection.
Biomechanics and Modeling in Mechanobiology , 19(5):1607–1628, jan 2020., 19(5):1607–1628, jan 2020.