An immersed boundary technique for simulating complex flows with rigid boundary

(1)

An immersed boundary technique for simulating complex ﬂows

with rigid boundary

Shen-Wei Su

a

, Ming-Chih Lai

b

, Chao-An Lin

a,*

a_{Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 300, Taiwan} b_{Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan}

Received 7 March 2005; received in revised form 8 July 2005; accepted 29 September 2005 Available online 19 January 2006

Abstract

A new immersed boundary (IB) technique for the simulation of flow interacting with solid boundary is presented. The present for-mulation employs a mixture of Eulerian and Lagrangian variables, where the solid boundary is represented by discrete Lagrangian mark-ers embedding in and exerting forces to the Eulerian fluid domain. The interactions between the Lagrangian markmark-ers and the fluid variables are linked by a simple discretized delta function. The numerical integration is based on a second-order fractional step method under the staggered grid spatial framework. Based on the direct momentum forcing on the Eulerian grids, a new force formulation on the Lagrangian marker is proposed, which ensures the satisfaction of the no-slip boundary condition on the immersed boundary in the inter-mediate time step. This forcing procedure involves solving a banded linear system of equations whose unknowns consist of the boundary forces on the Lagrangian markers; thus, the order of the unknowns is one-dimensional lower than the fluid variables. Numerical exper-iments show that the stability limit is not altered by the proposed force formulation, though the second-order accuracy of the adopted numerical scheme is degraded to 1.5 order. Four different test problems are simulated using the present technique (rotating ring flow, lid-driven cavity and flows over a stationary cylinder and an in-line oscillating cylinder), and the results are compared with previous experimental and numerical results. The numerical evidences show the accuracy and the capability of the proposed method for solving complex geometry flow problems both with stationary and moving boundaries.

1. Introduction

The fluid–solid interaction problems are frequently encountered in many engineering applications. When the immersed object is complex and moving, the problem poses two difficulties; namely, the accuracy of the fluid computa-tion and the treatment of complex boundary. This accuracy of fluid calculations can be improved by employing high resolution schemes[1–3]or using adaptive or moving mesh techniques such as in[4–7]to resolve the small scale struc-tures. The other difficulty arises from the handling of the interaction between the fluid and the immersed body which is the major issue of this work. Traditionally, the body-fitted or unstructured grid methods are adopted to simulate

the ﬂow with complex rigid boundary. However, the computational cost and memory requirements of these methods are generally high. Also, the modeling of complex time evolving moving boundary ﬂows requiring transient re-meshing strategies increases further the computational overhead of these methods.

Another simple alternative to tackle the complex fluid– structure problem is to use the Cartesian grid instead of body-fitted grid. Since the boundary now does not conform with the Cartesian grid, the main difficulties are to model the complex immersed boundary accurately and at the same time maintain the solution efficiency discretized on the Cartesian grids. One approach, the so-called Cartesian grid method [8], tracks the boundary by identifying the control volumes (cut cell) in the Cartesian grid that are cut by the immersed boundary and determines the intersec-tion of the boundary with the faces of these cut cells. More

*

Corresponding author.

E-mail address:[email protected](C.-A. Lin).

(2)

specifically, a control volume containing the boundary is reshaped by discarding the region that lies within the solid and this results in the trapezoidal shape control volume. Due to the irregular shapes of the cut cells, complex inter-polating procedures to approximate the fluxes must be introduced and this unavoidably affects the computational efficiency of the Cartesian solvers.

Instead of using the cut cell approach, one approach also within the Cartesian grid framework, the immersed interface method (IIM) has attracted attention recently for the simulation of ﬂuid ﬂows in irregular regions. This method was originally developed by LeVeque and Li [9]

for solving elliptic equations with discontinuous coeﬃ-cients and singular sources along an arbitrary interface. A derived jump conditions are employed near the bound-ary, so that the sharp solutions can be achieved. The imple-mentation of the IIM has been applied within the vorticity and stream function formulation [10,11]. Also, based on the work of Li and Lai [12] for the Navier–Stokes equa-tions with singular forces, the immersed interface method using the Navier Stokes equation with primitive variables

[13,14] has also been explored. A recent review on the immersed interface method can be found in[15]by Li.

The complex geometry within the Cartesian grid can be simulated by generating external force field to mimic the immersed boundary. The immersed boundary (IB) method proposed by Peskin[16,17], has been applied successfully to blood–valve interaction and other biological problems. The proposed immersed boundary formulation employs a mixture of Eulerian and Lagrangian variables, where the immersed boundary is represented by a set of discrete Lagrangian markers embedding in the Eulerian fluid domain. Those markers can be treated as force generators to the fluid, and meanwhile move along with the fluid. The interaction between the Lagrangian markers and the fluid variables defined on the fixed Eulerian grid is linked by a well-chosen discretized delta function. Although the IB method is developed to handle mostly the fluid problem with elastic structures, it has also been used to simulate the flow with rigid boundaries or structures. Lai and Peskin[18]

proposed a new formally second-order accurate immersed boundary method to simulate flows past a rigid cylinder. The rigid boundary is simply mimicked by a stiff spring con-necting the boundary points to their target positions and the boundary force is generated by the deviation of the markers from their target positions. Since this feedback force is com-puted and distributed into the fluid at the beginning of each time step, the problem becomes very stiff which has the con-sequence of the small time step. Besides, the force generates oscillations along the surface of the boundary.

Also adopting the mixture of Eulerian and Lagrangian variables, Goldstein et al. [19,20] proposed the virtual boundary method to simulate the flow with solid boundary within the spectral method framework. The solid boundary is treated as a force generator where the force field is calcu-lated by a feedback method based on the difference between the predicted velocity and the actual velocity of

the boundary. However, in order to prevent the generation of the spurious oscillations in simulating the start-up flow around a cylinder, the magnitude of the CFL number has been kept below O(102). The feedback force of Goldstein et al. was also adopted by Saiki and Biringen [21] using the finite difference formulation. The distribution of the boundary force to the grid points was achieved by the area-weighted average function. This same function is used to interpolate the fluid velocity to the boundary points. A very good agreement has been found between their numer-ical computations and experimental results for the simula-tion of flow around a cylinder. However, the time step restriction of the method is similar to the scheme proposed by Goldstein et al.[19]and an excessive number of surface points O(1000) were used to model the immersed boundary even for two-dimensional flows.

In [23], Silva et al. calculated the momentum forcing along the immersed boundary using the pressure and veloc-ity derivatives interpolated by a second-order Lagrange polynomial approximation. Once the momentum forcing at the boundary points has been calculated, they are dis-tributed over the Eulerian grid by the discrete delta func-tion. The whole numerical scheme is very similar to the immersed boundary method of Lai and Peskin except the evaluation of the forcing term. While Silva et al.’s approach is ideally simple, the calculations of momentum forcing at the boundary points are quite complex.

Without adopting the Lagrangian markers, Mohd-Yusof [24] proposed a direct forcing method within the spectral framework, where direct momentum forcing is applied to a set of points adjacent to the surface and inte-rior to the body, which is equivalent to the direct imposi-tion of the velocity boundary condiimposi-tions on the Eulerian grids to ensure the no-slip condition at the immersed boundary. Therefore, information regarding the locations of the Eulerian grids either external or internal to the immersed boundary must be determined. The major advantage of the method is that since the force field is directly computed from the momentum equations, so that the time step can be larger than the previous methods. Fadlun et al. [25] further extended the Mohd-Yusof approach to a finite-difference formulation on a stagger gird system, where direct forcing or velocity boundary con-dition is applied along the first Eulerian grid external to the immersed boundary. Although this velocity boundary con-dition is obtained by linear interpolation procedure, the choice of interpolation direction is arbitrary.

Based on the concept of direct forcing, Kim et al. [26]

introduced both the momentum forcing and mass source/ sink to properly represent the immersed body. The momen-tum forcing and the mass source/sink are applied only on the body surface or inside the body so that the no-slip boundary condition on the immersed boundary and the continuity for the cell containing the boundary are both sat-isﬁed. Since the immersed boundary in general does not coincide with the grid points, an interpolation scheme for computing the momentum forcing must be employed.

(3)

Tseng and Ferziger[27]extended the idea of Fadlun et al. via a ghost cell approach. The immersed boundary is repre-sented by piecewise linear segments and the ghost cells are defined to lie just inside the body but adjacent to computa-tional grids of the fluid domain. The values of fluid variables at those ghost cells are obtained by extrapolation using a local quadratic scheme which involves the neighboring flow nodes and the associated velocity boundary condition.

It can be summarized that the major drawback of the existing Eulerian Lagrangian based IB techniques to simulate solid boundary flows is the restriction of small CFL number due to the spring or feedback force formula-tions employed, despite the simple interpolation procedure used to link the Lagrangian marker and the Eulerian grids. On the other hand, though the existing direct forcing approaches can compute flows using a larger time step, the determinations of the forcing locations and their mag-nitudes may not be straightforward on the Eulerian grids, especially for time evolving moving boundary flows. Com-bining the merits of the above two approaches, a new immersed boundary (IB) technique for the simulation of flows interacting with solid boundary is proposed.

The present formulation employs a mixture of Eulerian and Lagrangian variables, where the solid boundary is represented by discrete Lagrangian markers embedding in and exerting forces to the Eulerian fluid domain. Based on the direct momentum forcing on the Eulerian grids, a new force formulation on the Lagrangian marker is pro-posed, which ensures the satisfaction of the no-slip bound-ary condition on the immersed boundbound-ary. The interactions between the Lagrangian markers and the fluid variables on the fixed Eulerian grid are linked by a simple discrete delta function. The boundary forces are first computed on the Lagrangian markers and then distributed to the Eulerian grids using discrete delta function. As will be shown in the next section that the present formulation does not need to employ any modifications due to the jump conditions such as in the immersed interface method, so the present implementation is simpler but can be potentially less accurate.

Four different test problems are simulated using the present technique (rotating ring flow, lid-driven cavity and flows over a stationary cylinder and an in-line oscillat-ing cylinder). These cases are used to examine the numeri-cal accuracy of the method and its capability to model complex flows including the time evolving moving bound-ary. The results are compared with previous experimental and numerical results to assess the accuracy and the capa-bility of the proposed method for solving complex geome-try flow problems.

2. The methodology of immersed boundary technique 2.1. Mathematical formulation

Here, consider a problem of a viscous incompressible ﬂuid in a two-dimensional square domain X = [0, L]·

[0, L] containing an immersed massless boundary in the form of a simple closed curve C, as shown in Fig. 1. The lowercase and uppercase letters are used to represent the variables deﬁned on the Cartesian grid and the Lagrangian markers, respectively. The immersed boundary is tracked by the parametric form, X(s), 0 6 s 6 Lb, where s is the

parameter of the reference configuration of the boundary. (For simplicity, the boundary configuration is assumed to be fixed in time.) As mentioned before, the influence of the immersed boundary to the fluid is represented by forces that exert to the fluid so that the fluid will move by the pre-scribed velocity U(X(s), t) (no-slip condition) at the immersed boundary. Thus, the governing equations of this fluid–structure interaction system are

ou otþ rðuuÞ ¼ rp þ 1 Rer 2_u_{þ f;} _ð1Þ r u ¼ 0; ð2Þ fðx; tÞ ¼ Z Lb 0 FðXðsÞ; tÞdðx XðsÞÞ ds; ð3Þ UðXðsÞ; tÞ ¼ Z X uðx; tÞdðx XðsÞÞ dx. ð4Þ

Here, x = (x, y), U(x, t) is the fluid velocity, p(x, t) is the fluid pressure, and the number Re is the non-dimensional Reynolds number. Note that, the Eulerian force field f(x, t) is defined on the whole fluid domain X while the Lagrangian force F(X(s), t) is on the immersed boundary C. Eqs. (1) and (2) are the familiar Navier–Stokes equa-tions of a viscous incompressible fluid. Eqs.(3) and (4) rep-resent the interaction between the immersed boundary and the fluid. In particular, Eq.(3)describes that the force act-ing on the fluid f is due to the boundary force F alone, and Eq.(4)represents the fluid moves with the same prescribed velocity of the immersed boundary. One can easily see that the present formulation employs a mixture of Eulerian (x) and Lagrangian (X) variables which are linked by the

Ω

Γ

(4)

two-dimensional Dirac delta function d(x) = d(x)d(y). To close the system, the initial and the physical boundary (oX) conditions of velocity should be given as well.

The main diﬃculty of the above mathematical formula-tion is that the forcing term f is not known a priori (since the boundary force F is unknown).The force ﬁeld F, how-ever, can be determined by enforcing the no-slip boundary condition of the immersed boundary.

2.2. Numerical scheme

The present numerical scheme is a mixed Eulerian– Lagrangian finite difference method for simulating the complex fluid flows interacting with an immersed bound-ary. Therefore, two distinct discretized grids are intro-duced; the regular lattice points to cover the whole fluid domain, and the marker points to discretize the immersed boundary. Throughout this work, the spatial discretization makes use of the staggered marker-and-cell (MAC) mesh introduced by Harlow and Welsh[28]. Thus, the fluid vari-ables are defined at different locations of the lattice grids. For instance, the pressure is defined on the grid points labeled as x = (xi, yj) = ((i 1/2)h, (j 1/2)h) for i, j =

1, 2, . . . , N, the velocity components u and v are deﬁned at (xi1/2, yj) = ((i 1)h, (j 1/2)h) and (xi, yj1/2) = ((i 1/

2)h, (j 1)h), respectively, where the spacing h = Dx = Dy= L/N. On the other hand, the immersed boundary is discretized by the M discrete Lagrangian markers Xk= (Xk, Yk), where the marker spacing is Ds = Lb/M.

In general, the markers Xk do not coincide with

the Eulerian grid points x, so a discrete delta function linkage between those two grids must be employed. Here, the following discrete version of Dirac delta function is used:

dhðx XkÞ ¼ dhðxi XkÞdhðyj YkÞ; ð5Þ

where dh is the hat function[22]deﬁned by

dhðrÞ ¼

ð1 jrj=hÞ=h; for jrj h 0; otherwise.

ð6Þ

The adopted delta function is equivalent to the bilinear interpolation scheme, which is consistent with the interpo-lation procedure employed in the Navier–Stokes solver.

Based on the adopted grid arrangement, the governing equations(1)–(4) can be discretized directly. The integrals in Eqs. (3) and (4) are approximated by their Riemann sums and the Dirac delta function d is replaced by its dis-crete version dh in5. On the other hand, Eqs. (1) and (2)

are integrated using the fractional-step method proposed by Choi and Moin[29], and the non-linear convection term and the diﬀusion term are treated by the Adams–Bashforth and the Crank–Nicholson methods, respectively. All the spatial derivatives are approximated by the second-order central diﬀerence method. At the beginning of the time step, the solutions un1, unmust be given in order to march to un+1.

The time advancement and spatial discretization can be done in the following steps:

~ u un Dt ¼ 3 2rhðuuÞ n þ1 2rhðuuÞ n1 rhpnþ 1 2Rer 2 hðu n_{þ ~}_uÞ; _ð7Þ fðxÞ ¼X M k¼1 FðXkÞdhðx XkÞDs; for all x¼ ðxi; yjÞ ð8Þ u_~_u Dt ¼ f _ð9Þ UðXkÞ ¼ X x uðxÞdhðx XkÞh2; k¼ 1; 2; . . . ; M. ð10Þ u_u Dt ¼ rhp n_; _ð11Þ r2 hp nþ1 _¼rh u Dt ; ð12Þ unþ1_u Dt ¼ rhp nþ1_; _ð13Þ

where ~u, u*_{, and u}**_{are the intermediate velocity}

compo-nents between the time step n and n + 1. The discrete oper-ators $h and r2h represent the regular centered diﬀerence

approximations for the gradient and Laplace operators. Eq. (8) describes that the Eulerian forces are obtained from distributing the boundary forces at those Lagrangian markers to their neighboring Cartesian grids. Eq.(9) repre-sents the direct momentum forcing to modify the velocities on the Eulerian grids using the forces obtained from Eq.

(8). Eq. (10) says that the velocity at the markers can be interpolated from the neighboring grids and should equal to the prescribed boundary velocity. Note that, the summa-tionP_xin(10)is only over four neighboring grid points of Xk, since the support of the discrete delta function is as

wide as the mesh width h. Also, the above approximation in (10) is nothing but the bi-linear interpolation of the velocity. The continuity condition is satisﬁed by solving the pressure Poisson equation(12).

It should be noted the discrete momentum forcing in the present algorithm seems to be calculated by the formula f Dtr2

hf

_{=ð2ReÞ rather than f}*_{. This is due to the implicit}

treatment of the intermediate velocity in the prediction step

(7). This confusion can be avoided by simply using the explicit Adams–Bashforth discretization for the diffusion term but it will lead to a restrictive time step constraint when Reynolds number is small. Despite the above momentum forcing discrepancy, the physical quantities such as the drag and lift coefficients (as shown in the numerical results) will not be affected. This is because, in the present formulation, the lift and drag are calculated by the discrete sum of momentum forces over the whole fluid grid points and the discrete Laplacian term makes no contribution in the summation. This can be derived eas-ily as follows. (Note that, f*_{is non-zero only in the}

(5)

X i X j r2 hf ¼X j X i f iþ1;j fi;j h2 f i;j fi1;j h2 þX i X j f_i;jþ1 f i;j h2 f_i;j f i;j1 h2 ¼ 0.

The above is verified numerically by computing flows over a cylinder employing the Adams–Bashforth and Crank– Nicholson schemes, respectively, and identical lift and drag coefficients are obtained.

In the present numerical setting, the support of f*_{is only}

in the neighborhood of the immersed boundary C. Also, the force ﬁeld is computed in the intermediate step, so that the solutions of coupling linear systems can be avoided. The most time-consuming part of the above scheme is to solve three Helmholtz-type equations (two for velocity in

(7) and one for the pressure in (12)) which can be solved by the eﬃcient Bi-CGSTAB method[30].

The remaining question is how to ﬁnd the boundary force ﬁeld F*_{properly in Eq.}₍₈₎_{, so that the resulting}

veloc-ity can be adjusted by the presence of the forces. More pre-cisely, the present forcing ﬁeld will bring the intermediate velocity u* _{to the desired velocity U}* _{at the immersed}

boundary. It is important to mention that the present scheme (7)–(13) is very similar to the direct momentum forcing of Mohd-Yusof used in[24,25] except in Eq. (8). Since the present formulation employs a mixture of Eule-rian and Lagrangian variables, the boundary forces are ﬁrst computed on the Lagrangian markers and then are distributed to the Eulerian grids using discrete delta func-tion, which depends on the distances between the Eulerian grids and the Lagrangian markers. Therefore, there is no need to determine the relative locations of the Eulerian grid either internal or external to the immersed boundary, and this further simpliﬁes the computational procedure. 2.3. The boundary force derivation

Here, description as how to ﬁnd the marker forces F*

properly is introduced, so that the intermediate velocity u*_{at Lagrangian markers will satisfy the desired boundary}

values U*_{. Firstly, the velocity ~}_{u calculated in Eq.} ₍₇₎ _is

interpolated to the markers to obtain the marker velocity e UðXkÞ, as e UðXkÞ ¼ X x ~ uðxÞdhðx XkÞh2. ð14Þ

If the above interpolation procedure is applied to Eq. (9)

(i.e. direct momentum forcing on Eulerian grids) directly, then X x fðxÞdhðx XkÞh2¼ UðXkÞ eUðXkÞ Dt ; ð15Þ where U*_(X

k) is the interpolating velocity of u*at the

mar-ker Xk. By simply letting the marker velocity U*(Xk) =

Un+1(Xk), the force ﬁeld f* will bring the velocity u* in

the intermediate step of the Navier–Stokes solver to the prescribed values at the immersed boundary.

Note that, it is the intermediate velocity U*_{rather than}

the ﬁnal velocity Un+1 satisfying the prescribed boundary velocity, so that the prediction step (Eq.7) can be decou-pled from the projection step (Eqs. (12) and (13)) com-pletely in the Navier–Stokes integration. It is possible to let the ﬁnal velocity satisfy the prescribed boundary condi-tion by doing a few iteracondi-tions at each time step, but this will complicate the present implementation.

Combining Eqs.(8) and (15), the Lagrangian boundary force can be obtained as

X x XM j¼1F _ðX jÞdhðx XjÞDs |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} fðxÞ dhðx XkÞh2 ¼U nþ1_ðX kÞ eUðXkÞ Dt . ð16Þ Therefore XM j¼1 X x dhðx XjÞdhðx XkÞDsh2 ! FðXjÞ ¼U nþ1_ðX kÞ eUðXkÞ Dt . ð17Þ

There are M equations for the M Lagrangian marker forces F*_(X

k), k = 1, 2, . . . , M. Since the support of the discrete

delta function is the same as the grid spacing h, the above coeﬃcient matrix in(17) is banded. (For a closed bound-ary, it is a periodic banded matrix.) The matrix can become explicit if dh(x Xj) = 0, when j 5 k. This only happens in

regions where the Eulerian grids inﬂuenced by Lagrangian marker force F*_(X

k) are not aﬀected by Lagrangian marker

forces F*_(X

j), where j = 1, M and j 5 k. For example, the

Lagrangian markers may collocate with the Eulerian grids or reside on the edges of the Eulerian grids. While the latter may seem plausible, the restriction is that the grid points connected to the edge should not be aﬀected by other mar-ker forces, which may not be trivial and straightforward.

For simplicity, the markers are generally uniformly dis-tributed along the immersed boundary. As will be shown later, the optimum marker spacing (Ds) scales with the grid spacing h. Therefore, the number of markers (M) required is one dimensional lower than the fluid variables and is the same order of magnitude to the Eulerian grid points sur-rounding the immersed boundary either on the fluid or the solid sides. Also, for two dimensional flows, the banded matrix is tria- or penta-diagonal. Therefore, the cost of solving Eq.17is O(M), which is negligibly small compared to the cost of solving the momentum and the pressure pois-son equations. It should also be noted that a singular matrix may occur if the markers are too close to each other. Numerical experiments show that if the ratio of the marker spacing to the grid width Ds/h is greater than 0.5, no singu-lar matrix will be found, and the optimum marker spacing (Ds) scales with the grid spacing h. This further excludes the

(6)

possibility of the singular matrix and reduces the sizes of the banded force matrix.

2.4. The full solution procedure

A whole numerical procedure for each time step of the present method is summarized as follows.

(1) Solve Eq. (7) in order to obtain the intermediate velocity component ~uðxÞ.

(2) Calculate the Lagrangian boundary force F*_(X k) for

all markers using Eq.(17)and the prescribed marker velocity Un+1(Xk).

(3) Distribute the Lagrangian force F*_(X

k) to the

Eule-rian grid f*_{(x) using Eq.}₍₈₎_.

(4) Correct the intermediate velocity u*_{(x) using the}

newly obtained force f*_{(x) through Eq.}₍₉₎_.

(5) Obtain the intermediate velocity u**_{(x) using Eq.}₍₁₁₎_.

(6) Compute the pressure by solving the pressure Poisson equation(12).

(7) Update the ﬂuid velocity un+1(x) by Eq.(13).

3. Numerical results

3.1. Flow induced by a rotating ring

Here, the flow considered is a ring rotating at a constant angular velocity locating at the center of a square domain 1 6 x, y 6 1. The radius of the ring is 0.3, and the angular velocity X of the ring is 2, where the corresponding Rey-nolds number is 18. Since the square domain boundary is stationary, the flow is induced by the rotating ring. The velocity boundary conditions are imposed along the ring using the immersed boundary method. Since the immersed boundary does not coincide with the Cartesian grids, the forcing interpolation procedure discussed in the previous section is required to obtain the velocity correction near the immersed boundary. At the steady state, the flow inside the ring would undergo a solid body rotation, if the flow remains laminar. Therefore, this can provide a means of examining the numerical accuracy of the present scheme.

Simulations are conducted using four sets of grids, i.e. 40· 40, 80 · 80, 160 · 160 and 320 · 320. The grid spacing ratio Ds/h is 0.85, and the maximum CFL number is approximately 0.5. Fig. 2 shows the velocity vector plots of the ﬂow, where the rotating ring is marked with the bold line. The solid body rotation can be clearly observed inside the rotating ring. The L1and L2error norms are shown in

Fig. 3. The spatial accuracy of the present scheme is around 1.5 order, though at higher grid density the order of accu-racy does increase. Since this case is a discontinuous prob-lem using the simple bi-linear interpolation procedure, the order of accuracy achieved is rather encouraging.

The time of the immersed boundary calculation relative to the overall cost is also examined here using the rotating ring ﬂow as the test bed. The immersed boundary method

includes four parts, namely the distribution of the markers, the calculation of the interpolation function, the determi-nation of the maker forces and the distributions of marker forces to the Eulerian grids. The measurements based on the 320· 320 grid indicates that the immersed boundary method occupies slightly less than 1% of the overall putational cost, and this is rather small. Most of the com-putational cost is on the solution of the pressure Poisson equation, which is about 90% of the overall time. For stationary cylinder, the distribution of the markers and the calculation of the interpolation function are calculated only once, and this further reduces the cost of the IBM operation. -0.5 -0.25 0 0.25 0.5 X -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Y

Fig. 2. Velocity vector of ﬂow induced by a rotating ring within a square domain. 100 200 300 400 number of grid 10-3 10-2 er ro r n o rm L1 L₂ 2nd-order

(7)

3.2. Lid-driven cavity

The lid-driven cavity flow is widely used to verify the accuracy of numerical method too. Fig. 4shows the geo-metric layout of the lid-driven cavity within a square domain [1, 1] · [1, 1], where the cavity (jxj þ jyj ¼ 1=pffiffiffi2) is obliqued at 45 with respect to the Eulerian grid. Thus, the width of the cavity is 1 and the lid locates along the edge xþ y ¼ 1=pffiffiffi2 with constant driven velocity one. This inclined lid-driven cavity is mimicked by the present immersed boundary method, where the respective velocity distributions are assigned at Lagrangian markers. No-slip conditions are prescribed along the computational bound-ary of the Eulerian grids.

Three diﬀerent uniform grids (N· N, N = 40, 80, 160) are used in the simulations with the associated number of Lagrangian markers (M = 80, 160, 320), respectively. The grid spacing ratio of marker points and the lattice points Ds/h is 1. The Reynolds number is chosen as Re = 100 and the maximum CFL number is 0.4.

Fig. 5shows the steady velocity component U inside the cavity but along the line y = x. The coordinate Y = 0 rep-resents the intersection of the lower left edge and y = x.

Fig. 6 shows the steady velocity component V inside the cavity but along the line y =x. The coordinate X = 0 rep-resents the intersection of the upper left edge and y =x. One can see that present numerical results converge to the result obtained by Ghia et al.[31]quite well.

To demonstrate the capability of the present scheme to mimic flows with multiple objects, the simulation is applied to flow within the lid-driven cavity with five circular obsta-cles. The size of the computational domain is the same as the previous one, and the Reynolds number is 100 based on the lid-driven velocity and the height of the cavity.

The grid spacing ratio Ds/h is 0.85 within a 160· 160 grid.

Fig. 7 shows the vector plot of the simulation, where the original lid-driven cavity ﬂow structures have been modi-ﬁed by the presence of the obstacles.

3.3. The ﬂow past a cylinder

The ﬂow past a stationary circular cylinder is a typical problem and has been widely investigated experimentally

[32,33] and numerically [8,18,23,26,34–36]. For Reynolds number below 47, the ﬂow structure remains symmetric with stationary recirculating vortices behind the cylinder.

-1 -0.5 0 0.5 1 x -1 -0.5 0 0.5 1 y Y X 1 1 0 U= 1

Fig. 4. The conﬁguration of the lid-driven cavity in the computational domain. 0 0.25 0.5 0.75 1 Y -0.2 0 0.2 0.4 0.6 0.8 1 U 40X40 80X80 160X160 Ghia et al.

Fig. 5. The steady velocity component U along the line (y = x) through the center of cavity with diﬀerent Eulerian grid sizes.

0 0.25 0.5 0.75 1 X -0.2 0 0.2 0.4 0.6 0.8 1 V 40X40 80X80 160X160 Ghia et al.

Fig. 6. The velocity component V along the central line (y =x) through the center of cavity with diﬀerent Eulerian grid sizes.

(8)

As the Reynolds number is elevated, the symmetry breaks down and the vortex starts to shed up and down alterna-tively. This shedding frequency and the intensity of the vor-tex also increase in tandem with the elevated level of the Reynolds number.

In this test, ﬂows with a broad range of Reynolds num-ber, Re = 20, 40, 80, 100 and 150 are examined. The geo-metric set up in the computational domain and the associate physical boundary conditions are shown in

Fig. 8. Based on the diameter D of the cylinder, the range of the computational domain chosen is (13.4D 6 x 6 16.5D,8.35D 6 y 6 8.35D). A non-uniform grid (250 · 160) is adopted to discretize the computational domain, within which a uniform grid 60· 60 is employed in the region D 6 x, y 6 D. Here, the maximum CFL number employed is 0.37. In the simulation, the boundary of the stationary cylinder is modeled by the present immersed boundary methods, where the Lagrangian markers are assigned the no-slip boundary conditions.

There are three diﬀerent quantities that are often used to make a comparison of the performance of numerical meth-ods; namely, the drag and lift coeﬃcients, and the Strouhal

number. The deﬁnitions and computing of these quantities are brieﬂy described below.

The drag coeﬃcient is deﬁned as

CD¼

FD

U2₁D=2 ð18Þ

where FDis the drag force. As in[18], the drag force can be

obtained easily from

FD¼ Z X f1ðxÞ dx ¼ X x f1ðxÞh2 ð19Þ

where f1(x) is the x component of the Eulerian force f(x).

Similarly, the lift coeﬃcient is deﬁned as

CL ¼

FL

U2₁D=2 ð20Þ

where FLis the lift force and it can be obtained by the

Eule-rian force as FL¼ Z X f2ðxÞ dx ¼ X x f2ðxÞh2 ð21Þ

where f2(x) is the y component of the Eulerian force f(x).

The other interesting quantity called the Strouhal num-ber (St) is the dimensionless frequency of vortex shedding. When the ﬂow ﬁeld becomes unstable, the originally sta-tionary vortices behind the cylinder start moving down-stream and shedding alternatively with frequency fq. The

Strouhal number is deﬁned as St = fqD/U1.

The influences of Ds/h on the predicted lift and drag coefficients and Strouhal number are shown in Table 1, where the Reynolds number is 100. It can be seen that the differences of CD and St resulting from the variations

of Ds/h ratios are within 3%. It is also observed that Ds/ h = 0.94 and 0.78 produce the similar results. Therefore, in subsequent computations for ﬂows past a cylinder, the Ds/h is chosen as 0.94.

Table 2 shows the comparison of drag coefficient with previous numerical predictions [8,18,23,26,34,35] and experimental measurements[33]at different Reynolds num-bers. In[34,35], the flows were assumed to be steady and symmetric, therefore only the results of Re = 20 and 40 are shown here. The present drag predictions at lower Rey-nolds numbers are slightly higher than others, but at higher Reynolds numbers they are comparable with other results.

Table 3 shows the comparison of lift coeﬃcient for

Re = 100. Here, CL is deﬁned as the maximum amplitude

-1 -0.5 0 0.5 1 X -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 Y 1 m/s

Fig. 7. The velocity vector plot of lid-driven cavity with multi-obstacles.

u=1 v=0 ∂p ∂x ∂p ∂u ∂y ∂y =0 ∂v ∂x ∂u ∂x ∂p ∂x ∂p ∂x =0 =0 =0 =0 v=0 ∂p ∂u ∂y=0 ∂y=0 v=0 ∂p ∂x =0

Fig. 8. The boundary condition and computational domain of ﬂow over cylinder.

Table 1

Inﬂuences of the marker spacing to grid width on the lift and drag coeﬃcients and Strouhal number for the case of Re = 100

Ds/h 0.94 0.78 0.67 0.58 0.52

CL 0.34 0.34 0.30 0.31 0.32

CD 1.40 1.41 1.38 1.42 1.42

(9)

of the lift coeﬃcient. Table 4 shows the comparison of Strouhal number with the previous numerical and experi-mental results for diﬀerent Reynolds numbers. One can see that we obtained the close results with those methods as mentioned in the introduction.

For the flow of Re = 40, the drag coefficient remains constant and the lift is zero. The vorticity contour lines of the flow is shown inFig. 9. It has confirmed the exper-imental observation that at this Reynolds number, there are two vortices attached behind the cylinder and the lift force is constantly zero. However, if the Reynolds number becomes larger, the symmetry breaks down and the two vortices shed alternatively.Fig. 10(a)–(b) shows the evolu-tion of the drag and lift coefficients, and the vorticity con-tour lines for Re = 100.

4. The ﬂow past an in-line oscillating cylinder

In order to explore the capability of the present tech-nique in computing moving boundary flows, simulations are further applied to the in-line oscillating cylinder in uni-form flow at Reynolds number 100. The present flow lay-out and numerical details are exactly the same as the flow past a stationary cylinder in the previous section, except the cylinder is now oscillating parallel to the free stream at a frequency fc= 2fq, i.e. two times the vortex shedding

frequency (fq) of the ﬁxed cylinder ﬂow. The motion of

the cylinder is prescribed by setting the horizontal velocities on the Lagrangian markers to U(Xk) = 0.14D cos(2pfct),

where the amplitude of the oscillation is 0.14 of the cylinder diameter D. This ﬂow has been examined numerically by Hurlbut et al. [37]using the non-inertia coordinate trans-formation technique, and the results are used to examine the accuracy of the present predictions using the immersed boundary method.

According to Tanida et al. [38] and Griﬃn, Ramberg

[39]and Hurlbut et al.[37], the in-line oscillation of the cyl-inder at a range of frequencies near twice the Strouhal shedding frequency for the stationary cylinder causes the synchronization, i.e. the phase-locking of the vortex shed-ding with the cylinder motion. Both the drag and lift coef-ﬁcients increase. Table 5 shows the comparisons of the present predictions with the numerical results by Hurlbut et al. [37]. It can be seen that the present results are com-patible with their results. Further examination of the vor-tex patterns is referred to the instantaneous vorticity contours over two oscillating periods of the cylinder as shown inFig. 11.

5. Conclusion

In this paper, a new immersed boundary technique is proposed for the simulation of two-dimensional viscous incompressible flow interacting with complex solid bound-ary. A mixture of Eulerian and Lagrangian variables is adopted, where the solid boundary is represented by dis-crete Lagrangian markers embedding in and exerting forces to the Eulerian fluid domain. The interactions between the Lagrangian markers and the fluid variables are linked by a simple discretized delta function. The numerical integra-tion is based on a second-order fracintegra-tional step method under the staggered grid spatial framework. Based on the direct momentum forcing on the Eulerian grids, a new

Table 4

The comparison of Strouhal number for diﬀerent Reynolds numbers Re Present Lai and

Peskin[18]

Silva et al.[23] Ye et al.[8] Williamson (exp)[32] 80 0.153 – 0.15 0.15 0.15 100 0.168 0.165 0.16 – 0.166 150 0.187 0.184 0.18 – 0.183 -0.5 0 0.5. 1 1.5 2 2.5 3 X -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 Y

Fig. 9. The instantaneous vorticity contours near the cylinder; dotted and solid lines denote negative and positive contours.

Table 2

The comparison of drag coeﬃcients for diﬀerent Reynolds numbers

Re Present Lai and Peskin[18] Kim et al.[26] Silva et al.[23] Ye et al.[8] Fornberg[34] Ingham and Tang[35] Tritton (exp)[33]

20 2.20 – – 2.04 2.03 2.000 1.995 2.22 40 1.63 – 1.51 1.54 1.52 1.498 – 1.48 80 1.43 – – 1.40 1.37 – – 1.29 100 1.40 1.44 1.33 1.39 – – – – 150 1.39 1.44 – 1.37 – – – – Table 3

The comparison of maximum lift coeﬃcient at Re = 100

Re = 100 Present Lai and Peskin[18] Kim et al.[26]

(10)

force formulation on the Lagrangian marker is proposed, where the boundary forces are first computed on the Lagrangian markers and then distributed to the Eulerian grids using discrete delta function to ensure the satisfaction of the no-slip boundary condition on the immersed bound-ary in the intermediate step. This forcing procedure involves solving a banded linear system of equations whose unknowns consist of the boundary force on the Lagrangian markers; thus, the order of the unknowns is one-dimen-sional lower than the fluid variables. Four different test problems are simulated using the present technique

(rotat-ing r(rotat-ing flow, lid-driven cavity and flows over a stationary cylinder and an in-line oscillating cylinder). Numerical experiments show that the stability limit is not altered by the proposed force formulation, though the second-order accuracy of the adopted numerical scheme is degraded to 1.5 order. Influences of the marker spacing on the solutions are also examined numerically, and the optimal value is found to scale with the grid width. The numerical evidences show the accuracy and the capability of the proposed method for solving complex geometry flow problems both with stationary and moving boundary.

Acknowledgements

This research work was supported by the National Sci-ence Council of Taiwan under grant NSC-93-2212-E-007-036 and the computational facilities were provided by the National Centre for High-Performance Computing of Taiwan which the authors gratefully acknowledge.

22 24 26 28 time -0.5 -0.25 0 0.25 0.5 0.75 1 1.25 1.5

Drag and Lift coefficient

Drag coef. Lift coef. -0.5 0 0.5 1 1.5 2 2.5 3 X -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 Y (a) (b)

Fig. 10. (a) The time evolution of drag and lift coeﬃcients of Re = 100 and (b) the instantaneous vorticity contours near the cylinder; dotted and solid lines denote negative and positive contours.

Table 5

The comparisons of lift and drag coeﬃcients of in-line oscillating cylinder at Reynolds number 100

Present Hurlbut et al.[37]

fc/fq 0 2 0 2

CD 1.4 1.70 1.41 1.68

(11)

References

[1] Osher S, Shu CW. High-order essentially non-oscillatory schemes for Hamilton-Jacobi equations. SIAM J Numer Anal 1991;28:907. [2] Jiang GS, Shu CW. Eﬃcient implementation of weighted ENO

schemes. J Comp Phys 1996;126:202.

[3] Lin CH, Lin CA. Simple high-order bounded convection scheme to model discontinuities. AIAA J 1997;35:563.

[4] Brackbill JU, Saltzman JS. Adaptive zoning for singular problems in two dimensions. J Comp Phys 1982;46:342.

[5] Dvinsky AS. Adaptive grid generation from Harmonic maps on Riemannian-Manifolds. J Comp Phys 1991;95:450.

[6] Li R, Tang T, Zhang PW. Moving mesh methods in multiple dimensions based on harmonic maps. J Comp Phys 2001;170:562. [7] Di Y, Li R, Tang T, Zhang PW. Moving mesh ﬁnite element methods

for the incompressible Navier–Stokes equations. SIAM J Sci Comput 2005;26:1036.

[8] Ye T, Mittal R, Udaykumar HS, Shyy W. An accurate Cartesian grid method for viscous incompressible ﬂows with complex immersed boundaries. J Comp Phys 1999;156:209.

[9] LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coeﬃcients and singular sources. SIAM J Numer Anal 1994;31:1019.

[10] Calhoun D. A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions. J Comp Phys 2002;176:231.

[11] Li Z, Wang C. A fast ﬁnite diﬀerence method for solving Navier– Stokes equations of irregular domains. Commun Math Sci 2003;1:180.

[12] Li Z, Lai M-C. The immersed interface method for the Navier–Stokes equations with singular forces. J Comp Phys 2001;171:822.

[13] Le DV, Khoo BC, Peraire J. An immersed interface method for the incompressible Navier–Stokes equations in irregular domains. Pro-ceedings of the third MIT conference on computational ﬂuid and solid mechanics, 710. Elsevier Science; 2005.

[14] Xu S, Wang ZJ. An immersed interface method for simulating the interaction of a ﬂuid with moving boundaries. preprint; 2005. [15] Li Z. An overview of the immersed interface method and its

applications. Taiwanese J Math 2003;7:1.

[16] Peskin CS. Flow patterns around heart valves: a numerical method. J Comp Phys 1972;10:252.

[17] Peskin CS. The immersed boundary method. Acta Numer 2002;1. [18] Lai M-C, Peskin CS. An immersed boundary method with formal

second order accuracy and reduced numerical viscosity. J Comp Phys 2000;160:705.

[19] Goldstein D, Handler R, Sirovich L. Modeling a no-slip ﬂow with an external force ﬁeld. J Comp Phys 1993;105:354.

[20] Goldstein D, Hadler R, Sirovich L. Direct numerical simulation of turbulent ﬂow over a modeled riblet covered surface. J Fluid Mech 1995;302:333.

[21] Saiki EM, Biringen S. Numerical simulation of a cylinder in uniform ﬂow: application of a virtual boundary method. J Comp Phys 1996;123:450.

[22] Beyer RP, LeVeque RL. Analysis of a one-dimensional model for the immersed boundary method. SIAM J Numer Anal 1992;29:332. [23] Lima E Silva ALF, Silveira-Neto A, Damasceno JJR. Numerical

simulation of two-dimensional ﬂows over a circular cylinder using the immersed boundary method. J Comp Phys 2003;189:351.

[24] Mohd-Yusof J. Combined immersed boundary/B-Spline method for simulations of ﬂows in complex geometries in complex geome-tries CTR annual research briefs. NASA Ames/Stanford University; 1997. -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (b) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (c) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (d) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (e) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (f) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (g) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (h) -0.5 0 0.5 1 1.5 2 2.5 3 X -0.5 0 0.5 Y (a)

Fig. 11. The instantaneous vorticity contours near the oscillating cylinder—Re = 100; dotted and solid lines denote negative and positive contours: (a) t = T/4, (b) t = T/2, (c) t = 3T/4, (d) t = T, (e) t = 5T/4, (f) t = 3T/2, (g) t = 7T/4 and (h) t = 2T (T is the oscillation period of the cylinder).

(12)

[25] Fadlun EA, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary methods for three dimensional complex ﬂow simulations. J Comp Phys 2000;161:30.

[26] Kim J, Kim D, Choi H. An immersed-boundary ﬁnite-volume method for simulations of ﬂow in complex geometries. J Comp Phys 2001;171:132.

[27] Tseng Y-H, Ferziger JH. A ghost-cell immersed boundary method for ﬂow in complex geometry. J Comp Phys 2003;192:593.

[28] Harlow FH, Welsh JE. Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with a free surface. Phys Fluids 1965;8:2181–9.

[29] Choi H, Moin P. Eﬀects of the computational time step on numerical solutions of turbulent ﬂow. J Comp Phys 1993;113:1.

[30] Van den Vorst HA, Sonneveld P. CGSTAB, a more smoothly converging variant of CGS, Technical Report 90-50, Delft University of Technology; 1990.

[31] Ghia U, Ghia KN, Shin CT. High-Re solutions for incompressible ﬂow using the Navier–Stokes Equations and a multi-grid method. J Comp Phys 1982;48:387.

[32] Williamson CHK. Vortex dynamics in the cylinder wake. Ann Ver Fluid Mech 1996;28:477.

[33] Triton DJ. Experiments on the ﬂow past a circular cylinder at low Reynolds number. J Fluid Mech 1959;6:547.

[34] Fornberg B. A numerical study of steady viscous ﬂow past a circular cylinder. J Fluid Mech 1980;98:819.

[35] Ingham DB, Tang T. A numerical investigation into the steady ﬂow past a rotating circular cylinder at low and intermediate Reynolds numbers. J Comp Phys 1990;87:91.

[36] Tang T, Ingham DB. On steady ﬂow past a rotating circular cylinder at Reynolds numbers 60 and 100. Comput Fluids 1991;19:217. [37] Hurlbut SE, Spaulding ML, White FM. Numerical solution for

laminar two dimensional ﬂow about a cylinder oscillating in a uniform stream. Trans ASME, J Fluids Eng 1982;104:214.

[38] Tanida Y, Okajima A, Watanabe Y. Stability of a circular cylinder oscillating in uniform ﬂow or in a wake. J Fluid Mech 1973;61:769. [39] Griﬃn OM, Ramberg SE. Vortex shedding from a cylinder

vibrat-ing in line with an incident uniform ﬂow. J Fluid Mech 1976;75: 257.