Numerical simulation of flows past multiple cylinders using the hybrid Local Domain Free Discretization and Immersed Boundary Method

Ocean Engineering xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Numerical simulation of flows past multiple cylinders using the hybrid Local Domain Free Discretization and Immersed Boundary Method Y.L. Wu

⁎

Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China

A R T I C L E I N F O

A BS T RAC T

Keywords: CFD Viscous flow Flow past multiple cylinders Cartesian mesh Immersed Boundary Method Local Domain-Free Discretization

Study of the flow around a group of cylinders (such as risers of the offshore platforms) has been the challenging subject because of its inherent complexity. The research on this problem is important because flow interference between cylinders is responsible for some significant changes in the characteristics of fluid loads. In this paper, the hybrid Local Domain-Free Discretization-Immersed Boundary Method (LDFD-IBM), which combines the strengths of LDFD method and Immersed Boundary Method (IBM), is adopted to solve this problem. The advantages of the LDFD-IBM over conventional IBM are its truly second-order accuracy and simple implementation using a regular fixed Cartesian mesh. In the present investigation, the solver has been applied to simulate the unsteady viscous flow past two cylinders and three cylinders at different configurations in a uniform cross flow. The streamlines and vorticity contours are plotted for the flow visualization. The lift and drag coefficients time histories are presented.

1. Introduction Flow past a system of circular cylindrical structures is a classic flow problem. Since its fundamental and practical significance, a great deal of research work has been carried out to understand the flow behavior of an isolated, single circular cylinder in cross flow. On the other hand, cylindrical structures are usually appeared in groups (or arrays) and find numerous real-life and industrial applications, such as heat exchanger tubes, cooling systems for nuclear power plants, transmission lines, buildings and offshore structures, etc. These structures are subject to shear and oscillatory flows, and therefore may experience flow-induced vibration, which could lead to structural failure or shutdown of industrial facilities under certain conditions. Due to mutual interference between the cylinders in proximity to each other, the flow pattern, the force coefficients, vortex shedding of each cylinder would significantly change. Therefore, it is worthwhile to study such flow problems in order to understand the details of wake behavior when multiple cylinders is placed in a fluid stream. Although the study of flow past multiple cylinders have drawn considerable attention over the last two decades (Zdravkovich, 1968, 1977, 1985; Williamson, 1985; Chang and Song, 1990; Rockwell, 1998; Sumner et al., 1999; Meneghini et al., 2001; Farrant et al., 2001; Lee et al., 2009; Mittal et al., 1997, etc), many of previous research on this topic have been done by experimental studies. This is because the geometrical complexity of this kind of problem poses a challenge for

⁎

numerical study. Traditionally, people used body-fitted mesh generation techniques to solve the flow problems with complex geometries. However, this involves the complicated transformation between physical domain and computational domain, and may not be suitable for complex flow field around multiple cylinders. Many people then turned to use unstructured grid or mesh-free method. Bao et al. (2010) and Thapa et al. (2015) used Finite Element Method (FEM) to study the flow past two and three cylinders on unstructured mesh. Rhie and Chow (1983) and Choi et al. (1993) are among the pioneers who adopted unstructured Finite Volume Method (FVM) to do the simulation for flow past cylinder problems. Harichandan and Roy (2010) used “consistent flux reconstruction” technique in Finite Volume Method (FVM) and simulated flow over multiple cylinders on a collocated unstructured mesh. Ding et al. (2007) applied a mesh free least squarebased Finite Difference Method (FDM) to study the flow field around two circular cylinders arranged in side-by-side and tandem configurations. Although the unstructured mesh method or mesh-free method is more flexible to handle the complex geometries problems, these methods surfer the drawback of high cost of either maintaining the unstructured mesh data base or finding the supporting points for a meshfree node. Another option for these problems is to use the non-conforming mesh methods. Immersed Boundary Method (IBM), first proposed by Peskin in the 1970s (1977), gained great popularity in the past decade. The major advantages of IBM are the simplicity and the easy

Correspondence to: Department of Modern Mechanics, School of Engineering Science, University of Science and Technology of China, Hefei 230027, China. E-mail addresses: [email protected], [email protected].

http://dx.doi.org/10.1016/j.oceaneng.2017.06.045 Received 25 February 2017; Received in revised form 20 May 2017; Accepted 18 June 2017 0029-8018/ © 2017 Elsevier Ltd. All rights reserved.

Please cite this article as: Wu, Y.L., Ocean Engineering (2017), http://dx.doi.org/10.1016/j.oceaneng.2017.06.045

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

implementation for flow problems with complex geometries. The basic idea of IBM is that, the rigid solid boundary, which is represented by a set of Lagrangian points, is treated as a deformable one with high stiffness. A small distortion of the boundary will lead to a force which tends to restore the boundary into its original position. As a result, an additional force fb is included in the momentum equation, as shown in Eq. (1),

∂u 1 + u⋅∇u = −∇p + Δu + fb. ∂t Re

problem, the authors of this paper proposed a hybrid method, named hybrid of Local Domain-Free Discretization and Immersed Boundary Method (LDFD-IBM) (Wu et al., 2012). LDFD-IBM is a delicate combination of LDFD method and Immersed Boundary Method (IBM), and enjoys the merits of both methods. LDFD-IBM uses the same concept of IBM, that is the numerical discretization of the governing equations is applied to the whole Euler mesh, either the mesh nodes in the fluid domain or in solid body. The flow variables at these mesh nodes has the updating values with the time marching. Therefore, the determination of the value of flow variables near the solid boundary can be done by interpolation process rather than less accurate extrapolation. In this way, the nonslip boundary condition is enforced more accurately therefore the penetration of streamlines into solid objects in the conventional IBM (due to inaccurate satisfaction of no-slip boundary conditions) can be avoided. On the other hand, the treatment of boundary condition for pressure at the solid boundary in the LDFD method and other sharp interface immersed boundary methods, which is not a trivial task, is no longer necessary when the governing equation is solved in the whole domain. Compared with conventional IBM, the force term is added in the governing equation implicitly other than explicitly in LDFD-IBM, therefore the procedure of solving the governing equation is similar to the normal N-S solver. In LDFD-IBM, only the updating the velocity at the grid nodes near the immersed boundary is necessary, which is more straightforward as well as efficient. In this paper, the LDFD-IBM is used for the prediction of the flow past two and three cylinders in different Re numbers and different configurations. The adaptive stencil refinement approach proposed by Ding and Shu (2006) is introduced to enhance the efficiency of LDFDIBM method. This paper provides sufficient technical details of LDFDIBM for an interested reader to program. To prove it is a truly secondorder solver, numerical accuracy study of LDFD-IBM has been carried out via simulating the Taylor-Green decaying vortex problem.

(1)

Note that Eq. (1) is solved on the Eulerian mesh, thus the conventional IBM requires a distribution of this body force from the Lagrangian points to the Eulerian mesh, including exterior and interior of the object. Using IBM, the governing equations can be solved by using a regular fixed Cartesian mesh for discretization of flow field. It is clear that the evaluation of the force term in Eq. (1) is crucial for the performance of IBM. Different ways of deriving the forcing term lead to different variants of IBM. Basically, they can be classified into two categories One is the ‘continuous forcing scheme’ (Peskin, 1977; Lai and Peskin, 2000; Griffith and Peskin, 2008), the other is ‘discrete forcing scheme’ (Udaykumar et al., 2001, Kim et al., 2001; Shu et al., 2007; Mittal et al., 2008). The ‘continuous forcing scheme’ uses continuous forcing function, such as a Dirac delta function which the original IBM method (Peskin, 1977) adopted in the information exchange between the Lagrangian and Eulerian points. This method can be implemented with relatively easy, however the no-slip boundary condition is only approximately satisfied and may result in the flow penetration, and less accuracy near the body boundary. The ‘Discrete forcing scheme’ did better job on precisely satisfaction of the no-slip boundary condition at the fluid/solid interface. Fadlun et al. (2000) showed that the force term in Eq. (1) can be imposed directly by using the interpolation formulas together with normal momentum equation. The boundary condition exactly enforced without any dynamical process. However, their interpolation procedure seems to be complex and some errors may introduced. In the work of Shu et al. (2007), a fractional step technique is used, the body force term is calculated implicitly by making a correction in the velocity field. In this way, the no-slip boundary condition is correctly satisfied. However, the numerical accuracy near the immersed boundary is still of first-order. Recently, a number of other immersed boundary methods, such as the sharp interface immersed boundary method (Mittal et al., 2008) and ghost-cell immersed boundary method (Luo et al., 2017) are developed to achieve the second-order or an even higher-order accuracy. These methods can be viewed as one of ‘Discrete forcing scheme’. In such methods, numerical discretization for the fluid grid nodes immediately adjacent to an immersed solid boundary may involve some grid nodes (called ghost nodes) inside the solid object. The values of velocity and pressure at these nodes can be determined via complex interpolation or extrapolation schemes embedded the nonslip boundary condition. Normally, they do not solve the N-S equations for a spurious flow inside the solid bodies therefore are not truly variants of immersed boundary method which proposed by Peskin. Local Domain-Free Discretization (LDFD) method, which was developed by Shu and Wu (2006), belongs to this category. The concept of LDFD method (Shu and Wu, 2006) is based on the mathematical fact that the solution inside the domain can be extended locally to the outside of the domain. With this idea, we can solve the problems with complex geometries on a non-conforming regular Cartesian mesh. The numerical discretization stencil for the fluid grid nodes immediately adjacent to an immersed solid boundary may involve some grid nodes inside the solid object. The values of velocity and pressure at these nodes can be determined by a local extrapolation process, in which the boundary conditions of the curved embedded boundaries are enforced. This treatment allows the LDFD method to flexibly handle flow problems with complex geometry. However, the extrapolation in LDFD may bring numerical errors and instabilities. To overcome this

2. The numerical algorithm based on LDFD-IBM Fig. 1 gives an example of implementing the LDFD-IBM for a curved solid boundary immersed in a fixed Cartesian mesh. In Fig. 1, the mesh nodes inside the physical domain (named interior node) are represented by solid circles (such as B1, C1), while the mesh nodes outside the domain (named exterior nodes) are represented by shaded circles (such as D1, D2). There are two types of mesh nodes in the immediate neighborhood of solid boundary, and we call them dependent points. If the dependent point is located in the flow domain, it is termed an interior dependent point (represented by empty square, such as A1, A2) and if the dependent point is located outside of flow domain, it is termed an exterior dependent point (represented by empty circle, such as E1, E2). The small solid squares denote the

Fig. 1. Configuration of LDFD-IBM.

2

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

up2 = Ub vp2 = Vb

intersection points of mesh lines and the immersed boundary, i.e., the Lagrangian nodes in conventional IBM (such as p1, p2, q1, q2). LDFD-IBM improves the performance of the LDFD with IBM concept. Before introducing the procedure of LDFD-IBM, we start from LDFD. In LDFD (Shu and Wu, 2006), the discretization of the PDE is always performed at the mesh nodes inside the physical domain, including interior nodes and interior dependent nodes (such as A1,B1,C1 in Fig. 1). For the mesh nodes outside the physical domain, the discretization is skipped (such as E2, D2 in Fig. 1). Therefore, the functional values at the nodes outside the domain (E2 and D2) are not defined. However, for the interior dependent points near the boundary such as A1, if the central-difference scheme on the uniform mesh is still used, the discrete form of the first and second order derivatives in the x direction at the position A1 may need the functional values at E2. In LDFD, the functional value at node E2 can be calculated by extrapolation through the local approximate form of the solution, such as using the three local points (p1, A1, B1) to construct a polynomial-based extrapolation at the exterior dependent point E2. On the other hand, to enforce the no-slip boundary condition for pressure at intersection point such as p1 (which is the local normal ∂p 1 ⎡ ∂ 2 (u ) ⎤ momentum equation ( ∂n ) = Re ⎢ 2n ⎥ in the normal direction to the P1 ⎣ ∂n ⎦P1 n ), one needs to use the value at points aa, bb, cc (which are boundary ⇀ the intersection points of the perpendicular line through p1 and the

where Ub, Vb are velocity components of the immersed boundary. The velocity of immersed boundary is supposed to be prescribed during the computation, and is zero for stationary solid boundaries. Note that in Eqs. (6) and (7), the no-slip boundary condition is directly implemented through the velocity at points p1 and p2. The approximation in Eqs. (4) and (5) have second-order of accuracy, and in the context of IBM, to change the local velocity to the values determined by Eqs. (4) and (5) is equivalent to solving the momentum equation with an implicit force. In LDFD-IBM, since E2, D1, D2 are inside the solid body, their velocities are simply assigned to the body velocities, i.e.:

, 2 . The functional values grid lines as shown in Fig. 1) to calculate ∂⇀ n ∂⇀ n at the points aa, bb, cc can be calculated by interpolation between nodes A1 and G (for point aa), nodes H and G (for point bb), and nodes H and F (for point cc). On the contrary, in the LDFD-IBM, the governing equations are discretized at all mesh nodes. Take the discretization of the secondorder derivatives at exterior dependent point E2 as an example. The discrete form of the second-order derivatives in the x and y directions at E2 can be approximated by ∂ 2uE 2 1 = 2 (uA1 − 2uE 2 + uD2 ) + O (h2 ), ∂x 2 h

(2)

∂ 2uE 2 1 = 2 (uA2 − 2uE 2 + uD1) + O (h2 ), ∂y 2 h

(3)

in which the velocities at A1, A2, E2, D1, D2 are needed for the discretization. To enforce the no-slip boundary condition, the velocity at A1 can be evaluated by an interpolation scheme along the x direction, which involves three points p1, B1, C1:

uA1 =

(xA1 − xB1)(x A1− xC1) u (xp1 − xB1)(xp1 − xC1) p1

+

+

(xA1 − xC1)(xA1 − xp1) u (xB1 − xC1)(xB1 − xp1) B1

(xA1 − xp1)(xA1 − xB1) u (xC1 − xp1)(xC1 − xB1) C1

. (4)

Here, B1 and C1 are the interior mesh nodes, and p1 is the intersection point between the horizontal mesh line and the immersed boundary. Similarly, the spatial discretization of the derivatives in the y direction for E2 needs the flow information at A2. The velocity at the point A2 can be updated by the following interpolation scheme that involves three points p2, B2, C2,

uA2 =

( yA2 − yC 2)( yA2 − yp2) ( yA2 − yB2)( y A2−yC 2) up2 + uB2 ( yp2 − yB2)( yp2 − yC 2) ( yB2 − yC 2)( yB2 − yp2)

+

( yA2 − yp2)( yA2 − yB2) ( yC 2 − yp2)( yC 2 − yB2)

(8)

uD1 = Ub vD1 = Vb

(9)

uD2 = Ub vD2 = Vb

(10)

∼ u n +1/3 − un 1 = − H ( u n) − G ( p n ) + L ( u n) Δt /3 Re

(11)

∼ u n +1/2 − un 1 = −H (∼ u n +1/3) − G ( p n ) + L (∼ u n +1/3) Δt /2 Re

(12)

∼ u n +1 − un 1 = −H (∼ u n +1/2) − G ( p n ) + L (∼ u n +1/2) Δt Re

(13)

u* − ∼ u n +1 = Gp n Δt

(14)

1 (D u*) Δt

(15)

Lp n+1 =

uC 2 (5)

un +1 = u* − ΔtG ( p n +1)

where B2, C2 are the interior mesh points, and p2 is the intersection point of a vertical mesh line with the immersed boundary. The velocity at the intersection points p1 and p2 is assigned as the immersed boundary velocity, i.e.,

up1 = Ub , vp1 = Vb

uE 2 = Ub vE 2 = Vb

Obviously, in the above process, the solid boundary condition is implemented accurately. Again in the context of IBM, it is equivalent to solving the momentum equation with an implicit force by changing the local velocity at E2, D1,D2 to the values determined by Eqs. (8)–(10). For the computation of momentum equation at each time step, the velocity at interior dependent points needs to be updated by interpolation at least once, either along the x direction or along the y direction. It should be noted that two distinct approximate solutions may be obtained at some interior dependent points. Taking the point A1 as an example, the velocity can be obtained by interpolation along the x direction, i.e. using B1, C1, p1 via Eq. (4); on the other hand, it can also be calculated by interpolation along the y direction, i.e. using interior nodes R, G, and the boundary point q1. To uniquely determine the value at these interior dependent points such as A1 and A2, we choose to use the average of approximate solutions from the x- and ydirections. In this paper, we adopted that the second-order finite difference scheme for the approximation of spatial derivatives and the explicit three-step formulation (Jiang and Kawahara, 1993) for the approximation of the time derivative. The fractional step method is employed to couple the momentum equations and the continuity equation in solving the N-S equations. The solution procedure is as shown in the following:

∂ 2u

∂p

(7)

(16)

where H denotes the discrete advection operator, G is the discrete gradient operator, L the discrete Laplacian operator, D is the divergence operator. Superscripts n, n+1/3, n+1/2, and n+1 denote the time levels. It should be noted that, when solving Eqs. (11)–(13), the additional force fb in Eq. (1) is added implicitly through updating the velocities at

(6) 3

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

the dependent points. In fact, the velocity at an interior dependent point is computed according to Eqs. (4) and (5), and the velocity at an exterior dependent point or exterior node is set to the wall velocity (Eqs. (8)–(10)). u n+1 has been calculated by Eq. (13), Eq. (14) is used to After ∼ compute the intermediate velocity u*, which is required to solve the pressure Poisson Eq. (15). Once the pressure p n+1 is obtained, the velocity un+1 can be determined by Eq. (16). It should be noted that through Eqs. (15) and (16), the pressure is obtained by solving the governing equation directly, and there is no special treatment of pressure boundary condition along the solid bodies involved. The important advantages of LDFD-IBM over LDFD are: 1) The evaluation of the functional value at some mesh nodes in LDFD-IBM became interpolation other than extrapolation. 2) The pressure is obtained by solving the governing equation directly, therefore, no tedious process is involved to handle with the Neumann boundary condition for the pressure along the normal direction of the solid boundary as that in LDFD. Fig. 2. Position of the embedded circle and the contours of vorticity at t = 0.3 for decaying vortex problem.

3. Numerical accuracy study of LDFD-IBM via decaying vortex problem It is well known that in the original IBM, the distribution of restoring force and the interpolation of flow velocity are through the Dirac delta function, which only has first order accuracy. Later in 2000, Lai and Peskin (2000) presented a formally second-order accurate IBM with adoption of a well-chosen Dirac delta function. However, due to the velocity does not have a continuous derivative across the boundary, therefore the method from Lai and Peskin (2000) is still first-order accurate in practice. In this section, the numerical accuracy of LDFD-IBM is examined. The test example is the Taylor-Green decaying vortex problem. This problem is often used to determine the accuracy of a numerical scheme since the analytical solution is available. Kim et al. (2001) calculated this problem in a quadrilateral embedded domain. Uhlmann (2005) solved this case in an embedded circular domain. We solve this problem in a square domain Ω = [−1.5, 1.5] × [−1.5, 1.5] in which a solid circle with radius of 0.5 is embedded. Reynolds number here is taken as 100. The analytical solutions of the 2-D decaying vortices are

u (x, y, t ) = − cos(πx )⋅ sin(πy)exp(−2π 2t /Re) v (x, y, t ) = sin(πx )⋅ cos(πy)exp(−2π 2t /Re) 1 p (x, y, t ) = − (cos(2πx ) + cos(2πy)exp(−4π 2t /Re) 4

Fig. 3. Convergence rate for decaying vortex problem.

The superscripts 'numerical' and 'exact' represent the values obtained from the numerical results and from analytical solution (Eq. (17)). In Fig. 3, the second order slope and the convergence rate of LDFD method are also included for comparison. To make the comparison a fair one, the mesh nodes outside of the flow domain, i.e, exterior nodes and exterior dependent nodes, are excluded because in LDFD, the functional values at those nodes are not defined. Therefore, N is the total number of interior mesh nodes between the square outer boundary and the circular inner boundary in the flow domain, i.e., the total number of the interior nodes and interior dependent nodes (but excluding the external boundary nodes because the errors at those nodes are zero since we have already assigned the analytical solution to the boundary of square domain). From Fig. 3, it is found that for relatively fine mesh size of 41 × 41–161 × 161, both LDFD and LDFD-IBM have essentially the second order accuracy (the slope of the convergence rate line is about 2.0 for both LDFD and LDFD-IBM). This demonstrates that they are truly second order accurate numerical methods. It is as expected because the second order numerical discretization is applied, and the second order interpolation/extrapolation scheme is used in both LDFD-IBM and LDFD. Decaying vortex problem is a tough problem for the simulation over relatively coarse mesh such as mesh size of 21ⅹ21. From Fig. 3, it is found that the performance of LDFD-IBM on the mesh size of 21 × 21 is better than LDFD method in terms of numerical error. This is because LDFD uses the extrapolation scheme when obtaining the approximate solution at exterior dependent points, while LDFD-IBM uses the interpolation scheme to do the functional value approximation at interior dependent points. It is well known that the accuracy of extrapolation is worse than interpolation when the nodal spacing is large.

(17)

In this test, the analytical solution (17) provides: 1) initial field at t = 0; 2) The time-dependent boundary conditions at the boundary of square domain; 3) The time-dependent variable values at the circumference of the embedded circle; 4) the exact solution which is used to compute the numerical error so that the overall accuracy of solution can be measured. Several uniform meshes (21 × 21, 41 × 41, 51 × 51, 81 × 81, 101 × 101, 161 × 161) are used in the simulations. The time is advanced from t = 0 to t = 0.3 using a time step of Δt = 0.00001. Fig. 2 shows the position of the embedded circle and contours of vorticity at t = 0.3 with the mesh size of 161 × 161. The solution at time level t = 0.3 is selected to do the comparison. Fig. 3 shows the convergence rate for spatial accuracy. The numerical error of u in terms of L 2 norm is defined as

L 2 error =

∑N (u numerical − uexact )2 N

(18) 4

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

sented by the variation of the drag and lift coefficients at the surface of body, which are defined by

CD =

FD , 1 2D ρU 2

CL =

FL 1 2D ρU 2

(24)

where D is the diameter of cylinder, ρ is the fluid density, and U is the free-stream velocity, FD is the drag force and FL is the lift force acting on the circular cylinder. In this study, the drag force FD and the lift force FL are calculated from the integral of the x- and y-component of the momentum equations over a rectangular control volume Ω enclosing each cylinder as follows: Fig. 4. Configuration of flow past a pair of cylinders.

⎧d ⎛ FD = −⎨ ⎜ ⎩ dt ⎝ ⎪ ⎪

4. Flows over a pair of circular cylinders

⎧d ⎛ FL = −⎨ ⎜ ⎩ dt ⎝ ⎪ ⎪

∫Ω u2 dΩ⎟⎠ + ∫

⎪ ⎪

In this arrangement, three transverse gaps of 1.5D, 3D and 4D are used for the simulations. The computational domain is taken as 50Dⅹ40D. Uniform Cartesian mesh of 201 × 161 is taken as the initial mesh. To get satisfactory resolution around the boundaries of two cylinders, we use stencil mesh refinement technique as in the paper of Ding and Shu (2006) to perform 10 levels of mesh refinement, and the final mesh contains 66,739 nodes, 70,621 nodes, 71,637 nodes for cases of T=1.5D, 3D, 4D, respectively. The mesh spacing near the cylinder is about 0.0078D. The non-dimensional time step used for the calculation is 0.002. Fig. 5(a) shows the local refined mesh around the two cylinders for T=3D. Fig. 6 and Fig. 8 present instantaneous vorticity contours and streamlines for the flow around two cylinders in side-by-side arrangement with a transverse gap of 1.5D at Re = 100 and Re = 200. As indicated by Zdravkovich (1985), T = 1.5D is within the range of intermediate critical gap (1.1D < T < 2.2D), where the occurrence of biased flow pattern can be observed in experiments. From our numerical results, it is found that this flow pattern is successfully reproduced. Similarly as reported by Bearman and Wadcock (1973), this biased flow pattern is bistable, i.e., the narrow and wide wakes, and the direction of the gap flow, switch from one cylinder to the other. Figs. 7 and 9 show the time evolution of drag and lift coefficients. "Time" here is non-dimensionalized by speed of incoming flow and reference length. It can be observed that drag and lift coefficients exhibit irregular variation with time, which implies that the irregular vortex shedding occurs behind the pair of cylinders. When T > 2D, there are two types of synchronized Karman vortex streets appearing in the experiments (Zdravkovich, 1977): one is symmetric (anti-phase) and the other is anti-symmetric (in-phase). Both phenomena are captured in our numerical solutions for the transverse gaps of 3D and 4D, as shown in Figs. 10, 12, 14, 16, where the plots of instantaneous vorticity contours and streamlines are presented. In most cases, the anti-phase type of vortex streets occurs in the flow behind two cylinders. The only exception is the case of T = 3D at Re = 100, where the wake appears as clearly in-phase type of vortex. This is consistent with observations reported from other researchers, such as Ding et al. (2007). The time evolutions of drag and lift coefficients are shown in Figs. 11, 13, 15, 17. As shown in these figures, drag and lift coefficients vary synchronously with time, which also confirm the synchronized behavior of vortex shedding from the

(19)

(20)

(21)

At the far field downstream boundary of the cylinder, the zero gradient boundary condition is applied by

(22)

At beginning of numerical simulation, the undisturbed velocity field at t=0 is assumed as the initial flow field, i.e.

⎧ u t =0 = 1.0 ⎨ ⎩ v t =0 = 0.0

⎞

⎡ 1 ⎛ ∂u2 ∂u ⎞ ⎤ ⎫ ⎢u2 u⋅n + pn2 − + i ⎟ ⎥ ds⎬ ⎜ Re ⎝ ∂xi ∂x2 ⎠ ⎦ ⎭ ∂Ω ⎣

4.1. Side-by-side arrangement

On the top and bottom boundary, the symmetric boundary condition is applied, that is,

⎧ ∂u = 0 ⎪ ∂x ⎨ ⎪ ∂v = 0 ⎩ ∂x

⎪

(26)

On the cylinder surface, no-slip boundary condition is implemented

∂u ⎧ ⎪ =0 ∂y ⎨ ⎪ ⎩v = 0

⎪

where ∂Ω is the control surface of Ω , n = (n1, n2 ) is the normal vector to the boundary of the control surface, the subscripts 1 and 2 denote the x-direction and y-direction, respectively.

by

⎧u = 0 ⎨ ⎩v = 0

1 ⎛ ∂u1 ∂u ⎞ ⎤ ⎫ + i ⎟ ⎥ ds⎬ ⎜ Re ⎝ ∂xi ∂x1 ⎠ ⎦ ⎭ (25)

The problem configuration is shown in Fig. 4, which is an incompressible, viscous fluid flow at a constant velocity U∞ = 1 past a pair of stationary cylinders with same diameter D. Three standard classifications of the configuration of a pair of cylinders are side-byside, tandem and staggered arrangements. There are three parameters which can significantly affect the flow structures: the Reynolds number Re, the transverse gap spacing T, and the longitudinal gap L, as shown U D in Fig. 4. Here Re is defined as Re= ∞ν , where D is the cylinder diameter, and ν is the kinematic viscosity. Experiments showed that there exists a range of flow regimes characterized by the critical gap. For example, for the side-by-side arrangements, Williamson (1985) found that three flow patterns can be observed for T < 2.2D, 2.2D < T < 5D, T > 5D,respectively. Zdravkovich (1977) identified three flow regimes for the tandem arrangement, i.e., L < 1.2–1.8D, 1.2–1.8D < L < 3.4–3.8D and L > 4D. These critical gaps may vary in a certain range according to the Reynolds number. In the present study, the numerical simulations are performed for the cases with different transverse and longitudinal gaps for the sideby-side and tandem arrangements at Re = 100 and Re = 200. The computational domain is taken as 50Dⅹ40D. The outer boundary is set in order not to cause any perturbation in the near-wake flow. The boundary conditions of the problem are: At the inlet, a free stream velocity profile is specified, that is,

⎧u = 1 ⎨ ⎩v = 0

⎡

⎞

∫Ω u1dΩ⎟⎠ + ∫∂Ω ⎢⎣u1u⋅n + pn1 −

(23)

The time-evolution of the flow field can be quantitatively repre5

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 5. Local refined mesh for simulation of flow past a pair of circular cylinders. (a) Side by side arrangement (T = 3D) (b) Tandem arrangement (L = 5.5D).

Fig. 6. Vorticity contours and streamlines for flow over a pair of side-by-side cylinders (T = 1.5D) at Re = 100.

Fig. 7. Time histories of drag (CD) and lift (CL) coefficients of flow past a pair of side-by-side cylinders (T = 1.5D, Re = 100).

Fig. 8. Vorticity contours and streamlines for side-by-side cylinders (T = 1.5D) at Re = 200.

upper and lower cylinders. The drag coefficient CD , lift coefficient CL and Strouhal number St (St = fD / U , where f is the shedding frequency) for the case of T = 1.5D,

3D and 4D at Re = 100, 200 are presented in Tables 1–3, together with some previous results available in the literature. Our numerical results on an isolated single cylinder with Re = 100 and 200 are also included

6

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 9. Time histories of CD and CL for side-by-side cylinders (T = 1.5D, Re = 200).

Fig. 10. Vorticity contours and streamlines for flow over a pair of side-by-side cylinders (T = 3D) at Re = 100.

Fig. 11. Drag and lift coefficients of flow past a pair of side-by-side cylinder (T = 3D) at Re = 100.

Fig. 12. Vorticity contours and streamlines for side-by-side cylinders (T = 3D) at Re = 200.

7

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 13. Drag and lift coefficients of side-by-side cylinders (T = 3D) at Re = 200.

Fig. 14. Vorticity contours and streamlines for side-by-side cylinders (T = 4D) at Re = 100.

Fig. 15. Drag and lift coefficients of side-by-side cylinder (T = 4D) at Re = 100.

Fig. 16. Vorticity contours and streamlines for side-by-side cylinders (T = 4D) at Re = 200.

8

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 17. Drag and lift coefficients of side-by-side cylinder (T = 4D) at Re = 200. Table 1 Comparison of drag coefficients, lift coefficients and Strouhal number for flow past two side-by-side cylinders T=3D for Re=100. Results Ding et al. (2007)

Upper cylinder Lower cylinder Upper cylinder Lower cylinder Upper cylinder Lower cylinder Upper cylinder Lower cylinder Single cylinder

Chang and Song (1990) Lee et al. (2009) Present results Wu et al. (2012)

CD

CL

St

1.560 ± 0.038 1.560 ± 0.038 1.533 ± 0.04 1.533 ± 0.04 1.48 ± 0.04 1.48 ± 0.04 1.457 ± 0.04 1.457 ± 0.04 1.330 ± 0.01

−0.131 ± 0.253 0.131 ± 0.253 −0.108 ± 0.31 0.108 ± 0.31 −0.135 ± 0.465 0.135 ± 0.465 −0.120 ± 0.31 0.120 ± 0.31 ± 0.341

0.182 0.182 0.18 0.18 – – 0.181 0.181 0.166

4.2. Tandem arrangement

Table 2 Comparison of drag coefficients, lift coefficients and Strouhal number for flow past two side-by-side cylinders T = 3D for Re = 200.

Ding et al. (2007) Present results Wu et al. (2012)

Upper cylinder Lower cylinder Upper cylinder Lower cylinder Single cylinder

CD

CL

St

1.548 ± 0.083 1.548 ± 0.083 1.536 ± 0.110 1.536 ± 0.110 1.30 ± 0.041

−0.104 ± 0.866 0.104 ± 0.866 −0.11 ± 0.911 0.11 ± 0.911 ± 0.69

0.215 0.215 0.221 0.221 0.200

In this arrangement, the longitudinal gaps of 2.5D, 5.5D are considered. The computational domain is again taken as 50Dⅹ40D. Uniform Cartesian mesh of 201 × 161 is taken as the initial mesh. 10 levels of mesh refinement are performed by stencil mesh refinement technique (Ding and Shu, 2006), and the final mesh contains 71,123 nodes, 72,607 nodes for cases of T = 2.5D, 5.5D, respectively. The mesh spacing near the cylinder is about 0.0078D. The non-dimensional time step used for the calculation is 0.002. Fig. 5(b) shows the local refined mesh around the two cylinders for T = 5.5D. Figs. 18 and 19 present instantaneous vorticity contours and streamlines for L = 2.5D at Re=100 and Re = 200. For the case of Re = 100, it can be observed that the flow maintains a steady state though the Reynolds number is greater than the critical value for one isolated cylinder (Re≈49 ). The shear layer separated from the upstream cylinder reattaches on the surface of the downstream cylinder. This phenomenon is called "quasi-steady reattachment" by Zdravkovich

in Tables 1 and 2. It can be found that the present results are consistent with previous data. From the figures and tables, it is also found that when the gap between cylinders increase, the values of the flow parameters for each cylinder tend towards the numerical results for a single cylinder. It indicates that the flow interference reduces as the two cylinders are placed further apart.

Table 3 Comparison of drag coefficients, lift coefficients and Strouhal number for flow past two side-by-side cylinders T = 4D for Re = 100.

Ding et al. (2007) Farrant et al. (2001) Lee et al. (2009) Bao et al. (2010) Present

Upper Lower Upper Lower Upper Lower Upper Lower Upper Lower

cylinder cylinder cylinder cylinder cylinder cylinder cylinder cylinder cylinder cylinder

CD

CL

St

1.514 ± 0.013 1.514 ± 0.013 1.44 ± 0.03 1.44 ± 0.03 1.44 ± 0.03 1.44 ± 0.03 1.429 1.429 1.435 ± 0.026 1.435 ± 0.026

0.077 ± 0.348 −0.077 ± 0.348 0.08 ± 0.32 −0.08 ± 0.32 0.079 ± 0.415 −0.079 ± 0.415 0.071 ± 0.288 −0.071 ± 0.288 0.07 ± 0.422 −0.07 ± 0.422

0.184 0.184 0.172 0.171 – – – – 0.183 0.183

9

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 18. Vorticity and streamlines for tandem cylinders (L = 2.5D) at Re = 100.

Fig. 19. Vorticity and streamlines for tandem cylinders (L = 2.5D) at Re = 200.

Fig. 20. Drag and lift coefficients of tandem cylinders (L = 2.5D) at Re = 200.

Fig. 21. Vorticity and streamlines for tandem cylinders (L = 5.5D) at Re = 100.

(1985). When Reynolds number increases, the flow field behind the upstream cylinder remains steady, i.e., quasi-steady attachment, but the flow behind the downstream cylinder becomes unsteady, as shown in Fig. 19. The flow behind the downstream cylinder at Re = 200 finally evolves into an oscillatory pattern with obvious periodicity. The corresponding drag and lift coefficients for Re = 200 are shown in Fig. 20. It can be seen that the drag coefficients of both upstream and downstream cylinders remain almost constant, but the lift coefficients vary with time. When the longitudinal gap is increased from 2.5D to 5.5D, the flow patterns are completely different from those for the cases of L=2.5D, as shown in Figs. 21 and 23. The separated shear layer from the upstream cylinder no longer reattaches to the front of downstream cylinder. Instead, a Karman street is formed behind the upstream cylinder. The temporal histories of drag and lift coefficients are shown in Figs. 22 and 24. The figures show that vortex shedding frequencies for the two

cylinders are identical and the shedding behind the cylinders is synchronized nearly in anti-phase. This “lock-in” phenomenon confirms the findings of experiments by Rockwell (1998). The drag coefficient CD , lift coefficient CL and Strouhal number St for the case of L = 2.5D, 5.5D at Re = 100, 200 are presented in Table 4. Some numerical and experimental results from the literature are also included in this table. The agreement here also is very good, showing the applicability of the present algorithm to study the problem of flow past multiple cylinders. 5. Flow over three circular cylinders From the literature survey, we found that some researchers (Harichandan and Roy, 2010; Kang, 2003) have reported that numerical results for three circular cylinders in side-by-side or tandem arrangements. Although three cylinders in equidistant triangular 10

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 22. Drag and lift coefficients of tandem cylinders (L = 5.5D) at Re = 100.

Fig. 23. Vorticity and streamlines for tandem cylinders (L = 5.5D) at Re = 200.

Fig. 24. Drag and lift coefficients of tandem cylinders (L = 5.5D) at Re = 200. Table 4 Comparison of drag coefficients, lift coefficients and Strouhal number for flows over a pair of tandem cylinders at Re = 100 and Re = 200. Results

L = 2.5D Re=100

L = 2.5D Re=200

L = 5.5D Re=100

Present Ding et al. (2007) Mittal et al. (1997) Wu (2010) Present Wu (2010) Ding et al. (2007) Present Ding et al. (2007) Mittal et al. (1997) Wu (2010)

St

CL

CD Upstream cylinder

Down- stream

Upstream cylinder

Down- stream

1.135 1.163 1.271 1.146 0.996 1.02 1.022 1.223 1.329 1.433 1.301

−0.0868 −0.0895 −0.075 −0.077 −0.181 −0.18 −0.22 0.759 ± 0.858 ± 0.952 ± 0.761 ±

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

± 0.015 ± 0.013 ± 0.015 ± 0.012

11

0.121 0.125 0.164 0.139

± ± ± ± ± ± ±

0.024 0.023 0.002 0.367 0.330 0.403 0.374

± ± ± ± ± ± ±

0.223 0.24 0.023 1.307 1.554 1.741 1.527

– – – – 0.194 0.19 0.19 0.156 0.16 0.168 0.154

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 26 shows the local refined mesh around the three cylinders for Types I, II, III. All the simulations are conducted at Re = 100 with undisturbed initial field by Eq. (23). The instantaneous vorticity contours and streamlines for three types of configuration are presented in Figs. 27, 29, 31, 32. The drag coefficient CD , lift coefficient CL for all the cases of flow past three cylinders are presented in Figs. 28, 30, 33, 34. Table 5 compares the drag and lift coefficients for all three types (where the fluctuation from the mean value is taken as the average of peak-to-peak amplitude to the mean force coefficient). For Type I arrangement, it is found from Fig. 27 that the vortex shedding developed behind Cylinders #1 and #2 is suppressed due to the presence of Cylinder #3. Due to the shear layer of upstream cylinders rolling into the wake vortex of downstream cylinder, a relatively complete vortex shedding is formed behind the downstream cylinder (Cylinder #3). The flow patterns behind upstream cylinders #1 and #2 show nearly symmetric feature with respect to the incoming stream direction. The symmetric flow feature of Type I arrangement is demonstrated by the force coefficients for Cylinder #1 and #2 in Fig. 28, while the mean lift coefficients of Cylinder #1 and #2 are contrary in sign while almost equal in magnitude. Lam and Cheung (1988) and Gu and Sun (2001) reported that this symmetric phenomenon is also exists at higher Re in their experiments. This symmetricity may be the reason for the mean value of the life coefficient of Cylinder #3 close to zero. The mean drag coefficient of downstream cylinder #3 is smaller than those of upstream Cylinder #1 and #2 (Fig. 28), due to the shielding effect from Cylinder #1 and #2. The fluctuating drag and lift coefficients of downstream cylinder are much larger than those of upstream cylinders, which means the impact which comes from the wake of upstream cylinders on the downstream cylinder is severe, but there is phase difference in vortex shedding behind upstream cylinders #1 and #2. CD and CL histories as shown in Fig. 28 indicates this is true. For arrangement Type II, the prominent feature of the flow pattern is an individual vortex shed behind Cylinder #1, indicating that the wake interference from Cylinder #2 and #3 are weak. Because Cylinders #2 and #3 are in the tandem arrangement, quasi-steady attachment occurs in the flow field behind Cylinder #2 as described in Section 4.2. The vortex shedding behind Cylinder #3 is evidently disturbed by the wake of Cylinder #1. The shielding effect from Cylinder #1 and #2 causes the mean drag coefficient of Cylinder #3 is much smaller than those of Cylinder #1 and #2 (Fig. 30). The proximity effect causes the drag force of Cylinder #1 is larger than an isolated single cylinder, and it also suppresses the fluctuation of lift coefficients of Cylinder #2 and #3. For arrangement Type III, unlike the finding of Bao et al. (2010) where only one type of synchronized flow is observed, both in-phase

Fig. 25. Configuration of different types of arrangement for flow past three cylinders.

arrangement is a fundamental arrangement, there are not much numerical study have been performed. Zdravkovich (1968) used a smoke-visualization technique to study the complex interaction of the vortex streets for Reynolds number between 60 and 300 experimentally. Lam and Cheung (1988) studied the flow patterns and vortexshedding frequency data of three cylinders arranged equidistantly at different angles of incidence and different spacing ratios at Reynolds numbers of 2.1 × 103 and 3.5 × 103 using a dye-injection technique in their experiments. Dong et al. (1999) investigated the frequency of vortex shedding experimentally from three cylinders of same diameter in equilateral triangular arrangement at a Reynolds number of 3.0 × 104. Gu and Sun (2001) investigated the flow pattern of this problem at Reynolds number of 5.5 × 104 and 1.4 × 104 in cross uniform flow with smoke-wire technique flow visualization. Bao et al. (2010) employed the second-order characteristic-based split scheme in FEM and investigated the low Reynolds number flow over three equal circular cylinders arranged in equilateral-triangle arrangement numerically. In this paper, LDFD-IBM is applied to numerical simulations over three circular cylinders in equidistant triangular arrangements at low Reynolds number. In this study, three configurations of three identical circular cylinders in equilateral triangular array are investigated, as shown in Fig. 25. The distance between circular cylinder centers l is 3.0D, where D is the diameter of the cylinder. The reason that l = 3.0D being selected as the focus of this study is because according to Lam and Cheung (1988), when spacing ratio l/D < 1.43, the group of three cylinders can essentially be treated as a single large cylinder. When spacing ratio l/D reached 4.65, the effects of flow interference became negligible thereafter. The computational domain is taken as 50Dⅹ40D. Uniform Cartesian mesh of 201×161 is taken as the initial mesh. 10 levels of mesh refinement are performed around the cylinders by stencil mesh refinement technique, and the final mesh contains 81,061 nodes for Type I, 81,125 nodes for Type II, 81,061 nodes for Type III, respectively. The mesh spacing near the cylinder is about 0.0078D. The non-dimensional time step used for the calculation is 0.002.

Fig. 26. Local refined mesh for simulation of flow past three circular cylinders.

12

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 27. Instantaneous vorticity contours and streamlines for Type I (Re = 100).

Fig. 28. Drag and lift coefficients of flow past three cylinders (Type I) at Re = 100.

Fig. 29. Instantaneous vorticity contours and streamlines for Type II (Re = 100).

Fig. 30. Drag and lift coefficients of flow past three cylinders (Type II) at Re = 100.

13

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Fig. 31. Instantaneous vorticity contours and streamlines for Type III (anti-phase).

Fig. 32. Instantaneous vorticity contours and streamlines for Type III (in-phase).

Fig. 33. History of lift coefficients of three cylinders (Type III) at Re = 100.

Fig. 34. Drag and lift coefficients of flow past three cylinders (Type III) at Re = 100 (In-phase). Table 5 Comparison of drag coefficients, lift coefficients for flows over three cylinders at Re = 100 with l/D = 3.0. Type

I II III

CD

CL

#1

#2

#3

#1

#2

#3

1.244 ± 0.036 1.489 ± 0.065 1.285 ± 0.033

1.245 ± 0.036 1.134 ± 0.040 1.045 ± 0.004

1.157 ± 0.143 0.308 ± 0.165 1.284 ± 0.033

−0.130 ± 0.084 −0.109 ± 0.492 0.056 ± 0.311

0.131 ± 0.085 0.077 ± 0.075 ± 0.031

−0.003 ± 0.527 0.086 ± 0.272 0.061 ± 0.311

14

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

into in-phase wake. The drag force and amplitude of the lift coefficients for the two-cylinder combination in Type III are larger than that in Type I. This indicates when the third cylinder is placed upstream, the two-cylinder combination (Cylinder #1 and #3 in Type III) experiences larger drag force and stronger life force oscillation in the flow than the situation (Cylinder #1 and #2 in Type I) when the third cylinder placed downstream. In general, from the numerical study of the flow past three cylinders in different equilateral arrangements at intermediate spacing (l/D = 3.0) under Re = 100, it is found that interference effect causes changed fluid loads acting on the cylinders. The vortex shedding from the upstream cylinder(s) always seems to be suppressed by the downstream cylinder(s). This results in the drag of the cylinders in combination is always less than the sum of the drags of the isolated cylinders, which is due to the proximity effect and wake effect of flow over multiple cylinders.

and anti-phase synchronized wake patterns behind the downstream cylinders are observed during the single simulation(see Figs. 31 and 32). At the initial stage of the simulation, the flow pattern is completely symmetric with respect to the free stream direction, as shown in Fig. 31 (at t = 110). An anti-phase synchronized vortex shedding can be observed behind downstream cylinders #1 and #3. It is also found that the flow behind the Cylinder #2 is in a steady state, and symmetric recirculation zone is attached behind it. However, the wake of Cylinder #2 sandwiched by the effect of downstream cylinders could lead to unstable flow structures. As shown in the history of lift coefficients (Fig. 33), one can see clearly that the flow around Cylinder #2 becomes more and more unstable. At around t = 220, it triggers the flow deflecting from the state of anti-phase vortex shedding. An in-phase synchronized flow behind the downstream cylinders can be observed, as shown in Fig. 32 (at t = 500). This is an interesting finding that both types of synchronized Karman vortex streets appear in present numerical study, and this needs to be verified by experiments. The temporal drag coefficient CD , lift coefficient CL are presented in Fig. 34 when the flow eventually settles into in-phase vortex shedding. The mean value of drag coefficients CD of all cylinders are smaller than that for a single cylinder, reflecting that the flow interference has significant effect on the reduction of the drag force acted on the cylinders. The upstream cylinder #2 experiences lower CD and CL than that for Cylinder #1 and #3, and a recirculation zone is attached behind upstream cylinder #2.

Acknowledgements Support from “the Fundamental Research Funds for the Central Universities” (WK2090000009) is gratefully acknowledged. References Bao, Y., Zhou, D., Huang, C., 2010. Numerical simulation of flow over three circular cylinders in equilateral arrangements at low Reynolds number by a second-order characteristic-based split finite element method. Comput. Fluids 39, 882–889. Bearman, P.W., Wadcock, A.J., 1973. The interaction between a pair of circular cylinders normal to a stream. J. Fluid Mech. 61, 499–511. Chang, K.S., Song, C.J., 1990. Interactive vortex shedding from a pair of circular cylinders in a transverse arrangement. Int. J. Numer. Methods Fluids 11, 317–329. Choi, S.K., Nam, H.Y., Cho, M., 1993. Use of the momentum interpolation method for numerical solution of incompressible flows in complex geometries: Choosing cell face velocities. Numer. Heat. Transf. Part B 23, 21–41. Ding, H., Shu, C., 2006. A stencil adaptive algorithm for finite difference solution of incompressible viscous flows. J. Comput. Phys. 214, 397–420. Ding, H., Shu, C., Yeo, K.S., Xu, D., 2007. Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods. Int. J. Numer. Methods Fluids 53, 305–332. Dong, L.Y., Xu, Y.H., Huang, D.Q., 1999. Vortex shedding frequency determined by surface pressure fluctuations of three equi-spaced cylinders. J. Hydrodyn. 4, 75–84. Fadlun, E.A., Verzicco, R., Orlandi, P., 2000. Combined Immersed-boundary finitedifference methods for three dimensional complex flow simulations. J. Comput. Phys. 161 (1), 35–60. Farrant, T., Tan, M., Price, W.G., 2001. A cell boundary element method applied to laminar vortex shedding from circular cylinders. Comput. Fluids 30, 211–236. Griffith, B.E., Peskin, C.S., 2008. On the order of accuracy of the immersed boundary method: higher order convergence rates for sufficiently smooth problems”. J. Comput. Phys. 208, 75–105. Gu, Z., Sun, T., 2001. Classification of flow pattern on three circular cylinders in equilateral-triangular arrangements. J. Wind Eng. Ind. Aerodyn. 89 (6), 553–568. Harichandan, A.B., Roy, A., 2010. Numerical investigation of low Reynolds number flow past two and three circular cylinders using unstructured grid CFR scheme. Int. J. Heat Fluid Flow 31, 154–171. Jiang, C.B., Kawahara, M., 1993. A three-step finite element method for unsteady incompressible flows. Comput. Mech. 11, 355–370. Kang, S., 2003. Characteristics of flow over two circular cylinders in a side-by-side arrangement at low Reynolds numbers. Phys. Fluids 15 (9), 2486–2498. Kim, J., Kim, D., Choi, H., 2001. An immersed-boundary finite-volume method for simulations of flow in complex geometries. J. Comput. Phys. 171, 132–150. Lai, M., Peskin, C.S., 2000. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J. Comput. Phys. 160, 705–734. Lam, K., Cheung, W.C., 1988. Phenomena of vortex shedding and flow interference of three cylinders in different equilateral arrangements. J. Fluid Mech. 196, 1–26. Lee, K., Yang, K.S., Yoon, D.H., 2009. Flow-induced forces on two circular cylinders in proximity. Comput. Fluids 38, 111–120. Luo, K., Mao, C., Zhuang, Z., Fan, J., Haugen, N.E.L., 2017. A ghost-cell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions Part-II: Complex geometries. International Journal of Heat and Mass Transfer 104, 98–111. Meneghini, J.R., Saltara, F., Siqueira, C.L.R., Ferrari, J.A., Jr., 2001. Numerical simulation of flow interference between two circular cylinders in tandem and sideby-side arrangements. J. Fluids Struct. 15, 327–350. Mittal, R., Dong, H., Bozkurttas, M., Najjar, F.M., Vergas, A., von Loebbecke, A., 2008. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J. Comput. Phys. 227 (10), 4825–4852. Mittal, S., Kumar, V., Raghuvanshi, A., 1997. Unsteady incompressible flows past two cylinders in tandem and staggered arrangements. Int. J. Numer. Methods Fluids 25,

6. Conclusions In this study, the LDFD-IBM method, which combines the merits of LDFD and IBM, is applied to simulate flows around multiple circular cylinders in different configurations and Reynolds numbers. As a Cartesian mesh method, the mesh generation in LDFD-IBM is straightforward. Our numerical test shows that the accuracy of LDFD-IBM method is truly second order. The good agreement between present numerical results and benchmark solutions shows the good accuracy and robustness of LDFD-IBM. This paper focused on flow around two and three cylinders within the low Reynolds number range, i.e. Re=100 or 200. For the flow around two cylinders in tandem and side-by-side arrangement, it has been found that the flow features significantly depend on gap spacing between cylinders and Re number. The flow patterns behind the cylinders can be divided into several categories depending on the gap between the cylinders, as in experimental observation by Zdravkovich (1977, 1985) and by Williamson (1985), have been confirmed in numerical investigation here. The flow characteristics which have been observed in the experiments, such as two types of synchronized Karman vortex streets for T > 2D in the side-by-side configuration and “reattachment” phenomenon for values of L < 3D in tandem arrangement, are successfully reproduced. The present numerical results for two cylinders in tandem and side-by-side arrangements compared well with the available data from the literature. Although three-cylinder arrangement can be considered as twocylinder combination with one more cylinder, the flow characteristic is changed significantly owing to the third cylinder. For example, the wake behind the configuration of two cylinders in side-by-side arrangement with T = 3D in Re = 100 is in-phase type of vortex street as shown in Fig. 10. However, when the third cylinder is placed downstream (Type I), the anti-symmetric (in-phase) vortex shedding does not exist anymore. On the contrary, the flow exhibits symmetric feature. Both drag and lift force for the upstream two-cylinder combination decrease compared with original two cylinders placed side-by side in the stream. In the configuration of one upstream cylinder sandwiched by the two downstream cylinders (Type III), the flow firstly experiences anti-phase flow pattern and eventually settles 15

Ocean Engineering xxx (xxxx) xxx–xxx

Y.L. Wu

Udaykumar, H.S., Mittal, R., Rampunggoon, P., Khanna, A., 2001. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J. Comput. Phys. 174, 345–380. Uhlmann, M., 2005. An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209, 448–476. Williamson, C.H.K., 1985. Evolution of a single wake behind a pair of bluff bodies. J. Fluid Mech. 159, 1–18. Wu, J., 2010. Development of a Novel Immersed Boundary-lattice Boltzmann Method and its Applications (Ph.D. Thesis). National University of Singapore. Wu, Y.L., Shu, C., Ding, H., 2012. Simulation of incompressible viscous flows by local DFD-immersed boundary method’. Adv. Appl. Math. Mech. 4 (3), 311–324. Zdravkovich, M.M., 1968. Smoke observation of the wake of a group of three cylinders at low Reynolds number. J. Fluid Mech. 32, 339–351. Zdravkovich, M.M., 1977. Review of flow interference between two circular cylinders in various arrangements. ASME J. Fluids Eng. 99, 618–633. Zdravkovich, M.M., 1985. Flow induced oscillations of two interfering circular cylinders. J. Sound Vib. 101, 511–521.

1315–1344. Peskin, C.S., 1977. Numerical analysis of blood flow in the heart. J. Comput. Phys. 25, 220. Rhie, C.M., Chow, W.L., 1983. Numerical study of turbulent flow past an airfoil with trailing edge separation. AIAA J. 21 (11). Rockwell, D., 1998. Vortex-body interactions. Ann. Rev. Fluids Mech. 30, 199–229. Shu, C., Liu, N.Y., Chew, Y.T., 2007. A novel immersed boundary velocity correctionlattice Boltzmann method and its application to simulate flow past a circular cylinder. J. Comput. Phys. 226, 1607–1622. Shu, C., Wu, Y.L., 2006. Adaptive mesh refinement-enhanced local DFD method and its application to solve Navier-Stokes equations. Int. J. Numer. Methods Fluids 51, 897–912. Sumner, D., Wong, S.S.T., Price, S.J., Paidoussis, M.P., 1999. Fluid behavior of side-byside circular cylinders in steady cross-flow. J. Fluids Struct. 13, 309–338. Thapa, J., Zhao, M., Leng, L., Zhou, T., 2015. Three dimensional simulations of flow past two circular cylinders in side-by-side arrangements at right and oblique attacks. J. Fluids Struct. 55, 64–83.

16