A discrete-forcing immersed boundary method for the fluid–structure interaction of an elastic slender body

Journal of Computational Physics 280 (2015) 529–546

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A discrete-forcing immersed boundary method for the fluid–structure interaction of an elastic slender body Injae Lee, Haecheon Choi ∗ Department of Mechanical & Aerospace Engineering, Seoul National University, Seoul, 151-744, Republic of Korea

a r t i c l e

i n f o

Article history: Received 11 January 2014 Received in revised form 12 July 2014 Accepted 24 September 2014 Available online 2 October 2014 Keywords: Immersed boundary method Fluid–structure interaction Elastic slender body Discrete forcing CFL number

a b s t r a c t We present an immersed boundary (IB) method for the simulation of flow around an elastic slender body. The present method is based on the discrete-forcing IB method for a stationary, rigid body proposed by Kim, Kim and Choi (2001) [25]. The discreteforcing approach is used to relieve the limitation on the computational time step size. The incompressible Navier–Stokes equations are implicitly coupled with the dynamic equation for an elastic slender body motion. The first is solved in the Eulerian coordinate and the latter is described in the Lagrangian coordinate. The elastic slender body is modeled as a thin and flexible solid and is segmented by finite number of thin blocks. Each block is moved by external and internal forces such as the hydrodynamic, elastic and buoyancy forces, where the hydrodynamic force is obtained directly from the discrete forcing used in the IB method. All the spatial derivative terms are discretized with the second-order central difference scheme. The present method is applied to three different fluid–structure interaction problems: flows around a flexible filament, a flapping flag in a free stream, and a flexible flapping wing in normal hovering, respectively. Computations are performed at maximum CFL numbers of 0.75–1. The results obtained agree very well with those from previous studies. © 2014 Elsevier Inc. All rights reserved.

1. Introduction There are many elastic slender bodies interacting with ambient fluid such as flags fluttering in the wind and flapping insect wings [1–3]. Motions of these elastic slender bodies are of interest because they may provide us new ideas for bio-mimetic engineering applications such as an energy harvesting eel or wings of micro air vehicles [4,5]. However, understanding the dynamics of fluid–structure interaction is quite complicated, because an elastic slender body, unlike rigid one, is significantly deformed by the fluid. Thus, for the movement and deformation of an elastic slender body, an efficient and accurate numerical method is required. There have been quite a few numerical studies about elastic slender bodies using a conventional body-fitted or unstructured mesh to handle fluid–structure interaction (see, for example, [6–9]). However, they suffer from complex mesh generation and frequent remeshing that may lead to losing accuracy during the procedure of projecting a solution from old mesh to new mesh. On the other hand, the immersed boundary (called IB hereafter) method initially developed by Peskin [10] uses Cartesian or cylindrical mesh and does not require mesh regeneration even for a deforming body. Therefore, the IB method is appropriate for the fluid–structure interaction of an elastic slender body in terms of its efficiency and convenience.

*

Corresponding author. Also at: Institute of Advanced Machines and Design, Seoul National University, Republic of Korea. E-mail address: [email protected] (H. Choi).

http://dx.doi.org/10.1016/j.jcp.2014.09.028 0021-9991/© 2014 Elsevier Inc. All rights reserved.

530

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Over the last few decades, many researchers have conducted simulations of flows around elastic slender bodies by developing various versions of IB methods. Among them, Peskin [10] firstly simulated the blood flow in human heart assuming a low Reynolds number and two-dimensional flow. For the fluid–structure interaction inside human heart, the IB was made of moving Lagrangian points linked by springs, and the discrete delta function was used to interpolate the fluid velocity for Lagrangian points and to spread out the momentum forcing on Lagrangian points to fluid grid points (so is called continuous-forcing IB method). It is known that the discrete delta function has a role of attenuating the spurious force oscillations [11,12]. This original IB method has been widely used to deal with various problems of elastic slender bodies such as animal locomotion, valveless pumping, and parachute [13–15]. Furthermore, there have been improved versions of this IB method for elastic slender bodies [16–19] by considering the mass of a solid body, but they bring about a smearing interface problem from mass handling owing to the use of additional discrete delta function or virtual boundary, which decreases the order of accuracy near the body. These IB methods [13–19] based on the original one [10] have a clear advantage of their straightforward implementation into existing solvers [20], and satisfy the strong coupling between fluid and structure inherently because they use the same numerical method for both fluid and structure equations (monolithic algorithm). However, they have a fundamental difficulty in handling the structure with mass. To handle the mass of structure more strictly, researchers employed fluid and structure solvers independently (partitioned algorithm). Huang et al. [21] adopted a feedback forcing approach [22] for the interaction between the fluid and flexible filament, where the inextensibility condition of the filament was satisfied by treating the tension force as the constraint of inextensibility. In this study, the discrete delta function was also used for interaction as done in the original IB method [10]. The two additional constants, α and β , used in the feedback forcing approach [22] had to be very large to obtain results especially for unsteady flow problems. These large values of constants provided a severe limitation on the size of computational time step. In another method [23], the momentum forcing was formulated from the inertial term of the structure equation to reduce the number of constants, but there still remained a free constant to be tuned. It is remarkable to note that these IB methods explicitly/semi-implicitly integrate the structure equation in time, which can cause the numerical instability. Nevertheless, this explicit/semi-implicit treatment of the structure equation may be acceptable because the time-step constraint from the IB method is less than or equal to that from the structure equation. In case stiff material is considered, an implicit treatment of the structure equation should be adopted. Depending on how to implement the momentum forcing, the IB method is classified as two groups, continuous forcing and discrete forcing [24]. An important advantage of the discrete-forcing IB method over the continuous-forcing IB method is that the first has moderate limitation on the size of the computational time step [24–26]. Although previous studies for the fluid–structure interaction of an elastic body are mainly based on the continuous-forcing IB method, there have been some studies using the discrete-forcing approach [27–33]. The vocal fold vibration during phonation was studied as an example of biological flows [27], and the vortex-induced vibration of a flexible splitter behind a circular/rectangular cylinder was studied by varying the material properties and Reynolds numbers [28,29]. Luo and co-workers [30–33] studied flapping flights with a larger size of the computational time step than those used in the continuous-forcing IB methods [15–18,23,34]. However, these discrete-forcing IB methods need additional interpolation schemes on the IB to obtain the hydrodynamic force or diminish the sharpness of the interface owing to the flow reconstruction remedy to reduce spurious force oscillations [12]. As mentioned above, there are many application areas associated with elastic slender bodies interacting with ambient fluid. For this type of bodies, there is a need to develop a numerical method that is more efficient than those developed for general bodies. Therefore, in the present study, we develop an IB method using the discrete-forcing approach [25] for the fluid–structure interaction of an elastic slender body and to relieve the limitation on the size of the computational time step. We obtain the hydrodynamic force directly from the Navier–Stokes equations without any interpolation on the IB and maintain the sharpness of the IB by avoiding any treatment on the fluid node. We test our method to three different fluid–structure interaction problems: flows around a flexible filament, a flapping flag in a free stream, and a flexible flapping wing in normal hovering. We show that the present numerical method accurately predicts these flows at moderate CFL numbers of 0.75–1. 2. Numerical method The present numerical method for the fluid–structure interaction of an elastic slender body is constructed by coupling fluid and structure systems. The incompressible Navier–Stokes equations are solved in the Eulerian coordinate, whereas the motion of elastic slender body is described in the Lagrangian coordinate. Construction of the elastic slender body surface and coupling of the fluid and elastic slender body are given in the framework of the discrete-forcing IB method. 2.1. Navier–Stokes equations for fluid motion in the Eulerian coordinate The present IB method is based on the discrete-forcing IB method developed for a stationary, rigid body [25]. In this method, the non-dimensional governing equations for unsteady incompressible viscous flow are given as

∂u 1 + ∇ · (uu ) = −∇ p + ∇ 2 u + f , ∂t Re ∇ · u − q = 0,

(1) (2)

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 1. Locations of the discrete momentum forcing and mass source/sink [25]. The solid and dashed lines denote the grid lines and IB, respectively. discrete forcing; , mass source/sink; , fluid velocity.

531

,

where u is the velocity, p is the pressure, f is the momentum forcing to satisfy the no-slip boundary condition on the IB, q is the mass source/sink to satisfy the mass conservation for the cell containing the IB, and Re = U L /ν is the Reynolds number (U and L are the characteristic velocity and length scales, and ν is the kinematic viscosity). In the staggered mesh system, f is applied directly at the grid cell faces, and q is applied at the cell center (Fig. 1). A fractional step method is applied to solve Eqs. (1) and (2), and a semi-implicit time advancement scheme is used (a third-order Runge–Kutta method (RK3) for the convection term and the Crank–Nicolson method for the diffusion term):

uˆ − uk−1 k

t ∇ 2 φk =

 k       = αk L uˆ + αk L uk−1 − 2αk ∇ pk−1 − γk N uk−1 − ρk N uk−2 + 2αk f k , 

1 2αk t

 k ∇ · uˆ − qk ,

(4)

k

uk = uˆ − 2αk t ∇φ k , p k = p k −1 + φ k −

αk t Re

(3)

(5)

∇ 2 φk ,

(6)

where L (u ) = ∇ 2 u / Re, N (u ) = ∇ · (uu ), uˆ is the intermediate velocity, φ is the pseudo-pressure, t is the computational time step, k is the RK3 substep index (k = 1, 2, 3), and αk , γk , and ρk are the coefficients of RK3 substeps (α1 = 4/15, γ1 = 8/15, ρ1 = 0; α2 = 1/15, γ2 = 5/12, ρ2 = −17/60; α3 = 1/6, γ3 = 3/4, ρ3 = −5/12). k

The momentum forcing f k in Eq. (3) is determined in advance such that uˆ satisfies the no-slip condition on the IB k k instead of uk , which does not affect the overall second-order temporal accuracy because u k = uˆ − 2αk t ∇φ k = uˆ + ϑ(t 2 ). k To determine f near the IB, Eq. (1) is provisionally discretized explicitly in time (RK3 for the convection term and forward Euler method for the diffusion term) only for the fluid cells nearby the IB ( f k = 0), as done in [25]:

u˜ − uk−1 k

t

      = 2αk L uk−1 − 2αk ∇ pk−1 − γk N uk−1 − ρk N uk−2 ,

(7)

k

where u˜ is the provisional velocity obtained for interpolation purpose. Now, U k is the velocity that we want to obtain at a forcing point just inside the IB by applying the momentum forcing. Its value is determined using the linear, bilinear k or trilinear interpolation from the velocities u˜ [25]. Then, f k just inside the IB is obtained as follows (with the same integration method as in Eq. (7)):

2αk f k =

U k − u k −1

t

      − 2αk L uk−1 + 2αk ∇ pk−1 + γk N uk−1 + ρk N uk−2 .

(8)

On the other hand, we obtain f k at a grid location further inside the IB from Eq. (8) with U k = ukb , where ukb is the solid-body velocity at that location but is constant within a segmented block (see Section 2.2). The procedure of obtaining qk in Eq. (4) is described in brief here and more details are found in [25]. In Fig. 2, uk1 and v k1

are the velocity components inside the body, uk2 and v k2 are those outside the body, and uks and v ks are those on the IB. For the triangular cell containing only fluid ( P B P C P D ) and the rectangular cell containing both the solid and fluid ( P A P B P C P D ), the continuity reads

532

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 2. Mass conservation for a cell containing the IB.

uk2  y + v k2 x = uks  y + v ks x, uk2  y

+

By replacing

uk1

qk =

v k2 x and

v k1



1

x y

=

uk1  y with

+

uˆ k1

(9)

v k1 x + qk x y .

and



vˆ k1

(10) k

without losing the overall temporal accuracy, the mass source/sink q is obtained as







uks − uˆ k1  y + v ks − vˆ k1 x .

(11)

Note that Eq. (11) results from a simple geometric configuration between the solid and fluid as depicted in Fig. 2. For general configurations, one may want to use a more complicated (but accurate) formula for qk and should consult Refs. [35,36]. In moving body problems, the solid body passes through fixed mesh, and consequently the locations for the momentum forcing f k and the mass source/sink qk change in time. Therefore, at every time step, the procedure to identify the fluid and solid nodes is required. In three-dimensional problems, the computational cost would be very expensive if this identification procedure is carried out for whole nodes. To avoid this unnecessary searching process, we search nodes only near the solid nodes in the previous time step. 2.2. Dynamic equation for the motion of an elastic slender body in the Lagrangian coordinate In our approach, an elastic slender body is segmented by finite number of blocks (Fig. 3). The total number of blocks is N b . Each block has a central point X in the Lagrangian coordinate (s1 , s2 ), a finite thickness of h, and Lagrangian lengths of s1 and s2 . The center of each block is moved by the external and internal forces such as the elastic ( F E ), buoyancy ( F G ), and hydrodynamic ( F H ) forces, and the dynamic equation for each block is described by the Newton second law:

ρs V

∂2 X = F E + F G + F H, ∂t2

(12)

where ρs is the density of elastic slender body, V = hs1 s2 . The elastic force is composed of the tension ( F T ), shearing ( F S ), bending ( F B ), and twisting ( F W ) forces (Fig. 3(b)). The principle of virtual work is used to obtain the elastic force as done in [34]: the elastic energy is expressed as

E(X) =

¨  2 



T 0 cab T ab − T ab

2

   B 0 2 + cab B ab − B ab ds1 ds2 ,

(13)

a,b=1

where T ab (= ∂ X /∂ sa · ∂ X /∂ sb ) is the tension or shearing effect, B ab (= ∂ 2 X /∂ sa ∂ sb · ∂ 2 X /∂ sa ∂ sb )1/2 is the bending or T B twisting effect, and cab and cab are the tension or shearing constant and bending or twisting constant, respectively. The 0 0 initial conditions, T ab and B ab , for thin flat body are given as



0 T ab

=

1, if a = b, 0, if a = b,

and

0 B ab = 0.

(14)

The elastic force F E is derived from the variational derivative of the elastic energy,







 2  2  ∂2 ∂ ∂X ∂X ∂X T 0 B ∂ X − FE = κ · − T ab κ A. ∂ sa ab ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ab ∂ sa ∂ sb a,b=1

(15)

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

533

Fig. 3. Elastic slender body: (a) overall shape; (b) forces on each block; (c) movement of the center of each block; (d) surface reconstruction. In (d), the distance between 2’s is h, and N is the pseudo-normal vector. T T B B Here, κab (= 4cab ) is the tension (a = b) or shearing (a = b) coefficient, κab (= 2cab ) is the bending (a = b) or twisting (a = b) coefficient, and A = s1 s2 . Then, the dynamic equation for each block of an elastic slender body is expressed as [7,34]

ρs V







 2  2  ∂2 ∂2 X ∂ ∂X ∂X ∂X T 0 B ∂ X − = κ · − T κ A + ρs V g + F H , ab ∂ sa ab ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ab ∂ sa ∂ sb ∂t2

(16)

a,b=1

where g is the gravitational acceleration, and F H is the hydrodynamic force acting on the body (see in the below). The non-dimensional equation corresponding to Eq. (16) is

ρ∗

∗

∗

2  2 ∗   F∗ ∂ ∂ X ∂ X∗ ∂X ∂2 ∂2 X∗ T 0 B ∂ X − + ρ ∗ ˆi g Fr + H∗ , = K · − T K ab ab ab ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ 2 ∂ sa ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ∂ sa ∂ sb V ∂t

(17)

a,b=1

T T B B where ρ ∗ = ρs /ρ f (ratio of solid to fluid densities), K ab = κab /ρ f U 2 h, K ab = κab /ρ f U 2 L 2 h, ˆi g = g /| g |, Fr = | g | L /U 2 (Froude ∗ 2 2 ∗ number), F H = F H /ρ f U L (non-dimensional hydrodynamic force), and V = s1 s2 h/ L 3 (non-dimensional volume of each block), and L is the characteristic length defined for each flow problem. The non-dimensional tension or shearing T coefficient K ab is chosen to be sufficiently large to satisfy the inextensibility condition [34]. From now on, we drop the superscript ∗ for convenience. Note that the structure model described by Eq. (17) is valid only for fibrous shells in which the stretching and bending stiffness are independent from each other [21,34]. This model is adequate for highly inextensible structures such that their displacement in the normal direction is much larger than that in the tangential direction. For general structures having complex geometries and large thickness, one may adopt other structure models such as the Saint Venant–Kirchhoff model [33].

534

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 4. Schematic diagram for the pseudo-normal vector at the center of each block: (a) X (!; center of a block) and C ("; geometric centroid of four adjacent centers of blocks); (b) outward normal vectors from the discretized surfaces.

In this study, an implicit second-order finite-difference scheme is used to discretize Eq. (17) in time as similar to [6,7]: for constant t,

ρ

X n+1 − 2X n + X n−1

t 2

=

2 1 

2

a,b=1

+



n +1

n +1

2 n+1  ∂2 ∂ X n +1 ∂ ∂X ∂X T 0 B ∂ X − K ab · − T ab K ab ∂ sa ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ∂ sa ∂ sb

2 1 

2



a,b=1

+ ρ ˆi g Fr +

n

n

2 n  ∂2 ∂ ∂ X ∂ Xn ∂X T 0 B ∂ X − K ab · − T ab K ab ∂ sa ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ∂ sa ∂ sb 1 1 2V



F nH+1 + F nH .

(18)

The center location of each block moves in time according to Eq. (18) as shown in Fig. 3(c). Note that Eq. (18) is obtained from a time integration of Eq. (17) using a constant time step, and thus it should be modified when the structure equation is coupled with the Navier–Stokes equations, because the numerical method used for the latter, Eq. (3), requires three substeps per time step and also allows variable time step sizes during time integration. For the details, see Section 2.4. 2.3. Construction of the surface of an elastic slender body Now, we explain how to construct the surface of each block. Fig. 4 shows the schematic diagram of the geometric components and the pseudo-normal vectors at discretized surfaces. Suppose the discretized surface is composed of centers of blocks ( X ) as shown in Fig. 4(a). Let us consider X I , J at which we want to obtain the pseudo-normal vector. This center of the block is surrounded by adjacent centers, X I −1, J −1 , X I −1, J , · · · , X I +1, J +1 . The geometric centroid C I , J of four adjacent centers is calculated using the arithmetic mean for the triangulation of the discretized surface:

C I, J =

X I −1 , J −1 + X I −1 , J + X I , J −1 + X I , J 4

.

(19)

For convenience, let us denote the center X I , J as O which is located at the origin, and the adjacent centers of blocks and geometric centroids as V i (= X and C ) (Fig. 4(b)). The face-normal vector N i of each triangle O V i V i +1 is given as

Ni =

V i × V i +1 , | V i || V i +1 | sin αi

(20)

where αi is the angle between V i and V i +1 (here, V 9 = V 1 ). Several algorithms to compute the pseudo-normal vector at a center have been suggested [37–39]. The common strategy to get a pseudo-normal vector is the estimation by weighted sum of adjacent face-normal vectors defined on neighboring triangles, but each algorithm is identified by different weighting factor. Among these algorithms, the method proposed by Max [39] is used here; i.e., with the weighting factor ωi (= sin αi /| V i || V i +1 |), the pseudo-normal vector at the center O is obtained as

N(O ) =

8  i =1

ωi N i =

8  V i × V i +1 . | V i |2 | V i + 1 |2

(21)

i =1

The Max algorithm is also suitable when the faces surrounding the center differ greatly in size [40]. Using this pseudonormal vector at the center of each block, we reconstruct the body surface with maintaining the thickness of h as shown in Fig. 3(d).

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

535

2.4. Coupling between the fluid and elastic slender body In our numerical method, the Navier–Stokes equations and the dynamic equation for the elastic slender body are solved separately in two different coordinate systems. Therefore, coupling between these two different systems is needed for the fluid–structure interaction, and this is accomplished by the imposition of the no-slip boundary condition on the IB and by the hydrodynamic force acting on the solid. The no-slip boundary condition on the IB is satisfied by applying the momentum forcing in the Navier–Stokes equations (see Section 2.1). Therefore, the procedure to obtain the hydrodynamic force used in Eq. (18) is explained hereafter. The hydrodynamic force acting on each block of solid body is

ˆ 

− pn +

FH =

  ∇ u + ∇ u T · n dS ,

1  Re

(22)

Sb

where n is the outward-facing normal vector on the wetted surface (S b ) of each block. Connell and Yue [6] directly calculated the right hand side of Eq. (22) for a moving thin body using a body-fitted mesh. However, in IB methods like the present one, additional interpolations are required to obtain the pressure and stresses on the IB, which inevitably introduce errors. Thus, to avoid additional interpolations, we perform the volume integration of the Navier–Stokes equations as done in [11],

ˆ



ˆ ˆ 1 ∂u −∇ p + ∇ 2 u dV + + ∇ · (uu ) dV = f dV , ∂t Re

Vb

Vb

(23)

Vb

where V b is the volume of each block. Using the divergence theorem, the first term in the right hand side of Eq. (23) becomes

ˆ

−∇ p +

1 Re



ˆ    1 −∇ p + ∇ · ∇ u + ∇ u T dV ∇ 2 u dV = Re

Vb

Vb

ˆ  − pn +

=

  ∇ u + ∇ u T · n dS

1  Re

Sb

= F H.

(24)

Then, from Eqs. (23) and (24), the hydrodynamic force becomes

ˆ

FH =

∂u + ∇ · (uu ) − f dV . ∂t

(25)

Vb

As shown in Eq. (25), the hydrodynamic force on each block is directly obtained from the Navier–Stokes equations as sum of unsteady, inertial and momentum forcing terms without any interpolation on the surface. Note that this scheme is only feasible for one-layer elastic slender bodies considered in the present study (Fig. 3). If an elastic solid body is composed of two or more layers or is a thick bluff body, the surface interpolation to obtain the right hand side of Eq. (22) should be inevitable. There are monolithic and partitioned algorithms for the purpose of coupling between the fluid and structure [41]. The former uses the same numerical method for both fluid and structure equations, and thus the coupling between both systems is satisfied inherently. The original IB method by Peskin [10] is a representative for monolithic algorithms. On the other hand, the latter makes use of different numerical methods for the fluid and structure equations and gives more flexibility in selecting suitable methods for each equation, respectively [41]. Although, the partitioned algorithm requires additional iteration strategy for strong coupling to avoid the time lagging between the fluid and structure, it is more proper than the monolithic one especially in handling the elastic body with mass. Therefore, the partitioned algorithm adopting the strong coupling scheme is used in the present study. As mentioned before, the time integration scheme for the structure equation, Eq. (18), should be modified when it is coupled with the Navier–Stokes equations, because of the three substeps per time step and variable time step size used in Eq. (3). Thus, the strong coupling procedure for Eqs. (3)–(8), (11) and (18) is

ρ

( X k,m+1 − X n )to − ( X n − X n−1 )t p t p to (t p + to )/2 k,m+1

k,m+1

2  2 k,m+1   ∂ ∂X ∂X ∂2 ∂ X k,m+1 1 T 0 B ∂ X − = K ab · − T ab K ab 2 ∂ sa ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ∂ sa ∂ sb a,b=1

+

2 1 

2

a,b=1



n

n

2 n   ∂2 1 1  k−1,m ∂ ∂ X ∂ Xn ∂X T 0 B ∂ X − + ρ ˆi g Fr + K ab · − T ab K ab FH + F nH , (26) ∂ sa ∂ sa ∂ sb ∂ sb ∂ sa ∂ sb ∂ sa ∂ sb 2V

536

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

2αk f k,m+1 =

U k,m+1 − uk−1,m+1

      − 2αk L uk−1,m+1 + 2αk ∇ pk−1,m+1 + γk N uk−1,m+1 + ρk N uk−2,m+1 , (27)

t  k,m+1        − uk−1,m+1 + αk L uk−1,m+1 − 2αk ∇ pk−1,m+1 − γk N uk−1,m+1 − ρk N uk−2,m+1 = αk L uˆ t + 2αk f k,m+1 ,     1  k,m+1 us − uˆ k1,m+1  y + v ks ,m+1 − vˆ k1,m+1 x , qk,m+1 = x y  1  k,m+1 ∇ · uˆ ∇ 2 φ k,m+1 = − qk,m+1 , 2αk t uˆ

k,m+1

uk,m+1 = uˆ

k,m+1

− 2αk t ∇φ k,m+1 , αk t 2 k,m+1 ∇ φ , pk,m+1 = pk−1,m+1 + φ k,m+1 − Re

ˆ uk,m+1 − uk−1,m+1 γk  k−1,m+1  ρk  k−2,m+1  k,m+1 k,m+1 + − f FH = + N u N u dV , 2αk t 2αk 2αk

(28) (29) (30) (31) (32) (33)

Vb

k

k ,0

where m is the overall iteration index, t p = i =1 2αi t, to = tn − tn−1 , uk,1 = un , pk,1 = pn , and F H = F nH . For the treatment of the non-linear term in Eq. (26), a linearization technique (let X k,m+1 = X n + δ X and neglect O (δ X 2 )) is used [7]. The spatial derivative terms in this equation are discretized with the second-order central difference scheme as fully described in [23], and the resulting equation is solved using the Gauss–Seidel method [23]. About 80 iterations are needed to obtain a converged solution (| X k,m+1,l+1 − X k,m+1,l | < 10−12 , where l is the Gauss–Seidel iteration index). Note that a structured mesh used here can easily model any type of thin solid geometries. The second-order central difference scheme is also used for all the spatial derivative terms in Eqs. (27)–(33). The whole iterative procedure is repeated until the convergence criterion | X k=3,m+1 − X k=3,m | < ε is satisfied, and then X n+1 = X k=3,m+1 , un+1 = uk=3,m+1 , and pn+1 = pk=3,m+1 . Typically, 3–5 overall iterations are required to reach the convergence criterion of ε = 10−8 . In the present study, the momentum forcing is implemented inside the body to satisfy the no-slip boundary condition (Fig. 1). Since the elastic slender body has thin but finite thickness, the Eulerian grid sizes near the path of elastic slender body should be smaller than its thickness h to enforce the momentum forcing inside the body, and they are equal to or smaller than the Lagrangian lengths s1 and s2 to satisfy the no-slip boundary condition on the IB and estimate the hydrodynamic force. We conducted simulations with the Eulerian grid sizes equal to or four times smaller than the Lagrangian lengths, and found that there were nearly no difference among the solutions obtained. 3. Numerical examples The present numerical method is applied to three different fluid–structure interaction problems of elastic slender bodies: flows around a flexible filament, a flapping flag in a free stream, and a flexible flapping wing in normal hovering, respectively. All the computations are carried at maximum CFL = 0.75∼1. For the first flow problem, we discuss the convergence and stability of the present numerical method. The results obtained agree very well with those from previous studies, indicating the accuracy and efficiency of the present IB method. 3.1. Two-dimensional flexible filament in a free stream We simulate a flexible-filament motion in a free stream that is initially inclined with respect to the free stream direction. The material parameters of flexible filament are nearly the same as those used in Huang et al. [21]. In [21], the inextensibility condition is satisfied by treating the tension force term in Eq. (17) as an implicit constraint, whereas a large tension T coefficient K 11 is given in the present study to satisfy the inextensibility condition. As shown in Fig. 5, a slender filament of length L is pinned (but not clamped) at the leading edge and freely moves elsewhere. The directions of the free stream velocity and gravity vectors are the same as the x-direction, and the initial inclined angle between the filament and x-axis is θ = 0.1π . The computational domain size is [−2L , 6L ] × [−4L , 4L ] and the number of grid points is 513 × 481 in the streamwise (x) and transverse ( y) directions. The smallest grid spacing is x =  y = 0.005L, and uniform grids with this smallest grid spacing are distributed on and near the path of the filament. The Dirichlet boundary condition (u = U , v = 0) is applied at the inflow and side boundaries, and a convective boundary condition (∂ u i /∂ t + c ∂ u i /∂ x = 0) is given at the outflow boundary, where c is the velocity averaged over the outflow boundary. The simply supported boundary condition is used at the pinned boundary ( X |s1 =0 = X 0 , ∂ 2 X /∂ s21 |s1 =0 = 0), and the free-end boundary condition is applied at the freely-moving boundary (∂ X /∂ s1 · ∂ X /∂ s1 |s1 = L = 1, ∂ 2 X /∂ s21 |s1 = L = 0, ∂ 3 X /∂ s31 |s1 =L = 0). The values of other dimensionless material and numerical parameters are given as ρ = 100 or 150, B T K 11 = 0.1 or 0.15, K 11 = 2500, Fr = 0.5, N b = 50, and h/ L = 0.01. Numerical simulations are conducted for two different Reynolds numbers of 200 and 300 based on the free stream velocity U and filament length L at the maximum CFL number

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

537

Fig. 5. Schematic diagram of the computational domain and boundary conditions for the simulation of flow around a flexible filament in a free stream.

Fig. 6. Time traces of the trailing-edge transverse location of a flexible filament of 0.02; , 0.03; , 0.04.

B ρ = 150 and K 11 = 0.15 at Re = 200:

, h/ L = 0.01;

,

of 1.0. Note that the use of this relatively high CFL number is attributed to the use of the discrete-forcing approach. We first conduct a convergence test by increasing the number of grid points (513 × 481, 609 × 581, 705 × 681, respectively) for the B case of ρ = 150 and K 11 = 0.15 (Re = 200) and obtain nearly the same solutions. Next, we conduct a convergence test for the number of blocks N b as similar to [21], and obtain nearly same results when N b ≥ 40. The dependence of numerical solution by the filament thickness h is also examined. As shown in Fig. 6, nearly no difference is observed for h/ L ≤ 0.03. Therefore, the filament thickness of h/ L = 0.01 is proper for the present simulation. Fig. 7 shows the time traces of the trailing-edge transverse location of a flexible filament, together with those in [21]. The flexible filament reaches flapping states for both cases. The present trailing-edge locations agree well with those of [21], although there are slight phase differences between two results. Note that the present grid spacing in the transverse direction with h/ L = 0.01 is twice finer than that of [21] because a few grids should be located inside the structure in this direction for the present IB method. A half grid resolution in the transverse direction with h/ L = 0.03, which is comparable to that of [21], results in nearly the same solution as that from the fine grid resolution (not shown here). Fig. 8 shows the trajectories of the trailing-edge location during t = 30∼40. The ‘figure of eight’ trajectory due to the flapping motion is manifest in this figure. Fig. 9 shows the contours of the instantaneous vorticity around the flexible filament for the case B of ρ = 150 and K 11 = 0.15 (Re = 200) at four different time instants during one cycle. As expected, the vortices are shed alternatively behind the flexible filament. Fig. 10 shows the time traces of the drag and lift coefficients, where the drag and lift coefficients are defined as C d = 2D/ρ U 2 L and C l = 2L/ρ U 2 L, respectively, D is the drag force, and L is the lift force. In this figure, we also show the results obtained from two different ways of computing the drag and lift forces, i.e., from Eqs. (22) and (25). As shown, serious spurious oscillations are found in the force coefficients obtained from Eq. (22). Lee et al. [11] identified the sources of these spurious oscillations, and showed that these oscillations are reduced with reducing the grid spacing and increasing

538

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 7. Time traces of the trailing-edge transverse location of a flexible filament: (a) Re = 300. , present study; 2, Huang et al. [21].

Fig. 8. Trajectories of the trailing-edge location at 30 ≤ t ≤ 40:

,

B B ρ = 150 and K 11 = 0.15 at Re = 200; (b) ρ = 100 and K 11 = 0.1 at

B ρ = 150 and K 11 = 0.15 at Re = 200;

,

B ρ = 100 and K 11 = 0.1 at Re = 300.

the computational time step size for moving rigid body problems. Accordingly, we conducted a simulation with half the grid size used, and the spurious force oscillations from Eq. (22) were indeed reduced with smaller grid spacing for the present flexible body problem. To reduce these spurious force oscillations more strictly, a fully implicit time integration scheme [26] may be used because it allows the use of larger computational time step size [11]. Other ways of reducing the spurious force oscillations are also found in [12,20,42,43]. A flexible filament sustains the flapping state or goes to the stable state depending on the values of the density ratio B and bending coefficient, and the boundary between two states is given as ρcr = 1.3 Re−0.5 /(h/ L ) + 4K 11 π 2 ([6]; Fig. 11)

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 9. Contours of the instantaneous vorticity around the flexible filament for the case of 10.5; (d) 10.75. T is the period.

Fig. 10. Time traces of the drag and lift coefficients for the case of

B ρ = 150 and K 11 = 0.15 at Re = 200: (a) t / T = 10; (b) 10.25; (c)

B ρ = 150 and K 11 = 0.15 at Re = 200:

Fig. 11. Critical density ratio of the flexible filament for Re = 1000 and h/ L = 0.01: 2, present study (stable state).

539

, Eq. (22);

, Eq. (25).

B ρcr = 1.3 Re−0.5 /(h/ L ) + 4K 11 π 2 [6]: ", present study (flapping state);

540

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 12. Time traces of the trailing-edge transverse locations of two flexible filaments (ρ = 5 and Re = 1000): Also shown in this figure are the contours of the instantaneous vorticity at t = 14.

B , K 11 = 0.01;

B , K 11 = 0.1.

Table 1 Convergence tests for flow over a flexible filament at Re = 1000 (O, converge; X, diverge). B K 11

ρ (CFLmax = 1)

0.1 1 5 10 100

(CFLmax = 0.6)

1.01

1.02

1.05

1.1

1.5

2

10

100

0.8

0.9

1.0

10

100

X X X X X

O X X X X

O X X X X

O O X X X

O O X X X

O O O O O

O O O O O

O O O O O

X X X X X

O O O O O

O O O O O

O O O O O

O O O O O

from a linear stability analysis. To test our numerical method, we conduct simulations for Re = 1000 and Fr = 0 with B different values of the parameters (ρ and K 11 ) corresponding to the flapping and stable states. The initial inclined angle between the filament and x-axis is θ = arctan(0.1). The computation is performed at the maximum CFL number of 1.0. In Fig. 11, we display the present results as square (stable state) and circle (flapping state) symbols. Fig. 12 shows the time B traces of the trailing-edge transverse locations of two cases of flexible filaments for flapping (ρ = 5, K 11 = 0.01) and stable B (ρ = 5, K 11 = 0.1) states, together with the instantaneous flow fields at t = 14. As shown, the flow becomes steady and the B trailing-edge location goes to y = 0 for the case of ρ = 5 and K 11 = 0.1, whereas the flexible filament sustains the flapping B state for the case of ρ = 5 and K 11 = 0.01. These results agree very well with the criterion suggested by [6]. Here, we want to discuss the numerical stability of the present numerical method for light (low ρ ) and stiff (high B K 11 ) structures. In Table 1, we provide the convergence/divergence of the solution according to the maximum CFL number B (CFLmax ) applied during the computations for different magnitudes of ρ and K 11 . As shown, with CFLmax = 1, the solution B B always converges for any K 11 ’s when ρ ≥ 2, but diverges for any K 11 ’s when ρ ≤ 1.01. In between these two values (i.e., B B 1.01 < ρ < 2), the present numerical scheme is stable for low K 11 (= 0.1 or 1) but unstable for higher K 11 ’s. On the other hand, with CFLmax = 0.6, the present numerical scheme is stable even for very low values of ρ = 0.9 but unstable for B ρ = 0.8, regardless of the magnitude of K 11 . The reason why the present numerical method is stable for the simulation of light and stiff structures is owing to the strong coupling between the fluid and structure equations used here [41,42]. However, as shown in Table 1, the numerical convergence depends on the maximum CFL number applied, indicating that the implicit treatment of the convection term in the Navier–Stokes equations may further improve the numerical convergence. 3.2. Three-dimensional flapping flag in a free stream We conduct a three-dimensional simulation of flow around a flapping flag in a free stream. Fig. 13 shows the schematic diagram of the computational domain and boundary conditions. The flag has the length L and width W , and its ratio is kept to be W / L = 1. The center of the leading edge of the flag is positioned at the origin. Initially, the flag shape is a flat plate inclined at θ = 0.1π from the xz-plane, where x and z are the free stream (streamwise) and vertical directions. The flag is pinned (not clamped) at the leading edge and freely moves elsewhere. The computational domain size is [− L , 7L ] ×

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

541

Fig. 13. Schematic diagram of the computational domain and boundary conditions for simulation of flow around a flapping flag in a free stream.

[−4L , 4L ] × [− L , L ] in x, y (spanwise) and z directions, and the corresponding number of grid points is 353 × 353 × 193. The smallest grid spacings are x =  y = 0.005L and  z = 0.01L. Uniform grids with these smallest grid spacings are distributed on and near the path of the flag. The Dirichlet boundary condition (u = U , v = w = 0) is applied at the inflow, top and bottom, and a convective boundary condition (∂ u i /∂ t + c ∂ u i /∂ x = 0) is given at the outflow boundary, where c is the velocity averaged over the outflow boundary. The periodic boundary condition is applied in the z direction. B T T T T The dimensionless material parameters are ρ = 100, K ab = 0.01 (for a, b = 1, 2), K 11 = K 22 = 2500 and K 12 = K 21 = 100. Here, the density ratio and bending and twisting coefficients are the same as those used in Huang and Sung [34], and the tension and shearing coefficients are taken to be sufficiently large to resist the extension or in-plane shear for satisfying the inextensibility condition. The ratio of thickness to the length of the flag is h/ L = 0.01. The number of blocks is N b = 50(s1 ) × 50(s2 ) (Fig. 3(a)). The gravity force is neglected (Fr = 0) for the comparison of the result with that of [34]. The simply supported boundary condition is used at the pinned boundary ( X |s1 =0 = X 0 , ∂ 2 X /∂ s21 |s1 =0 = 0), and the free-end boundary condition is used at the freely-moving boundaries (∂ X /∂ s1 · ∂ X /∂ s1 |s1 = L = 1, ∂ X /∂ s2 · ∂ X /∂ s2 |s2 =0, W = 1, ∂ X /∂ sa · ∂ X /∂ sb |s1 =L&s2 =0,W = 0 for a = b, ∂ 2 X /∂ s21 |s1 =L = 0, ∂ 3 X /∂ s31 |s1 =L = 0, ∂ 2 X /∂ s22 |s2 =0,W = 0, ∂ 3 X /∂ s32 |s2 =0,W = 0, ∂ 2 X /∂ sa ∂ sb |s1 =L&s2 =0,W = 0 for a = b). Computations are conducted for Re = 100 and 200 at the maximum CFL number of 1.0. Fig. 14 shows the time traces of the trailing-edge (i.e. at s1 = L) transverse locations of the flapping flag for Re = 100 and 200 in the case of zero gravity, together with those of [34]. The agreements between the results of present study and [34] are very good. At Re = 100, the amplitude of middle trailing edge (B) is larger than those of two corner trailing edges (A and C), and there is a phase difference between B and A (or C) due to the bendable property (Fig. 14(a)). The non-dimensional flapping frequency and amplitude of middle trailing edge (B) at this Reynolds number are f = 0.241 and A m = 0.309, respectively. Also, the non-dimensional frequency based on the flapping frequency and amplitude is St = 2 A m f = 0.149, which is also in good agreement with that of [34]. At a higher Reynolds number of Re = 200 (Fig. 14(b)), the transverse locations of A, B and C are almost identical, and the flapping frequency and amplitude are higher than those of Re = 100 ( f = 0.265, A m = 0.376 and St = 0.2 for Re = 200). This result is also in good agreement with that of [34]. Fig. 15 shows the instantaneous vortical structures around the flapping flag at four different time instants during one cycle for Re = 200, where the vortical structures are identified by the iso-surface of λ2 = −0.2 [44]. As shown, hairpin-like vortices are formed and shed from the trailing edge at each flapping. The time traces of the drag and lift coefficients for Re = 100 and 200 are shown in Fig. 16, where C d = 2D/ρ U 2 A f , C l = 2L/ρ U 2 A f , and A f is the flag area. The flapping frequency and amplitude at Re = 200 are higher and larger, respectively, than those at Re = 100. Fig. 17 shows the time traces of trailing-edge (i.e. at s1 = L) transverse locations of the flapping flag for Re = 200 in the case of Fr = 0.05 (under gravity). Unlike the case without gravity, phase differences are manifest among three trailing-edge locations (A, B and C in Fig. 13). As shown, the lower trailing edge fluctuates more than the upper one due to gravity. Fig. 18 shows the instantaneous vortical structures around the flapping flag at four different instants, together with the flag configurations. The flag shape changes significantly from the fluid–structure interaction under gravity. 3.3. Normal hovering of flexible flapping wing In this subsection, we conduct a simulation of normal hovering of a two-dimensional flexible flapping wing and compare the result with that of [31]. As shown in Fig. 19, a two-dimensional flexible wing with the chord length L moves back and forth without incoming fluid flow. Its translational and rotational motions are a priori given only at the leading edge

542

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 14. Time traces of the trailing-edge transverse locations of the flapping flag (Fr = 0): (a) Re = 100; (b) Re = 200. middle (B); , lower (C) trailing edges; 2 (A and C) and 1 (B), Huang and Sung [34].

, upper (A in Fig. 13);

,

Fig. 15. Vortical structures around the flapping flag for Re = 200 (Fr = 0): (a) t = 15.93; (b) 17.08; (c) 17.81; (d) 18.96. The blue-colored plate denotes the flag. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and are determined elsewhere through the fluid–structure interaction. The leading-edge translational motions are given as xl (t ) = ( A 0 /2) cos(2π f t ), yl (t ) = 0.5, and α (t ) = α0 + β sin(2π f t ), where A 0 = 2.5L, α0 = π /2 and β = −π /8. The computational domain size is [−10L , 10L ] × [−17.5L , 17.5L ] and the number of grid points is 481 × 609 in the horizontal (x) and vertical ( y) directions. The smallest grid spacing is x =  y = 0.01L and uniform grids with this smallest grid spacing are distributed on and near the path of the flexible wing. The Dirichlet boundary condition (u = v = 0) is applied at the whole outside boundaries. The values of dimensionless material and numerical parameters are given as ρ = 16.67,

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 16. Time traces of the drag and lift coefficients (Fr = 0):

, Re = 100;

Fig. 17. Time traces of the trailing-edge transverse locations of the flapping flag for Re = 200 (Fr = 0.05): (B); , lower (C) trailing edges.

543

, Re = 200.

, upper (A in Fig. 13);

, middle

Fig. 18. Vortical structures around the flapping flag for Re = 200 (Fr = 0.05): (a) t = 15.8; (b) 16.6; (c) 18.0; (d) 18.8. The blue-colored plate denotes the flag. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

544

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

Fig. 19. Schematic diagram of the computational domain and boundary conditions for simulation of normal hovering of a flexible flapping wing.

Fig. 20. Time trace of the trailing-edge angle:

, flexible wing (αt .e. );

Fig. 21. Time trace of the lift coefficient of flexible wing:

, rigid wing (α ). T is the flapping period.

, present study; 1, Yin and Luo [31]. T is the flapping period.

B T K 11 = 5.383, K 11 = 833, N b = 100, and h/ L = 0.06. The values of these parameters are determined to match those used in [31]. The prescribed boundary condition is used at the leading edge ( X |s1 =0 = (xl , yl ), ∂ X /∂ s1 |s1 =0 = (cos α (t ), sin α (t ))), and the free-end boundary condition is applied at the trailing edge (∂ X /∂ s1 · ∂ X /∂ s1 |s1 = L = 1, ∂ 2 X /∂ s21 |s1 = L = 0, ∂ 3 X /∂ s31 |s1 = L = 0). The Reynolds number is Re = ρ U m L /ν = 150 based on the maximum velocity of leading edge U m and chord length L. The maximum CFL number is taken to be 0.75. Fig. 20 shows the time trace of the trailing-edge angle of flexible wing (αt .e. ), together with that of rigid one (α ). The trailing-edge angle of flexible wing is much larger than that of rigid wing due to the bending properties, and the trailing-edge angle lags behind that of the rigid wing. Fig. 21 shows the time trace of the lift coefficient of flexible wing, 2 together with that of [31], where the lift coefficients are C l = 2L /ρ U m L and U m = π A 0 f . The present result is in an excellent agreement with that of [31].

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

545

4. Summary In the present study, an IB method was developed for the simulation of flow around an elastic slender body. The present method was based on the discrete-forcing IB method for a stationary body proposed by Kim, Kim and Choi [25] and was fully coupled with the dynamic equation for an elastic slender body motion. The discrete-forcing approach was used to relieve the limitation on the size of the computational time step. The incompressible Navier–Stokes equations were solved in the Eulerian coordinate using a semi-implicit fractional step method. On the other hand, the dynamic equation for an elastic slender body was solved in the Lagrangian coordinate. To describe the body motion, the elastic slender body was modeled as a thin flexible solid and was segmented by a finite number of thin blocks, and a pseudo-normal vector was calculated at the center of each block to construct the thin body surface. Each block was moved by external and internal forces such as the hydrodynamic, elastic and buoyancy forces. The hydrodynamic force was obtained directly by integrating the Navier–Stokes equations without any additional interpolation. All the spatial derivative terms were discretized with the second-order central difference scheme. In this paper, we suggested a discrete-forcing IB method together with strong coupling between the fluid and structure equations. This numerical method enabled us to simulate flow interacting with light and stiff structures at maximum CFL numbers of 0.75–1, but strong coupling required 3–5 overall iterations per time step. Although weak coupling does not require overall iteration, it may have serious numerical instability problems for light and stiff structures. Therefore, even if the present numerical method requires a few overall iterations per time step, it is numerically efficient owing to the use of high CFL number. With this method, we simulated three different fluid–structure interaction problems: flows around a flexible filament, a flapping flag in a free stream, and a flexible flapping wing in normal hovering, respectively. The results obtained agreed well with those from previous studies, indicating the accuracy and efficiency of the present immersed boundary method. Acknowledgements This research is supported by the NRF Programs (NRF-2011-0028032, NRF-2012M2A8A4055647) of MSIP, Korea. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

M. Shelley, N. Vandenberghe, J. Zhang, Heavy flags undergo spontaneous oscillations in flowing water, Phys. Rev. Lett. 94 (2005) 094302. T.L. Daniel, S.A. Combes, Flexible wings and fins: bending by inertial or fluid-dynamic forces? Integr. Comp. Biol. 42 (2002) 1044–1049. M.J. Shelley, J. Zhang, Flapping and bending bodies interacting with fluid flows, Annu. Rev. Fluid Mech. 43 (2011) 449–465. J.J. Allen, A.J. Smits, Energy harvesting eel, J. Fluids Struct. 15 (2001) 629–640. W. Shyy, Y. Lian, J. Tang, H. Liu, P. Trizila, B. Stanford, L. Bernal, C. Cesnik, P. Friedmann, P. Ifju, Computational aerodynamics of low Reynolds number plunging, pitching and flexible wings for MAV applications, Acta Mech. Sin. 24 (2008) 351–373. B.S.H. Connell, D.K.P. Yue, Flapping dynamics of a flag in a uniform stream, J. Fluid Mech. 581 (2007) 33–67. B.S.H. Connell, Numerical investigation of the flow-body interaction of thin flexible foils and ambient flow, Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA, 2006. T. Sawada, T. Hisada, Fluid–structure interaction analysis of the two-dimensional flag-in-wind problem by an interface-tracking ALE finite element method, Comput. Fluids 36 (2007) 136–146. K. Namkoong, H.G. Choi, J.Y. Yoo, Computation of dynamic fluid–structure interaction in two-dimensional laminar flows using combined formulation, J. Fluids Struct. 20 (2005) 51–69. C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252–271. J. Lee, J. Kim, H. Choi, K.-S. Yang, Sources of spurious force oscillations from an immersed boundary method for moving-body problems, J. Comput. Phys. 230 (2011) 2677–2695. J.H. Seo, R. Mittal, A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations, J. Comput. Phys. 230 (2011) 7347–7363. L.J. Fauci, C.S. Peskin, A computational model of aquatic animal locomotion, J. Comput. Phys. 77 (1988) 85–108. E. Jung, C.S. Peskin, Two-dimensional simulations of valveless pumping using the immersed boundary method, SIAM J. Sci. Comput. 23 (1) (2001) 19–45. Y. Kim, C.S. Peskin, 2-D parachute simulation by the immersed boundary method, SIAM J. Sci. Comput. 28 (6) (2006) 2294–2312. L. Zhu, C.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. Y. Kim, C.S. Peskin, Penalty immersed boundary method for an elastic boundary with mass, Phys. Fluids 19 (2007) 053103. Y. Kim, C.S. Peskin, Numerical study of incompressible fluid dynamics with nonuniform density by the immersed boundary method, Phys. Fluids 20 (2008) 062101. Y. Kim, M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty immersed boundary method, J. Comput. Phys. 229 (2010) 4840–4853. J. Yang, E. Balaras, An embedded-boundary formation for large-eddy simulation of turbulent flows interacting with moving boundaries, J. Comput. Phys. 215 (2006) 12–40. W.-X. Huang, S.J. Shin, H.J. Sung, Simulation of flexible filaments in a uniform flow by the immersed boundary method, J. Comput. Phys. 226 (2007) 2206–2228. D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys. 105 (1993) 354–366. W.-X. Huang, H.J. Sung, An immersed boundary method for fluid–flexible structure interaction, Comput. Methods Appl. Mech. Eng. 198 (2009) 2650–2661. R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, J. Comput. Phys. 171 (2001) 132–150.

546

I. Lee, H. Choi / Journal of Computational Physics 280 (2015) 529–546

[26] J. Lee, D. You, An implicit ghost-cell immersed boundary method for simulations of moving body problems with control of spurious force oscillations, J. Comput. Phys. 233 (2013) 295–314. [27] H. Luo, R. Mittal, X. Zheng, S.A. Bielamowicz, R.J. Walsh, J.K. Hahn, An immersed-boundary method for flow-structure interaction in biological systems with application to phonation, J. Comput. Phys. 227 (2008) 9303–9332. [28] R. Bhardwaj, R. Mittal, Benchmarking a coupled immersed-boundary-finite-element solver for large-scale flow-induced deformation, AIAA J. 50 (2012) 1638–1642. [29] J. Lee, D. You, Study of vortex-shedding-induced vibration of a flexible splitter plate behind a cylinder, Phys. Fluids 25 (2013) 110811. [30] H. Luo, B. Yin, H. Dai, J.F. Doyle, A 3D computational study of the flow-structure interaction in flapping flight, AIAA Paper 2010-556, 2010. [31] B. Yin, H. Luo, Effect of wing inertia on hovering performance of flexible flapping wings, Phys. Fluids 22 (2010) 111902. [32] H. Dai, H. Luo, J.F. Doyle, Dynamic pitching of an elastic rectangular wing in hovering motion, J. Fluid Mech. 693 (2012) 473–499. [33] F.-B. Tian, H. Dai, H. Luo, J.F. Doyle, B. Rousseau, Fluid–structure interaction involving large deformations: 3D simulations and applications to biological systems, J. Comput. Phys. 258 (2014) 451–469. [34] W.-X. Huang, H.J. Sung, Three-dimensional simulation of a flapping flag in a uniform flow, J. Fluid Mech. 653 (2010) 301–336. [35] S. Kang, G. Iaccarino, F. Ham, P. Moin, Prediction of wall-pressure fluctuation in turbulent flows with an immersed boundary method, J. Comput. Phys. 228 (2009) 6753–6772. [36] W.-X. Huang, H.J. Sung, Improvement of mass source/sink for an immersed boundary method, Int. J. Numer. Methods Fluids 53 (2007) 1659–1671. [37] H. Gouraud, Continuous shading of curved surfaces, IEEE Trans. Comput. C-20 (6) (1971) 623–629. [38] G. Thürrner, C.A. Wüthrich, Computing vertex normals from polygonal facets, J. Graph. Tools 3 (1) (1998) 43–46. [39] N. Max, Weights for computing vertex normals from facet normals, J. Graph. Tools 4 (2) (1999) 1–6. [40] J. Du, B. Fix, J. Glimm, X. Jia, X. Li, Y. Li, L. Wu, A simple package for front tracking, J. Comput. Phys. 213 (2006) 613–628. [41] J. Yang, S. Preidikman, E. Balaras, A strongly coupled, embedded-boundary method for fluid–structure interactions of elastically mounted rigid bodies, J. Fluids Struct. 24 (2008) 167–182. [42] J. Yang, F. Stern, A simple and efficient direct forcing immersed boundary framework for fluid–structure interactions, J. Comput. Phys. 231 (2012) 5029–5061. [43] D. Kim, H. Choi, Immersed boundary method for flow around an arbitrarily moving body, J. Comput. Phys. 212 (2006) 662–680. [44] J. Jeong, F. Hussain, On the identification of a vortex, J. Fluid Mech. 285 (1995) 69–94.