- Email: [email protected]

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Simulating ﬂows with moving rigid boundary using immersed-boundary method Chuan-Chieh Liao a, Yu-Wei Chang a, Chao-An Lin a,*, J.M. McDonough b a b

Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 30013, Taiwan Departments of Mechanical Engineering and Mathematics, University of Kentucky, Lexington, KY 40506, USA

a r t i c l e

i n f o

Article history: Received 10 March 2009 Received in revised form 10 July 2009 Accepted 23 July 2009 Available online 6 August 2009

a b s t r a c t The present study is to apply the immersed-boundary method to simulate 2- and 3-D viscous incompressible ﬂows interacting with moving solid boundaries. Previous studies indicated that for stationary-boundary problems, different treatments inside the solid body did not affect the external ﬂow. However, the relationship between internal treatment of the solid body and external ﬂow for movingboundary problems was not studied extensively and is investigated here. This is achieved via directmomentum forcing on a Cartesian grid by combining ‘‘solid-body forcing” at solid nodes and interpolation on neighboring ﬂuid nodes. The inﬂuence of the solid body forcing within the solid nodes is ﬁrst examined by computing ﬂow induced by an oscillating cylinder in a stationary square domain, where signiﬁcantly lower amplitude oscillations in computed lift and drag coefﬁcients are obtained compared with those without solid-body-forcing strategy. Grid-function convergence tests also indicate second-order accuracy of this implementation with respect to the L1 norm in time and the L2 norm in space. Further test problems are simulated to examine the validity of the present technique: 2-D ﬂows over an asymmetrically-placed cylinder in a channel, in-line oscillating cylinder in a ﬂuid at rest, in-line oscillating cylinder in a free stream, two cylinders moving with respect to one another, and 3-D simulation of a sphere settling under gravity in a static ﬂuid. All computed results are in generally good agreement with various experimental measurements and with previous numerical simulations. This indicates the capability of the present simple implementation in solving complex-geometry ﬂow problems and the importance of solid body forcing in computing ﬂows with moving solid objects. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Body-ﬁtted or unstructured-grid methods have been commonly adopted to simulate ﬂow in domains with complex rigid boundaries, such as in ﬂuid–solid interaction problems. However, the computational cost and memory requirements of these methods are generally high. Also, modeling complex time-evolving moving-boundary ﬂows requiring transient re-meshing strategies further increases the computational overhead and algorithmic complexity of these approaches. Thus, one would seek to adopt the Cartesian grid based non-boundary conforming methods to address complex ﬂuid–structure interaction problems. Since boundary points now do not conform to the underlying Cartesian grid, the main challenges are to model possibly complicated immersed boundaries accurately and at the same time maintain efﬁciency in solving the associated governing equations discretized on the Cartesian grid.

* Corresponding author. Address: Department of Power Mechanical Engineering, National Tsing Hua University, 101, Section 2, Kuang Fu Road, Hsinchu 30013, Taiwan. Tel.: +886 3 5742602; fax: +886 3 5722840. E-mail address: [email protected] (C.-A. Lin). 0045-7930/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compﬂuid.2009.07.011

There are many approaches available. For example, sharp interface method (Ye et al. [1]), tracks the boundary by identifying the control volumes (cut cells) in the Cartesian grid that are cut by the immersed boundary. Consequently, when simulations of ﬂow with moving boundaries are attempted (see, e.g., Udaykumar et al. [2,3] and Marella et al. [4]) the irregular shapes of cut cells result in a need for complicated interpolation procedures to approximate ﬂuxes, and this unavoidably affects computational efﬁciency of Cartesian-grid solvers. Another approach, also within the Cartesian-grid framework, the immersed interface method (IIM), originally developed by LeVeque and Li [5], employed jump conditions to produce sharp solutions near boundaries. Implementation of the IIM has been carried out by Calhoun [6] and Li and Wang [7] in the context of vorticity-stream function formulations. Based on work of Li and Lai [8] for the Navier–Stokes equations with singular forcing, the IIM has also been explored by Le et al. [9], Russell and Wang [10] and Xu and Wang [11] using the Navier–Stokes (N.–S.) equations in primitive-variable form. However, the implementation of the jump conditions may not be trivial. An alternative that not only retains the advantage of structured Cartesian grids but also provides an ability to handle complex geometries is the so-called immersed-boundary method

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

introduced by Peskin [12,13], where complex geometries within a Cartesian grid can be replaced by generating an external force ﬁeld which acts on the ﬂuid in the same manner as would a solid boundary. In general, immersed-boundary methods are categorized as either feedback-forcing or direct-forcing approaches [14]. In addition to the works of Peskin, feedback forcing has been studied by numerous other researchers, including Goldstein et al. [15,16], Saiki and Biringen [17], Lai and Peskin [18] and Silva et al. [19]. It is applied to the equations of ﬂuid motion discretized on a Cartesian grid, and it spreads forcing effects to neighboring grid points via a distribution function to account for the presence of any immersed boundaries. Feedback forcing is a particularly attractive immersed-boundary formulation for problems with elastic boundaries, and successful simulations of biological and multi-phase ﬂows have been accomplished within this framework. The main problem associated with feedback-forcing approaches is that smoothing of the forcing function prohibits sharp representation of the immersed boundary. This is often not acceptable at high Reynolds numbers and leads to smearing of interface information. To maintain the suitable resolution, Uzgoren et al. [20] adopted local adaptive grid reﬁnement techniques near the irregularly shaped boundaries. A detailed review on this topics applying to multiphase ﬂow can be referred to Uzgoren et al. [21]. However, the feedback-forcing approach’s another disadvantage is the restriction of the time step (CFL number < Oð102 Þ) [15,17]. On the other hand, the direct-forcing approach explicitly enforces boundary conditions at Eulerian nodes near the immersed boundary. Mohd-Yusof [22] proposed a direct-forcing method that introduces a body force such that the desired velocity distribution is satisﬁed at the boundary. Fadlun et al. [23] further extended the approach of [22] to a ﬁnite-difference formulation on a staggered grid, where direct forcing was applied at the ﬁrst Eulerian grid points external to the immersed boundary. Stability of the timeintegration scheme was not altered, and good quantitative agreement with experimental measurements was achieved. Other direct forcing formulations include Kim et al. [24], Tseng and Ferziger [25], Balaras [26], Gilmanov and Sotiropoulos [27], Kim and Choi [28], Yang and Balaras [29], Su et al. [30], Zhang and Zheng [31], Choi et al. [32] and Ghias et al. [33]. Moving objects within a ﬁxed grid are usually a source of error for non-boundary-conforming strategies on a Cartesian grid. In a feedback-forcing approach, such deﬁciencies are usually minimal due to the smooth transition from solid to ﬂuid. In the direct-forcing strategy, however, difﬁculties arise because of the role change of the ﬂuid–solid interface Eulerian grid points from time step to time step. In particular, a point in the solid can become a point in the ﬂuid (or, vice versa) at the next time step. This produces severe oscillations of the predicted lift and drag coefﬁcients. Miller and Peskin [34] recommends that lift and drag coefﬁcients were ﬁltered to remove the high frequency noise. Alternatively, several remedies have been proposed for treating this deﬁciency. In order to handle such ‘‘phase-change” points, Yang and Balaras [29] proposed a ﬁeld-extension procedure to treat points emerging from a moving solid into the ﬂuid. However, evaluations of lift and drag coefﬁcients may not be trivial. Kim and Choi [28] developed an immersed-boundary method in a non-inertial reference frame to simulate ﬂow around an arbitrarily-moving body in an inﬁnite domain. However, the non-inertial reference frame approach exhibited difﬁculties in simulating ﬂow interactions between two (or more) moving bodies. Another approach is to adopt a combination of Eulerian and Lagrangian variables, such that the solid boundary was represented by discrete Lagrangian markers embedded in, and exerting forces on, the Eulerian ﬂuid domain [30,31]. Based on direct-momentum forcing on the Eulerian grid, Su et al. [30] proposed a new implicit force formulation on the Lagrangian marker to ensure exact satis-

153

faction of the no-slip boundary condition on the immersed boundary. Interactions between the Lagrangian markers and the ﬂuid variables on the ﬁxed Eulerian grid are linked by a simple discrete delta function as in the spreading procedure of feedback forcing approaches. A similar strategy was applied in lattice Boltzmann simulations reported by Chen et al. [35] and Chang et al. [36]. Although stability of the time-integration scheme was not altered using this approach, smoothing of the forcing function may potentially prohibit sharp representation of the immersed boundary. On the other hand, the inﬂuence of the internal treatment of the solid body was investigated in [23,37], and it was indicated that treatments inside the body did not affect the external ﬂow when the body is stationary. Yang and Balaras [29] further indicated that the solid points do not inﬂuence the rest of the computation. However, the relationship between internal treatment of the solid body and external ﬂow for moving-boundary problems is not studied extensively. In the present study, an immersed-boundary technique is used for simulation of 2- and 3-D viscous incompressible ﬂow interacting with moving solid boundaries. The aim is to investigate the inﬂuence of solid body forcing within the solid nodes and provide a procedure which is easy to implement during non-stationary ﬂuid and solid interaction. This is based on direct-momentum forcing on a Cartesian grid utilizing a combination of interpolation at ﬂuid nodes adjacent to the solid body and solid-body forcing at nodes which had been part of the solid region in the current time step. The inﬂuence of the solid body forcing within the solid nodes, grid-function convergence tests and ﬂow with high Reynolds number are examined by computing ﬂow induced by an oscillating cylinder in a stationary square domain. Further test problems are simulated to examine the validity of the present technique: 2D ﬂows over an asymmetrically-placed cylinder in a channel, inline oscillating cylinder in a ﬂuid at rest, in-line oscillating cylinder in a free stream, two cylinders moving with respect to one another, and 3-D simulation of a sphere settling under gravity in a static ﬂuid. These cases are used to examine capability of the present method in solving complex-geometry ﬂow problems and the importance of solid body forcing in computing ﬂows with moving solid objects. 2. Methodology of the immersed-boundary technique In this section, the general mathematical formulation of the problems under consideration is ﬁrst introduced. This is followed by the numerical methods to be employed, and then details of constructing the forcing function needed to represent effects of the immersed boundary. This section is concluded with a fairly detailed description of the complete solution algorithm. 2.1. Mathematical formulation Consider the problem of a viscous incompressible ﬂuid in a 2-D rectangular domain X ¼ ½0; L ½0; H containing an immersed massless boundary in the form of a simple closed curve C, as shown in Fig. 1. The immersed boundary is tracked in parametric form, XðsÞ; 0 6 s 6 Lb , with s being an arclength parameter associated with the (immersed) boundary C, and Lb is the length of C. It is clear from this ﬁgure that few, if any, of the points comprising C coincide with grid points. As mentioned above, the inﬂuence of the immersed boundary on the ﬂuid is represented by forces exerted on the ﬂuid and causing it to move by the prescribed velocity (via the no-slip condition) at the immersed boundary. Thus, the governing equations of this ﬂuid–structure interaction system are the usual incompressible N.–S. equations,

154

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

H

y x L Fig. 1. Flow domain X with an immersed boundary C.

@u 1 þ r ðuuÞ ¼ rp þ Du þ f ; @t Re r u ¼ 0:

ð1Þ ð2Þ

Here, x ¼ ðx; yÞ; uðx; tÞ is the ﬂuid velocity having components u and v ; pðx; tÞ is the ﬂuid pressure, and Re is the dimensionless parameter, Reynolds number; r and D are the usual Cartesian-coordinate gradient and Laplacian. Note that the discrete-time momentum forcing f ðx; tÞ is applied to satisfy the no-slip condition on immersed boundaries, as in [22], and momentum forcing is applied only at nodes adjacent to the immersed boundary in an Eulerian grid system. We assume throughout that sufﬁcient initial and boundary data are provided to constitute well-posed mathematical problems for these equations. Speciﬁc details of these required conditions will differ from problem to problem in the sequel. 2.2. Numerical scheme The numerical procedure used herein consists of a ﬁnite-volume method discretized in Cartesian coordinates with a staggered-grid arrangement of dependent variables. This is based on integration of the transport equations over arbitrary control volumes, leading to conservation of mass and balance of momentum, and any scalar ﬂow property, over each such volume. Spatial derivatives are approximated using second-order centered differencing, and a fractional-step (projection) method implemented with a combination of Adams–Bashforth and Crank–Nicolson methods for advective and diffusive terms, respectively, is used for temporal discretization. Several variants of fractional-step methods are available, such as Kim and Moin [38], Bell et al. [39] and Choi and Moin [40]. Here, the time-advancement formulation of Choi and Moin [40] is employed for the momentum equations due to its storage economy and straightforward implementation of boundary conditions. The discrete form of this algorithm is embodied in the following sequence of equations:

e un 3 1 u ¼ rh ðuuÞn þ rh ðuuÞn1 2 2 Dt 1 e Þ þ f nþ1 ; rh pn þ Dh ðun þ u 2Re e u u ¼ rh pn ; Dt r h u Dh / ¼ ; Dt unþ1 u ¼ rh /; Dt 1 e; rh u pnþ1 ¼ / 2Re

ð3Þ

ð4Þ ð5Þ ð6Þ ð7Þ

e and u might be viewed as intermediate velocities occurwhere u ring between the time steps n and n þ 1 although they are, in fact, non-divergence-free approximations to unþ1 . The discrete operators rh and Dh represent the usual centered-difference approximations for the gradient and Laplace operators, respectively. To evaluate nþ1 must be determined in advance to satisfy Eq. (3), the forcing f the no-slip boundary condition on any immersed boundary at the advanced time level. This will be treated in the next subsection. The most time-consuming part of executing the above scheme is solving three (in 2-D) Helmholtz-type equations (two for velocity components in (3) and one for pressure in (5)) which can be done by the efﬁcient Bi-CGSTAB method of Van den Vorst and Sonneveld [41]. In particular, at least Oð1Þ iterations per time step for the momentum equations are typically required to lower iteration errors to within the local truncation error of the overall method, in addition to the Oð102 Þ iterations required by pressure Poisson equation solves. This requires signiﬁcantly more arithmetic than nþ1 , as will be apparent from what is needed for construction of f follows. 2.3. Forcing strategies The force applied to the computational node due to presence of a solid body must be determined at the advanced time level before the solution procedure presented in Eqs. (3)–(6) can be started. This force (per unit mass) is simply a Newtonian acceleration of any node near the solid body. Thus,

f

nþ1

¼

b uF u ; Dt

ð8Þ

b is an estimate of the time-level n þ 1 velocity ﬁeld in a where u neighborhood of the solid boundary in the absence of forcing, and uF is the velocity (at the same location) when forcing is taken into b is estimated with a simple, explicit Adams–Bashaccount. Here, u forth scheme where Eq. (1) is provisionally discretized explicitly in time to derive the momentum forcing value, as described in Kim et al. [24]:

( n 3 1 b ¼ u Dt rh ðuuÞ Dh u u 2 Re ) n1 1 1 n rh ðuuÞ Dh u þ rh p : 2 Re n

ð9Þ

The issue of computing the forcing function f of Eq. (8), and hence uF , when the interface does not coincide with grid nodes is discussed here. A combination of two strategies is adopted for calculations to be reported herein: viz., interpolation [23] (at forcing points in the ﬂuid domain) and interpolation plus a solid-bodyforcing procedure, the latter of which leads to improved accuracy (in particular, signiﬁcant reduction of non-physical temporal oscillations of dependent variables) when treating moving boundaries. Thus, two issues should be addressed here. First is the forcing location in the ﬂuid region, and the second is speciﬁc details of solid-body forcing. A schematic of forcing points in the ﬂuid region used to calculate uF is shown in Fig. 2. This is carried out via linear interpolation in the present study, but more accurate interpolation methods could be employed if needed. From the ﬁgure it can be seen that forcing nodes are Cartesian grid points outside the solid body, but immediately adjacent to it. The no-slip condition is imposed at points of the boundary, C, of the solid body at points xB lying on the same grid line as the forcing point xF . Then the velocity uF at this point is obtained by (linear) interpolation between points xB and xA , with the latter being the next point outward from xF on the grid line containing points xB and xF . Thus, for the u velocity component at the point xF in the inset to Fig. 2,

155

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

A

Fluid Domain forcin ng nodes A

So olid Dom main B

F

F2

F1

A

ambigu uous node e

B F ambiguo uous node

F3

B

Flu uid Dom main

Solid Domain

B

F

A

Solid Domain ambiguous node

B

F

Fluid Domain

A

Fig. 2. Location of nodes A; B and forcing point F; treatment of ambiguous interpolation points in 2-D and 3-D.

bA uF ¼ u

bA u bB u ðyA yF Þ; yA yB

ð10Þ

with analogous computations used, scalarly, for other velocity components (in both 2-D and 3-D, the latter of which is indicated in the lower part of the ﬁgure). nþ1 This leads to direct calculation of f via substitution of (9) and (10) into Eq. (8):

f

nþ1

¼

n uF un 3 1 þ rh ðuuÞ Dh u 2 Re Dt n1 1 1 rh ðuuÞ Dh u þ rpn : 2 Re

a suitable way to deal with the ambiguity. Otherwise, interpolated forcing points F 1 and F 2 can be used as the interpolated ﬂuid points to obtain the velocity of the forcing node F 3 . Also, the solid-body-forcing strategy is imposed within the solid region through the solid-body force calculated as

f

nþ1

¼

b U mov e u vðxnþ1 Þ; Dt

ð12Þ

where U mov e represents the velocity of the moving object. It is clear that within the solid region uF ¼ U mov e in Eqs. (8) and (11). The parameter v is merely a characteristic function for points of the solid, deﬁned in the usual way as

ð11Þ

In convex solid boundary, the preceding approach works well if the forcing point is associated with only one solid node, as is true for forcing nodes shown in the upper-left part of Fig. 2. On the other hand, if the forcing point is associated with two (or more in 3-D) solid nodes, then ambiguity occurs as indicated in the ﬁgure. Here, two or three possible interpolations can be performed; an arbitrary selection was used in [23]. To remove this ambiguity we propose utilizing simple linear interpolation between the diagonal node xA and the node xB on the boundary, as shown in Fig. 2. A similar, more elaborate interpolation procedure was employed in [29]. In concave solid boundary however, as shown in the upperright part of Fig. 2, the present interpolation procedure provides

(

vðxnþ1 Þ ¼

1; xnþ1 2 solid region; 0; otherwise:

Note that the solid-body-forcing procedure must be carried out at every time step for moving-boundary problems. To conclude this treatment of forcing required to account for immersed boundaries, we note that the various ﬁgures presented here are generic. Application of the technique in the context of staggered grids, as employed in the current study, is performed scalarly for each individual set of points associated with discrete locations of u; v (and w in 3-D). Clearly, there are other fairly natural possible approaches which we have not yet explored.

156

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

2.4. Complete solution procedure Fluid

The full numerical procedure for each time step of the proposed method is summarized in the following pseudo-language algorithm.

D

Algorithm 1 (Fractional-step/solid-body forcing immersed-boundary method). Suppose n time steps have been completed. To calculate the solution at time level n þ 1 carry out the following steps.

1. Determine the immersed-boundary location at step n þ 1. (a) If the trajectory of the moving object is prescribed, the location of the object can be easily determined. (b) If the object moves in accordance with the net imposed force on the solid object (for example: balance of gravitational force ðf g Þ, drag force ðf D Þ and P n n n n n f ¼ f g þ f D þ f B ), then buoyancy ðf B Þ, then f net ¼ nþ1 Þ can be comacceleration ðanobject Þ and velocity ðuobject n puted, respectively, as anobject ¼ f net =mobject and nþ1 n n uobject ¼ uobject þ aobject Dt, where mobject and Dt are the mass of the object and time step, respectively. Finally the location of the object can be simply computed as Dt. xnþ1 ¼ xn þ unþ1 object 2. 3. 4. 5. 6.

7. 8. 9. 10.

Identify ﬂuid and solid points on Eulerian grid. Establish locations of forcing points. b. Use Eq. (9) to determine u nþ1 from Eq. (11) with uF Obtain ﬂuid region forcing term f calculated via interpolation using Eq. (10). nþ1 (solid-body-forcing Obtain solid region forcing term f procedure) from Eq. (11) with uF ¼ U mov e , when dealing with moving-boundary problems. e implicitly by Compute non-divergence-free velocity ﬁeld u nþ1 . solving Eq. (3) with force f Obtain intermediate velocity u from Eq. (4). Compute pressure at time level n þ 1; pnþ1 , by solving pressure Poisson equation, Eq. (5). Compute divergence-free time level n þ 1 ﬂuid velocity unþ1 via projection, Eq. (6).

Two remarks regarding this algorithm are needed to aid in its implementation. First, in the case of stationary, non-deforming solid bodies, steps (1) and (2) need be performed only at the beginb is not divergence free, ning of the simulation. Second, since u step (3) is required at only sufﬁciently many grid points of the ﬂuid near the immersed boundary to permit the interpolations carried out in step (4) and of the solid region in step (5). 3. Numerical results In this section we begin by demonstrating accuracy of the overall algorithm on a 2-D, but nontrivial, problem. We then proceed to solve several different problems for which laboratory experiments have been conducted, or numerical simulations have been performed by other researchers, or both. These results, for both 2-D and 3-D problems, provide a quite rigorous test of the algorithm, and at the same time lead to qualitative and (at least) semi-quantitative assessments of solid-body-forcing accuracy. All the simulations are performed with CFL numbers within the range 0.2–0.55. 3.1. Accuracy test for moving-boundary problem Here, 2-D ﬂow is induced by an oscillating circular cylinder in a stationary square domain with side L ¼ 1 (in arbitrary units) such that 0:5 6 x; y 6 0:5 shown in Fig. 3. The diameter D of the cylin-

L

Solid

L Fig. 3. Conﬁguration of moving-cylinder problem used in accuracy test.

der is 0.5, and the motion is prescribed by setting the horizontal location of the center to XðtÞ ¼ 0:1 cosð2pftÞ, where the oscillation frequency f of the cylinder is equal to 1 Hz. No-slip boundary conditions are imposed on the walls of the box containing the cylinder, _ the Reynolds numbers for as indicated in Fig. 3. Based on D and X, these calculations are 100 and 800. Three kinds of direct-forcing strategies are considered in these simulations: extrapolation (as used in [24]), interpolation (forcing points in the ﬂuid domain) and interpolation plus the proposed solid-body-forcing procedure. The drag force F D and lift force F L generated by the cylinder can be obtained by integrating the local pressure and stress distributions along the cylinder wall. However, this method may not be straightforward due to the interpolation procedure involved. Another simple alternative is by applying the volume integral of the Navier–Stokes equation, and the momentum deﬁcits due to the presence of the cylinder are the drag and lift.

FD ¼

@ qu @ quu @ quv @p @ @u @ @u l l þ þ dv ol; þ @x @x @x @x @y @y @y v ol @t

Z

ð13Þ FL ¼

@ qv @ quv @ qvv @p @ @v @ @v l l þ þ þ dv ol: @y @x @x @y @y @t @x @y v ol

Z

ð14Þ In the present study, another alternative to compute the drag and lift incurred by the presence of the cylinder is used here [18,30],

FD ¼

Z

f1 ðxÞ dx ﬃ

X

X

2

f1 ðxÞ h ;

ð15Þ

x nþ1

where f1 ðxÞ is the x component of the force f in Eq. (3). Similarly, the lift force F L is calculated from

FL ¼

Z X

f2 ðxÞ dx ﬃ

X

2

f2 ðxÞ h ;

ð16Þ

x

with f2 ðxÞ being the y component of the force. To validate the accuracy of the present force calculation, a comparison is made with the force obtained by the volume integral of the Navier–Stokes equation. The ﬂow induced by an oscillating cylinder in a stationary square domain at Re = 100 is adopted here, and the forcing strategy is interpolation plus the proposed solidbody-forcing procedure. Table 1 shows the inﬂuences of grid densities and force computing methods on the predicted mean drag force. It is clear that the present method agrees perfectly with the volume integral method. Therefore, instead of resorting to

157

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167 Table 1 Mean drag force F D variations using different grid for the ﬂow induced by an oscillating cylinder in a stationary square domain at Re = 100. GRID

Mean drag Eq. (15) (present study)

Mean drag Eq. (13) (volume integral)

jF D1 F D2 j jF D1 ðref Þj

100 100 200 200 400 400 (ref)

0.00247627814571708 0.00253377798404919 0.00253345120147349

0.00247627773960265 0.00253377795814951 0.00253345119950619

0.00001603 0.00000102 0.00000008

the complex integration of the pressure and shear stress distributions along the surface, the present force calculation is a simple and yet accurate way of computing moving object ﬂows. If the detailed pressure and stress distributions along the surface are needed for other purposes, thus interpolation procedure similar to that adopted by Yang and Balaras [29] could be adopted. Fig. 4 shows evolution of drag force F D (Reynolds number being 100 and with grid density of 300 300) for the different forcing strategies after 10 periods, from which we see that differences among these are very signiﬁcant (as much as 40% discrepancy at some instances of time). It can be observed from Fig. 4, especially in the zoom-in ﬁgure on the right, that severe oscillations of predicted drag forces occur when using either interpolation alone or extrapolation, and it is interesting to note that extrapolation appears to perform somewhat better than interpolation for this particular problem. As indicated earlier, complications for moving-boundary problems arise from the fact that Eulerian forcing points near the interface change from time step to time step. Without modiﬁcation to the computational procedure, these points will be a source of non-physical oscillation in computed solutions. To quantify such oscillations for the current problem, single-sided (Fourier) amplitude spectra of the drag forces produced by each of the three methods were computed using MATLAB, with results shown in Fig. 5. It is clear that the dominant frequency of the drag forces is equal to the frequency of cylinder movement, i.e., 1 Hz, which is the same in each ﬁgure, implying that the low-frequency essentials of the drag force curves are (nearly) identical, but with a slightly lower amplitude for the fundamental frequency corresponding to interpolation. The remaining parts of the amplitude spectra for higher frequency reﬂect different oscillation details for each of the three approaches and, based on the simple form of cylinder translation used here, must be due to numerical errors. In parts (a) and (b) of Fig. 5, corresponding to extrapolation and interpolation results, one sees very distinct odd harmonics of the fundamental, with at least noticeable amplitude for frequencies as high as

19 Hz. Extrapolation results, displayed in Fig. 5a, are of somewhat lower amplitude at frequencies above the fundamental, but at the same time showing a somewhat slower decay of amplitude with frequency than in the interpolation case of part (b). Finally, interpolation combined with the solid-body-forcing strategy produces the lowest amplitude oscillations at essentially all frequencies above the fundamental. Indeed, except for a small peak at about 3 Hz (which is also seen in the other two power spectra, and is of signiﬁcant amplitude in the interpolation spectrum), there are no sharp peaks. It is quite clear that this approach produces a very effective ﬁltering procedure for these non-physical oscillations, suggesting its potential for use in problems involving ﬂows in, and around, complex geometries with moving boundaries. Grid-function convergence tests were conducted using four uniform grids of 200 200; 300 300; 400 400 and 800 800 cells. In the absence of analytical solutions for this moving-boundary problem, results from the ﬁnest grid were viewed as ‘‘exact” for purposes of grid-function convergence assessments. For all grids, the simulation was continued for a complete oscillation period to mitigate effects of the impulsive start. L1 - and L2 -error norms are used to measure the difference between velocity components computed on the coarser grids and the reference grid. For example, error norms for the u-component velocity are calculated as N 1 X L1 ¼ jui;c ui;f j; N jij>0

vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ u N u1 X and L2 ¼ t ðui;c ui;f Þ2 : N jij>0

3.5

Interpolation

3 Extrapolation Interp.+solid-body-forcing

drag force

drag force

2

0

2.5

2 -2 1.5 -4

3.6

3.8

4

4.2

4.4

4.6

time

4.8

5

5.2

5.4

ð17Þ

Here, N is the total number of cells in the discrete representation; ui;c and ui;f are numerical solutions from a coarser grid and the ﬁnest grid, respectively, and i is a generic multi index. Fig. 6 shows the relationship between velocity error norms and grid resolution. Dashed and solid lines indicate ﬁrst- and secondorder accuracy slopes, respectively. It is evident that the L2 error in velocity decreases asymptotically with Dx2 in a manner supporting the expected second-order spatial accuracy of the present method. This implies that the solid-body-forcing technique leads

Extrapolation Interpolation Interp.+solid-body-forcing

4

(%)

4.9

5

5.1

time

Fig. 4. Comparison of drag coefﬁcients for different forcing strategies.

5.2

158

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

25

(a) 20 15 10 5 0 0

2

4

6

8

10 12 Frequency (Hz)

14

16

18

20

25

(b) 20 15 10 5 0 0

2

4

6

8

10 12 Frequency (Hz)

14

16

18

20

25

(c) 20 15 10 5 0 0

2

4

6

8

10 12 Frequency (Hz)

14

16

18

20

Fig. 5. Amplitude spectra of drag forces for different forcing strategies: (a) extrapolation, (b) interpolation and (c) interpolation plus solid-body forcing.

Fig. 6. L1 - and L2 -error norms for velocity components.

159

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

to global second-order accuracy in space. The CFL number was held constant ð¼ 0:2Þ on all grids employed for these tests. This permits a straightforward transformation from space to time leading to an assessment of temporal accuracy shown in terms of the L1 norm. Again, second-order accuracy is indicated; we remark that with ﬁxed (and relatively large) Courant number as used here, this can also be (implicitly) deduced from the ﬁgure for L2 error with respect to Dx2 . The capability of the adopted approach is further investigated by computing ﬂow at higher Reynolds number, Re = 800 with 300 300 grid. The inﬂuence of the reduction of the viscosity can be clearly observed in Fig. 7, where thinner boundary layer is developed at the Re = 800 ﬂow in contrast to the thicker boundary layer prevailed at Re = 100 ﬂow. Also shown in Fig. 8, the evolution of the drag force for Re = 800 correlates well with the oscillation frequency of the cylinder with negligible oscillatory noise. 3.2. Flow over an asymmetrically-placed cylinder Schafer and Turek [42] reported a set of 2-D and 3-D numerical benchmark results using different boundary conforming numerical techniques for laminar ﬂows over a circular cylinder of diameter D that is asymmetrically placed inside a channel.

The distances from the center of the cylinder to the upper and lower walls are 2:1D and 2:0D, respectively, as shown in more detail in Fig. 9. A parabolic velocity proﬁle of maximum speed U max ¼ 1 is assigned at the channel inlet; the no-slip condition is applied along walls, and a convective boundary condition,

@u @u þc ¼ 0; @t @x is used at the outﬂow boundary. Here, c is space-averaged streamwise velocity at the exit. The Reynolds number, based on the average inlet velocity, is deﬁned as Re ¼ U av g D=m with U av g ¼ 2U max =3. A non-uniform grid of 300 200 cells is adopted to discretize the overall computational domain; within this grid is a reﬁned uniform grid of 60 60 cells employed in the region D 6 x; y 6 D, but with the additional grid lines extending into the remainder of the domain only in each separate direction, as shown in Fig. 9. The CFL number adopted is 0.3. At Re ¼ 100, used in the present calculations, the ﬂow becomes unsteady, and periodic vortex shedding is observed as indicated in [42]. Vorticity contours are shown in Fig. 10 for a single instant in time to provide a qualitative indication of the unsteady nature of the ﬂow, where vortex sheds downstream periodically away from the top and bottom part of the cylinder.

Fig. 7. Contours of u velocity for different Reynolds number.

3

Drag force

2 1 0 -1 -2 -3

0

2

4

6 Time

Fig. 8. Time evolution of drag force for Re = 800.

8

10

160

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

u= v= 0 inlet

Stationary solid cylinder u= v= 0

Lc

outlet

H

Fluid

D

Hc u= v= 0

L

zoom-in showing uniform grid in neighborhood of solid object

Fig. 9. 2-D ﬂow over asymmetically-placed cylinder; L ¼ 3:3; H ¼ 0:615; D ¼ 0:15; Lc ¼ Hc ¼ 0:225.

Fig. 10. Instantaneous vorticity contours near stationary cylinder, dotted and solid lines denote negative and positive contours.

4

3.28

(a) 3.26

drag coefficient

Drag coefficiences

Lift & Drag coefficiences

3

(b)

2 lift coefficient

1

3.24

3.22

0 3.2 -1 6

7

8

Time

9

10

3.18

6

7

8

Time

Fig. 11. Time evolution of (a) drag and lift coefﬁcients and (b) zoom-in drag coefﬁcient for Re ¼ 100.

9

10

161

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

The drag coefﬁcient is deﬁned as

CD ¼ 1 2

FD U 2av g D

ð18Þ

;

where F D is the drag force computed from Eq. (15). Similarly, the lift coefﬁcient is deﬁned as

CL ¼ 1

FL

U2 D 2 av g

ð19Þ

;

where F L is the lift force calculated from Eq. (16). Corresponding variations of lift and drag coefﬁcients are presented in Fig. 11 where effects of asymmetric geometric can be seen in the slightly non-symmetric drag coefﬁcient variations. As was true of the vorticity contours, these results provide a mainly qualitative assessment of performance of the method in the absence of comparisons with other corresponding results.

Table 2 Values of C D;max ; C L;max and St for 2-D ﬂow over a cylinder asymmetrically placed in a channel with Re ¼ 100.

Schafer and Turek [42] Interpolation Interpolation plus solid-body forcing

C D;max

C L;max

St

3.22–3.24 3.255 3.255

0.99–1.01 0.9913 0.9913

0.295–0.305 0.303 0.303

Pressure isolines

Table 2 summarizes values of C D;max ; C L;max and St at Re ¼ 100 for ﬂow over a 2-D cylinder asymmetrically placed as in the channel of Fig. 9 and presents a more quantitative veriﬁcation of solution accuracy. The present results are shown to be comparable with those from previous numerical studies [35,42]. The predicted Strouhal number (St ¼ Dfq =U av g , where fq is the shedding frequency) is 0.303. This agrees well with the range of Strouhal numbers given in [42] using different numerical procedures. For forcing points located in the ﬂuid domain, two different approaches are used for the simulations, viz., interpolation and interpolation with solid-body forcing. From Table 2, we see that the results of the two approaches are identical. In other words, it appears that use of solid-body forcing is not necessary for this stationary-boundary problem, as should be expected since no computational points are moving from solid to ﬂuid as time progresses. On the other hand, this problem provides an easy test of solid-body forcing in a trivial case, showing that its use does not cause errors in situations when motion of the solid body ceases.

3.3. In-line oscillating cylinder in ﬂuid at rest The interaction of an oscillating circular cylinder with a ﬂuid at rest is frequently reported in the literature [28,29,32,43]. The key dimensionless parameters in this ﬂow are the Reynolds number, Re ¼ U max D=m and the Keulegan–Carpenter number, KC ¼ U max =fD (U max and f are maximum velocity and characteristic frequency, respectively, of the cylinder.) The set up in the present work corre-

Vorticity isolines

0

o

96

o

192

o

288

o

Fig. 12. Pressure and Vorticity contours at four different phase angles, dotted and solid lines denote negative and positive contours.

162

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

sponds to the LDA experiment and numerical simulations with non-inertial reference frame reported by Dutsch et al. [43]. Motion of the cylinder is prescribed by

X CEN ðtÞ ¼

1 sinð2pftÞ; 2p

ð20Þ

where X CEN denotes location of the center of the cylinder. The Reynolds and Keulegan–Carpenter numbers are set to Re ¼ 100 and KC ¼ 5, respectively. The computational domain is 55D and 35D in the x and y directions, with the cylinder initially located in the center of the domain. A non-uniform grid of 250 200 cells is adopted to discretize the overall computational domain; within this grid is a reﬁned uniform grid of 80 80 cells employed in the region D 6 x; y 6 D, but with the additional grid lines extending into the remainder of the domain only in each separate direction. Here the CFL number is about 0.55.

Simulations are started in a static ﬂuid ﬁeld, and the cylinder moves according to Eq. (20) until periodic vortex shedding is established. Pressure and vorticity contours at four phase angles after vortex shedding begins are shown in Fig. 12. When the cylinder moves to the left, upper and lower boundary layers develop on the leading face the cylinder and separate at symmetric (top and bottom) positions downstream. Two counter-rotating vortices are produced by the separating ﬂow. These vortices persist until the cylinder reaches the left-maximum position and starts its rightward motion; then the same vortex structure is created on the opposite side of the cylinder. It is interesting to note that this reversed motion carries the cylinder back through already-stirred ﬂuid, resulting in splitting of vortex pairs produced during the preceding leftward motion. These results are consistent with the corresponding laboratory experiments reported in [43], implying that qualitative features of the vorticity ﬁeld dynamics can be properly captured by the present method.

(a) 0.2

0.2

0.1

-0.1

Y

0

Y

0.1

exp.(x=-0.6D) exp.(x=0D) exp.(x=0.6D) exp.(x=1.2D) (x=-0.6D) (x=0D) (x=0.6D) (x=1.2D)

-0.1

-0.2

-0.2 -1.5

0

-1

-0.5

0 U

0.5

1

-1.5

1.5

-1

-0.5

0 V

0.5

1

1.5

-1

-0.5

0 V

0.5

1

1.5

-1

-0.5

0 V

0.5

1

1.5

0.2

0.2

0.1

0.1

0

0

Y

Y

(b)

-0.1

-0.1

-0.2

-0.2 -1.5

-1

-0.5

0 U

0.5

1

-1.5

1.5

0.2

0.2

0.1

0.1

0

0

Y

Y

(c)

-0.1

-0.1

-0.2

-0.2 -1.5

-1

-0.5

0 U

0.5

1

1.5

-1.5

Fig. 13. Comparison of velocity proﬁles at four cross sections with constant x location; phase angles: (a) 180°, (b) 210° and (c) 330°. Experimental results taken from [43].

163

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

2 Present calculations Experimental data, Dutsch et al.[41]

Fx

1

0

-1

-2

0

0.2

0.4

t/T

0.6

0.8

1

Fig. 14. Comparison of x-direction force computed at Re ¼ 100 and KC ¼ 5.

To further compare the present numerical results with experimental data, velocity proﬁles at four x locations and three different phase angles (for each) are plotted in Fig. 13 along with corresponding experimental results from [43]. It is expected that the velocity distributions should be symmetric and anti-symmetric for the u and v velocities, respectively, along the y direction. There is generally good agreement, although measurements show slightly nonsymmetric velocity distribution in the y direction, indicating the level of measurement uncertainty. Fig. 14 shows time evolution of the drag force produced by both computation and experiment; it is clear that present predictions agree very well with those of [43], although force is generally a rather weak measure of accuracy. On the other hand, this is often one of the most important physical quantities required in analysis of bodies moving in a ﬂuid, and these results demonstrate that forces on the surface of the cylinder can be predicted accurately by the present solid-body forcing method with a relatively simple algorithm and quite acceptable computational cost. 3.4. In-line oscillating cylinder in a free stream In order to further explore the ability of the proposed technique for computing moving-boundary ﬂows, simulations are conducted with an in-line oscillating cylinder in uniform ﬂow at Re ¼ 100. According to Tanida et al. [44], Grifﬁn and Ramberg [45], and Hurlbut et al. [46], in-line cylinder oscillation in a range of frequencies near twice the Strouhal shedding frequency for a stationary cylinder causes synchronization, i.e., phase-locking of vortex shedding, with the (oscillatory) motion of the cylinder. Both drag and lift coefﬁcients increase in the presence of such phase-locked oscillations. The present ﬂow domain and boundary conditions are show in Fig. 15 with the cylinder now oscillating parallel to the free stream at a frequency fc ¼ 2f q . Motion of the cylinder is prescribed by setting its horizontal velocity to

cells employed to more accurately compute forces on the cylinder when using the interpolation plus solid-boundary-forcing technique, as was described earlier in discussions associated with Fig. 9. Results are summarized in the following table and ﬁgure. Table 3 shows a comparison of drag and lift coefﬁcients with those computed by Su et al. [30], and Hurlbut et al. [46] for the in-line oscillating cylinder at Re ¼ 100. Results from [24] for a non-oscillating cylinder are also provided. The lock-in effect can be observed by the vortex patterns over two periods of drag coefﬁcient oscillation shown in Fig. 16; here T is the period of the drag coefﬁcient, and the period of the lift coefﬁcient is 2T. Due to the lock-in effect, the contours indicate an essentially periodic behavior with period equaling that of the forced oscillations. As noted above, the cylinder is oscillating at a frequency at which lock-in occurs. Here, the predicted lift amplitude is approximately three times that obtained with a rigid cylinder at the same Reynolds number, and there is an obvious increase in mean drag over the rigid-cylinder value. This is consistent with results from previous experiments and numerical studies, demonstrating that the present method can capture the important practical ﬂow properties very well for moving-body problems.

dp/dy=0 du/dy=0 v=0

solid cylinder

dp/dx=0 H

u=1

D

fluid

dp/dx=0 du/dx=0 dv/dx=0

v=0

U c ðtÞ ¼ 2pfc A cosð2pfc tÞ; where A ¼ 0:028 is oscillation amplitude in the x direction. In the present calculations L ¼ 11; H ¼ 14; D ¼ 0:2 and Re ¼ 100. In Hurlbut et al. [46], the distance from the cylinder to the outer boundary is ﬁfty-ﬁve cylinder diameters. For these calculations the CFL number is about 0.45 and a non-uniform grid of 400 330 cells is employed to discretize the entire computational domain. Within this is a reﬁned uniform-grid region, 2 D 6 x; y 6 2 D, consisting of 160 160

dp/dy=0 du/dy=0 v=0 L Fig. 15. Computational domain and boundary conditions for ﬂow with in-line oscillating cylinder.

164

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

Table 3 Comparisons of lift and drag coefﬁcients of in-line oscillating cylinder at Reynolds number 100. Re ¼ 100

Present

fc =fq C D;mean C L;max

0 1.36 0.34

Su et al. [30] 2 1.71 0.95

0 1.40 0.34

Hurlbut et al. [46] 2 1.70 0.97

0 1.41 0.31

Kim et al. [24] 2 1.68 0.95

0 1.33 0.32

– – –

Y

0.5 0

-0.5

(a)

(e)

(b)

(f)

(c)

(g)

(d)

(h)

Y

0.5 0

-0.5

Y

0.5 0

-0.5

Y

0.5 0

-0.5

0

1

2 X

3

0

1

2

3

X

Fig. 16. Instantaneous vorticity contours near oscillating cylinder, Re=100; dotted and solid lines denote, respectively, 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 period of drag coefﬁcient).

3.5. Two cylinders moving with respect to each other This example was originally investigated numerically in [10], and subsequently via a reformulation in [11]. The initial (with respect to cylinder location) geometry is shown in Fig. 17. The far-

ﬁeld boundaries are rigid walls. To mitigate effects of an impulsive start, both cylinders are oscillated about their initial positions for two periods, and then toward the other, at speeds corresponding to Re ¼ 40. The motion of lower and upper cylinders is given by

Fig. 17. Geometry for ﬂow around two cylinders moving with respect to each other: L ¼ 32; H ¼ 16; Lc ¼ Hc ¼ 8; D ¼ 1; Hs ¼ 1:5, with all values being dimensionless.

165

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

(a) 5

y

0

-5

5

(b)

y

0

-5

-15

-10

-5

0 X

5

10

15

Fig. 18. Contours of vorticity (a) cylinders closest together, and (b) cylinders farthest apart, dotted and solid lines denote negative and positive contours.

( xlc ¼

1 4

sin p4t ;

t 16;

4.0

0 6 t 6 16; 16 6 t 6 32;

Xu and Wang [11]

and

3.0

16 p4 sin p4t ; 32 t;

0 6 t 6 16;

16 6 t 6 32:

A non-uniform grid of 1120 680 cells is employed to discretize the overall computational domain. Here the CFL number is about 0.3. Fig. 18 shows vorticity contours at time t ¼ 24, when the two cylinders are closest to each other, and at time t ¼ 32 corresponding to maximum distance between the two cylinders in parts (a) and (b). The present results are in good qualitative agreement with those in [10], suggesting that dynamics of the vorticity ﬁeld can be captured by the present immersed-boundary formulation. Further, more quantitative solution characterizations are presented in terms of temporal evolution of drag and lift coefﬁcients for the upper cylinder for t 2 ½23; 25, the interval during which the two cylinders are closest to each other. These results are in general agreement with those in [11], as shown in Fig. 19. 3.6. Simulation of a sphere settling under gravity In addition to the above 2-D moving-boundary problems, the present formulation will now be applied to a straightforward 3-D problem for which experimental data are available. Consider a sphere falling under gravity in an enclosure completely ﬁlled with ﬂuid, as depicted in Fig. 20. Dimensions of the box are 0:1 0:16 0:1 m3 , and the diameter of the sphere is 0.015 m. The sphere starts its motion at a height Hs ¼ 0:12 m from the bottom of the box. Experimental measurements using PIV have been performed by ten Cate et al. [47] in an enclosure having the foregoing geometry to obtain the trajectory and velocity of the falling sphere, and these investigators also reported corresponding lattice Boltzmann simulations. In a more complex physical situation, Feng

Drag

xuc ¼

Present calculations

2.0

1.0

0.0

1.0

0.5 Lift

(

0.0

−0.5 23

24 Time

25

Fig. 19. Drag and lift coefﬁcients versus time for upper cylinder in ﬂow around two cylinders moving with respect to each other at Re ¼ 40.

et al. [48] combined direct-forcing immersed-boundary and lattice Boltzmann methods to simulate particulate ﬂows, including sedimentation of single and multiple particles.

166

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

0

(a)

-0.04

Ds

u

H g

-0.08

Fluid

Re32.2, simulation Re32.2, exp. Simulation Experiment Re11.6, simulation Re11.6, exp. Re4.1, simulation Re4.1, exp. Re1.5, simulation Re1.5, exp.

Re

Hs

1.5 4.1 11.6 32.2

-0.12

y

z

x

W

(b) 8

W Fig. 20. Schematic of sphere settling in 3-D container of ﬂuid.

6

Case Case Case Case

E1 E2 E3 E4

H/D

Table 4 Fluid properties used in simulations. Re

q ðkg=m3 Þ

l ð103 Ns=m2 Þ

1.5 4.1 11.6 32.2

970 965 962 960

373 212 113 58

In simulations reported here, density of the falling sphere is 3 1120 kg=m ; ﬂuid density is in the range 960–970 kg=m3 , and ﬂuid dynamic viscosity is varied between 0.058 and 0:353 103 Ns=m2 . Speciﬁc ﬂuid properties are in accord with those employed in the experiments of [47], and are summarized in Table 4. In the present study, a non-uniform grid of ð90 90Þ cells is used to discretize the domain in x and z directions, within which a ﬁner uniform grid of 40 40 cells is employed in the region D 6 x; z 6 D with the coarser grid utilized in the remaining region analogous to griddings employed in preceding 2-D calculations (see Fig. 9, above). In the y direction, 214 grid cells with uniform spacing are used, and the CFL number is 0.2. Fig. 21a and b show comparisons of current calculations with experimental measurements of falling sphere velocities and trajectories (in terms of normalized height above the ﬂoor of the container), respectively, for different cases corresponding to varying ﬂuid properties reported in [47]. It is evident that the present results for both of these agree reasonably well with the experimental measurements reported in [47] (with some discrepancies similar to those reported for simulations of that reference and in [48]), indicating the accuracy and capability of the present immersedboundary method for solving 3-D moving-boundary problems.

4. Conclusion In the present study, an immersed-boundary method is used for the simulation of 2- and 3-D viscous incompressible ﬂows interacting with moving solid boundaries. This is achieved via direct-momentum forcing on a Cartesian grid by adopting a com-

4

2

0

0

1

2

t

3

4

Fig. 21. Comparisons between present simulation and experimental measurements [47] for (a) falling velocity and (b) sphere height.

bination of a ‘‘solid-body-forcing” strategy applied at solid nodes and interpolation at ﬂuid nodes. Numerical integration is based on a combination of second-order accurate Adams–Bashforth and Crank–Nicholson methods applied to a fractional-step technique implemented with a staggered-grid arrangement of dependent variables. The evaluation of the lift and drag forces is simple as proposed in [18], instead of resorting to the complex integration of the pressure and shear stress distributions along the surface. This provides a simple and yet efﬁcient ways of computing moving object ﬂows. The inﬂuence of the solid body forcing within the solid nodes is ﬁrst examined by computing ﬂow induced by an oscillating cylinder in a stationary square domain, where signiﬁcantly lower amplitude oscillations in computed lift and drag coefﬁcients are obtained compared with those without solidbody-forcing strategy. Grid-function convergence tests also indicate second-order accuracy of this implementation with respect to the L1 norm in time and the L2 norm in space. Further test problems are simulated to examine the validity of the present technique: 2-D ﬂows over an asymmetrically-placed cylinder in a channel, in-line oscillating cylinder in a ﬂuid at rest, in-line oscillating cylinder in a free stream, two cylinders moving with respect to one another, and 3-D simulation of a sphere settling

C.-C. Liao et al. / Computers & Fluids 39 (2010) 152–167

under gravity in a static ﬂuid. All computed results are in generally good agreement with various experimental measurements and with previous numerical simulations. This indicates the capability of the present simple implementation in solving complexgeometry ﬂow problems and the importance of solid body forcing in computing ﬂows with moving solid objects. However, it is also noted the simulated Reynolds numbers ranging from 1 to 800 are low. Therefore, the extension of the present technique to compute high Reynolds number complex ﬂows should be further pursued in future. Acknowledgments This research work was supported by the National Science 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. McDonough acknowledges the support by the National Science Council of Taiwan under grant NSC-96-2811-E007-003 during his stay in Taiwan. References [1] Ye T, Mittal R, Udaykumar HS, Shyy W. An accurate Cartesian gird method for viscous imcompressible ﬂows with complex immersed boundaries. J Comput Phys 1999;156:209–40. [2] Udaykumar HS, Mittal R, Shyy W. Computation of solid–liquid phase fronts in the sharp interface limit on ﬁxed grids. J Comput Phys 1999;153:535–74. [3] Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating ﬂows with complex moving boundaries. J Comput Phys 2001;174:345–80. [4] Marella S, Krishnan S, Liu H, Udaykumar HS. Sharp interface Cartesian grid method I: an easily implemented technique for 3D moving boundary computations. J Comput Phys 2005;210:1–31. [5] LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefﬁcients and singular sources. SIAM J Numer Anal 1994;31:1019–44. [6] Calhoun D. A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions. J Comput Phys 2002;176:231–75. [7] Li Z, Wang C. A fast ﬁnite difference method for solving Navier–Stokes equations of irregular domains. Commun Math Sci 2003;1:182–98. [8] Li Z, Lai MC. The immersed interface method for the Navier–Stokes equations with singular forces. J Comput Phys 2001;171:822–42. [9] Le DV, Khoo BC, Peraire J. An immersed interface method for the incompressible Navier–Stokes equations in irregular domains. Proceedings of the third MIT conference on computational ﬂuid and solid mechanics, vol. 710. Elsevier Science; 2005. [10] Russell D, Wang ZJ. A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous ﬂow. J Comput Phys 2003;191:177–205. [11] Xu S, Wang ZJ. A immersed interface method for simulating the interaction of a ﬂuid with moving boundaries. J Comput Phys 2006;216:454–93. [12] Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972;10:252–71. [13] Peskin CS. The immersed boundary method. Acta Numer 2002:459–517. [14] Mittal R, Iaccarino G. Immersed boundary methods. Ann Rev Fluid Mech 2005;37:239–61. [15] Goldstein D, Handler R, Sirovich L. Modeling a no-slip ﬂow with an external force ﬁeld. J Comput Phys 1993;105:354–66. [16] Goldstein D, Handler R, Sirovich L. Direct numerical simulation of turbulent ﬂow over a modeled riblet covered surface. J Fluid Mech 1995;302:333–76. [17] Saiki EM, Biringen S. Numerical simulation of a cylinder in uniform ﬂow: application of a virtual boundary method. J Comput Phys 1996;123:450–65. [18] Lai MC, Peskin CS. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J Comput Phys 2000;160:705–19.

167

[19] Lima E Silva ALF, Silveira-Neto A, Damasceno JJR. Numerical simulation of twodimensional ﬂows over a circular cylinder using the immersed boundary method. J Comput Phys 2003;189:351–70. [20] Uzgoren E, Singh R, Sim J, Shyy W. Computational medeling for multiphase ﬂows with spacecraft application. Prog Aerosp Sci 2007;43:138–92. [21] Uzgoren E, Sim J, Shyy W. Marker-based, 3-D adaptive cartesian grid method for multiphase ﬂow around irregular geometries. Commun Comput Phys 2009;5:1–41. [22] Mohd-Yusof J. Combined immersed boundary/B-Spline method for simulations of ﬂows in complex geometries CTR annual research briefs. NASA Ames/Stanford University; 1997. [23] Fadlun EA, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersedboundary methods for three dimensional complex ﬂow simulations. J Comput Phys 2000;161:35–60. [24] Kim J, Kim D, Choi H. An immersed-boundary ﬁnite-volume method for simulations of ﬂow in complex geometries. J Comput Phys 2001;171:132–50. [25] Tseng YH, Ferziger JH. A ghost-cell immersed boundary method for ﬂow in complex geometry. J Comput Phys 2003;192:593–623. [26] Balaras E. Modeling complex boundaries using an external force ﬁeld on ﬁxed Cartesian grids in large-eddy simulations. Comput Fluids 2004;33:375–404. [27] Gilmanov A, Sotiropoulos F. A hybrid Cartesian/immersed boundary method for simulating ﬂows with 3D, geometrically complex, moving bodies. J Comput Phys 2005;207:457–92. [28] Kim D, Choi H. Immersed boundary method for ﬂow around an arbitrarily moving body. J Comput Phys 2006;212:662–80. [29] Yang J, Balaras E. An embedded-boundary formulation for large-eddy simulation of turbulent ﬂows interacting with moving boundaries. J Comput Phys 2006;215:12–24. [30] Su SW, Lai MC, Lin CA. A simple immersed boundary technique for simulating complex ﬂows with rigid boundary. Comput Fluids 2007;36:313–24. [31] Zhang N, Zheng ZC. An improved direct-forcing immersed-boundary method for ﬁnite difference applications. J Comput Phys 2007;221:250–68. [32] Choi Jung-Il, Oberoi RC, Edwards JR, Rosati JA. An immersed boundary method for complex incompressible ﬂows. J Comput Phys 2007;224:757–84. [33] Ghias R, Mittal R, Dong H. A sharp interface immersed boundary method for compressible viscous ﬂows. J Comput Phys 2007;225:528–53. [34] Miller LA, Peskin CS. When vortices stick: an aerodynamic transition in tiny insect ﬂight. J Exp Biol 2004;207:3073–88. [35] Chen DJ, Lin KH, Lin CA. Immersed boundary method based lattice Boltzmann method to simulate 2d and 3d complex geometry ﬂows. Int J Mod Phys C 2007;18:585–94. [36] Chang C, Liu CH, Lin CA. Boundary conditions for lattice Boltzmann simulations with complex geometry ﬂows. Comput Math Appl 2009;58:940–9. [37] Iaccarino G, Verzico R. Immersed boundary technique for turbulent ﬂow simulations. Appl Mech Rev 2003;56:331–47. [38] Kim J, Moin P. Application of a fractional step method to incompressible Navier–Stokes equations. J Comput Phys 1985;59:308–23. [39] Bell JB, Colella P, Glaz HM. A second-order projection method for the incompressible navier–stokes equations. J Comput Phys 1989;85:257–83. [40] Choi H, Moin P. Effects of the computational time step on numerical solutions of turbulent ﬂow. J Comput Phys 1994;113:1–4. [41] Van den Vorst HA, Sonneveld P. CGSTAB, a more smoothly converging variant of CGS, Technical report 90-50. Delft University of Technology; 1990. [42] Schafer M, Turek S. The benchmark problem ﬂow around a cylinder. In: Hirschel EH, editor. Flow simulation with high-performance computer II. Notes in numerical ﬂuid mechanics, vol. 52. Braunschweig: Vieweg; 1996. p. 547–66. [43] Dutsch H, Durst F, Becker S, Lienhart H. Low-Reynolds-number ﬂow around an oscillating circular cylinder at low Keulegan–Carpenter numbers. J Fluid Mech 1998;360:249–71. [44] 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–84. [45] Grifﬁn OM, Ramberg SE. Vortex shedding from a cylinder vibrating in line with an incident uniform ﬂow. J Fluid Mech 1976;75:257–71. [46] 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–22. [47] ten Cate A, Nieuwstad CH, Derksen JJ, Van den Akker HEA. Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity. Phys Fluids 2002;14:4012–25. [48] Feng ZG, Michaelides EE. Proteus: a direct forcing method in the simulations of particulate ﬂows. J Comput Phys 2005;202:20–51.