by user

Category: Documents





Journal of Computational Physics 189 (2003) 351–370
Numerical simulation of two-dimensional flows over a
circular cylinder using the immersed boundary method
A.L.F. Lima E Silva a, A. Silveira-Neto
J.J.R. Damasceno
Mechanical Engineering College, Federal University of Uberl^
andia, Jo~
ao Naves de Avila Avenue,
Campus Santa M^
onica, Minas Gerais, Brazil
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
In this work, a virtual boundary method is applied to the numerical simulation of a uniform flow 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–fluid interface. These equations are discretized, using the finite differences method.
The immersed boundary is represented with a finite number of Lagrangian points, distributed over the solid–fluid
interface. A Cartesian grid is used to solve the fluid flow 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 flow 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 field, 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 coefficients and the Strouhal number,
calculated for an immersed cylinder, are compared with previous experimental and numerical results, for different
Reynolds numbers.
Ó 2003 Elsevier Science B.V. All rights reserved.
Keywords: Immersed-boundary method; Eulerian–Lagrangian formulation; Drag coefficient; Finite-difference method
1. Introduction
Immersed bodies present a very broad spectrum of practical problems due to their geometric complexity
and movement inside the flow. This class of problems is very difficult 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.
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
structure, which induces turbulence in the flow and turbulence induces vibration on the immersed structure.
This is the well-known fluid–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 [21].
Recently, some new types of methods based on the immersed boundary idea of Peskin [15], have been
proposed. The main idea is to use a regular Eulerian mesh for the fluid dynamics simulation, coupled with a
Lagrangian representation of the immersed boundary. The immersed boundary exerts a singular force on
the fluid. The Lagrangian mesh can move independent of the Eulerian mesh. Different models to calculate
this force field have been proposed recently.
Lai and Peskin [10] have proposed a second-order accurate immersed boundary method that has been
used to simulate flows over a circular cylinder. The interaction between the fluid 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 fixed to a
equilibrium point X e with a very stiff spring, whose stiffness constant is j 1.
Goldstein et al. [6,7] proposed a model called the Virtual Boundary Formulation, which was used to
simulate turbulent flow 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 fluid, applying a force field to
the fluid, 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 field 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 fluid velocity close to the interface velocity. The force intensity is calculated using
the expression
f ð~
xk ; tÞ ¼ a
Z th
~ ð~
~ ð~
xk ; tÞ dt þ b ~
xk ; tÞ U
xk ; tÞ ;
xk ; tÞ U
where a and b are constants that must be adjusted in order to obtain the expected physical behavior at the
flow. Saiki and Biringen [18] used this model to simulate a two-dimensional flow around stationary and
moving cylinders.
Fogelson and Peskin [4] solved the three-dimensional Stokes equations to simulate the sedimenting
particle flows 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 fluid–particle interface. The momentum equations are solved at all points
of the Eulerian grid. The particle points move with the local velocity fluid. This velocity is affected by
density force term that is calculated over the Lagrangian points and distributed over the neighboring
Eulerian points. The distributed force field has to model the non-linear and viscous effects between fluid
and particle surfaces. Knowing this Eulerian force field, the entire domain is viewed as one entity without
solid interfaces. The authors made simulations with one and two particles under the influence 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 flows with suspended
Ye et al. [26] proposed a method called the Cartesian Grid Method for simulating two-dimensional
unsteady, viscous, incompressible flows over complex geometries. In this method, a control volume near the
immersed boundary is re-formed into a body-fitted 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
Fadlum et al. [3] used discrete-time forcing, as suggested by Mohd-Yusof [12], and showed that it is more
efficient than the feedback forcing method for three-dimensional problems. The velocity at the first 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 first 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. [9] proposed a method based on a finite-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 field 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 flow and the flows around a circular cylinder, in order to validate the
methodology and the numerical procedure. Detailed informations about the flows over the cylinder, at
different Reynolds numbers are presented. These flow quantities are the mean drag and mean lift coefficients, the drag and lift fluctuations, 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 flows in a Cartesian square domain X containing an immersed boundary, as
illustrated in Fig. 1, can be modeled by the Navier–Stokes equations:
~P þ lr2 V~ þ ~
~ ¼ r
þ ðV~ r
~ V~ ¼ 0;
Fig. 1. Immersed boundary illustration; Eulerian mesh ð~
xÞ and Lagrangian mesh ð~
xk Þ.
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
where ~
F is given by
F ð~
xÞ ¼ ~
x ~
xk Þ d~
xk ;
f ð~
xk Þdð~
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
fluid flow, by injecting the force field on the fluid.
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 [15] and modified by Juric [8], represented by the following
Dij ð~
xk Þ ¼
f ½ðxk xi Þ=h
f ½ðyk yj Þ=h
< f1 ðrÞ
f ðrÞ ¼ 12 f1 ð2 krkÞ
if krk < 1;
if 1 < krk < 2;
if krk > 2;
f1 ðrÞ ¼
3 2krk þ
1 þ 4krk 4krk2
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 field, 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 Þ:
The different 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:
xk Þ;
xk Þ ¼ q
fa ð~
fi ð~
xk Þ ¼ qðV~ r
xk Þ;
xk Þ;
xk Þ ¼ lr2 V~ð~
fv ð~
~P ð~
xk Þ ¼ r
xk Þ:
fp ð~
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
Fig. 2. Scheme for the pressure and velocity interpolation.
The different terms described by Eqs. (9)–(12) must be evaluated over the interface using the velocity field
xÞ and pressure field P ð~
xÞ. These calculations must also take into account that, at the interface, the fluid
velocity must be equal to the interface velocity, which guarantees the no-slip condition. The velocity and
pressure spatial derivatives are calculated using the fluid 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 defined as being the region where the variables will be analyzed. In the
Fig. 3. Interpolation scheme for the horizontal velocity at point 3.
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 fields 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
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
different. 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 first and second derivatives in the x direction are approximated
ðxk x2 Þ
ðxk x1 Þ
ðxk x1 Þ þ ðxk x2 Þ
ðxk ; yk Þ ¼
/ þ
/ þ
ðx1 x2 Þðx1 xk Þ 1 ðx2 x1 Þðx2 xk Þ 2
ðxk x1 Þðxk x2 Þ k
o2 /
ðxk ; yk Þ ¼
ðx1 x2 Þðx1 xk Þ ðx2 x1 Þðx2 xk Þ ðxk x1 Þðxk x2 Þ
The derivatives in the y direction are given by:
ðyk y4 Þ
ðyk y3 Þ
ðyk y3 Þ þ ðyk y4 Þ
ðxk ; yk Þ ¼
/ þ
/ þ
ðy3 y4 Þðy3 yk Þ 3 ðy4 y3 Þðy4 yk Þ 4
ðyk y3 Þðyk y4 Þ k
o2 /
ðxk ; yk Þ ¼
oy 2
ðy3 y4 Þðy3 yk Þ ðy4 y3 Þðy4 yk Þ ðyk y3 Þðyk y4 Þ
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
flow properties.
The acceleration force qðoV~=otÞ, Eq. (9), is calculated taking into account that the fluid 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 fluid 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
Fij ¼
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.
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 finite central difference method in space and a
Euler first order in time. The velocity and pressure coupling was solved using a second-order pressure
correction method, as presented by Armfied and Street [1]. This method is described as follows:
1 opn
uinþ1 uni
N ðuni Þ þ Lðuni Þ þ Fin ;
q oxi
o2 unþ1
q o
Dt oxi
oxi oxi
uinþ1 ¼ uinþ1 Dt ounþ1
q oxi
pnþ1 ¼ pn þ unþ1 ;
Lðui Þ ¼
l o2 ui
q oxj oxj
N ðui Þ ¼
oðui uj Þ
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 field ~
f ð~
xk Þ, over the Lagrangian points ðxk ; yk Þ, using Eqs. (8)–(12) and the initial
2. Distribute the force ~
f ð~
xk Þ to the Eulerian grid using Eq. (17).
3. Estimate the fluid velocity field under the influence of the force field ~
F using Eq. (18).
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
4. Solve the linear system for the pressure correction given by Eq. (19).
5. Compute the actual pressure field using Eq. (21).
6. Update the fluid velocity using Eq. (20).
7. Verify the divergence using Eq. (3).
This procedure completes the time step loop. Normally, one loop is sufficient to obtain the mass conservation, verified 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.
The resulting linear system for the pressure correction u is solved using the interactive solver MSI
(Modified Strongly Implicit Procedure) of Schneider and Zedan [19]. Note that the interface force field
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 flow and the flow around a circular cylinder in order to validate the methodology and the
numerical procedure. For the two-dimensional Poiseuille flow simulation, the channel walls where modeled
using the PVM model. Simulations of the flow past a circular cylinder were done for different low Reynolds
numbers and consequently, the drag and lift coefficients, the length of the recirculation bubble and the
pressure coefficient on the cylinder surface were obtained. These results were compared with results present
by other authors.
4.1. Poiseuille flow
Fig. 7 shows the steady state of a vertical velocity profile 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 field was calculated at the channel walls. This
Fig. 7. (a) Force field at the channels walls. (b) Detail of the force at the wall. (c) Vertical velocity profile at y ¼ 2. The continuous line
represents the analytical result and the crosses the numerical result.
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 flow over a stationary cylinder. The boundary conditions
were imposed in such a way that the flow 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 effects on the flow development. A Newman boundary condition was used on the lateral
boundaries. A constant velocity profile U1 was specified at the domain entrance.
A grid refinement study was done to verify the result independence and the accuracy of the method. The
results are shown in Table 1 for five 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 flow 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 refinement study
Re ¼ 80
Dx ¼ Dy
Grid points
100 200
125 250
150 300
200 400
250 500
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
The velocity was imposed in such a way that the desired Reynolds number was obtained. The fluid
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 [18], who used a refined Lagrangian mesh in their simulations.
The Reynolds number in this flow was defined as
qU1 d
Re ¼
where d is the cylinder diameter. Finally, the dimensionless time was defined as
T ¼
U1 t
Once the velocity and pressure fields are calculated, the drag and lift coefficients and the Strouhal
number can also be calculated, using the force field directly.
Drag coefficient. As demonstrated by Lai and Peskin [10], 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
coefficient can be defined as
CD ¼
2 d
where FD is the drag force. Here, it can be calculated using the Eulerian force or the Lagrangian force as:
FD ¼ Fx dx ¼ fx ds;
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 coefficient. Similar to the drag coefficient, the lift coefficient can be defined as
CL ¼
2 d
where FL is the lift force, which can be calculated using the y component of the Eulerian force or the
Lagrangian force as
FL ¼ Fy dy ¼ fy ds;
where Fy and fy are the y components of the Eulerian and Lagrangian forces, respectively.
Strouhal number. The Strouhal number is defined as the dimensionless frequency with which the vortices
are shed behind the body
St ¼
where f is the vortex shedding frequency. This frequency can be obtained using the Fast Fourier Transform
of the lift coefficient time distribution.
A small time step ðDt ¼ 1 106 Þ was used at the first loop, in order to guarantee the code stability. It is
increased as the flow develops. When the force field is enough to satisfy the no-slip condition on the cylinder
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 [10]). 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 flow 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 different from zero also inside the interface. In fact, we have a velocity field 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 field illustration for Re ¼ 40: (a) Fx , (b) Fy .
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
not affect the Lagrangian force calculation because only the external fields are used during the interpolation
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. Specific 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).
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 fluid 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 flow 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 different from zero due to the fact that the fluid velocity at the interface is
recalculated at each time and, consequently, it assume a value different from zero.
All force components present the same oscillating behavior, with the same frequency, but with different
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 field. The fluid velocity at the
Lagrangian point
the solid/fluid interface is obtained by interpolating the nearest velocity field. We
define 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 field for six
Reynolds numbers and Fig. 16 shows the corresponding vorticity field. It can be visualized in both figures
that for Reynolds 20 and 40 the flow is symmetric. For Reynolds 80 the instabilities take place and the flow
is completely asymmetric.
The length of the bubble recirculation ðLw Þ is here defined 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 field. (a) Re ¼ 20, (b) Re ¼ 40, (c) Re ¼ 80, (d) Re ¼ 100, (e) Re ¼ 150 and (f) Re ¼ 300.
Fig. 16. Vorticity field. (a) Re ¼ 20, (b) Re ¼ 40, (c) Re ¼ 80, (d) Re ¼ 100, (e) Re ¼ 150 and (f) Re ¼ 300.
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 different 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
The drag coefficient as a function of the dimensionless time for different values of the Reynolds number
is shown in Fig. 20. These results were obtained using Eq. (24).
Fig. 21 shows the lift coefficient as a function of time for different Reynolds numbers. These results were
obtained using Eq. (28). It should be noted that for Re ¼ 45 the flow is still stable. As the Reynolds number
increases the flow 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.
Fig. 23 shows the pressure coefficient on the cylinder surface, defined as Cp ¼ pk =0:5qU1
[22], 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 coefficient, 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 coefficients 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 coefficients with time: (a) Re ¼ 45 and 47; (b) Re ¼ 50 and 80.
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 coefficient distribution, between the stagnation points ðh ¼ 0Þ and ðh ¼ 180Þ: (a) Re ¼ 10, 20, 40 (present study);
.: Re ¼ 10, 20, 40 [14]; (b) Re ¼ 50, 80, 100, 150 (present study);O: Re ¼ 80, 100 [14].
Table 2
Comparison of mean drag coefficient (CD ) with those of other authors
Present work
Park et al.
Sucker and Brauer
Dennis and Chang
Ye et al.
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
Fig. 24. Strouhal number vs. Reynolds number.
Fig. 24 shows the Strouhal number, defined 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 flows over an immersed circular cylinder using a Cartesian grid. Statistical parameters, including drag and lift coefficients, 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 sufficient to obtain the statistic parameters and it is compatible with
other methods such as Finite-Volume and the Lattice–Boltzman Methods.
~o CAPES for the financial 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
[1] S. Armfield, 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.
[2] S.C.R. Dennis, G. Chang, Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100, J. Fluid
Mech. 42 (1970) 471.
[3] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yousof, Combined immersed-boundary finite-difference methods for threedimensional complex flows simulations, J. Comp. Phys. 161 (2000) 35.
A.L.F. Lima E Silva et al. / Journal of Computational Physics 189 (2003) 351–370
[4] 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.
[5] B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, J. Fluid Mech. 98 (1980) 819.
[6] D. Goldstein, R. Hadler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comp. Phys. 105 (1993)
[7] D. Goldstein, R. Hadler, L. Sirovich, Direct numerical simulation of turbulent flow over a modelled riblet covered surface,
J. Comp. Phys. 302 (1995) 333.
[8] D. Juric, Computation of phase change, Ph.D. Thesis, Mechanical Engineering, University of Michigan, USA, 1996.
[9] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, J. Comp.
Phys. 171 (2001) 132.
[10] 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.
[11] R. Mittal, S. Balachandar, Inclusion of three-dimensional effects in simulations of two-dimensional bluff-body wake flows, ASME
Fluids Engineering Division Summer Meeting, 1997.
[12] J. Mohd-Yusof, Combined immersed boundaries/B-splines methods for simulations of flows in complex geometries, CTR Annual
Reserch Briefs, NASA Ames/Stanford University, 1997.
[13] 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.
[14] J. Park, K. Kwon, H. Choi, Numerical solutions of flow past a circular cylinder at Reynolds number up to 160, KSME Int. J. 12
(1998) 1200.
[15] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comp. Phys. 25 (1977) 220.
[16] A. Roshko, On the wake and drag of bluff bodies, J. Aeronaut. Sci. 22 (1955) 124.
[17] D. Sucker, H. Brauer, Fluiddynamik bei der angestromten Zilindern, W€arme Stoffubertragung 8 (1975) 149.
[18] E.M. Saiki, S. Biringen, Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method, J. Comp.
Phys. 123 (1996) 450.
[19] G.E. Schneider, M. Zedan, A modified strongly implicit procedure for the numerical solution of field problems, Numer. Heat
Transfer 4 (1981) 1.
[20] D.J. Triton, Experiments on the flow past a circular cylinder at low Reynolds number, J. Fluid Mech. 6 (1959) 547.
[21] Q. Wang, K.D. Squires, Int. J. Multiphase Flow 22 (1996).
[22] F.M. White, Viscous Fluid Flow, second ed., McGraw-Hill, New York, 1991.
[23] C. Wieselsberger, New data on the laws of fluid resistance, NACA TN 84 (1992).
[24] 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.
[25] C.H. Williamson, Vortex dynamics in the cylinder wake, Ann. Ver. Fluid Mech. 28 (1996) 477.
[26] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex
boundaries, J. Comp. Phys. 156 (1999) 209.
Fly UP