PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Boyce E. Griffith

arXiv:1703.09265v1cs.CE

TL;DR

The paper addresses how to simulate aortic-valve fluid-structure interaction with elastic leaflets and root mechanics under realistic physiological driving and loading. It develops an adaptive, staggered-grid immersed boundary model coupled to a reduced circulation model, and finds physiological, periodic flow rates emerging at physiological pressures over multiple cardiac cycles.

  • Problem

    Aortic valve function and dysfunction require fluid-structure interaction models that represent immersed elastic structures in viscous incompressible flow under physiological conditions.

  • Method

    An adaptive, staggered-grid immersed boundary method models elastic leaflets and a semi-rigid aortic root, with physiological pressure driving and reduced-model pressure loading.

  • Results

    Physiological flows emerge over multiple cardiac cycles with physiological pressure differences, including approximately 65 ml mean stroke volume and 12.2 mmHg maximum systolic pressure difference.

  • Takeaways & Limitations

    The coupled fluid-structure simulation produces realistic time-periodic flow rates without prescribing flow-rate boundary conditions.

  • Takeaways & Limitations

    Higher-resolution simulations are constrained by severe timestep restrictions, while model realism is also limited by simplified elasticity descriptions for the valve and root.

Abstract

from arXiv · show

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.

1. INTRODUCTION

The paper applies an adaptive, staggered-grid immersed boundary method to three-dimensional aortic valve fluid-structure interaction, motivated by the need to study valve function and dysfunction. Its formulation models leaflet and root mechanics while improving numerical treatment of large deformations and pressure discontinuities.

  • 1. INTRODUCTION: The immersed boundary method couples Lagrangian elastic-structure equations with Eulerian fluid equations through regularized Dirac delta interactions.It uses curvilinear meshes for Lagrangian equations and Cartesian grids for Eulerian equations.
  • 1. INTRODUCTION: Adaptive immersed-boundary discretizations avoid conforming fluid and structure meshes, simplifying grid generation and accommodating large elastic deformations.The Lagrangian mesh may cut through the Cartesian Eulerian grid arbitrarily, without dynamically generated body-fitted meshes.
  • 1. INTRODUCTION: The study addresses clinically relevant valve fluid mechanics, given the large number of valve procedures and the relationship between prosthetic-valve difficulties and fluid dynamics.The motivation is to enable study of fluid-mechanical mechanisms underlying valve function and dysfunction.
  • 1. INTRODUCTION: The model represents thin aortic leaflets as elastic fiber systems and the aortic root and ascending aorta as a thick, semi-rigid elastic structure.The leaflet fibers resist extension, compression, and bending, while the geometry is based on descriptions derived from imaging data.
  • 1. INTRODUCTION: The present method uses a staggered-grid discretization that improves volume conservation by one to two orders of magnitude over cell-centered discretization and better resolves pressure discontinuities.These discontinuities are especially important along closed leaflets supporting physiological pressure loads.

2. THE CONTINUOUS EQUATIONS OF MOTION

The continuous IB formulation couples Lagrangian elastic structures to an Eulerian viscous incompressible fluid through force and velocity interaction equations. The aortic-valve model represents leaflet and vessel elasticity, imposes pressure-based physiological boundary conditions, and couples downstream loading to a reduced circulation model.

  • IB formulation: The IB equations describe structure elasticity in Lagrangian form and fluid momentum, velocity, and incompressibility in Eulerian form.The interaction equations connect the two descriptions through force spreading and structure advection.
  • IB formulation: Dirac-delta interaction equations convert Lagrangian elastic force density into Eulerian force density and move material points with the local fluid velocity.This formulation enforces no slip by determining immersed-structure motion rather than imposing a separate fluid constraint.
  • Elastic structures: The aortic valve leaflets are thin elastic boundaries, while the vessel wall is a thick, semi-rigid elastic structure represented by fiber-based extension, compression, and bending resistance.The elastic force is derived from a strain-energy functional, with total energy combining stretching and bending contributions; bending resistance is used in the leaflets.
  • Boundary conditions: The model prescribes a time-dependent left-ventricular pressure waveform upstream to drive flow through the valve and uses a reduced circulation model downstream to provide dynamic pressure loading.The downstream model is a three-element Windkessel model with characteristic resistance Rc, peripheral resistance Rp, and arterial compliance C.
  • Boundary conditions: Pressure boundary conditions let the model determine the boundary velocity profile and preserve a realistic pressure difference across the valve during diastole.Prescribing flow rate instead would require an appropriate boundary velocity profile and prevent imposing a realistic diastolic pressure difference.
  • Boundary conditions: A zero-pressure open boundary allows vessel-volume changes and mismatched instantaneous inflow and outflow, while periodic steady state requires time-integrated inflow and outflow volumes to match.Because the vessel wall is semi-rigid, the instantaneous rates are approximately equal; physiological root compliance would generally produce instantaneous differences.

3. THE DISCRETE EQUATIONS OF MOTION

The discrete model combines fiber-aligned Lagrangian meshes with locally refined, staggered-grid Cartesian discretizations of the Eulerian equations. Adaptive refinement follows the moving immersed structure, while explicit bending forces impose a severe timestep restriction.

  • Lagrangian and Eulerian spatial discretizations: The Lagrangian equations use a fiber-aligned curvilinear mesh, while Eulerian equations use a block-structured, locally refined Cartesian grid adapted to the moving fibers.The Lagrangian mesh spacings are Δs1, Δs2, and Δs3, with mesh nodes carrying structure positions and elastic force densities.
  • Lagrangian and Eulerian spatial discretizations: Velocity components are stored at corresponding Cartesian cell-face centers, whereas pressure is stored at cell centers.The staggered arrangement also places the three body-force components at their corresponding cell-face centers.
  • Cartesian grid adaptive mesh refinement: The adaptive Cartesian grid is organized as nested levels with refinement ratio n, and simulations generally use two levels with ℓmax = 1 and n = 4.Grid patches satisfy proper nesting, with finer patch faces coincident with cells on the next coarser level.
  • Lagrangian and Eulerian spatial discretizations: The three-dimensional Eulerian discretization uses a locally refined staggered-grid finite-difference scheme for the incompressible Navier–Stokes equations.The coarsest grid uniformly discretizes a rectangular computational box.
  • Composite-grid operators: Composite-grid divergence, gradient, and Laplacian operators couple patchwise discretizations through restriction, interpolation, and prolongation procedures.The staggered-grid operators compute divergence at cell centers and pressure-gradient and Laplacian terms on cell faces.
  • Temporal stability: The explicit treatment of bending-resistant elastic forces imposes a severe timestep stability restriction.For structures with only extension- and compression-resistant elements, the restriction is reduced to a less severe order in Δt.

4. SOLUTION METHODOLOGY

The solution method advances the coupled fluid–structure system by solving an unsplit block system with FGMRES and a projection-based preconditioner. This preserves direct boundary-condition imposition while avoiding timestep-splitting error, although projection-based solver boundary conditions can be difficult to make accurate.

  • Krylov solution: The coupled velocity–pressure system is solved with FGMRES using prior iterates as initial approximations and the projection method as a preconditioner.The method solves for the updated structure position, velocity, and pressure within the timestep iteration.
  • Preconditioning: The projection-based preconditioner uses a pressure multigrid V-cycle and conjugate-gradient iterations for the velocity subsystem.The implementation avoids multigrid for velocity because the relevant linear system is well conditioned under the stated flow and timestep conditions.
  • Boundary conditions: Projection methods used as solvers can require carefully chosen artificial boundary conditions, and high-order accuracy may be difficult for some outflow-boundary problems.The passage identifies stable and accurate boundary-condition discretization as difficult in practice.
  • Preconditioning: Using projection as a preconditioner allows the Krylov method to eliminate timestep-splitting error from the basic projection method.The overall solve targets an unsplit discretization of the incompressible Stokes equations.
  • Boundary conditions: The unsplit coupled solve permits direct imposition of the true velocity and pressure boundary conditions.Artificial boundary conditions required by the basic projection method therefore do not determine the accuracy of the overall solver.

5. IMPLEMENTATION

The adaptive immersed boundary method is implemented in IBAMR, a freely available C++ framework supporting distributed-memory parallelism and Cartesian adaptive mesh refinement.

  • Software framework: IBAMR provides the implementation framework for adaptive immersed boundary simulations of fluid–structure interaction.IBAMR is a freely available C++ library designed for models using the immersed boundary method.
  • Software framework: The framework supports MPI-based distributed-memory parallelism and Cartesian-grid adaptive mesh refinement.It relies on SAMRAI, PETSc, and hypre for substantial functionality.

6. COMPUTATIONAL RESULTS

The adaptive immersed-boundary model simulates aortic-valve dynamics under physiological driving and loading conditions, including leaflet opening and closure over multiple cardiac cycles. It rapidly approaches periodic behavior with physiological stroke volume and transvalvular pressures, while adaptive-grid performance depends on the refinement hierarchy.

  • Model configuration: The three-dimensional model represents valve leaflets as elastic fiber systems and the aortic root and ascending aorta as a thick, semi-rigid elastic structure.Leaflet geometry and fiber architecture follow the theory of Peskin and McQueen; vessel geometry is based on prior anatomical descriptions and measurements.
  • Physiological conditions: A time-periodic ventricular pressure drives the simulation, while dynamic loading conditions are supplied by a three-element Windkessel circulation model fit to clinical data.Pressure boundary conditions are imposed at both upstream and downstream vessel boundaries.
  • Valve dynamics: Approximately 65 ml mean stroke volume and 420 ml s−1 peak flow rate are obtained during the second and third beats as the model rapidly approaches a periodic steady state.The flow rate is not prescribed; it emerges from the fluid-structure interaction simulation.
  • Valve dynamics: 12.2 mmHg and 12.0 mmHg maximum systolic pressure differences are obtained in the second and third beats, respectively, compared with an experimentally reported value of 12.8 mmHg.Mean transvalvular pressure differences are 5.5 mmHg and 5.1 mmHg in the two beats.
  • Valve dynamics: During closure, the model permits minor regurgitation but no further leak once the valve is closed.The closing sequence is shown during the second cardiac cycle.
  • Performance analysis: More than two adaptive Cartesian-grid levels do not improve efficiency at the present effective fine-grid resolution, and higher-resolution bending models may require efficient implicit time stepping.The optimal number of refinement levels is problem dependent.

7. CONCLUSIONS

The study applies an adaptive, staggered-grid immersed boundary method to multibeat aortic valve simulations with realistic driving and loading conditions. Physiological flows and pressures emerge from the coupled model, while spatial resolution and constitutive realism remain limitations.

  • Physiological flows are obtained with physiological pressure differences over multiple cardiac cycles, although only upstream and downstream pressure boundary conditions are prescribed.
  • Realistic, time-periodic flow rates emerge from the coupled fluid-structure interaction model rather than being specified directly.
  • The present simulations combine a physiological driving pressure waveform, multiple cardiac cycles, and an adaptive staggered-grid IB method.
  • Fully resolved three-dimensional aortic-valve fluid dynamics remain an unresolved goal because semi-implicit timestepping severely restricts the timestep size and higher spatial resolution.The authors identify efficient parallel implicit IB methods as a possible route beyond this stability restriction.
  • Model realism is also limited by simple elasticity descriptions of the aortic valve and root, motivating experimentally based constitutive models and more general finite-element-compatible IB extensions.

A. COMPOSITE-GRID DISCRETIZATION

The composite-grid discretization uses coarse and fine levels within an adaptive staggered-grid framework, with distinct indexing and storage conventions for velocity and pressure.

  • Coarse grid cells at level ℓ use indices (I, J, K), while cells at the next finer level ℓ+1 use indices (i, j, k).
  • The fine and coarse grid indices are related by (i, j, k) = (nI, nJ, nK), so corresponding cells share a vertex.
  • On the coarse level, u = (u, v, w) and p denote stored velocity and pressure values, with analogous notation on the fine level.

A.1. Restriction

The adaptive composite-grid scheme defines restriction operators by relating one coarse cell to the overlying n × n × n fine-grid cells.

  • Restriction procedures are defined for a coarse cell (I, J, K) and its overlying n × n × n fine-grid cells on level ℓ+1.The overlying fine cells are indexed as (i + α, j + β, k + γ) for α, β, γ = 0, …, n − 1.

A.1.1. Cell-centered cubic restriction

Cell-centered cubic restriction computes coarse-grid quantities from nearby fine-grid values, with linear interpolation used for the two-to-one refinement case.

  • The cell-centered cubic restriction procedure requires even n ≥ 4 and uses the closest overlying 4 × 4 × 4 fine-grid values.
  • When n = 2, the implementation reverts to linear interpolation, reducing the formal order of accuracy at coarse-fine interfaces.
  • Odd n is disallowed because the recursive divergence- and curl-preserving interpolation scheme requires n to be a power of two.

A.1.2. Face-centered conservative restriction

The face-centered conservative restriction procedure defines coarse-level face-centered quantities from overlying fine-grid values and preserves conservation in the finite-volume sense.

  • Coarse-level face-centered quantities are computed from values stored in the overlying n × n × n fine-grid cells.
  • The restriction procedure is conservative in the sense of a finite-volume scheme.

A.1.3. Face-centered cubic restriction

Face-centered cubic restriction uses nearby fine-grid values to define coarse-level quantities when the refinement factor is sufficiently large and even, with a linear fallback for n = 2.

  • The face-centered cubic restriction procedure requires an even refinement factor n of at least four.
  • For n = 2, the implementation reverts to linear interpolation, reducing the formal order of accuracy at coarse-fine interfaces.
  • Odd refinement factors are disallowed because the recursive divergence- and curl-preserving regridding interpolation requires n to be a power of two.

A.2. Interpolation at coarse-fine interfaces

The coarse-fine interface scheme computes ghost-cell values through staged interpolation tailored to cell-centered quantities and velocity components, then extends the construction to three dimensions with tensor-product rules.

  • The scheme is described in two dimensions because its three-dimensional extension is straightforward but cumbersome to visualize and describe.
  • For an even refinement factor, coarse-grid values are first interpolated tangentially to align with valid fine-grid locations, then ghost cells are interpolated normally.
  • Cell-centered interpolation: Cell-centered quantities use quadratic interpolation in both tangential and normal directions, combining coarse-grid and adjacent fine-grid values.
  • Face-centered velocity interpolation: Velocity components normal to the interface use quadratic interpolation tangentially and normally, while tangential components use cubic interpolation tangentially followed by quadratic interpolation normally.
  • Three-dimensional extension: In three dimensions, the tangential coarse-grid interpolation becomes two-dimensional and uses tensor-product rules combining the described schemes.
  • AMR differential operators: Discrete divergence, gradient, and Laplacian approximations on the AMR hierarchy combine restriction procedures with finite-difference formulas.
Loading 1703.09265v1…