On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems

Abstract

The immersed boundary method is both a mathematical formulation and a numerical scheme for problems involving the interaction of a viscous incompressible fluid and a (visco-)elastic structure. In [M.-C. Lai, Simulations of the flow past an array of circular cylinders as a test of the immersed boundary method, Ph.D. thesis, Courant Institute of Mathematical Sciences, New York University, 1998; M.-C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comput. Phys. 160 (2000) 705–719], Lai and Peskin introduced a formally second order accurate immersed boundary method, but the convergence properties of their algorithm have only been examined computationally for problems with nonsmooth solutions. Consequently, in practice only first order convergence rates have been observed. In the present work, we describe a new formally second order accurate immersed boundary method and demonstrate its performance for a prototypical fluid–structure interaction problem, involving an immersed viscoelastic shell of finite thickness, studied over a broad range of Reynolds numbers. We consider two sets of material properties for the viscoelastic structure, including a case where the material properties of the coupled system are discontinuous at the fluid–structure interface. For both sets of material properties, the true solutions appear to possess sufficient smoothness for the method to converge at a second order rate for fully resolved computations.

Introduction

Many problems of interest in biofluid mechanics involve the dynamic interaction of a viscous incompressible fluid and an elastic or viscoelastic structure. One approach to modeling and simulating such interactions is provided by the immersed boundary method [3], [4], [5], [6], [7], [8], [9]. In the immersed boundary formulation of such problems, the configuration of the elastic structure is described by Lagrangian variables (i.e., variables indexed by a coordinate system attached to the elastic structure), whereas the momentum, velocity, and incompressibility of the coupled fluid–structure system are described by Eulerian variables (i.e., in reference to fixed physical coordinates). In the continuous equations of motion, these two descriptions are connected by making use of the Dirac delta function, whereas a smoothed approximation to the delta function is used to link the Lagrangian and Eulerian descriptions when the continuous equations are discretely approximated for computer simulation.

A formally second order accurate version of the immersed boundary method was introduced in the Ph.D. thesis of Lai [1], [2] – prior to this work, computations performed with the immersed boundary method typically employed a variety of genuinely first order accurate schemes. The second order accuracy of the method of Lai and Peskin is formal in the sense that second order convergence rates are expected only for problems where the true solution is sufficiently smooth. However, the rate of convergence of the immersed boundary method has almost always been assessed in situations where the true solutions do not possess enough regularity for the formal convergence rate to be attained.

In the present work, we introduce a new formally second order accurate version of the immersed boundary method and demonstrate actual second order numerical convergence rates for a prototypical fluid–structure interaction problem. In our computational convergence study, we consider the interaction between a viscous incompressible fluid and a viscoelastic shell (i.e., a body which, although thin, is not infinitely thin). The numerical performance of the method is examined over a broad range of Reynolds numbers for shells with two sets of elastic properties. For the first set of elastic properties, the stiffness of the shell tapers to zero at its edges, so that there is a continuous transition in material properties between the fluid and the structure. We also consider the case where the stiffness of the shell is constant, so that there is a sharp discontinuity in the material properties of the coupled system at the fluid–structure interface. At least at low and moderate Reynolds numbers, in each situation the true solution appears to be sufficiently regular for the numerical method to converge at its formal order of accuracy as the computational grids are refined.

To our knowledge this is the first time that convergence rates in excess of first order have been documented using the immersed boundary method, although higher order convergence rates have been observed for related methods [10], [11], [12], [13]. Unlike the problem considered in the present work, previous convergence studies for the immersed boundary method have typically considered the case of an infinitely thin elastic membrane (i.e., an elastic boundary or interface). The analytic solutions to such interface problems possess discontinuities in the pressure and in derivatives of the velocity, and the immersed boundary method does not accurately capture these discontinuities. By considering the interaction between a viscoelastic shell and a viscous incompressible fluid, we avoid these difficulties and are able to obtain second order convergence rates.

The numerical scheme we present is essentially a refinement of the formally second order method of Lai and Peskin [2]. Several modifications are made to the method detailed in [2] in an attempt to reduce the occurrence of nonphysical oscillations in the computed dynamics. The simplest of these modifications is our use of a strong stability-preserving Runge–Kutta method [14] for the time integration of the Lagrangian equations of motion (i.e., the equations that specify the evolution of the configuration of the elastic structure).

We more drastically depart from [2] in our treatment of the Eulerian equations of motion, namely the incompressible Navier–Stokes equations. In the present work, the solution of these equations is by a projection method that makes use of an implicit L-stable discretization of the viscous terms [15], [16] and a second order Godunov method for the explicit treatment of the nonlinear advection terms [17], [18], [19]. Generally speaking, projection methods [20], [21], [22] are a class of fractional step algorithms for incompressible flow problems that update the velocity by first solving the momentum equation over a time interval without imposing the constraint of incompressibility. Doing so yields an “intermediate” velocity field that is generally not divergence free. The “true” updated velocity is then obtained by solving a Poisson problem to enforce the incompressibility constraint. Mathematically speaking, this process projects the intermediate velocity onto the space of divergence free vector fields.

When an “exact” projection method is used, the discrete divergence of the updated velocity is identically zero (to within roundoff error). In the present work, we employ a projection method that is not exact but rather is “approximate” in the sense that the discrete divergence of the velocity only converges to zero at a second order rate as the computational grid is refined. When such methods are used with the immersed boundary method, we have found that it is beneficial to determine the updated velocity and pressure in terms of the solutions to two different approximate projection equations at each timestep. This so-called hybrid approach was originally proposed by Almgren et al. for inviscid flow [23], and our approach is essentially an extension of their algorithm (“version 5”) to the viscous case. A more traditional projection method would employ only a single projection at each timestep. Consequently, when compared to more traditional projection algorithms, hybrid methods require the solution of additional systems of linear equations at each timestep, although this additional expense can be made modest. In the present context, we have found that the use of a more traditional projection method can result in spurious oscillations in the computed pressure. These oscillations, sometimes considered to be characteristic of the immersed boundary method [13], can be dramatically reduced by making use of the hybrid approach we present. Notably, this is an improvement that holds for thick structures (such as shells) as well as the more challenging thin interface case.

The remainder of the paper begins with a presentation of the immersed boundary formulation of the continuous fluid–structure interaction equations in Section 2. Formally second order accurate spatial and temporal discretizations of the continuous equations are then described in Section 3, although some important numerical details regarding our treatment of the nonlinear advective terms are postponed to an appendix. In Section 4, we verify that the scheme attains second order rates for a prototypical fluid–structure interaction problem, using two different sets of elastic material properties and several smoothed delta functions, and in Section 5, we demonstrate in the context of a thin interface problem that the hybrid projection method we employ reduces the magnitude of nonphysical pressure oscillations when compared to a more standard projection method. Conclusions and directions for future work are discussed in Section 6.