{'fraction_non_alphanumeric': 0.03361633120289993, 'fraction_numerical': 0.011046456167127731, 'mean_word_length': 4.610189446644546, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Immersed Boundary Method with Divergence-Free Velocity Interpolation and Force Spreading

The Immersed Boundary (IB) method is a mathematical framework for constructing robust numerical methods to study fluid-structure interaction in problems involving an elastic structure immersed in a viscous fluid. The IB formulation uses an Eulerian representation of the fluid and a Lagrangian representation of the structure. The Lagrangian and Eulerian frames are coupled by integral transforms with delta function kernels. The discretized IB equations use approximations to these transforms with regularized delta function kernels to interpolate the fluid velocity to the structure, and to spread structural forces to the fluid. It is well-known that the conventional IB method can suffer from poor volume conservation since the interpolated Lagrangian velocity field is not generally divergence-free, and so this can cause spurious volume changes. In practice, the lack of volume conservation is especially pronounced for cases where there are large pressure differences across thin structural boundaries. The aim of this paper is to greatly reduce the volume error of the IB method by introducing velocity-interpolation and force-spreading schemes with the properties that the interpolated velocity field in which the structure moves is at least C1 and satisfies a continuous divergence-free condition, and that the force-spreading operator is the adjoint of the velocity-interpolation operator. We confirm through numerical experiments in two and three spatial dimensions that this new IB method is able to achieve substantial improvement in volume conservation compared to other existing IB methods, at the expense of a modest increase in the computational cost. Further, the new method provides smoother Lagrangian forces (tractions) than traditional IB methods. The method presented here is restricted to periodic computational domains. Its generalization to non-periodic domains is important future work.

three spatial dimensions that this new IB method is able to achieve substantial improvement in volume conservation compared to other existing IB methods, at the expense of a modest increase in the computational cost. Further, the new method provides smoother Lagrangian forces (tractions) than traditional IB methods. The method presented here is restricted to periodic computational domains. Its generalization to non-periodic domains is important future work.

Introduction
The Immersed Boundary (IB) method [34] is a general mathematical framework for the numerical solution of fluid-structure interaction problems arising in biological and engineering applications. The IB method was introduced to simulate flow patterns around the heart valves [32,33], and since its success in modeling cardiac fluid dynamics [31,15,16,14], it has been extended and applied to various other applications, including but not limited to motion of biological swimmers [3,30], dynamics of red-blood cells [9] and dry foam [21,22], and rigid body motion [20,43].
The essence of the IB method as a numerical scheme lies in its simple way of coupling an Eulerian representation of the fluid and a Lagrangian representation of the structure. The force spreading linear operator S that spreads forces (stresses) from the structure to the fluid and the velocity interpolation linear operator S that interpolates velocities from the fluid to the structure are carried out via a regularized delta function δ h . One effective way to construct δ h is to require the regularized delta function to satisfy a set of moment conditions to achieve approximate grid translation-invariance and desired interpolation accuracy [1,2], thereby avoiding special grid treatment near the fluid-structure interface. In spite of its wide applicability and ease of implementation, the conventional IB method with a collocatedgrid discretization (referred to herein as IBCollocated) has two well-known shortcomings in accuracy: it achieves only first-order convergence for problems that possess sharp-interface solutions [26,19], and it can be relatively poor of volume conservation [35]. Much research effort has been put into improving the convergence rate of the IB method to second order or even higher order for problems with singular forcing at the sharp interface. Notable examples include, the Immersed Interface Method (IIM) [28,27], and more recently, a new method known as Immersed Boundary Smooth Extension [37,38]. Our focus here, however, is on improving the volume conservation properties of the IB method.
As an immediate consequence of fluid incompressibility, which is one of the basic assumptions of the IB formulation, the volume enclosed by the immersed structure is exactly conserved as it deforms and moves with the fluid in the continuum setting. Thus, a desirable feature of an IB method is to conserve volume as nearly as possible. In practice, however, it is observed that, even in the simplest case of a quasi-static pressurized membrane [17], the conventional IB method (regardless of collocated-or staggered-grid discretization) produces volume error that persistently grows in time, as if fluid "leaks" through the boundary. An intuitive explanation for this "leak" is that fluid is "squeezing" between the marker points used to discretize the boundary in a conventional IB method; however, this is not the full story, because refining the Lagrangian discretization does not improve the volume conservation of the method for a fixed Eulerian discretization.
In the conventional IB method, we can extend the notion of velocity interpolation to any point in the domain (not restricted to the immersed structure), denoted here with an italic X. The continuous interpolated velocity field can be written as U (X) = (J u)(X), where J denotes the continuous interpolation operator that interpolates the velocity at X from the discrete fluid velocity u. If a closed surface moves with velocity that is continuously divergence-free with respect to the continuum divergence operator, i.e., (∇ · U )(X) = 0, then the volume enclosed by the (deformed) surface will be exactly conserved. However, in the discrete setting, even if the interpolated velocity field is continuously divergence-free, exact volume conservation is generally not achieved because of the time-stepping error from the temporal integrator. Another source of error comes from discretizing the surface itself.
In the IB method, only a discrete collection of points on the surface, i.e., the Lagrangian markers, move according to the interpolated velocity field. A closed discretized surface can be constructed by simply connecting the Lagrangian markers defining a facet, and the resulting faceted surface by this construction does not enclose a constant volume. In the absence of temporal integration errors, this kind of volume-conservation error will approach zero as the discretization of the surface is refined. Peskin and Printz realized that the major cause of poor volume conservation of IBCollocated is that the continuous interpolated velocity field given by the conventional IB interpolation operator (denoted by J IB ) is not continuously divergence-free [35], despite that the discrete fluid velocity is enforced to be discretely divergence-free with respect to the discrete divergence operator by the fluid solver.
To improve the volume conservation of the conventional IB method, Peskin and Printz proposed a modified finite-difference approximation to the discrete divergence operator to ensure that the average of the continuous divergence of the interpolated velocity is equal to zero in a small control volume with size of a grid cell [35]. Their IB method with modified finite-difference operators (herein referred to as IBModified) was applied to a two-dimensional model of the heart, and it achieved improvement in volume conservation by one-to-two orders of magnitude compared to IBCollocated. Nevertheless, a major drawback of IBModified that limits its use in applications is its complex, non-standard finite-difference operators that uses coefficients derived from the regularized delta function (but see [21,22] for applications).
To address the issue of spurious currents across immersed structure supporting extremely large pressure differences, Guy and Strychalski [39] developed a different extension of the IB method that uses non-uniform Fast Fourier Transform [8,12] (NUFFT) to generate "spectral" approximations to the delta function, which also has superior volume conservation.
Over the past two decades, the staggered-grid (MAC) discretization has been widely adopted by the IB community [15,14,3,20,43,7]. In addition to its most celebrated feature of avoiding the odd-even decoupling in the Poisson solver that can otherwise occur with collocated-grid discretization, which leads to "checkerboard" instability in the solutions, Griffith [17] concluded from his numerical studies that the improvement in volume conservation of the IB method with staggered-grid discretization (IBMAC) is essentially the same as that of IBModified. In practice, IBMAC is more practical than IBModified in that the improvement in volume conservation directly comes as a byproduct of grid discretization without any modification to the finite-difference operators, and it is relatively straightforward to extend IBMAC to include adaptive mesh refinement [15,36] and physical boundary conditions [13]. However, we emphasize that the nature of Lagrangian velocity interpolation of IBMAC remains the same as that of IBCollocated, and, hence, there is much room for further improvement in volume conservation by ensuring that the interpolated velocity is constructed to be nearly or exactly divergence-free. We note that the methods designed to improve the convergence rate of IB methods, such as IIM [28,27] and the Blob-Projection method [5], also improve volume conservation, because the solution near the interface is computed more accurately. These methods, however, are somewhat more complex and less generalizable than the conventional IB method. This paper is concerned with further improving volume conservation of IBMAC by constructing a continuous velocity-interpolation operator J that is divergence-free in the continuous sense. The discrete IB interpolation operator S is simply the restriction of J to the Lagrangian markers. The key idea introduced in this paper is first to construct a discrete vector potential that lives on an edge-centered staggered grid from the discretely divergencefree fluid velocity, and then to apply the conventional IB interpolation scheme to obtain a continuum vector potential, from which the interpolated velocity field is obtained by applying the continuum curl operator. Note that the existence of the discrete vector potential relies on the fact that the discrete velocity field is discretely divergence-free. The interpolated velocity field obtained in this manner is guaranteed to be continuously divergence-free, since the divergence of the curl of any vector field is zero. We also propose a new force-spreading operator S that is defined to be the new adjoint of the interpolation operator S , so that Lagrangian-Eulerian interaction conserves energy. The Eulerian force density that is the result of applying this force-spreading operator to a Lagrangian force field turns out to be discretely divergence-free, so we refer to this new force-spreading operation as divergence-free force spreading. We name the IB method equipped with the new interpolation and spreading operators as the Divergence-Free Immersed Boundary (DFIB) method. As presented here, the DFIB method is limited to periodic domains.
In contrast to the local nature of interpolation and spreading in the conventional IB method, the spreading and interpolation operators of the DFIB method turn out to be nonlocal in that their construction requires the solution of discrete Poisson equations, although these operators can be evaluated efficiently using the Fast Fourier Transform (FFT) or multigrid methods. Another new feature of our method is that transferring information between the Eulerian grid and the Lagrangian mesh involves derivatives of the regularized delta function ∇δ h instead of only δ h . We confirm through various numerical tests in both two and three spatial dimensions that the DFIB method is able to reduce volume error by several orders of magnitude compared to IBMAC and IBModified at the expense of only a modest increase in the computational cost. Moreover, we confirm that the volume error for DFIB decreases as the Lagrangian mesh is refined with the Eulerian grid size held fixed, which is not the case in the conventional IB method [35]. In addition to the substantial improvement in volume conservation, the DFIB method is quite straightforward to realize from an existing modular IB code with staggered-grid discretization, that is, by simply switching to the new velocity-interpolation and force-spreading schemes while leaving the fluid solver and time-stepping scheme unchanged.
The rest of the paper is organized as follows. In Sec. 2, we begin by giving a brief description of the continuum equations of motion in the IB framework. Then we define the staggered grid on which the fluid variables live and introduce the spatial discretization of the equations of motion. Sec. 3 introduces the two main contributions of this paper: divergencefree velocity interpolation and force spreading. In Sec. 4, we present a formally second-order time-stepping scheme that is used to evolve the spatially-discretized equations, followed by a cost comparison of DFIB and IBMAC. Numerical examples of applying DFIB to problems in two and three spatial dimensions are presented in Sec. 5, where the volume-conserving characteristics of the new scheme are assessed.

Equations of motion
This section provides a brief description of the continuum equations of motion in the IB framework [34]. We assume a neutrally-buoyant elastic structure Γ that is described by the Lagrangian variables s, immersed in a viscous incompressible fluid occupying the whole fluid domain Ω ⊂ R 3 that is described by the Eulerian variables x. Eqs. (2.1) and (2.2) are the incompressible Navier-Stokes equations describing mass and momentum conservation of the fluid, in which u(x, t) denotes the fluid velocity, p(x, t) is the pressure, and f (x, t) is the Eulerian force density (force per unit volume) exerted by the structure on the fluid. In this formulation, we assume that the density ρ and the viscosity µ of the fluid are constant. The fluid-structure coupled equations are:

Spatial discretization
Throughout the paper, we assume the fluid occupies a periodic domain Ω = [0, L] 3 that is discretized by a uniform N × N × N Cartesian grid with meshwidth h = L N . Each grid cell is indexed by (i, j, k) for i, j, k = 0, . . . , N − 1. For the Eulerian fluid equations, we use the staggered-grid discretization, in which the pressure p is defined on the cell-centered grid (Fig. 1a), denoted by C, i.e., at positions x i,j,k = ((i + 1 2 )h, (j + 1 2 )h, (k + 1 2 )h). The discrete fluid velocity u is defined on the face-centered grid (Fig. 1b)  i.e., at positions , k for each component respectively. In Sec. 3, we will use these half-shifted staggered grids to construct divergence-free velocity interpolation and force spreading.
To discretize the differential operators in Eqs. (2.1) and (2.2), we introduce the central difference operators corresponding to the partial derivatives ∂/∂x α , where ϕ is a scalar grid function and {e 1 , e 2 , e 3 } is the standard basis of R 3 . We can use D h α to define the discrete gradient, divergence and curl operators: where v is a vector grid function, ijk is the totally antisymmetric tensor, and the Einstein summation convention is used here. The discrete differential operators may be defined on different pairs of domain and range (half-shifted staggered grids), and therefore, in a slight abuse of notation, we will use the same notation to denote the different operators, Although the curl operator does not appear in the equations of motion explicitly, we define it here for use in Sec. 3. The discrete scalar Laplacian operator can be defined by which yields the familiar compact second-order approximation to ∇ 2 : Note that the range and domain of L h are a set of grid functions defined on the same grid, and that grid can be C or N or any of the three subgrids of E or F on which the different components of vector-valued functions are defined. We will use the notation L h to denote the discrete vector Laplacian operator that applies (the appropriately shifted) L h to each component of a vector grid function.

Advection
We follow the same treatment of discretization of the advection term as in earlier presentations of the IB method [7]. From the incompressibility of the fluid flow ∇ · u = 0, we can write the advection term in the skew-symmetric form (2.14) Let N (u) denote the discretization of Eq. (2.14), and we define whereũ denotes an averaged collocated advective velocity whose components all live on the same grid as u α . The advective velocityũ in [7] is obtained by using the same interpolation scheme as the one used for moving the immersed structure. In our work, we simply take the average of u on the grid. For example, the three components ofũ in the x-component equation areũ Note that in the y-and z-component equations, we need different averages of u to construct u. We choose to use the wide-stencil operators in Eq. (2.15) so that the resulting grid functions are all defined on the same grid as u α . A more compact discretization of the advection term has been previously described in [17,13,42].

Fluid-Structure Interaction
The immersed structure Γ is discretized by a Lagrangian mesh of M points or markers, denoted here by a non-italic X = {X m } M m=1 , and the discrete Lagrangian force densities defined on the Lagrangian markers are F = {F m } M m=1 . As discussed in the introduction, we can extend the notion of velocity interpolation to any point X in the domain, not just restricted to the Lagrangian markers X, and define a continuous interpolated velocity field U (X) = (J u)(X). In the conventional IB method, the continuous velocity-interpolation operator J IB can be defined as (2. 16) We note that the interpolated velocity field given by Eq. (2.16) is not generally divergencefree with respect to the continuum divergence operator 1 , i.e., generally even if u is discretely divergence-free with respect to the discrete divergence operator.
The restriction of J to the collection of Lagrangian markers X defines the discrete IB interpolation operator We will also develop a new force-spreading operator S[X] that is the adjoint of the new velocity-interpolation operator S [X]. Here we use the notation [X] to emphasize that these 1 The interpolated velocity given by Eq. (2.16) has the same regularity as the regularized delta function δ h which are generally at least C 1 in the IB method, and therefore, the divergence of U is well-defined.
linear operators are parametrized by the position of the markers, as will be important when discussing temporal integration.

The discretization of the interaction equations (Eqs. (2.3) and (2.4)) can be concisely
written in the form where f is the discrete Eulerian force density defined on the appropriate subgrid of F for each component, and U = {U m } M m=1 denotes the interpolated velocities at the Lagrangian markers X.
In the conventional IB method, the force-spreading operator S IB and the velocity-interpolation operator S IB are simply discrete approximations of the surface and volume integrals in Eqs. (2.3) and (2.4), i.e., is replaced by a regularized delta function δ h to facilitate the coupling between the Eulerian and Lagrangian grids, which is taken to be of the tensor-product form where φ(r) denotes the one-dimensional immersed-boundary kernel that is constructed from a set of moment conditions to achieve approximate grid translation invariance [34,2].
In the following section, we will develop a new velocity-interpolation operator J that produces a continuously divergence-free interpolated velocity field constructed from a discretely divergence-free discrete fluid velocity.
In summary, the spatially-discretized equations of motion are

Divergence-free velocity interpolation and force spreading
This section presents the two main contributions of this paper: divergence-free velocity interpolation and force spreading. Familiarity with discrete differential operators on staggered grids and with some discrete vector identities, reviewed and summarized in Appendix A, will facilitate the reading of this section.

Divergence-free velocity interpolation
Here we introduce a new recipe for constructing an interpolated velocity field U (X) = (J u)(X) that is continuously divergence-free with respect to the continuum divergence operator, i.e., (∇ · U )(X) = 0 for all X. For now we drop the dependence on time and emphasize again that X is an arbitrary position in the domain Ω ⊂ R 3 , not just on the Lagrangian structure Γ. The main idea is first to construct a discrete vector potential a(x) that is defined on the edge-centered staggered grid E, and then to apply the conventional IB interpolation to a(x) to obtain a continuum vector potential A(X), so that the Lagrangian velocity defined by U (X) = (∇ × A)(X) is automatically divergence-free.
Suppose the discrete velocity field u(x) is defined on F and is discretely divergence-free, i.e., D h · u = 0. Let u 0 be the mean of u(x), where V = x∈F h 3 is the volume of the domain. Using the Helmholtz decomposition, we construct a discrete velocity potential a(x) for x ∈ E that satisfies where the requirement that a(x) is discretely divergence-free is an arbitrary gauge condition that makes a(x) uniquely defined up to a constant. If the gauge condition of a(x) is omitted in Eq. (3.2), then the discrete velocity potential a(x) is only uniquely defined up to G h ψ, where ψ is some unknown scalar grid function defined on N. Note that D h · a is a scalar field defined on N. In Appendix B, we prove that the discrete vector potential a(x) defined by Eq. The next step is to interpolate the discrete vector potential a(x) to obtain the continuum vector potential Lastly, we take the continuum curl of A(X) with respect to X, and our new interpolation is completed by adding the mean flow u 0 , that is, We note that the interpolation Eq. is C k , we note that the interpolated velocity field given by Eq. (3.6) is a globally-defined function that is C k−1 . We can think of the regularized delta function concentrated at X as being defined everywhere with zero outside a cube of fixed edge length (e.g. 6h for the C 3 6-point kernel [2]). Moreover, the continuity of derivatives of δ h also applies globally, including at the edges for the cube. Since the continuum vector potential defined by Eq. (3.4) is a finite sum of such C k functions, and the interpolated velocity field U (X) is obtained by differentiating A(X) once, then the resulting interpolated velocity field must have k − 1 continuous derivatives. Note that if we use an IB kernel that is C 1 , then the interpolated velocity U is C 0 , and ∇ · U is defined in only a piecewise manner. This naturally brings into question whether the volume of a closed surface is strictly conserved as the surface passes over the discontinuity of the velocity derivatives. Indeed, we observe numerically that the DFIB method offers only marginal improvement in volume conservation for C 1 kernel functions, such as the standard 4-point kernel [34], unless the Lagrangian mesh is discretized with impractically high resolution (8 markers per fluid meshwidth, see Fig. 5). By contrast, we will show that with only a moderate Lagrangian mesh size (1 to 2 markers per fluid meshwidth), the DFIB method offers a substantial improvement in volume conservation for kernels of higher smoothness, which gives a continuously differentiable interpolated velocity U . Further, we observe that volume conservation of the DFIB method improves with the smoothness of the interpolated velocity field.
In addition to the standard 4-point kernel (denoted by φ 4h ), the IB kernels considered in this paper include the C 3 5-point and 6-point kernels [1,2] (denoted by φ new 5h and φ new 6h respectively), and the C 2 4-point B-spline kernel [40], and the C 4 6-point B-spline kernel, These B-spline kernels are members of a sequence of functions obtained by recursively con-volving each successive kernel function against a rectangular pulse (also known as the window function), starting from the window function itself [40]. The limiting function in this sequence is a Gaussian [41], which is exactly translation-invarant and isotropic. The family of IB kernels with nonzero even moment conditions, such as φ 4h and φ new 6h , also have a Gaussian-like shape, but it is not currently known whether this sequence of functions also converges to a Gaussian.

The force-spreading operator
The force-spreading operator S is constructed to be adjoint to the velocity-interpolation opeartor S * so that energy is conserved by the Lagrangian-Eulerian interaction, where (·, ·) x and (·, ·) X denote the corresponding discrete inner products on the Eulerian and Lagrangian grids. In other words, the power generated by the elastic body forces is transferred to the fluid without loss, 2 where U m is the Lagrangian marker velocity at X m , and F m ∆s is the Lagrangian force where the average of f (x) over the domain is Note that we have used the summation-by-parts identity Eq. (A.5) to transfer the discrete curl operator D h × from a(x) to f (x), and thus, the grid on which the summation is performed in Eq. (3.11) is E not F. On the the right-hand side of Eq. (3.10), we substitute for U m by using the divergence-free velocity interpolation Eq. (3.6), Since u 0 and a(x) are arbitrary (except for D h · a = 0), the power identity Eq. and where ϕ is an arbitrary scalar field that lives on the node-centered grid N. Note that we have the freedom to add the term G h ϕ in Eq. (3.15), since from the identity Eq. (A.4) and Indeed, we are required to include this term since the left-hand side of Eq. (3.15) is discretely divergence-free but there is no reason to expect the first term on the right-hand side of Eq. (3.15) is also divergence-free. Note that it is not required to find ϕ in order to determine f (x), because we can eliminate ϕ by taking the discrete curl on both sides of Eq. (3.15), By imposing the gauge condition we obtain a vector Poisson equation for f (x), Note again that ∇δ h is computed on E, so that the cross-product with F m is face-centered, which agrees with the left-hand side of Eq. (3.18). Note that the solution of Eq. (3.18) can be uniquely determined by the choice of f 0 . Like our velocity interpolation scheme, the new force-spreading scheme is also non-local because it requires solving discrete Poisson equations. We remark that the new force-spreading scheme is also constructed so that the resulting force density f (x) is discretely divergence-free. This means that f (x) includes the pressure gradient that is generated by the Lagrangian forces. We do not see a straightforward way to separate the pressure gradient from f (x) in case it is needed for output purposes.

Time-stepping scheme
In this section, we present a second-order time-stepping scheme, similar to the ones developed previously [14], that evolves the spatially-discretized system Eqs. (2.24) to (2.26).
Let u n , X n denote the approximations of the fluid velocity and Lagrangian marker velocities at time t n = n∆t. To advance the solutions to u n+1 and X n+1 , we perform the following steps: Step 1. First, update the Lagrangian markers to the intermediate time step n + 1 2 using the interpolated velocity, Step 2. Evaluate the intermediate Lagrangian force density at X n+ 1 2 from the force density functional or the energy functional, and spread it to the Eulerian grid using the force-spreading scheme to get Step 3. Solve the fluid equations on the periodic grid [7], where the second-order Adams-Bashforth (AB2) method is applied to approximate the nonlinear advection term 4) and N n = N (u n ).
Step 4. In the last step, update the Lagrangian markers X n+1 by using the mid-point ap- Note that the time-stepping scheme described above requires two starting values because of the treatment of the nonlinear advection term using the AB2. To get the starting value at t = ∆t, we can use the second-order Runge-Kutta (RK2) scheme described in [34,7].
In Table 1  Therefore, the DFIB method is about two times slower than IBMAC per time step in 3D.
We point out that if the RK2 scheme [34,7] is employed rather than the scheme above, then we can save one interpolation step per time step, but the fluid equations need to be solved twice.

Numerical Results
This section presents numerical results of the DFIB method for various benchmark problems in 2D and 3D. We first consider in 2D a thin elastic membrane subject to surface tension of the membrane only. The continuum solution of this simple 2D problem has the special feature that the tangential component of the elastic force vanishes, and therefore, the normal derivative of the tangential fluid velocity does not suffer any jump across the  immersed boundary. This has the effect that second-order convergence in the fluid velocity u and the Lagrangian deformation map X can be achieved [19]. In the second set of tests, we compare volume conservation in 2D, i.e., area conservation of DFIB and IBMAC by applying them to a circular membrane under tension, and we discuss the connection between area conservation and the choice of Lagrangian marker spacing relative to the Eulerian grid size. In the third set of computations, we apply the DFIB method to a problem in which a 2D elastic membrane actively evolves in a parametrically-forced system. In the last set of numerical experiments, we extend the surface tension problem to 3D, and compare volume conservation of DFIB with that of IBMAC.

A thin elastic membrane with surface tension in 2D
It is well-known that the solutions to problems involving an infinitely thin massless membrane interacting with a viscous incompressible fluid possess jump discontinuities across the interface in the pressure and in the normal derivative of the velocity due to singular forcing at the interface [27,25]. These sharp jump discontinuities cannot be fully resolved by the conventional IB method because of the use of the regularized delta function at the interface.
Consequently, the numerical convergence rate for the Lagrangian deformation map X is gen-erally only first order even if the discretization is carried out with second-order accuracy. To achieve the expected rate of convergence, we consider problems with solutions that possess sufficient smoothness.
As a simple benchmark problem with a sufficiently smooth continuum solution we consider a thin elastic membrane that deforms in response to surface tension only. Suppose that the elastic interface Γ is discretized by a collection of Lagrangian markers X = {X 1 , . . . , X M }.
The discrete elastic energy functional associated with the surface tension of the membrane is the total (polygonal) arc-length of the interface [21], where X 0 = X M and γ is the surface tension constant (energy per unit length). The Lagrangian force generated by the energy functional at the marker X m is In our tests, we set the initial configuration of the membrane to be the ellipse To avoid artifacts in the error-norm computation because of Lagrangian markers getting too clustered during the simulation, we reparametrize the interface (for the purpose of the error computation only) from the computed markers using periodic cubic splines after each time step, and compute the l p -norm error of X based on a collection of M uniformly sampled markers X from the reparametrized interface, that is, where M does not change with N . We emphasize that the resampled markers are only used to compute the error norm and are discarded after each time step. In Fig. 2 Fig. 2) confirms that second-order convergence in u and X is achieved.

Area conservation and IB marker spacing
As an immediate consequence of fluid incompressibility, the volume enclosed by a closed immersed boundary should be exactly conserved as it deforms and moves with the fluid.
However, it is observed that even in the simplest scenario of a pressurized membrane in its circular equilibrium configuration [17], the volume error of an IB method with conventional interpolation and spreading systematically grows at a rate proportional to the pressure jump across the elastic interface [35]. In this set of tests, we demonstrate that, the "volume" or area enclosed by a 2D membrane is well-conserved by the DFIB method when the Lagrangian interface is sufficiently resolved.
We follow the same problem setup as in the test described in [17]. A thin elastic membrane Lagrangian force density on the interface is described by in which κ is the uniform stiffness coefficient. The elastic membrane is discretized by a uniform Lagrangian mesh of M points in the variable s. We approximate the Lagrangian force density by to form the regularized delta function δ h .
In Fig. 3 we compare the computational results of DFIB with those of IBMAC for different h s = 4h, 2h, h and h wiggly pattern in the passive tracers. As the the Lagrangian mesh is refined gradually from h s = 4h to h 2 (from left to right in Fig. 3), we see that u ∞ decreases from 10 −3 to 10 −7 in the DFIB method, whereas u ∞ stops improving around 10 −4 in IBMAC. Moreover, in the columns where h s = 2h, h, h 2 , we see a clear global pattern in the spurious velocity field in IBMAC, while the spurious velocity field of DFIB appears to be much smaller in magnitude and random in pattern.
We define the normalized area error with respect to the initial configuration For the coarse Lagrangian marker spacings, for example, when h s = 2h, 4h, the area errors for IBMAC and DFIB have similar orders of magnitude (compare Fig. 4a, 4b to Fig. 4d, 4e).
As the Lagrangian marker spacing is reduced from 2h to h, we see a decrease in ∆A(t; X) for IBMAC by approximately a factor of 10 (see Fig. 4b, 4c) for all the IB kernels we consider in this set of tests. In contrast, the area errors for DFIB improve by at least a factor of 10 3 for the IB kernels that are at least C 2 (see Fig. 4e, 4f), and in the best scenario, ∆A(t; X) for φ B 6h decreases from 10 −4 to 10 −9 . Moreover, as the Lagrangian mesh is refined from h to h 8 , area errors for DFIB keep improving, even approaching the machine epsilon in double precision Fig. 5e, 5f). For a moderate Lagrangian marker spacing, such as h s = h and h 2 , area errors for DFIB are several orders of magnitude smaller than those of IBMAC. On the other hand, area errors for IBMAC stop improving around 10 −5 for h s ≤ h, no matter how densely the Lagrangian mesh is refined (see Fig. 4c,   5c). We remark that the smoothness of the IB kernel appears to play an important role in volume conservation of DFIB. In this study DFIB achieves the best volume conservation result for h s ≤ h with φ B 6h , and this kernel also has the highest regularity of the kernel functions considered in this work.
The area errors of DFIB shown in Fig. 4 and Fig. 5 can be attributed to two sources of error. The first source of error is the time-stepping error from the temporal integrator, which is relatively small in the quasi-static circle test. The second source of area loss comes from discretizing the continuous curve (circle) as a polygon whose vertices are the IB markers.
This kind of error can be substantially reduced by using a high-order representation of the interface, such as a periodic cubic spline. We define a normalized area error with respect to the true initial area of the interface A true by using the tracers, where we compute A(t; X tracer ) via exact integration of the cubic spline interpolant. Fig. 6 shows that the area enclosed by the passive tracers using the cubic spline approximation is far more accurately preserved than the polygonal approximation. Furthermore, the area error approaches zero as the discretization of the tracer interface is refined, even if the marker and grid spacings are held fixed, as shown in Fig. 6. Indeed, the area of the spline interpolant through the tracers of spacing h s /20 = h/20 is conserved to the machine epsilon in double precision.
It is well-known that the traditional IB method produces non-smooth surface tractions, and a number of improvements have been proposed [11,44,45,39,29]. Somewhat unexpectedly, the divergence-free force spreading used in our DFIB method proposed here offers smoother and more accurate tractions without any post-processing such as filtering [11]. This is inherently linked to the reduced spurious flows compared to traditional methods [39].
In Fig. 7, we compare the errors of the tangential and normal components of F(s, t = 1) for DFIB and IBMAC for the quasi-static circle problem. We observe that the DFIB method dramatically improves the accuracy of Lagrangian forces by only refining the Lagrangian mesh, keeping the Eulerian grid fixed. By contrast, in the IBMAC method, the tractions do not improve as the Lagrangian grid is refined.

A thin elastic membrane with parametric resonance in 2D
In many biological applications, the immersed structure is an active material, interacting dynamically with the surrounding fluid and generating time-dependent motion. It has been reported that the simulation of active fluid-structure interactions using the conventional IB method may suffer from significant loss in the volume enclosed by the structure [35]. A , and with Lagrangian marker spacings h s ∈ {4h, 2h, h} indicated above each figure panel.    simple prototype problem for active fluid-structure interaction is a thin elastic membrane that dynamically evolves in a fluid in response to elastic forcing with periodic variation in the stiffness parameter [6,23], that is, where κ(t) is a periodic time-dependent stiffness coefficient of the form It is quite remarkable that such a purely temporal parameter variation can result in the emergence of spatial patterns, but that is indeed the case. We assume that the immersed structure is initially in a configuration that has a small-amplitude perturbation from a circle of radius R, X (s, 0) = R(1 + 0 cos(ps))r(s), (5.13) wherer(s) denotes the position vector pointing radially from the origin. For certain choices of parameters, the perturbed mode in the initial configuration may resonate with the driving frequency ω 0 in the periodic forcing, leading to large-amplitude oscillatory motion in the membrane. The stability of the parametric resonance has been studied in the IB framework using Floquet linear stability analysis for a thin elastic membrane in 2D [6,23], and recently for an elastic shell in 3D [24]. Motivated by the linear stability analysis of [6,23], we consider two sets of parameters listed in Table 2 for our simulations. The first set of parameters with τ = 0.4 leads to a stable configuration in which the membrane undergoes damped oscillations (Fig. 8a), and the second set with τ = 0.5 leads to an unstable configuration in which the membrane oscillates with growing amplitude (Fig. 8b). The computational domain Ω = [0, L] 2 is discretized by a 128 × 128 uniform Cartesian grid with meshwidth h = L 128 . The number of Lagrangian markers is determined so that the distance between the Lagrangian markers is h s ≈ h 2 in the initial configuration. The discretization of the Lagrangian force density Eq. (5.11) is constructed in the same way as Eq. (5.8). The time step size is ∆t = h 10 to ensure the stability of computation. On the left panel of Fig. 8 we show snapshots of the membrane configuration for each case, and on the right panel we plot the time-dependent amplitude (t) of the ansatz X (s, t) = R(1 + (t) cos(ps))r(s) (5.14) by applying the FFT to the Lagrangian marker positions X. In the case of growing oscillation (Fig. 8b), the amplitude of the perturbed mode increases from 0.05 to 0.3 until nonlinearities eventually stabilize the growing mode and the membrane starts to oscillate at a fixed amplitude.
We next give a direct comparison of area conservation of IBModified (with φ cos 4h [35]), IBMAC, and DFIB with the IB kernels 6h ∈ C 4 . In this test, the area enclosed the Lagrangian markers is computed by the cubic spline approximation discussed in Sec. 5.2. In Fig. 9 and Fig. 10 we show time-dependent area errors enclosed by the parametric membrane for the damped-oscillation and the growingamplitude cases respectively. For the damped-oscillation case (Fig. 9), we see that area errors for DFIB are at least two orders of magnitude smaller than those of IBMAC and IBModified for IB kernels that are at least C 2 . The volume conservation of IBModified and IBMAC was not directly compared in the previous work [17], but it was anticipated that they are similar.
In our comparison, we find that IBModified is only slightly better than IBMAC in volume conservation, yet IBMAC is much simpler to use in practice.
In this set of tests, the choice of IB kernel also plays a role in affecting area conservation.
In particular, the area errors for DFIB using φ new 6h and φ B 6h are smaller than those of φ B 4h and φ new 5h by approximately one order of magnitude. Additionally, the error curves of DFIB with φ new 6h and φ B 6h remain oscillating below 10 −7 while apparent growth of error in time is observed with φ B 4h and φ new 5h , and in the other IB methods. Unlike the quasi-static circle test in which the time-stepping error is negligible compared to the area loss due to moving a finite collection of Lagrangian markers, the time-stepping error in this example can be observed by considering a dense collection of tracers with different time step sizes. With N tracer = 4M , we first confirm that the area error of the tracer interface cannot be reduced by further including more tracers, but as we reduce the time step size ∆t ∈ h 10 , h 20 , h 40 , we observe an improvement in the area error, as shown in the bottom panel of Fig. 9. The improvements in area conservation of DFIB is consistently more than 10 4 times over IBCollocated and about 10 3 times over IBMAC. Similar results are obtained for the growing-amplitude case (see Fig. 10) except that the parametrically-unstable membrane has experienced some area loss due to the growing-amplitude oscillation before its motion is stabilized by the nonlinearities.

A 3D thin elastic membrane with surface tension
In our final test problem, we examine volume conservation of the DFIB method by extending the surface tension problem to 3D. We consider in 3D a thin elastic membrane that is initially in its spherical equilibrium configuration. The spherical surface of the membrane is discretized by a triangulation consisting of approximately equilateral triangles with edge length approximately equal to h s , constructed from successive refinement of a regular icosahedron by splitting each facet into four smaller equilateral triangles and projecting the vertices onto the sphere to form the refined mesh (see Fig. 11 for the first two levels of refinement). We use {X 1 , X 2 , . . . , X M } and {T 1 , T 2 , . . . , T P } to denote the vertices (Lagrangian markers) and the triangular facets of the mesh respectively. The generalization of discrete elastic energy functional of surface tension in 3D is the product of surface tension constant γ (energy per unit area) and the total surface area of the triangular mesh [22], that is, where |T p | is the area of the p th triangle. The Lagrangian force F k ∆s at the k th vertex is minus the partial derivative of E[X 1 , . . . , X M ] with respect to X k ,   Fig. 8b) are plotted on the semi-log scale, as done in Fig. 9 for damped motion. The computations are performed using IBMAC and DFIB with the IB kernels: where nbor(k) denotes the set of indices of triangles that share X k as a vertex 3 . Each component of the term ∂|T l |/∂X k in Eq. (5.16) can be computed analytically [22], where X k , X k , X k denote the three vertices of the triangle T l ordered in the counterclockwise direction andn is the unit outward normal vector of T l .
The computation is performed in the periodic box Ω = [0, 1] 3 with Eulerian meshwidth h = 1 128 using DFIB with φ new 6h . For the quasi-static test, the initial fluid velocity is set to be zero, and for the dynamic test, we set u(x, 0) = (0, sin(4πx), 0). In the computational results shown in Fig. 12, the spherical membrane is discretized by triangulation (as shown in Fig. 11) with 5 successive levels of refinement from the regular icosahedron (Fig. 11a) using IBMAC and DFIB with h s = h, h 2 , h 4 , which correspond to triangular meshes with 4,5,6 levels of refinement from the regular icosahedron respectively. For the quasi-static case ( Fig. 13a), volume errors for DFIB are at least 2 orders of magnitude smaller than those of IBMAC. Further, volume errors for DFIB keep decreasing as the Lagrangian mesh is refined from h s = h to h 4 . For the dynamic case (Fig. 13b), both methods suffer a significant amount of volume loss arising from the rapid deformation at the beginning of simulation. The volume error of DFIB with h s = h is similar to those of IBMAC in magnitude, but the volume error of DFIB decreases as the Lagrangian mesh is refined for h s = h 2 , h 4 . It appears that the behavior of volume error changes in nature from h s = h to h 2 , which coincides with the conventional recommendation that the best choice of Lagrangian mesh spacing in the IB method is h s = h 2 in practice. Similar to the two sources of error that contribute to the area loss in 2D, the volume error observed in Fig. 13 can also be explained by contribution from the timestepping error, and the volume loss due to only moving the vertices (Lagrangian markers) that constitute the triangular mesh. This kind of error in volume conservation decreases as the discretization of the surface is refined, even on a fixed Eulerian grid (as shown in Fig. 13). Finally, we remark that the improvement in volume conservation does not seem to be as substantial as the improvement in area conservation in 2D. We suspect that this may be attributed to the larger approximation error in computing the volume using the tetrahedral approximation (after the triangular mesh is deformed), whereas in two dimensions we use a higher-order representation of the interface (cubic splines). Nevertheless, the reduction in volume error from the Lagrangian mesh-refinement experiments indeed confirms that the DFIB method can generally achieve better volume conservation if the immersed boundary is sufficiently resolved (h s ≤ h 2 ).

Conclusions
In this paper, we introduce an IB method with divergence-free velocity interpolation and force spreading. Our IB method makes use of staggered-grid discretization to define an edge-centered discrete vector potential. By interpolating the discrete vector potential in the conventional IB fashion, we obtain a continuum vector potential whose curl directly yields a continuum Lagrangian velocity field that is exactly divergence-free by default. The corre-   using ∇δ h on the edge-centered staggered grid E. We have found that volume conservation of DFIB improves with the smoothness of the IB kernel used to construct δ h , and we have numerically tested that IB kernels that are at least C 2 are good candidate kernels that can be used to construct the regularized delta function in the DFIB method.
We have incorporated the divergence-free interpolation and spreading operators in a second-order time-stepping scheme, and applied it to several benchmark problems in two and three spatial dimensions. First, we have tested that our method achieves second-order convergence in both the fluid velocity and the Lagrangian deformation map for the 2D surface tension problem, which is admittedly a special case, since its continuum solution has a continuous normal derivative of the tangential velocity across the immersed boundary.
The highlight of the DFIB is its capability of substantially reducing volume error in the immersed structure as it moves and deforms in the process of fluid-structure interaction.  [19,10,18,4]. We also remark that the current version of the DFIB method is accompanied with a single-fluid Navier-Stokes fluid solver. This is not a fundamental limitation in our approach, and as a direction of future research, the DFIB method may be extended to work with variable-viscosity and variable-density fluid solvers [9].
Unlike other improved IB methods that either use non-standard finite-difference operators (IBModified [35]) that complicate the implementation of the fluid solver, or rely on analytically-computed correction terms (IIM [28,27] or Blob-Projection method [5]) that may not be readily accessible in many applications, the DFIB method is generally applicable, and it is straightforward to implement in both 2D and 3D from an existing IB code that is based on the staggered-grid discretization. Moreover, the additional costs of performing the new interpolation and spreading do not increase the overall complexity of computation and are modest compared to the existing IB methods.
We point out two limitations of our present work. A first limitation of the current version of DFIB method is based on the assumption of periodic boundary conditions. Extending the method to include physical boundary conditions at the boundaries of the computational domain is one possible direction of future work, but there are several challenges to overcome.
First, a special treatment of spreading and interpolation is required near the boundaries since the support of the IB kernel can extend outside of the physical domain [16,20]. Second, with non-periodic BCs, instead of FFTs, the resulting linear system needs to be solved by geometric or algebraic multigrid method to achieve high performance. For unbounded domains, a lattice Green's function technique was recently proposed as an alternative approach [29]. Third, for domains with physical boundaries, the use of projection-based fluid solvers to eliminate pressure introduces splitting errors near the physical boundaries, and instead one ought to solve a coupled velocity-pressure system at every time step [13]. In the DFIB method, it is also nontrivial to specify boundary conditions for the vector potential a(x) at physical boundaries, which would lead to a Poisson equation with non-periodic BCs. It may be that volume conservation as the structure passes near a boundary requires solving a coupled velocity-potential system. A second limitation of our DFIB method is that the pressure gradient generated by the Lagrangian forces is part of the resulting Eulerian force density because force spreading is also constructed to be discretely divergence-free. However, this may also be an important advantage of our method from the standpoint of accuracy, since it means that jumps in pressure across the interface do not require any explicit representation.
We do not yet see an obvious way to extract the pressure from the Eulerian force density in case it is needed for output purposes, or for the purposes of imposing physical boundary conditions involving tractions or avoiding splitting errors near boundaries [13].
on the periodic staggered grid just as in the continuum case, Thus, Since D h · u = 0 and D h × u = 0 by hypothesis, the left-hand side of Eq. (B.1) is also zero.
But this implies u i is constant for i = 1, 2, 3.

Lemma 2.
If ψ is a scalar grid function that lives on one of the staggered grids, such that Now we are ready to state the theorem that guarantees the existence of a discrete vector potential a(x) for x ∈ E given a discretely divergence-free velocity field u(x) for x ∈ F. The constant must be zero, however, since D h × a has zero sum by "summation by parts", and u has zero sum by assumption. This completes the proof of the existence of a vector potential satisfying u = D h × a.