Hybrid finite difference/finite element immersed boundary method

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.


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 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][5][6][7][8][9][10][11][12][13][14][15][16] including three-dimensional simulations of cardiac fluid dynamics. [17][18][19][20][21][22][23][24][25][26][27][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 (ie, 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 coworkers developed the immersed finite element (IFE) method, [30][31][32][33][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 coworkers developed a fully variational IB method that avoids regularized delta functions altogether. [37][38][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 Navier-Stokes 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][31][32][33][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 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. *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.
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][31][32][33][34][37][38][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.

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 Section 2.2. Let x = (x 1 , x 2 , … ) ∈ Ω ⊂ R d , d = 2 or 3 denote Cartesian physical coordinates, with Ω denoting the physical region that is occupied by the coupled fluid-structure system; let X = (X 1 , X 2 , … ) ∈ U ⊂ R d 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 Figure 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 a 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 fluid-structure system. In the present formulation, we assume that 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, in which p(x, t) is the pressure, I is the identify tensor, 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 P e (X, t), which is defined in terms of e via in which the deformation gradient tensor associated with the deformation ∶ ( and J(X, t) = det(F(X, t)) is the Jacobian determinant of the deformation gradient. Although P e is only defined within the solid region, it is convenient to extend e (x, t) by zero in the fluid region. 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, 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

Strong formulation
As shown by Boffi et al, 37 the strong form of the equations of motion is in which is the mass density of the coupled fluid-structure system, Du Dt (x, t) = u t (x, t) + u(x, t) · ∇u(x, t) is the material derivative, f e (x, t) is the Eulerian elastic force density, and (x) = ∏ d i=1 (x i ) is the d-dimensional delta function. Two different types of Lagrangian elastic force densities appear in these equations. The Lagrangian internal force density, ∇ X · P e , is a volumetric force density that is distributed throughout the elastic body, whereas the Lagrangian transmission force density, −P e N, is a surface force density that is distributed along the fluid-solid interface. Equation 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 P e (X, t) is sufficiently smooth, the internal force acts as a regular (ie, 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 Equation 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 Equation 9 is equivalent to 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.

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 uses a weak form of the Lagrangian equations. Because we use a finite difference scheme to approximate the incompressible Navier-Stokes 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 · P e by requiring to hold for arbitrary smooth V(X). Integrating by parts, it is clear that Equation 11 implies that for all smooth V(X). Assuming sufficient regularity, F = ∇ X · P e pointwise. It is also convenient to define the transmission force density, T(X, t), pointwise along U via T(X, t) = −P e (X, t) N(X).
With F and T so defined, the equations of motion become ∇ · u(x, t) = 0, 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 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 for all V(X). Thus, G(X, t) can be a continuous function only if P e 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 U. It is possible to use Equation 19 with standard finite element methods; however, when doing so, the transmission surface force density is effectively projected (in an L 2 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: in which g(x, t) is the Eulerian total elastic force density. This weak form of the equations of motion is essentially the formulation used in the IFE method, 30-34 the energy-based formulation of Devendran and Peskin, 35,36 and the fully variational IB method. [37][38][39] To our knowledge, the partitioned formulation described here has not been widely used in previous work.

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: 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 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.

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.

Eulerian discretization
We use a staggered-grid finite difference scheme to discretize the incompressible Navier-Stokes equations in space. Compared to collocated discretizations (ie, 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 Δx 1 = Δx 2 = h = 1 N . 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 = (u 1 , u 2 ) are approximated at the centers of the x 1 -and x 2 -edges of the Cartesian grid cells, ie, at positions ) and , respectively. A staggered scheme is also used for the Eulerian body force f = ( f 1 , f 2 ). We use the 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 p i,j . See Figure 2.
Let ∇ h · u ≈ ∇ · u, ∇ h p ≈ ∇p, and ∇ 2 h 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 ∇ 2 h 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 uses Cartesian grid adaptive mesh refinement (AMR) following the discretization approach described by Griffith. 26

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 L 2 inner product on Ω by Minor adjustments to this definition are required when nonperiodic physical boundary conditions are used 51 or near coarse-fine interfaces in locally refined grids. 26

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

Lagrangian discretization
An approximation to the deformation gradient is given by

Immersed elastic structures
For an immersed elastic structure, we use the FE approximation to the deformation gradient tensor F h (X, t) to compute P e h (X, t) and T h (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 The nodal values {F l } M l=1 and {G l } M l=1 must be determined from P e h (X, t) via discretizations of Equation 11 and Equation 19. We use a standard Galerkin approach, so that after rearranging terms, Equation 11 becomes for each m = 1, … , M. Similarly, Equation 19 becomes for each m = 1, … , M. In practice, these integrals are approximated via Gaussian quadrature.

Immersed rigid structures
For a fixed, rigid immersed structure, we directly evaluate the discretized Lagrangian penalty force F h (X, t) via Because h (X, t) and F h (X, t) are defined in terms of the same basis functions, F h (X, t) is given by in which

Lagrangian inner products
Letting [F] denote the vector of nodal coefficients of F h , we write Equation 35 as in which [] is the mass matrix that has entries of the form ∫ U l (X) m (X) dX. Equation 36 may be rewritten similarly. The mass matrix [] can also be used to evaluate the L 2 inner product of Lagrangian functions on U. In particular, for any Different choices of mass matrices (eg, 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.

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 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 = ( f 1 , f 2 ) on the Cartesian grid, we construct for each element U e ∈  h a Gaussian quadrature rule with N e quadrature points X e Q ∈ U e and weights e Q , Q = 1, … , N e . We then compute f 1 and f 2 on the edges of the Cartesian grid cells via with F(X, t) = ( F 1 (X, t), F 2 (X, t)). We use the notation in which  ( (·, t)) is the force-prolongation operator implicitly defined by Equations 42 and 43. Notice that in which ΔX is proportional to the Lagrangian mesh spacing and O(ΔX q ) corresponds to quadrature error that may be controlled by the choice of numerical quadrature rules. A corresponding velocity-restriction operator ( (·, t)) determines the motion of the nodes of the Lagrangian mesh from the Cartesian grid velocity field via There are many possible ways to construct  ; however, we have found that an effective approach is to require t =  u to satisfy the discrete power identity, which implies that the semidiscrete unified formulation conserves energy during Lagrangian-Eulerian interaction. 3 This power identity can be rewritten as ie,  =  * . To construct  explicitly, it is convenient to use matrix notation. Identifying [] and [ ] with the matrix representations of  and  , we have that Equation 49 then becomes In our time-stepping scheme, which is stated below in Appendix A.1, notice that we need only to apply [ ] to discrete velocity fields defined on the Cartesian grid. Specifically, we do not need to compute [ ] explicitly. It is straightforward to show that this construction of  implies that t (X, t) is an approximation to the L 2 projection of the Lagrangian vector field U IB (X, t) = (U IB 1 (X, t), U IB 2 (X, t)), with Because the components of U IB (X, t) are not generally linear combinations of the FE basis functions, generally t ≠ U IB . For the semidiscretized partitioned formulation, f is computed on the Cartesian grid via f =  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 =  U T to denote this operation. We use the same regularized delta function h (x) to construct both  and  U . For simplicity, we use the same velocity-restriction operator for both formulations. This choice ensures that the 2 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 in which f IB denotes the Eulerian force determined by the standard IB force-spreading operator and IB l denotes the volume associated with Lagrangian node l. In this approach, each nodal force F l 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 restriction used by such methods would be to set d l dt = U IB (X l , 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 2 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 Figure 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  and  that provide, on average, at least a 3 × 3 net of quadrature points per Cartesian grid cell.

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

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. The IBAMR provides support for distributed-memory parallelism and AMR. The IBAMR relies upon the SAMRAI, 55-57 PETSc, 58-60 hypre, 61,62 and libMesh 63,64 libraries for much of its functionality.

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 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 ≠ 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-by-N Cartesian grid. The Lagrangian discretization is constructed so that the nodes of the Lagrangian mesh are physically separated by a distance of approximately M fac Δx. Specifically, we discretize U using a 28M-by-M mesh of bilinear (Q 1 ) elements, with M = 1 M fac N 16 . Representative numerical results using N = 128 are shown in Figures 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 P e 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.

Anisotropic shell
We first consider an idealized anisotropic material model defined in terms of so that 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 s 1 direction, P e 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. Setting = 0, so that the structure is in equilibrium, and requiring ∫ Ω p(x, t) dx = 0, it can be shown 37 that . We set = 1, = 1, and e = 1, and we consider the time interval 0 ⩽ t ⩽ 3. Figure 5 summarizes the error data at time t = 3 for N = 64, 128, 256, 512, and 1024, using M fac = 1, 2, and 4, with Δt = 0.25Δx. Second-order convergence rates are observed in the L 1 , L 2 , and L ∞ norms for the velocity field. Second-order convergence rates are also observed for the pressure in the L 1 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 L 2 norm. 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 L 1 norm, at first-order in the L ∞ norm, and at order 1.5 in the L 2 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 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 Figure 6 summarizes the error data at time t = 0.75 for N = 64, 128, 256, and 512 and M fac = 1, 2, and 4, with Δt = 0.25Δx. Essentially second-order convergence rates are observed in the L 1 , L 2 , and L ∞ norms for the velocity field. The pressure converges at a second-order rate in the L 1 norm, at a first-order rate in the L ∞ norm, and at an intermediate rate (approximately 1.5) in the L 2 norm. Convergence rates for the deformation are somewhat less regular, with nearly second-order convergence rates being observed in the L 1 and L 2 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 M fac 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.

Orthotropic shell
The second case that we consider uses a neo-Hookean material model, with C = F T F and I 1 (C) = tr(C), so that 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 The velocity converges at second-order accuracy in all norms, whereas the pressure converges at second-order in the L 1 norm, at first-order in the L ∞ norm, and at order 1.5 in the L 2 norm. The displacement converges at essentially second-order rates in the L 1 and L 2 norms, and at a slightly lower rate in the L ∞ norm 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. Setting = 0, so that the structure is in equilibrium, and requiring ∫ Ω p(x, t) dx = 0, it can be shown 37 that with r = ||x − (0.5, 0.5)|| and p 0 = . We set = 1, = 1, and e = 1, and we consider the time interval 0 ⩽ t ⩽ 3. Figure 8 summarizes the error data at time t = 3 for N = 64, 128, 256, 512, and 1024, using M fac = 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 p in the L 1 norm. Because p possesses discontinuities at fluid-structure interfaces for this problem; however, the present method yields convergence rates of 0.5 in the L 2 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 Figure 9 summarizes the error data at time t = 0.75 for N = 64, 128, 256, and 512 and M fac = 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 L 1 norm.
For this problem, we find that the uni fied and partitioned formulations yield similar accuracy in most cases for u and  Figure 4, the coarse and fine structural meshes yield very similar kinematics . 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 Section 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 (x 1 = 0), right (x 1 = 1), and bottom (x 2 = 0) boundaries of Ω, and u ≡ (1, 0) on the top (x 2 = 1) boundary of Ω. We use an isotropic neo-Hookean model, and we consider both p 0 = 0 and p 0 = e . Because generally P e N ≢ 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 Figure 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 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 M fac 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 M fac Δ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. Figure 11  Setting p 0 = 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 p 0 = e , the maximum volume change is less than 0.4% in all cases considered. At smaller values of M fac , 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 volume-conservation 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% 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 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 L 1 and L 2 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.

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] × [− 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 (x 1 = −15), we set a uniform inflow velocity, u ≡ (1, 0). Along the outflow boundary (x 1 = 45), we set the normal traction and tangential velocity to zero, whereas along the top (x 2 = 30) and bottom (x 2 = −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 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 2 spatial dimensions) conservation than the unified (unsplit) formulation. Volume conservation is substantially improved for the unified formulation by setting p 0 = e . The volume change is less than 0.4% for all cases considered except for the unified formulation with p 0 = 0 and M fac = 4. discretized using an adaptively refined Cartesian grid with six nested grid levels and a refinement ratio of two between levels.
The Cartesian grid spacing on the finest grid level is Δx finest = 2 −5 Δx coarsest , with Δx coarsest = 60 N . The cylinder is discretized using a mesh of one-dimensional quadratic (P 2 ) elements with a node spacing of approximately M fac Δx finest . The time step size is Δt = 0.1Δx finest , yielding a maximum CFL (Courant-Friedrichs-Lewy) number of approximately 0.1 to 0.2. Rigidity constraints are approximately imposed using tether forces as in Equation 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 Figure 12.  Table 1 compares results obtained using the present method with the standard four-point kernel, 3 N = 128, and M fac = 2 with corresponding experimental and computational results from prior studies. We compute the lift coefficient, F y ) is the net force on the cylinder and f s is the shedding frequency. The present results are seen to be in excellent quantitative agreement with these earlier results. Figure 13 shows the lift and drag coefficients as functions of time for N = 32, 64, and 128 and for M fac = 1, 2, and 4 using the four-point regularized kernel function. 3 Similar results are shown in Figure 14 for the similar three-point kernel, 71 whereas Figure 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 (ie, 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 M fac . For the four-and six-point kernels, however, the best accuracy is obtained for larger values of M fac . In particular, by using relatively coarser structural meshes, the scheme yields more accurate forces (C L and C D ) and vortex shedding dynamics (characterized by, eg, St ). In fact, for the coarsest Eulerian discretization considered (N = 32), the six-point kernel yields erratic results for the finest structural discretization (M fac = 1) that do not occur for the coarser structural discretizations (M fac = 2 and 4). In contrast, with the three-point kernel, comparable results are obtained for all values of M fac 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.

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

FIGURE 13
Lift (C L ) and drag (C D ) 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 M fac Δx finest . Under simultaneous Lagrangian and Eulerian grid refinement, the scheme converges to the same dynamics for all values of M fac . Notice, however, that the best accuracy is obtained for larger values of M fac . Using a relatively finer structural mesh spacing results in larger lift and drag coefficients and lower vortex shedding frequencies. Thus, smaller values of M fac increase the effective numerical size of the cylinder 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) , with length measured in millimeters. The endocardial surface is defined by t = 0, u ∈ [− , − arccos 5 17 ], and v ∈ [− , ], and the epicardial surface is defined by t = 1, u ∈ [− , − arccos 5 20 ] and v ∈ [− , ]. The model ventricle is truncated at the base, which is taken to correspond to z = 5mm. See Figure 16A. 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 with is the Green-Lagrange strain tensor and C, b f , b t , b fs are material parameters. The structure is discretized using trilinear (Q 1 ) 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 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. Figures 16C,F and 17B show the results from a simulation of active contraction using this model. In this case, we use a fiber-reinforced material model, in which the fiber direction in the reference configuration E f is given by

Active contraction of a fiber-reinforced left ventricle model
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 accounted for in the first Piola-Kirchhoff stress via with T a corresponding to the active stress acting in the fiber direction. We use C = 2 kPa, b f = 8, b t = 2, b fs = 4, and T a = 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 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 Figure 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 Figure 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.

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][31][32][33][34] as well as other similar approaches. [35][36][37][38][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 (ie,  and  ) couple the volumetric structural variables to the background grid, and corresponding surface operators (eg,  U and  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 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 mesh width. 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 M fac = 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 U IB 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][84][85][86][87][88][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.

A.1 Basic time-stepping scheme
Let n , u n , and p n− 1 2 denote the approximations to the values of and u at time t n and to the value of p at time t n− 1 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 at time t n+ 1 2 via Then, we solve ( ∇ h · u n+1 = 0, (A3) for n+1 , u n+1 , and p n+ 1 2 , in which A n+ 1 2 = 3 2 u n ·∇ h u n − 1 2 u n−1 ·∇ h u n−1 is computed via a PPM-type approximation to the nonlinear advection term. 50,51 The time-stepping scheme for the unified formulation is analogous. Notice that solving Equations A2-A6 for n+1 , u n+1 , and p n+ 1 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, 90 using u n and p n− 1 2 as initial approximations to u n+1 and p n+ 1 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 Section 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 − u n Δt + A n ) = −∇ hp for̃n +1 ,ũ n+1 , andp n+ 1 2 , with A n = u n · ∇ h u n . Because we do not have an initial value for the pressure, we use p = 0 as an initial guess for p n+ 1 2 . Second, we set n+ 1 2 =̃n +1 + n 2 and solve Equations A2 -A6 for n+1 , u n+1 , and p n+ 1 2 , except that we use A n+ 1 2 = u n+ 1 2 · ∇ h u n+ 1 2 with u n+ 1 2 = 1 2 (ũ n+1 + u n ) .