Immersed boundary model of aortic heart valve dynamics with physiological driving and loading conditions

TL;DR

This paper presents a three-dimensional immersed boundary model of the aortic heart valve that simulates realistic cardiac cycle dynamics under physiological conditions, using adaptive numerical methods and clinical data-driven loading.

Contribution

It introduces an adaptive, staggered-grid IB method for simulating aortic valve dynamics with realistic physiological driving and loading conditions, including a thick aortic root wall.

Findings

Model approaches steady periodic state rapidly

Flow rates emerge naturally at physiological levels

Simulation captures realistic valve and flow dynamics

Abstract

The immersed boundary (IB) method is a mathematical and numerical framework for problems of fluid-structure interaction, treating the particular case in which an elastic structure is immersed in a viscous incompressible fluid. The IB approach to such problems is to describe the elasticity of the immersed structure in Lagrangian form, and to describe the momentum, viscosity, and incompressibility of the coupled fluid-structure system in Eulerian form. Interaction between Lagrangian and Eulerian variables is mediated by integral equations with Dirac delta function kernels. The IB method provides a unified formulation for fluid-structure interaction models involving both thin elastic boundaries and also thick viscoelastic bodies. In this work, we describe the application of an adaptive, staggered-grid version of the IB method to the three-dimensional simulation of the fluid dynamics of the aortic heart valve. Our model describes the thin leaflets of the aortic valve as immersed elastic boundaries, and describes the wall of the aortic root as a thick, semi-rigid elastic structure. A physiological left-ventricular pressure waveform is used to drive flow through the model valve, and dynamic pressure loading conditions are provided by a reduced (zero-dimensional) circulation model that has been fit to clinical data. We use this model and method to simulate aortic valve dynamics over multiple cardiac cycles. The model is shown to approach rapidly a periodic steady state in which physiological cardiac output is obtained at physiological pressures. These realistic flow rates are not specified in the model, however. Instead, they emerge from the fluid-structure interaction simulation.

Tables

Table 1. Table 1: Average wallclock time (in seconds) per timestep during the first 128 timesteps of the simulation for different numerical and computational parameters.

	number of cores
	24	48	96
uniform grid	13.02	7.62	7.24
2-level grid	9.09	6.21	4.54
3-level grid	9.19	6.24	4.54

Equations

$$\rho \left(\frac{\partial \mathrm{u}}{\partial t}(\mathrm{x},t)+\mathrm{u}(\mathrm{x},t)\cdot \nabla \mathrm{u}(\mathrm{x},t)\right)$$

$$\nabla \cdot \mathrm{u}(\mathrm{x},t)$$

$$\mathrm{f}(\mathrm{x},t)$$

$$\frac{\partial \mathrm{X}}{\partial t}(\mathrm{s},t)$$

$$\mathrm{F}(\mathrm{s},t)$$

$$\mathrm{F}=-\frac{\delta E}{\delta \mathrm{X}},$$

$$\delta E[\mathrm{X}(\cdot ,t)]=-{\int }_U\mathrm{F}(\mathrm{s},t)\cdot \delta \mathrm{X}(\mathrm{s},t)\mathrm{d}\mathrm{s}.$$

$$E_{\text{s}}={\int }_{\Omega }{\mathcal{E}}_{\text{s}}\left(\left|\frac{\partial \mathrm{X}}{\partial s_3}\right|;\mathrm{s}\right)\mathrm{d}\mathrm{s},$$

$${\mathrm{F}}_{\text{s}}=\frac{\partial }{\partial s_3}\left({\mathcal{E}}_{\text{s}}^{\prime }\left(\left|\frac{\partial \mathrm{X}}{\partial s_3}\right|;\mathrm{s}\right)\frac{\partial \mathrm{X}/\partial s_3}{\left|\partial \mathrm{X}/\partial s_3\right|}\right),$$

$${\mathrm{F}}_{\text{s}}=\frac{\partial }{\partial s_3}\left(T\tau \right).$$

$$E_{\text{b}}=\frac{1}{2}{\int }_{\Omega }{c}_{\text{b}}(\mathrm{s}){\left|\frac{{\partial }^2\mathrm{X}}{\partial s_3^2}-\frac{{\partial }^2\bar{\mathrm{X}}}{\partial s_3^2}\right|}^2\mathrm{d}\mathrm{s},$$

$${\mathrm{F}}_{\text{b}}=\frac{{\partial }^2}{\partial s_3^2}\left(c_{\text{b}}(\mathrm{s})\left(\frac{{\partial }^2\bar{\mathrm{X}}}{\partial s_3^2}-\frac{{\partial }^2\mathrm{X}}{\partial s_3^2}\right)\right).$$

$$\sigma =-p\mathbb{I}+\mu \left[\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^T\right].$$

$$\mathrm{n}\cdot \sigma \cdot \mathrm{n}=-p+2\mu \frac{\partial }{\partial n}(\mathrm{u}\cdot \mathrm{n}).$$

$$\mathrm{n}\cdot \sigma \cdot \mathrm{n}=-p.$$

$$p(\mathrm{x},t)=P^{\text{LV}}(t),\mathrm{x}\in \partial {\Omega }^{\text{LV}}.$$

$$Q^{\text{LV}}(t)=-{\int }_{\partial {\Omega }^{\text{LV}}}\mathrm{u}(\mathrm{x},t)\cdot \mathrm{n}\mathrm{d}a,$$

$$p(\mathrm{x},t)=P^{\text{Ao}}(t),\mathrm{x}\in \partial {\Omega }^{\text{Ao}}.$$

$$Q^{\text{Ao}}(t)={\int }_{\partial {\Omega }^{\text{Ao}}}\mathrm{u}(\mathrm{x},t)\cdot \mathrm{n}\mathrm{d}a.$$

$$\frac{1}{R_{\text{c}}}\left(P^{\text{Ao}}(t)-P^{\text{Wk}}(t)\right)$$

$$C\frac{\text{d}{P}^{\text{Wk}}}{\text{d}t}(t)+\frac{1}{R_{\text{p}}}{P}^{\text{Wk}}(t)$$

$${\left({\nabla }_h\cdot \text{mbox}\mathrm{u}\right)}_{i,j,k}=$$

$$\frac{{\left(u_3\right)}_{i,j,k+\frac{1}{2}}-{\left(u_3\right)}_{i,j,k-\frac{1}{2}}}{h^{\ell }}.$$

$$({\nabla }_{h}p{)}_{i-\frac{1}{2},j,k}$$

$$({\nabla }_{h}p{)}_{i,j-\frac{1}{2},k}$$

$$({\nabla }_{h}p{)}_{i,j,k-\frac{1}{2}}$$

$$({\nabla }_h^2{u}_1{)}_{i-\frac{1}{2},j,k}=$$

$$\frac{{\left(u_1\right)}_{i-\frac{1}{2},j+1,k}-2{\left(u_1\right)}_{i-\frac{1}{2},j,k}+{\left(u_1\right)}_{i-\frac{1}{2},j-1,k}}{{\left(h^{\ell }\right)}^2}+\text{mbox}$$

$$\frac{{\left(u_1\right)}_{i-\frac{1}{2},j,k+1}-2{\left(u_1\right)}_{i-\frac{1}{2},j,k}+{\left(u_1\right)}_{i-\frac{1}{2},j,k-1}}{{\left(h^{\ell }\right)}^2}.$$

$$(f_1{)}_{i-\frac{1}{2},j,k}$$

$$(f_2{)}_{i,j-\frac{1}{2},k}$$

$$(f_3{)}_{i,j,k-\frac{1}{2}}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}(X_1{)}_{l,m,n}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}(X_2{)}_{l,m,n}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}(X_3{)}_{l,m,n}$$

$$\rho \left(\frac{{\mathrm{u}}^{n+1,k+1}-{\mathrm{u}}^n}{\Delta t}+{\mathrm{N}}^{n+\frac{1}{2},k}\right)$$

$${\nabla }_h\cdot \text{mbox}{\mathrm{u}}^{n+1,k+1}$$

$${\mathrm{f}}^{n+\frac{1}{2},k}$$

Full text

1 Introduction

The immersed boundary (IB) method is a mathematical and numerical approach to problems of fluid-structure interaction that was introduced by Peskin to model the fluid dynamics of heart valves [1, 2]. The IB methodology has subsequently been used to model diverse problems in biological fluid dynamics [3] and other problems in which a rigid or elastic structure is immersed in a fluid flow [3, 4, 5]. The IB method for fluid-structure interaction treats problems in which an elastic structure is immersed in a viscous incompressible fluid, describing the elasticity of the immersed structure in Lagrangian form, and describing the momentum, viscosity, and incompressibility of the coupled fluid-structure system in Eulerian form. Integral equations with Dirac delta function kernels couple the Lagrangian and Eulerian variables. When discretized for computer simulation, the IB method approximates the Lagrangian equations on a curvilinear mesh, approximates the Eulerian equations on a Cartesian grid, and approximates the Lagrangian-Eulerian interaction equations by replacing the singular delta function with a regularized version of the delta function.

A key strength of the IB approach to fluid-structure interaction is that it does not require conforming Lagrangian and Eulerian discretizations. Specifically, the IB method permits the Lagrangian mesh to cut through the background Eulerian grid in an arbitrary manner and does not require dynamically generated body-fitted meshes. This attribute of the method greatly simplifies the task of grid generation and facilitates simulations involving large deformations of the elastic structure. An additional feature of the IB formulation is that it provides a unified approach to constructing models involving both thin elastic boundaries (i.e., immersed structures that are of codimension one with respect to the fluid) and also thick elastic bodies (i.e., immersed structures that are of codimension zero with respect to the fluid) [3].

In this work, we describe an adaptive version of the IB method and the application of this method to the simulation of the fluid dynamics of the aortic heart valve. Each year, approximately 250,000 procedures are performed to repair or replace damaged or destroyed heart valves [6]. Severe aortic valve disease is generally treated by replacement with either a mechanical or a bioprosthetic valve [7], and approximately 50,000 aortic valve replacements are performed annually to treat severe aortic stenosis [8, 9, 10]. Because many of the difficulties of prosthetic heart valves are related to the fluid dynamics of the replacement valve [6], mathematical and computational models that enable the study of the fluid-mechanical mechanisms of valve function and dysfunction may ultimately aid in improving treatment outcomes for the many patients suffering from valvular heart diseases.

The aortic valve model employed herein is similar, but not identical, to that described by Griffith et al. [11]. We model the thin leaflets of the aortic valve as immersed boundaries comprised of systems of elastic fibers that resist extension, compression, and bending, and we model the aortic root and ascending aorta as a thick, semi-rigid elastic structure. To construct the model valve leaflets, we use the mathematical theory of Peskin and McQueen [12], which describes the architecture of the systems of collagen fibers within the valve leaflets that allow the closed valve to support a significant pressure load. The geometry of the model aortic root is based on the idealized description of Reul et al. [13], which was derived from imaging data collected from healthy patients, and we use dimensions that are based on measurements by Swanson and Clark [14] of human aortic roots harvested after autopsy. A Windkessel model fit to human data by Stergiopulos et al. [15] provides physiological loading conditions for the model valve. Two limitations of our earlier model [11], which are overcome in the present work, are that it used only a highly idealized left-ventricular driving pressure waveform, and that it considered only a single cardiac cycle. In this work, we use a physiological driving pressure waveform that is based on human clinical data [16], and we perform multibeat simulations of the fluid dynamics of the aortic heart valve. We emphasize that we do not prescribe the flow rate at either the upstream or downstream boundaries of the model vessel. Instead, we impose a realistic, periodic left-ventricular driving pressure at the upstream boundary along with a dynamic circulatory loading pressure at the downstream boundary. With the driving and loading conditions used in the present work, our model rapidly approaches a periodic steady state in which physiological cardiac output is obtained at physiological pressures.

There are also important differences between the numerical methods used in the present study and those of our earlier simulations of aortic valve dynamics [11]. Although both use an adaptive version of the IB method for fluid-structure interaction, our earlier study used a cell-centered IB method [17], whereas herein we use a staggered-grid discretization. This is notable because we have recently demonstrated that staggered-grid IB methods yield substantially improved accuracy when compared to cell-centered discretizations [18]. Specifically, we have found that using a staggered-grid Eulerian discretization improves the volume-conservation properties of the IB method by one to two orders of magnitude in comparison to a cell-centered discretization [18]. Staggered-grid IB methods also yield improved resolution of pressure discontinuities [18]. In the present application, such discontinuities occur along the thin heart valve leaflets and are especially pronounced when the valve is closed and supporting a significant, physiological pressure load.

The three-dimensional adaptive IB method used in our simulations is similar to the two-dimensional adaptive IB method of Roma et al. [19]. Specifically, both schemes use a globally second-order accurate staggered-grid (i.e., marker-and-cell or MAC [20]) discretization of the incompressible Navier-Stokes equations on block-structured adaptively refined Cartesian grids, and both schemes implement formally second-order accurate versions of the IB method (i.e., schemes that yield second-order convergence rates for problems with sufficiently smooth solutions [21, 22]). There are also important differences between the present scheme and the scheme of Roma et al. For instance, the method of Roma et al. uses centered differencing to approximate the nonlinear advection terms of the incompressible Navier-Stokes equations, whereas we use a staggered-grid version [23] of the xsPPM7 variant [24] of the piecewise parabolic method (PPM) [25] that enables the application of the present method to high Reynolds number flows. Our scheme also uses the projection method not as a fractional-step solver for the incompressible Navier-Stokes equations, but rather as a preconditioner for an iterative Krylov method applied to an unsplit discretization of those equations [23]. Our approach eliminates the timestep-splitting error associated with standard projection methods. It also greatly simplifies the specification of physical boundary conditions along the outer boundaries of the computational domain. In this work, such physical boundary conditions couple the detailed, three-dimensional fluid-structure interaction model to the reduced circulation models that provide realistic driving and loading conditions.

2 The continuous equations of motion

The IB formulation of the equations of motion for a coupled fluid-structure system describes the elasticity of the immersed structure in Lagrangian form and describes the momentum, velocity, and incompressibility of the fluid-structure system in Eulerian form. Let $\bm{\mathrm{x}}=(x_{1},x_{2},x_{3})\in\Omega$ denote Cartesian physical coordinates, with $\Omega\subset\mathbb{R}^{3}$ denoting the physical region that is occupied by the fluid-structure system; let $\bm{\mathrm{s}}=(s_{1},s_{2},s_{3})\in U$ denote Lagrangian material coordinates that are attached to the immersed elastic structure, with $U\subset\mathbb{R}^{3}$ denoting the Lagrangian coordinate domain; and let $\bm{\mathrm{X}}(\bm{\mathrm{s}},t)\in\Omega$ denote the physical position of material point $\bm{\mathrm{s}}$ at time $t$. We consider the case in which the fluid possesses a uniform mass density $\rho$ and dynamic viscosity $\mu$, and we assume that the structure is neutrally buoyant and has the same viscous properties as the fluid in which it is immersed. These assumptions are not essential to the method, however, and generalizations of the IB method have been developed to permit the mass density of the structure to differ from that of the fluid [3, 26, 27, 28, 29, 30]. Work is also underway to develop new extensions of the IB method that permit the viscosity of the structure to differ from that of the fluid.

The IB formulation of the equations of fluid-structure interaction is [3]:

[TABLE]

in which $\bm{\mathrm{u}}(\bm{\mathrm{x}},t)=(u_{1}(\bm{\mathrm{x}},t),u_{2}(\bm{\mathrm{x}},t),u_{3}(\bm{\mathrm{x}},t))$ is the Eulerian velocity field, $p(\bm{\mathrm{x}},t)$ is the Eulerian pressure, $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)=(f_{1}(\bm{\mathrm{x}},t),f_{2}(\bm{\mathrm{x}},t),f_{3}(\bm{\mathrm{x}},t))$ is the Eulerian elastic force density (i.e., the elastic force density with respect to the physical coordinate system, so that $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)\,{\rm d}\bm{\mathrm{x}}$ has units of force), $\bm{\mathrm{F}}(\bm{\mathrm{s}},t)=(F_{1}(\bm{\mathrm{s}},t),F_{2}(\bm{\mathrm{s}},t),F_{3}(\bm{\mathrm{s}},t))$ is the Lagrangian elastic force density (i.e., the elastic force density with respect to the material coordinate system, so that $\bm{\mathrm{F}}(\bm{\mathrm{s}},t)\,{\rm d}\bm{\mathrm{s}}$ has units of force), $\bm{\mathrm{\mathcal{F}}}:\bm{\mathrm{X}}\mapsto\bm{\mathrm{F}}$ is a functional that specifies the Lagrangian elastic force density in terms of the deformation of the immersed structure, and $\delta(\bm{\mathrm{x}})=\delta(x_{1})\,\delta(x_{2})\,\delta(x_{3})$ is the three-dimensional Dirac delta function. In this formulation, eqs. (3) and (4) are the interaction equations that couple the Lagrangian and Eulerian variables. Eq. (3) converts the Lagrangian force density $\bm{\mathrm{F}}(\bm{\mathrm{s}},t)$ into the equivalent Eulerian force density $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)$. Eq. (4) states that the physical position of each Lagrangian material point $\bm{\mathrm{s}}$ moves with velocity $\bm{\mathrm{u}}(\bm{\mathrm{X}}(\bm{\mathrm{s}},t),t)$, thereby implying that there is no fluid slip at fluid-structure interfaces. Notice, however, that the no-slip condition of a viscous fluid does not appear in the equations as a constraint on the fluid motion. Instead, the no-slip condition determines the motion of the immersed structure. See, e.g., Peskin [3] for further discussion of these equations.

We next describe the form of the Lagrangian elastic force density functional $\bm{\mathrm{\mathcal{F}}}:\bm{\mathrm{X}}\mapsto\bm{\mathrm{F}}$ used in our model. Like our earlier study of cardiac valve dynamics [11], we model the flexible leaflets of the aortic valve as thin elastic boundaries, and we model the vessel wall as a thick, semi-rigid elastic structure. The elasticity of these structures is described in terms of families of fibers that resist extension, compression, and bending. We identify the model fibers of the valve leaflets with the collagen fibers that enable the real valve to support a significant pressure load when closed. In the case of the vessel wall, we do not identify the model fibers with particular physiological features; instead, these fibers are used simply to fix the geometry of the vessel.

As is frequently done in IB models [3], we define the fiber elasticity in terms of a strain-energy functional $E=E[\bm{\mathrm{X}}(\cdot,t)]$. The corresponding Lagrangian elastic force density may be expressed in terms of the Fréchet derivative of $E$. Specifically, $\bm{\mathrm{F}}$ is defined by

[TABLE]

which is shorthand for

[TABLE]

Notice that in eqs. (6) and (7), $\delta$ denotes the perturbation operator, not the Dirac delta function. To specify $E$, it is convenient to choose the Lagrangian coordinates $\bm{\mathrm{s}}=(s_{1},s_{2},s_{3})\in U$ so that each fixed value of $(s_{1},s_{2})$ labels a particular fiber. This implies that the mapping $s_{3}\mapsto\bm{\mathrm{X}}(s_{1}^{0},s_{2}^{0},s_{3})$ is a parametric representation of the fiber labeled by $(s_{1},s_{2})=(s_{1}^{0},s_{2}^{0})$. The curvilinear coordinate $s_{3}$ need not correspond to arc length along the fiber, however, and even if $s_{3}$ were to correspond to arc length in an initial or reference configuration, notice that it generally will not remain arc length as the structure deforms.

As we have done previously [11], we describe the total elastic energy functional $E$ as the sum of a stretching energy $E_{\text{s}}$ and a bending energy $E_{\text{b}}$, so that $E=E_{\text{s}}+E_{\text{b}}$. In turn, these elastic energy functionals determine a stretching force density $\bm{\mathrm{F}}_{\text{s}}$ and a bending force density $\bm{\mathrm{F}}_{\text{b}}$, so that $\bm{\mathrm{F}}=\bm{\mathrm{F}}_{\text{s}}+\bm{\mathrm{F}}_{\text{b}}$. The stretching energy is

[TABLE]

in which ${\mathcal{E}}_{\text{s}}$ is a local stretching energy. The corresponding stretching force is given by [3]

[TABLE]

in which ${\mathcal{E}}_{\text{s}}^{\prime}$ indicates the derivative of ${\mathcal{E}}_{\text{s}}$ with respect to its first argument. By identifying $T={\mathcal{E}}_{\text{s}}^{\prime}\left(\left|\partial\bm{\mathrm{X}}/\partial s_{3}\right|;\bm{\mathrm{s}}\right)$ as the fiber tension and $\bm{\mathrm{\tau}}=\partial\bm{\mathrm{X}}/\partial s_{3}/\left|\partial\bm{\mathrm{X}}/\partial s_{3}\right|$ as the fiber-aligned unit tangent vector, we may rewrite $\bm{\mathrm{F}}_{\text{s}}$ as

[TABLE]

The bending energy used in our model is

[TABLE]

in which $c_{\text{b}}=c_{\text{b}}(\bm{\mathrm{s}})$ is the spatially inhomogeneous bending stiffness, and $\bar{\bm{\mathrm{X}}}=\bar{\bm{\mathrm{X}}}(\bm{\mathrm{s}})$ is the reference configuration of the structure. The corresponding bending-resistant force is given by [11]

[TABLE]

We take the reference configuration to be the initial configuration, i.e., $\bar{\bm{\mathrm{X}}}(\bm{\mathrm{s}})=\bm{\mathrm{X}}(\bm{\mathrm{s}},0)$. We remark that we use bending-resistant forces only within the model valve leaflets, i.e., $c_{\text{b}}(\bm{\mathrm{s}})\neq 0$ only for those fibers that comprise the valve leaflets. Because we model the valve leaflets as thin elastic surfaces immersed in fluid, including bending-resistant forces allows the model leaflets to account for the thickness of real valve leaflets, which are thin but, of course, not infinitely thin. In our model, we increase $c_{\text{b}}$ near the tips of the free edges of the valve leaflets to account for the fibrous noduli arantii.

Next, we specify the boundary conditions imposed along the outer boundary of the physical domain $\Omega$. We take $\Omega$ to be a $\mbox{10 cm}\times\mbox{10 cm}\times\mbox{15 cm}$ rectangular box, and we employ a combination of solid-wall and prescribed-pressure boundary conditions along $\partial\Omega$. Solid-wall boundary conditions are simply homogeneous Dirichlet conditions for the velocity field $\bm{\mathrm{u}}(\bm{\mathrm{x}},t)$. At solid-wall boundaries, a boundary condition for the pressure is neither needed nor permitted. By prescribed-pressure boundary conditions, we mean a combination of normal-traction and zero-tangential-slip boundary conditions. For a viscous incompressible fluid, it is easy to show that combining normal-traction and zero-tangential-slip boundary conditions along a flat boundary allows for the pointwise specification of the pressure $p(\bm{\mathrm{x}},t)$ on that boundary. To see this, recall that the Cauchy stress tensor of a viscous incompressible fluid is

[TABLE]

Let the outward unit normal at a position $\bm{\mathrm{x}}\in\partial\Omega$ be denoted by $\bm{\mathrm{n}}=\bm{\mathrm{n}}(\bm{\mathrm{x}})$, and let a unit tangent vector at a position $\bm{\mathrm{x}}\in\partial\Omega$ be denoted by $\bm{\mathrm{t}}=\bm{\mathrm{t}}(\bm{\mathrm{x}})$. By prescribing the normal traction at the boundary, we are prescribing the value of the normal component of the normal stress, i.e.,

[TABLE]

The zero-tangential-slip condition imposed on $\bm{\mathrm{u}}$ implies that $\bm{\mathrm{u}}\cdot\bm{\mathrm{t}}\equiv 0$ along $\partial\Omega$, and combining this condition with the incompressibility constraint implies that $\frac{\partial}{\partial n}(\bm{\mathrm{u}}\cdot\bm{\mathrm{n}})\equiv 0$ along $\partial\Omega$. Therefore, along $\partial\Omega$, the normal component of the normal stress reduces to

[TABLE]

Thus, combining normal-traction and zero-tangential-slip boundary conditions allows us to prescribe the value of the pressure pointwise along the boundary.

The model vessel attaches directly to $\partial\Omega$, the outermost boundary of the physical domain. A schematic diagram is provided in fig. 1. Along $\partial\Omega^{\text{LV}}$, the upstream boundary of the vessel, a time-dependent left-ventricular pressure waveform $P^{\text{LV}}(t)$ is prescribed to drive flow through the model valve, so that

[TABLE]

The specific left-ventricular pressure waveform used in the present model is adapted from the study of Murgo et al. [16]. The rate of flow entering the model vessel via the upstream boundary is

[TABLE]

in which $\bm{\mathrm{n}}$ is the outward unit normal along $\partial\Omega$, and ${\rm d}a={\rm d}a(\bm{\mathrm{x}})$ is the area element in the Cartesian coordinate system.

On $\partial\Omega^{\text{Ao}}$, the downstream boundary of the vessel, we use a reduced (i.e., ordinary differential equation) circulation model to determine the pressure $P^{\text{Ao}}(t)$ that provides dynamic loading conditions for the model valve, so that

[TABLE]

The reduced circulation model used in this work is a three-element Windkessel model with characteristic resistance $R_{\text{c}}$, peripheral resistance $R_{\text{p}}$, and arterial compliance $C$; see fig. 1. The rate of flow leaving the model vessel via the downstream boundary is

[TABLE]

Because we specify the pressure along $\partial\Omega^{\text{Ao}}$, the value of $Q^{\text{Ao}}(t)$ is not known in advance; instead, it must be determined by the coupled model. The flow leaving the fluid-structure interaction model via $\partial\Omega^{\text{Ao}}$ is exactly the flow through the circulation model, so that

[TABLE]

in which $P^{\text{Wk}}(t)$ is the stored pressure in the Windkessel model. Notice that the value of $P^{\text{Ao}}(t)$ is completely determined by $P^{\text{Wk}}(t)$ and $Q^{\text{Ao}}(t)$. In our simulations, we set $R_{\text{c}}=0.033$ (mmHg ml*-1* s), $R_{\text{p}}=0.79$ (mmHg ml*-1* s), and $C=1.75$ (ml mmHg*-1*), corresponding to the human “Type A” beat characterized by Stergiopulos et al. [15].

Although we have found that coupling the detailed and reduced models via prescribed-pressure boundary conditions works well in practice, other choices of boundary conditions are possible. For instance, it is straightforward to devise boundary conditions for the incompressible Navier-Stokes equations that prescribe the mean pressure along a portion of a flat boundary; see ch. 3 sec. 8 of Gresho and Sani [31] for details. Alternatively, one may wish to couple the detailed and reduced models by prescribing boundary conditions for the normal component of the velocity. A drawback of this approach is that it would require determining an appropriate velocity profile at the boundary. A more serious limitation of this alternative approach is that prescribing the flow rate as a boundary condition, at either the upstream or the downstream boundary, makes it impossible to impose a realistic pressure difference across the model valve during the diastolic phase of the cardiac cycle. By using pressure boundary conditions at both the upstream and downstream boundaries of the vessel, we allow the normal component of the velocity profile at the boundary to be determined by the model, and we are able to impose realistic pressure loads on the model valve throughout the cardiac cycle.

Along $\partial\Omega^{\text{O}}$, the outermost portion of $\partial\Omega$ exterior to the model vessel, we set the pressure to equal zero. This external boundary condition provides an open boundary that acts to couple the fluid-structure interaction model to a zero-pressure fluid reservoir. This constant-pressure reservoir allows the vessel to change volume during the course of the simulation, i.e., it permits a mismatch between the instantaneous flow rates at the inflow and outflow boundaries of the vessel. Because we model the vessel wall as a semi-rigid elastic structure, we have that $Q^{\text{LV}}(t)\approx Q^{\text{Ao}}(t)$. Of course, once the model reaches periodic steady state, the time-integrated inflow and outflow volumes must match. The real aortic root is a flexible structure with significant compliance, however, and a model vessel that accounts for the physiological compliance of the aortic root would generally result in instantaneous differences between the inflow rate through $\partial\Omega^{\text{LV}}$ and the outflow rate through $\partial\Omega^{\text{Ao}}$.

All that remains is to specify the initial conditions. At time $t=0$, we set $\bm{\mathrm{u}}(\bm{\mathrm{x}},t)=0$ along with $P^{\text{LV}}(t)=P^{\text{Wk}}(t)=0$, so that all prescribed normal traction along the boundary are equal to zero. During a brief initialization period lasting 12.8 ms, we increase the left-ventricular driving pressure to a value of approximately 10 mmHg, and we increase the stored pressure in the Windkessel model to 85 mmHg, thereby establishing a realistic pressure load on the closed valve. During this initialization period, $P^{\text{Wk}}(t)$ is treated as a boundary condition and not as a state variable. That is to say, $P^{\text{Wk}}(t)$ does not satisfy eq. (21) for $t\leq\mbox{12.8 ms}$; rather, the value of $P^{\text{Wk}}(t)$ is prescribed. Once the model is initialized, however, $P^{\text{Wk}}(t)$ is treated as a state variable, the dynamics of which are determined by eq. (21).

3 The discrete equations of motion

3.1 Lagrangian and Eulerian spatial discretizations

As in our earlier simulation studies of cardiac fluid dynamics [11, 17], we discretize the Lagrangian equations on a fiber-aligned curvilinear mesh, and we discretize the Eulerian equations on a block-structured locally refined Cartesian grid that is adaptively generated to conform to the moving fiber mesh. The curvilinear mesh spacings are $\Delta s_{1}$, $\Delta s_{2}$, and $\Delta s_{3}$, and we use the indices $(l,m,n)$ to label the nodes of the Lagrangian mesh, so that $\bm{\mathrm{X}}_{l,m,n}$ and $\bm{\mathrm{F}}_{l,m,n}$ are the position and Lagrangian elastic force density associated with curvilinear mesh node $(l,m,n)$. The nodal values of $\bm{\mathrm{F}}_{l,m,n}$ are computed from the physical positions of the nodes of the curvilinear mesh via standard second-order accurate finite difference approximations to $\frac{\partial}{\partial s_{3}}$ and to $\frac{\partial^{2}}{\partial s_{3}^{2}}$. This approach is equivalent to describing the elasticity of the discretized model in terms of systems of springs and beams.

The locally refined Cartesian grid is organized as a hierarchy of nested grid levels that are labeled $\ell=0,\ldots,\ell_{\text{max}}$, with $\ell=0$ denoting the coarsest level of the hierarchical grid and with $\ell=\ell_{\text{max}}$ denoting the finest level. Each grid level $\ell$ is comprised of the union of rectangular Cartesian grid patches. All grid patches in a given level $\ell$ of the grid hierarchy share the same uniform grid spacing $h^{\ell}$, and the grid spacings are chosen so that $h^{\ell}=h^{\ell-1}/{\mathfrak{n}}$, in which ${\mathfrak{n}}>\mbox{1}$ is an integer refinement ratio. The patch levels are constructed to satisfy the proper nesting condition [32], which generally requires that the union of the level $\ell+1$ grid patches be strictly contained within the union of the level $\ell$ grid patches. The proper nesting condition is relaxed at the outermost boundary of the physical domain, thereby allowing high Eulerian spatial resolution all the way up to $\partial\Omega$ in cases in which such resolution is needed. The patch levels are generated so that the faces of the grid patches that comprise level $\ell>0$ are coincident with the faces of the Cartesian grid cells that comprise level $\ell-1$, the next coarser level of the grid. This construction simplifies the development of composite-grid discretization methods that couple the levels of the locally refined grid. Except where noted, in our simulations, we employ a two-level adaptive grid, so that $\ell_{\text{max}}=1$, and we use ${\mathfrak{n}}=4$. An example two-dimensional grid is shown in fig. 2.

To discretize the Eulerian incompressible Navier-Stokes equations in space, we employ a locally refined version of a three-dimensional staggered-grid finite difference scheme; see fig. 3. The computational domain $\Omega$ is a rectangular box, $\Omega=[0,L_{1}]\times[0,L_{2}]\times[0,L_{3}]$, and the coarsest level of the locally refined Cartesian grid is a uniform discretization of $\Omega$, so that the union of the level $\ell=0$ grid patches form a regular $N_{1}\times N_{2}\times N_{3}$ Cartesian grid with grid spacings $\Delta x_{1}=\frac{L_{1}}{N_{1}}$, $\Delta x_{2}=\frac{L_{2}}{N_{2}}$, and $\Delta x_{3}=\frac{L_{3}}{N_{3}}$. For simplicity, we assume that $\Delta x_{1}=\Delta x_{2}=\Delta x_{3}=h^{0}$. On each level $\ell$ of the locally refined grid, $(i,j,k)$ labels a particular Cartesian grid cell, and \bm{\mathrm{x}}_{i,j,k}=\mbox{\left(\left(i+\frac{1}{2}\right)h^{\ell},\left(j+\frac{1}{2}\right)h^{\ell},\left(k+\frac{1}{2}\right)h^{\ell}\right)} denotes the physical location of the center of that cell. The physical region covered by Cartesian grid cell $(i,j,k)$ on level $\ell$ is denoted by $\bm{\mathrm{c}}_{i,j,k}^{\ell}$, and the set of Cartesian grid cell indices associated with level $\ell$ is denoted by $\mathcal{G}^{\ell}$. The components of the Eulerian velocity field $\bm{\mathrm{u}}=(u_{1},u_{2},u_{3})$ are respectively approximated at the centers of the $x_{1}$, $x_{2}$, and $x_{3}$ faces of the Cartesian grid cells, i.e., at positions \bm{\mathrm{x}}_{i-\frac{1}{2},j,k}=\mbox{\left(ih^{\ell},\left(j+\frac{1}{2}\right)h^{\ell},\left(k+\frac{1}{2}\right)h^{\ell}\right)}, \bm{\mathrm{x}}_{i,j-\frac{1}{2},k}=\mbox{\left(\left(i+\frac{1}{2}\right)h^{\ell},jh^{\ell},\left(k+\frac{1}{2}\right)h^{\ell}\right)}, and \bm{\mathrm{x}}_{i,j,k-\frac{1}{2}}=\mbox{\left(\left(i+\frac{1}{2}\right)h^{\ell},\left(j+\frac{1}{2}\right)h^{\ell},kh^{\ell}\right)}. The pressure $p$ is approximated at the centers of the Cartesian grid cells. We use $(u_{1})_{i-\frac{1}{2},j,k}$, $(u_{2})_{i,j-\frac{1}{2},k}$, $(u_{3})_{i,j,k-\frac{1}{2}}$, and $p_{i,j,k}$ to denote the values of $\bm{\mathrm{u}}$ and $p$ that are stored on the grid. A staggered discretization is also used for the Eulerian body force $\bm{\mathrm{f}}=(f_{1},f_{2},f_{3})$, so that $f_{1}$, $f_{2}$, and $f_{3}$ are respectively approximated at the centers of the $x_{1}$, $x_{2}$, and $x_{3}$ faces of the Cartesian grid cells.

Let $\Omega^{\ell}\subseteq\Omega$ denote the physical region covered by the union of the level $\ell$ grid patches. By construction, $\Omega^{0}=\Omega$, and $\Omega^{\ell+1}\subseteq\Omega^{\ell}$. Moreover, away from $\partial\Omega$, $\Omega^{\ell+1}$ is strictly contained within $\Omega^{\ell}$. Notice that $\Omega^{\ell}=\cup_{(i,j,k)\in\mathcal{G}^{\ell}}\bm{\mathrm{c}}_{i,j,k}^{\ell}$. The coarse-fine interface between levels $\ell$ and $\ell+1$ is $\partial\Omega^{\ell+1}\setminus\partial\Omega^{\ell}$. Because the grid levels are constructed to satisfy the proper nesting condition, however, $\partial\Omega^{\ell+1}\cap\partial\Omega^{\ell}\subset\partial\Omega$, so that the coarse-fine interface between levels $\ell$ and $\ell+1$ is $\partial\Omega^{\ell+1}\setminus\partial\Omega$.

The refined region of level $\ell<\ell_{\text{max}}$, denoted by $\Omega^{\ell,\text{ref}}$, consists of the portion of $\Omega^{\ell}$ that is covered by $\Omega^{\ell+1}$, i.e., $\Omega^{\ell,\text{ref}}=\Omega^{\ell}\cap\Omega^{\ell+1}=\Omega^{\ell+1}$. The grid values of any quantity stored on the Cartesian grid that are physically located in $\Omega^{\ell,\text{ref}}$ are referred to as invalid values. The remaining values, i.e., those that are located in $\Omega^{\ell}\setminus\Omega^{\ell,\text{ref}}$, are referred to as valid values. The invalid values of level $\ell<\ell_{\text{max}}$ are constrained to be the restriction of the overlying level $\ell+1$ values. Notice that these overlying values could be either valid or invalid values, depending on the configuration of the grid hierarchy. The simplest restriction procedure is to define the coarse-grid invalid values to be the averages of the overlying fine-grid values; however, other restriction procedures are possible and, in fact, necessary to obtain higher-order accuracy at coarse-fine interfaces.

Each Cartesian grid patch is augmented by a layer of ghost cells that provide the values at the patch boundary that are needed to evaluate the approximations to the Eulerian spatial differential operators in the patch interior. We consider three types of ghost cells for the level $\ell$ grid patches:

(1) ghost cells $(i,j,k)$ that overlap the interior of another level $\ell$ grid patch, so that $\bm{\mathrm{c}}_{i,j,k}^{\ell}\subset\Omega^{\ell}$;

(2) ghost cells $(i,j,k)$ that overlap the interior of some level $\ell-1$ grid patch but not the interior of any of the level $\ell$ grid patches, so that $\bm{\mathrm{c}}_{i,j,k}^{\ell}\subset\Omega^{\ell-1}$ but $\bm{\mathrm{c}}_{i,j,k}^{\ell}\not\subset\Omega^{\ell}$; and

(3) ghost cells that are exterior to $\Omega$, so that $\bm{\mathrm{c}}_{i,j,k}^{\ell}\not\subset\Omega$.

For level $\ell$ ghost cells that overlap the interior of a neighboring level $\ell$-grid patch, the values in that ghost cell are simply copies of the neighboring interior values. If a level $\ell$ ghost cell overlaps the interior of a neighboring level $\ell-1$ grid patch but does not overlap the interior of any level $\ell$ grid patch, that ghost cell is said to be located on the coarse side of a coarse-fine interface. Ghost values located on the coarse side of a coarse-fine interface are defined by an interpolation procedure that uses both fine- and coarse-grid values. As we have done previously [11], we use a specialized quadratic interpolation procedure [33, 34, 35] to compute cell-centered ghost values at coarse-fine interfaces. For face-centered quantities, we use a generalization of this procedure that employs a combination of quadratic and cubic interpolation at coarse-fine interfaces. (Details of these coarse-fine interface discretizations are provided in the appendix.) Finally, ghost values that are exterior to the physical domain are determined via the physical boundary conditions in a manner that is analogous to the uniform-grid scheme for the incompressible Navier-Stokes equations of Griffith [23].

We denote by $\nabla_{h}\cdot\mbox{}\bm{\mathrm{u}}\approx\nabla\cdot\bm{\mathrm{u}}$, $\nabla_{h}p\approx\nabla p$, and $\nabla_{h}^{2}\bm{\mathrm{u}}\approx\nabla^{2}\bm{\mathrm{u}}$ composite-grid finite difference approximations to the divergence, gradient, and Laplace operators, respectively. We employ a standard staggered-grid (MAC) discretization approach, in which $\nabla_{h}\cdot\mbox{}\bm{\mathrm{u}}$ is computed at cell centers, whereas both $\nabla_{h}p$ and $\nabla_{h}^{2}\bm{\mathrm{u}}$ are computed on cell faces. Briefly stated, we use standard second-order accurate staggered-grid approximations to these operators within each Cartesian grid patch [23]. These uniform-grid patch-based discretizations are coupled to form composite-grid discretizations via restriction and prolongation operators. Although our composite-grid finite difference discretization of the incompressible Navier-Stokes equations is globally second-order accurate, this discretization does not retain the pointwise second-order accuracy of the basic uniform-grid approximation because of localized reductions in accuracy at coarse-fine interfaces.

To compute $\nabla_{h}\cdot\mbox{}\bm{\mathrm{u}}$ on the composite grid, we first use simple averaging to restrict $\bm{\mathrm{u}}$ from finer levels of the grid to coarser levels of the grid. We then compute for each level $\ell$ grid cell

[TABLE]

To compute $\nabla_{h}p$ on the composite grid, we first use cubic interpolation to restrict $p$ from finer levels of the grid to coarser levels of the grid, and we then compute the composite-grid cell-centered interpolation of $p$ in the ghost cells along the coarse side of the coarse-fine interface. This approach ensures that the ghost values are at least third-order accurate interpolations of $p$. With the coarse-fine interface ghost values so determined, we then compute for each level $\ell$ cell face

[TABLE]

Finally, to compute $\nabla_{h}^{2}\bm{\mathrm{u}}$, we first use cubic interpolation to restrict $\bm{\mathrm{u}}$ from finer levels of the grid to coarser levels of the grid, and we then compute coarse-fine interface ghost values of $\bm{\mathrm{u}}$ using a composite-grid face-centered interpolation scheme, again ensuring that the ghost values are at least third-order accurate. We then compute for each $x_{1}$ face on level $\ell$

[TABLE]

Similar formulae are used to evaluate $\nabla_{h}^{2}u_{2}$ and $\nabla_{h}^{2}u_{3}$ on the $x_{2}$ and $x_{3}$ faces of the composite grid.

In the case of a uniform-grid discretization, notice that $\nabla_{h}\cdot\mbox{}\nabla_{h}$ is the usual 7-point cell-centered finite difference approximation to the Laplacian. This property facilities the construction of efficient preconditioners based on the projection method [23] and approximate Schur complement methods [36, 37]. In the presence of local mesh refinement, $\nabla_{h}\cdot\mbox{}\nabla_{h}$ is a relatively standard composite-grid generalization of the 7-point Laplacian [11, 17, 33, 34, 35] for which efficient solution algorithms, such as FAC (the fast adaptive composite-grid method) [38, 39, 40], have been developed.

3.2 Lagrangian-Eulerian interaction

To approximate the Lagrangian-Eulerian interaction equations, eqs. (3) and (4), we replace the delta function $\delta(\bm{\mathrm{x}})$ with a regularized version of the delta function $\delta_{h}(\bm{\mathrm{x}})$. The regularized delta function that we use is of the tensor-product form $\delta_{h}(\bm{\mathrm{x}})=\delta_{h}(x_{1})\,\delta_{h}(x_{2})\,\delta_{h}(x_{3})$, and we use a one-dimensional regularized delta function that is of the form $\delta_{h}(x)=\frac{1}{h}\phi\left(\frac{x}{h}\right)$. We take $\phi=\phi(r)$ to be the four-point delta function of Peskin [3].

We have found that fluid-structure interfaces generally require the highest available spatial resolution in an IB simulation. Therefore, we construct the Cartesian grid hierarchy so that the curvilinear mesh is embedded in the finest level of the grid. Additionally, because we use a regularized delta function with a support of four Cartesian meshwidths in each coordinate direction, we ensure that the locally refined grid is generated so that each curvilinear mesh point is physically located at least two level $\ell_{\text{max}}$ grid cells away from $\partial\Omega^{\ell_{\text{max}}}\setminus\partial\Omega$, the coarse-fine interface between levels $\ell_{\text{max}}-1$ and $\ell_{\text{max}}$. With these constraints on the configuration of the locally refined Cartesian grid, the nodes of the curvilinear mesh are directly coupled only to valid values defined on the finest level of the locally refined grid. This construction therefore allows us to discretize the Lagrangian-Eulerian interaction equations as if we were using a uniformly fine grid with resolution $h^{\ell_{\text{max}}}$.

With $\bm{\mathrm{f}}=(f_{1},f_{2},f_{3})$ and $\bm{\mathrm{F}}=(F_{1},F_{2},F_{3})$, we approximate eq. (3) componentwise by

[TABLE]

for $(i,j,k)\in\mathcal{G}^{\ell_{\text{max}}}$. The valid values of $\bm{\mathrm{f}}$ are set to equal zero on all coarser levels of the grid hierarchy. Similarly, with $\bm{\mathrm{u}}=(u_{1},u_{2},u_{3})$ and with $\bm{\mathrm{X}}=(X_{1},X_{2},X_{3})$, we approximate eq. (4) by

[TABLE]

in which we again consider only Cartesian grid cells on the finest level of the hierarchical grid. For those curvilinear mesh nodes that are in the vicinity of physical boundaries, we use the modified regularized delta function formulation of Griffith et al. [11]. This approach ensures that force and torque are conserved during Lagrangian-Eulerian interaction, even near $\partial\Omega$. To simplify the description of our timestepping algorithm, we use the shorthand $\bm{\mathrm{f}}=\bm{\mathrm{\mathcal{S}}}[\bm{\mathrm{X}}]\,\bm{\mathrm{F}}$ and $\frac{{\rm d}}{{\rm d}t}\bm{\mathrm{X}}=\bm{\mathrm{\mathcal{S}}}^{*}[\bm{\mathrm{X}}]\,\bm{\mathrm{u}}$, in which the force-spreading and velocity-interpolation operators, $\bm{\mathrm{\mathcal{S}}}[\bm{\mathrm{X}}]$ and $\bm{\mathrm{\mathcal{S}}}^{*}[\bm{\mathrm{X}}]$, are implicitly defined by eqs. (27)–(29) and (30)–(32), respectively.

3.3 Temporal discretization

We employ a simple timestep-splitting scheme to discretize the equations in time. Briefly, during each timestep, we first solve the fluid-structure interaction equations, treating the values of the upstream and downstream pressure boundary conditions as fixed. We then update the state variables of the circulation model, treating the fluid-structure interaction model state variables as fixed. This approach amounts to a first-order timestep splitting of the equations.

We discretize the fluid-structure interaction equations in time using truncated fixed-point iteration. We treat the linear terms in the incompressible Navier-Stokes equations implicitly, and we treat all other terms explicitly. Let $\bm{\mathrm{X}}^{n+1,k}$, $\bm{\mathrm{u}}^{n+1,k}$, and $p^{n+\frac{1}{2},k}$ denote the approximations to the values of $\bm{\mathrm{X}}$ and $\bm{\mathrm{u}}$ at time $t^{n+1}=(n+1)\Delta t$ and to the value of $p$ at time $t^{n+\frac{1}{2}}=(n+\frac{1}{2})\Delta t$ obtained after $k$ steps of fixed-point iteration, with $\bm{\mathrm{X}}^{n+1,0}=\bm{\mathrm{X}}^{n}$, $\bm{\mathrm{u}}^{n+1,0}=\bm{\mathrm{u}}^{n}$, and $p^{n+\frac{1}{2},0}=p^{n-\frac{1}{2}}$. Letting $\bm{\mathrm{X}}^{n+\frac{1}{2},k}=\frac{1}{2}\left(\bm{\mathrm{X}}^{n}+\bm{\mathrm{X}}^{n+1,k}\right)$ and $\bm{\mathrm{u}}^{n+\frac{1}{2},k}=\frac{1}{2}\left(\bm{\mathrm{u}}^{n}+\bm{\mathrm{u}}^{n+1,k}\right)$, we obtain $\bm{\mathrm{X}}^{n+1,k+1}$, $\bm{\mathrm{u}}^{n+1,k+1}$, and $p^{n+\frac{1}{2},k+1}$ by solving the linear system of equations

[TABLE]

in which $\bm{\mathrm{N}}^{n+\frac{1}{2},k}\approx\left[\bm{\mathrm{u}}\cdot\nabla\bm{\mathrm{u}}\right]^{n+\frac{1}{2}}$ is an explicit approximation to the advection term that uses the xsPPM7 variant [24] of the piecewise parabolic method (PPM) [25] to discretize the nonlinear advection terms; see Griffith [23] for details. We use two cycles of fixed-point iteration per timestep to obtain a second-order accurate timestepping scheme. When solving for $\bm{\mathrm{X}}^{n+1,k+1}$, $\bm{\mathrm{u}}^{n+1,k+1}$, and $p^{n+\frac{1}{2},k+1}$, we fix the pressure at the upstream boundary $\partial\Omega^{\text{LV}}$ to be $p=P^{\text{LV}}(t^{n})$, and we fix the pressure at the downstream boundary $\partial\Omega^{\text{Ao}}$ to be $p=P^{\text{Ao}}(t^{n})$.

Next, having computed $\bm{\mathrm{X}}^{n+1}$, $\bm{\mathrm{u}}^{n+1}$, and $p^{n+\frac{1}{2}}$, we use $\bm{\mathrm{u}}^{n+1}$ to compute $Q^{\text{Ao}}\left(t^{n+1}\right)$, the instantaneous rate of flow leaving the vessel through $\partial\Omega^{\text{Ao}}$ at time $t^{n+1}$. We then update the value of $P^{\text{Wk}}$ via a second-order accurate explicit Runge-Kutta method, so that

[TABLE]

We then set

[TABLE]

which serves as the downstream pressure boundary condition for the fluid-structure interaction model during the subsequent timestep. For the initial timestep, we set $P^{\text{LV}}(0)=P^{\text{Ao}}(0)=0$, values that are consistent with the initial conditions of the continuous system.

3.4 Cartesian grid adaptive mesh refinement

The locally refined grid is constructed in a recursive fashion: First, level 0 is constructed to cover the entire physical domain $\Omega$. Next, having constructed levels $0,\ldots,\ell<\ell_{\text{max}}$, level $\ell+1$ is generated by

(1) tagging cells on level $\ell$ for refinement,

(2) covering the tagged level $\ell$ grid cells by rectangular boxes generated by the Berger-Rigoutsos point-clustering algorithm [41], and

(3) refining the generated boxes by the integer refinement ratio ${\mathfrak{n}}$ to form the level $\ell+1$ grid patches.

Our cell-tagging criteria are simple rules that ensure that the immersed structure remains covered throughout the simulation by the grid cells that comprise the finest level of the hierarchical grid, and that attempt to ensure that flow features requiring enhanced resolution, such as vortices shed from the free edges of the valve leaflets, remain covered by grid cells of an appropriate resolution. Specifically, we tag grid cell $(i,j,k)$ on level $\ell<\ell_{\text{max}}$ for refinement whenever there exists some curvilinear mesh node $(l,m,n)$ such that $\bm{\mathrm{X}}_{l,m,n}\in\bm{\mathrm{c}}_{i,j,k}^{\ell}$, or whenever the local magnitude of the vorticity $\|\bm{\mathrm{\omega}}\|_{i,j,k}=\|\nabla_{h}\times\bm{\mathrm{u}}\|_{i,j,k}$ exceeds a relative threshold. Additional cells are added to the finest level of the grid to ensure that the coarse-fine interface between levels $\ell_{\text{max}}-1$ and $\ell_{\text{max}}$ is sufficiently far away from each of the curvilinear mesh nodes to prevent complicating the discretization of the Lagrangian-Eulerian interaction equations, as discussed previously. We emphasize that the positions of the curvilinear mesh nodes are not constrained to conform in any way to the locally refined Cartesian grid. Instead, the Cartesian grid patch hierarchy is adaptively updated to conform to the evolving configuration of the immersed elastic structure.

To prevent the immersed structure from “escaping” from the finest level of the grid, it is necessary to regenerate the locally refined grid at an interval that is determined by the CFL number of the flow. In our simulations, we choose the timestep size $\Delta t$ to satisfy a CFL condition of the form

[TABLE]

This condition implies that each curvilinear mesh point moves at most $\frac{1}{5}$ fractional meshwidths per timestep. Therefore, to ensure that the immersed structure remains covered by $\Omega^{\ell_{\text{max}}}$, we must adaptively regenerate the grid hierarchy at least every five timesteps. In our simulations, we actually regenerate the grid every four timesteps. In practice, we could generally postpone regridding because the actual timestep size is generally smaller than that required by eq. (41) as a consequence of an additional stability restriction on $\Delta t$ that is of the form $\Delta t=O\left(\left(h^{\ell_{\text{max}}}\right)^{4}\right)$. This severe stability restriction results from our time-explicit treatment of the bending-resistant elastic force. For a model that includes only extension- and compression-resistant elastic elements, the stability restriction is reduced to $\Delta t=O\left(\left(h^{\ell_{\text{max}}}\right)^{2}\right)$.

Each time that the locally refined Cartesian grid is regenerated, Eulerian quantities must be transferred from the old grid hierarchy to the new one. In newly refined regions of the physical domain, the velocity field is prolonged from coarser levels of the old grid via a specialized conservative interpolation scheme that preserves the discrete divergence and curl of the staggered-grid velocity field [42]. (The basic divergence- and curl-preserving interpolation scheme [42] considers only the case in which ${\mathfrak{n}}=2$; however, this procedure is easily generalized to cases in which ${\mathfrak{n}}$ is a power of two via recursion.) The pressure, which is not a state variable of the system and which is used in the subsequent timestep only as an initial approximation to the updated pressure computed during that timestep (see below), is prolonged by simple linear interpolation. In newly coarsened regions of the domain, the values of the velocity and pressure are set to be the averages of the overlying fine-grid values from the old grid hierarchy.

4 Solution methodology

Solving for $\bm{\mathrm{X}}^{n+1,k+1}$, $\bm{\mathrm{u}}^{n+1,k+1}$, and $p^{n+\frac{1}{2},k+1}$ in eqs. (33)–(37) requires the solution of the linear system of equations,

[TABLE]

We solve eq. (42) via the FGMRES algorithm [43], using $\bm{\mathrm{u}}^{n+1,k}$ and $p^{n+\frac{1}{2},k}$ as initial approximations to $\bm{\mathrm{u}}^{n+1,k+1}$ and $p^{n+\frac{1}{2},k+1}$, and using the projection method as a preconditioner. (Recall that we define $\bm{\mathrm{u}}^{n+1,0}=\bm{\mathrm{u}}^{n}$ and $p^{n+\frac{1}{2},0}=p^{n-\frac{1}{2}}$.) Letting $A$ denote the matrix corresponding to this block system, i.e.,

[TABLE]

we write the corresponding projection method-based preconditioner matrix $B$ as [23]

[TABLE]

Because we use a preconditioned Krylov method to solve for $\bm{\mathrm{u}}^{n+1,k+1}$ and $p^{n+\frac{1}{2},k+1}$, it is not necessary to form explicitly the matrices corresponding to $A$ and $B$. Instead, we need only to be able to compute the application of these operators to arbitrary velocity- and pressure-like quantities. Computing the action of $A$ requires implementations of the finite difference approximations to the divergence, gradient, and Laplace operators described previously. Computing the action of $B$ additionally requires solvers for cell-centered and face-centered Poisson-type problems. It is not necessary, however, to employ exact solvers for these subdomain problems. In fact, these subdomain solvers may be quite approximate. At least in the present application, we have found that the performance of our implementation is optimized by using a single multigrid V-cycle of a cell-centered FAC preconditioner for the pressure subsystem, corresponding to the block $(\nabla_{h}\cdot\mbox{}\nabla_{h})^{-1}$ of the preconditioner $B$, and by applying two iterations of the conjugate gradient method for the velocity subsystem, corresponding to the block $\left(\frac{\rho}{\Delta t}I-\frac{\mu}{2}\nabla_{h}^{2}\right)^{-1}$ of the preconditioner. We are able to avoid using a multigrid method for the velocity subsystem because the Reynolds number of the flow is relatively large and the timestep size is relatively small, so that the linear system $\left(\frac{\rho}{\Delta t}I-\frac{\mu}{2}\nabla_{h}^{2}\right)$ is well conditioned.

We remark that although the projection method can be an extremely effective preconditioner [23], it does not appear to be widely used in this manner in practice. Instead, it is generally used as an approximate solver for the incompressible Navier-Stokes equations [34, 35, 44, 45, 46]. As a solver for the incompressible Navier-Stokes equations, the projection method is a fractional-step scheme that first solves the momentum equation over a time interval $[t^{n},t^{n+1}]$ for an “intermediate” velocity field without imposing the constraint of incompressibility, and then projects that intermediate velocity field onto the space of discretely divergence-free vector fields to obtain an approximation to the incompressible velocity field at time $t^{n+1}$. These two steps correspond to the two subdomain solves of our projection preconditioner, and each step requires the imposition of “artificial” physical boundary conditions. When the projection method is used as a solver, the artificial boundary conditions must be chosen carefully to yield a stable and accurate approximation to the true boundary conditions that are to be imposed on the coupled equations [46, 47]. Obtaining high-order accuracy may not be possible in all cases, such as for problems involving outflow boundaries [48], and constructing discretizations that are both stable and accurate can be difficult in practice [49].

We prefer to use the projection method as a preconditioner rather than as a solver. First, doing so permits the Krylov method to eliminate the timestep-splitting error of the basic projection method. Second, solving eq. (42), an unsplit discretization of the incompressible Stokes equations, permits us to impose directly the true boundary conditions on $\bm{\mathrm{u}}$ and $p$ in a coupled manner. Specifically, the form of the artificial boundary conditions required of the basic projection method does not affect the accuracy of the overall solver. This greatly simplifies the development of higher-order accurate discretiations of various types of boundary conditions. See Griffith [23] for details on the construction of accurate discretizations of various combinations of normal and tangential velocity and traction boundary conditions.

5 Implementation

Our adaptive IB method is implemented in the IBAMR software framework, a freely available C++ library for developing fluid-structure interaction models that use the IB method [50]. IBAMR provides support for distributed-memory parallelism via MPI and Cartesian grid adaptive mesh refinement. IBAMR relies upon the SAMRAI [51, 52, 53], PETSc [54, 55, 56], and hypre [57, 58] libraries for much of its functionality.

6 Computational results

6.1 Model results

Using the methods described herein, we have simulated the fluid dynamics of the aortic valve over multiple cardiac cycles, using a time-periodic left-ventricular driving pressure and dynamic loading conditions. In these simulations, the physical domain $\Omega$ is a $\mbox{10 cm}\times\mbox{10 cm}\times\mbox{15 cm}$ rectangular box that we discretize using a two-level adaptively refined Cartesian grid with refinement ratio ${\mathfrak{n}}=4$ between grid levels. The coarse-grid resolution is $h^{0}=\mbox{3.125 mm}$, which corresponds to that of a $32\times 32\times 48$ uniform discretization of $\Omega$, and the fine-grid resolution is $h^{1}=h^{\ell_{\text{max}}}=\mbox{0.78125 mm}$, which corresponds to that of a $128\times 128\times 192$ uniform discretization of $\Omega$. We model the valve leaflets as thin elastic surfaces because real aortic valve leaflets are only 0.2 mm thick [59], which is somewhat thinner than the Cartesian grid spacing on the finest level of the locally refined Cartesian grid. In contrast, we model the vessel as a thick elastic structure because the thickness of the aortic wall is 1 mm [60], which is relatively thick at the scale of the Cartesian grid. In our simulations, we use a uniform timestep size of \Delta t=\mbox{6.25 \mus}, thereby requiring 128,000 timesteps per 0.8 s cardiac cycle. This value of $\Delta t$ was empirically determined to be approximately the largest stable timestep permitted by the present model and semi-explicit numerical method.

As discussed earlier, the elasticity of the valve leaflets and vessel wall is modeled using systems of fibers. In our model, each valve leaflet is spanned by two families of fibers. The first family of fibers runs from commissure to commissure, and the second family runs orthogonal to the commissural fibers. We view the commissural fibers as corresponding to the collagen fibers that allow the real valve leaflets to support a significant pressure load, and we use the mathematical theory of Peskin and McQueen [12] to determine the architecture of these fibers. We remark that the construction of Peskin and McQueen yields discretized fibers that form an orthogonal net, thereby facilitating the construction of the second family of fibers. Because this mathematical theory of the valve fiber architecture describes the geometry of the closed, loaded valve, we assume that the commissural fibers are initialized in a condition of 10% strain, and we choose the commissural fiber stiffness so that the closed valve supports a pressure load of approximately 80 mmHg. The second family of fibers is 10% as stiff as the commissural fibers, thereby approximating the anisotropic material properties of real aortic valve leaflets [61]. Notice that although the valve leaflets are initialized in a strained configuration, the boundary conditions at time $t=0$ do not provide the closed valve with any pressure load. This initial strain thereby induces oscillations in the valve leaflets, as if the valve had been struck like a drumhead at time $t=0$. These oscillations are rapidly damped by the viscosity of the fluid but nonetheless render the results of the first simulated beat somewhat atypical. The initial configuration of the model vessel and valve leaflets are shown in figs. 4 and 5.

The geometry of the aortic root and ascending aorta is based on the geometric description of Reul et al. [13], and the dimensions of the model vessel are based on measurements by Swanson and Clark of human aortic roots harvested post autopsy that were pressurized to 120 mmHg [14]. The stiffnesses of the fibers that comprise the vessel wall are empirically determined to keep the vessel essentially fixed in place. Our model therefore neglects the significant compliance of the real aortic root, which increases in volume by approximately 35% during ejection [62]. The incorporation of a realistic description of the elasticity of the aortic root and ascending aorta into our fluid-structure interaction model remains important future work.