A 3D immersed interface method for fluid–solid interaction

Available online at www.sciencedirect.com

Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086 www.elsevier.com/locate/cma

A 3D immersed interface method for fluid–solid interaction Sheng Xu a

a,*

, Z. Jane Wang

b

Department of Mathematics, Southern Methodist University, P.O. Box 750156, Dallas, TX 75275-0156, United States b Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14853, United States Received 20 February 2007; received in revised form 4 June 2007; accepted 17 June 2007 Available online 14 August 2007 This paper is written in honor of Professor Charles Peskin’s 60th birthday

Abstract In immersed interface methods, solids in a fluid are represented by singular forces in the Navier–Stokes equations, and flow jump conditions induced by the singular forces directly enter into numerical schemes. This paper focuses on the implementation of an immersed interface method for simulating fluid–solid interaction in 3D. The method employs the MAC scheme for the spatial discretization, the RK4 scheme for the time integration, and an FFT-based Poisson solver for the pressure Poisson equation. A fluid–solid interface is tracked by Lagrangian markers. Intersections of the interface with MAC grid lines identify finite difference stencils on which jump contributions to finite difference schemes are needed. To find the intersections and to interpolate jump conditions from the Lagrangian markers to the intersections, parametric triangulation of the interface is used. The velocity of the Lagrangian markers is interpolated directly from surrounding MAC grid nodes with interpolation schemes accounting for jump conditions. Numerical examples demonstrate that (1) the method has near second-order accuracy in the infinity norm for velocity, and the accuracy for pressure is between first and second order; (2) the method conserves the volume enclosed by a no-penetration boundary; and (3) the method can efficiently handle multiple moving solids with ease. Published by Elsevier B.V. Keywords: Immersed interface method; Immersed boundary method; Fluid–solid interaction; Singular force; Jump conditions

1. Introduction Immersed interface methods are offspring of the immersed boundary method. The original immersed boundary method was proposed by Peskin to simulate blood flow in the human heart [20,21]. It treats heart walls and heart valves as fiber-reinforced fluid. The immersed boundary method is therefore a mathematical formulation, in which the effects of solid boundaries are formulated as forces in the Navier–Stokes equations. The forces are determined from boundary configurations according to constitutive laws. They involve the form of the Dirac d function and are thus called singular forces. The immersed *

Corresponding author. Tel.: +1 214 768 2985; fax: +1 214 768 2355. E-mail addresses: [email protected] (S. Xu), [email protected] (Z.J. Wang). 0045-7825/$ - see front matter Published by Elsevier B.V. doi:10.1016/j.cma.2007.06.012

boundary method has been applied to a wide variety of problems, especially biological flows, as summarized in [23]. When applied to flow simulation, an immersed interface method shares the same mathematical formulation as the immersed boundary method, which reads ov 1 þ r  ðvvÞ ¼ rp þ Dv ot Re Z þ

fða1 ; a2 ; tÞdðx  Xða1 ; a2 ; tÞÞ da1 da2 ;

B

r  v ¼ 0;

ð1Þ ð2Þ

where v is velocity, p is pressure, t is time, Re is the Reynolds number, B is the boundary of a solid, dðx  Xða1 ; a2 ; tÞÞ is the 3D Dirac d function, x is Cartesian coordinates, Xða1 ; a2 ; tÞ is the coordinates of the boundary,

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

they were extended to 1D nonlinear parabolic equations [30], Poisson equations with Neumann boundary conditions [8], elliptic equations with variable coefficients [31,6,2], and the 2D Navier–Stokes equations [17,14,19, 33,13]. These various methods are summarized in the recent book by Li and Ito [18]. To extend immersed interface methods to the 3D Navier–Stokes equations, necessary jump conditions have been systematically derived by Xu and Wang [32]. A list of these jump conditions is given in Section 3. The incorporation of jump conditions into finite difference schemes is based on the following generalized Taylor expansion [32]:

and f is the density of a singular force in the parameter space which is formed by two Lagrangian parameters a1 and a2 parameterizing the boundary B as shown in Fig. 1. In the above formulation, only one boundary, the boundary B, is considered. This formulation is used hereafter for the presentation of this paper. If multiple boundaries are considered, they can be easily included in the same manner. In the immersed boundary method, the boundary of an immersed solid is tracked by Lagrangian markers that are convected by a fluid. Numerically, the communication between the solid and the fluid is obtained by spreading the singular forces from the Lagrangian markers to nearby Cartesian grid nodes and interpolating the velocity from nearby Cartesian grid nodes to the Lagrangian markers with the use of discrete Dirac d functions. Many research efforts have been devoted to analyze and improve the accuracy, stability, conservation, and robustness of the immersed boundary method [3,28,22,25,27,5,12,34]. Motivated to improve the accuracy of the immersed boundary method from first order to second order, LeVeque and Li [15,16] proposed immersed interface methods. To avoid the use of discrete Dirac d functions, an immersed interface method directly incorporates singularity-induced jump conditions of flow quantities into finite difference schemes, which gives it second-order or higher accuracy, sharp fluid– solid interfaces, and very good conservation of mass enclosed by no-penetration boundaries. Immersed interface methods were initially proposed for elliptic equations [15] and the Stokes equations [16]. Later,

gðs mþ1 Þ ¼

n

z

α2

2 gðs Þ  gðsþ dgðs 1 X ½gðnÞ ðnÞ i1 Þ i Þ ¼ iþ1 ðsi1  nÞn þ ds 2h n¼0 n! 2h ! 2 X ½gðnÞ ðgÞ n  ðsiþ1  gÞ þ Oðh2 Þ; n! n¼0

ð4Þ

3 þ gðs d2 gðs 1 X ½gðnÞ ðnÞ iþ1 Þ  2gðsi Þ þ gðsi1 Þ i Þ ¼  ðsi1  nÞn ds2 n! h2 h2 n¼0 ! 3 X ½gðnÞ ðgÞ n ð5Þ þ ðsiþ1  gÞ þ Oðh2 Þ: n! n¼0

An interpolation scheme also needs to account for jump conditions if its interpolation stencil contains discontinuities. The following second-order interpolation scheme applies to the case shown in Fig. 2b    gðsþ 1 ogðnÞ i1 Þ þ gðsiþ1 Þ 2 gðsi Þ ¼ þ Oðh Þ þ ðsi1  nÞ 2 os 2   1 ogðgÞ  ð6Þ ðsiþ1  gÞ: 2 os

α1

X y

x

The jump conditions derived in [32] have been employed in an immersed interface method to simulate the interaction

Fig. 1. A parametrized boundary in a Cartesian coordinate system.

a

ð3Þ

where g(s) is a non-smooth and discontinuous function as shown in Fig. 2a, and ½gðnÞ ðsl Þ denotes jump conditions ðnÞ  along the s-axis, i.e. ½gðnÞ ðsl Þ ¼ gðnÞ ðsþ l Þ  g ðsl Þ. Secondorder central finite difference schemes with discontinuities at n and g on its stencil shown in Fig. 2b can be modified as follows to keep their second-order accuracy:

τ

B

1 m X 1 X X gðnÞ ðsþ ½gðnÞ ðsl Þ 0Þ ðsmþ1  s0 Þ þ n! n! n¼0 l¼1 n¼0

 ðsmþ1  sl Þ ;

n

b

2069

b

g(s) +

-

-

+

+ + s0 s1

s2

......

+ - +

......

sm–1sm

h s

s i–1

s3

h

ξ

si η

s i+1

sm+1 s

Fig. 2. Examples for generalized Taylor expansion and finite differences: (a) a non-smooth and discontinuous function, (b) a finite difference stencil with discontinuities.

2070

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

of a fluid with moving boundaries in 2D [33,1]. Simulation results indicate that the 2D immersed interface method (1) achieves near second-order accuracy in the infinity norm for both velocity and pressure, (2) introduces relatively insignificant cost with the addition of a solid in a simulation, and (3) conserves volumes enclosed by non-penetration boundaries. In this paper, the derived jump conditions are used in an immersed interface method to simulate fluid–solid interaction in 3D. Compared with the existing 3D immersed boundary method, the current method has the improvements on spatial accuracy and resolution. It achieves near second-order accuracy in the infinity norm for the velocity, the accuracy for the pressure is between first and second order, and it does not smear sharp fluid–solid interfaces. Compared with body-fitted grid methods, the current method has the advantage in efficiency for moving boundary problems. Because of the use of a fixed Cartesian grid for fluids and Lagrangian markers for moving boundaries, the method does not need costly 3D grid regeneration, which is required in body-fitted grid methods. As shown in Section 6, the cost count of the current method in each time step is OðN ln N Þ þ OðM ln MÞ þ OðN Þ þ OðMÞ, where N is the total number of Cartesian grid nodes and M is the total number of Lagrangian markers. This paper is organized as follows. In Section 2, an overview of the method is given, which summarizes the major components needed by the method. Each major component is then presented in following sections. In Section 3, linear systems to determine jump conditions are listed. These jump conditions are the necessary ones to be incorporated into finite difference schemes. In Section 4, parametric triangulation of an interface is introduced. The parametric triangulation is used to identify finite difference stencils which pass across a fluid–solid interface. In Section 5, the interpolation of the velocity on staggered grid nodes and Lagrangian markers is presented. In Section 6, the major procedures of the current method are listed along with their cost counts. In Section 7, numerical examples are given to demonstrate the accuracy, conservation, and efficiency of the method. Last, Section 8 concludes the paper.

Terms with the divergence D are kept in Eq. (7) to better is discretized enforce the divergence-free condition, and oD ot by assuming D = 0 at the next time level. The current method solves the momentum equation, Eq. (1), and the pressure Poisson equation, Eq. (7), using the MAC scheme for the spatial discretization the RK4 scheme for the temporal integration, and an FFT-based Poisson solver for the pressure Poisson equation. A MAC grid is a staggered Cartesian grid, on which the pressure p and the velocity components u, v, and w are arranged as in Fig. 3. Define the central finite difference operators dx ; dy ; dz ; dxx ; dyy , and dzz as dx ðÞi;j;k ¼ dy ðÞi;j;k ¼

ðÞiþ1;j;k  ðÞi1;j;k 2

2

Dx ðÞi;jþ1;k  ðÞi;j1;k 2

2

Dy ðÞi;j;kþ1  ðÞi;j;k1

þ cx ðÞi;j;k ;

ð9Þ

þ cy ðÞi;j;k ;

ð10Þ

2 þ cz ðÞi;j;k ; Dz ðÞiþ1;j;k  2ðÞi;j;k þ ðÞi1;j;k dxx ðÞi;j;k ¼ þ cxx ðÞi;j;k ; Dx2 ðÞi;jþ1;k  2ðÞi;j;k þ ðÞi;j1;k þ cyy ðÞi;j;k ; dyy ðÞi;j;k ¼ Dy 2 ðÞi;j;kþ1  2ðÞi;j;k þ ðÞi;j;k1 þ czz ðÞi;j;k ; dzz ðÞi;j;k ¼ Dz2

dz ðÞi;j;k ¼

ð11Þ

2

ð12Þ ð13Þ ð14Þ

where Dx; Dy; Dz are the spatial steps as shown in Fig. 3, and cx ; cy ; cz ; cxx ; cyy , and czz are jump contributions. If the stencils of the above finite difference operators do not cross any fluid–solid interface, the jump contributions are zero, and usual central finite difference schemes are recovered. If the stencils of the above finite difference operators cross a fluid–solid interface, the jump contributions are non-zero, and they can be calculated according to Eqs. (4) and (5). With these central finite difference operators, the momentum equation, Eq. (1), and the pressure Poisson equation, Eq. (7), are spatially discretized as follows. The

2. Overview of the method

w

Δy

(i, j, k+1/2)

Taking the divergence of the momentum equation, Eq. (1), the pressure Poisson equation is obtained, which reads   oD 1 Dp ¼  þ r  ð2vDÞ  DD þ sp ot Re Z  þr fða1 ; a2 ; tÞdðx  Xða1 ; a2 ; tÞÞ da1 da2 ; ð7Þ B

where D = $ Æ v is the divergence of the velocity, and sp is   ou ov ou ov ou ow ou ow ov ow ov ow   þ  sp ¼ 2 þ : ox oy oy ox ox oz oz ox oy oz oz oy ð8Þ

v p

(i, j+1/2, k)

(i, j, k)

u

(i+1/2, j, k)

Δz z

y x Δx

Fig. 3. Arrangements of the velocity components and the pressure on a MAC grid.

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

spatially discretized equation for the velocity  momentum  component u at i þ 12 ; j; k can be written as ou 1 ¼ dx ðuuÞ  dy ðvuÞ  dz ðwuÞ  dx p þ ðdxx þ dyy þ dzz Þu; ot Re ð15Þ  1  in the operawhere the subscript i þ 2 ; j; k is neglected  1 tors. The1similar equations for v at i; j þ 2 ; k and w at i; j; k þ 2 can be obtained. The spatially discretized pressure Poisson equation at (i, j, k) can be written as oD  2ðdx ðuDÞ þ dy ðvDÞ ot 1 þ dz ðwDÞÞ þ ðdxx þ dyy þ dzz ÞD þ sp ; Re ð16Þ

ðdxx þ dyy þ dzz Þp ¼ 

where sp is calculated at (i, j, k) as sp ¼ 2ðdx udy v  dy udx v þ dx udz w  dz udx w þ dy vdz w  dz vdy wÞ: ð17Þ Eqs. (15)–(17) and the discretized equations for v and w needs the values of the velocity components at the grid nodes with subscripts listed in Table 1. They can be interpolated from uiþ12;j;k ; vi;jþ12;k , and wi;j;kþ12 . The interpolation schemes are given in Section 5. The RK4 temporal integration is used to march Eq. (15) and the spatially discretized equations for v and w in time. The reason to choose an explicit scheme is to be consistent with the explicit treatment of the motions of fluid–solid interfaces and the singular forces on the interfaces, which are functions of interface configurations. The reason to choose a high order scheme is to ensure numerical stability in the flow regime of moderate Reynolds numbers considered here. As pointed out by Johnston and Liu [9,10] and Weinan and Liu [7], high order explicit schemes are appropriate for flows of moderate to high Reynolds numbers, where viscous time step constraint is less restrictive than the convective one. The stability region of the RK4 scheme includes a portion of the imaginary axis, which ensures numerical stability of the background flow solver for the current flow regime. Eq. (16) is a discretized Poisson equation as oD is approxot imated by assuming D = 0 at the next time level, and it is solved in each substep of the RK4 temporal integration, as shown in [33]. Since the MAC grid is uniform in the current method, an FFT-based Poisson solver is adopted. The FFT-based Poisson solver can handle periodic boundary

Table 1 Subscripts of MAC grid nodes where velocity components need to be interpolated u v w

i, j, k i, j, k i, j, k

i þ 12 ; j þ 12 ; k i þ 12 ; j þ 12 ; k i þ 12 ; j; k þ 12

i þ 12 ; j; k þ 12 i; j þ 12 ; k þ 12 i; j þ 12 ; k þ 12

i; j þ 12 ; k i þ 12 ; j; k i þ 12 ; j; k

i; j; k þ 12 i; j; k þ 12 i; j þ 12 ; k

2071

conditions and inhomogeneous Dirichlet, Neumann, and mixed boundary conditions by using FFT, sine, cosine, and quarter wave transformations, respectively [24]. A summary of the major components required by the current method can be given below. • As indicated by Eqs. (4)–(6), necessary jump conditions are needed to obtain jump contributions in finite difference and interpolation schemes for Eqs. (15)–(17) and the spatially discretized equations for v and w. • Jump contributions are non-zero only if the stencils of finite difference and interpolation schemes cross fluid– solid interfaces. In order to distinguish these stencils, the intersections between MAC grid lines and the interfaces need to be identified, including the coordinates of the intersections and the necessary jump conditions at the intersections. • A fluid–solid interface follows the motion of the surrounding fluid. The fluid velocity is solved on MAC grid nodes, but the location of the interface is updated using the velocity of Lagrangian markers distributed on the interface. The velocity of Lagrangian markers needs to be interpolated from surrounding MAC grid nodes. 3. Computing the jump conditions The derivation of jump conditions listed in this section can be found in [32]. The formulas for the jump conditions in this section are different from those in [32], but they are equivalent mathematically. The formulas given in this section are more amenable to numerical implementation. The tangent vectors s and b, and the normal vector n shown in Fig. 1 appear in the expressions for the jump conditions below. They are defined as follows:   oX oX oY oZ s ¼ ðs1 ; s2 ; s3 Þ ¼ ¼ ; ; ; ð18Þ oa1 oa1 oa1 oa1   oX oX oY oZ b ¼ ðb1 ; b2 ; b3 Þ ¼ ¼ ; ; ; ð19Þ oa2 oa2 oa2 oa2 n ¼ ðn1 ; n2 ; n3 Þ ¼ s  b:

ð20Þ

The parameters a1 and a2 are chosen such that the vector n points to outside a solid. In addition, the following definitions are used: J ¼ knk; n n ¼ ; J f F¼ ; J F n ¼ F  n ;

ð21Þ ð22Þ ð23Þ ð24Þ 

Fs ¼ F  F n n ;

ð25Þ

where n* is the unit normal vector, and F is the density of singular force in the Cartesian space. In the numerical

2072

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

examples presented in Section 7, the force density F is calculated based on force models that relate the force density F to the configuration of the Lagrangian markers. The force models are given in each numerical example in Section 7. The jump conditions for the velocity and the pressure are ½v ¼ 0;

ð26Þ

½p ¼ F n :

ð27Þ

The jump conditions for the first derivatives of the velocity satisfy 2 3

0

ov ox

6 7 B 6 ov 7 B C 1 6 oy 7¼ 4 5 @ ov oz

0

1

0

C C; A

ð28Þ

1

s2

s3

b2

C b3 A :

n1

n2

n3

ð29Þ

The jump conditions for the first derivatives of the pressure satisfy 2 op 3 0 oF n 1 oa1

ox

6 7 B 6 7 B C 1 6 op 7¼B 4 oy 5 @ op oz

oF n oa2 of~ 1 oa1

þ

of~ 2 oa2

C C C; A

ð30Þ

where the contravariant components f~ 1 and f~ 2 in the parameter space are calculated by f~ 1 ¼ ðb  n Þ  F;

ð31Þ

f~ 2 ¼ ðn  sÞ  F:

ð32Þ



s1 s2 þ s2 s1 s1 s3 þ s3 s1 s2 s2 s2 s3 þ s3 s2 s3 s3

The jump conditions for the second derivatives of the velocity satisfy 2 2 3 1 0 2 o v o X o2 Y o2 Z 1 0 ox ox 0 oa oa oa oa oa oa 1 1 1 1 7 6 C B 1 1 6 o2 v 7 B C C B o2 X 2 2 7 B 6 o Y o Z C C B 0 6 ox oy 7 B oa2 oa2 oa2 oa2 oa2 oa2 C2 ov 3 C B 7 6 C ox C B 6 o2 v 7 B C B o2 X 0 6 ox oz 7 B o2 Y o2 Z C C6 ov 7 C B B 7 6 7 C  B oa1 oa2 oa1 oa2 oa1 oa2 C6 C26 2 7 ¼ B oy 7; C6 C B B oJ F 6 o v 7 B Re oa s C B on1 4 5 on3 C on2 1 6 oy oy 7 B C B oa1 oa1 oa1 C ov 7 B 6 C oz B 6 2 7 B Re oJ Fs C C C B on 6 ov 7 @ on3 C on2 oa2 A 1 B 6 oy oz 7 A @ oa oa oa 2 2 2 5 4 Re½rp o2 v 0 0 0 oz oz ð33Þ

1

C b1 b2 þ b2 b1 b1 b3 þ b3 b1 b2 b2 b2 b3 þ b3 b2 b3 b3 C C C s 1 b2 þ s 2 b1 s 1 b3 þ s 3 b1 s 2 b 2 s 2 b3 þ s 3 b2 s 3 b3 C C C; C s 1 n2 þ s 2 n1 s 1 n3 þ s 3 n1 s 2 n 2 s 2 n3 þ s 3 n2 s 3 n3 C C C b1 n2 þ b2 n1 b1 n3 þ b3 n1 b2 n2 b2 n3 þ b3 n2 b3 n3 C A 0 0 1 0 1 

on1 oa1

2 ; on oa1

3 ; on oa1



and ¼ and numerically according to

ReJ Fs

s1 B C 1 ¼ @ b1

C2 ¼ 0 s1 s1 B Bb b B 1 1 B Bs b B 1 1 B B B s 1 n1 B B B b1 n1 @ 1

on oa1

where ½ denotes jump conditions along the direction of n, and the coefficient matrix C1 is 0

where the coefficient matrix C2 is

on oa2

¼



on1 oa2

2 ; on oa2

3 ; on oa2

on o2 X o2 X ¼ bþs ; oa1 oa1 oa1 oa1 oa2 on o2 X o2 X ¼ bþs : oa2 oa1 oa2 oa2 oa2



ð34Þ are calculated

ð35Þ ð36Þ

The jump conditions for the second derivatives of the pressure satisfy 1 2 2 3 0 o2 F n 1 0 2 o p o X o2 Y o2 Z oa1 oa1 ox ox 7 B C 6 oa1 oa1 oa1 oa1 oa1 oa1 C B C 6 2 7 B o2 F n C B o2 X C 6op7 B 2 2 o Y o Z C C 6 ox oy 7 B B oa2 oa2 7 B C B oa2 oa2 oa2 oa2 oa2 oa2 C2 op 3 6 C C ox 6 2 7 B B o2 F n C B o2 X 6op7 B 2 Y o2 Z C C B oa oa oao oa C6 op 7 6 ox oz 7 B oa1 oa2 7 7 B C B 1 2 1 2 oa1 oa2 C6 C26 oy 7: C6 6 o2 p 7 ¼ B of~ of~ C  B 4 5 1 2 C 7 B o 6 B on on3 C on2 6 oy oy 7 B oa1 oa1 þ oa2 C B oa11 oa1 oa1 C op 7 B C B C oz 6 C B C 6 o2 p 7 B  ~ on3 C on2 7 B o of 1 of~2 C B on1 6 6 oy oz 7 B oa2 oa1 þ oa2 C @ oa2 oa2 oa2 A 5 @ A 4 o2 p 0 0 0 ½r  ðv  rvÞ oz oz ð37Þ The calculation of surface derivatives with respect to a1 and a2 is presented in Section 4. The non-singularity (except at the two poles of a spherical interface) of the coefficient matrices C1 and C2 is proved in [32]. The small linear systems given by Eqs. (28), (30), (33) and (37) can be solved analytically, as shown in Appendix . A flow quantity at a fixed point in space may have a jump in terms of time when a fluid–solid interface crosses the point, and a temporal jump condition can be related to a corresponding spatial jump condition [32,33]. In simulation of viscous flow, the incorporation of temporal jump conditions in temporal discretization has negligible effect on simulation results [33]. In the current method, temporal jump conditions are not included. 4. Parametric triangulation of an interface The current method considers immersed solids whose surfaces are smooth, orientable, and topologically equivalent to a sphere or a torus.

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

4.1. Interface parametrization The parametrization of an ellipsoidal shape is X s ¼ a sinða1 Þ cosða2 Þ; Y s ¼ b sinða1 Þ sinða2 Þ;

ð38Þ ð39Þ

Z s ¼ c cosða1 Þ;

ð40Þ

where the coordinates ðX s ; Y s ; Z s Þ are used to express the shape, a and b are the equatorial radii along the x- and y-axes, c is the polar radius along the z-axis, and a1 2 ½0; p and a2 2 ½0; 2p. The coordinates ðX ; Y ; ZÞ of an ellipsoidal fluid–solid interface are related to the shape coordinates ðX s ; Y s ; Z s Þ through translation and rotation transformations. The interface is spherical if a = b = c. A Lagrangian interface mesh as illustrated in Fig. 4a is generated with the following parameter discretization: Da1 þ m1 Da1 ; m1 ¼ 0; 1; . . . ; M 1 ; 2 ¼ m2 Da2 ; m2 ¼ 0; 1; . . . ; M 2 ;

a1m1 ¼

ð41Þ

a2m2

ð42Þ

where Da1 ¼ M 1pþ1 ; Da2 ¼ M2p2 , and M2 is an even integer. The interface is tracked by Lagrangian markers located at the nodes of the interface mesh. The integers M1 and M2 are chosen such that the maximum distance between two neighboring Lagrangian markers is about the spatial step of the background MAC grid. The jump conditions in Section 3 are calculated at the Lagrangian markers. The total

2073

number of the Lagrangian markers on the interface is M = M1(M2  1). No Lagrangian markers locate at the two poles corresponding to a1 = 0 and a1 = p. At the two poles, J = 0, and many equations in Section 3 have J in denominators of fractions. Realizing that the force density in the Cartesian space, F ¼ Jf , is finite at the two poles, it can be shown that these fractions are also finite at the two poles. The current choice of the Lagrangian interface mesh excludes the two poles, so the values of these fractions at the two poles are not needed. Otherwise numerical extrapolation has to be applied. The ellipsoidal interface is periodic in a2. The first derivatives with respect to a2 in Section 3 are calculated numerically using periodic cubic splines. The second derivatives with respect to a2 are calculated from their corresponding first derivatives also using periodic cubic splines. In order to use periodic cubic splines to calculate surface derivatives with respect to a1, a closed smooth curve on the interface for each a2m2 2 ½0; p is composed by two branches corresponding to values of a2m2 and p þ a2m2 , as demonstrated in Fig. 5a. So M2 has to be even. It can be shown that ~ ~ the functions JFs, f~ 1 , and ooaf 11 þ ooaf 22 are continuous at the two poles and smooth away from the two poles on this closed curve. All other functions to be differentiated with respect to a1 in Section 3 are smooth on the curve. On the branch corresponding to a2 ¼ a2m2 , a1 ¼ a1 ; on on ¼ n ; n oa1 oa1

ð43Þ ð44Þ

where a1 is defined in Fig. 5a, and n ¼ 0; 1; 2; . . . On the branch corresponding to a2 ¼ p þ a2m2 , a1 ¼ 2p  a1 ; 2n

Fig. 4. Interface parametrization and interface meshes: (a) a sphere, (b) a torus.

ð45Þ

2n

o o ¼ 2n ; oa2n oa 1 1

ð46Þ

o2nþ1 o2nþ1 ¼  : oa12nþ1 oa12nþ1

ð47Þ

Thus, periodic cubic splines with respect to a1 can be used n to calculate oao n with the above transformations. The cost 1

a

b α2

∗

α1

2π

p α2 π+α2m 2

α1

α1

π α2m 2

o 0

q π

2π α1

Fig. 5. Periodicity of an ellipsoidal interface: (a) formation of a closed curve past the two poles, (b) composition of periodic data.

2074

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

count to compute all the surface derivatives in Section 3 is OðMÞ. Because of the differentiation along the interface, it is important to maintain the smoothness of the interface for numerical stability. The current method employs Fourier filtering to smooth the interpolated velocity of Lagrangian markers on the interface. The interpolation approach is presented in Section 5. To filter in the direction of a2, original data defined on the parameter regions ‘‘o’’ and ‘‘p’’ in Fig. 5b are directly used because of their periodicity in a2. To filter in the direction of a1, data defined on the parameter regions ‘‘o’’ and ‘‘q’’ in Fig. 5b are used, where data on the region ‘‘q’’ are obtained by re-organize the data on the region ‘‘p’’, as illustrated in Fig. 5b. It is required that a1 start at Da2 1 and end at p  Da2 1 . The cost count of Fourier filtering is OðM ln MÞ. The surface derivatives in Section 3 can also be computed using Fourier transformations with more cost than periodic cubic splines. The parametrization of a torus is given by X s ¼ r cosða1 Þ;

ð48Þ

Y s ¼ ðR þ r sinða1 ÞÞ cosða2 Þ;

ð49Þ

Z s ¼ ðR þ r sinða1 ÞÞ sinða2 Þ;

ð50Þ

where R is the distance from the center of the torus tube to the center of the torus, r is the radius of the tube, and a1 2 ½0; 2p and a2 2 ½0; 2p. A Lagrangian interface mesh as illustrated in Fig. 4b is generated with the following parameter discretization: a1m1 ¼ m1 Da1 ;

m1 ¼ 0; 1; . . . ; M 1 ;

ð51Þ

a2m2 ¼ m2 Da2 ;

m2 ¼ 0; 1; . . . ; M 2 ;

ð52Þ

where Da1 ¼ M2p1 and Da2 ¼ M2p2 . The total number of the Lagrangian markers on the torus is M ¼ ðM 1  1ÞðM 2  1Þ. The torus has periodicity in the both directions of a1 and a2. So implementation of periodic cubic splines to calculate surface derivatives and Fourier filtering to smooth the torus is straightforward. 4.2. Interface triangulation As illustrated in Fig. 6, an interface are be approximated by small triangular patches formed from neighboring Lagrangian markers. In Fig. 6, two triangular patches, DP 1 P 2 P 4 and DP 2 P 3 P 4 , are formed from four neighboring Lagrangian markers (nodes of the interface mesh): P 1 ða1m1 ; a2m2 Þ; P 2 ða1m1 þ Da1 ; a2m2 Þ; P 3 ða1m1 þ Da1 ; a2m2 þ D a2 Þ, and P 4 ða1m1 ; a2m2 þ Da2 Þ. The parametrization of an ellipsoidal interface leaves two holes at the two poles. The two holes are covered by triangular patches as illustrated for the hole at a1 = 0 in Fig. 7. The intersections between an interface and MAC grid lines are found by projecting triangular patches along the x-, y-, and z-axes. Here the triangular patch DP1P2P4 in Fig. 8 is taken as an example. To find the intersections between this triangular patch and MAC grid lines parallel to the z-axis, DP1P2P4 is projected to the x–y plane along

α2m2+Δα 2 α2m2

P4

P1

P3

α1m1

α1m1+Δα 1

P2

Fig. 6. Parametric triangulation of an interface.

(0,2)

(0,M 2 - 1) (m 1 ,m2 )=(0,0) = (0,M 2 )

(0,1)

Fig. 7. Triangular patches covering the hole at a pole.

the z-axis to obtain the projection DQ1Q2Q4, which is contained inside a rectangle I II III IV. The rectangle is used to determine MAC grid lines which are parallel to the z-axis and may intersect DP1P2P4. If the projection QX(xI, yJ) of such a MAC grid line l, where I = i or i þ 12 and J = j or j þ 12, falls inside DQ1Q2Q4 (as the case in Fig. 8), on the edges of DQ1Q2Q4, or on the vertices of DQ1Q2Q4, the MAC grid line l intersects the triangular patch DP1P2P4 at the corresponding locations. The intersection is denoted as the point PX in Fig. 8. The x- and y-coordinates of the intersection PX are (xI, yJ). The z-coordinate and the necessary jump conditions at the intersection PX are interpolated from the three vertices P1, P2, and P4, where the values of the three Cartesian coordinates, the two Lagrangian parameters, and the necessary jump conditions are all known. If DP1P2P4 is not a triangular patch which covers the hole at a pole of an ellipsoidal interface, the following linear interpolation is used: gða1 ; a2 Þ ¼ cg0 þ cg1 a1 þ cg2 a2 ;

ð53Þ

where g can be a Cartesian coordinate of an interface or a jump condition across an interface, and cg0, cg1, and cg2 are constants. The constants cg0, cg1, and cg2 are determined from the vertices in the following linear system: gða1m1 ; a2m2 Þ ¼ cg0 þ cg1 a1m1 þ cg2 a2m2 ; gða1m1 þ Da1 ; a2m2 Þ ¼ cg0 þ cg1 ða1m1 þ Da1 Þ þ cg2 a2m2 ;

ð54Þ ð55Þ

gða1m1 ; a2m2 þ Da1 Þ ¼ cg0 þ cg1 a1m1 þ cg2 ða2m2 þ Da2 Þ:

ð56Þ

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

l

P4

a n*

n*

z

n*

Px

P1

l

b

l

z

z

2075

n* Fig. 9. Intersection on an edge: (a) recorded, (b) discarded.

P2

y IV

yj ( yj+1/2 )

Q1

r1

I

Q4

III

r4 Qx r2

a

b

l

z

II

x i (x i+1/2 )

z y

y

Q2

l

x

Fig. 8. Projection of a triangular patch for finding intersections between an interface and grid lines.

x

x

Fig. 10. Intersection on a vertex: (a) recorded, (b) discarded.

The values of the two parameters, a1 and a2, at the intersection PX are determined from xI ¼ cg0 þ cg1 a1 þ cg2 a2 ;

ð57Þ

y J ¼ cg0 þ cg1 a1 þ cg2 a2 :

ð58Þ

This interpolation is of second-order accuracy in terms of Da1 and Da2. If DP1P2P4 is a triangular patch which covers the hole at a pole of an ellipsoidal interface, the following area-based interpolation formula is used instead: gðP X Þ ¼

ag1 gðP 1 Þ þ ag2 gðP 2 Þ þ ag4 gðP 4 Þ ; ag1 þ ag2 þ ag4

ð59Þ

where ag1, ag2, and ag4 are net areas calculated as follows: ag1 ¼ ðr2  r4 Þ  ez ;

ð60Þ

ag2 ¼ ðr4  r1 Þ  ez ; ag4 ¼ ðr1  r2 Þ  ez ;

ð61Þ ð62Þ

where the vectors r1, r2, and r4 are defined in Fig. 8, and ez is the unit basis vector of the z-axis. After all the coordinates of the point PX and all the necessary jump conditions at the point PX are known, finite difference stencils which contain the point PX are identified, and jump contributions to usual finite differences on these stencils are calculated. The above considers the situation in which MAC grid lines are parallel to the z-axis. A similar consideration applies to the other two directions. Before ending this section, some special situations which need special care are described below. When an intersection falls on an edge (as in Fig. 9) or a vertex (as in Fig. 10) of a triangular patch, it is important to ensure that the intersection is not missed or repeated. First, the intersections falls on the edges represented by

dashed lines or the vertices represented by open circles in Fig. 6 are not counted to avoid repetition. Second, if a MAC grid line is aligned with a triangular patch (the area of its projection is zero), the MAC grid line is regarded either parallel or tangential to the patch, and no intersection is recorded. Third, when an intersection falls on an edge but not a vertex of a triangular patch, it is counted only if the patch is oriented with respect to the corresponding MAC grid line in the same way as the adjacent patch, as illustrated in Fig. 9, where the orientation of the patches in Fig. 9 is defined as the sign of the z-component of the normal direction n*. Fourth, when the intersection is on a vertex, it is kept only if the projection of the vertex falls inside the polygon formed by the projection of all the triangular patches that share this common vertex, as illustrated in Fig. 10. The equation of a projected edge can have different forms. For example the equation of the edge Q1Q2 in Fig. 8 can be written in the two forms as follows: ðy  Y 1 ÞðX 2  X 1 Þ ¼ ðx  X 1 ÞðY 2  Y 1 Þ; ðy  Y 2 ÞðX 1  X 2 Þ ¼ ðx  X 2 ÞðY 1  Y 2 Þ;

ð63Þ ð64Þ

where (X1, Y1) and (X2, Y2) are the coordinates of the points Q1 and Q2, respectively. When determining whether the point QX falls inside DQ1Q2Q4, on the edge Q1Q2, or inside the neighboring triangle sharing the same edge Q1Q2, the same form of the equation of the edge Q1Q2 has to be used for the two adjacent triangles. Otherwise, the intersection PX may be missed or repeated in counting due to rounding errors if the point QX falls on or is very close to the edge Q1Q2. This is a very trivial case, but it

2076

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

happens in practice. Finally, a crude check on intersection finding is the total number of intersections. It should be even. An intersection on a triangular patch may coincide with a MAC grid node. If this occurs, the intersection is regarded at either one side or the other of the grid node along the x-, y-, and z-axes. When computing jump contributions to finite difference or interpolation schemes along each axis, a consistent choice can be made for this axis by infinitesimally detaching the triangular patch away from the grid node either toward the normal direction of the triangular patch or opposite to it.

The staggered arrangement of the velocity components u, v and w and the pressure p, as illustrated in Fig. 3, necessitates the interpolation of the velocity components at the MAC grid nodes with subscripts listed in Table 1 from the defined velocity components uiþ12;j;k ; vi;jþ12;k , and wi;j;kþ12 . Define the interpolation operators ei, ej, ek as

ej ðÞi;j;k ¼ ek ðÞi;j;k ¼

ðÞiþ1;j;k þ ðÞi1;j;k 2

2

2 ðÞi;jþ1;k þ ðÞi;j1;k 2

2

2 ðÞi;j;kþ1 þ ðÞi;j;k1 2

2

2

þ ci ðÞi;j;k ;

ð65Þ

þ cj ðÞi;j;k ;

ð66Þ

þ ck ðÞi;j;k ;

ð67Þ

ð68Þ ð69Þ

uiþ12;j;kþ12 ¼ ek uiþ12;j;kþ12 ;

ð70Þ

a

N+

ui;j;kþ12 ¼ ek ui;j;kþ12 ;

ð72Þ

 where the normal derivative onov ¼ ReFs according to Eq. (28). Trilinear interpolation is used to interpolate the velocity of each supplemental points from surrounding MAC grid points via transitional points in three separate steps, as shown in Fig. 11b, where the Cartesian grid points are the vertices of the cell, the transitional points lie on the edges and the faces of the cell, and a supplemental point locates inside the cell. The order of each interpolation step is marked in Fig. 11b. The distance Dnffi * is chosen to be a little qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi larger than ðDxÞ2 þ ðDyÞ2 þ ðDzÞ2 such that all the interpolation cell in Fig. 11b is not cut by any interface. Thus a standard interpolation scheme can be applied in each interpolation step in Fig. 11b. Because of the staggered arrangement of the velocity components, the interpolation cell shown in Fig. 11b is different for the different velocity components. In Fig. 11b, if the cell is for the velocity component u, the transitional point II can be interpolated from the MAC grid points I and III as the following:

where ci, cj, and ck are jump contributions. If the stencils of the above interpolation operators do not cross any fluid– solid interface, the jump contributions are zero, and usual interpolation schemes are recovered. If the stencils of the above interpolation operators cross a fluid–solid interface, the jump contributions are non-zero, and they can be calculated according to Eq. (6). With these interpolation operators, the interpolation for the first row in Table 1 can be written as follows in the listing order: ui;j;k ¼ ei ui;j;k ; uiþ12;jþ12;k ¼ ej uiþ12;jþ12;k ;

ð71Þ

The interpolation for the second and the third rows in Table 1 can be written similarly. A fluid–solid interface moves with the fluid. To update the interface location, the velocity of Lagrangian markers at the interface is interpolated from surrounding fluid velocity. The current method takes the interpolation strategy illustrated in Fig. 11. As shown in Fig. 11a, the velocity of the Lagrangian marker L on the interface is interpolated from two supplemental points N+ and N along the normal direction n* at the marker, and the two supplemental points are at the different sides of the interface with the equal distance Dn* away from the interface. According to Eq. (6), the interpolation scheme at this step is   vðN þ Þ þ vðN  Þ 1 ovðLÞ 2  vðLÞ ¼ ð73Þ Dn þ OððDn Þ Þ; 2 2 on

5. Interpolation of the velocity

ei ðÞi;j;k ¼

ui;jþ12;k ¼ ej ui;jþ12;k ;

uðIIÞ ¼

zðIIIÞ  zðIIIÞ zðIIÞ  zðIÞ uðIÞ þ uðIIÞ þ OððDzÞ2 Þ: zðIIIÞ  zðIÞ zðIIIÞ  zðIÞ ð74Þ

III

b

n 4th

N +( N - )

L 1st

B N-

z

3rd

y x

Fig. 11. Interpolation of the velocity of a Lagrangian marker.

II 2nd I

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

More transit points and Cartesian grid nodes than in Fig. 11 can be used to achieve higher order interpolation as long as no interfaces cut interpolation cells. There are other options available to interpolate the velocity of Lagrangian markers from surrounding fluid velocity. The following two options are avoided in the current practice. In the first one, the velocity of a Lagrangian marker is interpolated from the three closest intersections using a area-based formula similar to Eq. (59), and the velocity of the intersections is previously interpolated from MAC grid nodes both inside and outside of the interface. It turns out this option only gives first-order accuracy in the infinity norm for the velocity. In the second option, the velocity of a Lagrangian marker is extrapolated from MAC grid nodes either inside or outside of the interface. Accuracy of high order can be achieved, but the method suffers from a small numerical stability region with relatively low Reynolds numbers.

6. Summary of the method The major procedures of the current method can now be summarized as follows. The computational cost associated with each step is also given, with N denoting the total number of MAC grid nodes for the pressure and M denoting the total number of Lagrangian markers. (Since the integers M1 and M2 are chosen such that the maximum distance between two neighboring Lagrangian markers is about the spatial step of the background MAC grid, the total number of intersections between MAC grid lines and fluid–solid interfaces is of order M too.) 1. Initializing the flow field and the interfaces ðOðN Þ þ OðMÞÞ; 2. calculating surface derivatives of geometric quantities ðOðMÞÞ; 3. modeling singular forces ðOðMÞÞ; 4. calculating surface derivatives of forcing quantities ðOðMÞÞ; 5. calculating jump conditions ðOðMÞÞ; 6. finding the intersections ðOðMÞÞ; 7. calculating jump contributions to finite difference and interpolation schemes ðOðMÞÞ; 8. interpolating and smoothing the velocity of the Lagrangian markers ðOðMÞ þ OðM ln MÞÞ; 9. updating the interface configurations ðOðMÞÞ; 10. solving the pressure field ðOðN ln N ÞÞ; 11. updating the velocity field ðOðN ÞÞ.

2077

7.1. Flow inside a rotating object This first example considers the steady flow inside an object which rotates with a constant angular velocity X around a unit vector R. The object is in the middle of a ½1; 1  ½1; 1  ½1; 1 cuboid. Rigid-wall boundary conditions are applied at the six sides of the cuboid. The analytical solution of the flow inside the rotating object is vðxÞ ¼ XR  x; 1 pðxÞ ¼ X2 kR  xk22 þ p0 ; 2

ð75Þ ð76Þ

where p0 is an arbitrary constant. To enforce the prescribed rotation, each Lagrangian marker on the object is connected to its prescribed position by a linear spring, and the density of the singular force from the spring model is given by the following equation: F ¼ K s ðXe  XÞ;

ð77Þ

where Ks is the spring stiffness, and Xe is the prescribed position of a Lagrangian marker. 7.1.1. Spatial convergence Spatial convergence of the method is indicated by the change of simulation errors with grid refinement. The simulation errors are based on the analytical solution, Eqs. (75) and (76). Table 2 presents the results of spatial convergence analysis for the steady flow inside a rotating sphere with the diameter equal to 1 at Re = 10. The angular velocity of the rotation = 1, and the rotation direction is pffiffi pis pffiffi ffiffi X given by R ¼ 33 ; 33 ; 33 . The value of the spring stiffness in the spring model is Ks = 1000. In Table 2, Nx, Ny, and Nz are the numbers of MAC cells for the pressure p along the x-, y-, and z-axes, respectively. The order of accuracy is calculated by the following formula: order ¼

lnðkecurrent k1 =keprevious k1 Þ ; lnðDcurrent =Dprevious Þ

ð78Þ

7. Numerical examples

where ecurrent and eprevious denote the errors at the current and the previous rows in Table 2, respectively, and Dcurrent and Dprevious denote the corresponding spatial resolutions. The results in Table 2 indicate that the accuracy for the three velocity components u, v, and w is near second-order, and the accuracy for the pressure p is between first and second order. Table 3 are the results for the steady flow inside a rotating torus with R = 0.5 and r = 0.25 at Re = 10. The angular velocity of the rotation is X = 1, and the rotation direction is R ¼ ð1; 0; 0Þ. The value of the spring stiffness is Ks = 1000. The results in Table 3 also indicate that the accuracy for the velocity components u, v and w is near second-order, and the accuracy for the pressure p is between first and second order.

In this section, numerical examples simulated by the current method are given to test its accuracy, conservation, and efficiency.

7.1.2. Effect of spring stiffness The effect of spring stiffness is investigated by simulating the steady flow inside a rotating sphere with different values

2078

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

Table 2 Spatial convergence analysis for the flow inside a rotating sphere Nx · Ny · Nz, M1 · M2 25 · 25 · 25, 33 · 33 · 33, 49 · 49 · 49, 65 · 65 · 65, 97 · 97 · 97,

24 · 48 32 · 64 48 · 96 64 · 128 96 · 192

keuk1

Order 2

kevk1

Order 2

3.51 · 10 2.23 · 102 1.12 · 102 6.84 · 103 3.29 · 103

1.63 1.74 1.75 1.83

kewk1

Order 2

3.45 · 10 2.21 · 102 1.14 · 102 6.88 · 103 3.32 · 103

Order 2

3.41 · 10 2.19 · 102 1.12 · 102 7.51 · 103 3.29 · 103

1.60 1.70 1.41 2.06

Order

kewk1

1.14 1.28 3.01 1.89

1.15 · 101 8.18 · 102 4.81 · 102 2.08 · 102 9.68 · 103

1.60 1.67 1.79 1.82

kepk1 2.33 · 10 1.15 · 102 6.83 · 103 5.40 · 103 2.94 · 103

2.45 1.28 0.82 1.52

Order

kepk1

Order

1.23 1.34 2.97 1.91

1.42 · 101 9.23 · 102 5.98 · 102 3.56 · 102 2.27 · 102

1.50 1.07 1.80 1.11

Table 3 Spatial convergence analysis for the flow inside a rotating torus Nx · Ny · Nz, M1 · M2

keuk1

25 · 25 · 25, 33 · 33 · 33, 49 · 49 · 49, 65 · 65 · 65, 97 · 97 · 97,

1.58 · 102 1.51 · 102 8.16 · 103 5.30 · 103 2.97 · 103

24 · 48 32 · 64 48 · 96 64 · 128 96 · 192

Order

kevk1

0.16 1.56 1.53 1.45

1.12 · 101 8.17 · 102 4.92 · 102 2.10 · 103 9.85 · 103

of spring stiffness. The spatial resolution of the simulation corresponds to Nx · Ny · Nz = 33 · 33 · 33 and M1 · M2 = 32 · 64. Table 4 lists the change of simulation errors against the spring stiffness. For this particular flow, the spring stiffness in the considered range has very small effect on the infinity norm of the simulation errors for both the velocity and the pressure. The errors for the positions of Lagrangian markers are plotted for Ks = 10 and Ks = 5000 in Fig. 12. The amplitudes of the position errors are much larger for much smaller spring stiffness, as expected. The spring model introduces a vibration time scale into the flow. The effect of this time scale on a time-dependent flow has been investigated in [33]. How to choose the values of the spring stiffness has also been discussed in [33]. 7.2. Flow induced by a relaxing balloon In the second example, a 3D pressurized balloon immersed in an incompressible fluid relaxes to its spherical equilibrium shape from the initial distortion. The initial velocity and the pressure are set zero, and the coupled motion of the balloon and the fluid is driven only by balloon tension. The simulation domain is a square box with dimensions ½0:8; 0:8  ½0:8; 0:8  ½0:8; 0:8. The initial shape of the balloon is an ellipsoid with a = 0.64, b = 0.4, and c = 0.25. The center coordinates of the ellipsoid are ð0; 0; 0Þ. The density of singular force is modeled by F ¼ EH n ;

ð79Þ

where H is the mean curvature of the balloon surface, and E = 0.2 is a constant. At equilibrium, the balloon should be pffiffiffiffiffiffiffi a sphere with radius ro ¼ 3 abc ¼ 0:4 and center coordinates ð0; 0; 0Þ, the velocity should be zero everywhere, and the pressure should be piecewise constants with a jump, ½p ¼  rEo ¼ 0:5, across the balloon surface. Shown in Fig. 13 is the simulated evolution of balloon shape at Re = 100. The spatial resolution of this simulation is Nx · Ny · Nz = 33 · 33 · 33 and M1 · M2 = 32 · 64, and the time step of this simulation is fixed with Dt = 0.01. At this Reynolds number, the balloon undergoes a couple of oscillations before it settles down to equilibrium, as indicated in Fig. 14a. In Fig. 14a, the distances ra, rb, and rc to the domain  center   ð0; 0; 0Þ from three Lagrangian points ðm1 ; m2 Þ ¼ M21 ; 0 ; M21 ; M42 , and (0, 0), are plotted against time at four different spatial simulation resolutions. The volume conservation of the balloon in the relaxation process is checked in Fig. 14b. The volume is well conserved during the process with its calculated values very close to the analytical one, 4pabc  0:2681. The volume errors during 3 the relaxation process are of the same amplitudes as their initial values (at the time t = 0) which are due to the approximation of the volume by a Riemann sum. The relative volume errors do not exceed 0.37% at the four resolutions. The equilibrium pressure at two slices shown in Fig. 15 indicates that it is piecewise constants with a jump across the balloon surface. The volume conservation is checked in Fig. 16 for Re = 1. The spatial resolution of this simulation is Nx · Ny · Nz = 33 · 33 · 33 and M1 · M2 = 32 · 64, and

Table 4 Effect of spring stiffness on the simulation accuracy for the flow inside a rotating sphere Ks

10

100

500

1000

2000

5000

keuk1 kevk1 kewk1 kepk1

2.35 · 102 2.35 · 102 2.39 · 102 1.26 · 102

2.25 · 102 2.26 · 102 2.22 · 102 1.61 · 102

2.23 · 102 2.25 · 102 2.20 · 102 1.28 · 102

2.23 · 102 2.21 · 102 2.19 · 102 1.15 · 102

2.23 · 102 2.60 · 102 2.19 · 102 1.98 · 102

2.25 · 102 2.25 · 102 2.19 · 102 1.44 · 102

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

2079

Fig. 12. Errors for the positions of Lagrangian markers: (a) Ks = 10, (b) Ks = 5000.

the time step Dt = 0.005 is also fixed. No oscillations are observed at this low Reynolds number, and it takes a very long time for the balloon to reach equilibrium. Fig. 17 plots the distances ra, rb, and rc, and the volume against time for Re = 10 at different CFL numbers. The convective and viscous CFL numbers CFLc and CFLl are defined as   umax vmax wmax CFLc ¼ Dt þ þ ; ð80Þ Dx Dy Dz ! Dt 1 1 1 þ þ ; ð81Þ CFLl ¼ Re ðDxÞ2 ðDyÞ2 ðDzÞ2 where umax, vmax, and wmax are the maximum velocity components in the flow field. Again, the volume of the balloon is well preserved in the relaxation process, with the relative errors less than 0.52% for the four values of the CFL numbers.

ple. The spring model given by Eq. (77) is used to calculate the density of singular force with Ks = 1000, 500, 200, and 100 for Re = 10, 20, 100, and 200, respectively. The sphere locates in a box domain. Its diameter equals to 1, and its center coordinates are (0, 0, 0). The spatial resolution of the simulation is given by Nx · Ny · Nz = 256 · 129 · 129 and M1 · M2 = 32 · 64. For Re = 10, 20, and 100, the domain sizes are ½4; 16  ½4; 4  ½4; 4 along the x-, y-, and z-axes, and the fixed time step Dt = 0.005 is used. For Re = 200, the domain sizes are ½2; 8  ½2; 2 ½2; 2, and the time step is controlled by CFLc = CFLl = 0.2. A free stream enters the domain in the direction of the x-axis. At the four sides of the domain, symmetric boundary conditions are used. At the domain outlet, the following outflow boundary conditions are used:

7.3. Flow past a stationary sphere

ov ¼ 0; ox op 1 o2 u ¼ : ox Re ox2

Flow past a stationary sphere at varying Reynolds numbers, Re = 10, 20, 100, and 200, is simulated in this exam-

At Re = 10, 20, and 100, the initial flow field is set uniform with u = 1 as the far field, and the sphere impulsively stops

Fig. 13. Evolution of the balloon shape at Re = 100.

ð82Þ ð83Þ

2080

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

a

b

0.65

r

0.2695

r

0.269

r

0.2685

a

0.6

b

0.55

0.27

c

0.5

0.268 0.2675

0.45

0.267

0.4

0.2665 0.35 0.266 0.3

0.2655

0.25 0

5

10

15

0.265 0

5

10

15

t

t

Fig. 14. Volume conservation at Re = 100 with different spatial resolutions: (a) temporal change of distances from three Lagrangian markers to the domain center, (b) temporal change of volume. Nx · Ny · Nz and M1 · M2: 97 · 97 · 97 and 96 · 192 for lines, 65 · 65 · 65 and 64 · 128 for lines with ‘‘+’’ marks, 49 · 49 · 49 and 48 · 96 for lines with ‘‘o’’ marks, and 33 · 33 · 33 and 32 · 64 for lines with ‘‘·’’ marks.

is used to ramp the far-field uniform flow from u = 0 to u = 1 to avoid impulsive stop of the sphere     tn1 tn  exp  ; n ¼ 1; 2; . . . ; un ¼ un þ exp  tc tc

0.8 0.4

z

ð84Þ 0 –0.4 –0.8 0.8 0.4

0.8 0.4

0

0

–0.4

y

–0.8 –0.8

–0.4

x

Fig. 15. Equilibrium pressure p at the two slices x = 0 and y = 0.

from the uniform velocity. At Re = 200, the initial flow field is set quiescent with u = 0, and the following formula

a 0.65

where the subscript n denotes a discrete time level, t0 = 0 corresponds to the initial time, and tc = 1 is a characteristic time of the ramping process. Correspondingly, the bound ary condition for u at the domain inlet is u ¼ 1  exp  ttnc instead of u = 1 at the other Reynolds numbers. In general, the fluid force G applied by a fluid to an object can be calculated by  þ  Z  1 ov G¼ pþ n þ dS Re on S    Z  Z 1 ov p n þ dS; ð85Þ ¼  f da1 da2 þ Re on S S

b 0.2672 ra

0.6

0.2671

rb

0.55

0.267

rc

0.5

0.2669 0.45 0.2668 0.4 0.2667

0.35

0.2666

0.3 0.25 0

5

10

15

20

t

25

30

35

40

0.2665 0

5

10

15

20

25

30

35

40

t

Fig. 16. Volume conservation at Re = 1: (a) temporal change of distances from three Lagrangian markers to the domain center, (b) temporal change of volume.

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

a 0.65

2081

b 0.2672 r

0.6

a

0.2671

r

b

0.55

rc

0.5

0.267

0.45

0.2669

0.4

0.2668

0.35 0.2667

0.3 0.25 0

1

2

3

4

5

6

7

8

0.2666 0

1

2

3

4

5

6

7

8

t

t

Fig. 17. Volume conservation at Re = 10 at different CFL numbers: (a) temporal change of distances from three Lagrangian markers to the domain center, (b) temporal change of volume. lines: CFLc = CFLl = 0.05, lines with ‘‘+’’ marks: CFLc = CFLl = 0.1, lines with ‘‘o’’ marks: CFLc = CFLl = 0.3, and lines with ‘‘·’’ marks: CFLc = CFLl = 0.6.

where the subscript ‘‘+’’ denotes the outer side of the object surface (fluid–solid interface) and the subscript ‘‘’’ deR notes the inner side. The term S f da1 da2 is just the resultant external force on the object generated by a force model (a collection of springs in this example). For a centrally symmetric object, its motion can be regarded as the superposition of translation and rotation at its geometric center, and the following relations apply: Z dUt ; ð86Þ ðp n ÞdS ¼ V dt S Z   ov dS ¼ 0; ð87Þ on S

6

Re=10

which is simply the application of Newton’s second law to the translational motion of the fluid contained inside the object. The sphere is stationary in the current example, so the calculation of the fluid force is given by Z G ¼  f da1 da2 : ð89Þ S

Fig. 18 shows the time history of the drag coefficient C1 for flow past the sphere. The force coefficients C1, C2, and C3 are defined as ðC 1 ; C 2 ; C 3 Þ ¼ C ¼ 2G=Af , where Af ¼ p4 is the projected frontal area for the sphere. A constant value of the drag coefficient is achieved for each of the Reynolds numbers as the flow becomes steady after the transit. Its value is compared with previous computational results [26,11] in Table 5, and the agreement is very good. Due to the spring model, oscillations in the drag coefficient are observed in the transit. With the ramping process de-

Re=100

Re=200

C1

4

3

2

1

0

where Ut is the translational velocity of the object, and V is the object volume. Thus the calculation of the fluid force can be simplified to Z dUt ; ð88Þ G ¼  f da1 da2 þ V dt S

Re=20

5

0

10

20

30

40

50

60

70

t Fig. 18. Time history of the drag coefficient for flow past a stationary sphere at different Reynolds numbers.

Table 5 Drag coefficient and wake structure of flow past a sphere

C1 Hs, LTB (xTB, zTB)

Re

10

20

100

200

Current Previous Current Previous Current previous

4.42 4.4 – – – –

3.73 3.8 – – – –

1.15 1.2 130, 0.89 128, 0.90 (0.76, 0.29) (0.76, 0.30)

0.88 0.8 119, 1.44 116, 1.45 (0.90, 0.36) (0.90, 0.36)

scribed by Eq. (84), the amplitudes of the oscillations are significantly reduced, as indicated in Fig. 19. The distributions of the streamwise velocity u along the x-axis are plotted in Fig. 20. Negative values of u at Re = 100 and 200 indicate recirculating flow behind the sphere. No recirculating flow exists at Re = 10 and 20. Streamline plots in Fig. 21 confirm that recirculating flow exists at Re = 100 but not at Re = 10.

2082

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086 100

a

b

80

C

60

C

40

C3

1

3 C

1

2.5

C

2

2

2

C3

1.5

20

1

0 –20

0.5

–40 0 –60 –0.5

–80 –100 0

1

2

3

4

5

–1 0

5

10

15

20

t

t

Fig. 19. Transient behavior of force coefficients for flow past a stationary sphere: (a) Re = 10 with no ramping process, (b) Re = 200 with a ramping process.

subject to a constant. It is axisymmetric around the x-axis with the highest value at the stagnation point at this Reynolds number.

1

0.5

u

7.4. Flow around a flapper

Re=10

0

Re=20 Re=100 Re=200 –0.5 –2

–1

–0

1

2

3

4

5

6

x

Fig. 20. The distribution of the streamwise velocity u along the x-axis.

The separation angle Hs, the length of a separation bubble LTB, and the coordinates of a separation bubble center in the x–z plane (xTB, zTB) are illustrated in Fig. 21b. Their values from the current simulation are compared with previous computational results [26,11] in Table 5, and agree with the previous results very well. In Fig. 22, the surface pressure ps on the sphere is shown for Re = 200. The surface pressure is calculated by ps = Fn

In this example, flow around a hovering flapper is simulated. The spring model given by Eq. (77) with Ks = 100 is used to calculate the density of singular force. The flapper is an ellipsoid with a = 0.4, b = 0.5, and c = 0.2, and its two poles are located in the middle of flat surfaces. The flapper is contained in a rigid box of dimensions ½2; 4  ½1; 1  ½2; 4. The simulation resolution is Nx · Ny · Nz = 129 · 65 · 129 and M1 · M2 = 32 · 64. The time step is fixed in this simulation, and it is Tf Dt ¼ 1:96  103  4000 , where Tf is the flapping period. The motion of the flapper is formulated as p ; ð90Þ xc ¼ 1:25ðcosð0:8tÞ þ 1Þ cos 3 y c ¼ 0; ð91Þ p zc ¼ 1:25ðcosð0:8tÞ þ 1Þ sin ; ð92Þ 3   3p p t þ sinð0:8tÞ 1  exp  ; ð93Þ h¼ 4 4 tc

Fig. 21. Streamlines on the top of contours of the y-component of vorticity at the x–z plane: (a) Re = 10, (b) Re = 100.

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

2083

Fig. 22. Surface pressure distribution at Re = 200 in (a) Cartesian space, (b) parameter space.

where ðxc ; y c ; zc Þ are the coordinates of the flapper center, h is the rotation angle of the flapper with respect to the x-axis, and tc = 1 is a characteristic time of a ramping process to avoid the impulsive start of the flapper rotation. The rotation is around the flapper center, and the rotation direction is given by R = (0, 1, 0). The Reynolds number 2p of the flow is Re = 196. The flapping period is T f ¼ 0:8 .A

2D flapper with the similar kinematics has been simulated in [29,33]. Fig. 23 shows four snapshots of vortex structures during one flapping period. As suggested by Chong et al. [4], vortex structures can be identified by the isosurface of the unit value of Q, where Q is the second invariant of the rate-ofdeformation tensor $v and can be calculated by Q ¼

Fig. 23. Snapshots of vortex structures generated by a flapper.

2084

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086 6

C

Table 6 Relative computational cost for different number of spheres

C2

Number of objects

0

1

2

3

4

C

Relative computational cost

1

1.51

1.94

2.46

2.94

1

5 4

3

3 2

The computational time corresponding to the case with no sphere is 0.202 h.

1 0 –1 –2 –3 –4 0

5

10

15

20

25

30

35

40

t Fig. 24. Time history of force coefficients for flow around a flapper.

 12 rv : rv for incompressible flow. Fig. 23 indicates that a vortex ring is shed downward in the upstroke of the flapper. Fig. 24 shows the time history of the force coefficients C1, C2, and C3. The force coefficients C1, C2, and C3 are defined as ðC 1 ; C 2 ; C 3 Þ ¼ C ¼ 2G=A, where A = pab is the projected area of flat surfaces. In the current simulation setup, C1 is the drag coefficient, C3 is the lift coefficient, and C2 is the side force coefficient. The time average values of C1, C2, and C3 in the first 5 periods are 0.23, 0.0069, and 0.53, respectively. 7.5. Flow induced by multiple spinning spheres In this last example, flow induced by different number of spinning spheres is simulated to examine the efficiency of the current method in handling multiple solids. Each sphere spins around an axis through its center, and the spinning direction is given by R = (0, 1, 0). Again, the model of singular force is given by Eq. (77) with Ks = 1000. Each sphere is represented by M1 · M2 =

32 · 64 Lagrangian markers. The simulation domain is shown in Fig. 25, and it is discretized with Nx · Ny · Nz = 65 · 33 · 65. Also shown in Fig. 25 is contours of u in the x–z plane for the case with four spheres. Table 6 lists the relative computational cost spent on 5000 time steps using the same computer for different number of spheres, and indicates a linear relation with the slope equal to about 0.5. In this example, the ratio of Nx · Ny · Nz to M1 · M2 is about 65, which is relatively small. In the previous example of flow around a flapper, the ratio is about 4 times larger. If this ratio is larger, the slope is expected to become smaller, and the method is relatively more efficient in handling multiple moving solids. 8. Conclusions This paper presents the detailed numerical implementation of a 3D immersed interface method with the jump conditions derived in [32]. The method is tested in simulating (a) flow inside a rotating object, (b) flow induced by a relaxing balloon, (c) flow past a stationary sphere, (d) flow around a flapper, and (d) flow induced by multiple spinning spheres. The test results suggest that (1) the method has near second-order accuracy in the infinity norm for velocity, and the accuracy for pressure is between first and second order; (2) the method conserves the volume enclosed by a no-penetration boundary very well; and (3) the method can efficiently handle multiple moving solids with ease. Future work on the method includes the analysis and improvement of its numerical stability, the use of different interface tracking schemes, and the application of the method to relatively high Reynolds numbers through implicit treatment of solid motion and adaptive mesh refinement. Acknowledgements S. Xu thanks Professor Charles Peskin’s help with his career development. The authors want to thank Professor Randall LeVeque and Professor Zhilin Li for helpful discussions. This work is supported by the AFOSR. Appendix. Analytical solutions to the linear systems in Section 3

Fig. 25. Computational domain for flow involving multiple spinning spheres.

The coefficient matrices C1 and C2 in the small linear systems given by Eqs. (28), (30), (33), and (37) are formed from three independent vectors, s, b, and n. They are

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

non-singular (except at the two poles of a spherical interface). The inverse of C1 is 0 C 1 1 ¼

t1

1B @ t2 J2 t3

b1

n1

1

b2

C n2 A ;

b3

n3

ð94Þ

where ðt1 ; t2 ; t3 Þ :¼ t ¼ b  n, and (b1, b2, b3) :¼ b = n · s. The linear systems with the form of C2q = r can be expanded, split, and rewritten as below 0

s1

B B0 B 0 1B C1 0 0 B 0 B B CB @ 0 C 1 0 A B b1 B 0 0 1 B B0 B B0 @ 1 0 1 0 s1 r1 B C B B r3 C B0 B C B B C B B r4 C B0 B C B B C B ¼ B r 3 C ) B b1 B C B Br C B0 B 2C B B C B Br C B0 @ 5A @ r6 1

0

C 1 1 B ¼B 0 @ 0

0 C 1 1 0

s2

s3

0

0

s1

0

s2

s3

0

s1

0

s2

b2

b3

0

0

b1

0

b2

b3

0

b1

0

b2

0

0

1

0

s2

s3

0

0

s1

0

s2

s3

0

s1

0

s2

b2

b3

0

0

b1

0

b2

b3

0

b1

0

b2

0

1

0 1 C q1 0C C CB q2 C CB B s3 CB C CB q C 3C C C 0 CB B C B q4 C C C 0C CB B C @ q5 C A C b3 A q6 1 1 0 0 1 C q1 0C C CB q2 C CB B s3 CB C CB q C 3C C C 0 CB B C B q4 C C C C 0 CB C CB q A @ 5 b3 C A q6 1

ð95Þ T

where q ¼ ðq1 ; q2 ; q3 ; q4 ; q5 ; q6 Þ and r ¼ ðr1 ; r2 ; r3 ; r4 ; r5 ; r6 Þ . As s is non-zero, at least one component of s is nonzero. If s1 5 0, then 0

1 B B s1 B B B b1 B B B0 B B B0 B B0 @ 0

0

0

1

0

s2

s3

0

0

b2

b3

0

0

s1

0

s2

s3

b1

0

b2

b3

0

s1

0

s2

0

b1

0

b2

1

1

0

1

d7 0 1 C q1 B C 0 CB C B d 1 C CB q C B C 2 C C C B 0 CB B d4 C B C B q3 C B C C C C C¼B 0 CB d 2 C: C B CB C B C B q4 C C C B 0 CB d B 5 C B C CB C q A @ 5 Bd C s3 C A @ 3A q6 b3 d6

then

If s3 5 0, 0 1 0 B B s3 s2 B B b3 b2 B B B 0 s3 B B 0 b3 B B @0 0

then

0

0 0 1 0 0 1 r1 B C B r3 C C 1B B C 0 B r4 C C CB B C T :¼ ðd 1 ; d 2 ; d 3 ; d 4 ; d 5 ; d 6 ; d 7 Þ ; r B 0C 3 AB C C B C 1 B r2 C B C Br C @ 5A r6 T

If s2 5 0, 0 1 0 Bs s B 2 3 B B b2 b3 B B B 0 s2 B B 0 b2 B B @0 0 0 0

Eqs. 0 1 B Ba B Bx B B B0 B B0 B B @0 0

0

0 s1

1 0

0 0

b1 0

0 s3

0 s1

0

b3

b1

s2 b2

0 0

s3 b3

0

1

0

s1 b1

0 0

0 0

0 0

s2 b2

s1 b1

s3

0

s2

b3

0

b2

1 0 1 d7 1 0 1 q4 C B 0 CB C B d 2 C C C B q5 C B C C B 0C d C 5C CB q2 C B CB B C ¼ B d3 C 0 CB C: C CB C q B 6C B C B 0 CB C B d 6 C C C @ q3 A B C s1 A @ d1 A q1 b1 d4 0 1 d7 1 q 6 C B C 0 CB C B d 3 C C B q5 C B C 0C C B d6 C CB C q3 C B CB C C¼B : 0 CB d B 2 C B C CB q B 4C B C C 0C d C B C@ q A B 5 C C B C 2 s1 A @ d1 A q1 b1 d4 1

1

2085

ð97Þ

0

ð98Þ

(96)–(98) can be denoted with the following form: 1 0 1 0 0 1 0 1 0 1 R1 C v1 B C b c 0 0 0 C B C B R2 C C B v2 C B C y z 0 0 0C C B R3 C CB C v3 C B CB B C B a 0 b c 0 CB C ð99Þ ¼ R C; B 4 C B v4 C C C B C C B x 0 y z 0 C B C B R5 C C @ v5 A B C 0 a 0 b cA @ R6 A v6 0 x 0 y z R7

where a 5 0. With Gaussian elimination, the coefficient matrix in Eq. (99) can be transformed to 1 0 1 0 0 1 0 1 C B b c 0 C B0 a 0 C B B0 0 a 0 b c C C B C B ð100Þ B 0 0 0 s3 2bc s2 C; C B C B0 0 0 0 e f C B C B 0 g h A @0 0 0 0

0

0

0

m3

m2

where s2 ¼ a2 þ c2 > 0; s3 ¼ a2 þ b2 > 0; m2 ¼ cx  az; m3 ¼ ay  bx; e ¼ 2bcðax þ byÞ  s3 ðcy þ bzÞ;

ð96Þ

f ¼ s2 ðax þ byÞ  s3 ðax þ czÞ; g ¼ s3 m2  2bcm3 ; h ¼ s2 m3 : Since a and n are non-zero except at the poles of an ellipsoidal interface, one of m2 and m3 must be non-zero [32].

2086

S. Xu, Z.J. Wang / Comput. Methods Appl. Mech. Engrg. 197 (2008) 2068–2086

The right lower corner matrix in (100) can be transformed to 0 1 0 m3 f þ m2 e B C ð101Þ @ 0 m3 h þ m2 g A; m3

m2

if m3 5 0, or 0 m3 f þ m2 e B @ m3 h þ m2 g m3

0

1

C 0 A; m2

ð102Þ

if m2 5 0. It can be shown that m3h + m2g < 0 [32]. So far, Eq. (95) has been solved analytically. References [1] Attila J. Bergou, Sheng Xu, Z. Jane Wang, Passive wing rotation in dragonfly flight, J. Fluid Mech., accepted for publication. [2] Petter Andreas Berthelsen, A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions, J. Comput. Phys. 197 (2004) 364–386. [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] M.S. Chong, A.E. Perry, B.J. Cantwell, A general classification of three-dimensional flow fields, Phys. Fluids A 2 (5) (1990) 765–777. [5] R. Cortez, M. Minion, The blob projection method for immersed boundary problems, J. Comput. Phys. 161 (2000) 428–453. [6] Shaozhong Deng, Kazufumi Ito, Zhilin Li, Three-dimensional elliptic solvers for interface problems and applications, J. Comput. Phys. 184 (2003) 215–243. [7] W. E, Jian-Guo Liu, Vorticity boundary condition and related issues for finite difference schemes, J. Comput. Phys. 124 (1996) 368–382. [8] Aaron L. Fogelson, James P. Keener, Immersed interface method for Neumann and related problems in two and three dimensions, SIAM J. Sci. Comput. 22 (5) (2000) 1630–1654. [9] Hans Johnston, Jian-Guo Liu, Finite difference schemes for incompressible flow based on local pressure boundary conditions, J. Comput. Phys. 180 (2002) 120–154. [10] Hans Johnston, Jian-Guo Liu, Accurate stable and efficient Navier– Stokes solvers based on explicit treatment of the pressure term, J. Comput. Phys. 199 (2004) 221–259. [11] T.A. Johnson, V.C. Patel, Flow past sphere up to a Reynolds number of 300, J. Fluid Mech. 378 (1999) 19–70. [12] Ming-Chih Lai, Charles S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comput. Phys. 160 (2000) 705–719. [13] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys. 220 (2006) 109–138. [14] Long Lee, Randall J. LeVeque, An Immersed Interface Method for Incompressible Navier–Stokes Equations, SIAM J. Sci. Comput. 25 (3) (2003) 832–856.

[15] Randall J. LeVeque, Zhilin Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (4) (1994) 1019–1044. [16] Randall J. LeVeque, Zhilin Li, Immersed interface methods for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput. 18 (3) (1997) 709–735. [17] Zhilin Li, Ming-Chih Lai, The immersed interface method for the Navier–Stokes equations with singular forces, J. Comput. Phys. 171 (2001) 822–842. [18] Zhilin Li, Kazufumi Ito, The Immersed Interface Method – Numerical Solutions of PDEs Involving Interfaces and Irregular Domains, SIAM Frontiers in Applied Mathematics, vol. 33, SIAM Publisher, 2006, ISBN: 0-89971-609-8. [19] Mark N. Linnick, Hermann F. Fasel, A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains, J. Comput. Phys. 204 (2005) 157–192. [20] Charles S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252–271. [21] Charles S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [22] Charles S. Peskin, Beth Feller Printz, Improved volume conservation in the computation of flows with immersed elastic boundaries, J. Comput. Phys. 105 (1993) 33–46. [23] Charles S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [24] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Numerical Recipes in Fortran 77, second ed., The Art of Scientific Computing, 1999, pp. 848–852. [25] A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (1999) 509– 534. [26] G.J. Sheard, K. Hourigan, M.C. Thompson, Computations of the drag coefficients for low-Reynolds-number flow past rings, J. Fluid Mech. 526 (2005) 257–275. [27] John M. Stockie, Brian R. Wetton, Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes, J. Comput. Phys. 154 (1999) 41–64. [28] C. Tu, C.S. Peskin, Stability and instability in the computation of flows with moving immersed boundaries, SIAM J. Statist. Comput. 13 (1992) 70–83. [29] Z. Jane Wang, Two dimensional mechanism for insect hovering, Phys. Rev. Lett. 85 (10) (2000). [30] Andreas Wiegmann, Kenneth P. Bube, The immersed interface method for nonlinear differential equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 35 (1) (1998) 177–200. [31] Andreas Wiegmann, Kenneth P. Bube, The explicit-jump immersed interface method: finite difference methods for PDEs with piecewise smooth solutions, SIAM J. Numer. Anal. 37 (3) (2000) 827– 862. [32] Sheng Xu, Z. Jane Wang, Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation, SIAM J. Sci. Comput. 27 (6) (2006) 948–1980. [33] Sheng Xu, Z. Jane Wang, An immersed interface method for simulating the interaction of a fluid with moving boundaries, J. Comput. Phys. 216 (2) (2006) 454–493. [34] Luoding Zhu, Charles S. Peskin, Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method, J. Comput. Phys. 179 (2002) 452–468.