Immersed Boundary Smooth Extension (IBSE): A high-order method for solving incompressible flows in arbitrary smooth domains

# Immersed Boundary Smooth Extension (IBSE): A high-order method for solving incompressible flows in arbitrary smooth domains

David B. Stein Robert D. Guy Becca Thomases Department of Mathematics, University of California, Davis, Davis, CA 95616-5270, USA
###### Abstract

The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in general domains, yet for fluid problems it only achieves first-order spatial accuracy near embedded boundaries for the velocity field and fails to converge pointwise for elements of the stress tensor. In a previous work we introduced the Immersed Boundary Smooth Extension (IBSE) method, a variation of the IB method that achieves high-order accuracy for elliptic PDE by smoothly extending the unknown solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of simple Cartesian-grid discretizations. In this work, we extend the IBSE method to allow for the imposition of a divergence constraint, and demonstrate high-order convergence for the Stokes and incompressible Navier-Stokes equations: up to third-order pointwise convergence for the velocity field, and second-order pointwise convergence for all elements of the stress tensor. The method is flexible to the underlying discretization: we demonstrate solutions produced using both a Fourier spectral discretization and a standard second-order finite-difference discretization.

###### keywords:
Embedded boundary, Immersed Boundary, Incompressible Navier-Stokes, Fourier spectral method, Complex geometry, High-order
\newdefinition

remarkRemark

## 1 Introduction

The Immersed Boundary (IB) method was originally developed for the study of moving, deformable structures immersed in a fluid, and it has been widely applied to such problems since its introduction [1, 2, 3]. Recently, the method has been adapted to more general fluid-structure problems, including the motion of rigid bodies immersed in a fluid [4] and fluid flow through a domain with either stationary boundaries or boundaries with prescribed motion [5, 6]. In this broadened context, we use the term Immersed Boundary method to refer only to methods in which (i) the boundary is treated as a Lagrangian structure embedded in a geometrically simple domain, (ii) the background PDE (e.g. the Navier-Stokes equations) are solved on a Cartesian grid everywhere in that domain, and (iii) all communication between the Lagrangian structure and the underlying PDE is mediated by convolutions with regularized -functions. These methods have many desirable properties: they make use of robust and efficient Cartesian-grid methods for solving the underlying PDE, are flexible to a wide range of problems, and are simple to implement, requiring minimal geometric information and processing to describe the boundary.

The IB method belongs to the broad category of methods known as embedded boundary (EB) methods, including the Immersed Interface [7], Ghost Fluid [8], and Volume Penalty methods [9]. These methods share a common feature: they enable solutions to PDE on nontrivial domains to be computed using efficient and robust structured-grid discretizations; yet these methods differ largely in how boundary conditions are enforced and whether or not the solution is produced in the entirety of a simple domain. Methods which compute the solution everywhere in a -dimensional rectangle admit the simplest discretizations and enable the use of high-order discretizations such as Fourier spectral methods. Unfortunately, this simplicity comes coupled with a fundamental difficulty: the analytic solution to these problems is rarely globally smooth on the entire domain. Consider the one-dimensional Poisson problem on the periodic interval with Dirichlet boundary conditions for . Even if , the solution will typically display jumps in its derivative at the values and . The lack of regularity in the analytic problem leads to low-order convergence in many numerical schemes, including the Immersed Boundary method: the addition of (regularized) singular forces supported at the boundary causes the solution to be , and solutions are accurate only to first-order in the grid spacing [4, 5].

The advantages of EB methods are substantial enough that significant effort has been expended on improving their accuracy [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Two different approaches are generally taken. The first approach involves locally altering the discretization of the PDE in the vicinity of the boundary to accommodate the lack of smoothness in the solution. One example of this approach is the Immersed Interface method [7]. Such approaches are particularly useful for interface problems where the solution is required on both sides of the embedded boundary. When the solution is only needed on one side of the interface, a second approach may be taken in which variables are redefined outside of the domain of interest to obtain higher global regularity. Improved convergence rates are achieved as a natural consequence of the properties of the discretization scheme when applied to smooth problems. Variations of this basic idea have been used by the Fourier Continuation (FC) [24, 25, 23] and the Active Penalty (AP) [26] methods to provide high order solutions to PDE on general domains.

In [27], we introduce another approach: the Immersed Boundary Smooth Extension (IBSE) method for the solution of the Poisson and related problems (e.g. the heat equation and Burgers equation). The IBSE method builds off of the basic framework of the direct forcing IB method. Drawing on ideas from the AP and FC methods, the IBSE method remedies the slow convergence of the direct forcing IB method by solving an auxiliary problem to extend the unknown solution from the physical domain into the entire computational domain. This extension is used to define the body forcing in the non-physical portion of the domain, forcing the solution to be globally . High-order accuracy is then achieved naturally, up to or the maximum accuracy supported by the underlying discretization for a smooth problem.

In this paper, we extend the work from [27] to apply to PDE requiring the imposition of a divergence constraint, including the Stokes and Navier-Stokes equations. This is a non-trivial difficulty for a high-order embedded boundary scheme based on smooth extensions. Prior work has either been restricted to the compressible Navier-Stokes equations [23], or has dealt with the pressure in an ad-hoc manner, limiting the overall accuracy of the numerical scheme to second-order [26]. These methods either extend the forcing function or solution from a prior timestep, and are thus unable to provide solutions to the Stokes problem. For time-dependent problems, they require either explicit timestepping or the use of Alternating-Direction-Implicit (ADI) schemes to take implicit timesteps, complicating the ability to achieve high-order accuracy in time.

In contrast, we work directly with the Stokes equations, rather than employing dimensional splitting or decoupling the problem into a set of Poisson equations. The smooth extensions to the solution of the Stokes equation define forcing functions in the non-physical domain in a coupled manner, and the Lagrange multipliers supported on the boundary that define these extensions are fully coupled through the Stokes problem, eliminating any need to impose artificial boundary conditions for a pressure Poisson equation. The IBSE method smoothly extends the unknown solution, allowing for efficient direct solutions of the steady incompressible Stokes equation and simple high-order in time implicit-explicit time discretization of the Navier-Stokes equations. Remarkably, the fully coupled problem for the solution to the incompressible Stokes equations (combined with the equations that define the smooth extensions to the velocity and pressure fields) can be reduced to the solution of a relatively small111With size a small multiple of the number of points used to discretize the boundary dense system of equations, two Stokes solves on a simple domain ignoring the internal boundary, and a handful of fast Fourier transforms. The dense system depends only on the boundary and the discretization, and so it can be formed and prefactored to allow for efficient time-stepping.

The IBSE method retains the essential robustness and simplicity of the original IB method. All communication between the Lagrangian boundary and the underlying Cartesian grid is achieved by convolution with regularized -functions or normal derivatives of those -functions. This allows an absolute minimum of geometric information to be used. In the traditional IB method, only the position of the Lagrangian structure must be known; the IBSE method additionally requires normals to that structure and an indicator variable denoting whether Cartesian grid points lie inside or outside of the physical domain where the PDE is defined. As with the IB method, the IBSE method is flexible to the choice of the underlying PDE discretization, so long as the mesh is uniform in the vicinity of the immersed boundary, facilitating simple coupling to existing code. In this paper, we demonstrate the IBSE solver with an underlying Fourier spectral discretization as well as a second-order staggered-grid finite difference discretization.

This paper is organized as follows. In Section 2, we introduce the methodology, ignoring the details of the numerical implementation. This section includes a review of the relevant methodology introduced in [27]. In Section 3, we discuss the numerical implementation of the method, not including the choice of the underlying discretization of the PDE. In Section 4, we study a steady Stokes problem with an analytic solution, using a Fourier spectral discretization. We analyze the convergence of the velocity and pressure fields, as well as derivatives of the velocity fields. We demonstrate up to third-order pointwise convergence of the velocity fields, as well as second-order pointwise convergence of all elements of the stress tensor. In Section 5, we analyze the flow of a viscous incompressible fluid around a cylinder in a confined channel at zero Reynolds number. We solve the problem using both a Fourier spectral discretization and a second-order finite difference discretization, comparing the accuracy and convergence of the solution, as well as the convergence of the scalar drag coefficient to known values from the literature. Finally, in Section 6, we solve the incompressible Navier-Stokes equations in both steady and unsteady cases, for the flow around a cylinder in a confined channel, demonstrating rapid convergence and quantitative agreement with a number of benchmarks from the literature.

## 2 Methods

We begin by considering the time dependent, incompressible Navier-Stokes equations: \cref@addtoresetequationparentequation

 ∂tu+u⋅∇u−νΔu+∇p =fu in Ω, (1a) ∇⋅u =0 in Ω, (1b) u =b on Γ, (1c)

where the boundary of the domain is denoted by . One way to discretize these equations in time is to treat the non-linearity explicitly, and to treat the Stokes operator implicitly. The simplest such time discretization, using forward Euler for the nonlinear terms and backward Euler for the remaining terms, is \cref@addtoresetequationparentequation

 (I−νΔtΔ)ut+Δt+Δt∇pt+Δt =ut+Δt(ft+Δtu−ut⋅∇ut) in Ω, (2a) ∇⋅ut+Δt =0 in Ω, (2b) ut+Δt =bt+Δt on Γ. (2c)

Although the Navier-Stokes equations are often discretized in time using projection methods [28], these discretizations produce large splitting errors near the boundary at low Reynolds number. Implicit-Explicit time-discretizations of the form of Equation 2 eliminate these errors, and are simpler to generalize to higher-order in time schemes (see Section 6). The unsplit Stokes system may be inverted efficiently using GMRES with a projection preconditioner [29]. We find this inversion strategy to be very efficient (convergence to an algebraic tolerance of achieved in 3-4 iterations) for moderate to high Reynolds numbers (), and moderately efficient (convergence achieved in approximately 35 iterations) for the Stokes problem where projection methods cannot be used.

Motivated by the form of Equation 2, in this section we focus primarily on developing the methodology required to solve the general Stokes problem: \cref@addtoresetequationparentequation

 (αI−Δ)u+∇p =fu in Ω, (3a) ∇⋅u =fp in Ω, (3b) u =b on Γ. (3c)

When and , Equation 3 reduces to the steady incompressible Stokes equations. The case where arises in the study of porous media [30] and from discretizing the Navier-Stokes equations in time using an implicit-explicit discretization, as in Equation 2. The non-zero divergence constraint () adds no additional complexity and is included for completeness. Before discussing the IBSE method as applied to Equation 3, we first review relevant material from [27], in order to introduce the central ideas in the simplified context of the Poisson problem.

### 2.1 Review of [27]: the IBSE method for the Poisson problem

Suppose that we wish to solve the Poisson problem \cref@addtoresetequationparentequation

 Δu =f in Ω, (4a) u =g on Γ, (4b)

in an arbitrary smooth domain . The IBSE method works by smoothly extending the unknown solution of a PDE from to a simplier domain . The domain will be referred to as the physical domain and will be referred to as the computational domain. It is assumed that . The interior boundary will be assumed to be smooth, not self-intersecting, and must separate into the two disjoint regions and , referred to as the extension domain. Neither nor need be connected. An additional simple computational domain is required to solve the auxiliary equation that defines the extension of the unknown solution. Although may coincide with , it need only contain . The boundary of the computational domain will be denoted by . Two typical domains are shown in Figure 1. In Figure 1, is periodic, and no boundary conditions need to be specified on . In Figure 1, Dirichlet conditions are imposed on . To prevent confusion with the interior boundary , the boundary and its respective boundary condition will only be explicitly included in equations when required to prevent confusion, e.g. in Sections 5.2 and 3.4.

We will denote the extension of the unknown solution by . This extension of the solution is then used to define a volumetric forcing in the region . With this forcing, an extended problem in all of the simple domain may be solved: \cref@addtoresetequationparentequation

 Δue−χEFe =χΩf in C, (5a) ue =g on Γ. (5b)

The solution gives the desired solution in and is equal to in . Because was chosen to be a smooth extension to , the function is globally smooth in .

The extension to the unknown function is defined as a solution to a high-order PDE which takes for its boundary conditions matching criteria of the form . This allows the extension to be defined by a small number of unknowns (proportional to the number of points used to discretize the boundary). The extension PDE for is solved efficiently in the simple domain using an Immersed Boundary type method.

In order to succinctly describe the methodology we will require some notation. We define the spread operator:

 (S(j)F)(x)=(−1)j∫ΓFj(s)∂jδ(x−X(s))∂njdX(s) (6)

and the interpolation operator:

 (S∗(j)ξ)(s)=(−1)j∫Cξ(x)∂jδ(x−X(s))∂njdx (7)

for the normal derivative, where is a parametrization for with in the parameter interval . The integral on the right hand side of Equation 7 is notation for the action of the distribution on the smooth function . The traditional Immersed Boundary spread and interpolation operators, and are denoted by and respectively. The interpolation operator maps function values (or the normal derivatives of the function) from to , while the spread operator maps singular (and hyper-singular) forces supported on to . To simplify equations presented in the forthcoming material, we also define the operators , , and by: \cref@addtoresetequationparentequation

 Tk =k∑j=0S(j), (8a) T∗k (8b) R∗k =(S∗(1)⋯S∗(k))⊺. (8c)

The operator provides an interpolation of a function and its first normal derivatives to the boundary; provides an interpolation of the first normal derivatives to the boundary, but excludes the values; the spread operator represents a set of singular forces (-like) and hyper-singular forces (like the first normal derivatives of the -function) on the boundary.

The central challenge of the IBSE method is to compute the smooth extension to an unknown solution. We first discuss how to compute an extension to a given function. Let be given. To compute a extension to , we solve the following high-order PDE in the region : \cref@addtoresetequationparentequation

 Hkξ =0 in E, (9a) ∂jξ∂nj =∂jv∂nj on Γ, 0≤j≤k. (9b)

Here is an appropriate differential operator such as the polyharmonic operator ; we discuss the choice of this operator in more detail, from a numerical perspective, in Section 3.2. This problem may be solved on the simpler domain using methodology directly analogous to the direct forcing Immersed Boundary method. The boundary conditions given in Equation 9b that force to share its first normal derivatives with along are enforced by the addition of unknown singular and hyper-singular forces supported on the boundary: \cref@addtoresetequationparentequation

 Hkξ(x)+(TkF)(x) =0 for x∈CE, (10a) (S∗(j)ξ)(s) =∂jv∂nj(s) for s∈IΓ, 0≤j≤k. (10b)

Future equations similar to Equation 10b will be shortened simply to , and assumed to hold for all . Notice that is not actually an extension to : that is, in . We will only be interested in the function in , and so need not form its literal extension (which is .

To solve the Poisson problem given by Equation 4 using the IBSE method, we instead solve the extended problem given in Equation 5. The forcing function that is specified in must be chosen so that it forces the extended solution to be . Let smoothly extend , that is, we ask that is globally smooth in and that it satisfies the constraints

 R∗kξ=R∗ku (11)

at the interface . These constraints require that the first normal derivatives of agree with the first normal derivatives of on the boundary. Notice that it is not enforced that the values of agree with the values of on ; we will see momentarily that this is unnecessary. The forcing function is then defined as:

 Fe=Δξ. (12)

Coupling these equations together, we obtain the IBSE formulation for the Poisson problem given by Equation 4: \cref@addtoresetequationparentequation

 Δue−χEΔξ =χΩf in C, (13a) Hk+TkF =0 in CE, (13b) R∗kξ =R∗kue, (13c) S∗ue =0. (13d)

These equations are (13a) the extended Poisson problem, (13b) the extension PDE, (13c) the interface constraints that force and to share their first normal derivatives, and (13d) the physical boundary condition to the Poisson problem. Both the interface constraints and the physical boundary condition are enforced by the unknown forces in the extension PDE. We will subsequently drop the notation and refer to the solution of the IBSE system of equations simply as .

In , satisfies the physical Poisson equation; while in , satisfies: \cref@addtoresetequationparentequation

 Δu =Δξ in E, (14a) ∂u∂n =∂ξ∂n on Γ, (14b)

and thus in up to a constant. Since implies that , then for , the forcing , and so , and is thus at least (for spatial dimensions )222The notation denotes the space of measurable functions with the property that the function itself, as well as its first two weak derivatives, are square integrable. The Sobolev embedding theorem implies that any function in is equal almost everywhere to a continuous function, at least in spatial dimensions and [31].. This continuity, along with the smoothness constraints guaranteed by the constraint that imply that is globally .

In [27], we verify that the IBSE formulation of the problem produces solutions that converge at a rate of in the grid-spacing , for the Poisson problem, as well as the heat equation, Burgers equation, and the Fitzhugh-Nagumo equations.

{remark}

In [27], both a volumetric force and a singular force distribution are added to the Poisson equation, giving a slightly different formulation for the IBSE method:

 Δue−χEΔξ+SG =χΩf in C, (15) Hk+TkF =0 in CE, (16) T∗kξ =T∗kue, (17) S∗ue =0. (18)

The singular force distribution acts to enforce that the boundary values of equal those of , this is an unnecessary requirement for the IBSE method. We will proceed with the simpler formulation given in Equation 13 when developing the methodology for the Stokes equation.

### 2.2 Solution of the general Stokes problem

We now focus on the modifications to the IBSE method required to solve the general Stokes problem: \cref@addtoresetequationparentequation

 Lu+∇p =fu in Ω, (19a) ∇⋅u =fp in Ω, (19b) B(u,p) =b on Γ. (19c)

We assume that and are sufficiently smooth functions defined on (for ) and (for ). Here denotes the Helmholtz operator and denotes a boundary operator. We discuss the boundary operator associated with the Dirichlet problem (); imposition of other boundary conditions is similar and discussed in [27]. Note that this equation reduces to the incompressible Stokes problem when and are zero.

In the direct forcing Immersed Boundary (IB) method, this problem is solved in a simple computational domain by adding singular forces supported on the boundary that act as Lagrange multipliers which enforce the boundary condition. For Dirichlet problems this can be represented as \cref@addtoresetequationparentequation

 Lu+∇p+SG =fu in C, (20a) ∇⋅u =fp in C, (20b) S∗u =b. (20c)

The singular forces supported on induce jumps in the normal derivatives of at the boundary; the velocities produced by the IB method are generically while the pressure field typically has jump discontinuities. For Eulerian discretizations of the underlying PDE that are ignorant of the boundary, this leads to the slow convergence rate for the velocity field in . The pressure field fails to converge near to the boundary, but does converge at the rate in .

In order to improve this convergence, the Immersed Boundary Smooth Extension (IBSE) method extends the unknown solutions and to be and , respectively. Suppose that the solution vector to Equation 19 is known in . Define by \cref@addtoresetequationparentequation

 =0 in CE, (21a) R∗kξu =R∗ku, (21b)

and by \cref@addtoresetequationparentequation

 Hk−1ξp+Tk−1Fp =0 in CE, (22a) T∗k−1ξp =T∗k−1p. (22b)

Recall that is a high-order differential operator such as , which will be defined precisely in Section 3.2. We make two remarks concerning the extension of the pressure. {remark} The pressure function is typically thought of as a Lagrange multiplier: it equals whatever is required to enforce the divergence constraint that . One may suspect therefore that it is unnecessary to construct an explicit extension of the pressure function as we do in Equation 22. Unfortunately, this is not true. Consider the problem: \cref@addtoresetequationparentequation

 −Δu+∇p =fu in Ω, (23a) ∇⋅u =0 in Ω, (23b) u =b on Γ, (23c)

along with the modified problem where is replaced by for any scalar function . The pressure will be changed: . But the solution will be left unchanged, and thus its extension functions will be unchanged. If the extension of the pressure function depends only locally on , it cannot extend both and smoothly. {remark} The choice to extend the pressure function to be rather than yields improved numerical stability. When , the structure of the Stokes equations implies that . Applying the operator to the function with that global regularity yields an inconsistent estimate of , i.e. , and thus enforcing that does not actually enforce that .

From these extensions, we may define a set of forces : \cref@addtoresetequationparentequation

 Fu =Lξu+∇ξp, (24a) Fp =∇⋅ξu. (24b)

Analogous to Equation 5, to smoothly extend the Stokes problem given by Equation 19 from to , we add to Equation 19a and to Equation 19b in the extension region to obtain \cref@addtoresetequationparentequation

 Lu+∇p =χΩfu+χELξu+χE∇ξp in C, (25a) ∇⋅u =χΩfp+χE∇⋅ξu in C, (25b) B(u,p) =b on Γ. (25c)

In the physical domain , we have that \cref@addtoresetequationparentequation

 Lu+∇p =fu in Ω, (26a) ∇⋅u =fp in Ω, (26b) B(u,p) =b on Γ, (26c)

and thus solves the original Stokes problem given in Equation 19. In the extension domain , we have that \cref@addtoresetequationparentequation

 L(u−ξu)+∇(p−ξp) =0 in E, (27a) ∇⋅(u−ξu) =0 in E, (27b) ∂(u−ξu)∂n =0 on Γ, (27c)

Thus and , both up to a constant (in ). Since is enforced to match at the boundary by Equation 22, they are equal at the boundary, and thus in , implying that is globally smooth in , while the global regularity for is implied by the same argument given in Section 2.2. Numerical approximations of these solutions may converge rapidly, even when the PDEs are discretized using an underlying Cartesian grid discretization that is ignorant of the boundary. The difficulty in this formulation comes from the fact that the extensions and depend on the unknown solution .

Combining Equations 19, 22 and 21, we obtain a coupled system for the solution together with the smooth extensions to that solution:
\cref@addtoreset equationparentequation (28a) (28b) (28c) (28d) (28e) (28f) (28g)

We will refer to this system of equations as the IBSE- equations. The remainder of this paper is dedicated to demonstrating that these equations provide accurate solutions to the Stokes and Navier-Stokes equations, that they allow a flexible choice of discretization, and that they admit an efficient inversion strategy. We organize the presentation of information as follows:

1. In Section 3, we discuss the details of numerical implementation that do not depend upon the underlying choice of discretization.

2. In Section 4, we solve a Stokes problem with an analytic solution that allows us to carefully verify the convergence rates of the solution (up to third-order for and second-order for , in ), as well as the convergence of elements of the fluid stress tensor (up to second-order in ).

3. In Section 5, we describe both a Fourier spectral and a finite-difference discretization to the confined channel flow around a cylinder problem for incompressible Stokes flow. We validate the convergence rates for solutions produced by the IBSE method and compare to known benchmarks.

4. In Section 6, we discuss time-stepping in the IBSE framework, solve both steady and unsteady Navier-Stokes problems, and demonstrate rapid convergence of the solutions and agreement with known benchmarks.

## 3 Numerical implementation

In this section, we discretize the coupled IBSE- system given in Equation 28. We assume that the functions are discretized using an underlying Cartesian mesh with uniform333For the discretization described here, the mesh needs to be uniform only in a small neighborhood of the boundary, with the same width as the regularized -function used to discretize the operators and . In fact, the mesh need not even be uniform, but this would add some complexity to the discretization [32]. grid-spacing but otherwise make no assumptions regarding the underlying discretization of the functions and differential operators. These final details will be treated in Section 4 for a Fourier spectral discretization, and in Section 5.2 for a standard second-order finite difference scheme. Little varies in the details.

1. In Section 3.1, we discretize the spread and interpolation operators for the normal derivative that were introduced in Equations 7 and 6. Note that the discretization of these operators automatically induces a discretization for the composite operators , , and , defined in Equation 8.

2. In Section 3.2, we discuss how we choose the extension operator introduced in Equation 10.

3. In Section 3.3, we describe an efficient inversion strategy for the IBSE- system given by Equation 28.

4. In Section 3.4, we provide a brief summary of the algorithm.

### 3.1 Discretization of singular integrals

Let the boundary be parametrized by the function . In all examples in this manuscript, we work with a two-dimensional fluid, and the boundary is one-dimensional and closed, with the single parameter defined on the periodic interval . The discrete version of the spread operator for the normal derivative, defined in Equation 6 requires a regularized -function and a discretization of the integral over . The regularized -function that we use is analytically three times differentiable, satisfies four discrete moment conditions, and has a support width of . The -function and its properties are introduced in [27, 33]; this discrete function will be represented by . Multivariate -functions are computed as Cartesian products of the univariate and also denoted by . Normal derivatives of are computed by the formula [34]

 ∂j~δ∂nj=ni1⋯ni2∂j~δ∂xi1⋯∂xij, (29)

where the Einstein summation convention has been used to indicate sums over repeated indices and
is computed as Cartesian products of the appropriate derivatives of the one dimensional . For example, in two dimensions is computed as

 ∂~δ∂n=nx~δ′⊗~δ+ny~δ⊗~δ′. (30)

Discretization of the integral over is made by choosing quadrature nodes , equally spaced in the parameter interval so that and . Quadrature weights are computed at each quadrature node to be ; this is a spectrally accurate quadrature rule for the integral of smooth periodic functions on . The discrete spread operator maps functions sampled at points in to by

 (S(j)F)(x)=nbdy∑i=1F(si)∂j~δ(x−Xi)∂njΔsi. (31)

We do not adopt explicit notation to distinguish between the analytic and discrete operators. The number of points in the quadrature is chosen so that . This choice of node-spacing is wider than that recommended for the traditional IB method [1] but has been observed empirically to be the optimal choice in other studies of direct-forcing IB methods [4]. We define the interpolation operator by the adjoint property , but note that the discrete interpolation operator produces a discrete function

 (S∗(j)u)(sk)=∫Cu(x)∂j~δ(x−Xk)∂njdx. (32)

Discrete integrals over are straightforward sums computed over the underlying uniform Cartesian mesh, and may be efficiently evaluated due to the finite support of . When acting on vector functions, all operators defined in this section are assumed to operate elementwise.

### 3.2 Solution of the extension equations and the choice of extension operator Hk

Due to the difficulty in accurately and efficiently inverting the high-order differential operator needed to define the smooth extensions, we solve these equations utilizing a Fourier spectral discretization regardless of the discretization used for the Stokes equations. Let be the largest wave-number associated with the discrete Fourier transform (DFT) on the discretized domain . We choose the extension operator introduced in Equation 9 to be:

 Hk=Δk+1+(−1)k+1Θ(k,m). (33)

Here is a positive scalar function that depends on the smoothness of the extension () and the largest wave-number (). The function is chosen to mitigate the numerical condition number of the operator :

 κ=1+m2(k+1)Θ. (34)

Minimizing the condition number (by taking to be large) must be balanced against the need to resolve the intrinsic length scale introduced to the problem. We set the length scale as a function of , in effect allowing the extension to decay to zero over some number of grid cells. We thus choose

 Θ∗=(1NΔx)2(k+1), (35)

where is proportional to the number of gridpoints over which the extension decays. The choice of can impact the accuracy of the scheme: if is too large, the condition number of may be high with respect to the precision of the computer being used, and numerical instability causes a degradation of the quality of the solution. If is too small, short length scales that are not resolvable by the discretization are introduced to the problem. For all examples presented in this paper, is chosen to be . In Section 4, we investigate the effect of the choice of on a specific problem.

### 3.3 Inversion of the IBSE-k system, Equation 28

In this section, we describe an efficient inversion strategy for Equation 28. We will assume that the discrete Stokes operator is invertible444The null-space in the pressure function should be fixed by the discrete operator; e.g. by enforcing that .. The details for when it is not invertible (e.g. when solving on a periodic domain with ) are analogous to the Poisson problem, and the solution is presented in [27]. In matrix form, Equation 28 is given by:

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝L∇−χEL−χE∇∇⋅−χE∇⋅HkTkHk−1Tk−1S∗R∗k−R∗kT∗k−1−T∗k−1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝upξuξpFuFp⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝χΩfuχΩfp00b00⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (36)

This system is large, but the number of constraints (the last three equations) is relatively small. By block Gaussian elimination, we can find a Schur-complement for this matrix, which we label :

 SC=⎛⎜⎝S∗R∗k−R∗kT∗k−1−T∗k−1⎞⎟⎠⎛⎜ ⎜ ⎜ ⎜⎝L∇−χEL−χE∇∇⋅−χE∇⋅HkHk−1⎞⎟ ⎟ ⎟ ⎟⎠−1⎛⎜ ⎜ ⎜⎝0000Tk00Tk−1⎞⎟ ⎟ ⎟⎠ (37)

along with an associated system of equations for the Lagrange multipliers and ,

 (38)

The size of is comparatively small, only square. This equation is the key to the efficiency of the algorithm, as it allows the Lagrange multipliers and to be computed rapidly without first solving for , , , or . Once and are determined, and may be found by solving Equations 28d and 28c. With and known, and may be computed simply by solving the Stokes equation in , ignoring the boundary (that is, by solving Equations 28b and 28a).

Rapid inversion of the system of equations for and given in Equation 38 is not a trivial task. Because we have restricted to problems set on stationary domains and the size of is small, it is feasible to proceed using dense linear algebra. We form by repeatedly applying it to basis vectors. This operation is expensive: for two dimensions it is in the total number of unknowns . The Schur-complement is then factored by the LU algorithm provided by LAPACK [35]. Once this factorization is computed Equation 38 can be solved rapidly. This Schur-complement depends only on the domain and the discretization, so the LU-decomposition can be reused to solve multiple problems on the same domain or in each timestep when solving time-dependent problems.

### 3.4 Outline of the methodology

We provide a brief outline of the steps required for solving the general Stokes problem given by Equation 19 using the IBSE- method. The algorithm proceeds as follows:

1. Form the right hand side of the Schur-complement equation, as defined in Equation 38. This may be done by solving \cref@addtoresetequationparentequation

 Lu0+∇p0 =χΩfu in C, (39a) ∇⋅u0 =χΩfp in C, (39b) B(u0,p0) =bC on ∂C, (39c)

for and then computing , and . Note that the internal boundary is ignored when computing the solution to Equation 39.

2. Compute and by solving the Schur-complement system given by Equation 38, given the right hand side computed above. We assume that the Schur complement has been formed and its LU decomposition found as a pre-computation.

3. Compute and from and by solving Equations 28d and 28c. We take to always have periodic boundary conditions, and these equations may be simply inverted using the fast Fourier transform.

4. With and known, find and by solving Equations 28b and 28a. To be precise, we solve \cref@addtoresetequationparentequation

 Lu+∇p =χELξu+χE∇ξp+χΩfu in C, (40a) ∇⋅u =χE∇⋅ξu+χΩfp in C, (40b) B(u,p) =bC on ∂C. (40c)

Note again that the internal boundary is ignored when computing this solution.

Step two of this algorithm requires that the Schur-complement be formed and factored as a pre-processing step. The Schur complement operator maps a set of singular (and hyper-singular) forces to the constraints given by equations Equations 28g, 28f and 28e. We form the matrix by applying to basis vectors. To apply to any set of singular forces, first compute and from and by solving Equations 28d and 28c. Then and may be found by solving \cref@addtoresetequationparentequation

 Lu+∇p =χELξu+χE∇ξp in C, (41a) ∇⋅u =χE∇⋅ξu in C, (41b) B(u,p) =0 on ∂C, (41c)

and finally with , , and known, then , , and can be evaluated. Note that , and do not enter into this computation. The Schur complement depends only on the geometry and the discrete operators , , .

## 4 Results: Stokes equation with analytic solution

To demonstrate the improved convergence properties of the IBSE method as compared to the IB method, we construct a simple example with an analytic solution and compare the accuracy of the solutions with those produced by the Immersed Boundary method. For the numerical comparisons with the IB method provided in this section, we show two sets of values: those produced using a regularized -function commonly used in IB simulations [10]:

 8δ(x)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩3−2|x|+√1+4|x|−4|x|20≤x≤1,5−2|x|−√−7+12|x|−4|x|21

as well as those produced using the smoother and more accurate function introduced in [27] and used in all IBSE simulations in this paper. These two methods will be denoted by and , respectively. We solve the general Stokes problem with and Dirichlet boundary conditions: \cref@addtoresetequationparentequation

 u−Δu+∇p =fu in Ω, (43a) ∇⋅u =0 in Ω, (43b) u =b on Γ, (43c)

where the solution is given by \cref@addtoresetequationparentequation

 u(x,y) =esinxcosy, (44a) v(x,y) =−cosxesinxsiny, (44b) p(x,y) =ecos2x, (44c)

the forcing function is defined as

 fu={u−Δu+∇pin Ω,0in E, (45)

and the boundary condition is found by evaluating on . The physical domain is taken to be , where denotes the 2-torus identified with the periodic rectangle , and where denotes the ball of radius centered at . The computational domain is taken to be . Plots of the analytic solution in the physical domain are shown in Figure 2.

For this problem, it is convenient to use a Fourier spectral discretization for and to choose . We discretize the domain using a regular Cartesian mesh with points discretizing each dimension: and for . Differential operators are discretized in the usual way.

A refinement study for the and errors in the solutions for and is shown in Figure 3. The refinement path for is omitted; it is similar to what is shown for . In , the velocity field converges at first-, second-, and third-order accuracy in for the IB, IBSE- and IBSE- methods, respectively, while the pressure field converges at first- and second-order accuracy for IBSE- and IBSE-, but at the slow rate of for the IB method. In , the rates of convergence for the velocity remain unchanged in all methods, while for the pressure the difference is starker: the IBSE- and IBSE- methods maintain first- and second-order convergence, while the IB method fails to converge pointwise.

{remark}

It may be surprising that the velocity fields produced by the Immersed Boundary method converge at only first-order in . Indeed, for a given and , the velocity fields produced by solving \cref@addtoresetequationparentequation

 −Δu+∇p+SF =0, (46a) ∇⋅u =0, (46b) S∗u =b, (46c)

converge at , , and in the , , and norms, respectively [36]. However, is not given, it is computed by inverting a Schur-complement equation whose elements include terms of the form where , and thus contain errors.

In the refinement study shown in Figure 3, the pressure function converges more slowly than the velocity function; in the case of the IB method, it fails to converge pointwise. This will be true in general for elements of the stress tensor, which are one derivative less smooth than the velocity field. In Figure 4, we show a refinement study for the and errors in the derivatives and (the convergence paths for and are omitted, but essentially the same). A similar pattern to the convergence of the pressure function is observed: for all of the derivatives, the IB method converges slowly (at ) in , but fails to converge pointwise; the IBSE- and IBSE- method provide first- and second-order convergence in both and . We emphasize that with pointwise convergence of the pressure function as well as pointwise convergence of these derivatives, the IBSE method is able to capture the stress tensor pointwise up to the boundary .

Finally, we use this test problem to numerically investigate the choice of the free parameter in the definition of the extension operator . This parameter controls the rate at which the extension functions and decay to away from the boundary. We compute solutions with , using extensions across a wide range of choices of (from to ). The error for , and the condition number of the Schur complement, as a function of , are shown in Figures 5 and 4(a). For very small values of , the error is dominated by the inability of the discretization to resolve the small length scale introduced to the problem. As is increased, the error produced by the IBSE method improves, but the condition number of the Schur-complement worsens. As grows large, the error again begins to degrade, dominated by the error produced when inverting the ill-conditioned Schur-complement. This presents a balance between resolution and conditioning. For this test, the minimum error is somewhere between and , a result that has been fairly robust across the range of test problems we have examined. Other than for this test, all numerical results in this paper are produced using with to construct . In Figures 5 and 5, we show the extension function for the solution to Equation 43 for and .