INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING
Int. J. Numer. Meth. Biomed. Engng. 2017; 00:1–32
Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/cnm

Hybrid finite difference/finite element immersed boundary method

arXiv:1612.05916v2 [math.NA] 23 Feb 2017

Boyce E. Griffith1∗ , Xiaoyu Luo2
1 Departments

of Mathematics and Biomedical Engineering, Carolina Center for Interdisciplinary Applied
Mathematics, and McAllister Heart Institute, University of North Carolina, Chapel Hill, NC, USA
2 School of Mathematics and Statistics, University of Glasgow, Glasgow, UK

SUMMARY
The immersed boundary method is an approach to fluid-structure interaction that uses a Lagrangian
description of the structural deformations, stresses, and forces along with an Eulerian description of the
momentum, viscosity, and incompressibility of the fluid-structure system. The original immersed boundary
methods described immersed elastic structures using systems of flexible fibers, and even now, most
immersed boundary methods still require Lagrangian meshes that are finer than the Eulerian grid. This
work introduces a coupling scheme for the immersed boundary method to link the Lagrangian and Eulerian
variables that facilitates independent spatial discretizations for the structure and background grid. This
approach employs a finite element discretization of the structure while retaining a finite difference scheme
for the Eulerian variables. We apply this method to benchmark problems involving elastic, rigid, and actively
contracting structures, including an idealized model of the left ventricle of the heart. Our tests include cases
in which, for a fixed Eulerian grid spacing, coarser Lagrangian structural meshes yield discretization errors
that are as much as several orders of magnitude smaller than errors obtained using finer structural meshes.
The Lagrangian-Eulerian coupling approach developed in this work enables the effective use of these coarse
structural meshes with the immersed boundary method. This work also contrasts two different weak forms
of the equations, one of which is demonstrated to be more effective for the coarse structural discretizations
facilitated by our coupling approach. Copyright c 2017 John Wiley & Sons, Ltd.
Received . . .

KEY WORDS: immersed boundary method; fluid-structure interaction; finite element method; finite
difference method; incompressible elasticity; incompressible flow

1. INTRODUCTION
Since its introduction [1, 2], the immersed boundary (IB) method has been widely used to
simulate biological fluid dynamics and other problems in which a structure is immersed in
a fluid flow [3]. The IB formulation of such problems uses a Lagrangian description of the
deformations, stresses, and forces of the structure and an Eulerian description of the momentum,
viscosity, and incompressibility of the fluid-structure system. Lagrangian and Eulerian variables
are coupled by integral transforms with delta function kernels. When the continuous equations
are discretized, the Lagrangian equations are approximated on a curvilinear mesh, the Eulerian
equations are approximated on a Cartesian grid, and the Lagrangian-Eulerian interaction equations
are approximated by replacing the singular kernels by regularized delta functions. A major
advantage of this approach is that it permits nonconforming discretizations of the fluid and the
immersed structure. Specifically, the IB method does not require dynamically generated body-fitted
∗ Correspondence to: Phillips Hall, Campus Box 3250, University of North Carolina, Chapel Hill, NC 27599-3250, USA.
E-mail: boyceg@email.unc.edu

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls [Version: 2010/03/27 v2.00]

2

B. E. GRIFFITH AND X. Y. LUO

meshes, a property that is especially useful for problems involving large structural deformations or
displacements, or contact between structures.
In many applications of the IB method, the elasticity of the immersed structure is described by
systems of fibers that resist extension, compression, or bending [3]. Such descriptions can be well
suited for highly anisotropic materials commonly encountered in biological applications, and have
facilitated significant work in biofluid dynamics [4–16], including three-dimensional simulations of
cardiac fluid dynamics [17–28]. Fiber models are also convenient to use in practice because they
permit an especially simple discretization as collections of points that are connected by springs
or beams. The fiber-based approach to elasticity modeling also presents certain challenges. For
instance, it can be difficult to incorporate realistic shear properties into spring network models [29].
Further, fiber models often must use extremely dense collections of Lagrangian points to avoid leaks
if the models are subjected to very large deformations [14].
The fiber-based elasticity models often used with the conventional IB method are a special case
of finite-deformation structural models in which the structural response depends only on strains in a
single material direction (i.e., strains in the fiber direction). The mathematical framework of the IB
method is not restricted to fiber-based material models, however, and several distinct extensions of
this methodology have sought to use more general finite-deformation structural mechanics models.
For instance, Liu, Wang, Zhang, and co-workers developed the immersed finite element (IFE)
method [30–34], which is a version of the IB method in which finite element (FE) approximations
are used for both the Lagrangian and the Eulerian equations. Like the IB method, the IFE method
couples Lagrangian and Eulerian variables by discretized integral transforms with regularized delta
function kernels, although because the IFE method is designed to use unstructured discretizations
of the Eulerian momentum equation, the IFE method uses different families of smoothed kernel
functions than those typically used with the IB method. Devendran and Peskin proposed an energy
functional-based version of the conventional IB method that obtains a nodal approximation to the
elastic forces generated by an immersed hyperelastic material via an FE-type approximation to
the deformation of the material [35, 36], and Boffi, Costanzo, Gastaldi, Heltai, and co-workers
developed a fully variational IB method that avoids regularized delta functions altogether [37–39].
Other work led to the development of the immersed structural potential method, which uses a
meshless method to describe the mechanics of hyperelastic structures immersed in fluid [40, 41].
In this paper, we describe an alternative approach to using FE structural discretizations with the
IB method that combines a Cartesian grid finite difference method for the incompressible NavierStokes equations with a nodal FE method for the structural mechanics. We consider both flexible
structures with a hyperelastic material response, and also rigid immersed structures in which rigidity
constraints are approximately imposed via a simple penalty method. The primary contribution
of this work is its treatment of the equations of Lagrangian-Eulerian interaction. Conventionally,
structural forces are spread directly from the nodes of the Lagrangian mesh using the regularized
delta function, and velocities are interpolated directly from the Eulerian grid to the Lagrangian mesh
nodes. This discrete Lagrangian-Eulerian coupling approach is adopted by the classical IB method
as well as the IFE method [30–34], the energy-based method of Devendran and Peskin [35, 36],
and the immersed structural potential method [40, 41]. A significant limitation of this approach
is that if the physical spacing of the Lagrangian nodes is too large in comparison to the spacing
of the background Eulerian grid, severe leaks will develop at fluid-structure interfaces. A widely
used empirical rule that generally prevents such leaks is to require the Lagrangian mesh to be
approximately twice as fine as the Eulerian grid [3], potentially requiring very dense structural
meshes, especially for applications that involve large structural deformations.
The Lagrangian-Eulerian coupling scheme introduced in this paper overcomes this longstanding
limitation of the classical IB method. Specifically, rather than spreading forces from the nodes of the
Lagrangian mesh and interpolating velocities to those mesh nodes, we instead spread forces from
and interpolate velocities to dynamically selected interaction points located within the Lagrangian
structural elements. Here, the interaction points are constructed by quadrature rules defined on the
structural elements. In contrast to both the conventional IB method and the IFE method, in our
method, the nodes of the structural mesh act as control points that determine the positions of the

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

3

interaction points, but the mesh nodes are not required to be locations of direct Lagrangian-Eulerian
interaction.† Our approach thereby takes full advantage of the kinematic information provided by the
FE description of the structural deformation, which yields approximations not only to the positions
of the nodes of the Lagrangian mesh, but also to the positions of all material points of the structure.
A related IB method was introduced by Shankar et al. [42], who use a radial basis function
(RBF) approach to describe the deforming geometry of two-dimensional models of circulating
platelets, which are described as elliptical shells with linear elasticity models. In that approach,
velocities are interpolated to data sites, which are analogous to the control points of the present
FE scheme, but forces are spread from a fixed collection of sample sites that are distributed along
the surface of the platelet. Unlike the RBF scheme, the interaction operators of the present FE
approach satisfy a discrete adjoint property that implies that the method satisfies energy conservation
during Lagrangian-Eulerian interaction. As demonstrated by Shankar et al., satisfying the discrete
adjoint property is not required to obtain a convergent scheme, but this property is crucial for
some applications, such as the design of efficient IB methods for rigid bodies [43, 44]. The
importance of satisfying the adjoint property also seems to depend on the choice of regularized
kernel functions. Lagrangian-Eulerian coupling schemes that do not satisfy the adjoint property
seem to be prone to induce numerical instabilities when used with kernel functions developed by
Peskin that satisfy moment conditions along with a “sum-of-squares” condition [3], like those used
in this work. Shankar et al. consider a kernel function that may be less sensitive to these aspects of
the discretization.
We apply the present IB method to benchmark fluid-structure interaction (FSI) problems
involving elastic, rigid, and actively contracting structures, including an idealized model of the
left ventricle of the heart [45]. For elastic structures, we consider two weak formulations of the
equations of motion suitable for standard nodal (C 0 ) FE methods for structural mechanics. One of
these formulations, referred to herein as the unified weak form, is similar to those used by earlier
IB-like methods [30–34,37–39]. This formulation uses a single volumetric force density to describe
the mechanical response of the immersed structure. The other formulation, referred to herein as the
partitioned weak form, uses both an internal volumetric force density, which is supported throughout
the immersed elastic structure, and also a transmission surface force density, which is supported only
on the surface of the immersed structure. In a numerical scheme, the unified formulation effectively
regularizes the transmission surface force density by projecting it onto the volumetric FE shape
functions. The partitioned formulation described herein, which does not appear to have been used
in previous versions of the IB or IFE methods, avoids this additional regularization step by treating
the surface and volume force densities separately. We present numerical tests that demonstrate that
this partitioned scheme can yield accurate results even for Lagrangian meshes that are significantly
coarser than the background Eulerian grid.
An important feature of our discretization approach is that obtaining a “watertight” structure
simply requires using a dense enough collection of interaction points so as to prevent leaks.
Moreover, because it is straightforward to use dynamic quadrature schemes that account for highly
deformed elements, this approach can ensure that the structure remains watertight even in the
presence of very large structural deformations. One benefit of our approach is that the Lagrangian
resolution may be determined primarily by accuracy requirements for the structural model rather
than by the requirements of the Lagrangian-Eulerian coupling scheme. Further, for problems
involving immersed rigid structures, we demonstrate that using coarser Lagrangian meshes can
reduce discretization errors by an order of magnitude compared to finer Lagrangian meshes for
a given Eulerian grid spacing. The present method also yields improved volume conservation in
comparison to the IFE method [33] for both coarse and fine structural meshes. To our knowledge,
the present IB method is the first IFE-type method to explicitly enable the effective use of such
relatively coarse Lagrangian meshes.
† In fact, the structural mesh nodes can be both control points and interaction points if the quadrature rules that generate
the interaction points include the mesh nodes as quadrature points, as in Gauss-Lobatto rules.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

4

B. E. GRIFFITH AND X. Y. LUO

Ω
χ(U, t)

U
χ : (U, t) 7→ Ω

Figure 1. Lagrangian and Eulerian coordinate systems. The Lagrangian material coordinate domain is U ,
and the Eulerian physical coordinate domain is Ω. The physical position of material point X at time t is
χ(X, t), the physical region occupied by the structure is χ(U, t), and the physical region occupied by the
fluid is Ω \ χ(U, t).

2. CONTINUOUS FORMULATIONS
2.1. Immersed elastic bodies
In the IB formulation of problems involving an immersed elastic body, the momentum, velocity,
and incompressibility of the coupled fluid-structure system are described in Eulerian form, whereas
the deformation and elastic response of the immersed structure are described in Lagrangian form.
A similar formulation is used in the case of an immersed rigid structure; see sec. 2.2. Let
x = (x1 , x2 , . . .) ∈ Ω ⊂ Rd , d = 2 or 3, denote Cartesian physical coordinates, with Ω denoting
the physical region that is occupied by the coupled fluid-structure system; let X = (X1 , X2 , . . .) ∈
U ⊂ Rd denote Lagrangian material coordinates that are attached to the structure, with U denoting
the Lagrangian coordinate domain; and let χ(X, t) ∈ Ω denote the physical position of material
point X at time t. The physical region occupied by the structure at time t is χ(U, t) ⊆ Ω, and the
physical region occupied by the fluid at time t is Ω \ χ(U, t). See fig. 1. We do not assume that the
Lagrangian coordinates are the initial coordinates of the elastic structure, nor, more generally, do
we require that U ⊆ Ω.
To use an Eulerian description of the fluid and a Lagrangian description of the elasticity of
the immersed structure, it is necessary to describe the stress of the fluid-structure system in both
Eulerian and Lagrangian forms. Let σ(x, t) denote the Cauchy stress tensor of the coupled fluidstructure system. In the present formulation, we assume that
(
σ e (x, t) for x ∈ χ(U, t),
σ(x, t) = σ f (x, t) +
(1)
0
otherwise,
in which σ f (x, t) is the stress tensor of a viscous incompressible fluid and σ e (x, t) is the stress
tensor that describes the elastic response of the immersed structure. The fluid stress tensor is the
usual one for a viscous incompressible fluid,
h
i
σ f (x, t) = −p(x, t) I + µ ∇u(x, t) + (∇u(x, t))T ,
(2)

in which p(x, t) is the pressure, µ is the dynamic viscosity, and u(x, t) is the Eulerian velocity field.
To describe the elasticity of the structure with respect to the Lagrangian material coordinate
system, we use the first Piola-Kirchhoff elastic stress tensor Pe (X, t), which is defined in terms
of σ e via
Pe (X, t) = J(X, t) σ e (χ(X, t), t) F−T (X, t),
(3)
in which the deformation gradient tensor associated with the deformation χ : (U, t) 7→ Ω is
F(X, t) =

∂χ
(X, t),
∂X

(4)

and J(X, t) = det(F(X, t)) is the Jacobian determinant of the deformation gradient. Although Pe is
only defined within the solid region, it is convenient to extend σ e (x, t) by zero in the fluid region.
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

5

For simplicity, we primarily restrict our attention to hyperelastic constitutive models, which may
be characterized by a strain-energy functional W e (F). For such constitutive laws,
Pe =

∂W e
.
∂F

(5)

Our formulation does not rely on the existence of such an energy functional, however, and it permits
a material description defined only in terms of a stress response. For instance, separate work using
the present methodology relies on this feature to treat active tension generation in dynamic models
of muscle contraction [45–48].
2.1.1. Strong formulation As shown by Boffi et al. [37], the strong form of the equations of motion
is
Du
(x, t) = −∇p(x, t) + µ∇2 u(x, t) + f e (x, t),
Dt
∇ · u(x, t) = 0,
Z
e
f (x, t) =
∇X · Pe (X, t) δ(x − χ(X, t)) dX
U
Z
−
Pe (X, t) N(X) δ(x − χ(X, t)) dA(X),
Z ∂U
∂χ
(X, t) =
u(x, t) δ(x − χ(X, t)) dx,
∂t
Ω
ρ

(6)
(7)
(8)

(9)

∂u
in which ρ is the mass density of the coupled fluid-structure system, Du
Dt (x, t) = ∂t (x, t) + u(x, t) ·
e
∇u(x, t) is the material derivative, f (x, t) is the Eulerian elastic force density, and δ(x) =
Qd
i=1 δ(xi ) is the d-dimensional delta function.
Two different types of Lagrangian elastic force densities appear in these equations. The
Lagrangian internal force density, ∇X · Pe , is a volumetric force density that is distributed
throughout the elastic body, whereas the Lagrangian transmission force density, −Pe N, is a surface
force density that is distributed along the fluid-solid interface. Eq. (8) generates a corresponding
Eulerian description of the elastic forces from the volumetric and surface forces. Notice that f e (x, t)
will generally have a δ -layer of force along the fluid-solid interface. It is possible to show that f e
is variationally equivalent to ∇ · σ e (for instance, by integrating against a test function). Doing so,
it is clear that f e is singular along the fluid-solid interface wherever σ e n is discontinuous. Indeed,
because elastic stresses are present only within the structure, σ e n is generally discontinuous at the
fluid-structure interface. These discontinuities must be exactly balanced by discontinuities in σ f n to
ensure that the total stress, σ = σ f + σ e , has a continuous traction vector. If Pe (X, t) is sufficiently
smooth, the internal force acts as a regular (i.e., nonsingular) body force and may be treated with
higher-order accuracy by the IB method [37, 49]. However, the transmission force always acts as a
singular force layer, and although this force will induce jumps in the pressure and shear stress along
∂χ(U, t), such force layers may also be readily treated by the IB method.
An integral transform is also used in eq. (9) to determine the velocity of the immersed elastic
structure from the material velocity field u(x, t). The defining property of δ(x) implies that eq. (9)
is equivalent to
∂χ
(X, t) = u(χ(X, t), t),
(10)
∂t
which may be interpreted as the no-slip and no-penetration conditions of a viscous incompressible
fluid. Notice, however, that the no-slip and no-penetration conditions do not appear as constraints
on the fluid motion. Instead, these conditions determine the motion of the immersed structure.

2.1.2. Weak formulations To use standard nodal (C 0 ) FE methods for nonlinear elasticity with
the IB formulation, it is necessary to introduce a weak formulation of the equations of motion.
Here we consider two different formulations that each employ a weak form of the Lagrangian
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

6

B. E. GRIFFITH AND X. Y. LUO

equations. Because we use a finite difference scheme to approximate the incompressible NavierStokes equations, we do not develop a weak formulation for the Eulerian equations or the equations
of Lagrangian-Eulerian interaction.
We begin with the partitioned weak formulation of the equations. To do so, we first define the
internal elastic force density, F(X, t), that is variationally equivalent to ∇X · Pe by requiring
Z
Z
F(X, t) · V(X) dX = −
Pe (X, t) : ∇X V(X) dX
(11)
U
U
Z
+
(Pe (X, t) N(X)) · V(X) dA(X)
∂U

to hold for arbitrary smooth V(X). Integrating by parts, it is clear that eq. (11) implies that
Z
Z
F(X, t) · V(X) dX =
(∇X · Pe (X, t)) · V(X) dX,
U

(12)

U

for all smooth V(X). Assuming sufficient regularity, F = ∇X · Pe pointwise. It is also convenient
to define the transmission force density, T(X, t), pointwise along ∂U via
T(X, t) = −Pe (X, t) N(X).

(13)

With F and T so defined, the equations of motion become
Du
(x, t) = −∇p(x, t) + µ∇2 u(x, t) + f (x, t) + t(x, t),
Dt
∇ · u(x, t) = 0,
Z
f (x, t) =
F(X, t) δ(x − χ(X, t)) dX,
ZU
t(x, t) =
T(X, t) δ(x − χ(X, t)) dA(X),
Z∂U
∂χ
(X, t) =
u(x, t) δ(x − χ(X, t)) dx,
∂t
Ω
ρ

(14)
(15)
(16)
(17)
(18)

in which f (x, t) is the Eulerian internal elastic force density and t(x, t) is the Eulerian transmission
elastic force density. It is important to notice that, under relatively mild regularity requirements, F
and T are both smooth functions on their domains of definition, and f is piecewise smooth. Under
these conditions, the only singular function in this formulation is t.
Alternative weak definitions of the Lagrangian elastic force density are possible. The formulation
typically used with the IFE method defines a total elastic force per unit volume, G(X, t), by
requiring
Z
Z
G(X, t) · V(X) dX = −
Pe (X, t) : ∇X V(X) dX
(19)
U

U

for arbitrary V(X). It is possible to see that G(X, t) accounts for the effects of both the internal
and transmission force densities of the strong form of the equations. To do so, we again integrate by
parts to find that
Z
Z
G(X, t) · V(X) dX =
(∇ · Pe (X, t)) · V(X) dX
(20)
U
U
Z
−
(Pe (X, t) N(X)) · V(X) dA(X),
∂U

for all V(X). Thus, G(X, t) can be a continuous function only if Pe N ≡ 0. In general, G is in fact
a distribution that accounts for both the internal force per unit volume and the transmission force
per unit area, in which the transmission force gives rise to a singular force layer concentrated along
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

7

∂U . It is possible to use eq. (19) with standard finite element methods; however, when doing so, the
transmission surface force density is effectively projected (in an L2 sense) onto the volumetric basis
functions. Specifically, the FE basis functions serve to regularize the transmission force.
Using this definition of G, we may state a unified weak formulation that includes only a single,
unified body forcing term:
Du
(x, t) = −∇p(x, t) + µ∇2 u(x, t) + g(x, t),
Dt
∇ · u(x, t) = 0,
Z
g(x, t) =
G(X, t) δ(x − χ(X, t)) dX,
ZU
∂χ
(X, t) =
u(x, t) δ(x − χ(X, t)) dx,
∂t
Ω
ρ

(21)
(22)
(23)
(24)

in which g(x, t) is the Eulerian total elastic force density. This weak form of the equations of motion
is essentially the formulation employed in the IFE method [30–34], the energy-based formulation
of Devendran and Peskin [35, 36], and the fully variational IB method [37–39]. To our knowledge,
the partitioned formulation described here has not been widely used in previous work.
2.2. Immersed rigid structures
The equations of motion for a fixed, rigid structure are similar to those used in the case of an
immersed elastic structure:
Du
(x, t) = −∇p(x, t) + µ∇2 u(x, t) + f (x, t),
Dt
∇ · u(x, t) = 0,
Z
f (x, t) =
F(X, t) δ(x − χ(X, t)) dX,
ZU
∂χ
(X, t) =
u(x, t) δ(x − χ(X, t)) dx = 0,
∂t
Ω
ρ

(25)
(26)
(27)
(28)

in which here, F(X, t) is a Lagrange multiplier for the constraint ∂χ
∂t ≡ 0. Thus, this fully
constrained formulation takes the form of an extended saddle-point problem with two Lagrange
multipliers, p(x, t) for the incompressibility constraint and F(X, t) for the rigidity constraint.
Solving this system effectively requires specialized techniques that are the subject of active research
[43, 44]. In this work, we consider instead a penalty formulation for an immersed stationary
structure, in which the Lagrange multiplier force is approximated by
F(X, t) = κ (χ(X, 0) − χ(X, t)) − η

∂χ
(X, t),
∂t

(29)

in which κ ≥ 0 is a stiffness penalty parameter and η ≥ 0 is a damping penalty parameter. As
κ → ∞, χ(X, t) → χ(X, 0) and ∂χ
∂t → 0, so that this formulation is equivalent to the constrained
formulation. In principle, it is not necessary to include the damping parameter η ; however, we have
found that using damping reduces numerical oscillations, especially at moderate-to-high Reynolds
numbers.

3. SPATIAL DISCRETIZATION
For simplicity, we describe the numerical scheme in two spatial dimensions. The extension of the
numerical scheme to the case d = 3 is straightforward.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

8

B. E. GRIFFITH AND X. Y. LUO

1)
(i, j + 2

(i − 1
, j)
2

(i, j)

1 , j)
(i + 2

(i, j − 1
)
2

Figure 2. Schematic of the staggered-grid layout of Eulerian degrees of freedom in two spatial dimensions.
The pressures are approximated at cell centers, indicated by (i, j), the x1 -components of the velocity and
force are approximated on the x1 -edges, (i − 12 , j) or (i + 21 , j), and the x2 -components of the velocity and
force are approximated on the x2 -edges, (i, j − 21 ) or (i, j + 12 ).

3.1. Eulerian discretization
We use a staggered-grid finite difference scheme to discretize the incompressible Navier-Stokes
equations in space. Compared to collocated discretizations (i.e., purely cell- or node-centered
schemes), such Eulerian discretization approaches yield superior accuracy when used with the
conventional IB method [50]. To simplify the exposition, assume that Ω is the unit square and
is discretized using a regular N × N Cartesian grid with grid spacings ∆x1 = ∆x2 = h = N1 .
Let (i, j) label the individual Cartesian grid cells for integer values of i and j , 0 ≤ i, j < N .
The components of the Eulerian velocity field u = (u1 , u2 ) are approximated at the centers
  of
the x1 - and x2 -edges of the Cartesian grid cells, i.e., at positions xi− 21 ,j = ih, j + 21 h and


xi,j− 21 = i + 21 h, jh , respectively. A staggered scheme is also used for the Eulerian body
force f = (f1 , f2 ). We use the notation (u1 )i− 21 ,j , (u2 )i,j− 21 , (f1 )i− 12 ,j , and (f2 )i,j− 21 to denote
the discrete values of u and f . The pressure p is approximated at the centers of the Cartesian grid
cells, and its discrete values are denoted pi,j . See fig. 2.
Let ∇h · u ≈ ∇ · u, ∇h p ≈ ∇p, and ∇2h u ≈ ∇2 u respectively denote standard second-order
accurate finite difference approximations to the divergence, gradient, and Laplace operators [51]. In
this approach, ∇h · u is defined at cell centers, whereas both ∇h p and ∇2h u are defined at cell edges.
We use a staggered-grid version [50, 51] of the xsPPM7 variant [52] of the piecewise parabolic
method (PPM) [53] to discretize the nonlinear advection terms. Where needed, physical boundary
conditions are discretized and imposed along ∂Ω as described by Griffith [51]. In some numerical
examples, we use a locally refined Eulerian discretization that employs Cartesian grid adaptive mesh
refinement (AMR) following the discretization approach described by Griffith [26].
3.1.1. Eulerian inner products If u and v are discrete staggered-grid vector fields, we denote by
[u] and [v] the corresponding vectors of grid values. If Ω has periodic boundaries, we define the
discrete L2 inner product on Ω by
(u, v)x = [u]T [v] h2 .
(30)
Minor adjustments to this definition are required when non-periodic physical boundary conditions
are used [51] or near coarse-fine interfaces in locally refined grids [26].
3.2. Lagrangian discretization
Let Th = ∪e U e be a triangulation of U composed of elements U e , in which e indexes the elements
M
of the mesh. We denote by {Xl }M
l=1 the nodes of the mesh, and by {φl (X)}l=1 nodal (Lagrangian)
FE basis functions. The time-dependent physical positions of the nodes of the Lagrangian mesh are
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

9

{χl (t)}M
l=1 , and, using the FE basis functions, we define an approximation to χ(X, t) by
χh (X, t) =

M
X

χl (t) φl (X).

(31)

l=1

An approximation to the deformation gradient is given by
M

X
∂χh
∂φl
Fh (X, t) =
(X, t) =
(X).
χl (t)
∂X
∂X

(32)

l=1

3.2.1. Immersed elastic structures For an immersed elastic structure, we use the FE approximation
to the deformation gradient tensor Fh (X, t) to compute Peh (X, t) and Th (X, t), which are
approximations to the first Piola-Kirchhoff stress tensor and the Lagrangian transmission force
density, respectively. We approximate the Lagrangian force densities F(X, t) and G(X, t) by
Fh (X, t) =

M
X

Fl (t) φl (X), and

(33)

Gl (t) φl (X).

(34)

l=1

Gh (X, t) =

M
X
l=1

M
e
The nodal values {Fl }M
l=1 and {Gl }l=1 must be determined from Ph (X, t) via discretizations of
eq. (11) and eq. (19). We use a standard Galerkin approach, so that after rearranging terms, eq. (11)
becomes

Z
M Z
X
φl (X) φm (X) dX Fl (t) = −
Peh (X, t) ∇X φm (X) dX
l=1

U

U

+

Z

Peh (X, t) N(X) φm (X) dA(X),

(35)

∂U

for each m = 1, . . . , M . Similarly, eq. (19) becomes
M Z
X
l=1

U


Z
φl (X) φm (X) dX Gl (t) = −
Peh (X, t) ∇X φm (X) dX,

(36)

U

for each m = 1, . . . , M . In practice, these integrals are approximated via Gaussian quadrature.
3.2.2. Immersed rigid structures For a fixed, rigid immersed structure, we directly evaluate the
discretized Lagrangian penalty force Fh (X, t) via
Fh (X, t) = κ (χh (X, 0) − χh (X, t)) − η

∂χh
(X, t).
∂t

(37)

Because χh (X, t) and Fh (X, t) are defined in terms of the same basis functions, Fh (X, t) is given
by
X
Fh (X, t) =
Fl (t) φl (X)
(38)
l

in which

Fl (t) = κ (χl (0) − χl (t)) − η

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

dχl
(t).
dt

(39)

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

10

B. E. GRIFFITH AND X. Y. LUO

3.2.3. Lagrangian inner products Letting [F] denote the vector of nodal coefficients of Fh , we
write eq. (35) as
[M] [F] = [B],
(40)
R
in which [M] is the mass matrix that has entries of the form U φl (X) φm (X) dX. Eq. (36) may
2
be rewritten similarly. The mass matrix [M] can also be used to
Pevaluate the L inner product
of
P Lagrangian functions on U . In particular, for any Uh (X, t) = l Ul (t) φl (X) and Vh (X, t) =
l V l (t) φl (s),
(Uh , Vh )X = [U]T [M] [V].
(41)
Different choices of mass matrices (e.g., lumped mass matrices) induce different discrete inner
products on U .
To simplify notation, in the remainder of this paper we drop the subscript “h” from our numerical
approximations to the Lagrangian variables.
3.3. Lagrangian-Eulerian interaction
We next describe Lagrangian-Eulerian coupling operators that take advantage of the kinematic
information encoded in the FE approximation to the deformation of the immersed structure. As
in the conventional IB method, we approximate the singular delta function kernel appearing in the
Lagrangian-Eulerian interaction equations by a smoothed d-dimensional Dirac delta function δh (x)
Qd
that is of the tensor-product form δh (x) = i=1 δh (xi ). Except where otherwise noted, in this work,
we take the one-dimensional smoothed delta function δh (x) to be the four-point delta function of
Peskin [3].
To compute an approximation to f = (f1 , f2 ) on the Cartesian grid, we construct for each
e
element U e ∈ Th a Gaussian quadrature rule with N e quadrature points XeQ ∈ U e and weights ωQ
,
e
Q = 1, . . . , N . We then compute f1 and f2 on the edges of the Cartesian grid cells via
e

N
X X

(f1 )i− 21 ,j =

U e ∈T

h

e
F1 (XeQ , t) δh (xi− 21 ,j − χ(XeQ , t)) ωQ
, and

(42)

e
F2 (XeQ , t) δh (xi,j− 21 − χ(XeQ , t)) ωQ
,

(43)

Q=1
e

(f2 )i,j− 21 =

N
X X

U e ∈Th Q=1

with F(X, t) = (F1 (X, t), F2 (X, t)). We use the notation
f (x, t) = S (χ(·, t)) F(X, t),

(44)

in which S (χ(·, t)) is the force-prolongation operator implicitly defined by eqs. (42) and (43).
Notice that
X Z
(f1 )i− 21 ,j =
(45)
F1 (X, t) δh (xi− 12 ,j − χ(X, t)) dX + O(∆Xq ), and
U e ∈Th

(f2 )i,j− 21 =

X Z

U e ∈T

h

Ue

Ue

F2 (X, t) δh (xi,j− 12 − χ(X, t)) dX + O(∆Xq ),

(46)

in which ∆X is proportional to the Lagrangian mesh spacing and O(∆Xq ) corresponds to
quadrature error that may be controlled by the choice of numerical quadrature rules.
A corresponding velocity-restriction operator J (χ(·, t)) determines the motion of the nodes of
the Lagrangian mesh from the Cartesian grid velocity field via
∂χ
(X, t) = J (χ(·, t)) u(x, t).
∂t

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

(47)

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

11

There are many possible ways to construct J ; however, we have found that an effective approach is
to require ∂χ
∂t = J u to satisfy the discrete power identity,


∂χ
F,
= (f , u)x ,
(48)
∂t X
which implies that the semi-discrete unified formulation conserves energy during LagrangianEulerian interaction [3]. This power identity can be rewritten as
(F, J u)X = (S F, u)x ,

(49)

i.e., J = S ∗ .
To construct J explicitly, it is convenient to use matrix notation. Identifying [S] and [J ] with the
matrix representations of S and J , we have that
[f ] = [S] [F], and

(50)

d[χ]
= [J ] [u].
dt

(51)

Eq. (49) then becomes
T

[F]T [M] [J ] [u] = ([S] [F]) [u] h2 .

(52)

If eq. (52) is to hold for any [F] and [u], then [J ] must be defined via
[J ] = [M]−1 [S]T h2 .

(53)

In our time-stepping scheme, which is stated below in appendix A, notice that we need only to apply
[J ] to discrete velocity fields defined on the Cartesian grid. Specifically, we do not need to compute
[J ] explicitly.
It is straightforward to show that this construction of J implies that ∂χ
∂t (X, t) is an approximation
to the L2 projection of the Lagrangian vector field UIB (X, t) = (U1IB (X, t), U2IB (X, t)), with
X
(u1 )i− 1 ,j δh (xi− 12 ,j − χ(X, t)) h2 , and
U1IB (X, t) =
(54)
2

i,j

U2IB (X, t) =

X

(u2 )i,j− 1 δh (xi,j− 12 − χ(X, t)) h2 .
2

(55)

i,j

Because the components of UIB (X, t) are not generally linear combinations of the FE basis
IB
functions, generally ∂χ
∂t 6= U .
For the semi-discretized partitioned formulation, f is computed on the Cartesian grid via f = S F.
The Eulerian transmission force density t is computed in a similar manner, but in this case, the
numerical integration occurs only over those element boundaries that coincide with ∂U . We use the
notation t = S ∂U T to denote this operation. We use the same regularized delta function δh (x) to
construct both S and S ∂U . For simplicity, we use the same velocity-restriction operator for both
formulations. This choice ensures that the two formulations coincide whenever T ≡ 0.
The Lagrangian-Eulerian interaction operators introduced in this work are different from
analogous operators generally used in standard IB methods. Standard IB methods and schemes such
as the IFE method use regularized delta functions to apply nodal forces directly to the Cartesian grid
and to interpolate Cartesian grid velocities directly to the Lagrangian nodes [3]. In such schemes,
f (x, t) would be approximated on the Eulerian grid by expressions similar to
(f1IB )i− 12 ,j

=

M
X

(F1 )l (t) δh (xi− 12 ,j − χl (t)) ωlIB , and

(56)

(F2 )l (t) δh (xi,j− 12 − χl (t)) ωlIB ,

(57)

l=1

(f2IB )i,j− 21

=

M
X
l=1

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

12

B. E. GRIFFITH AND X. Y. LUO

A.

B.

C.

Figure 3. Prolonging the elastic force density from the Lagrangian mesh onto the Eulerian grid. Starting with
an approximation to the force density at the nodes of the Lagrangian mesh (A), we use the interpolatory FE
basis functions to determine the force density at interaction points that are defined by a quadrature rule (B),
and then spread those forces from the interaction points to the background Eulerian grid using the smoothed
delta function δh (x) (C). This approach permits Lagrangian meshes that are significantly coarser than the
Eulerian grid so long as the “net” of interaction points is sufficiently dense. Denser nets of interaction points
can be obtained, for instance, by increasing the order of the numerical quadrature scheme.

in which f IB denotes the Eulerian force determined by the standard IB force-spreading operator and
ωlIB denotes the volume associated with Lagrangian node l. In this approach, each nodal force Fl is
applied only to Cartesian grid cells within the support of δh (x − χl ), and the Lagrangian mesh must
therefore be finer than the Cartesian grid to avoid leaks [3]. The corresponding approach to velocity
IB
l
restriction used by such methods would be to set dχ
dt = U (Xl , t).
Our Lagrangian-Eulerian interaction operators communicate between a collection of Lagrangian
control points (the nodes of the structural mesh) and the Cartesian grid via a collection of interaction
points (the Lagrangian quadrature points). The force-prolongation operator can be seen as the
composition of two operations: first, the Lagrangian force density is evaluated at the interaction
points in terms of data defined at the control points; then, the standard IB delta function δh (x)
spreads volume- or area-weighted force densities from the interaction points to the Cartesian grid.
See fig. 3. Velocity restriction is similar: First, the Cartesian velocity field is interpolated to the
interaction points using δh (x); then, these velocities are accumulated to form the right-hand-side of
a system of equations that determines ∂χ
∂t at the control points. Our approach is similar to methods
used in RBF-based IB methods [42].
In general, it is necessary that the same interaction points are used in the discrete force-spreading
and velocity-interpolation operators if those operators are to satisfy a discrete adjoint property. It is
possible to construct an interpolation operator that uses the control points as the interaction points,
but then satisfying the discrete adjoint property requires that the control points are also used as the
interaction points in the spreading operator.
Our numerical tests indicate that in our scheme, the Lagrangian structure appears watertight to the
fluid so long as the net of interaction points is sufficiently dense. Denser nets of interaction points
can be obtained by increasing the order of the quadrature rule, and this may be done in an adaptive
manner as the simulation progresses. In our numerical tests, we use dynamically adapted Gaussian
quadrature rules to construct S and J that provide, on average, at least a 3 × 3 net of quadrature
points per Cartesian grid cell.

4. IMPLEMENTATION
This version of the IB method is implemented in the open-source IBAMR software [54], which
is a C++ framework for developing FSI models that use the IB method. IBAMR provides support
for distributed-memory parallelism and adaptive mesh refinement (AMR). IBAMR relies upon the
SAMRAI [55–57], PETSc [58–60], hypre [61, 62], and libMesh [63, 64] libraries for much of its
functionality.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

13

5. NUMERICAL RESULTS
5.1. Thick elliptical shell
This section presents results from tests that use a thick elastic shell [23, 37, 49] to demonstrate the
convergence properties of our method for different types of material models. In these computations,
the physical domain is Ω = [0, 1] × [0, 1] with periodic boundary conditions, and, following Boffi et
al. [37], the Lagrangian coordinate domain U is defined using curvilinear coordinates s = (s1 , s2 ) ∈
U instead of reference coordinates, with U = [0, 2πR] × [0, w] for R = 0.25 and w = 0.0625, and
with periodic boundary conditions in the s1 direction. The coordinate mapping χ : (U, t) 7→ Ω is
given at time t = 0 by
χ(s, 0) = (cos(s1 /R)(R + s2 ) + 0.5, sin(s1 /R)(R + γ + s2 ) + 0.5) .

(58)

For γ = 0, the initial configuration of the shell is a circular annulus with inner radius R and
thickness w, which corresponds to an equilibrium configuration of the structure. For γ 6= 0, the
initial configuration is an elliptical annulus that is out of equilibrium. In our tests, we use γ = 0 for
static problems and γ = 0.15 for dynamic problems. In either case, we discretize Ω using an N -byN Cartesian grid. The Lagrangian discretization is constructed so that the nodes of the Lagrangian
mesh are physically separated by a distance of approximately Mfac ∆x. Specifically, we discretize
N
U using a 28M -by-M mesh of bilinear (Q1 ) elements, with M = M1fac 16
. Representative numerical
results using N = 128 are shown in figs. 4 and 7.
Although this is a relatively simple benchmark problem, the static version (γ = 0) is one of the
only test problems available for the IB method that we know of that permits a simple analytic
solution. Moreover, because certain choices of Pe allow the IB method to attain higher-order
convergence rates, this test case allows us to verify that our implementation attains its designed
order of accuracy.
5.1.1. Anisotropic shell We first consider an idealized anisotropic material model defined in terms
of
2
µe
µe ∂χ
Fα1 Fα1 ,
(59)
=
W e (F) =
2w ∂s1
2w
so that
∂W e
µe
Pe =
=
∂F
w

∂χ1
∂s1
∂χ2
∂s1

!

0
µe F11
=
0
w F21


0
.
0

(60)

This model corresponds to an idealized elastic material composed of a continuous family of
extension-resistant fibers that wrap the thick shell. Because U is periodic in the s1 direction,
Pe N ≡ 0 along ∂U . If we view the structure as a fiber-reinforced material, none of the fibers
terminate along the boundary of the structure. Because the transmission force vanishes in this case,
the unified and partitioned weak formulations are identical.
R
Setting γ = 0, so that the structure is in equilibrium, and requiring Ω p(x, t) dx = 0, it can be
shown [37] that

µe

r ≤ R,
p0 + Re
µ 1
p(x, t) = p0 + w R (R + w − r) R < r ≤ R + w,
(61)

p
R + w < r,
0


3
e
R2 − (R+w)
. We set ρ = 1, µ = 1, and µe = 1, and we
with r = kx − (0.5, 0.5)k and p0 = πµ
3w
R
consider the time interval 0 ≤ t ≤ 3. Fig. 5 summarizes the error data at time t = 3 for N = 64, 128,
256, 512, and 1024, using Mfac = 1, 2, and 4, with ∆t = 0.25∆x. Second-order convergence rates
are observed in the L1 , L2 , and L∞ norms for the velocity field. Second-order convergence rates are
also observed for the pressure in the L1 norm; however, because the pressure field is C 0 but not C 1 ,
only first-order convergence rates are observed for the pressure in the L∞ norm, and intermediate
convergence rates (approximately order 1.5) are observed in the L2 norm.
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

14

B. E. GRIFFITH AND X. Y. LUO

A.

B.

C.

Figure 4. Representative results from the dynamic (γ = 0.15) version of the idealized anisotropic shell
model of sec. 5.1.1 for N = 128 over the time interval 0 ≤ t ≤ 0.75. The computed pressure and structure
deformation for Mfac = 4 are shown in A and B. The computed deformations obtained with Mfac = 1 and
Mfac = 4 are compared in C. The coarse and fine structural meshes yield essentially identical kinematics.

We also consider the case in which γ = 0.15, so that the initial configuration of the shell is not
in equilibrium. We set ρ = 1, µ = 0.01, and µe = 1, yielding a Reynolds number of approximately
50. We consider the time interval 0 ≤ t ≤ 0.75, which corresponds to approximately one damped
oscillation of the shell. Because an exact solution is not available in this case, we use a Richardson
extrapolation approach, as described in detail in previous work [49]. Fig. 6 summarizes the error
data at time t = 0.75 for N = 64, 128, 256, and 512 and Mfac = 1, 2, and 4, with ∆t = 0.25∆x.
Essentially second-order convergence rates are observed in the L1 , L2 , and L∞ norms for the
velocity field. The pressure converges at a second-order rate in the L1 norm, at a first-order rate
in the L∞ norm, and at an intermediate rate (approximately 1.5) in the L2 norm. Convergence rates
for the deformation are somewhat less regular, with nearly second-order convergence rates being
observed in the L1 and L2 norms and between first- and second-order convergence rates observed in
the L∞ norm. The robustness of the convergence rate in the deformation can be improved by using
higher-order structural elements. Because the overall accuracy of the discretization is also limited
by the Eulerian discretization and the form of the regularized kernel function, however, the use of
higher-order elements does not in itself increase the overall order of accuracy of the method.
Notice that in both the static and dynamic test cases, virtually identical errors are attained for all of
the values of Mfac considered. This indicates that for these tests, the method is able to use relatively
coarse structural meshes without appreciable loss in accuracy. In particular, these results suggest
that the scheme does not allow leaks at fluid-structure interfaces, even for Lagrangian meshes that
are quite coarse compared to the Eulerian grid.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

15

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD
L1
10

L2

-4

10

L∞

-4

u

10 -4

10

-6

10

-6

10 -6

10

-2

10

-2

10

-2

10 0
10

-2

p

10

-2
M

FAC

=4

M FAC = 2

10

M

-4

10
10

-2

10

∆x
1

FAC

=1

1st-order
2nd-order

-2

∆x
2

-2

10

-2

∆x

∞

Figure 5. Errors in u and p in L , L , and L norms for the static (γ = 0) version of the idealized anisotropic
shell model of sec. 5.1.1. Reference lines with slopes of -1 and -2 are also shown. The velocity converges
at second-order accuracy in all norms, whereas the pressure converges at second-order in the L1 norm, at
first-order in the L∞ norm, and at order 1.5 in the L2 norm.

5.1.2. Orthotropic shell The second case that we consider uses a neo-Hookean material model,
W e (F) =

µe
I1 (C),
2w

(62)

µe
F.
w

(63)

with C = FT F and I1 (C) = tr(C), so that
Pe =

Because of the form of the mapping from curvilinear coordinates to initial coordinates, this material
behaves as an orthotropic fiber-reinforced solid rather than as an isotropic material. One way of
viewing the elastic response is that the body is composed of two continuous families of fibers.
The first family of fibers wraps the elliptical shell circumferentially, and the second family is
composed of radial fibers that are orthogonal to the circumferential fibers. Because one family of
fibers terminates along the fluid-structure interfaces, there are singular force layers along ∂χ(U, t)
that must be balanced by discontinuities in the pressure and viscous stress. Therefore, in this case
the discretized unified and partitioned formulations yield different results.
R
Setting γ = 0, so that the structure is in equilibrium, and requiring Ω p(x, t) dx = 0, it can be
shown [37] that


1
e 1

r ≤ R,
p0 + µe R − R+w

1
R
(64)
p(x, t) = p0 + µw R
R < r ≤ R + w,
(R + w − r) + R+w

p
R + w < r,
0


e
(R+w)3
2
with r = kx − (0.5, 0.5)k and p0 = πµ
. We set ρ = 1, µ = 1, and µe = 1,
3w 3wR + R −
R
and we consider the time interval 0 ≤ t ≤ 3. Fig. 8 summarizes the error data at time t = 3
for N = 64, 128, 256, 512, and 1024, using Mfac = 1, 2, and 4, with ∆t = 0.25∆x. First-order
convergence rates are observed for u in all norms. First-order convergence rates are also observed for
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

16

B. E. GRIFFITH AND X. Y. LUO
L1
-2

10

-3

10

-4

u

10

L2

10

L∞

-2

10

10 -3
10
10

-2

-4

-2

10

-2

10

p

10

10

-2

10
10

-2

10

-2

-2

10 -2

10

χ

0

-2

10

10

-2

-2

-4

10

M

10 -3

-4

M

10

-2

10

∆x

∆x
1

2

-2

10

FAC

=4

M FAC = 2
FAC

=1

1st-order
2nd-order

-4

10

-2

∆x

∞

Figure 6. Errors in u, p, and χ in L , L , and L norms for the dynamic (γ = 0.15) version of the idealized
anisotropic shell model of sec. 5.1.1. Reference lines with slopes of -1 and -2 are also shown. The velocity
converges at second-order accuracy in all norms, whereas the pressure converges at second-order in the
L1 norm, at first-order in the L∞ norm, and at order 1.5 in the L2 norm. The displacement converges at
essentially second-order rates in the L1 and L2 norms, and at a slightly lower rate in the L∞ norm.

p in the L1 norm. Because p possesses discontinuities at fluid-structure interfaces for this problem,
however, the present method yields convergence rates of 0.5 in the L2 norm and does not converge
in the L∞ norm.
We also consider the case in which γ = 0.15, so that the initial configuration of the shell is not
in equilibrium. We set ρ = 1, µ = 0.01, and µe = 1, yielding a Reynolds number of approximately
100. We consider the time interval 0 ≤ t ≤ 1.25, which corresponds to approximately one damped
oscillation of the shell. Again, an exact solution is not available, and so convergence rates are
estimated using Richardson extrapolation [49]. Fig. 9 summarize the error data at time t = 0.75
for N = 64, 128, 256, and 512 and Mfac = 1 and 4, with ∆t = 0.25∆x. Essentially first-order
convergence rates are observed for u and χ in all norms, whereas p exhibits first-order convergence
in only the L1 norm.
For this problem, we find that the unified and partitioned formulations yield similar accuracy in
most cases for u and χ. By contrast, the partitioned formulation offers significantly better accuracy
for the pressure for relatively coarse Lagrangian meshes. This property appears also to result in
improvements in volume conservation; see sec. 5.2.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

17

A.

B.

C.

Figure 7. Similar to fig. 4, but here showing results obtained using the orthotropic shell model of sec. 5.1.2
for N = 128 and the partitioned (split) weak formulation over the time interval 0 ≤ t ≤ 1.25. As in fig. 4,
the coarse and fine structural meshes yield very similar kinematics.

5.2. Soft elastic disc in lid driven cavity
This section presents results from tests that use a soft elastic structure in a lid-driven cavity flow to
demonstrate the volume-conservation properties of our method.
In these computations, the physical domain is Ω = [0, 1] × [0, 1], and the immersed structure is a
disc of radius 0.2 that is initially centered about x = (0.6, 0.5). The boundary conditions imposed
along ∂Ω are u ≡ 0 on the left (x1 = 0), right (x1 = 1), and bottom (x2 = 0) boundaries of Ω, and
u ≡ (1, 0) on the top (x2 = 1) boundary of Ω. We use an isotropic neo-Hookean model,
Pe = µe F − p0 F−T ,

(65)

and we consider both p0 = 0 and p0 = µe . Because generally Pe N 6≡ 0, the solution has
discontinuities in the pressure and viscous stress at fluid-structure interfaces. In such cases, we
expect the IB method to yield no better than first-order convergence rates.‡ The flow induced by the
driven lid brings the structure nearly into contact with the moving upper boundary of the domain;
see fig. 10. This near contact is handled automatically by the IB formulation using a version of a
modified kernel function approach introduced by Griffith et al. [24] and enhanced by Kallemov et
al. [43]. No additional specialized methods are required by the present scheme to handle this case.
As in previous studies of this test case [33,66], we set µ = 0.01 and ρ = 1. We consider µe = 0.2,
which is a relatively “soft” case. The initial velocity is u ≡ 0, and the reference coordinates X are
‡ This flow also possesses well-known corner singularities that reduce the convergence rate of the incompressible flow
solver. Although it is possible to devise numerical schemes that accurately treat the corner singularities present in the
classical lid-driven cavity flow [65], we do not employ such a method in this work.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

18

B. E. GRIFFITH AND X. Y. LUO
L1

L2

L∞

u

M FAC = 4

10 -4

10

10

-3

10

-4

M FAC = 1
1st-order

-4

10 -2

10 -2

10 -2

∆x
0

p

10

10

0

10 -1

M

FAC

=4

M FAC = 1
0.5-order
1st-order

10

-2

10

∆x

-2

∆x
1

2

∞

Figure 8. Errors in u and p in L , L , and L norms for the static (γ = 0) version of the orthotropic shell
model of sec. 5.1.2. Errors for the unified formulation appear as solid lines, and errors for the partitioned
formulation appear as dashed lines. Reference lines with slope -1 are provided for the u error data. For p,
reference lines with slopes -1 and -0.5 are provided for the L1 and L2 norm data, respectively. Because this
test includes discontinuities in the pressure at fluid-structure interfaces, the present method does not yield
convergence in p in the L∞ norm. The partitioned formulation generally yields improved accuracy compared
to the unified formulation, especially for the pressure for relatively coarse structural mesh spacings.

taken to be the initial coordinates, so that χ(X, 0) ≡ X. The physical domain is discretized using
an N × N Cartesian grid. The Lagrangian coordinate domain is discretized using an unstructured
mesh of quadratic triangular (P 2 ) elements constructed so that the elements are approximately a
factor of Mfac coarser than the background Eulerian grid, so that in the reference configuration, the
nodes of the Lagrangian mesh are physically separated by a distance of approximately Mfac ∆x.
The time step size is ∆t = 0.125∆x. We consider the time interval 0 ≤ t ≤ 10, during which
the disc is subjected to slightly more than one rotation within the cavity. The structure becomes
entrained in the shearing flow along the cavity lid from t ≈ 4 until t ≈ 6, and during this time
is subjected to very large deformations. Fig. 11 shows the percent change in disc volume for
different values of N and Mfac . With p0 = 0, the maximum volume change yielded by the unified
formulation is approximately 2.3% for N = 64 and Mfac = 4, 0.4% for N = 64 and Mfac = 2, and
0.2% for N = 64 and Mfac = 1. The split formulation yields substantially improved accuracy: the
maximum volume change is approximately 0.12% for N = 64 and Mfac = 4, 0.2% for N = 64
and Mfac = 2, and 0.15% for N = 64 and Mfac = 1. Setting p0 = µe improves the accuracy of the
unified formulation substantially, especially at coarser relative mesh spacings, whereas it results in
slightly poorer volume conservation for the split formulation. With p0 = µe , the maximum volume
change is less than 0.4% in all cases considered. At smaller values of Mfac , there is essentially no
difference in the volume changes produced by the unified and partitioned formulations. In all cases,
the maximum volume change converges to zero at a first-order rate.
In general, the partitioned formulation appears to yield superior volume conservation for
this problem, especially at relatively coarse Lagrangian mesh spacings. Moreover, the volumeconservation properties of the partitioned scheme seem to be relatively insensitive to the relative
coarseness of the Lagrangian mesh. Using either formulation, volume errors converge to zero at
essentially a first-order rate. These results compare very favorably to those obtained by the IFE
method, which, even for relatively fine Lagrangian meshes, can yield volume losses of up to 20%
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD
L1

L2

19

L∞
M FAC = 4

u

M

10

10

-1

10

-2

FAC

=1

1st-order

10 -2

-2

10

-2

10

-2

10

-2

10

-2

M FAC = 4

10 0

M

=1

p

10

FAC

0.5-order
1st-order

0

10

10 -1

-2

10 -2

10 -2

10

-2

M FAC = 4
M

FAC

=1

χ

1st-order

10

-4

10

10

-4

-3

10 -4
10

-2

10

∆x

∆x
1

2

-2

∆x

∞

Figure 9. Errors in u, p, and χ in L , L , and L norms for the dynamic (γ = 0.15) version of the
orthotropic shell model of sec. 5.1.2. Errors for the unified formulation appear as solid lines, and errors
for the partitioned formulation appear as dashed lines. Reference lines with slope -1 are provided for the
u and χ error data. For p, reference lines with slopes -1 and -0.5 are provided for the L1 and L2 norm
data, respectively. The unified formulation generally yields modest improvements in the accuracy for u and
χ, whereas the partitioned formulation generally yields improved accuracy for p, especially for relatively
coarse structural mesh spacings.

Figure 10. A soft neo-Hookean disc in a lid-driven cavity flow using the partitioned (split) weak formulation
with N = 128 and Mfac = 4 over the time interval 3 ≤ t ≤ 7.

when applied to the same test without using a volume-conservation algorithm. A modification of the
IFE method to improve its volume conservation still yields volume losses of approximately 2.5%
for this test [33].
5.3. Flow past a cylinder
This section presents results from tests using the widely used benchmark of viscous flow past
a stationary circular cylinder. In these computations, the physical domain is Ω = [−15, 45] ×
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

20

B. E. GRIFFITH AND X. Y. LUO

A.

p0 = 0

B.

p 0 = µe

Figure 11. Percent change in structure volume for the soft elastic disc benchmark of sec. 5.2 as a function
of time using eq. (65) with A. p0 = 0 and B. p0 = µe . Results are shown for Cartesian grids of size N = 64,
128, and 256 with relative Lagrangian mesh spacings Mfac = 4, 2, and 1. The amount of spurious volume
change converges to zero at first order. At coarser relative mesh spacings, the partitioned (split) formulation
yields substantially better volume (area in two spatial dimensions) conservation than the unified (unsplit)
formulation. Volume conservation is substantially improved for the unified formulation by setting p0 = µe .
The volume change is less than 0.4% for all cases considered except for the unified formulation with p0 = 0
and Mfac = 4.

[−30, 30], and the immersed structure is a thin circular interface of radius 0.5 centered about x =
(0, 0). This domain size and cylinder placement corresponds to ‘Case C’ considered by Taira and
Colonius [67]. Along the inflow boundary (x1 = −15), we set a uniform inflow velocity, u ≡ (1, 0).
Along the outflow boundary (x1 = 45), we set the normal traction and tangential velocity to zero,
whereas along the top (x2 = 30) and bottom (x2 = −30) boundaries, we set the normal velocity and
tangential traction to zero. The boundary condition treatment is the same as that described by Griffith
[51]. We set ρ = 1 and µ = 0.005. Using the inflow velocity as the characteristic velocity u∞ and
the cylinder diameter d as the characteristic length, the Reynolds number is Re = ρuµ∞ d = 200.
The computational domain is discretized using an adaptively refined Cartesian grid with six nested
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

A.

21

B.

Figure 12. Vortices shed from a stationary circular cylinder at Re = 200. This computation uses a six-level
locally refined grid with a refinement ratio of two between levels and an N × N coarse grid with N = 128.
Regions of local mesh refinement are indicated by rectangular boxes, with lighter grey boxes indicating
coarser levels of refinement and black boxes indicating the finest regions of the locally refined grid. A. The
full computational domain and B. a magnified view of the flow field near the cylinder.

CL
CD
St
—
—
0.190
±0.69 1.34 ± 0.044 0.197
±0.69 1.31 ± 0.049 0.192
±0.68 1.34 ± 0.047 0.195
—
—
0.190
—
—
0.197
±0.70 1.36 ± 0.046 0.195
Table I. Comparison of experimental and computational values of lift (CL ) and drag (CD ) coefficients and
Strouhal numbers (St) for flow past a circular cylinder at Re = 200.

Lai and Peskin [68]
Linnick and Fasel [69]
Liu et al. [70]
Taira and Colonius [67] (Case C)
Roshko (experimental) [70]
Williamson (experimental) [70]
present (four-point kernel, N = 128, Mfac = 2)

grid levels and a refinement ratio of two between levels. The Cartesian grid spacing on the finest
grid level is ∆xfinest = 2−5 ∆xcoarsest , with ∆xcoarsest = 60
N . The cylinder is discretized using a mesh
of one-dimensional quadratic (P 2 ) elements with a node spacing of approximately Mfac ∆xfinest .
The time step size is ∆t = 0.1∆xfinest , yielding a maximum CFL number of approximately 0.1–
0.2. Rigidity constraints are approximately imposed using tether forces as in eq. (29). For each grid
spacing considered, we use approximately the largest stable values of the penalty parameters κ and η
as permitted by the time step size. These values are empirically determined by a simple optimization
procedure. Representative results are shown in fig. 12.
Table I compares results obtained using the present method with the standard four-point kernel
[3], N = 128, and Mfac = 2 with corresponding experimental and computational results from prior
studies. We compute the lift coefficient, CL = Fy /(ρu2∞ d/2), drag coefficient, CD = Fx /(ρu2∞ d/2),
and Strouhal number, St = fs d/u∞ , in which F = (Fx , Fy ) is the net force on the cylinder and fs is
the shedding frequency. The present results are seen to be in excellent quantitative agreement with
these earlier results.
Fig. 13 shows the lift and drag coefficients as functions of time for N = 32, 64, and 128 and for
Mfac = 1, 2, and 4 using the four-point regularized kernel function [3]. Similar results are shown in
fig. 14 for the similar three-point kernel [71], whereas fig. 15 shows results for a recently developed
six-point kernel with higher continuity order [72]. By construction, the structural meshes obtained
for a fixed value of N are nested versions of each other (i.e. they can be viewed as obtained via
Lagrangian mesh refinement). Under simultaneous Lagrangian and Eulerian grid refinement, the
scheme converges to the same dynamics for all values of Mfac . For the four- and six-point kernels,
however, the best accuracy is obtained for larger values of Mfac . In particular, by using relatively
coarser structural meshes, the scheme yields more accurate forces (CL and CD ) and vortex shedding
dynamics (characterized by, e.g., St). In fact, for the coarsest Eulerian discretization considered
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

22

B. E. GRIFFITH AND X. Y. LUO
CL

N = 32

1

CD

2

0

1

M FAC = 4
M FAC = 2
M FAC = 1

N = 64

-1

10

20

30

40

50

1

2

0

1

-1

N = 128

0

0

0

10

20

30

40

50

0

1

2

0

1

-1

0

10

20

30

40

50

0

0

10

20

30

40

50

0

10

20

30

40

50

0

10

20

30

40

50

time

time

Figure 13. Lift (CL ) and drag (CD ) coefficients for flow past a stationary cylinder at Re = 200. The
computational domain Ω is discretized using a six-level locally refined grid with a refinement ratio of
two between levels and an N × N coarse grid. The spacing between the nodes of the immersed structure
is approximately Mfac ∆xfinest . Under simultaneous Lagrangian and Eulerian grid refinement, the scheme
converges to the same dynamics for all values of Mfac . Notice, however, that the best accuracy is obtained
for larger values of Mfac . Using a relatively finer structural mesh spacing results in larger lift and drag
coefficients and lower vortex shedding frequencies. Thus, smaller values of Mfac increase the effective
numerical size of the cylinder.
CL

N = 32

1

CD

2

0

1

M FAC = 4
M FAC = 2
M FAC = 1

N = 64

-1

10

20

30

40

50

1

2

0

1

-1

N = 128

0

0

0

10

20

30

40

50

0

1

2

0

1

-1

0

10

20

30

time

40

50

0

0

10

20

30

40

50

0

10

20

30

40

50

0

10

20

30

40

50

time

Figure 14. Similar to figs. 13 and 15, but here using a three-point kernel function [71]. Unlike the case of the
four- and six-point kernels, comparable accuracy is obtained for all values of Mfac considered.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

23

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD
CL

N = 32

1

CD

2

0

1

M FAC = 4
M FAC = 2
M FAC = 1

N = 64

-1

10

20

30

40

50

1

2

0

1

-1

N = 128

0

0

0

10

20

30

40

50

0

1

2

0

1

-1

0

10

20

30

40

50

0

0

10

20

30

40

50

0

10

20

30

40

50

0

10

20

30

40

50

time

time

Figure 15. Similar to figs. 13 and 14, but here using a six-point kernel function with higher continuity
order [72]. As with the four-point kernel, the best accuracy is obtained for larger values of Mfac . In this
case, however, using relatively fine structural mesh spacings yields erratic results on coarser Eulerian grids
(e.g. Mfac = 1 for N = 32).

(N = 32), the six-point kernel yields erratic results for the finest structural discretization (Mfac = 1)
that do not occur for the coarser structural discretizations (Mfac = 2 and 4). In contrast, with the
three-point kernel, comparable results are obtained for all values of Mfac considered here. Although
not shown here, results obtained using a two-point piecewise-linear kernel are similar to those
obtained using the three-point kernel.
A striking feature of these results is that the use of Lagrangian mesh refinement alone generally
lowers the accuracy of the computation for a fixed Eulerian grid. A complete theoretical explanation
for this behavior is lacking at present.
5.4. Idealized model of left ventricular mechanics
To demonstrate the capabilities of the methodology to treat more complex geometries and structural
models, this section briefly presents results from tests using an idealized model of the left ventricle
of the heart, as used in a previous benchmarking study of cardiac mechanics codes by Land et
al. [45].
For these tests, the undeformed geometry of the idealized left ventricle is described as a truncated
ellipsoid. Using parametric coordinates (t, u, v), the ventricular geometry is defined by
X(t, u, v) = ((7 + 3t) sin u cos v, (7 + 3t) sin u sin v, (17 + 3t) cos u) ,

(66)

with length measured in millimeters. The endocardial surface is defined by t = 0, u ∈
5
], and v ∈ [−π, π], and the epicardial surface is defined by t = 1, u ∈
[−π, − arccos 17
5
[−π, − arccos 20 ], and v ∈ [−π, π]. The model ventricle is truncated at the base, which is taken
to correspond to z = 5 mm. See fig. 16(A). The passive elasticity of the heart wall is described
using a transversely isotropic hyperelastic constitutive law by Guccione et al. [73], which is defined
with respect to a local fiber-aligned coordinate system via
We =
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls


C Q
e −1 ,
2

(67)
Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

24

B. E. GRIFFITH AND X. Y. LUO

A.

B.

C.

D.

E.

F.

Figure 16. Inflation and active contraction of an idealized model of the left ventricle of the heart. A. The
reference configuration of the left ventricle model is a truncated ellipsoid. B. Passive inflation of an isotropic
version of the left ventricular model. C. Active contraction of a transversely isotropic left ventricular model
that includes transmural fiber rotation. Notice the torsion induced by the active contractile forces. D–F. Slices
of the configurations shown in panels A–C.

with


2
2
2
2
2
2
2
2
2
,
+ E31
+ bfs E12
+ E21
+ E13
+ E23
+ E32
Q = bf E11
+ bt E22
+ E33

(68)

in which E = 12 (FT F − I) = (Eij ) is the Green-Lagrange strain tensor and C , bf , bt , bfs are material
parameters. The structure is discretized using trilinear (Q1 ) hexahedral elements. The computational
domain Ω is a 4 cm × 4 cm × 4 cm region. Except where noted, Ω is discretized using a uniform
64 × 64 × 64 Cartesian grid. The model is run to steady state under steady loading conditions,
corresponding to a quasi-static analysis, as in earlier work [74].
5.4.1. Passive inflation of an isotropic left ventricle model Figs. 16(B,E) and 17(A) show results
of passively inflating an isotropic version of the left ventricle model. We use C = 10 kPa and
bf = bt = bfs = 1 along with a uniform pressure load of 10 kPa. This corresponds to ‘Problem
2’ from Land et al. [45]. Notice that our formulation permits very large structural deformations.
The computed longitudinal, circumferential, and radial strains in the inflated configuration are in
excellent agreement (within 1%) with values obtained by other structural mechanics codes using
the same model [45]. This test problem was judged to be relatively easy, and all of the methods
considered in the benchmarking study yielded virtually identical results.
5.4.2. Active contraction of a fiber-reinforced left ventricle model Figs. 16(C,F) and 17(B) show
the results from a simulation of active contraction using this model. In this case, we use a fiberreinforced material model, in which the fiber direction in the reference configuration Ef is given
by


∂X
∂X
∂X
∂X
Ef = sin α
+ cos α
, α = 90 − 180t,
(69)
∂u
∂u
∂v
∂v

yielding a 180◦ transmural fiber rotation with a linear relationship between fiber angle and tissue
depth. In addition to passive elasticity, the test also includes active contractile forces, which are

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

A.

circ

80

long

endo

trans

-25

60

25

-30

60

-35
55

40

-40
-45

20

2

4

6

8

10

50

2

4

6

-50

8

50

epi
mid

6

8

10

2

4

6

8

10

4

6

8

10

-25
40

-30

30
20

30

10

20

-35
-40

2

4

6

8

10

2

4

6

-45

8

60

55

-25

50

50

-30

40

45

-35

30

40

-40

20

35

2

4

6

8

10

2

position

4

6

-45

8

2

position

circ

0

position

long

-16

40

-18
-10

trans

50

-17

-5

endo

4

-20

40

B.

2

30

-19
20

-15

2

4

6

8

10

15

epi

2

4

6

8

-6

20

2

4

6

8

10

2

4

6

8

10

4

6

8

10

22
20

-8

10

18
-10

5

16
-12

0
-5

mid

-20

2

4

6

8

10

-14

5

-6

0

-8

-5

-10

-10

-12

-15

-14

-20

-16

14
2

4

6

12

8

35
30
25

2

4

6

8

10

20

2

position

4

6

15

8

2

position

position

Figure 17. Circumferential, longitudinal, and transverse strains along the endocardium, epicardium, and
midwall in the idealized left ventricular model for (A) passive inflation and (B) active contraction. Strains
are plotted at selected points running from apex to base. See Land et al. [45] for additional details on the
positions of the sample points. The computed strains are in excellent agreement (within 1%) with consensus
results obtained by other structural mechanics codes using the same model [45].

accounted for in the first Piola-Kirchhoff stress via
Pe =

∂W
+ T a F Ef ⊗ Ef ,
∂F

(70)

with Ta corresponding to the active stress acting in the fiber direction. We use C = 2 kPa, bf = 8,
bt = 2, bfs = 4, and Ta = 60 kPa. A uniform pressure load of 15 kPa is also applied along the
endocardial surface. This corresponds to ‘Problem 3’ from Land et al. [45]. No special treatment
Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

26

B. E. GRIFFITH AND X. Y. LUO

A.

B.

Figure 18. Difference in the displacement of the actively contracting fiber-reinforced left ventricle model
using different N × N × N Cartesian grids. A. Distribution of difference in the structural displacement
obtained using N = 48 and N = 64. B. Distribution of difference in the structural displacement obtained
using N = 64 and N = 96.

Figure 19. Fiber strain distribution in the actively contracting fiber-reinforced left ventricle model obtained
using a 64 × 64 × 64 Cartesian grid.

is used for the fiber singularity at u = −π . To demonstrate convergence, the differences in the
structural displacement obtained using uniform N × N × N grids for N = 48, 64, and 96 are shown
in fig. 18. The maximum difference in the displacement between N = 48 and N = 64 is 2.5%; this
difference is 1.8% between N = 64 and N = 96, indicating approximately first-order convergence.
This test problem was judged to be relatively difficult when compared to Problem 2 [45], but as in the
case of passive deformations, the actively contracting model is in excellent agreement (within 1%)
with grid-converged consensus results obtained by other structural codes using the same idealized
left ventricular model [45]. The distribution of fiber strain within the model is shown in fig. 19.
It is interesting to see that although the ventricular wall is mostly compressed, as indicated by
the negative values of the fiber strain, there are parts of the wall that are very slightly stretched,
primarily in the midwall. This is because the fibers rotate across the wall, and the middle layer
experiences circumferential stretch even though the whole left ventricle is substantially compressed
in the longitudinal direction. The comparatively small stretches near the apex reflect the different
fiber distribution in that region.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

27

6. DISCUSSION AND CONCLUSIONS
This paper describes a version of the IB method for fluid-structure interaction that uses either a
structured or unstructured FE discretization of the immersed structure while retaining a Cartesian
grid finite difference scheme for the Eulerian equations. This method has already proved useful
in a variety of biological and biomedical applications [46–48, 74–81], but its performance has not
previously been examined using standard benchmark cases such as those that are the focus of this
work. This paper also demonstrates the suitability of this method to simulate cardiac mechanics
using an idealized model of the left ventricle of the heart [45].
A feature of this method is that it uses standard discretization technology for both the Lagrangian
and Eulerian equations. In practice, it should be feasible to use this approach to couple existing
structural analysis and incompressible flow codes by passing forces from the structural code to the
fluid solver, and by passing velocities from the fluid solver back to the structural code. Although we
focus on cases in which the immersed structure is hyperelastic, our numerical scheme does not rely
on the availability of a hyperelastic energy functional, and the present approach is also demonstrated
to treat immersed rigid bodies and structures with active contractile stresses.
A key contribution of this paper is that it introduces an approach to Lagrangian-Eulerian
interaction that overcomes a longstanding limitation of the IB method, namely the requirement
that the Lagrangian mesh must be relatively fine compared to the background Eulerian grid to avoid
leaks at fluid-structure interfaces. This restriction results in structural meshes that are excessively
dense, potentially leading to reduced efficiency and even reduced accuracy. Numerical examples
demonstrate that our scheme permits the use of Lagrangian meshes that are at least four times
as coarse as the background Eulerian grid without leak, and we expect that there are cases in
which even coarser structural meshes would yield good accuracy. Even with very coarse structural
discretizations, the present methodology can yield improved volume conservation when compared
to the IFE method [33]. Of course, fine Lagrangian meshes are still needed in practice to resolve
fine-scale geometrical features, and an approach such as Lagrangian adaptive mesh refinement may
be crucial for accurately treating extremely large structural deformations. Nonetheless, the present
approach enables the effective use of spatial discretizations that are tailored to the requirements of
the structural model rather than dictated by the background Eulerian discretization.
Our method is based on a formulation of the continuous IB equations introduced by Boffi et
al. [37]. We consider two different statements of these equations that each use a weak formulation
of the Lagrangian equations of nonlinear elasticity. One of these formulations treats the internal
and boundary stresses within the immersed elastic structure using a single, unified, volumetric
elastic force density. This form of the equations of motion is essentially that used in the IFE
method [30–34] as well as other similar approaches [35–39]. The other formulation considered
in this work partitions these stresses into a internal elastic force density supported on the interior
of the structure, and a transmission elastic force density supported on fluid-structure interfaces;
this approach does not appear to have been used previously to develop a numerical scheme. Both
formulations are demonstrated to yield similar convergence rates, but the partitioned formulation is
seen to yield higher accuracy for Lagrangian meshes that are relatively coarse in comparison to the
background Eulerian grid, especially in terms of volume conservation.
A limitation of the partitioned formulation is that it does not satisfy a discrete power identity
that implies that energy is conserved during Lagrangian-Eulerian interaction. Such a power identity
is satisfied by the unified formulation, and may be necessary to obtain an unconditionally stable
implicit time-stepping scheme [82]. Developing a symmetric partitioned formulation will likely
require the introduction of additional boundary degrees of freedom, so that volumetric operators
(i.e. S and J ) couple the volumetric structural variables to the background grid, and corresponding
surface operators (e.g. S ∂U and J ∂U ) couple the surface degrees of freedom to the Eulerian
variables. Lagrange multipliers or penalty methods could be used to ensure that the surface
discretization moves with the immersed body.
For immersed rigid bodies, tests suggest that relatively coarse structural discretizations can yield
superior accuracy when compared to relatively fine discretizations. For the standard test case of

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

28

B. E. GRIFFITH AND X. Y. LUO

viscous flow past a circular cylinder, we demonstrate that relatively fine structural discretizations
can produce spurious drag, suggesting that the effective numerical size of the immersed body is
determined in part by the spacing of the control points. Although not studied here in detail, there
appears to be a relatively sharp transition, and the results obtained by the method for structural
meshes that are coarser than a minimum spacing appear to be relatively insensitive to the choice of
Lagrangian meshwidth. The threshold spacing appears to depend on the choice of regularized kernel
function. For instance, when using a three-point kernel function, the method yields good results for
relative mesh spacing values of Mfac = 1, whereas this same relative structural mesh spacing yields
poor results with four- and six-point kernels.
In a constrained formulation, using a mesh of control points that is too fine with respect to the
background grid will result in an ill-conditioned system of equations, and for a sufficiently fine
mesh of control points, the constrained formulation becomes singular. It is interesting to note that in
our formulation, we obtain high accuracy even with extremely dense meshes of interaction points
so long as the control points are sufficiently far apart from each other. Projecting the standard IB
velocity field UIB onto the FE shape functions filters velocity fluctuations at length scales that are
smaller than structural mesh width. For relatively coarse structural meshes, this additional filtering
appears to improve the accuracy of the IB method, dramatically reducing errors in the lift and drag
forces. Such errors have been previously described to be a fundamental aspect of the IB method, but
we believe that our results show that such numerical artifacts are not intrinsic to this methodology,
but rather result from the use of overly dense structural meshes. The present work provides a
framework for coupling relatively course structural models without leak at fluid-structure interfaces.
In closing, we remark that the partitioned formulation developed in this work could be useful in
constructing higher-order versions of the IB method. For instance, this formulation is well suited for
developing a hybrid approach in which the IB method is used to treat the volumetric internal force
density, and another method is used to treat the singular transmission force density. Because the
transmission force density of the partitioned formulation is defined on a closed surface, it is feasible
to treat it with higher-order accuracy by a version of the immersed interface method [83–89].
Such an extension of this method could yield a fully second-order accurate generalization of the
IB method for cases, like those considered in this work, in which the immersed structure is of
codimension zero with respect to the fluid.

ACKNOWLEDGEMENTS
This work was sponsored in part by an award from the Royal Society of Edinburgh.
B.E.G. acknowledges research support from the National Science Foundation (NSF awards ACI
1047734/1460334, ACI 1450327, and DMS 1016554/1460368) and the National Institutes of
Health (NIH award HL117063). He also gratefully acknowledges discussions of this work with
Charles Peskin and David McQueen. X.Y.L. acknowledges research support from the Engineering
and Physical Sciences Research Council (UK EPSRC awards EP/I029990 and EP/N014642/1)
and from the Leverhulme Trust (RF-2015-510). We also acknowledge Hao Gao (U. Glasgow)
and Viatcheslav Gurev (IBM T. J. Watson Research Center) for their assistance with the cardiac
mechanics benchmark.

A. TEMPORAL DISCRETIZATION
A.1. Basic time-stepping scheme
1

Let χn , un , and pn− 2 denote the approximations to the values of χ and u at time tn , and to the value
1
of p at time tn− 2 , respectively. We advance the solution values forward in time by the increment ∆t
as follows. First, we determine a preliminary approximation to the deformed structure configuration

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

29

1

at time tn+ 2 via

1

χn+ 2 − χn
= J (χn ) un .
∆t/2

Then, we solve


 n+1
 n+1
1
1
u
− un
u
+ un
n+ 12
n+ 12
2
+A
= −∇h p
+ f n+ 2 + tn+ 2 ,
ρ
+ µ∇h
∆t
2
∇h · un+1 = 0,

 


1
1
f
= S χn+ 2 F χn+ 2 ,

 

1
1
1
tn+ 2 = S ∂U χn+ 2 T χn+ 2 ,
  un+1 + un 

χn+1 − χn
1
= J χn+ 2
,
∆t
2
n+ 21

(71)

(72)
(73)
(74)
(75)
(76)

1

1

for χn+1 , un+1 , and pn+ 2 , in which An+ 2 = 23 un · ∇h un − 21 un−1 · ∇h un−1 is computed via a
PPM-type approximation to the nonlinear advection term [50, 51]. The time-stepping scheme for
1
the unified formulation is analogous. Notice that solving eqs. (72)–(76) for χn+1 , un+1 , and pn+ 2
requires the solution of a Crank-Nicolson-type discretization of the time-dependent incompressible
Stokes equations. We solve this system of equations via the flexible GMRES (FGMRES) algorithm
1
1
[90], using un and pn− 2 as initial approximations to un+1 and pn+ 2 , and using a pressure-free
projection method with inexact multigrid subdomain solvers as a preconditioner [51].
A.2. Initial time step
Because time step-lagged values of u and p are used by the time-stepping scheme of sec. A.1,
we cannot use that scheme for the initial time step. Instead, we use a two-step predictor-corrector
method: First, we solve
 n+1


 n+1
1
ũ
+ un
ũ
− un
+ An = −∇h p̃n+ 2 + µ∇2h
+ f n + tn ,
(77)
ρ
∆t
2
∇h · ũn+1 = 0,
f n = S (χn ) F (χn ) ,

tn = S ∂U (χn ) T (χn ) ,
χ̃

n+1

(78)
(79)

(80)

n

−χ
= J (χn ) un ,
∆t

(81)

1

for χ̃n+1 , ũn+1 , and p̃n+ 2 , with An = un · ∇h un . Because we do not have an initial value for
n+1
n
1
1
the pressure, we use p = 0 as an initial guess for pn+ 2 . Second, we set χn+ 2 = χ̃ 2+χ and
1
1
1
1
solve eqs. (72)–(76) forχn+1 , un+1 , and pn+ 2 , except that we use An+ 2 = un+ 2 · ∇h un+ 2 with
1
n+1
1
un+ 2 = 2 ũ
+ un .
REFERENCES
Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972; 10(2):252–271.
Peskin CS. Numerical analysis of blood flow in the heart. J Comput Phys 1977; 25(3):220–252.
Peskin CS. The immersed boundary method. Acta Numer 2002; 11:479–517.
Miller LA, Peskin CS. When vortices stick: an aerodynamic transition in tiny insect flight. J Exp Biol 2004;
207(17):3073–3088.
5. Miller LA, Peskin CS. A computational fluid dynamics of ‘clap and fling’ in the smallest insects. J Exp Biol 2005;
208(2):195–212.
6. Fogelson AL, Guy RD. Immersed-boundary-type models of intravascular platelet aggregation. Comput Meth Appl
Mech Eng 2008; 197(25–28):2087–2104.

1.
2.
3.
4.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

30

B. E. GRIFFITH AND X. Y. LUO

7. Yang X, Dillon RH, Fauci LJ. An integrative computational model of multiciliary beating. Bull Math Biol 2008;
70(4):1192–1215.
8. Hsu CY, Dillon R. A 3D motile rod-shaped monotrichous bacterial model. Bull Math Biol 2009; 71(5):1228–1263.
9. Miller LA, Peskin CS. Flexible clap and fling in tiny insect flight. J Exp Biol 2009; 212(19):3076–3090.
10. Tytell ED, Hsu CY, Williams TL, Cohen AH, Fauci LJ. Interactions between internal forces, body stiffness, and fluid
environment in a neuromechanical model of lamprey swimming. Proc Natl Acad Sci U S A 2010; 107(46):19 832–
19 837.
11. Skorczewski T, Griffith BE, Fogelson AL. Multi-bond models for platelet adhesion and cohesion. Biological Fluid
Dynamics: Modeling, Computation, and Applications, Olson SD, Layton AT (eds.). Contemporary Mathematics,
American Mathematical Society: Providence, RI, USA, 2014; 149–172.
12. Hoover A, Miller L. A numerical study of the benefits of driving jellyfish bells at their natural frequency. J Theor
Biol 2015; 374:13–25.
13. Jones SK, Laurenza R, Hedrick TL, Griffith BE, Miller LA. Lift vs. drag based mechanisms for vertical force
production in the smallest flying insects. J Theor Biol 2015; 384:105–120.
14. Kou W, Bhalla APS, Griffith BE, Pandolfino JE, Kahrilas PJ, Patankar NA. A fully resolved active musculomechanical model for esophageal transport. J Comput Phys 2015; 298:446–465.
15. Kou W, Pandolfino JE, Kahrilas PJ, Patankar NA. Simulation studies of circular muscle contraction, longitudinal
muscle shortening, and their coordination in esophageal transport. Am J Physiol Gastrointest Liver Physiol 2015;
309(4):G238–G247.
16. Tytell ED, Leftwich MC, Hsu CY, Griffith BE, Cohen AH, Smits AJ, Hamlet C, Fauci LJ. The role of body stiffness
in undulatory swimming: Insights from robotic and computational models. Phys Rev Fluids 2016; 1:073 202 (17
pages).
17. Peskin CS, McQueen DM. Fluid dynamics of the heart and its valves. Case Studies in Mathematical Modeling:
Ecology, Physiology, and Cell Biology, Othmer HG, Adler FR, Lewis MA, Dallon JC (eds.). Prentice-Hall:
Englewood Cliffs, NJ, USA, 1996; 309–337.
18. 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.
19. Lemmon JD, Yoganathan AP. Three-dimensional computational model of left heart diastolic function with fluidstructure interaction. J Biomech Eng 2000; 122(2):109–117.
20. Lemmon JD, Yoganathan AP. Computational modeling of left heart diastolic function: Examination of ventricular
dysfunction. J Biomech Eng 2000; 122(4):297–303.
21. 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.
22. McQueen DM, Peskin CS. Heart simulation by an immersed boundary method with formal second-order accuracy
and reduced numerical viscosity. Mechanics for a New Millennium, Proceedings of the 20th International
Conference on Theoretical and Applied Mechanics (ICTAM 2000), Aref H, Phillips JW (eds.), Kluwer Academic
Publishers, 2001.
23. 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.
24. Griffith BE, Luo XY, 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(1):137–177.
25. Griffith BE, Hornung RD, McQueen DM, Peskin CS. Parallel and adaptive simulation of cardiac fluid dynamics.
Advanced Computational Infrastructures for Parallel and Distributed Adaptive Applications, Parashar M, Li X
(eds.). John Wiley and Sons: Hoboken, NJ, USA, 2009; 105–130.
26. Griffith BE. Immersed boundary model of aortic heart valve dynamics with physiological driving and loading
conditions. Int J Numer Meth Biomed Eng 2012; 28(3):317–345.
27. 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.
28. 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.
29. Hammer PE, Sacks MS, del Nido PJ, Howe RD. Mass-spring model for simulation of heart valve tissue mechanical
behavior. Ann Biomed Eng 2011; 39(6):1668–1679.
30. Zhang L, Gerstenberger A, Wang X, Liu WK. Immersed finite element method. Comput Meth Appl Mech Eng
2004; 193(21–22):2051–2067.
31. Liu WK, Liu Y, Farrell D, Zhang L, Wang XS, Fukui Y, Patankar N, Zhang Y, Bajaj C, Lee J, et al.. Immersed finite
element method and its applications to biological systems. Comput Meth Appl Mech Eng 2006; 195(13–16):1722–
1749.
32. Zhang LT, Gay M. Immersed finite element method for fluid-structure interactions. J Fluid Struct 2007; 23(6):839–
857.
33. Wang X, Zhang LT. Interpolation functions in the immersed boundary and finite element methods. Comput Mech
2010; 45(4):321–334.
34. Wang X, Wang C, Zhang LT. Semi-implicit formulation of the immersed finite element method. Comput Mech
2012; 49(4):421–430.
35. Peskin CS. The Immersed Boundary Method. III. Energy Functions for the Representation of Immersed
Elastic Boundaries and Materials 2007. Handwritten lecture notes available from http://math.nyu.edu/
faculty/peskin/ib_lecture_notes/index.html.
36. Devendran D, Peksin CS. An energy-based immersed boundary method for incompressible viscoelasticity.
J Comput Phys 2012; 231(14):4613–4642.
37. Boffi D, Gastaldi L, Heltai L, Peskin CS. On the hyper-elastic formulation of the immersed boundary method.
Comput Meth Appl Mech Eng 2008; 197(25–28):2210–2231.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

FINITE DIFFERENCE/FINITE ELEMENT IMMERSED BOUNDARY METHOD

31

38. Heltai L, Costanzo F. Variational implementation of immersed finite element methods. Comput Meth Appl Mech
Eng 2012; 229–232:110–127.
39. Roy S, Heltai L, Costanzo F. Benchmarking the immersed finite element method for fluid-structure interaction
problems. Comput Math Appl 2015; 69(10):1167–1188.
40. Gil AJ, Carreño AA, Bonet J, Hassan O. The Immersed Structural Potential Method for haemodynamic
applications. J Comput Phys 2010; 229(22):8613–8641.
41. Gil AJ, Carreno AA, Bonet J, Hassan O. An enhanced Immersed Structural Potential Method for fluid-structure
interaction. J Comput Phys 2013; 250:178–205.
42. Hankar V, Wright GB, Kirby RM, Fogelson AL. Augmenting the immersed boundary method with Radial Basis
Functions (RBFs) for the modeling of platelets in hemodynamic flows. Int J Numer Meth Fluid 2015; 79(10):536–
557.
43. Kallemov B, Bhalla APS, Griffith BE, Donev A. An immersed boundary method for rigid bodies.
Comm Appl Math Comput Sci 2016; 11(1):79–141.
44. Usabiaga FB, Kallemov B, Delmotte B, Bhalla APS, Griffith BE, Donev A. Hydrodynamics of suspensions of
passive and active rigid particles: A rigid multiblob approach. Comm Appl Math Comput Sci 2016; 11(2):217–296.
45. Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, et al..
Verification of cardiac mechanics software: Benchmark problems and solutions for testing active and passive
material behaviour. Proc R Soc A 2015; 471(2184):20150 641 (20 pages).
46. Gao H, Carrick D, Berry C, Griffith BE, Luo XY. Dynamic finite-strain modelling of the human left ventricle in
health and disease using an immersed boundary-finite element method. IMA J Appl Math 2014; 79(5):978–1010.
47. Chen WW, Gao H, Luo XY, Hill NA. Study of cardiovascular function using a coupled left ventricle and systemic
circulation model. J Biomech 2016; 49(12):2445–2454.
48. Gao H, Feng L, Qi N, Berry C, Griffith BE, Luo XY. A coupled mitral valve-left ventricle model with fluid-structure
interaction ; Submitted.
49. Griffith BE, Peskin CS. On the order of accuracy of the immersed boundary method: Higher order convergence
rates for sufficiently smooth problems. J Comput Phys 2005; 208(1):75–105.
50. Griffith BE. On the volume conservation of the immersed boundary method. Commun Comput Phys 2012;
12(2):401–432.
51. 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.
52. Rider WJ, Greenough JA, Kamm JR. Accurate monotonicity- and extrema-preserving methods through adaptive
nonlinear hybridizations. J Comput Phys 2007; 225(2):1827–1848.
53. Colella P, Woodward PR. The piecewise parabolic method (PPM) for gas-dynamical simulations. J Comput Phys
1984; 54(1):174–201.
54. IBAMR: An adaptive and distributed-memory parallel implementation of the immersed boundary method. https:
//github.com/IBAMR/IBAMR.
55. SAMRAI: Structured Adaptive Mesh Refinement Application Infrastructure. http://www.llnl.gov/CASC/
SAMRAI.
56. Hornung RD, Kohn SR. Managing application complexity in the SAMRAI object-oriented framework.
Concurrency Comput Pract Ex 2002; 14(5):347–368.
57. Hornung RD, Wissink AM, Kohn SR. Managing complex data and geometry in parallel structured AMR
applications. Eng Comput 2006; 22(3–4):181–195.
58. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Eijkhout V, Gropp WD, Kaushik D,
et al.. PETSc Web page. http://www.mcs.anl.gov/petsc 2016. URL http://www.mcs.anl.gov/
petsc.
59. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Eijkhout V, Gropp WD, Kaushik
D, et al.. PETSc users manual. Technical Report ANL-95/11 - Revision 3.7, Argonne National Laboratory 2016.
URL http://www.mcs.anl.gov/petsc.
60. Balay S, Gropp WD, McInnes LC, Smith BF. Efficient management of parallelism in object oriented numerical
software libraries. Modern Software Tools in Scientific Computing, Arge E, Bruaset AM, Langtangen HP (eds.),
Birkhäuser Press, 1997; 163–202.
61. hypre: High performance preconditioners. http://www.llnl.gov/CASC/hypre.
62. Falgout RD, Yang UM. hypre: a library of high performance preconditioners. Computational Science - ICCS 2002
Part III, Lecture Notes in Computer Science, vol. 2331, Sloot PMA, Tan CJK, Dongarra JJ, Hoekstra AG (eds.).
Springer-Verlag, 2002; 632–641. Also available as LLNL Technical Report UCRL-JC-146175.
63. libMesh: A C++ Finite Element Library. http://libmesh.github.io.
64. Kirk B, Peterson JW, Stogner RH, Carey GF. libMesh: A C++ library for parallel adaptive mesh
refinement/coarsening simulations. Eng Comput 2006; 22(3–4):237–254.
65. Botella O, Peyret R. Benchmark spectral results on the lid-driven cavity flow. Comput Fluid 1998; 27(4):421–433.
66. Zhao H, Freund JB, Moser RD. A fixed-mesh method for incompressible flow-structure systems with finite solid
deformations. J Comput Phys 2008; 227(6):3114–3140.
67. Taira K, Colonius T. The immersed boundary method: A projection approach. J Comput Phys 2007; 225(2):2118–
2137.
68. Lai MC, Peskin CS. An immersed boundary method with formal second-order accuracy and reduced numerical
viscosity. J Comput Phys 2000; 160(2):705–719.
69. Linnick MN, Fasel HF. A high-order immersed interface method for simulating unsteady incompressible flows on
irregular domains. J Comput Phys 2005; 204(1):157–192.
70. Liu C, Zheng X, Sung CH. Preconditioned multigrid methods for unsteady incompressible flows. J Comput Phys
1998; 139(1):35–57.
71. Roma AM, Peskin CS, Berger MJ. An adaptive version of the immersed boundary method. J Comput Phys 1999;
153(2):509–534.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

32

B. E. GRIFFITH AND X. Y. LUO

72. Bao Y, Kaye J, Peksin CS. A Gaussian-like immersed boundary kernel with three continuous derivatives and
improved translational invariance. J Comput Phys 2016; 316:139–144.
73. Guccione JM, Costa KD, McCullock AD. Finite element stress analysis of left ventricular mechanics in the beating
dog heart. J Biomech 1995; 28(10):1167–1177.
74. 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.
75. Griffith BE, Flamini V, DeAnda A, Scotten L. Simulating the dynamics of an aortic valve prosthesis in a pulse
duplicator: Numerical methods and initial experience. J Med Dev 2013; 7(4):040 912 (2 pages).
76. Gao H, Ma XS, Qi N, Berry C, Griffith BE, Luo XY. A finite strain model of the human mitral valve with fluidstructure interaction. Int J Numer Meth Biomed Eng 2014; 30(12):1597–1613.
77. Flamini V, DeAnda A, Griffith BE. Immersed boundary-finite element model of fluid-structure interaction in the
aortic root. Theor Comput Fluid Dynam 2016; 30(1):139–164.
78. Jones SK, Yun YJ, Hedrick TL, Griffith BE, Miller LA. Bristles reduce the force required to ‘fling’ wings apart in
the smallest insects. J Exp Biol 2016; 219:3759–3772.
79. Hoover AP, Griffith BE, Miller LA. Quantifying performance in the medusan mechanospace with an actively
swimming three-dimensional jellyfish model. J Fluid Mech 2017; 813:1112–1155.
80. Hasan A, Kolahdouz EM, Enquobahrie A, Caranasos TG, Vavalle JP, Griffith BE. Image-based immersed boundary
model of the aortic root ; Submitted.
81. Kou W, Griffith BE, Pandolfino JE, Kahrilas PJ, Patankar NA. A continuum mechanics-based musculo-mechanical
model for esophageal transport ; Submitted.
82. Newren EP, Fogelson AL, Guy RD, Kirby RM. Unconditionally stable discretizations of the immersed boundary
equations. J Comput Phys 2007; 222(2):702–719.
83. LeVeque RJ, Li Z. Immersed interface methods for Stokes flow with elastic boundaries or surface tension.
SIAM J Sci Comput 1997; 18(3):709–735.
84. Li ZL, Lai MC. The immersed interface method for the Navier-Stokes equations with singular forces.
J Comput Phys 2001; 171(2):822–842.
85. Lee L, LeVeque RJ. An immersed interface method for incompressible Navier-Stokes equations. SIAM J Sci Comput 2003; 25(3):832–856.
86. Xu S, Wang ZJ. An immersed interface method for simulating the interaction of a fluid with moving boundaries. J
Comput Phys 2006; 216(2):454–493.
87. Xu S, Wang ZJ. Systematic derivation of jump conditions for the immersed interface method in three-dimensional
flow simulation. SIAM J Sci Comput 2006; 27(6):1948–1980.
88. Xu S, Wang ZJ. A 3D immersed interface method for fluid-solid interaction. Comput Meth Appl Mech Eng 2008;
197:2068–2086.
89. Xu S. A boundary condition capturing immersed interface method for 3D rigid objects in a flow. J Comput Phys
2011; 230(19):7176–7190.
90. Saad Y. A flexible inner-outer preconditioned GMRES algorithm. SIAM J Sci Comput 1993; 14(2):461–469.

Copyright c 2017 John Wiley & Sons, Ltd.
Prepared using cnmauth.cls

Int. J. Numer. Meth. Biomed. Engng. (2017)
DOI: 10.1002/cnm

