The immersed boundary method: A projection approach

Journal of Computational Physics 225 (2007) 2118–2137 www.elsevier.com/locate/jcp The immersed boundary method: A projection approach Kunihiko Taira *, Tim Colonius Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA Received 11 April 2006; received in revised form 11 December 2006; accepted 12 March 2007 Available online 19 March 2007 Abstract A new formulation of the immersed boundary method with a structure algebraically identical to the traditional fractional step method is presented for incompressible flow over bodies with prescribed surface motion. Like previous methods, a boundary force is applied at the immersed surface to satisfy the no-slip constraint. This extra constraint can be added to the incompressible Navier–Stokes equations by introducing regularization and interpolation operators. The current method gives prominence to the role of the boundary force acting as a Lagrange multiplier to satisfy the no-slip condition. This role is analogous to the effect of pressure on the momentum equation to satisfy the divergence-free constraint. The current immersed boundary method removes slip and non-divergence-free components of the velocity field through a projection. The boundary force is determined implicitly without any constitutive relations allowing the present formulation to use larger CFL numbers compared to some past methods. Symmetry and positive-definiteness of the system are preserved such that the conjugate gradient method can be used to solve for the flow field. Examples show that the current formulation achieves second-order temporal accuracy and better than first-order spatial accuracy in L2-norms for oneand two-dimensional test problems. Results from two-dimensional simulations of flows over stationary and moving cylinders are in good agreement with those from previous experimental and numerical studies.  2007 Elsevier Inc. All rights reserved. MSC: 76D05; 76M12 PACS: 47.11.+j Keywords: Immersed boundary method; Fractional step method; Projection method; Staggered grid; Finite-volume method; Incompressible viscous flow 1. Introduction Immersed boundary methods (IBMs) have gained popularity for their ability to handle moving or deforming bodies with complex surface geometry [33,27]. Peskin [31] first introduced the method by describing the flow field with an Eulerian discretization and representing the immersed surface with a set of Lagrangian * Corresponding author. Tel.: +1 626 395 4766. E-mail addresses: [email protected] (K. Taira), [email protected] (T. Colonius). 0021-9991/$ - see front matter  2007 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2007.03.005 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2119 points. The Eulerian grid is not required to conform to the body geometry as the no-slip boundary condition is enforced at the Lagrangian points by adding appropriate boundary forces. The boundary forces that exist as singular functions along the surface in the continuous equations are described by discrete delta functions that smear (regularize) the forcing effect over the neighboring Eulerian cells. Peskin originally used the IBM to simulate blood flow inside a heart with flexible valves, where the forcing function was computed by Hooke’s law [31,32]. This technique was later extended to rigid bodies by taking the spring constant to be a large value [3,22]. Goldstein et al. [17] applied the concept of feedback control to compute the force on the rigid immersed surface. The difference between the velocity solution and the boundary velocity is used in a proportional-integral controller. For the aforementioned techniques to model flow over rigid bodies, the choice of gain (stiffness) remains ad hoc and large gain results in stiff equations. Our intention is to remove all tuning parameters and formulate the IBM in a general framework for rigid bodies (as well as bodies with prescribed surface motion). In our formulation, we treat the boundary forces in a manner analogous to the discretized pressure. For the incompressible Navier–Stokes equations, pressure may be viewed as a Lagrange multiplier required to satisfy the divergence-free constraint. Similarly, boundary forces can be regarded as Lagrange multipliers that satisfy the no-slip constraint [16]. By introducing regularization and interpolation operators and grouping the pressure and force unknowns together, the discretized incompressible Navier–Stokes equations can in fact be formulated with a structure algebraically identical to the traditional fractional step method. Although previous research has implemented immersed boundary techniques with the traditional fractional step algorithm, the entire IBM itself has not been regarded as a fractional step (projection) method, as reported here. We follow the algebraic approach of Perot [30], where the fractional step method is written as a block-LU decomposition. In the next section, we review the traditional fractional step method as it is the fundamental basis for our IBM. In Section 3, we introduce the immersed boundary projection method. This formulation is compared to previous methods in Section 4; namely the original IBM [31], the direct forcing method [28], the immersed interface method (IIM) [23], and the distributed Lagrange multiplier (DLM) method [16]. In Section 5, numerical examples are considered to assess the temporal and spatial accuracy of the current method. Flows over stationary and moving cylinders are simulated and results are compared to previous experimental and numerical studies. Section 6 summarizes the current formulation. Further details on the discretization of the immersed boundary projection method are placed in the Appendix. 2. Fractional step method We consider the incompressible Navier–Stokes equations ou 1 þ u  ru ¼ rp þ r2 u; ot Re r  u ¼ 0; ð1Þ ð2Þ where u, p, and Re are the suitably non-dimensionalized velocity vector, pressure, and the Reynolds number, respectively. Following Refs. [8,37,19,30,6], the equations are discretized with a staggered-mesh finite-volume formulation using the implicit Crank–Nicolson (CN) integration for the viscous terms and the explicit secondorder Adams-Bashforth (AB2) scheme for the convective terms. This produces an algebraic system of equations,   nþ1   n    A G r bc1 q ¼ ; ð3Þ þ D 0 bc2 0 / where qnþ1 and / are the discretized velocity flux and pressure vectors. The discrete velocity can be recovered by unþ1 ¼ R1 qnþ1 , where R is a diagonal matrix that transforms the discrete velocity into the velocity flux. Sub-matrices G and D correspond to the discrete gradient and divergence operators, respectively. The operator resulting from the implicit velocity term is A ¼ Dt1 M  12 L, where M is the (diagonal) mass matrix and L is the discrete (vector) Laplacian. We construct the Laplacian to be symmetric, hence making A symmetric as well. The right-hand side of Eq. (3) consists of the explicit terms from the momentum equation, rn, and 2120 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 inhomogeneous terms from the boundary condition, bc1 and bc2. Details on the discretization of Eqs. (1) and (2) can be found in the Appendix and [30,6]. It is interesting to note that G ¼ DT for the staggered grid formulation. The traditional fractional step method by Chorin [8] and Témam [37] was introduced to solve Eq. (3) in an efficient manner by using an approximation for A1. In the present analysis, we adopt the observation made by Perot [30] that the fractional step method can be regarded as an LU decomposition of Eq. (3):  A 0 GT GT BN G  I BN G 0 I  qnþ1 /  ¼  rn 0  þ  bc1 bc2  N N þ  Dt2N ðLM 1 Þ G/ 0 ! ; ð4Þ where BN is the Nth order Taylor series expansion of A1: A1 ffi BN ¼ DtM 1 þ N X Dt2 DtN Dtj N 1 j1 ðM 1 LÞM 1 þ    þ N 1 ðM 1 LÞ M 1 ¼ ðM 1 LÞ M 1 : j1 2 2 2 j¼1 ð5Þ The last term in Eq. (4) is the leading order error resulting from the truncation in BN. Let us note that BN is symmetric and can be made positive-definite with appropriate choices of Dt and N [30]. In the current situation, there also exists a second-order temporal discretization error from the AB2 and CN methods. As discussed in [30], the fractional step error can be absorbed by the discrete pressure if LM 1 and G are commutative (for example, in the case of periodic domains); otherwise there remains an Nth order error. Eq. (4) is more commonly written in three steps: Aq ¼ rn þ bc1 T T  N G B G/ ¼ G q þ bc2 q nþ1  N ¼ q  B G/ ðSolve for intermediate velocityÞ; ðSolve the Poisson equationÞ; ðProjection stepÞ: ð6Þ ð7Þ ð8Þ Since both A and GT BN G are symmetric positive-definite matrices, the conjugate gradient method can be utilized to solve the above momentum and Poisson equations in an efficient manner. In general, for non-symmetric matrices, various other Krylov solvers can be employed. Here the discrete pressure is denoted by / without any superscript for its time level, as we regard pressure as a Lagrange multiplier [6]. There has been extensive discussion on the exact time level of the discrete pressure variable for various treatments of pressure in fractional step methods [36,5]. For the present method, / is a first-order accurate approximation of pressure in time, vis. pnþ1=2 . Since the first-order accuracy of / does not affect the temporal accuracy of the velocity variable [30], we use / as a simple representation of the pressure variable. If a second-order accurate pressure is desired, Brown et al. [5] should be referred to for further modifications to the fractional step method. Although a detailed stability analysis is not offered in this paper, we demonstrate that the present method described in the next section can stably solve for the flow field for CFL numbers up to 1, as shown in Section 5. We mention that fractional step methods for incompressible flow can suffer numerical instability if Dt is decreased arbitrarily [18]. The time step is limited by a lower bound of Dt P cDxlþ1 if equal orders of interpolation are used for velocity and pressure, as in the present case (c is a constant and l is the interpolation order of velocity, here l = 2). While remedies are offered in [18,9], we have not utilized them here since a much larger Dt is usually selected based on the CFL constraint. We note in passing that the form of Eq. (3) is known as the Karush–Kuhn–Tucker (KKT) system that appears in constrained optimization problems [29]. This system minimizes a term similar to the kinetic energy:   1 T T min ðqnþ1 Þ Aqnþ1  ðqnþ1 Þ ðrn þ bc1 Þ subject to Dqnþ1 ¼ 0 þ bc2 : ð9Þ qnþ1 2 It is interesting that the discrete pressure / does not play a direct role in time advancement, but acts as a set of Lagrange multipliers to minimize the system energy and satisfy the kinematic constraint of divergence-free velocity field. K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2121 3. Immersed boundary projection method 3.1. The discretized Navier–Stokes equations with boundary force Since the discretized Navier–Stokes equations, Eq. (3), are observed to be a KKT system with pressure acting as a set of Lagrange multipliers to satisfy the continuity constraint, one can imagine appending additional algebraic constraints by increasing the number of Lagrange multipliers. Hence we incorporate the no-slip constraint from the IBM into the fractional step framework. The IBM introduces a set of Lagrangian points, nk , that represent the surface, oB, of an immersed object, B, within a computational domain, D, whose geometry need not conform to the underlying spatial grid. At the Lagrangian points, appropriate surface forces, fk, are applied to enforce the no-slip condition along oB. Fig. 1 illustrates the setup of the spatial discretization. Since the location of the Lagrangian boundary points does not necessarily coincide with the underlying spatial discretization, two operators are required: one that passes information from the boundary points to the neighboring staggered grid points and another one that conveys information in the opposite direction. We consider the continuous version of the incompressible Navier–Stokes equations and explain how the IBM can be discretized into a KKT system and solved with a fractional step/projection algorithm. The incompressible Navier–Stokes equations with a boundary force, f, and the no-slip condition can be considered as the continuous analog of the IBM Z ou 1 2 þ u  ru ¼ rp þ r u þ fðnðs; tÞÞdðn  xÞ ds; ð10Þ ot Re s r  u ¼ 0; ð11Þ Z uðnðs; tÞÞ ¼ uðxÞdðx  nÞ dx ¼ uB ðnðs; tÞÞ; ð12Þ x where x 2 D and nðs; tÞ 2 oB. The boundary oB, parametrized by s, is allowed to move at a velocity uB ðnðs; tÞÞ. Convolutions with the Dirac delta function d are used to allow the exchange of information from oB to D and vice versa in Eqs. (10) and (12), respectively. The discretization of the above system results in 1 30 nþ1 1 0 n 1 0 2 r A G H bc1 q C B C 7B 6 C B 0 5@ / A ¼ @ 0 A þ @ bc2 A; ð13Þ 4D 0 nþ1 uB 0 E 0 0 f Fig. 1. Staggered grid discretization of a two-dimensional computational domain D and immersed boundary formulation for a body B depicted by a shaded object. The horizontal and vertical arrows (!; ") represent the discrete ui and vi velocity locations, respectively. Pressure pj is positioned at the center of each cell (·). Lagrangian points, nk ¼ ðnk ; gk Þ, along oB are shown with filled squares (n) where boundary forces f k ¼ ðfx;k ; fy;k Þ are applied (); *). 2122 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 T where Hf corresponds to the last term in Eq. (10) with f ¼ ðfx ; fy Þ . Similar to the discrete pressure, we do not place a superscript for time level on f to emphasize its role as a Lagrange multiplier. The no-slip condition, Eq. (12), is enforced using the constraint, Eqnþ1 ¼ uBnþ1 . Here A, G, and D are the implicit operator for the velocity flux, the discrete gradient operator, and the discrete divergence operator, respectively, and rn, bc1, and bc2 are the explicit terms in the momentum equation, the boundary condition vector resulting from the Laplacian operator, and the boundary condition vector generated from the divergence operator, respectively. Note that these sub-matrices and vectors (A, G, D, rn, bc1, and bc2) are identical to those that appear in the traditional discretization, Eq. (3). The two additional sub-matrices H and E are introduced to regularize (smear) the singular boundary force over a few cells and interpolate velocity values defined on the staggered grid onto the Lagrangian points, respectively. We will refer to these sub-matrices as regularization (H) and interpolation (E) operators. The precise expressions of these operators are discussed below and details on the overall discretization are provided in the Appendix. 3.2. Interpolation and regularization operators The operators H and E are constructed from the regularized discrete delta function. Among the variety of discrete delta functions available, we choose to use the one by Roma et al. [34] which is specifically designed for its use on staggered grids (where even/odd de-coupling does not occur). This function has the form: 8 " ffi# rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 > > jrj jrj 1 > 5  3 Dr  3 1  Dr þ 1 for 0:5Dr 6 jrj 6 1:5Dr; > 6Dr > > <  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dðrÞ ¼ ð14Þ 2 > 1 > 1 þ 3 Drr þ 1 for jrj 6 0:5Dr; > 3Dr > > > : 0 otherwise; where Dr is the cell width of the staggered grid in the r-direction. This discrete delta function is supported over only three cells, which is an advantage for computational efficiency. We have not found significant differences in the results for the current formulation with alternative discrete delta functions. Refs. [33,3] may be consulted for a variety of delta functions. As observed by Peskin [31] and Beyer and LeVeque [3], discrete delta functions can be used both for regularization and interpolation. The interpolation operator can be derived from discretizing the convolution of u and d, Z uðnÞ ¼ uðxÞdðx  nÞ dx ð15Þ x yielding uk ¼ DxDy X i ui dðxi  nk Þdðy i  gk Þ ð16Þ for the two-dimensional case, where ui is the discrete velocity vector defined on the staggered grid ðxi ; y i Þ and uk is the discrete boundary velocity at the kth Lagrangian point ðnk ; gk Þ. For the three-dimensional case an extra factor of Dzdðzi  fk Þ is needed. Letting a denote the factor preceding the summation, the interpolation operator for Eq. (16) can be written as ^ k;i ¼ adðxi  nk Þdðy i  gk Þ; E ð17Þ so that the no-slip condition is represented by ^ k;i unþ1 ¼ Ek;i qnþ1 ¼ unþ1 ; E i i Bk ð18Þ K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2123 ^ 1 to allow the use of the flux, qnþ1 ¼ Runþ1 , from the fractional step formulation. The hat is used where E  ER to represent the original operator and is removed once a transformation (e.g. R1) is applied. Similarly, the regularization operator is a discrete version of the convolution operator in Eq. (10) that passes information from the Lagrangian points, nk, to the neighboring staggered grid points, xi. Defining H in a manner similar to E, we obtain ^ iE ^ i dðnk  xi Þdðgk  y i Þ ¼ b M ^T ; H i;k ¼ bM k;i a ð19Þ ^ is included for where b is the numerical integration factor proportional to ds. Note that a diagonal matrix M consistency with the fractional step formulation. It should be observed that E and H are symmetric up to a ^ are absent. constant if the diagonal matrices R1 and M ^ in Eq. Next, let us achieve symmetry between the (3, 1) and (1, 3) block entries in the presence of R1 and M (13). We absorb the offset in scaling into the unknown boundary force by introducing a transformed forcing function f~ that satisfies Hf ¼ ET f~ : T~ ð20Þ The original boundary force can be retrieved by f ¼ invðEH ÞEE f . In the case of using a uniform Cartesian grid with Dx ¼ Dy, the relation simplifies to f ¼  Dx1 2 ba f~ . The discrete delta function of Eq. (14) currently requires the use of a uniform grid in the vicinity of oB to satisfy a set of properties (e.g. moment conditions) [34]. Since the range and domain of E and H are only limited to the neighborhood of oB, non-uniform discretization can still be applied away from the body. Although it is not pursued here, it would be interesting to generate discrete delta functions that are suitable for a nonuniform spatial discretization around the immersed body. Note that symmetry between E and H is not necessary for discretization, but it allows us to solve the overall system in an efficient manner. There are unexplored possibilities using different discrete delta functions for interpolation and regularization operators. Beyer and LeVeque [3] consider such cases in a one-dimensional model problem. 3.3. Immersed boundary method via projection Now that we have formulated the sub-matrices G and D such that D ¼ GT and introduced a transformed force, f~ , the overall system of equations, Eq. (13), becomes 2 A 6 T 4G E G 0 0 30 nþ1 1 0 n 1 0 1 r bc1 q ET C B C B 7B C 0 5@ / A ¼ @ 0 A þ @ bc2 A: uBnþ1 0 f~ 0 ð21Þ As previously discussed, both the discrete pressure and boundary forcing functions are Lagrange multipliers and, algebraically speaking, it is no longer necessary to make a distinction between the two. Thus organizing the sub-matrices and vectors in Eq. (21) in the following fashion:     bc2 / T n ; ð22Þ Q  ½G; E ; k  ~ ; r1  r þ bc1 ; r2  unþ1 f B Eq. (21) can be simplified to a KKT system  A Q QT 0  qnþ1 k  ¼  r1 r2  ; ð23Þ which is now in a form identical to Eq. (3), providing motivation to apply the same fractional step technique in solving the overall system as in Section 2. Performing an LU decomposition of Eq. (23), 2124 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137  A QT 0 T N Q B Q  I BN Q 0 I  qnþ1 k  ¼  r1 r2  þ ! N N  Dt2N ðLM 1 Þ Qk : 0 ð24Þ As in the original fractional step method, there is an Nth order splitting error. Note that this error cannot be absorbed by the Lagrange multiplier, k, because LM 1 and Q do not commute (even for periodic domains). Hence, a third-order expansion for BN is recommended, as discussed in [30] and Section 5. Thus, the immersed boundary projection method consists of the same three steps as Eqs. (6)–(8) but with k replacing / and Q replacing G: Aq ¼ r1 T N T  Q B Qk ¼ Q q  r2 q nþ1  N ¼ q  B Qk ðSolve for intermediate velocityÞ; ðSolve the modified Poisson equationÞ; ðProjection stepÞ: ð25Þ ð26Þ ð27Þ The main differences between the present and the traditional fractional step methods are in the Poisson equation and the projection step. Here, the pressure and boundary force are determined implicitly from the modified Poisson equation. The projection step removes the non-divergence-free and slip components of the velocity from the intermediate velocity field in one step. The numerical constraint of no-slip exists only at the Lagrangian points, hence making the dimensions of H and f~ considerably smaller than those of G and /. Thus it is encouraging that there is no significant increase in size of QT BN Q in the modified Poisson equation from GT BN G in the classical fractional step method. We can still solve Eqs. (25) and (26) with the conjugate gradient method as both left-hand side operators are symmetric and positive-definite. Some care must be taken to make QT BN Q positive-definite and well-conditioned. First, as in the traditional fractional step method, one of the discrete pressure values must be pinned to a certain value to remove the zero eigenvalue.1 Second, no repeating Lagrangian points are allowed to avoid QT BN Q from becoming singular. Also, to achieve a reasonable condition number and to prevent penetration of streamlines caused by a lack of Lagrangian points, the distance between adjacent Lagrangian points, Ds, is set approximately to the Cartesian grid spacing. In the case of moving immersed bodies, the location of the Lagrangian points must be updated at each time and so must E, i.e. nþ1 Ek;i ¼ Ek;i ¼ Eðnk ðtnþ1 Þ; xi Þ ð28Þ and similarly for H. These operators can be pre-computed at each time step by knowing the location of the Lagrangian points a priori. The current technique is not limited to rigid bodies and can model flexible moving bodies if we are provided with the location of oB at time level n + 1. For deforming bodies, the volume of the body must be isochoric to satisfy the incompressibility constraint. The current formulation treats the density of the body and the outer fluid to be equal to each other. 4. Comparison with other immersed boundary methods Let us compare our current formulation with a few other IBMs, in particular the original IBM [31], the direct forcing approach [28,15], the IIM [23], and the DLM method [16] to clarify the fundamental differences. Since we only select a few IBMs that are most similar to the current formulation, Refs. [33,27] should be consulted for additional IBMs. The same notation introduced earlier is used in this section. Because the comparison of fundamental mechanisms for satisfying the no-slip condition along the immersed boundary is of interest here, we consider methods for simulating both rigid and elastic bodies. Some details such as the time integration schemes, the updating algorithms for the Lagrangian points, and the constitutive relations for the 1 There are alternatives to pinning the solution of the modified Poisson equation. Bochev and Lehoucq [4] discuss such techniques in detail for the Poisson equation with a Neumann boundary condition. Although the current staggered grid formulation does not require any explicit pressure boundary conditions, their analysis provides insight into the algebraic properties of the discretized Poisson equation. K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2125 boundary forces are omitted for clarity of discussion. The discrete spatial operators and the temporal treatment of the discrete pressure variable may not be identical to our version but remain conceptually similar. 4.1. The original immersed boundary method (IBM) The original IBM [31] is a modification to the traditional fractional step method, Eqs. (6)–(8), to simulate flow over a flexible body. An explicit boundary force term Hf n computed with Hooke’s law is added to the right-hand side of the momentum equation Aq ¼ rn þ bc1 þ Hf n ; T T  N G B G/ ¼ G q þ bc2 ; q nþ1  N ¼ q  B G/: ð29Þ ð30Þ ð31Þ At every time step, the location of the Lagrangian points on the elastic surface is updated. Although it is not considered here, a source/sink can be added to the pressure Poisson equation to apply a correction to the continuity equation [20]. Let us discuss how the original IBM may conceptually be related to our method. Hooke’s law can be written as: f ¼ jðne  nÞ, where j is the spring constant and ne is the equilibrium position for the boundary surface. If we are to differentiate and discretize this relation, we obtain: f nþ1  f n ¼ jðunþ1  Eqnþ1 Þ; B Dt ð32Þ using the implicit Euler time discretization. Adding the boundary force to the momentum equation, we observe that the overall system has the form: 1 2 30 nþ1 1 0 n r þ bc1 A G H q C B C B 6D 0 bc2 0 7 ð33Þ A: 4 5@ / A ¼ @ 1 1 n nþ1 nþ1 uB þ jDt f E 0 jDt I f For rigid body simulations, j  1 is chosen to reduce the effect from the (3, 3) sub-matrix [3,22]. In the limit of j ! 1, we recover our current formulation, Eq. (13). The above formulation, Eq. (33), has a structure identical to the artificial compressibility method [7] that approximately satisfies the continuity equation with: 1 op þ r  u ¼ 0, where a is an artificial speed of sound. This artificial parameter is typically set to a large value a2 ot similarly R to the spring constant, j, in Eq. (33). Instead of Hooke’s law, a feedback controller t (f ¼ j1 0 uðn; sÞ ds  j2 uðn; tÞ) with large gains (j1  1 and j2  1) has also been used to compute the boundary force [17], which results in an identical structure to Eq. (33). However, large gains used in such constitutive relations add stiffness to the governing system, thus prohibiting the use of high CFL numbers. For instance, CFL numbers used in [22,17] are Oð103 Þ to Oð101 Þ for simulations of flow over a rigid circular cylinder. It is possible to use higher CFL numbers by lowering the gains at the expense of relaxing the no-slip condition. In contrast, the current projection method solves for the boundary force implicitly with no constitutive relations and behaves similarly to the traditional fractional step method in terms of temporal stability. Hence simulations can be performed with CFL numbers as high as 1, which is reported later in Section 5. In previous methods, it is not clear how the gains or the magnitude of the forcing function relate to how well the no-slip condition is satisfied. On the other hand, our method satisfies the continuity equation and the no-slip condition exactly to machine precision or, if desired, to a prescribed tolerance. 4.2. The direct forcing method The direct forcing method [28] approximates the boundary force for rigid bodies with an intermediate velocity field q*. The force is not actually computed but implemented directly into the momentum equation by substituting the regularized no-slip condition near the immersed boundary. Conceptually speaking, the momentum equation, Eq. (25), is modified to yield 2126 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 1 e  HEÞðrn þ bc1 Þ þ 1 Hunþ1 ; HEq ¼ ð M Dt Dt B T N T  G B G/ ¼ G q þ bc2 ; e  HEÞAq þ ðM q nþ1  N ¼ q  B G/: ð34Þ ð35Þ ð36Þ Here HE interpolates and then regularizes a vector, which acts as a filtering operator to extract the velocity e is placed for scaling such that M e  HE 0 near oB. Factors of field near oB. A diagonal mass matrix M 1=Dt are inserted in Eq. (34) to keep the order with respect to Dt consistent (note that A ¼ Oð1=DtÞ). Conceptually, the above equation becomes Eq ¼ uBnþ1 near the immersed boundary and reduces to Aq ¼ ðrn þ bc1 Þ away from the body. The difference between the modified momentum equation, Eq. (34), and the momentum equation from the traditional fractional step method, Eq. (6), can be expressed as the boundary force for the direct forcing method: f nþ1 ¼ uBnþ1  Eq þ EAq  Eðrn þ bc1 Þ: Dt ð37Þ Note that this method enforces the no-slip condition on q* but not on qnþ1 . A projection step is applied later to project the intermediate velocity, q*, onto the solenoidal solution space. In order to satisfy the no-slip condition exactly, iterations over the entire fractional step algorithm is required for each time level. Although slip in qnþ1 is reported to be small [15], the magnitude of the error cannot be estimated in a deductive manner. 4.3. The immersed interface method (IIM) Next, we consider representing the IIM [23] for elastic membranes in an algebraic form. In the IIM, the boundary force is decomposed into tangential and normal components (fs and fn, respectively). A regularized tangential component of the force, Hf ns , is included in the momentum equation as an explicit term and the explicit normal boundary force is implemented into the pressure Poisson equation in terms of a pressure jump condition across the interface. The overall method can be described as Aq ¼ rn þ bc1 þ Hf ns ; T T  N T G B G/ ¼ G q þ bc2 þ G B q nþ1  N ¼ q  B ðG/  bðfnn ðsptÞÞ bðfnn ÞÞ; N bðfnn Þ; ð38Þ ð39Þ ð40Þ where b ¼ is a corrective term to calculate the pressure gradient (G/  b) taking the jump condition, spt, into consideration. Since the normal component of the boundary force is implemented directly into the pressure Poisson equation rather than in the momentum equation, a sharp velocity solution in the vicinity of the interface can be achieved resulting in second-order spatial convergence for some test problems. However, the construction of the correction term b requires explicit knowledge of the boundary force, and is not easily made implicit as desired in our formulation. We note in passing that Linnick and Fasel [24] recently developed a high order IIM that employs one-sided finite differences to obtain jump conditions for higher-order derivatives. Their results along with other numerical and experimental studies for flow over a stationary cylinder are compared to our results in Section 5. 4.4. The distributed Lagrange multiplier (DLM) method The most similar method to our formulation is the DLM method by Glowinski et al. [16] used in a variational principle (finite element) framework. Their work is closely related to ours as they introduce Lagrange multipliers (i.e. body force) on the immersed rigid body to satisfy the no-slip condition, essentially through projection. The main difference between our formulation and the DLM method lies in how the projection is applied to the velocity field. Conceptually speaking, we consider the DLM method as a different operator splitting applied to Eq. (13). Their overall system is solved with the Marchuk–Yanenko fractional step scheme [39,26] that decomposes the overall operations into three operators related to: (i) the divergence-free condition and pressure, (ii) the con- K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2127 vective and diffusive operators, and (iii) the no-slip condition and boundary force. Because the projection operators that remove the non-divergence-free and no-slip conditions are applied separately at different sub-time levels, these two constraints cannot be simultaneously satisfied by the velocity field. In our formulation, there is only one projection step that simultaneously removes both the non-divergencefree and slip component from the velocity field. We also note that our formulation achieves second-order accuracy in time by choosing a suitable approximation for A1. 4.5. Short summary on the comparisons In the first three approaches, the presence of an immersed object is treated as a corrective term to account for the no-slip condition. The fundamental difference between the aforementioned methods and our formulation is the implicit treatment of both the pressure and boundary force as a single set of Lagrange multipliers in the modified Poisson equation. Once the pressure and the force are determined, the continuity equation and the no-slip condition are satisfied through a projection at the same time level in our formulation. The DLM method is found to be the most similar method but differs in how the projections are applied. Our overall IBM is viewed as a projection method to allow further generalization and numerical investigation from an algebraic point of view. 5. Results We numerically investigate the temporal and spatial convergence of the current method in one- and twodimensional model problems; namely the Stokes’ problem and flow inside two concentric cylinders, respectively. Also, flow over a circular cylinder is considered to validate the current method in steady-state and transient flow. At last, a moving body example of an impulsively started circular cylinder is considered. Since the present method is a combination of the immersed boundary and the fractional step methods, we expect convergence analyses from both methods to carry over to the current formulation. The temporal accuracy of the immersed boundary projection method should follow the analysis from the fractional step algorithm as shown in Eq. (24). In all of the problems below, second-order finite-volume discretization (except for H and E) is applied. For the problems of flow over a cylinder, a non-uniform grid is employed, making the scheme formally first-order accurate. However, we suppress the first-order spatial error by using a very smooth grid stretching, effectively keeping the overall error to second-order. In the vicinity of the body, the spatial grid is kept uniform with its finest resolution and Dxmin ¼ Dy min Ds. Unless stated otherwise, N = 3 is chosen for approximating A1. 5.1. One-dimensional Stokes’ problem We first assess the accuracy of the current method using a one-dimensional Stokes’ problem where an infinitely long flat plate is impulsively set into motion with uwall ¼ 1 in an initially quiescent viscous fluid with m = 1 (Fig. 2). The initial condition for the simulation is set to the exact solution to the Stokes’ problem after Fig. 2. Setup for the one-dimensional Stokes’ problem. 2128 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 a finite time of t0 ¼ 0:1 has elapsed in order to avoid the temporal discontinuity due to the impulsive start from interfering with the convergence study. Simulations are performed in a periodic computational domain in both x- and y-directions with uniform grid discretization. The top and bottom boundaries are placed far enough to avoid periodicity from interfering with the velocity profile near the translating plate. Spatial and temporal convergence is analyzed in terms of the L1 and L2 norms pffiffiffiffiffiof ffi the horizontal velocity error, ej ¼ uðy j Þ  uj , over the domain y j 2 ½0; 1 (in non-dimensional length: y j = mt0 2 ½0; 3:162). Fig. 3a assesses the temporal L1 error for various sizes of non-dimensional time steps, mDt=Dy 2 . The error was computed by comparing the solution to a temporally refined reference solution at fixed grid resolution to isolate the spatial discretization error. We calculate the error at t ¼ 0:11 with Dy ¼ 102 . The three convergence curves on the plot result from the use of different orders of expansion N for BN ( A1). Note that the splitting error from Eq. (24) is larger in magnitude than the underlying second-order error resulting from the time integration schemes. Hence this splitting error directly influences the temporal accuracy for the range of Dt considered. As discussed in Perot [30], the splitting error cannot be absorbed by k because LM 1 and Q are not commutative even for a periodic domain. Next we perform simulations with a very fine time step (Dt ¼ 106 ) and compare the results to the exact solution at t ¼ 0:101 for varying Dy. The velocity profile in the vicinity of the plate is influenced by the regularization of the Dirac delta function. This alters the velocity derivative at the immersed boundary causing the first-order accuracy of the L1 norm as shown in Fig. 3b. Fortunately, this smearing effect is dominant only in close proximity of the plate and the underlying second-order convergence is achieved in the L2 sense. 5.2. Flow inside two concentric cylinders For a two-dimensional test problem, we simulate flow between two concentric hollow cylinders with radii r1 ¼ 1=2 and r2 ¼ 1 as well as the flow inside the smaller cylinder as shown in Fig. 4. The outer cylinder is held stationary while the inner cylinder is rotated with angular velocity X,   uh ðr1 Þ t  0:2 X¼ ¼ 1 þ tanh ; ð41Þ r1 0:05 moving the initially quiescent fluid at t = 0. We take a periodic computational domain of size ½1:05; 1:05 ½1:05; 1:05 with uniform spatial resolution and compute the azimuthal velocity error, ej ¼ uh ðrj Þ  uh;j , over rj 2 ½0; r2  (including flow inside the inner cylinder) reporting the L1 and L2 norms. We study the impact of the splitting error from Eq. (24) on the temporal convergence by comparing our results to a reference solution obtained with a very fine time step, Dt ¼ 5 106 , and spatial resolution, Dx ¼ Dy ¼ 2:1 102 . The spatial resolution is kept constant and viscosity is set to m = 1. Fig. 5a shows that the order of expansion N for A1 again influences the behavior of convergence in a fashion similar to the onedimensional case. As it can be seen from the N = 3 case, the second-order time integration error starts to affect Fig. 3. Error norms from the one-dimensional Stokes’ problem. (a) Temporal L1 norm errors with different orders of expansion, N, for A1; N = 1: s, N = 2: h, and N = 3: n. (b) The L1: s and L2: h spatial velocity error norms. K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2129 Fig. 4. Setup for the problem of two concentric cylinders (inner cylinder rotates with angular velocity X). Fig. 5. Error norms from the problem of two concentric cylinders. (a) Temporal L1 norm errors with different orders of expansion, N, for A1; N = 1: s, N = 2: h, and N = 3: n. (b) The L1: s and L2: h spatial velocity error norms. (c) The L1: s and L2: h spatial pressure error norms. the total error at the smallest shown time step. Based on both the one- and two-dimensional test problems, we recommend the use of third-order expansion N for practical problems. There also is an advantage in choosing N = 3 for achieving positive-definiteness of the modified Poisson equation with larger choice of Dt [30]. Next we consider the spatial accuracy of our method at steady-state by comparing our results to the exact solution. The viscosity is reduced to m ¼ 0:01 to use a fine Dx while satisfying mDt=Dx2 K 1 to keep BN positivedefinite. Fig. 5b shows the rate of decay for the spatial errors to be 1 and about 1.5 in the L1 and the L2 sense. Although the first-order convergence is expected from the use of discrete delta functions, further investigation is required to explain why second-order accuracy from the underlying spatial discretization cannot be achieved in an L2 measure. The spatial accuracy of the pressure is also studied by comparing the current solution to the exact solution at steady-state. Because the pressure based on the current scheme only solves up to a constant (since we pin the pressure to remove the zero eigenvalue), we compare the solutions by matching the pressure at r = 0 for all cases and compute the error norms along the x-axis from 0 to r1. The infinity and L2 error norms are plotted against the grid size in Fig. 5c for the same problem considered in assessing the spatial accuracy of velocity. As expected, the spatial accuracy follows the same trend as the velocity shown in Fig. 5b. Due to the presence of the discrete Delta function along the immersed boundary, the pressure distribution is affected limiting the spatial accuracy to orders of one and about 1.5 for the infinity and L2 norms, respectively. 5.3. Flow over a stationary cylinder We consider flow over a circular cylinder as another test problem because the dimensions of the recirculation zone and the force on the cylinder at various Reynolds numbers are readily available from previous exper- 2130 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 imental and numerical studies. For the numerical studies, we list results from the IBM of Lai and Peskin [22] and the IIM of Linnick and Fasel [24] among others when the data are available. Our two-dimensional simulations are performed by introducing a cylinder of diameter d = 1 in a large computational domain D with initially uniform flow, u ¼ u1 ¼ 1. Reynolds numbers of Re ¼ u1 d=m ¼ 20, 40, and 200 are chosen for validating the current method at steady-state and periodic vortex shedding conditions (m is the kinematic viscosity). The computational domain is discretized non-uniformly in both x- and y-directions, while the grid spacing is kept uniform with its finest size (Dxmin ) in the vicinity of the cylinder. Table 1 summarizes the parameters used in the simulations, where nx and ny are the number of cells in the x- and y-directions and nB is the number of Lagrangian points on the surface of the cylinder with Ds Dxmin ¼ Dy min . Computations are performed with different sizes of D to ensure that the boundary conditions along oD do not influence our solution. Left (inflow) and lateral boundary conditions along oD are set to uniform flow of ðu; vÞ ¼ ðu1 ; 0Þ and are placed far away from the cylinder. At the outlet, the convective boundary condition (ou=ot þ u1 ou=ox ¼ 0) is applied to allow vorticity to exit the domain freely. Various spatial and temporal resolutions are chosen to ensure that reliable solutions are obtained. We record the maximum CFL number (CFLmax ¼ umax Dt=Dxmin ) in Table 1 from cases of Re = 40 and Re = 200. Note that the current method yields a stable solution even with CFLmax ¼ 0:81. For comparison, we compute the force on the body applied by the flow in terms of the drag and lift coefficients: C D  F x = 12 qu21 d and C L  F y = 12 qu21 d, respectively, where qu21 d ¼ 1. The force on the cylinder, F, can be obtained simply by   Z Z X F x ðtÞ FðtÞ ¼ H i;k fk DxDy fðnðs; tÞÞdðnðs; tÞ  xÞ ds dx  ¼ ð42Þ F y ðtÞ x s i using the regularization operator and the boundary forcing function. Summation over i is implied to take place separately for each direction of the force vector. First, simulations are performed for Re = 20 and 40 to validate the steady-state characteristics. The resulting wake dimensions and drag coefficients are compared to values reported in the literature. The size of the wake is characterized by l, a, b, and h (appropriately non-dimensionalized by the diameter) defined in Fig. 6 following the notation used in Coutanceau and Bouard [12]. The parameters, l, a, and b represent the length of the recirculation zone, distance from the cylinder to the center of the wake vortex, and the gap between the centers of the wake vortices, respectively. The separation angle is denoted by h measured from the x-axis. The steady-state vorticity contours and streamlines from Case B are shown in Fig. 7 for Re = 20 and 40. The flow profiles are in close agreement with those reported in the literature. The wake properties from Cases A and B are compared against previous experimental and numerical studies in Table 2 and are also found to be in accord. Next, we consider flow over a cylinder at a Reynolds number of 200 to reproduce periodic vortex shedding. A short time after simulations are initiated from uniform flow, a perturbation in a form of an asymmetric body force is added to trigger the shedding instability. Numerical results replicate the periodic shedding of vortices to form the Kármán vortex street as shown in the vorticity contour of Fig. 8. The resulting lift and drag coefficients and the Strouhal number, St  fs d=u1 , where fs is the shedding frequency, are compared to previous studies in Table 3. Results obtained from Cases B–D are found to be in good agreement with previous findings. Table 1 Parameters for spatial and temporal discretization used in the simulations nx Case Case Case Case A B C D 150 300 300 300 ny 150 300 300 300 D ½30; 30 ½30; 30 ½15; 45 ½10; 10 ½30; 30 ½30; 30 ½30; 30 ½30; 30 Dxmin Dt CFLmax nB 0.04 0.02 0.0333 0.0333 0.005 0.005 0.0125 0.0125 0.22* 0.46* 0.81 0.75 78 157 94 94 The maximum CFL numbers are reported from Re = 40 (*) and Re = 200 ( ) cases. K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2131 a θ b l Fig. 6. Definition of the characteristic dimensions of the wake structure. Fig. 7. Vorticity contours (top) for steady-state flow over a cylinder, where contour levels are set from 3 to 3 in increments of 0.4, and corresponding streamlines (bottom). For left and right plots, Re = 20 and 40, respectively. Results from Case D compared to Cases B and C suggest that the placement of the outflow boundary is not too critical. As a pair of positive and negative vortices convect downstream, their effect on the cylinder become less important since their far-field induced velocity would appear to cancel. On the other hand, we have observed pronounced interference from the lateral boundary conditions when the height of the computational domain is shortened. 2132 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 Table 2 Comparison of experimental and numerical studies of steady-state wake dimensions and drag coefficient from flow over a cylinder for Re = 20 and 40 l=d a=d b=d h CD Re ¼ 20 Coutanceau and Bouard Tritton [38]* Dennis and Chang [14] Linnick and Fasel [24] Present (Case A) Present (Case B) [12]* 0.93 – 0.94 0.93 0.97 0.94 0.33 – – 0.36 0.39 0.37 0.46 – – 0.43 0.43 0.43 45.0 – 43.7 43.5 44.1 43.3 – 2.09 2.05 2.06 2.07 2.06 Re ¼ 40 Coutanceau and Bouard [12]* Tritton [38]* Dennis and Chang [14] Linnick and Fasel [24] Present (Case A) Present (Case B) 2.13 – 2.35 2.28 2.33 2.30 0.76 – – 0.72 0.75 0.73 0.59 – – 0.60 0.60 0.60 53.8 – 53.8 53.6 54.1 53.7 – 1.59 1.52 1.54 1.55 1.54 Experimental studies are listed with (*). Fig. 8. Snapshot of the vorticity field with contour levels from 3 to 3 in increments of 0.4 for Re = 200. Table 3 Comparison of Strouhal number and coefficients of drag and lift for flow over cylinder from experimental and numerical studies at Re = 200 Re ¼ 200 Belov et al. [2] Liu et al. [25] Lai and Peskin [22] Roshko [35]* Linnick and Fasel [24] Present (Case B) Present (Case C) Present (Case D) St CD CL 0.193 0.192 0.190 0.19 0.197 0.196 0.195 0.197 1.19 ± 0.042 1.31 ± 0.049 – – 1.34 ± 0.044 1.35 ± 0.048 1.34 ± 0.047 1.36 ± 0.043 ±0.64 ±0.69 – – ±0.69 ±0.68 ±0.68 ±0.69 Experimental study is listed with (*). 5.4. Flow around a moving cylinder As our last test problem, we simulate flow around a circular cylinder in impulsive translation to validate the present method for moving bodies. The simulation is performed by moving the Lagrangian body points at K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2133 each time step. As these points shift their positions in time, the regularization and interpolation operators are updated according to Eq. (28). We initially position the cylinder with unit diameter (d = 1) at the origin and impulsively set it into motion to the left with a constant velocity of u0 ¼ 1. Results are presented for Reynolds numbers of Re ¼ ju0 jd=m ¼ 40 and 200. The computational domain D is taken to be ½16:5; 13:5 ½15; 15 with no-slip boundary condition applied along oD. Non-uniform grid is used with uniform grid in the near field having a resolution of Dxmin ¼ 0:02 resulting in a grid size of 425 250. A constant time step of Dt ¼ Dxmin =2 is chosen such that the maximum CFL numbers are limited to 0.98 and 0.81, respectively for Re = 40 and 200 during the simulation from a non-dimensional time of t  ju0 jt=d ¼ 0 to 3.5. Quiescent flow is used for the initial condition. We present snapshots of the flow field at non-dimensional time of t ¼ 1, 2.5, and 3.5 in Fig. 9. Left and right figures illustrate the vorticity field for Re = 40 and 200, respectively. The flow fields are in agreement with those in [13,21] for Re = 40. For Re = 200, the flow exhibits a generation of stronger vortex pair in the wake of Fig. 9. Snapshots of the vorticity field around an impulsively moving circular cylinder for Re = 40 and 200 at non-dimensional time of t ¼ 1, 2.5, and 3.5. Contour levels from 3 to 3 in increments of 0.4 are chosen. 2134 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 Fig. 10. History of the drag coefficient of the body for Re = 40 and 200 (—) compared with numerical solutions from Koumoutsakos and Leonard [21] (Re = 40, - - -) and Cottet et al. [11] (Re = 200, – - –) and analytical solution by Bar-Lev and Yang [1] (  ) valid for early time. the cylinder. In the two cases, the solutions are resolved well even near the boundary and the difference in the effect of viscous diffusion is nicely captured. The drag coefficients for the two cases are also computed by Eq. (42) during the simulation and are plotted in Fig. 10. Computational results based on vortex methods from [21,11] along with the analytical series solution [1] valid for early time are superposed onpthe ffiffiffiffi current results. The current scheme reveals the singular behavior of the drag at the start up time ðOð1= t ÞÞ experienced by the cylinder due to the impulsive motion [1]. Our drag coefficients are about 4–5% larger than those from the vortex method. Additional simulations were performed with smaller grid spacings and larger computational domains. However, there were no noticable changes in our solutions to account for the difference. We also measure the length of the recirculation zone, previously defined as l=d in Fig. 6, in the frame of reference of the cylinder ðu  u0 ; vÞ for validation over time. In Fig. 11, these lengths are compared with the reported curves from a numerical study of [10] and experimental findings of [13] and are found to be in excellent agreement shown by the overlaps for both Reynolds numbers. Fig. 11. Length of the recirculation zone, l=d, in the frame of reference of the moving cylinder as a function time, t*, for (a) Re = 40 and (b) Re = 200 compared with previous studies. Present results: —; experimental measurements of Coutanceau and Bouard [13] (Re = 40, - - -); and numerical study of Collins and Dennis [10] (Re = 200, – - –). K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2135 6. Conclusion We presented a new formulation of the immersed boundary method that is algebraically identical to the traditional fractional step algorithm. The current method introduces regularization and interpolation operators and regards both the pressure and boundary force as Lagrange multipliers required to satisfy the kinematic constraints of divergence-free and no-slip. The no-slip condition along the immersed surface is satisfied by a projection in a manner analogous to the removal of non-divergence-free component of the velocity field in the classical projection methods. The overall method is constructed to preserve symmetry and positive-definiteness to efficiently solve for the flow field. The boundary force is determined implicitly without any constitutive relations for the rigid body formulation. This in turn allowed us to use a CFL number as high as 1 in our simulations. Deforming bodies whose motion is known a priori can also be treated. The current scheme is numerically found to be third-order temporally accurate for most practical sizes of time steps and secondorder accurate in time for small time steps. The spatial accuracy is observed to be better than first-order in the L2 norm for one- and two-dimensional problems. Results from simulations of flow over both stationary and moving cylinders show excellent agreement with previous experimental and numerical studies. Acknowledgments The authors thank Prof. Blair Perot for the enlightening discussions on the fractional step method. This work was supported by the United States Air Force Office of Scientific Research (AFOSR/MURI FA955005-1-0369). Appendix. Derivation of Eq. (13) This appendix describes the details on how the incompressible Navier–Stokes equations Eqs. (10)–(12) are discretized on a finite-volume staggered mesh to reach the form of Eq. (13). A two-dimensional case is presented, although an extension to three dimensions is straightforward. Here, the underlying spatial discretization is taken to be a non-uniform Cartesian staggered mesh ðxi ; y j Þ with the immersed surface represented by a set of Lagrangian points ðnk ; gk Þ as shown previously in Fig. 1. Readers can also consult the works of [30,6] for details on the fractional step method for staggered grid formulation. The incompressible Navier–Stokes equations with a boundary forcing function, Eqs. (10)–(12), can be discretized with the AB2 and CN methods for the convective and viscous terms, respectively: unþ1  un 3 ^ n 1 ^ n1 ^1þH ^ þ 1 Lðu ^ nþ1 þ un Þ þ bc ^f; ðu Þ ¼ G/ þ N ðu Þ  N 2 2 2 Dt ^ nþ1 ¼ bc2 ; Du ^ nþ1 ¼ unþ1 ; Eu B nþ1 ð43Þ ð44Þ ð45Þ where u , /, and f are the discrete velocity, pressure, and boundary force. We order the discrete velocity and force vectors, ðu; vÞT and ðfx ; fy ÞT , respectively. The spatial operators introduced above are listed side-by-side with their continuous analog in Table 4. The discrete Laplacian and divergence operators generate inhomo^ 1 and bc2 resulting from the boundary conditions along oD. Note that bc ^ 1 depends on time geneous terms bc ^ and H ^ are levels n and n + 1 (for CN method). For Eq. (44), bc2 is a function of time level n + 1. Details on E ^ or more precisely the Lagrangian over D, operator E provided in Section 3.2. If oB moves with velocity unþ1 B points are functions of time level n + 1. As stated in [30], a staggered grid formulation with velocity boundary conditions requires no pressure boundary condition. Operators and vectors with hats are later transformed with a diagonal matrices for scaling purposes. Collecting the unknowns on the left-hand side, Eqs. (43)–(45) can be written as 2 30 nþ1 1 0 n 1 0 1 ^1 ^ G ^ H ^ ^r u A bc 6^ 7B C B C C B ð46Þ 4D 0 0 5@ / A ¼ @ 0 A þ @ bc2 A; nþ1 ^ u f E 0 0 0 B 2136 K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 Table 4 Nomenclature of the discrete operators and their continuous analogs Operator Discrete Continuous Divergence Gradient Interpolation ^ D ^ G ^ n Eu Regularization ^ Hf ^ L ^ ðun Þ N r  ðÞ rðÞ R uðtn ; ðxÞÞdðx  nÞ dx Rx s fðnðs; tÞÞdðn  xÞ ds Laplacian Convection Re1 r2 ðÞ r  ðuðtn Þuðtn ÞÞ where   ^  1 I  1L ^ and ^rn  1 I þ 1 L ^ un  3 N ^ ðun Þ þ 1 N ^ ðun1 Þ: ð47Þ A Dt 2 Dt 2 2 2 ^ can be made symmetric and positive-definite quite easily. Although we omit the details, A In order to solve the above system efficiently, symmetry among the sub-matrices are desired. First, let us make the gradient and divergence operators transpose of each other by a simple transformation. Both oper^ ators can be scaled appropriately so that the entries consist solely of 1 by introducing R and M: " #   1 ðDxi þ Dxi1 Þ 0 Dy j 0 ^  2 R and M ; ð48Þ 1 0 ðDy þ Dy j1 Þ 0 Dxi j 2 ^ T . For details on the construction of ^ 1 ¼ ðM ^ GÞ where the non-zero sub-matrices are diagonal. As a result DR ^ and D, ^ refer to [6]. Note that R transforms velocity unþ1 to velocity flux qnþ1  Runþ1 . Using these transforms G in Eq. (46), we find 2 30 nþ1 1 0 n 1 0 1 ^H ^ r bc1 q A G M C B C B C 6 7B ð49Þ 4 D 0 0 5@ / A ¼ @ 0 A þ @ bc2 A; nþ1 1 ^ uB 0 f ER 0 0 where ^ 1: ^ 1 ; G  M ^ ^ AR ^ G; ^ 1 ¼ GT ; rn  M^ ^ rn ; bc1  M ^ bc AM D  DR ð50Þ Also, for ease of discussion in Sections 2 and 3 we define the mass matrix and the transformed Laplacian by ^ 1 and L  M ^ LR ^ 1 such that A ¼ 1 M  1 L. We note that A is symmetric and positive-definite by M  MR Dt 2 construction. All steps presented up to this point in this Appendix are for the fractional step method and nothing special has been performed for the immersed boundary portion of our formulation. We recover Eq. (3) if we remove f and the no-slip constraint from Eq. (49). Finally, we re-define the interpolation and regularization operators by combining the diagonal matrices, R ^ and M: ^ 1 E  ER and ^H ^: H M Combining Eqs. (49) and (51), we obtain Eq. (13): 1 2 30 nþ1 1 0 n 1 0 r bc1 A G H q C B C B 6 7B C 0 5@ / A ¼ @ 0 A þ @ bc2 A; 4D 0 uBnþ1 E 0 0 0 f ð51Þ ð52Þ Before closing, we note again that G ¼ DT and A ¼ AT by construction. References [1] M. Bar-Lev, H.T. Yang, Initial flow field over an impulsively started circular cylinder, J. Fluid Mech. 72 (4) (1997) 625–647. [2] A. Belov, L. Martinelli, A. Jameson, A new implicit algorithm with multigrid for unsteady incompressible flow calculations, AIAA Paper 95-0049, 33rd Aerospace Sciences Meeting and Exhibit, Reno, NV, January 9–12, 1995. K. Taira, T. Colonius / Journal of Computational Physics 225 (2007) 2118–2137 2137 [3] R.P. Beyer, R.J. LeVeque, Analysis of a one-dimensional model for the immersed boundary method, SIAM J. Numer. Anal. 29 (2) (1992) 332–364. [4] P. Bochev, R.B. Lehoucq, On the finite element solution of the pure Neumann problem, SIAM Rev. 47 (1) (2005) 50–66. [5] D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2001) 464–499. [6] W. Chang, F. Giraldo, B. Perot, Analysis of an exact fractional step method, J. Comput. Phys. 180 (2002) 183–199. [7] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1) (1967) 12–26. [8] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [9] R. Codina, Pressure stability in fractional step finite element methods for incompressible flows, J. Comput. Phys. 170 (2001) 112–140. [10] W.M. Collins, S.C.R. Dennis, Flow past an impulsively started circular cylinder, J. Fluid Mech. 60 (1) (1973) 105–127. [11] G. Cottet, P. Koumoutsakos, M.L.O. Salihi, Vortex methods with spatially varying cores, J. Comput. Phys. 162 (2000) 164–185. [12] M. Coutanceau, R. Bouard, Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow, J. Fluid Mech. 79 (2) (1977) 231–256. [13] M. Coutanceau, R. Bouard, Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 2. Unsteady flow, J. Fluid Mech. 79 (2) (1977) 257–272. [14] S.C.R. Dennis, G. Chang, Numerical solutions for steady flow past a circular cylinder at Reynolds number up to 100, J. Fluid Mech. 42 (3) (1970) 471–489. [15] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (2000) 35–60. [16] R. Glowinski, T.W. Pan, J. Périaux, Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies, Comput. Methods Appl. Mech. Engrg. 151 (1998) 181–194. [17] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys. 105 (1993) 354–366. [18] J.-L. Guermond, L. Quartapelle, On stability and convergence of projection methods based on pressure Poisson equation, Int. J. Numer. Methods Fluids 26 (1998) 1039–1053. [19] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (1985) 308–323. [20] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, J. Comput. Phys. 171 (1) (2001) 132–150. [21] P. Koumoutsakos, A. Leonard, High-resolution simulations of the flow around an impulsively started cylinder using vortex methods, J. Fluid Mech. 296 (1995) 1–38. [22] M. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comput. Phys. 160 (2000) 705–719. [23] L. Lee, R.J. LeVeque, An immersed interface method for incompressible Navier–Stokes equations, SIAM J. Sci. Comput. 25 (3) (2003) 832–856. [24] M.N. Linnick, H.F. Fasel, A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains, J. Comput. Phys. 204 (2005) 157–192. [25] C. Liu, X. Zheng, C.H. Sung, Preconditioned multigrid methods for unsteady incompressible flows, J. Comput. Phys. 139 (1998) 35– 57. [26] G.I. Marchuk, Methods of Numerical Mathematics, Springer-Verlag, 1975. [27] R. Mittal, G. Iaccarino, Immersed boundary methods, Ann. Rev. Fluid Mech. 37 (2005) 39–261. [28] J. Mohd-Yusof, Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries, Center for Turbulence Research, Annual Research Briefs, 1997, pp. 317–327. [29] J. Nocedal, S.J. Wright, Numerical Optimization, Springer, 1999. [30] J.B. Perot, An analysis of the fractional step method, J. Comput. Phys. 108 (1993) 51–58. [31] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252–271. [32] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [33] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [34] A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (1999) 509–534. [35] A. Roshko, Report 1191: On the development of turbulent wakes from vortex streets, NACA, 1954. [36] J.C. Strikwerda, Y.S. Lee, The accuracy of the fractional step method, SIAM J. Numer. Anal. 37 (1) (1999) 37–47. [37] R. Témam, Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires (I), Arch. Rat. Mech. Anal. 32 (2) (1969) 135–153. [38] D.J. Tritton, Experiments on the flow past a circular cylinder at low Reynolds number, J. Fluid Mech. 6 (1959) 547–567. [39] N.N. Yanenko, The Method of Fractional Steps: The Solution of Problems of Mathematical Physics in Several Variables, SpringerVerlag, 1971.