Removing the Stiffness of Elastic Force from the Immersed Boundary Method for the 2D Stokes Equations

Removing the Stiffness of Elastic Force from the Immersed Boundary Method for the 2D Stokes Equations

Thomas Y. Hou Applied and Comput. Math, Caltech, Pasadena, CA 91125. Email:    Zuoqiang Shi Applied and Comput. Math, Caltech, Pasadena, CA 91125 and Zhou Pei-Yuan Center for Applied Mathematics, Tsinghua University, Beijing 100084, China. Email:
June 16, 2007

The Immersed Boundary method has evolved into one of the most useful computational methods in studying fluid structure interaction. On the other hand, the Immersed Boundary method is also known to suffer from a severe timestep stability restriction when using an explicit time discretization. In this paper, we propose several efficient semi-implicit schemes to remove this stiffness from the Immersed Boundary method for the two-dimensional Stokes flow. First, we obtain a novel unconditionally stable semi-implicit discretization for the immersed boundary problem. Using this unconditionally stable discretization as a building block, we derive several efficient semi-implicit schemes for the immersed boundary problem by applying the Small Scale Decomposition to this unconditionally stable discretization. Our stability analysis and extensive numerical experiments show that our semi-implicit schemes offer much better stability property than the explicit scheme. Unlike other implicit or semi-implicit schemes proposed in the literature, our semi-implicit schemes can be solved explicitly in the spectral space. Thus the computational cost of our semi-implicit schemes is comparable to that of an explicit scheme, but with a much better stability property.

1 Introduction

The Immersed Boundary method was originally introduced by Peskin in the 1970’s to model the flow around heart valves. Now it has evolved into a general useful method in studying the motion of one or more massless, elastic surface immersed in an incompressible, viscous fluid, particularly in biofluid dynamics problems where complex geometries and immersed elastic membranes are present. The method has been successfully applied to a variety of problems including blood flow in the heart [25, 16, 17, 18, 26, 19, 20], vibrations of the cochlear basilar membrane [2, 8], platelet aggregation during clotting [7, 34], aquatic locomotion [5, 6, 11, 35, 3], flow with suspended particles [5, 31], and inset flight [21, 22], We refer to [27] for an extensive list of applications.

The Immersed Boundary method employs a uniform Eulerian grid over the entire domain to describe the velocity field of the fluid and a Lagrangian description for the immersed elastic structure. The force generated by the elastic structure drives the fluid and the fluid moves the elastic structure. This interaction is expressed in terms of the spreading and interpolation operations by use of smoothed Delta functions.

One of the main difficulties that the Immersed Boundary method encounters is that it suffers from a severe timestep restriction in order to keep the stability [27, 32, 30]. This has been the major limitation of the Immersed Boundary method. This restriction is typically much more severe than the one that would be imposed from using an explicit discretization for the convection term in the Navier-Stokes equation. The instability is known to arise from large boundary force and small viscosity [32]. Much effort has been made to remove this restriction. Some implicit and semi-implicit methods have been proposed in the literature [33, 23, 15]. Despite of these efforts, the timestep restriction has not been resolved satisfactorily. The computational cost of using an implicit or semi-implicit scheme is still too high to be effective in a practical computation. To date, almost all practical computations using the immersed boundary method have been performed using an explicit discretization.

In this paper, we develop several efficient semi-implicit schemes to compute the motion of an elastic interface immersed in a two-dimensional, incompressible Stokes flow. There are several important ingredients in deriving our semi-implicit schemes. The first one is to use the arclength and tangent angle formulation to describe the dynamics of the immersed interface [9]. We remark that Ceniceros and Roma have also used the arclength and tangent angle formulation to alleviate the stiffness of the viscous vortex sheet with surface tension in [4]. The second one is to obtain an unconditionally stable semi-implicit discretization of the immersed boundary problem. Throughout this paper, we use the term “stability” to mean that the energy norm of the solution can be bounded in terms of the energy norm of the initial data, which is a weaker result than proving that the difference between two solutions in the energy norm can be bounded in terms of the energy norm of their difference at time zero. The third ingredient is to perform Small Scale Decomposition to the unconditionally stable discretization to obtain our efficient semi-implicit schemes. An important feature of our small scale decomposition is that the leading order term, which is to be discretized implicitly, can be expressed as a convolution operator. This property enables us to solve for the implicit solution explicitly using the Fourier transformation. Thus, the computational cost of our semi-implicit schemes is comparable to that of an explicit method. This offers a significant computational saving in using the Immersed Boundary method.

The Small Scale Decomposition was first developed by Hou, Lowengrub and Shelley [9, 10]. They applied this method to remove the stiffness from interfacial flow with surface tension, which has proved to be very successful. Due to the coupling between the elastic boundary with the fluid, it is more difficult to remove the stiffness induced by the elastic force in the Immersed Boundary method. To remove the stiffness in the Immersed Boundary method, we need to decouple the stiffness induced by the elastic force from the fluid flow in such a way that the resulting semi-implicit discretization is still unconditionally stable. This is accomplished by using a semi-implicit discretization which preserves certain important solution structures which exist at the continuous level. Without obtaining this unconditionally stable semi-implicit discretization, a straightforward application of the Small Scale Decomposition to the Immersed Boundary method would not provide an efficient semi-implicit scheme with the desirable stability property. Very recently, Newren et al. have obtained an unconditionally stable discretization for linear force in [24]. However, they did not perform Small Scale Decomposition to their unconditionally stable discretization. As we will demonstrate in this paper, the unconditionally stable semi-implicit discretization without using the Small Scale Decomposition is still very expensive and the gain over the explicit discretization is quite limited.

We develop several efficient semi-implicit schemes for both the steady Stokes flow and the unsteady Stokes flow respectively. In both cases, our semi-implicit schemes work very well. In the steady Stokes flow, we also develop a fourth order semi-implicit scheme by using the integral factor method. For the unsteady Stokes flow, we develop a second order semi-implicit method by combining our Small Scale Decomposition with a well known second order temporal discretization [13, 27]. To illustrate the stability properties of our semi-implicit schemes, we apply our methods to several prototype problems and test our schemes for a range of elastic coefficients and viscosity coefficients. Our numerical results confirm that the semi-implicit schemes remove the high order stability constraint induced by the elastic force. In the case of unsteady Stokes equation, we also confirm the second order accuracy of our semi-implicit scheme.

This paper is organized as follows. First, we review the classical formulation of the Immersed Boundary method in Section 2. Then, we introduce the arclength and tangent angle formulation in Section 3. In Section 4, we describe the spatial discretization of the Immersed Boundary method. In Section 5-6, we develop the numerical schemes for steady Stokes flow and unsteady Stokes flow respectively. The numerical results are presented in Section 7. Our numerical studies will focus on the stability restriction and computational cost of our methods. Some concluding remarks are given in Section 8.

2 Review of the Immersed Boundary method

For simplicity, we just consider a viscous incompressible fluid in a two dimensional domain , containing an immersed massless elastic boundary in the form of a closed simple curve . The configuration of the boundary is given in a parametric form: , , tracks a material point of the boundary. We consider only the Stokes equations in this paper and would neglect the convection term. Then the governing equations are given as follows:


where is the fluid velocity, is the pressure, and are constant fluid density and viscosity, is the force density, which is not zero only on the boundary and which is infinite there. The force density can be expressed as below


denotes the two-dimensional Dirac delta function and


The choice of function in this paper is computed by Hook’s law


where is the elastic coefficient of the boundary, and is the unit tangent vector along the boundary, which is defined as


This choice of force density has been used widely in the literature in both computational and theoretical studies [12],[29],[33].

We can rewrite (3) in the following way:


Next, we introduce the spreading and interpolation operations. The spreading and interpolation operators are defined as follows:


It is easy to show that and are adjoint operators:


where the inner product are defined as follows:


Equations (1),(2) are the familiar Stokes equations of viscous incompressible fluid. Equations (3),(4) represent the interaction of the fluid and the elastic boundary. The elastic boundary applies the force to the fluid, the fluid carries the immersed boundary, and the force density is determined by the configuration of the boundary.

3 The arclength-tangent angle formulation

In studying the evolution of a curve, it is useful to represent the curve by its tangent angle and local arclength derivative . Previously, Hou, Lowengrub and Shelley [9] exploited this formulation and combined it with a so-called ”Small Scale Decomposition” reformulation to remove the stiffness induced by surface tension.

Consider the evolution of a simply closed curve with known normal and tangent velocity fields, . Assume the curve is represented by . We define the arclength derivative, , and the tangent vector, , as follows


The closed curve evolves according to


where and are the unit tangent and normal vectors of the curve respectively. According to the Frenet formula, we have , here is the arclength variable. It is easy to see that and satisfy the following evolution equations [9]:


Given and , the curve can be reconstructed up to a translation by integrating (16). However, we also need a point on the boundary to provide the constant of integration.

Using the formulation, we can reformulate the immersed boundary problem as follows:




4 Spatial Discretization

We use the spectral method to discretize both the Stokes equations and the immersed boundary equations in space since we work on periodic domains. We first discuss the discretization of the Stokes equations in a regular Cartesian grid with a uniform meshsize . Let and . The discrete Fourier transform and inverse Fourier transform are defined as follows:


Now we introduce the discrete differential operator using the discrete Fourier transform defined above. For a function defined in the fluid domain , we approximate its spatial derivatives as follows:


Denote . The differential operators are discretized in terms of :


Next, we describe the discretization of the immersed boundary. We employ a Lagrangian grid with grid space . The number of grid points along the boundary is . For a function defined on the interface , we define the discrete Fourier transform and its inverse as follows:


When the interface is a closed curve, we can approximate the derivative operator along the interface by the spectral derivative:


When the solution is not periodic, we can also use a finite difference method to discretize the derivative, we refer to [27] for more details.

Now we discuss the discretization of the spreading and interpolation operators. These two operators both involve the use of a discrete delta function. The discrete delta function we use is introduced by Peskin in [27]:




Using the above discrete delta function, we can discretize the spreading and interpolation operator as follows


The summation above is over grid points in in (41) and over grid points in in (42). Operator and are still adjoint using the following discrete inner product:


Using the inner product defined above, we have:


As we will see later, this discrete self-adjoint property is crucial in obtaining our unconditional stable semi-discrete scheme for the immersed boundary problem.

5 Steady Stokes flow

5.1 Formulation

For simplicity, we study the steady Stokes flow first. The governing equations for the steady Stokes flow are given as follows:


In this simple case, we can use a boundary integral method for the two dimension Stokes flow (page 60 of [28]) to solve equations (46)-(47) to get the velocity on the boundary:


where and


5.2 Small Scale Decomposition

As we can see from (52)-(53), the velocity field can be expressed as a singular integral with a kernel . However, the singular velocity integral is nonlinear and nonlocal. It is difficult to solve for the implicit solution if we treat the velocity integral fully implicitly. The main idea of the Small Scale Decomposition technique introduced in [9] is to decompose the singular velocity integral into the sum of a linear singular operator which is a convolution operator independent of time and the configuration of the curve, and a remainder operator which is regular. Since the remaining operator, which is nonlinear and nonlocal, is regular, the simplified convolution integral operator captures accurately the high frequency spectral property of the original velocity integral. Thus, if we treat only the leading order convolution operator implicitly, but keep the regular remainder operator explicitly, we can effectively remove the stiffness of the original velocity field which comes mainly from the high frequency modes of the solution. In this subsection, we will show how to perform such Small Scale Decomposition for the Immersed Boundary method applied to the Stokes flow.

Observe that in the integral representation of the velocity field, (52)-(53), the only singular part of the kernel is . The other part of the kernel is smooth. Thus to the leading order contribution of the velocity field can be expressed as follows:


Next, we perform a Taylor expansion for , and as a function of around . By keeping only the leading order term, we have


Substituting the above Taylor expansions to (5.2), we get


Integrating by part, we obtain


where is the Hilbert transform


Using the same method, we can get the leading order contributions of and as follows:


Note that if is a smooth function, then the commutator is a smoothing operator for . Thus we can factor a smooth function from the Hilbert transform without changing its leading order spectral property. Suppose that is smooth, then we obtain to the leading order that


Applying the same analysis to the Eqs (115)-(116) gives


Note that the leading order operator is linear. This suggests a natural semi-implicit discretization of the immersed boundary problem.

Since we are dealing with a closed immersed boundary, it is natural to work in the Fourier space. Furthermore, the Hilbert operator has a very simple kernel under the Fourier transformation. Notice that is not a periodic function of . Its value increases every time increases . Nevertheless, if we let


then is periodic. It is more convenient to work with than . Substituting (69) into (68) and taking the Fourier transform on both sides of (67),(68), we obtain


where . We have also used the fact that with being the signature function. The first term on the right hand side captures the leading order high frequency contribution of the terms from the right hand side. An important property of this leading order term is that it is linear in and and has constant coefficient in space. This provides a straightforward application of the implicit time discretization.

Since our small scale decomposition is exact near the equilibrium, we can use this result to get the stability constraint of the explicit scheme by using a frozen coefficient analysis. The stability constraint is given by


As we can see, the time step needs to be very small if is large and is small. For example, if , and , then the stability would require that .

5.3 Semi-implicit schemes

Based on the small scale decomposition presented in the previous subsection, we propose two types of semi-implicit schemes in this section. The first implicit time discretization uses the backward Euler method to discretize the leading order term while keeping the lower order term explicit. This gives rise to the following semi-implicit scheme:


We call the above discretization the semi-implicit method. Near equilibrium, the stability constraint of this numerical method is , independent of the meshsize . Since the small scale decomposition only captures the leading order contribution from the high frequency components, this method can not eliminate the effect of and completely. The coefficients and can still affect the time stability through the low frequency components of the solution, which comes from the second term of the right hand side. In order to obtain a semi-implicit discretization with better stability property, we can incorporate the low frequency contribution from the second term in our implicit discretization. This scheme can be found in the appendix . We call it the semi-implicit scheme of second kind.

The accuracy of the semi-implicit schemes presented above is just first order. In order to get a high order time discretization, we can use the integral factor method. The integral factor method factors out the leading order linear term prior to time discretization. They usually provide stable and high order time integration methods for stiff problems. To use the integral factor method, we rewrite (70),(71) as




Now it is straightforward to discretize this system to high order. In particular, we can apply the classical fourth order Runge-Kutta method to discretize the above system to obtain a fourth order semi-implicit scheme.

We remark that although the fourth order semi-implicit scheme based on the integral factor approach is much more accurate than the first order semi-implicit discretization, the stability of the fourth order method is weaker than the first semi-implicit scheme based on the backward Euler discretization. The fact that the higher order discretization gives a weaker stability property is a phenomenon which has been observed for almost all time integration methods. It is not a restriction of our semi-implicit schemes for the immersed boundary problem.

The semi-implicit schemes we describe above only update the and variables. We also need to reconstruct the boundary at from and . For this purpose, we need to update a reference point of the boundary. This will be done by using an explicit time integration method. The simplest one is the forward Euler method:


where and are evaluated at the reference point. A higher order integration method can be also used. In the explicit update of the reference point, we can use the values of and obtained using the semi-implicit discretization from the previous time steps to extrapolate the values of and in the intermediate time steps in our explicit update of the reference point. Once we have updated the reference point, we can obtain the configuration of the boundary from by integrating (16)


We can use more than one reference point, then average them to get the last configuration. This can improve the stability constraint significantly. Actually, in our computation, we use two reference points , then take the average to determine the position of the interface at next time step. Since we update only two reference points, the extra cost in updating the reference point is small compared with the overall computational cost.

6 Unsteady Stokes flow

6.1 Formulation

In this section, we will extend the semi-implicit discretization developed for the steady Stokes flow to the unsteady Stokes flow. The governing equations of the immersed boundary method for the unsteady Stokes flow are as follows:


It is much more difficult to solve the fluid velocity analytically from (84)-(85). As for the steady Stokes flow, we will first derive an unconditionally stable time discretization which will be given in next section and then apply the Small Scale Decomposition to the unconditionally stable time discretization to obtain our efficient semi-implicit schemes.

6.2 An unconditionally stable semi-implicit discretization

In this section, we will describe our unconditionally stable semi-implicit discretization of the Immersed Boundary method for the incompressible unsteady Stokes equations and prove its unconditional stability in the sense of total energy is non-increasing.

The unconditionally stable semi-implicit discretization is consisted of two steps. In the first step, we update from to , then we get in the second step.

Step 1: Update of and .


where , , , and are discrete derivative operators for the Eulerian grid and the Lagrangian grid respectively, and


Step 2: Update of . After we have obtained , and , we update at using the following semi-implicit scheme:




It is important to note that the above discretization is not fully implicit. In fact, both the spreading and interpolation operators are evaluated at the interface from the previous time step. Moreover, when solve the and , in (90) - (94), we use instead of to evaluate the force density. This makes our semi-implicit discretization linear with respect to the implicit solution variables, , , and . The above semi-implicit discretization essentially decouples the stiffness induced by the elastic force from the fluid equations. This enables us to remove the stiffness of the Immersed Boundary method effectively by applying the Small Scale Decomposition and arclength/tangent angle formulation as was done in [9].

In the following, we will prove that this semi-implicit discretization is unconditionally stable in the energy norm.

By using a discrete summation by parts, we can show that


First, we define the total energy of the physical system. The total energy includes the kinetic energy and the potential energy , which are defined below:


The total energy is then defined as


Below we will prove the unconditional stability of our semi-implicit discretization. To simplify the presentation, we still denote the discrete spectral derivative of a function as .

Taking the discrete inner product defined by (44) of (90) with and using (102), we obtain