Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Theor Comput Fluid Dyn. Author manuscript; available in PMC 2017 April 1.
Published in final edited form as:
PMCID: PMC4778980

Immersed boundary-finite element model of fluid–structure interaction in the aortic root


It has long been recognized that aortic root elasticity helps to ensure efficient aortic valve closure, but our understanding of the functional importance of the elasticity and geometry of the aortic root continues to e-volve as increasingly detailed in vivo imaging data become available. Herein, we describe a fluid-structure interaction model of the aortic root, including the aortic valve leaflets, the sinsuses of Valsalva, the aortic annulus, and the sinotubular junction, that employs a version of Peskin’s immersed boundary (IB) method with a finite element (FE) description of the structural elasticity. As in earlier work, we use a fiber-based model of the valve leaflets, but this study extends earlier IB models of the aortic root by employing an incompressible hyperelastic model of the mechanics of the sinuses and ascending aorta using a constitutive law fit to experimental data from human aortic root tissue. In vivo pressure loading is accounted for by a backward displacement method that determines the unloaded configurations of the root model. Our model yields realistic cardiac output at physiological pressures, with low transvalvular pressure differences during forward flow, minimal regurgitation during valve closure, and realistic pressure loads when the valve is closed during diastole. Further, results from high-resolution computations indicate that although the detailed leaflet and root kinematics show some grid sensitivity, our IB model of the aortic root nonetheless produces essentially grid-converged flow rates and pressures at practical grid spacings for the high-Reynolds number flows of the aortic root. These results thereby clarify minimum grid resolutions required by such models when used as stand-alone models of the aortic valve as well as when used to provide models of the outflow valves in models of left ventricular fluid dynamics.

Keywords: aortic valve, fluid-structure interaction, immersed boundary method, incompressible flow, hyperelasticity, finite element method, finite difference method

1 Introduction

The aortic root consists of the three semilunar aortic valve cusps, the three sinuses of Valsalva, which are bulbous cavities positioned behind each valve leaflet, the aortic annulus, where the aortic valve leaflets attach to the aorta, and the sinotubular junction, where the sinuses merge into the ascending aorta. The healthy aortic valve acts to ensure unidirectional flow of oxygenated blood from the heart to the tissues of the body, including to the heart itself via the coronary circulation. Diseases of the aortic valve can result in stenosis or regurgitation, and severe aortic valve disease is treated by repairing or by replacing the diseased valve with either a mechanical or a bioprosthetic valve [11]. Approximately 50,000 such procedures are performed each year in the United States [6, 26, 67]. Continuing advances in noninvasive in vivo imaging of blood flow and tissue deformation, mechanical tests specifically targeted to biological tissues, and computer modeling and simulation enable integrative studies of the dynamics and mechanics of the aortic root [41]. Such work has the potential to improve both medical devices and surgical procedures used to treat patients with valvular heart disease, and also clinical approaches to patient risk assessment and treatment planning.

Aortic compliance, and specifically the compliance of the aortic sinuses, has a fundamental role in the function of the aortic root. The sinuses act as reservoirs, storing blood during systole and then releasing it in diastole to facilitate flow in the coronary arteries [71]. The compliance of the sinuses also helps to ensure the proper closure of the aortic valve [4]. Further, the shape of the aortic root causes blood to circulate in the sinuses, generating vortices that act both on the valve leaflets and on the sinus walls. These vortices generate forces that facilitate effective valve closure and act as a regulatory mechanism that synchronizes closure [4]. Although this mechanism of efficient valve closure was first postulated by Leonardo da Vinci, there remains some debate in the surgical community about whether it is important to recreate the anatomic geometry of the native aortic root in aortic valve and root replacement procedures [12].

More recently, our understanding of the functional importance of the elasticity and anatomical geometry of the aortic root has evolved as increasingly detailed data have been acquired on the in vivo deformations that occur during the cardiac cycle. For instance, it has been shown that the aortic root undergoes a multi-modal series of conformational changes even before the leaflets open [18]. It has been argued that these deformations act to reduce shear stresses on the valve leaflets and thereby prolong the life of the native valve leaflets [12]. Annular and commissural flexibility may be a key component in this interaction, and annular flexibility is lost with the implantation of a stented artificial valve, possibly reducing the lifetime of bioprosthetic valve leaflets. The annular commissures also participate in load sharing, which reduces peak stresses on the aortic cusps immediately after valve closure. Devices such as the Medtronic 3f Aortic Bioprosthesis attempt to account for both the annular and the commissural flexibility of the native root, but despite the theoretical advantages offered by such designs, more traditional devices with rigid annuluses still dominate clinical practice.

A challenge in modeling and simulating the mechanical response of arteries, including the aorta, is that the arterial walls are continuously subjected to substantial intraluminal pressures in their in vivo state. Thus, arterial geometries that are derived from in vivo imaging generally correspond to a loaded configuration. In addition, as first demonstrated by Fung [27] and by Vaishnav and Vossoughi [72], even the unloaded configuration of the artery is not stress free. Instead, because of tissue growth and remodeling, the arterial wall includes residual stresses and is subject to axial tethering [10, 28, 41]. These factors all conspire to make determining the unloaded arterial geometry challenging. Indeed, both the deformation due to intraluminal pressure loading and the residual stresses cannot be measured directly in vivo, and therefore these quantities must be estimated. Residual stresses may be determined from ex vivo experiments, whereas the zero-pressure configuration of blood vessels can be determined by solving the inverse elastostatic problem.

Different approaches have been developed to solve the inverse elastostatic problem for complex arterial geometries derived from medical imaging data. Most of these strategies can be grouped in two broad categories. One group of approaches focuses on retrieving the initial deformation field and the initial stress field of the artery while keeping the imaging-derived geometry unaltered. An example of this type of approach is the modified updated Lagrangian formulation (MULF) introduced by Gee et al. [31], which is an incremental prestressing method that recovers the equilibrium configuration rather than the zero-pressure geometry. A second group of approaches determines an unloaded configuration of the vessel that, when subjected to intraluminal pressure loading, will inflate to match the imaging-derived geometry. Among the first to develop a stationary method to solve the inverse elastostatic problem were Govindjee and Mihalic [32]. This method was extended by Lu et al. [46] to arteries, but this approach does not appear to have been widely used in practice, possibly because it requires a modification in the finite element (FE) solution scheme that renders this method somewhat difficult to implement within commercial FE analysis software. In comparison, the backward incremental method developed by de Putter et al. [20] requires updating the deformation gradient and thereby can be implemented with fewer modifications to the finite element code. Perhaps the most straightforward approach to retrieving the zero-pressure configuration is the backward displacement method proposed by Bols et al. [5], which iteratively modifies the coordinates of the reference configuration rather than the deformation gradient tensor. We use the backward displacement method in this work to recover the unloaded configuration of an idealized model of the aortic root.

While the methods to obtain an estimation of the zero-pressure configuration of a blood vessel can be applied to any geometry, determining residual stresses is a local, vessel-specific task that presently requires harvesting and then destroying the blood vessel. Computational approaches to accounting for residual stresses are usually based on an assumed ‘open’ configuration, which is subsequently ‘closed,’ thereby reversing the process that releases residual stresses [2, 15]. This approach is difficult to apply to noncylindrical geometries and, in particular, cannot be readily applied to imaging-derived geometries [15]. A more generic approach to including residual stress in arteries is described in Alastrué et al. [2] and is based on a Kröner-Lee decomposition of the deformation gradient tensor. For simplicity, we do not consider residual stresses in the present work, although we do account for initial strains resulting from in vivo pressure loading.

Most the methods surveyed above have been developed and applied to study aneurysmal arteries [31, 46], although applications to normal vessels as well as other pathophysiological conditions are possible. Because arterial wall stress distributions are considered good predictors of aneurysm rupture, improvements in the accuracy of computed wall stress distributions are expected to facilitate improvements in aneurysm rupture prediction. In contrast, although the important role of aortic root compliance in valve closure is well documented [16, 17], many dynamic analyses of the aortic valve tend to focus on the mechanical behavior of the valve leaflets and often model the sinuses and the ascending aorta as rigid [34, 38, 74, 78] or as linear elastic materials [14, 49], but experimental data [3] show that the sinuses and the ascending aorta exhibit a nonlinear hyperelastic behavior. Examples of nonlinear hyperelastic constitutive models used in fluid-structure interaction simulations of the aortic root can be found in studies by De Hart et al. [19] and by Weinberg and Mofrad [75], which both adopt an arbitrary Lagrangian-Eulerian (ALE) approach. The application of ALE fluid-structure interaction schemes to modeling aortic valve dynamics is somewhat limited by the remeshing required by such methods as the structure deforms [68]. Attempts to tackle this limitation have included restricting the problem to a two-dimensional analysis [23] and the implementation of complex remeshing algorithms in three-dimensional analyses [13].

An alternative approach to modeling interactions between the blood and the aortic valve and root is offered by Peskin’s immersed boundary (IB) method [58], which was introduced to simulate cardiac valve dynamics [57] and was subsequently extended to model fluid-structure interaction in the heart and its valves [3436,38,44,47,48,5153,60]. This IB approach to fluid-structure interaction is to describe structural stresses and deformations in Lagrangian form, and to describe the momentum, viscosity, and incompressibility of the coupled fluid-structure system in Eulerian form. Integral transforms with Dirac delta function kernels mediate interaction between Lagrangian and Eulerian variables. When discretized, the IB formulation of the equations of motion replaces these singular kernel functions with regularized kernels that are designed to ensure conservation of physical quantities such as force and torque when converting between Lagrangian and Eulerian forms. This numerical approach allows the Lagrangian structural mesh to overlay the background Eulerian grid in an arbitrary manner, thereby avoiding the need to deploy dynamic body-fitted grids. In addition, because the immersed structures move according to a common interpolated velocity field, the IB method offers an implicit contact model that prevents the valve leaflets from interpenetrating, even when subjected to substantial diastolic pressure loads [34, 38]. Related sharp-interface IB methods have also been developed and used by Borazjani, Ge, and Sotiropoulos to simulate valvular dynamics for models of both rigid and flexible aortic valve prostheses [7, 8]. Previous work demonstrated that three-dimensional IB models of the aortic valve can yield physiological cardiac output at realistic pressures [34, 38]. However, we are aware of no previous study that demonstrates that Peskin’s IB method yields reasonably well-resolved simulation results for flexible aortic valve models in three spatial dimensions at practical grid spacings.

The primary aim of this study is to develop a new fluid-structure interaction model of the aortic root that substantially extends earlier IB models of aortic valve dynamics [34,38] by including descriptions of the elasticity, distensibility, and dynamics of the aortic sinuses and the ascending aorta. Herein, the aortic sinuses and proximal ascending aorta are modeled as an incompressible, isotropic, hyperelastic material with an exponential neo-Hookean strain-energy functional [21, 41] that we fit to experimental data obtained by Azadani et al. [3] from human aortic root tissue samples. A separate finite element analysis is employed to estimate the unloaded aortic geometry using a method based on the backward displacement method of Bols et al. [5]. Our fluid-structure interaction simulations employ hybrid models of the aortic root that use a fiber-based description of the thin aortic valve leaflets, as done in previous studies [34, 38], along with a finite element-based description [39] of the comparatively thick walls of the aortic sinuses and proximal ascending aorta.

In our dynamic analysis, we pace the aortic root model to an essentially periodic steady state using a prescribed, periodic left-ventricular pressure waveform. A human Windkessel model fit to clinical data [69] provides downstream loading for the aortic root, and a pressure waveform, taken from the clinical data set used to parameterize this reduced circulation model [55], serves to drive flow through the model. In the overall model, flow rates are not prescribed, but rather emerge from the computation. Physiological cardiac output is obtained at physiological driving and loading pressures, with low transvalvular pressure differences during forward flow, minimal regurgitation during valve closure, and realistic transvalvular pressure loads when the valve is closed during diastole.

The aortic root model considered in this work captures the complex interactions between the flow, the thin flexible aortic valve cusps, and the deformable walls of the aortic sinuses and the ascending aorta. It employs an idealized anatomical geometry, in which the initial and reference configurations of the aortic sinuses and valve cusps are assumed to exhibit three-fold symmetry, but can be extended to more realistic anatomic geometries. This model is fully three dimensional in the sense that there are no symmetry conditions imposed on the subsequent fluid dynamics or structural kinematics. Using this somewhat idealized model, we perform a refinement study to demonstrate that the present methods are able to yield essentially grid-converged flow rates and pressures at practical grid spacings. Specifically, under grid refinement, only relatively small differences are observed in stroke volume, maximum flow velocity, and vessel distensibility and stress distribution. Relatively large differences remain in the fine details of the flow in the vicinity of the valve cusps, and in the details of the kinematics of the valve cusps and the aortic wall. Higher spatial resolution is likely needed to resolve fully the dynamics of the valve during systole. Nonetheless, a key contribution of this study is that it demonstrates that bulk flow properties are reasonably resolved by IB models of aortic valve dynamics at practical spatial resolutions. These results are practically useful for determining minimum grid resolutions when using IB models of the outflow valves.

2 Methods

2.1 Immersed Boundary Formulation

The immersed boundary formulation used herein describes fluid-structure systems in which an elastic structure is immersed in a viscous incompressible fluid. It employs a Lagrangian description of the structural deformations and the stresses generated by those deformations, and an Eulerian description of the momentum, incompressibility, and viscosity of the coupled fluid-structure system. In the present model, we employ a fiber-based description of the thin valve leaflets, and we describe the aortic wall as an incompressible hyperelastic solid. Let q [set membership] U [subset or is implied by] R2 indicate Lagrangian curvilinear coordinates attached to the valve leaflets, with q = (q, r), let X [set membership] V [subset or is implied by] R3 indicate the Lagrangian material coordinate system of the aortic wall, with X = (X, Y, Z), and let x [set membership] Ω indicate the Eulerian physical coordinates of the physical domain, with x = (x, y, z). The physical position of fiber point q at time t is given by ϕ(q, t) [set membership] Ω, and the physical position of aortic wall material point X at time t is given by χ(X, t) [set membership] Ω. We use Ωl(t) [subset or is implied by] Ω to indicate the (codimension 1) physical region occupied by the valve leaflets at time t, Ωw(t) [subset or is implied by] Ω to indicate the (codimension 0) physical region occupied by the solid-body model of the vessel wall, and Ωf(t) = Ω \ Ωw(t) to indicate the physical region occupied by the fluid at time t.

The equations of motion for the coupled fluid-structure system are:







in which ρ is the mass density, μ is the dynamic viscosity, u(x, t) is the Eulerian velocity field of the fluid-structure system, p(x, t) is the Eulerian pressure field of the fluid-structure system, f (x, t) is the Eulerian elastic force density generated by deformations to the fiber model of the valve leaflets, g(x, t) is the Eulerian elastic force density generated by deformations to the solid-body model of the vessel wall, F(q, t) is the Lagrangian elastic force density of the fiber model, P(X, t) is the first Piola-Kirchhoff elastic stress tensor of the solid model, N(X) is the outward unit normal along [partial differential]V, δ(x) = δ(x) δ(y) δ(z) is the three-dimensional Dirac delta function, and D(·)Dt=(·)t+u·(·) is the material derivative. In the present work, we assume uniform mass density ρ and dynamic viscosity μ. These assumptions are not essential, however, and versions of the IB method have been developed that permit the use of spatially varying mass densities [25, 42, 43, 54, 62, 79] and viscosities [25, 62].

Eqs. (1) and (2) are the Eulerian incompressible Navier-Stokes equations. Here, the momentum equation (1) is augmented by two Eulerian body force densities. The first of these, f (x, t), is determined by eq. (3) to be the Eulerian force density that is equivalent to the Lagrangian fiber force density F(q, t). The second body force in (1), g(x, t), is determined by eq. (4) to be the Eulerian force density equivalent to the Lagrangian description of the forces generated by the solid-body model of the vessel wall, which are expressed in terms of P(X, t). Notice that f (x, t) is a singular force density supported on Ωl(t). By contrast, g(x, t) is supported on Ωw(t). Away from [partial differential]Ωw(t), g(x, t) is nonsingular, although g(x, t) is singular along [partial differential]Ωw(t) where PN0; see eq. (4).

Eqs. (5) and (6) determine the motions of the immersed structures from the Eulerian material velocity field u(x, t). Because of the presence of viscosity throughout the fluid-solid system, u(x, t) is continuous, and thus (5) and (6) are equivalent to



Because [nabla] · u(x, t) = 0, the immersed solid is automatically treated as incompressible in this formulation. Specifically, at least in the continuum formulation, there is no need to include terms in the elastic strain-energy functional to penalize compressible deformations. Notice that this formulation describes the inertia of both the fluid and the solid, and the immersed solid bodies are not assumed or required to be in (quasi-static) equilibrium.

In practice, we use a standard C0 finite element method to approximate the deformations and stresses of the vessel wall. To do so, it is useful to adopt a weak formulation for the structural equations, so that





in which G(X, t) is the Lagrangian elastic force density of the aortic wall, U(X, t) is the Lagrangian velocity field of the aortic wall, and V(X) is an arbitrary Lagrangian test function that is not assumed to vanish on [partial differential]V. In the continuum setting, eqs. (4) and (9)(10) are equivalent definitions for g(x, t); however, these formulations lead to different numerical schemes when discretized, and only eqs. (9)(10) lead directly to a standard nodal finite element structural discretization. Further, in the continuum equations, eqs. (11)(12) are equivalent to eq. (6) and to eq. (8). Further details of the numerical methods used in this work are described in a brief Appendix.

2.2 Valve Leaflets

2.2.1 Geometry

As in earlier work [34, 38], the leaflet geometry is determined by the mathematical theory of the fiber architecture of aortic valve leaflets developed by Peskin and McQueen [59]. This theoretical model of valve geometry and fiber architecture derives the shape of the leaflets from their function, which is to support a uniform pressure load during diastole, when the valve is closed. Each valve leaflet is composed of two families of fibers that are orthogonal to each other. One family of fibers runs from commissure to commissure, and the other runs from the bottom scalloped edge of the aortic sinus to the free edge of the valve leaflet. The leaflets are defined in a closed and loaded configuration; see fig. 1. In this configuration, the radius of each leaflet, measured from the tip of the valve to the end of the belly region, is 1.45 cm, and the coapting portion of each leaflet is 0.97 cm tall. The scale of the leaflets is based on measurements of human aortic roots [61,70]. The height of the coapting portion of the leaflet is chosen to be similar to that of the real valve while enabling the model valve to support a physiological pressure load when closed.

Fig. 1
Geometrically idealized model of the aortic root. We use Peskin and McQueen’s theoretical model of the collagen fiber architecture of aortic valve leaflets [59] to determine the shape of the leaflets, panel (a). We use measurements of human aortic ...

2.2.2 Mechanical response

The mechanical behavior of the leaflets is described by a strain-energy functional E = E[ϕ(·, t)] described previously [34, 38]. Briefly, let the curvilinear coordinates q = (q, r) be chosen so that a fixed value of q labels an individual fiber, and so that r runs along the fibers. Consequently, ([partial differential]ϕ/[partial differential]r)/‖[partial differential]ϕ/[partial differential]r‖ is the unit fiber tangent vector. The total elastic energy E is the sum of a stretching energy Es and a bending energy Eb,




Eq. (14) accounts for the total stretching energy associated with the fibers, in which x2130s is a spatially inhomogeneous local stretching energy with a quadratic length-tension relationship, and eq. (15) accounts for the total bending energy associated with the fibers, in which cb is a spatially inhomogeneous bending stiffness and ϕ¯ is the reference configuration, which is taken to be the initial configuration. Because the valve leaflets are modeled as thin elastic surfaces, the bending resistant energy allows the model to account for the thickness of the real valve leaflets. Larger values of cb model a thick, stenotic valve, whereas smaller values of cb model a thin, flexible valve. The resulting Lagrangian body force density is the Fréchet derivative of E.

Valve leaflet parameters are empirically determined for a fixed leaflet discretization by choosing the leaflet stiffness so that the leaflet supports a diastolic load of approximately 80 mmHg with 10% strain in the family of fibers running from commissure to commissure. This leads to a commissural fiber stiffness of 7.5e6 dyne/cm. We further assume that these commissural fibers are a factor of 10 stiffer than the radial fibers that run orthogonal to the commissural fibers [64]. The bending stiffness is chosen to be approximately the smallest value that yields valve leaflets that successfully coapt at the end of systole. These material parameters are chosen only to yield a functional valve, and are not represented as corresponding to those of an actual valve cusp. Nonetheless, this model of the valve leaflets does account for key structures of actual valve cusps, including the prominent collagen fibers that run from commissure to commissure and the 10-to-1 anisotropy of real valve leaflets [64]. See refs. 38 and 34 for further details, and see also the discussion of the limitations of this model in sec. 5. In preliminary studies, a sensitivity analysis was conducted for the bending stiffness of the aortic leaflets, and it was found that variations on the order of ± 25% in the bending stiffness have little effect on the physiological predictions of the models (not shown).

2.3 Aortic Wall

2.3.1 Geometry

An idealized model of the aortic root is constructed, as in previous work [34,38]. The dimensions of this model are based on measurements by Swanson and Clark [70] of human aortic roots collected at autopsy inflated to a diastolic pressure of 80 mmHg, and the geometry of the model is based on measurements by Reul et al. [61] from angiograms. The diameter of the aortic portion of the model is 3 cm, whereas the diameter of the sinus region, measured from a commissura to the center of a sinus, is 3.5 cm. The overall length of the model is of 10 cm, and the distance between the annulus and the aortic flow outlet is 7.75 cm. The thickness is constant throughout the model and is 2 mm [66]. In our simulations, we employ a semi-rigid model of the left-ventricular outflow tract while allowing for a fully flexible description of the aortic sinuses and the ascending aorta. The aortic annulus, which is the scalloped line of attachment between the valve leaflets and the aortic sinuses, was considered to be the lower boundary of the sinuses; see fig. 1. The aortic model used in this study differs from previous work [34,38] in that here we treat the aortic sinuses and ascending aorta as an incompressible hyperelastic material, whereas in earlier studies it was treated as a rigid structure.

Because the geometry of the valve leaflets is defined in a closed and loaded configuration, it is also necessary to specify the geometry of the vessel wall in a loaded configuration, and to compute the corresponding unloaded configuration. This configuration will depend on both the prescribed loaded configuration and the constitutive model.1 Because we use a constitutive model fit to data obtained from aortic tissue samples with a “J-shaped” stress-strain relationship (see fig. 2), accounting for initial deformations is also needed to obtain realistic deformations at systemic pressure loads. Specifically, if we did not account for these initial deformations, the sinuses and ascending aorta would experience unrealistically large deformations under even diastolic pressures. To determine the unloaded configuration, we employ an iterative method described by Bols et al. [5] detailed in sec. 2.3.3.

Fig. 2
Biaxial tensile test data obtained from human aortic root tissue samples by Azadani et al. [3] and our fit of the isotropic hyperelastic constitutive model (16) to these data.

2.3.2 Mechanical response

We describe the elasticity of the aortic sinuses and ascending aorta using a hyperelastic constitutive model fit to biaxial tensile test data collected by Azadani et al. [3] from tissue samples from human aortic sinuses. Because the data reported by Azadani et al. indicate a nearly isotropic material response (see fig. 2), in this work we model the material response of the aortic sinuses and ascending aorta using an isotropic strain-energy functional W with an exponential stress-strain relation,


in which I1 = I1(C) = tr(C) is the first invariant of the right Cauchy-Green strain tensor C = FT F, and F = [partial differential]χ/[partial differential]X is the deformation gradient tensor associated with the deformation mapping χ : (V, t) [mapsto] Ω. Material parameters c and b were obtained by least-squares fits to experimental data using MATLAB (Mathworks, Inc., Natick, MA, USA). The first Piola-Kirchhoff stress tensor P is determined from W via


in which ps is a structural pressure-like term that is chosen to improve the accuracy of the stress predictions of the method [29]. As in earlier work [29], we compute ps via


in which I3 = I3(C) = det(C) is the third invariant of C and βs = 2.5e4 kPa. Thus, P = 0 for F = I, and the structural model provides an energetic penalty for any compressible deformations.

2.3.3 Backward displacement method

To determine the unloaded configuration of the vessel wall, we use an iterative backward displacement method described by Bols et al. [5] implemented using a custom MATLAB script that interfaces with the ABAQUS (Dassault Systèmes Simulia Corp., Providence, RI, USA) finite element analysis software. Because the constitutive model (16) is not provided by ABAQUS, this model was implemented using a UHYPER subroutine.

Let xm denote positions of the finite element mesh nodes in the prescribed loaded configuration, let Xmi denote the corresponding nodal positions in the computed reference configuration after i iterations of the backward displacement algorithm, and let χmi=χm[Xi] denote nodal positions in the computed loaded configuration when Xi is used as the unloaded reference configuration. Notice that xm, Xmi, and χmi are all defined at the nodes of the finite element mesh.

The algorithm is straightforward: First, we initialize X0 := x (i.e., we use the deformed coordinates as an initial guess for the unloaded configuration). Then, given Xi, Xi+1 is determined by


Notice that if we define the displacement from the computed reference configuration to the deformed coordinates at step i by Di = χiXi, then (19) is equivalent to


Thus, the reference coordinates are determined via a backward displacement from the deformed configuration. If this process converges, χix and we obtain reference coordinates X such that χ[X] = x.

The convergence of this iterative process is assessed in terms of the discrete L2 norm of the residual xχi, i.e.


The routine was judged to have converged once the residual was less than 0.1% of the average vessel radius and the change in the residual between two consecutive iterations was less than 0.015% of the average radius. In general, obtaining convergence may require the use of damping, or of gradually increasing the loading pressure; however, in practice, we have found this algorithm to be robust, and neither damping nor incremental loading were required for the analyses performed herein.

2.4 Driving and Loading Conditions

A left-ventricular pressure waveform adapted from human clinical data of Murgo et al. [55] is prescribed at the upstream inlet of the left-ventricular outflow tract to drive flow through the model aortic root, and downstream loading conditions are provided by a three-element Windkessel model by Stergiopolus et al. [69] fit to the clinical data of Murgo et al. [55]. As in earlier work [34,38], zero pressure boundary conditions are imposed on the boundary of the fluid-filled region exterior to the aortic root model. Coupling between the detailed IB model and the reduced circulation models is also done as in earlier work [34,38]. Along the proximal and distal ends of the vessel, zero-displacement boundary conditions are imposed.

Because the normal traction boundary conditions are used at both the inlet and the outlet are subject to the development of flow instabilities [24], a small amount of flow stabilization is applied in the fluid region within a distance of approximately 2h from the boundary. Specifically, the net flow rate through the upstream or downstream boundary is first determined, and then any local flow that is opposite from the direction of net flow is penalized by including additional forcing that acts to avoid localized regions of flow reversal on the boundary. The same basic approach is used both for inflow and outflow, and it permits the boundary conditions to switch from inflow to outflow conditions based on the flow within the interior of the detailed model, and on the pressure differences across that model. This stabilization scheme is similar in its physical mechanism to that described by Esmaily Moghadam [24], but here we employ interior forcing rather than directly modifying the formulation of the physical boundary condition.

2.5 Discretization

In our dynamic fluid-structure interaction simulations, the computational domain is taken to be a 10 cm × 10 cm × 10 cm cubic region discretized using a three-level locally refined grid. For the purposes of a mesh convergence study, two different fine grid spacings are used: a coarser one corresponding to h = 0.78125 mm, which is equivalent to a uniform 128 × 128 × 128 Eulerian discretization, and a finer one corresponding to h = 0.390625 mm, which is equivalent to a uniform 256 × 256 × 256 Eulerian discretization. The curvilinear mesh used to discretize the valve leaflets has a physical grid spacing of approximately h/4 for the coarser Eulerian discretization, and to h/2 for the finer discretization. The volumetric mesh used to discretize the vessel wall uses 23,400 8-node (trilinear) hexahedral elements for the aortic root model, which results in an average mesh aspect ratio of 1.93 ± 0.25. A preliminary convergence study in ABAQUS was performed to ensure mesh convergence of the structural model for both the backward displacement method and the forward FE analysis. We use the standard four-point delta function of Peskin [58] for the coarser simulations, and a broadened 8-point version of this kernel function for the higher resolution simulations, so that the physical extent of the regularized delta functions is the same for both spatial resolutions.2 Consequently, the valve leaflets have the same effective thickness with respect to the fluid dynamics part of the simulation for both spatial resolutions. We remark that the structural response of the leaflets is primarily determined by the structural model and not by the effective leaflet thickness “seen” by the fluid model.

A uniform time step size of Δt = 9.94369e-6 s was used with the coarser Eulerian discretization, whereas a uniform time step size of Δt = 1.40625e-5 was used for the finer Eulerian discretization. These time step sizes are empirically determined to be within a factor of 2 of the largest time step that satisfies the stability restriction of our explicit time integration method. CFL numbers are on the order of 0.03 during systole and 0.003 during diastole. Thus, the time step size is generally substantially smaller than the maximum stable time step size allowed by our fluid solver, for which the CFL number is order 1. For this model, the stability is primarily determined by the explicit coupling between the solid and the Eulerian momentum equation. Asymptotically, we expect that for this model, Δt ~ h4 because of the presence of bending-resistant elastic elements in the model valve leaflets.

3 Results

3.1 Constitutive Parameters

The constitutive model (16) used to describe the mechanical response of the aortic sinuses and ascending aorta is fit to biaxial tensile tests of human aortic root tissue samples reported by Azadani et al. [3], yielding model parameter values of c = 12.8 kPa and b = 6.9; see fig. 2. Recall that these material parameters are used in the aortic sinuses and ascending aorta, whereas the outflow tract is treated as an essentially rigid structure.

3.2 Unloaded Geometry

The unloaded configuration is determined by assuming that the loaded configuration corresponds to a pressure load of 80 mmHg. In fig. 3, panel (a) shows the prescribed loaded geometry, and panel (b) shows the computed zero-pressure configuration. Fig. 4 shows the convergence history of the backward displacement method. When the computed zero-pressure geometry is subjected to a 80 mmHg pressure load, the deformed geometry agrees with the prescribed geometry to within approximately 15 µm; see fig. 5. The discrepancy between the prescribed loaded geometry and the computed loaded geometry is greatest at the commissures and smallest along the straight portion of the vessel.

Fig. 3
Views of the (a) loaded and (b) unloaded geometry of the aortic root.
Fig. 4
Convergence history of the inverse displacement procedure.
Fig. 5
Contour plots showing the norm of the distance between the prescribed and computed loaded geometries. Panels show a side view (a) and a section view (b) of the model used in the study. Recall that the unloaded geometry is determined from the loaded geometry. The ...

3.3 Fluid-Structure Interaction Analysis

A dynamic fluid-structure interaction analysis is performed that includes four cardiac cycles, including an initialization cycle and three cycles in which the model has attained periodic steady state. (Data are reported only for the beats in which the model has reached steady state.) Fig. 6 shows the prescribed left-ventricular driving pressure along with the computed aortic loading pressure determined by the coupled Windkessel model and the computed flow rate through the symmetric aortic root. We emphasize that the flow rate is not imposed in the model, but rather is predicted by the fluid-structure interaction analysis.

Fig. 6
Results of a dynamic fluid-structure interaction analysis of the aortic root over the final three simulated cardiac cycles using a coarser Cartesian grid (effective grid resolution of N = 128) and a finer Cartesian grid (effective grid resolution of ...

For the simulations using the coarser Cartesian spacing, h = 0.78125 mm, stroke volume is 96.2 ml in average, peak flow rate is 591.5 ml/s, and cardiac output is 6.5 l/min. For the finer Cartesian grid spacing, h = 0.3906 mm, stroke volume is 100.2 ml in average, peak flow rate is 643.8 ml/s, and cardiac output is 6.8 l/min. Both sets of simulation results are within 15% of the hemodynamic parameters of the subject-specific clinical data used to determine the driving pressure and the Windkessel model parameters. Specifically, the driving pressure and the Windkessel model parameters used in this work were based on clinical measurements from a subject with stroke volume of 100 ml, peak flow rate of 560 ml/s, and cardiac output of 6.8 l/min, as reported by Murgo et al. [55]. The simulations also capture the incisura in the flow profile corresponding to the backward flow generated by the leaflets. The backward flow volume is 1.77 ± 0.11 ml, or 1.8% of the stroke volume, for the coarser grid spacing, and 1.52 ± 0.059 ml, or 1.5% of the stroke volume, for the finer grid spacing. The net forward flow volume (i.e., into the aorta) following leaflet coaptation is 0.02 ± 0.01 ml. This forward flow results from interactions between the elastic aortic root and the dynamic downstream loading pressure. As the downstream pressure falls during diastole, the load on the aortic wall and leaflets decreases and allows the aortic root to gradually push fluid into the distal aorta.

The model captures the opening and closing dynamics of the valve at both Cartesian grid spacings. Fig. 7 shows the leaflets being pushed apart along with the formation of the vortices that help the leaflets to close at the end of systole. Fig. 8 shows that although there is some net flow towards the left ventricle as the valve closes, once the valve closes, no further flow is allowed between the aorta and the left ventricle.

Fig. 7
Valve leaflet configurations during valve opening shown at equally spaced times during the fourth simulated cardiac cycle for (a) the coarser Cartesian grid (effective grid resolution of N = 128) and (b) the finer Cartesian grid (effective grid resolution ...
Fig. 8
Similar to fig. 7, but here showing the dynamics of the closure of the aortic valve during the fourth cycle. In this case, the solid lines indicate the instants shown in the figure panels.

Fig. 9 shows the distribution of the structural displacements in the aortic root, and fig. 10 shows displacement as a function of time at four selected locations: the center of the sinus; the sinotubular junction; the commissure; and the ascending aorta. Notice that after the first cardiac cycle, the displacements are essentially periodic in time. Cartesian grid resolution has a minimal effect on the range of displacements except for the displacements measured at the commissure. The largest amplitude is observed in the ascending aorta (0.58 ± 0.05 mm) and at the sinotubular junction (0.52±0.02 mm). A measure of aortic stiffness is the radial pulsation, i.e., the ratio of the difference between the systolic radius and the diastolic radius to the average radius [56], which in the present analysis was approximately 3.9% for the coarser Cartesian grid spacing and 2.9% for the finer Cartesian grid spacing. Notice that radial pulsation is a highly sensitive measure of deformation, and the different values obtained for the coarse and fine simulations indicate a relatively small difference in the displacements of approximately 0.175 mm, which is less than 0.5h on the highest resolution Cartesian grid used in this study. These values are similar to the value of 2.5% reported in refs. 41 and 56 for the ascending aorta. The distribution of the displacements was also in good agreement with literature data, with the displacements in the aortic root greater in the sinotubular junction, and at the commissures than in the sinuses, as reported also by Lansac [45].

Fig. 9
Displacement contours along the aortic root model. Panel (a) shows displacement contours along the axial section of the aortic wall during valve opening, shown as the same time instants as in fig. 7, for the coarser Cartesian grid spacing. Panel (b) is similar, ...
Fig. 10
Displacements of the deformable aortic root model at four locations: (a) the straight portion of the aorta; (b) just above one of the commissures; (c) the sinotubular junction; and (d) the middle of the sinus.

The distribution of the maximum principal stress in the aortic root follows a distribution similar to that of the displacement, with the maximum stresses concentrated in the sinotubular junction, just above the sinuses, as shown in fig. 11. Stress in the aortic root follows the pulsatile waveform of the pressure; see fig. 12. Both Cartesian grid spacings yield similar stress distributions; compare fig. 11 panels (a) and (b), but fig. 12 shows some differences along the sinotubular junction and within the sinuses. The maximum principal stress in the sinuses is 20% of that observed in the rest of the aortic root. The distribution, and the range of the maximum principal stress predicted by our model, is in agreement with other analyses available in literature, in particular with the 220 kPa peak of the maximum principal stress in the sinotubular junction reported by Conti et al. [14].

Fig. 11
Maximum principal stress contours along the aortic root model. Panel (a) shows displacement contours along the axial section of the aortic wall during valve opening, shown as the same time instants as in fig. 7, for the coarser Cartesian grid spacing. ...
Fig. 12
Maximum principal stress of the deformable aortic root model at four locations: (a) the straight portion of the aorta; (b) just above one of the commissures; (c) the sinotubular junction; and (d) the middle of the sinus.

4 Discussion

In this work, previous IB models of aortic valve dynamics [34,38] are substantially extended by incorporating a deformable model of the aortic sinuses and ascending aorta. This extended model is thereby able to more closely replicate the complex in vivo dynamics of the real aortic root. As such, it promises to serve as a powerful model for simulations of aortic root behavior and valve dynamics in mechanobiology studies or in the design, testing, or selection of medical devices. The dynamic analysis of the aortic valve included four complete cardiac cycles, including a single initialization cycle, which are shown to be sufficient to yield essentially periodic results. As in previous IB models of aortic valve dynamics [34, 38], key hemodynamic parameters, including aortic pressure, maximum flow rate, cardiac output, and stroke volume, are in reasonable agreement with physiological ranges [22, 76] and are also in good agreement with the patient-specific data [55] used to fit the Windkessel model [69] employed in this study. Large oscillations in the flow rate, as well as in the deformations, indicate the reverberations of the aortic valve leaflets upon closure, and the resonances of the elastic model of the sinuses and ascending aorta; i.e., these are the so-called heart sounds that are generated by the aortic valve. Specifically, these dynamics are not the result of the explicit time coupling scheme used in the numerical method. Further, grid convergence results demonstrate for the first time that the IB method is able to yield essentially grid-resolved bulk flow parameters and vessel deformations at practical grid spacings.

Many strategies have been developed to determine the zero-pressure geometry of blood vessels [1, 5, 20, 31, 65, 73]. Indeed, medical images of arteries always provide the vessel geometry in a loaded configuration. It is known, however, that using the imaged geometry as the reference geometry of the vessel underestimates the amount of strain present in the blood vessel [20]. Among the approaches developed to estimate the zero-pressure geometry of blood vessels, we chose an iterative backward displacement approach, as described by Bols et al. [5] and by Sellier [65]. Iterative approaches based on geometric considerations are straightforward to implement and can be applied to a variety of problems and solvers because they require minimal implementation effort. In this work, the backward displacement analysis was implemented in ABAQUS and converged to a zero-pressure geometry that, when loaded, yielded a deformed configuration that is in good agreement with the initially prescribed loaded geometry; see figs. 35. However, the deflated geometry is quite different from the configuration assumed by a deflated aortic root harvested from a body; see fig. 3. This is a drawback of these methods, as shown also in a recent study by Wittek et al. [77], in which the unloaded configuration of the aorta is different from an unloaded aorta harvested from a body. This mismatch could be the result of many factors, such as the choice of the boundary conditions or the need for adding the residual stress to balance the level of stress within the aortic wall [28]. We expect that the addition of residual stress estimates, planned in future work, will reduce the mismatch between the deflated geometry computed and the deflated aortic root observed in vivo.

Interestingly, the addition of a compliant aortic wall did not substantially alter the global hemodynamic parameters previously determined by a noncompliant IB aortic root model [34]; see fig. 13, which shows that the peak pressure attained by the noncompliant model is slightly higher than that obtained in the compliant model, whereas the flow rates are similar in both models. This finding might suggest that our vessel wall model is overly stiff despite our use of a constitutive model fit to biaxial tensile test data from human aortic tissues samples. Another reason for the reduced motility of the aortic root model is the fact that the outflow tract of the left ventricle is fixed. In physiological conditions, due to ventricular contraction and relaxation, the aortic root moves upwards and sideways with every heartbeat. In this work, we constrained this motion, thus restricting the dynamics of the aortic root model. This constraint reduced greatly the accelerations the aortic root is subjected to as compared to the physiological conditions.

Fig. 13
Results for the first simulated cardiac cycle are compared to results from previous work that used a rigid model of the aortic root [34] (shown with grey dashed lines). We remark that the noncompliant aortic root model of ref. [34] employed a different ...

In contrast, we observed that using a compliant aorta affected more clearly the dynamics of valve leaflet closure. Figs. 7(b) and 8(b) show the flow inversion and the distension in the sinus, which collects fluid during systolic ejection, and the elastic recoil of the sinus, which pushes down on the leaflet during closure. The addition of this mechanism made the model substantially more robust. Specifically, in previous work [34] we found that the rigid aortic root model required the valve leaflet properties to be very carefully tuned in order for the model to remain competent throughout the cardiac cycle. This is in clear contrast to the present compliant aortic root model, which yields a competent valve for a relatively broad range of leaflet parameters (data not shown). A compliant root provides elastic recoil to “push” the leaflets closed and also makes it easier for the nearly-closed leaflets to deform to coapt fully and to prevent regurgitation. This is a trait that makes our model more relevant physiologically, as aortic valves are continuously remodeled and keep functioning in vivo for varying mechanical properties and configurations. Future work will address the impact that geometrical characteristics of the ascending aorta, such as length, curvature, and tortuosity, have on the FSI simulations. We expect that these changes will not substantially alter the bulk flow parameters calculated in this study, but we anticipate that changes in the shape of the ascending aorta will affect blood velocity distributions and will generate flow pattern qualitatively closer to those seen in vivo.

Similar model predictions were obtained for both Cartesian grid spacings used in this study. We observed hemodynamic values that were within 5% of the clinical values reported by Murgo et al. [55] except for the peak flow rate, which was within 5% of the values reported by Murgo et al. on the coarser Cartesian grid and within 15% on the finer grid. Cardiac output and stroke volume were both in very good agreement with the clinical data, within 4% on the coarser Cartesian grid, and within 1% on the finer grid. In addition, the pulsatile radius, a quantity used to measure aortic stiffness, was also shown to be in relatively good agreement with literature values in both cases [41].

5 Limitations and Future Work

Work is underway to address some of the limitations of these models. As described previously, our present models do not account for the curved geometry of the aortic root, or for the residual stresses present in real arterial vessels. We specifically aim to incorporate a description of residual stresses and realistic medical imaging-derived anatomic geometries into our IB models. Further, the present models use an isotropic model of the aortic sinuses and ascending aorta. Although the experimental data that were used to fit this model show an essentially isotropic material response, it is known that real aortic tissues have a well-defined collagen fiber structure with an anisotropic material response. We also aim to incorporate a more realistic hyperelastic model of the aortic wall that accounts for families of collagen fibers [30] as well as experimentally constrained models of the elasticity of the aortic valve leaflets [50]. Further, it is now well known that the valve leaflets are in fact layered structures [63]. These details are not accounted for in the present model, and we anticipate that including them in future work will require an important extension of the methods presented in this study. Specifically, because the valve leaflets are thin elastic structures, we expect that using realistic hyperelastic continuum models to describe the mechanics of the valve leaflets will require adopting a nonlinear shell formulation. To our knowledge, continuum shell models have not yet been widely used within the framework of the IB method, and we anticipate additional computational research will be required to implement these types of formulations into our IB simulation framework.

Another limitation of this work is the long run times required by these simulations. Because of the small time step sizes required by the present algorithm, each cardiac cycle required 50,000 or more time steps, requiring on the order of one day per cycle on 32–64 cores of a modern high-performance computing (HPC) cluster for the coarse simulations and on the order of one week per cycle for the fine-grid simulations. Higher-resolution simulations using the present methods and implementation are computationally infeasible at this time. We are currently working to develop effective multigrid solvers for stable implicit time stepping schemes for the IB method [40]. These methods promise to allow for substantially larger time step sizes as well as shorter run times at current grid resolutions. They could also enable higher resolution simulations that promise to resolve the fine-scale details eliminate the remaining grid-sensitivity observed in the kinematics of the valve leaflets and the distensibility of the aortic root.

6 Conclusions

This work has presented a new fluid-structure interaction model of the aortic root that extends earlier IB models of the aortic valve by incorporating a finite-strain model of the aortic root that uses a constitutive model fit to tensile test data obtained from human aortic root tissue specimens. This model captures both the complex fluid dynamics of the flow and the finite deformations of the vessel wall, and it yields results that are in good agreement with the clinical data that were used to fit the circulation model used in this study, as well as to physiological data in the research literature. Finally, a grid convergence study demonstrates that nearly grid converged results are obtained at practical spatial resolutions, although yet higher spatial resolution, or higher accuracy at present resolutions, is needed to obtain fully resolved simulation results.


This work has been funded in part by American Heart Association award 10SDG4320049, National Institutes of Health awards GM071558 and HL117063, and National Science Foundation awards DMS 1016554/1460368 (to NYU School of Medicine and UNC-Chapel Hill) and ACI 1047734/1460334 (to NYU School of Medicine and UNC-Chapel Hill). VF was supported in part by the Tandon School of Engineering of New York University. Computations were performed at New York University using computer facilities funded in large part by a generous donation by St. Jude Medical, Inc.


A Numerical Methods

We summarize the key discretization methods used in this study in this appendix. We use a locally-refined staggered-grid discretization of the Eulerian equations along with a finite difference-based discretization of the fiber model of the valve leaflets and a finite element-based description of the continuum vessel wall model. Our approach is similar to the spatially adaptive IB scheme described previously [34], except that here we also employ a finite element-based description of the aortic wall, using an approach introduced by Griffith and Luo [39].

A.1 Eulerian and Lagrangian spatial discretizations

The physical domain Ω is discretized using a locally-refined Cartesian grid, but for simplicity we describe only the uniform-grid version of this method; details on the adaptively refined discretization are provided by ref. 34. Let (i, j, k) index the cells of the Cartesian grid, and let xi,j,k=((i+12)h,(j+12)h,(k+12)h) indicate the position of the center of grid cell (i, j, k), in which h is the Cartesian grid meshwidth and Δx = h3 is the Cartesian grid cell volume. The Eulerian velocity field u = (u, υ, w) is approximated at the center of each face of the Cartesian grid cells in terms of the velocity component that is normal to that face, so that u is approximated at the locations xi12,j,k, υ is approximated at the locations xi,j12,k, and w is approximated at xi,j,k12. The Eulerian force densities f and g are approximated in the same staggered-grid fashion. The Eulerian pressure p is approximated at the centers of the Cartesian grid cells.

The deformations and forces associated with the valve leaflets are approximated on a fiber-aligned curvilinear mesh. Let l index the nodes of this mesh, and let ϕl and Fl indicate the current position and Lagrangian force density of node l, respectively, and let Δql indicate the area fraction (quadrature weight) associated with node l. Fl is computed from ϕl using a finite difference approximation to the fiber force density; see refs. 38 and 34 for details, including relevant finite difference formulae.

The deformations, stresses, and resultant forces associated with the aortic wall are approximated on a volumetric Lagrangian mesh. Let e index the elements of this mesh, with Ve indicating the volume associated with element e, so that V = [union or logical sum]eVe, let m index the mesh nodes, and let [var phi]m (X) indicate the interpolatory Lagrangian finite element basis function associated with node m. The structural deformation and elastic force density are approximated in a standard manner via



in which χm(t) are the nodal positions and Gm(t) are the nodal force densities. The deformation gradient F is computed by directly differentiating the approximation to χ(X, t), and P is evaluated from the approximation to F. The nodal values Gm(t) are determined so that


for all n. This leads to a linear system of equations that defines the nodal force densities in terms of the nodal deformations. In practice, these integrals are evaluated element-by-element using Gaussian quadrature rules that are exact for the left-hand side but are generally only approximate for the right-hand side.

A.2 Lagrangian-Eulerian interaction

To couple the Lagrangian and Eulerian discretizations, we employ approximations to the integral transforms of the continuum equations that replace the singular Dirac delta function δ(x) = δ(x) δ(y) δ(z) by a regularized delta function δh(x) = δh(x) δh(y) δh(z). In our computations, we construct the three-dimensional regularized delta function either by using the four-point one-dimensional regularized delta function of Peskin [58], or by using a broadened version of this function that has a spatial extent of 8 meshwidths. The same regularized delta function is used for both the thin model of the valve leaflets and the volumetric model of the vessel wall.

For the leaflet model, we employ a Lagrangian-Eulerian coupling approach frequently used with the IB method. The Lagrangian forces F = (Fx, Fy, Fz) associated with the leaflets are converted into equivalent Eulerian forces f = (fx, fy, fz) via




This amounts to using the trapezoidal rule to discretize the integral transform (3). We employ the compact notation


to express this force-spreading operation. The Eulerian velocity u = (u, υ, w) is used to determine the dynamics of the physical positions ϕ = (ϕx, ϕy, ϕz) of the material points of the valve leaflets via




We use the notation


to express this velocity-restriction operation. Notice that the force-spreading and velocity-interpolation operators are adjoints, i.e., Rl[ϕ] = Sl[ϕ]*.

The Lagrangian-Eulerian coupling scheme used for the vessel model is based on an approach developed by Griffith and Luo [39]. This approach is first to construct a force-spreading operator, and then to determine a velocity-restriction operator that is the adjoint of the force-spreading operator. Briefly, for each element e, quadrature points Q(e) are determined, and for each quadrature point q [set membership] Q(e), Xqe indicates the reference coordinates of the quadrature point, and ωqe indicates the quadrature weight associated. The Lagrangian force density G = (Gx, Gy, Gz) is converted into the equivalent Eulerian force density g = (gx, gy, gz) by discretizing eq. (9), with δ(x) replaced by δh(x), using this quadrature rule, i.e.




Notice that these formulae take advantage of the fact that we may evaluate the approximations to G and χ at arbitrary Lagrangian locations. Specifically, we are not restricted to evaluating G and χ at the nodes of the finite element mesh. In practice, the quadrature rules are dynamically determined to ensure that the physical spacing of the quadrature points is on average one-half of a Cartesian meshwidth, even under the presence of extremely large structural deformations. We use the notation


to define the force-spreading operator. To determine the corresponding velocity-restriction operator, we introduce a Lagrangian velocity field U = (U, V, W), which is defined in terms of nodal velocities Um via


This velocity field is required to satisfy, for all n,




This yields a system of linear equations to be solved for the nodal values of U. Notice that eqs. (38)(40) correspond to a discrete component-wise approximation to eq. (12). We use the notation


Notice that Rw[χ] is a nonlocal operator in the sense that it requires the solution of a system of linear equations. Although not shown here, it is the case that Rw[χ] = Sw[χ]*. See Griffith and Luo [39] for further discussion.

For the models of the leaflets and for the vessel wall, the Lagrangian-Eulerian coupling scheme is designed to conserve force and torque [39,58]. Moreover, because adjoint operators are used to spread force and to interpolate velocity, the scheme conserves energy at least in its semi-discrete (continuous time but discrete space) formulation.

A.3 Time stepping

The time stepping scheme is similar to that of Griffith and Luo [39]. Let Δt indicate the (uniform) time step size, and let [tn, tn+1] = [nΔt, (n + 1)Δt] indicate the nth time interval. In this subsection, superscript indices always indicate the time step number. The state variables u, ϕ, and χ are defined at integer time steps. The pressure p, which is not a state variable of the system, is defined at half-integer time steps. Approximations to the state variables and to derived quantities such as the Lagrangian forces are also evaluated at half-integer time steps.

The time stepping scheme proceeds as follows. First, predicted intermediate approximations to the structural deformations ϕ and χ at time tn+12 are determined via



This is an application of the forward Euler time stepping scheme. Lagrangian force densities Fn+12 and Gn+12 are determined from the predicted structure configurations, and these Lagrangian force densities are spread to the Cartesian grid via



Next, we solve the discretized incompressible Navier-Stokes equations for un+1 and pn+12 via a Crank-Nicolson-Adams-Bashforth scheme,





In the uniform grid case, the discrete operators [nabla]h, [nabla]h ·, and h2 are standard second-order accurate staggered-grid finite difference approximations to the gradient, divergence, and Laplace operators, respectively [33]. We compute u · [nabla]hu via a staggered-grid version of the piecewise parabolic method (PPM) [33]. On locally-refined Cartesian grids, these discrete operators are straightforward extensions of the uniform grid operators [34]. Solving eqs. (46) and (47) for un+1 and pn+12 requires the solution of the time-dependent Stokes equations. The Stokes equations are solved using an iterative Krylov method with a multigrid preconditioner derived from the classical projection method [9, 33].

Finally, we determine approximations to the structural deformations at the end of the time step via





This is an application of the explicit midpoint rule to the structural dynamics.

In the initial time step, we do not have a value for un−1 as required to evaluate the Adams-Bashforth approximation to the convective term. Consequently, in the initial time step, we employ a two-step Runge-Kutta scheme.

Within each time step, coupling between the detailed fluid-structure interaction model and the reduced Windkessel model of the circulation is done as described previously [34,38]. In this approach, we perform only a single Stokes solve per time step, except in the initial time step, when we perform two Stokes solves within the context of a two-step Runge-Kutta scheme.


1Keeping the loaded configuration fixed, if we assume a larger initial pressure, then the result will be larger deformations from the computed unloaded configuration to the prescribed loaded configuration, thereby leading to a stiffer model at the prescribed loaded configuration. If we assume a smaller initial pressure, then there will be smaller strains from the unloaded to the loaded configuration, and thus a more compliant model at the target configuration.

2Similar convergence studies, although for substantially different applications, are described in refs. 37 and 29.

Contributor Information

Vittoria Flamini, Department of Mechanical and Aerospace Engineering, New York University Tandon School of Engineering, Brooklyn, New York, USA.

Abe DeAnda, Division of Cardiothoracic Surgery, Department of Surgery, University of Texas Medical Branch, Galveston, TX.

Boyce E. Griffith, Departments of Mathematics and Biomedical Engineering and McAllister Heart Institute, Phillips Hall, Campus Box 3250, University of North Carolina, Chapel Hill, North Carolina, USA, Phone: (919) 962-7110, ude.cnu.liame@gecyob.


1. Alastrué V, Garía A, Peña E, Rodríguez J, Martínez M, Doblaré M. Numerical framework for patient-specific computational modelling of vascular tissue. Int J Numer Method Biomed Eng. 2010;26(1):35–51.
2. Alastrué V, Peña E, Martínez MÁ, Doblaré M. Assessing the use of the opening angle method to enforce residual stresses in patient-specific arteries. Ann Biomed Eng. 2007;35(10):1821–1837. [PubMed]
3. Azadani AN, Chitsaz S, Matthews PB, Jaussaud N, Leung J, Tsinman T, Ge L, Tseng EE. Comparison of mechanical properties of human ascending aorta and aortic sinuses. Ann Thorac Surg. 2012;93(1):87–94. [PubMed]
4. Bellhouse BJ, Bellhouse FH. Mechanism of closure of the aortic valve. Nature. 1968;217:86–87. [PubMed]
5. Bols J, Degroote J, Trachet B, Verhegghe B, Segers P, Vierendeels J. A computational method to assess in vivo stresses and unloaded configuration of patient-specific blood vessels. J Comput Appl Math. 2012
6. Bonow RO, Carabello BA, Chatterjee K, de Leon AC, Faxon DP, Freed MD, Gaasch WH, Lytle BW, Nishimura RA, O’Gara PT, O’Rourke RA, Otto CM, Shah PM, Shanewise JS. ACC/AHA 2006 guidelines for the management of patients with valvular heart disease. Circulation. 2006;114(5):E84–E231. [PubMed]
7. Borazjani I. Fluidstructure interaction, immersed boundary-finite element method simulations of bio-prosthetic heart valves. Comput Meth Appl Mech Eng. 2013;257:103–116.
8. Borazjani I, Ge L, Sotiropoulos F. Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. J Comput Phys. 2008;227(16):7587–7620. [PMC free article] [PubMed]
9. Cai M, Nonaka A, Bell J, Griffith B, Donev A. Efficient variable-coefficient finite-volume Stokes solvers. Submitted.
10. Cardamone L, Valentín A, Eberth JF, Humphrey JD. Origin of axial prestretch and residual stress in arteries. Biomech Model Mechanobiol. 2009;8(6):431–446. [PMC free article] [PubMed]
11. Carr JA, Savage EB. Aortic valve repair for aortic insufficiency in adults: A contemporary review and comparison with replacement techniques. Eur J Cardio Thorac Surg. 2005;25(1):6–15. [PubMed]
12. Cheng A, Dagum P, Miller DC. Aortic root dynamics and surgery: from craft to science. Phil Trans R Soc B. 2007;362(1484):1407–1419. [PMC free article] [PubMed]
13. Cheng R, Lai YG, Chandran KB. Three-dimensional fluid-structure interaction simulation of bileaflet mechanical heart valve flow dynamics. Ann Biomed Eng. 2004;32(11):1471–1483. [PMC free article] [PubMed]
14. Conti CA, Votta E, Della Corte A, Del Viscovo L, Bancone C, Cotrufo M, Redaelli A. Dynamic finite element analysis of the aortic root from mri-derived parameters. Med Eng Phys. 2010;32(2):212–221. [PubMed]
15. Creane A, Kelly DJ, Lally C. Patient specific computational modeling in cardiovascular mechanics. In: Lopez BC, Peña E, editors. Patient-Specific Computational Modeling. Springer; 2012. pp. 61–79.
16. Croft LR, Mofrad MRK. Computational modeling of aortic heart valves. In: De S, Guilak F, Mofrad MRK, editors. Computational Modeling in Biomechanics. Springer; 2010. pp. 221–252.
17. Crosetto P, Reymond P, Deparis S, Kontaxakis D, Stergiopulos N, Quarteroni A. Fluid–structure interaction simulation of aortic blood flow. Comput Fluid. 2011;43(1):46–57.
18. Dagum P, Green GR, Nistal FJ, Daughteres GT, Timek TA, Foppiano LE, Bolger AF, Ingels NB, Miller DC. Deformational dynamics of the aortic root: modes and physiologic determinants. Circulation. 1999;100(2):II54–II62. [PubMed]
19. De Hart J, Baaijens FPT, Peters GWM, Schreurs PJG. A computational fluid-structure interaction analysis of a fiber-reinforced stentless aortic valve. J Biomech. 2003;36(5):699–712. [PubMed]
20. de Putter S, Wolters BJBM, Rutten MCM, Breeuwer M, Gerritsen FA, van de Vosse FN. Patient-specific initial wall stress in abdominal aortic aneurysms with a backward incremental method. J Biomech. 2007;40(5):1081–1090. [PubMed]
21. Delfino A, Stergiopulos N, Moore JE, Meister JJ. Residual strain effects on the stress field in a thick wall finite element model of the human carotid bifurcation. J Biomech. 1997;30(8):777–786. [PubMed]
22. Driscol TE, Eckstein RW. Systolic pressure gradients across the aortic valve and in the ascending aorta. Am J Physiol. 1965;209(3):557–563. [PubMed]
23. Dumont K, Stijnen JMA, Vierendeels J, Van De Vosse FN, Verdonck PR. Validation of a fluid–structure interaction model of a heart valve using the dynamic mesh method in fluent. Comput Meth Biomech Biomed Eng. 2004;7(3):139–146. [PubMed]
24. Esmaily Moghadam M, Bazilevs Y, Hsia TY, Vignon-Clementel IE, Marsden AL. A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations. Comput Mech. 2011;48(3):277–291.
25. Fai TG, Griffith BE, Mori Y, Peskin CS. Immersed boundary method for variable viscosity and variable density problems using fast linear solvers. I: Numerical method and results. SIAM J Sci Comput. 2013;35(5):B1132–B1161.
26. Freeman RV, Otto CM. Spectrum of calcific aortic valve disease: Pathogenesis, disease progression, and treatment strategies. Circulation. 2005;111(24):3316–3326. [PubMed]
27. Fung YC. What principle governs the stress distribution in living organs? In: Fung YC, Fukada E, Wang J, editors. Biomechanics in China. Japan and USA: Science Press; 1983.
28. Fung YC. What are the residual stresses doing in our blood vessels? Ann Biomed Eng. 1991;19(3):237–249. [PubMed]
29. Gao H, Wang HM, Berry C, Luo XY, Griffith BE. Quasi-static image-based immersed boundary-finite element model of left ventricle under diastolic loading. Int J Numer Meth Biomed Eng. 2014;30(11):1199–1222. [PMC free article] [PubMed]
30. Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J R Soc Interface. 2006;3(6):15–35. [PMC free article] [PubMed]
31. Gee MW, Reeps C, Eckstein HH, Wall WA. Prestressing in finite deformation abdominal aortic aneurysm simulation. J Biomech. 2009;42(11):1732–1739. [PubMed]
32. Govindjee S, Mihalic PA. Computational methods for inverse deformations in quasi-incompressible finite elasticity. Int J Numer Method Eng. 1998;43(5):821–838.
33. Griffith BE. An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as a preconditioner. J Comput Phys. 2009;228(20):7565–7595.
34. Griffith BE. Immersed boundary model of aortic heart valve dynamics with physiological driving and loading conditions. Int J Numer Method Biomed Eng. 2012;28(3):317–345. [PubMed]
35. Griffith BE, Hornung RD, McQueen DM, Peskin CS. An adaptive, formally second order accurate version of the immersed boundary method. J Comput Phys. 2007;223(1):10–49.
36. Griffith BE, Hornung RD, McQueen DM, Peskin CS. Parallel and adaptive simulation of cardiac fluid dynamics. In: Parashar M, Li X, editors. Advanced Computational Infrastructures for Parallel and Distributed Adaptive Applications. Hoboken, NJ, USA: John Wiley and Sons; 2009.
37. Griffith BE, Lim S. Simulating an elastic ring with bend and twist by an adaptive generalized immersed boundary method. Commun Comput Phys. 2012;12(2):433–461.
38. Griffith BE, Luo X, McQueen DM, Peskin CS. Simulating the fluid dynamics of natural and prosthetic heart valves using the immersed boundary method. Int J Appl Mech. 2009;1(01):137–177.
39. Griffith BE, Luo XY. Hybrid finite difference/finite element version of the immersed boundary method. (submitted)
40. Guy RD, Phillip B, Griffith BE. Geometric multigrid for an implicit-time immersed boundary method. Adv Comput Math. 2015;41(3):635–662.
41. Humphrey JD. Cardiovascular Solid Mechanics: Cells, Tissues, and Organs. Springer Verlag; 2002.
42. Kim Y, Peskin CS. Penalty immersed boundary method for an elastic boundary with mass. Phys Fluid. 2007;19:18 pages. 053,103.
43. Kim Y, Zhu L, Wang X, Peskin CS. On various techniques for computer simulation of boundaries with mass. In: Bathe KJ, editor. Proceedings of the Second MIT Conference on Computational Fluid and Solid Mechanics; Elsevier; 2003. pp. 1746–1750.
44. Kovács SJ, McQueen DM, Peskin CS. Modelling cardiac fluid dynamics and diastolic function. Phil Trans R Soc Lond A. 2001;359(1783):1299–1314.
45. Lansac E, Lim HS, Shomura Y, Lim KH, Rice NT, Goetz W, Acar C, Duran CMG. A four-dimensional study of the aortic root dynamics. Eur J Cardio Thorac Surg. 2002;22(4):497–503. [PubMed]
46. Lu J, Zhou X, Raghavan ML. Inverse elastostatic stress analysis in pre-deformed biological structures: demonstration using abdominal aortic aneurysms. J Biomech. 2007;40(3):693–696. [PubMed]
47. Luo XY, Griffith BE, Ma XS, Yin M, Wang TJ, Liang CL, Watton PN, Bernacca GM. Effect of bending rigidity in a dynamic model of a polyurethane prosthetic mitral valve. Biomech Model Mechanobiology. 2012;11(6):815–827. [PubMed]
48. Ma XS, Gao H, Griffith BE, Berry C, Luo XY. Image-based fluid-structure interaction model of the human mitral valve. Comput Fluid. 2013;71:417–425.
49. Marom G, Haj-Ali R, Raanani E, Schäfers HJ, Rosenfeld M. A fluid–structure interaction model of the aortic valve with coaptation and compliant aortic root. Med Biol Eng Comput. 2012;50(2):173–182. [PubMed]
50. May-Newman K, Lam C, Yin FCP. A hyperelastic constitutive law for aortic valve tissue. J Biomech Eng. 2009;131(8):7 pages. 081,009. [PubMed]
51. McQueen DM, Peskin CS. Shared-memory parallel vector implementation of the immersed boundary method for the computation of blood flow in the beating mammalian heart. J Supercomput. 1997;11(3):213–236.
52. McQueen DM, Peskin CS. A three-dimensional computer model of the human heart for studying cardiac fluid dynamics. Comput Graph. 2000;34(1):56–60.
53. McQueen DM, Peskin CS. Heart simulation by an immersed boundary method with formal second-order accuracy and reduced numerical viscosity. In: Aref H, Phillips JW, editors. Mechanics for a New Millennium, Proceedings of the 20th International Conference on Theoretical and Applied Mechanics (ICTAM 2000); Kluwer Academic Publishers.2001.
54. Mori Y, Peskin CS. Implicit second order immersed boundary methods with boundary mass. Comput Meth Appl Mech Eng. 2008;197(25–28):2049–2067.
55. Murgo JP, Westerhof N, Giolma JP, Altobelli SA. Aortic input impedance in normal man: relationship to pressure wave forms. Circulation. 1980;62(1):105–116. [PubMed]
56. Nichols WW, O’Rourke MF. McDonald’s Blood Flow in Arteries: Theoretical, Experimental, and Clinical Principles. CRC Press; 2011.
57. Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys. 1972;10(2):252–271.
58. Peskin CS. The immersed boundary method. Acta Numer. 2002;11(0):479–517.
59. Peskin CS, McQueen DM. Mechanical equilibrium determines the fractal fiber architecture of aortic heart valve leaflets. Am J Physiol Heart Circ Physiol. 1994;266(1):H319–H328. [PubMed]
60. Peskin CS, McQueen DM. Fluid dynamics of the heart and its valves. In: Othmer HG, Adler FR, Lewis MA, Dallon JC, editors. Case Studies in Mathematical Modeling: Ecology, Physiology, and Cell Biology. Englewood Cliffs, NJ, USA: Prentice-Hall; 1996. pp. 309–337.
61. Reul H, Vahlbruch A, Giersiepen M, Schmitz-Rode TH, Hirtz V, Effert S. The geometry of the aortic root in health, at valve disease and after valve replacement. J Biomech. 1990;23(2):181–191. [PubMed]
62. Roy S, Heltai L, Costanzo F. Benchmarking the immersed finite element method for fluid-structure interaction problems ArXiv preprint arXiv:1306.0936
63. Sacks MS, Yoganathan AP. Heart valve function: a biomechanical perspective. Philosophical Transactions of the Royal Society B: Biological Sciences. 2007;362(1484):1369–1391. [PMC free article] [PubMed]
64. Sauren AAHJ. Ph.D. thesis. Technische Universiteit Eindhoven; 1981. The mechanical behaviour of the aortic valve.
65. Sellier M. An iterative method for the inverse elasto-static problem. J Fluid Struct. 2011;27(8):1461–1470.
66. Shunk KA, Garot J, Atalar E, Lima JA. Transesophageal magnetic resonance imaging of the aortic arch and descending thoracic aorta in patients with aortic atherosclerosis. Journal of the American College of Cardiology. 2001;37(8):2031–2035. [PubMed]
67. Singh IM, Shishehbor MH, Christofferson RD, Tuzcu EM, Kapadia SR. Percutaneous treatment of aortic valve stenosis. Cleve Clin J Med. 2008;75(11):805–812. [PubMed]
68. Sotiropoulos F, Borazjani I. A review of state-of-the-art numerical methods for simulating flow through mechanical heart valves. Med Biol Eng Comput. 2009;47(3):245–256. [PMC free article] [PubMed]
69. Stergiopulos N, Westerhof BE, Westerhof N. Total arterial inertance as the fourth element of the windkessel model. Am J Physiol Heart Circ Physiol. 1999;276(1):H81–H88. [PubMed]
70. Swanson WM, Clark RE. Dimensions and geometric relationships of the human aortic valve as a function of pressure. Circ Res. 1974;35(6):871–882. [PubMed]
71. Thubrikar M. The Aortic Valve. CRC Press; 1989.
72. Vaishnav RN, Vossoughi J. Estimation of residual strain in aortic segment. In: Hall CW, editor. Biomedical Engineering II: Recent Developments. Pergamon Press; 1983.
73. Vavourakis V, Papaharilaou Y, Ekaterinaris JA. Coupled fluid–structure interaction hemodynamics in a zero-pressure state corrected arterial geometry. J Biomech. 2011;44(13):2453–2460. [PubMed]
74. Viscardi F, Vergara C, Antiga L, Merelli S, Veneziani A, Puppini G, Faggian G, Mazzucco A, Luciani GB. Comparative finite element model analysis of ascending aortic flow in bicuspid and tricuspid aortic valve. Artif Organs. 2010;34(12):1114–1120. [PubMed]
75. Weinberg EJ, Kaazempur Mofrad MR. A multiscale computational comparison of the bicuspid and tricuspid aortic valves in relation to calcific aortic stenosis. J Biomech. 2008;41(16):3482–3487. [PubMed]
76. Westerhof N, Stergiopulos N, Noble MIM. Snapshots of Hemodynamics: An Aid for Clinical Research and Graduate Education. Springer; 2010.
77. Wittek A, Karatolios K, Bihari P, Schmitz-Rixen T, Moosdorf R, Vogt S, Blase C. In vivo determination of elastic properties of the human aorta based on 4d ultrasound data. Journal of the mechanical behavior of biomedical materials. 2013;27:167–183. [PubMed]
78. Yao J, Liu GR, Narmoneva DA, Hinton RB, Zhang ZQ. Immersed smoothed finite element method for fluid–structure interaction simulation of aortic valves. Comput Mech. 2012;50(6):789–804.
79. Zhu L, Peskin CS. Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method. J Comput Phys. 2002;179(2):452–468.