Journal of Computational Physics 189 (2003) 351–370 www.elsevier.com/locate/jcp Numerical simulation of two-dimensional ﬂows over a circular cylinder using the immersed boundary method A.L.F. Lima E Silva a, A. Silveira-Neto a a,* , J.J.R. Damasceno b Mechanical Engineering College, Federal University of Uberl^ andia, Jo~ ao Naves de Avila Avenue, Campus Santa M^ onica, Minas Gerais, Brazil b Chemical Engineering College, Federal University of Uberl^ andia, Jo~ ao Naves de Avila Avenue, Campus Santa M^ onica, Minas Gerais, Brazil Received 9 May 2002; received in revised form 14 November 2002; accepted 23 March 2003 Abstract In this work, a virtual boundary method is applied to the numerical simulation of a uniform ﬂow over a cylinder. The force source term, added to the two-dimensional Navier–Stokes equations, guarantees the imposition of the no-slip boundary condition over the body–ﬂuid interface. These equations are discretized, using the ﬁnite diﬀerences method. The immersed boundary is represented with a ﬁnite number of Lagrangian points, distributed over the solid–ﬂuid interface. A Cartesian grid is used to solve the ﬂuid ﬂow equations. The key idea is to propose a method to calculate the interfacial force without ad hoc constants that should usually be adjusted for the type of ﬂow and the type of the numerical method, when this kind of model is used. In the present work, this force is calculated using the Navier–Stokes equations applied to the Lagrangian points and then distributed over the Eulerian grid. The main advantage of this approach is that it enables calculation of this force ﬁeld, even if the interface is moving or deforming. It is unnecessary to locate the Eulerian grid points near this immersed boundary. The lift and drag coeﬃcients and the Strouhal number, calculated for an immersed cylinder, are compared with previous experimental and numerical results, for diﬀerent Reynolds numbers. Ó 2003 Elsevier Science B.V. All rights reserved. Keywords: Immersed-boundary method; Eulerian–Lagrangian formulation; Drag coeﬃcient; Finite-diﬀerence method 1. Introduction Immersed bodies present a very broad spectrum of practical problems due to their geometric complexity and movement inside the ﬂow. This class of problems is very diﬃcult to analyze using classical methods, because of these two characteristics. The most common problems involve vibration of the immersed * Corresponding author. Tel.: +55-34-323-9-4148. E-mail addresses: [email protected] (A.L.F. Lima E Silva), [email protected] (A. Silveira-Neto), [email protected] (J.J.R. Damasceno). 0021-9991/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0021-9991(03)00214-6 352 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 structure, which induces turbulence in the ﬂow and turbulence induces vibration on the immersed structure. This is the well-known ﬂuid–structure interaction problem. The movement of one or several particles, which are in a deposition process inside a reservoir, presents another example of this class of problems . Recently, some new types of methods based on the immersed boundary idea of Peskin , have been proposed. The main idea is to use a regular Eulerian mesh for the ﬂuid dynamics simulation, coupled with a Lagrangian representation of the immersed boundary. The immersed boundary exerts a singular force on the ﬂuid. The Lagrangian mesh can move independent of the Eulerian mesh. Diﬀerent models to calculate this force ﬁeld have been proposed recently. Lai and Peskin  have proposed a second-order accurate immersed boundary method that has been used to simulate ﬂows over a circular cylinder. The interaction between the ﬂuid and the immersed boundary is modeled using a well-chosen discretized approximation to the Dirac delta function. The singular Lagrangian force intensity is calculated using a Hook law, in which the boundary point X is ﬁxed to a equilibrium point X e with a very stiﬀ spring, whose stiﬀness constant is j 1. Goldstein et al. [6,7] proposed a model called the Virtual Boundary Formulation, which was used to simulate turbulent ﬂow over a riblet-covered surface and to address other similar problems. Their method is similar to the immersed boundary method but is appropriate only to rigid immersed surfaces. The main idea of the virtual boundary method is to treat the embedded boundary in the ﬂuid, applying a force ﬁeld to the ﬂuid, so that it takes the same velocity as the surface. This is a virtual model for the classic no-slip boundary condition. Since the force ﬁeld is not known a priori, it must be calculated in a feedback way such that the velocity on the boundary is used to determine the desired force distribution. This model involves two imposed constants, a and b, which are chosen to be both negative and large enough in magnitude to bring the ﬂuid velocity close to the interface velocity. The force intensity is calculated using the expression ~ f ð~ xk ; tÞ ¼ a Z th i h i ~ ð~ ~ ð~ ~ xk ; tÞ dt þ b ~ uð~ xk ; tÞ U xk ; tÞ ; uð~ xk ; tÞ U ð1Þ 0 where a and b are constants that must be adjusted in order to obtain the expected physical behavior at the ﬂow. Saiki and Biringen  used this model to simulate a two-dimensional ﬂow around stationary and moving cylinders. Fogelson and Peskin  solved the three-dimensional Stokes equations to simulate the sedimenting particle ﬂows in the Stokes regime. In this work, the particles are slightly deformable. A regular Lagrangian grid is used to represent the overall domain. Each particle is represented by a group of Lagrangian points, placed over the ﬂuid–particle interface. The momentum equations are solved at all points of the Eulerian grid. The particle points move with the local velocity ﬂuid. This velocity is aﬀected by density force term that is calculated over the Lagrangian points and distributed over the neighboring Eulerian points. The distributed force ﬁeld has to model the non-linear and viscous eﬀects between ﬂuid and particle surfaces. Knowing this Eulerian force ﬁeld, the entire domain is viewed as one entity without solid interfaces. The authors made simulations with one and two particles under the inﬂuence of gravity. They observed favorable behavior as compared with existing theories for multi particle sedimentation. They concluded that the proposed method is a powerful numerical tool to study ﬂows with suspended particles. Ye et al.  proposed a method called the Cartesian Grid Method for simulating two-dimensional unsteady, viscous, incompressible ﬂows over complex geometries. In this method, a control volume near the immersed boundary is re-formed into a body-ﬁtted trapezoidal shape by discarding the solid part of the cell and adding the neighboring cells. They used an interpolation scheme near the immersed boundary that gives second-order spatial accuracy. The concept of momentum forcing, that is, the force added on the Navier–Stokes equations, is not used. A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 353 Fadlum et al.  used discrete-time forcing, as suggested by Mohd-Yusof , and showed that it is more eﬃcient than the feedback forcing method for three-dimensional problems. The velocity at the ﬁrst grid point, external to the interface, is obtained by a linear interpolation of the velocity at the second grid point and the velocity of the body surface. The force is calculated at the ﬁrst external grid point using the Navier– Stokes Equations. In this model, there are no adjustable constants. This method requires a special algorithm to identify the grid points over which the force must be calculated. Kim et al.  proposed a method based on a ﬁnite-volume approach, using a staggered mesh. Flows over immersed complex geometries were simulated. The authors introduced the use of a mass source or sink and a momentum forcing to guarantee the no-slip boundary condition on the immersed boundary, and to satisfy the continuity for the cell containing the immersed boundary. In the present work, a new model is proposed. It is based on the calculation of the force ﬁeld over a sequence of Lagrangian points, which represent the interface, using the Navier–Stokes Equations. There are no ad hoc constants in this model and a special algorithm to capture the neighboring grid points of the immersed interface is not required. This model is being named here Physical Virtual Model (PVM), due to the fact that it is based only upon the laws of conservation. We simulated an internal channel ﬂow and the ﬂows around a circular cylinder, in order to validate the methodology and the numerical procedure. Detailed informations about the ﬂows over the cylinder, at diﬀerent Reynolds numbers are presented. These ﬂow quantities are the mean drag and mean lift coeﬃcients, the drag and lift ﬂuctuations, the wall pressure and the recirculation length behind the cylinder and the Strouhal number, for low Reynolds numbers. The results were compared with experimental and with other numerical studies. 2. Mathematical formulation 2.1. General aspects of the formulation Viscous and incompressible ﬂows in a Cartesian square domain X containing an immersed boundary, as illustrated in Fig. 1, can be modeled by the Navier–Stokes equations: " # oV~ ~ÞV ~P þ lr2 V~ þ ~ ~ ¼ r q F; ð2Þ þ ðV~ r ot ~ V~ ¼ 0; r ð3Þ Fig. 1. Immersed boundary illustration; Eulerian mesh ð~ xÞ and Lagrangian mesh ð~ xk Þ. 354 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 where ~ F is given by Z ~ F ð~ xÞ ¼ ~ x ~ xk Þ d~ xk ; f ð~ xk Þdð~ ð4Þ X and dð~ x ~ xk Þ is a Dirac delta function, ~ xk are the Lagrangian points placed over the immersed boundary (see Fig. 1), ~ f ð~ xk Þ is the Lagrangian force density, and ~ F ð~ xÞ is the Eulerian force that is not equal to zero only over the immersed boundary. Eq. (4) models the interaction between the immersed boundary and the ﬂuid ﬂow, by injecting the force ﬁeld on the ﬂuid. In order to discretize the Dirac delta function that appears in Eq. (4), it must be replaced by the distribution function, Dij , proposed by Peskin  and modiﬁed by Juric , represented by the following equations: Dij ð~ xk Þ ¼ f ½ðxk xi Þ=h f ½ðyk yj Þ=h ; h2 8 < f1 ðrÞ f ðrÞ ¼ 12 f1 ð2 krkÞ : 0 ð5Þ if krk < 1; if 1 < krk < 2; if krk > 2; ð6Þ where f1 ðrÞ ¼ 3 2krk þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ 4krk 4krk2 8 ; ð7Þ r being ðxk xi Þ=h or ðyk yj Þ=h, h is the Eulerian grid size and ðxi ; yj Þ are the Eulerian points. An important property of this function is that its integral is equal to one. 2.2. Immersed boundary force – Physical virtual model In the present work, the PVM is proposed to calculate the Lagrangian force ﬁeld, based only upon the momentum equation. All the Navier–Stokes terms are calculated over the Lagrangian points. Therefore the force ~ f ð~ xk Þ should be expressed by ~ f ð~ xk Þ ¼ ~ fa ð~ xk Þ þ ~ fi ð~ xk Þ þ ~ fv ð~ xk Þ þ ~ fp ð~ xk Þ: ð8Þ The diﬀerent terms that compose Eq. (8) are here referred to as acceleration force ~ fa , inertial force ~ fi , ~ ~ viscous force fv and pressure force fp . These force components are represented by: oV~ ~ ð~ xk Þ; xk Þ ¼ q fa ð~ ot ð9Þ ~ ~ÞV~ð~ fi ð~ xk Þ ¼ qðV~ r xk Þ; ð10Þ ~ xk Þ; xk Þ ¼ lr2 V~ð~ fv ð~ ð11Þ ~ ~P ð~ xk Þ ¼ r xk Þ: fp ð~ ð12Þ A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 355 Fig. 2. Scheme for the pressure and velocity interpolation. The diﬀerent terms described by Eqs. (9)–(12) must be evaluated over the interface using the velocity ﬁeld V~ð~ xÞ and pressure ﬁeld P ð~ xÞ. These calculations must also take into account that, at the interface, the ﬂuid velocity must be equal to the interface velocity, which guarantees the no-slip condition. The velocity and pressure spatial derivatives are calculated using the ﬂuid velocity and pressure obtained by Eqs. (2) and (3). ~ð~ One of the possible ways to do this is to interpolate V xÞ and P ð~ xÞ over appropriate points near the interface, as illustrated by Fig. 2. This interpolation method is described in the following section. 2.3. Interpolation and derivatives calculation procedure The interpolation scheme for the horizontal velocity, over point number 3, for example, is shown in Fig. 3. A square dotted domain is deﬁned as being the region where the variables will be analyzed. In the Fig. 3. Interpolation scheme for the horizontal velocity at point 3. 356 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 4. Interpolation scheme for the vertical velocity at point 3. shaded area that denotes the immersed body, the velocity and pressure ﬁelds do not contribute to the interpolation process. The horizontal velocity components, shown in Fig. 3, are interpolated over point 3 using the same model given by Eqs. (5)–(7). Only those components, which are at a distance less than or equal to 2Dx, from the point number 3, will be interpolated. This procedure is executed automatically by the distribution function model. Furthermore, the dotted box, shown in Figs. 4 and 5, is used to reduce CPU time calculation. Fig. 5. (a) Interpolation scheme for the velocities and (b) for the pressure at point k. A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 357 In Fig. 4, the vertical velocity components are shown for the interpolation scheme over point 3. The same procedure is used to obtain the interpolated pressure value. The interpolations over the other points (1, 2 and 4), are also done the same way. For the points over the interface (k points), the procedure is diﬀerent. The velocity components are interpolated using the internal and external grids (Fig. 5(a)). The pressure is taken as the value at a Dx distance from point k, in the normal direction, as shown in Fig. 5(b). Point p belongs to the nearest external Eulerian cell from point k. The pressure and velocity derivatives, that appear in Eqs. (9)–(12), are calculated using a secondorder Lagrange polynomial approximation. Generalizing the vertical and horizontal velocity components and the pressure by /, the ﬁrst and second derivatives in the x direction are approximated by: o/ ðxk x2 Þ ðxk x1 Þ ðxk x1 Þ þ ðxk x2 Þ ðxk ; yk Þ ¼ / þ / þ /; ox ðx1 x2 Þðx1 xk Þ 1 ðx2 x1 Þðx2 xk Þ 2 ðxk x1 Þðxk x2 Þ k ð13Þ o2 / 2/1 2/2 2/k þ þ : ðxk ; yk Þ ¼ ox2 ðx1 x2 Þðx1 xk Þ ðx2 x1 Þðx2 xk Þ ðxk x1 Þðxk x2 Þ ð14Þ The derivatives in the y direction are given by: o/ ðyk y4 Þ ðyk y3 Þ ðyk y3 Þ þ ðyk y4 Þ ðxk ; yk Þ ¼ / þ / þ /; oy ðy3 y4 Þðy3 yk Þ 3 ðy4 y3 Þðy4 yk Þ 4 ðyk y3 Þðyk y4 Þ k ð15Þ o2 / 2/3 2/4 2/k þ þ ; ðxk ; yk Þ ¼ oy 2 ðy3 y4 Þðy3 yk Þ ðy4 y3 Þðy4 yk Þ ðyk y3 Þðyk y4 Þ ð16Þ where /1 , /2 , /3 and /4 are obtained by the interpolation described in Figs. 3–5. The interface velocity components ðuk ; vk Þ, are also used for the derivative calculation. These velocity components are equal to zero for the purposes of the present work, for which the interface is stationary. The pairs ðxk ; yk Þ, ðx1 ; y1 Þ, ðx2 ; y2 Þ, ðx3 ; y3 Þ and ðx4 ; y4 Þ are the coordinates of the points k, 1, 2, 3 and 4, respectively, as shown in Fig. 2. We emphasize that the derivatives expressed by Eqs. (13)–(16) are calculated over the interface at points ðxk ; yk Þ. The distances between the points k and 1, k and 3, 1 and 2, and 3 and 4 are Dx (the mesh size). Their distribution is always at the outside part of the immersed boundary. Using the signal of the normal vector components, it is possible to conveniently locate these points, outside the interface. Therefore, the calculation of the force terms, Eqs. (10)–(12), is independent of the internal ﬂow properties. The acceleration force qðoV~=otÞ, Eq. (9), is calculated taking into account that the ﬂuid velocity over the interface must have the same value as the interface velocity. Therefore, this acceleration term is approxi~k ¼ ðuk ; vk Þ is the interface velocity and V ~fk ¼ ðufk ; vfk Þ is the ﬂuid velocity mated by ðqðV~k V~fk ÞÞ=Dt, where V at the same position on the interface. Once calculated, the Lagrangian force, given by Eq. (8), is distributed over the neighboring Eulerian grid, as illustrated in Fig. 6. The Dirac delta function, that appears in Eq. (4), is replaced by the distribution function, in order to calculate the Eulerian force, in a discrete form. Therefore Eq. (4) is replaced by X ~ Fij ¼ Dij~ ð17Þ fk Ds2 ; where Dij is given by Eqs. (5)–(7) and Ds is the distance between two Lagrangian points. The presented mathematical model is solved using the numerical method described in the following section. 358 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 6. Illustration of the force distribution process. 3. Numerical method Eqs. (2) and (3) were discretized using the second-order ﬁnite central diﬀerence method in space and a Euler ﬁrst order in time. The velocity and pressure coupling was solved using a second-order pressure correction method, as presented by Armﬁed and Street . This method is described as follows: 1 opn uinþ1 uni ¼ N ðuni Þ þ Lðuni Þ þ Fin ; q oxi Dt ð18Þ o2 unþ1 q o uinþ1 ¼ ; Dt oxi oxi oxi ð19Þ uinþ1 ¼ uinþ1 Dt ounþ1 ; q oxi pnþ1 ¼ pn þ unþ1 ; ð20Þ ð21Þ where Lðui Þ ¼ l o2 ui ; q oxj oxj ð22Þ N ðui Þ ¼ oðui uj Þ ; oxj ð23Þ where ui is the estimated velocity component, u is the pressure correction and Dt and n are the computational time step and the sub-step indices, respectively. The numerical solution of these equations is derived as follows: 1. Calculate the force ﬁeld ~ f ð~ xk Þ, over the Lagrangian points ðxk ; yk Þ, using Eqs. (8)–(12) and the initial conditions. 2. Distribute the force ~ f ð~ xk Þ to the Eulerian grid using Eq. (17). 3. Estimate the ﬂuid velocity ﬁeld under the inﬂuence of the force ﬁeld ~ F using Eq. (18). A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 359 4. Solve the linear system for the pressure correction given by Eq. (19). 5. Compute the actual pressure ﬁeld using Eq. (21). 6. Update the ﬂuid velocity using Eq. (20). 7. Verify the divergence using Eq. (3). This procedure completes the time step loop. Normally, one loop is suﬃcient to obtain the mass conservation, veriﬁed in step (7). The convergence criterion used for all the simulations was that ~ V ~g 6 108 , where the max value is calculated over the entire domain. maxfr The resulting linear system for the pressure correction u is solved using the interactive solver MSI (Modiﬁed Strongly Implicit Procedure) of Schneider and Zedan . Note that the interface force ﬁeld calculation and the momentum equation solution are performed in an explicit way. 4. Results The immersed boundary method, using the Physical Virtual Model, has been applied to simulate an internal channel ﬂow and the ﬂow around a circular cylinder in order to validate the methodology and the numerical procedure. For the two-dimensional Poiseuille ﬂow simulation, the channel walls where modeled using the PVM model. Simulations of the ﬂow past a circular cylinder were done for diﬀerent low Reynolds numbers and consequently, the drag and lift coeﬃcients, the length of the recirculation bubble and the pressure coeﬃcient on the cylinder surface were obtained. These results were compared with results present by other authors. 4.1. Poiseuille ﬂow Fig. 7 shows the steady state of a vertical velocity proﬁle in a vertical channel of width h and length L. This simulation was performed using a 100 200 grid in the x and y directions, respectively. A constant pressure gradient was imposed on the channel and the force ﬁeld was calculated at the channel walls. This Fig. 7. (a) Force ﬁeld at the channels walls. (b) Detail of the force at the wall. (c) Vertical velocity proﬁle at y ¼ 2. The continuous line represents the analytical result and the crosses the numerical result. 360 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 result was compared with the result obtained by the analytical solution. We would like to stress that boundary conditions over the walls were modeled using the immersed boundary method, instead of the classic no-slip velocity condition. The maximum velocity error obtained comparing the present result with analytical solution, was 1.0%. 4.2. Flow over a circular cylinder A rectangular domain was used to simulate the ﬂow over a stationary cylinder. The boundary conditions were imposed in such a way that the ﬂow was from the bottom up toward the top of the domain. A circular cylinder was placed inside the domain so that its center had coordinates x ¼ 7:5d and y ¼ 16:5d. The domain had a length of 30d and a width of 15d. These dimensions were chosen in order to minimize the boundary eﬀects on the ﬂow development. A Newman boundary condition was used on the lateral boundaries. A constant velocity proﬁle U1 was speciﬁed at the domain entrance. A grid reﬁnement study was done to verify the result independence and the accuracy of the method. The results are shown in Table 1 for ﬁve grid sizes. They became grid independent for a grid of 150 300 points and all results presented in this paper were obtained for a grid of 250 500 points. Fig. 8 shows the 250 500 uniform mesh used in the simulations of ﬂow past a stationary circular cylinder of diameter d ¼ 0:1. Simulations were performed at Re ¼ 10, 20, 40, 50, 60, 80, 100, 150, and 300. Table 1 A grid reﬁnement study Re ¼ 80 Dx ¼ Dy Grid points CD 0.0150 0.0120 0.0100 0.0075 0.0060 100 200 125 250 150 300 200 400 250 500 1.360 1.409 1.394 1.398 1.396 Fig. 8. (a) Complete view of the domain with the circular cylinder placed inside a vertical channel and (b) zoom view of the Cartesian and Lagrangian grids. A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 361 The velocity was imposed in such a way that the desired Reynolds number was obtained. The ﬂuid properties q and l were considered constants and were the same for all simulations. The criterion for the determination of the number of Lagrangian points was such as Ds=Dx 6 0:9. This implies that, in the present method, we can use a reduced number of Lagrangian points, as compared with Saiki and Biringen , who used a reﬁned Lagrangian mesh in their simulations. The Reynolds number in this ﬂow was deﬁned as qU1 d ; l Re ¼ ð24Þ where d is the cylinder diameter. Finally, the dimensionless time was deﬁned as T ¼ U1 t : d ð25Þ Once the velocity and pressure ﬁelds are calculated, the drag and lift coeﬃcients and the Strouhal number can also be calculated, using the force ﬁeld directly. Drag coeﬃcient. As demonstrated by Lai and Peskin , the drag force on an immersed body in a stream, arises from two sources: the shear stress and the pressure distribution along the body. The drag coeﬃcient can be deﬁned as CD ¼ FD ; 2 d ð1=2ÞqU1 ð26Þ where FD is the drag force. Here, it can be calculated using the Eulerian force or the Lagrangian force as: Z Z L FD ¼ Fx dx ¼ fx ds; ð27Þ X 0 where Fx and fx are the x components of the Eulerian and Lagrangian force, respectively, L is the length of the interface and X is a domain which contains the interface. Lift coeﬃcient. Similar to the drag coeﬃcient, the lift coeﬃcient can be deﬁned as CL ¼ FL ; 2 d ð1=2ÞqU1 ð28Þ where FL is the lift force, which can be calculated using the y component of the Eulerian force or the Lagrangian force as Z Z L FL ¼ Fy dy ¼ fy ds; ð29Þ X 0 where Fy and fy are the y components of the Eulerian and Lagrangian forces, respectively. Strouhal number. The Strouhal number is deﬁned as the dimensionless frequency with which the vortices are shed behind the body St ¼ fd ; U1 ð30Þ where f is the vortex shedding frequency. This frequency can be obtained using the Fast Fourier Transform of the lift coeﬃcient time distribution. A small time step ðDt ¼ 1 106 Þ was used at the ﬁrst loop, in order to guarantee the code stability. It is increased as the ﬂow develops. When the force ﬁeld is enough to satisfy the no-slip condition on the cylinder 362 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 surface, the time step is of 1 103 . This time step is of the order of that used by others authors (i.e., Lai and Peskin ). The results obtained with this time step, for Re ¼ 80, compare very well with others authors results. Although, by security measure the authors have used a time step of Dt ¼ 1:104 for all calculation presented in the present paper. For Reynolds number equal to 10, 20 and 40 the wake formed behind the cylinder attains a steady symmetric state, which is in agreement with the well-established result, by the linear stability theory. The cylinder wake instabilities rises for Re P 47. Fig. 9 shows the recirculation bubble behind the cylinder at Re ¼ 10, 20 and 40, respectively, in the steady ﬂow regime. The bubble increases with the Reynolds number. The force components ðFx ; Fy Þ distributed near the solid interface are shown in Fig. 10. We observe that it is diﬀerent from zero also inside the interface. In fact, we have a velocity ﬁeld inside this interface but it does Fig. 9. Flow past a circular cylinder. Streamline visualization: (a) Re ¼ 10, (b) Re ¼ 20 and (c) Re ¼ 40. Fig. 10. Force ﬁeld illustration for Re ¼ 40: (a) Fx , (b) Fy . A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 363 not aﬀect the Lagrangian force calculation because only the external ﬁelds are used during the interpolation process. The Lagrangian force terms, i.e., ~ fa (the acceleration force), ~ fi (the inertial force), ~ fv (the viscous force) and ~ fp (the pressure force) are given by Eqs. (9)–(12). Their temporal evolution is shown in Figs. 11–13. In Fig. 11(a) the points A and B represent the positions where this terms are analyzed. Fig. 11. Speciﬁc points used to analyze the Lagrangian force components (a); temporal variation of the vertical acceleration force at point B and vertical component of the velocity, inside the bubble recirculation, at point C (b). Fig. 12. Temporal variation of the Lagrangian force components at Point A: horizontal component and vertical component (b). Fig. 13. Temporal variation of the Lagrangian force components at Point B: horizontal component and vertical component (b). 364 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 11(b) shows the variation of the vertical component of the acceleration force ðfay Þ, at point B, and the temporal variation of the vertical component of the ﬂuid velocity ðvÞ, at point C. The components (x and y) of all force terms extracted at the point A, on the left side of the cylinder, can be visualized in Fig. 12. The horizontal and the vertical force components are shown in Figs. 12(a) and (b), respectively, as a function of the non-dimensional time. Fig. 13 shows the same force components at point B. We see that the vertical velocity component (Fig. 11(b)) presents an oscillating behavior, corresponding to the vortices shedding process behind the cylinder. All force components (Figs. 11–13) presents the same behavior as the vertical velocity component, showing that the force calculation takes into account the neighboring ﬂow dynamic. The vortices shedding process determines the frequencies and the amplitude of the force oscillations. The acceleration force ~ fa ¼ ðfax ; fay Þ might tend to zero if the Lagrangian force had been calculated implicitly in such a way that at each time the iterative process could bring this term ð~ fa Þ to zero, before to advance to the next time. As the calculation were performed in a explicit way, it cannot be zero because the force used at time t þ Dt is calculated at the time t. It is very important to stress that at each time the acceleration force assumes a value diﬀerent from zero due to the fact that the ﬂuid velocity at the interface is recalculated at each time and, consequently, it assume a value diﬀerent from zero. All force components present the same oscillating behavior, with the same frequency, but with diﬀerent amplitude levels. The acceleration force assumes a very important role to compose the total force. The pressure force is the second more important, the viscous force is third more import and the inertial force is negligible as compared with the others terms. The no-slipping boundary condition must be guaranteed by this force ﬁeld. The ﬂuid velocity at the Lagrangian point the solid/ﬂuid interface is obtained by interpolating the nearest velocity ﬁeld. We qover ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ PN 2 2 deﬁne here l2 ¼ k¼1 ðuk þ vk Þ=N as a parameter that shows how good is the model in the sense that no slipping boundary condition must be assured. Fig. 14 shows the time evolution for the parameter l2 , which tends to zero as time pass. For this simulation l2 is of the order 104 . Fig. 15 shows the pressure ﬁeld for six Reynolds numbers and Fig. 16 shows the corresponding vorticity ﬁeld. It can be visualized in both ﬁgures that for Reynolds 20 and 40 the ﬂow is symmetric. For Reynolds 80 the instabilities take place and the ﬂow is completely asymmetric. The length of the bubble recirculation ðLw Þ is here deﬁned as the distance between two stagnation points downstream of the cylinder, as illustrated by Fig. 17. It was determined using the vertical velocity Fig. 14. Time evolution of the parameter l2 , for Re ¼ 10. A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 15. Pressure ﬁeld. (a) Re ¼ 20, (b) Re ¼ 40, (c) Re ¼ 80, (d) Re ¼ 100, (e) Re ¼ 150 and (f) Re ¼ 300. Fig. 16. Vorticity ﬁeld. (a) Re ¼ 20, (b) Re ¼ 40, (c) Re ¼ 80, (d) Re ¼ 100, (e) Re ¼ 150 and (f) Re ¼ 300. 365 366 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 distribution, plotted over the central streamline downstream of the cylinder. Fig. 18 shows these distributions as a function of y, for diﬀerent values of Reynolds number. The second point, over this streamline, where the velocity is zero gives the Lw measure. Fig. 19 shows the Lw distribution as a function of the Reynolds number. Comparisons of the present results with those of other authors [13,14] demonstrate good agreement. The drag coeﬃcient as a function of the dimensionless time for diﬀerent values of the Reynolds number is shown in Fig. 20. These results were obtained using Eq. (24). Fig. 21 shows the lift coeﬃcient as a function of time for diﬀerent Reynolds numbers. These results were obtained using Eq. (28). It should be noted that for Re ¼ 45 the ﬂow is still stable. As the Reynolds number increases the ﬂow becomes unstable. The critical Reynolds number is approximately Re ¼ 47, as predicted by the linear theory of stability. The h orientation angle, illustrated by Fig. 22, varies from 0 to 180°, between the two stagnation points. 2 Fig. 23 shows the pressure coeﬃcient on the cylinder surface, deﬁned as Cp ¼ pk =0:5qU1 , as a function of the angle h, as illustrated by Fig. 22. The pressure is taken on the Lagrangian points ~ xk . U1 is the free stream velocity. The comparison of the drag coeﬃcient, obtained in the present work, with other numerical and experimental results, is presented in Table 2. Very good agreement has been obtained [5,11,23,26]. Fig. 17. Length of the bubble recirculation. Fig. 18. Vertical velocity component ðvÞ as a function of y. Lw determination. A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 19. Bubble recirculation length ðLw =dÞ as a function of the Reynolds number. Fig. 20. Drag coeﬃcients as a function of the dimensionless time: (a) Re ¼ 10, 20, 40 and 50; (b) Re ¼ 80, 100, 150 and 300. Fig. 21. Variation of lift coeﬃcients with time: (a) Re ¼ 45 and 47; (b) Re ¼ 50 and 80. 367 368 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 Fig. 22. Schematic illustration of the angle h. Fig. 23. Pressure coeﬃcient distribution, between the stagnation points ðh ¼ 0Þ and ðh ¼ 180Þ: (a) Re ¼ 10, 20, 40 (present study); .: Re ¼ 10, 20, 40 ; (b) Re ¼ 50, 80, 100, 150 (present study);O: Re ¼ 80, 100 . Table 2 Comparison of mean drag coeﬃcient (CD ) with those of other authors Re Present work Park et al.  Sucker and Brauer  Dennis and Chang  Ye et al.  Tritton  10 20 40 47 50 80 100 150 300 2.81 2.04 1.54 1.46 1.46 1.40 1.39 1.37 1.27 2.78 2.01 1.51 – – 1.35 1.33 – 1.37 2.67 2.08 1.73 – 1.65 1.51 1.45 1.36 1.22 2.05 1.52 – – – – – – 2.03 1.52 – – 1.37 – – 1.38 2.22 1.48 – – 1.29 – – – A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370 369 Fig. 24. Strouhal number vs. Reynolds number. Fig. 24 shows the Strouhal number, deﬁned by Eq. (30), as compared with the Strouhal number obtained by other authors [16,25,26]. It should be noted that for low Reynolds numbers the results are in good agreement [9,14,24]. As the Reynolds number increases the error factor of the numerical simulation increases. The present work presents consistent results even for Re ¼ 300. 5. Conclusions The Physical Virtual Model (PVM), proposed in the present work, permits the simulation of unsteady, viscous incompressible ﬂows over an immersed circular cylinder using a Cartesian grid. Statistical parameters, including drag and lift coeﬃcients, the Strouhal number, wall pressure and the length of the bubble recirculation, were reported for several values of the Reynolds number. The results show consistency agreement with previous numerical and experimental results. The proposed model seems to be a promising tool and can be also applied to other geometries, to mobile geometries and to higher Reynolds numbers. The CPU time used to compute 10 s of simulation, for Re ¼ 80, was approximately 36 h in a Pentium IV, 2 GHz. This time was suﬃcient to obtain the statistic parameters and it is compatible with other methods such as Finite-Volume and the Lattice–Boltzman Methods. Acknowledgements ~o CAPES for the ﬁnancial support. The computations were perThe authors acknowledge Fundacßa formed on a Pentium IV 2.0 GHz, at the Mechanical Engineering College of the Federal University of Uberl^ andia. References  S. Armﬁeld, R. Street, The fractional-step method for the Navier–Stokes equations on staggered grids: the accuracy of three variations, J. Comp. Phys. 153 (1999) 660.  S.C.R. Dennis, G. Chang, Numerical solutions for steady ﬂow past a circular cylinder at Reynolds numbers up to 100, J. Fluid Mech. 42 (1970) 471.  E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yousof, Combined immersed-boundary ﬁnite-diﬀerence methods for threedimensional complex ﬂows simulations, J. Comp. Phys. 161 (2000) 35. 370 A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370  A.L. Fogelson, C.S. Peskin, A fast numerical method for solving three-dimensional Stokes equation in the presence of suspended particles, J. Comp. Phys. 79 (1988) 50.  B. Fornberg, A numerical study of steady viscous ﬂow past a circular cylinder, J. Fluid Mech. 98 (1980) 819.  D. Goldstein, R. Hadler, L. Sirovich, Modeling a no-slip ﬂow boundary with an external force ﬁeld, J. Comp. Phys. 105 (1993) 354.  D. Goldstein, R. Hadler, L. Sirovich, Direct numerical simulation of turbulent ﬂow over a modelled riblet covered surface, J. Comp. Phys. 302 (1995) 333.  D. Juric, Computation of phase change, Ph.D. Thesis, Mechanical Engineering, University of Michigan, USA, 1996.  J. Kim, D. Kim, H. Choi, An immersed-boundary ﬁnite-volume method for simulations of ﬂow in complex geometries, J. Comp. Phys. 171 (2001) 132.  M.C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comp. Phys. 160 (2000) 705.  R. Mittal, S. Balachandar, Inclusion of three-dimensional eﬀects in simulations of two-dimensional bluﬀ-body wake ﬂows, ASME Fluids Engineering Division Summer Meeting, 1997.  J. Mohd-Yusof, Combined immersed boundaries/B-splines methods for simulations of ﬂows in complex geometries, CTR Annual Reserch Briefs, NASA Ames/Stanford University, 1997.  M. Nishioka, H. Sato, Mechanism of determination of the shedding frequency of vortices behind a cylinder at low Reynolds numbers, J. Fluid Mech. 89 (1978) 49.  J. Park, K. Kwon, H. Choi, Numerical solutions of ﬂow past a circular cylinder at Reynolds number up to 160, KSME Int. J. 12 (1998) 1200.  C.S. Peskin, Numerical analysis of blood ﬂow in the heart, J. Comp. Phys. 25 (1977) 220.  A. Roshko, On the wake and drag of bluﬀ bodies, J. Aeronaut. Sci. 22 (1955) 124.  D. Sucker, H. Brauer, Fluiddynamik bei der angestromten Zilindern, W€arme Stoﬀubertragung 8 (1975) 149.  E.M. Saiki, S. Biringen, Numerical simulation of a cylinder in uniform ﬂow: application of a virtual boundary method, J. Comp. Phys. 123 (1996) 450.  G.E. Schneider, M. Zedan, A modiﬁed strongly implicit procedure for the numerical solution of ﬁeld problems, Numer. Heat Transfer 4 (1981) 1.  D.J. Triton, Experiments on the ﬂow past a circular cylinder at low Reynolds number, J. Fluid Mech. 6 (1959) 547.  Q. Wang, K.D. Squires, Int. J. Multiphase Flow 22 (1996).  F.M. White, Viscous Fluid Flow, second ed., McGraw-Hill, New York, 1991.  C. Wieselsberger, New data on the laws of ﬂuid resistance, NACA TN 84 (1992).  C.H. Williamson, Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds number, J. Fluid Mech. 206 (1989) 579.  C.H. Williamson, Vortex dynamics in the cylinder wake, Ann. Ver. Fluid Mech. 28 (1996) 477.  T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible ﬂows with complex boundaries, J. Comp. Phys. 156 (1999) 209.