# LimaSilva2003.pdf

by user

on
Category: Documents
1

views

Report

#### Transcript

LimaSilva2003.pdf
```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.
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).
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 [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 ﬂ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 [10] 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 [18] used this model to simulate a two-dimensional ﬂow around stationary and
moving cylinders.
Fogelson and Peskin [4] 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. [26] 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. [3] used discrete-time forcing, as suggested by Mohd-Yusof [12], 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. [9] 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 [15] and modiﬁed by Juric [8], 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 [1]. 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 [19]. 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 [18], 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 [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
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 [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 ﬂ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
[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 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 [14]; (b) Re ¼ 50, 80, 100, 150 (present study);O: Re ¼ 80, 100 [14].
Table 2
Comparison of mean drag coeﬃcient (CD ) with those of other authors
Re
Present work
Park et al.
[14]
Sucker and Brauer
[17]
Dennis and Chang
[2]
Ye et al.
[26]
Tritton
[20]
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
[1] 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.
[2] 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.
[3] 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
[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 ﬂow past a circular cylinder, J. Fluid Mech. 98 (1980) 819.
[6] D. Goldstein, R. Hadler, L. Sirovich, Modeling a no-slip ﬂow boundary with an external force ﬁeld, J. Comp. Phys. 105 (1993)
354.
[7] D. Goldstein, R. Hadler, L. Sirovich, Direct numerical simulation of turbulent ﬂow 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 ﬁnite-volume method for simulations of ﬂow 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 eﬀects in simulations of two-dimensional bluﬀ-body wake ﬂows, ASME
Fluids Engineering Division Summer Meeting, 1997.
[12] 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.
[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 ﬂow past a circular cylinder at Reynolds number up to 160, KSME Int. J. 12
(1998) 1200.
[15] C.S. Peskin, Numerical analysis of blood ﬂow in the heart, J. Comp. Phys. 25 (1977) 220.
[16] A. Roshko, On the wake and drag of bluﬀ bodies, J. Aeronaut. Sci. 22 (1955) 124.
[17] D. Sucker, H. Brauer, Fluiddynamik bei der angestromten Zilindern, W€arme Stoﬀubertragung 8 (1975) 149.
[18] 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.
[19] G.E. Schneider, M. Zedan, A modiﬁed strongly implicit procedure for the numerical solution of ﬁeld problems, Numer. Heat
Transfer 4 (1981) 1.
[20] D.J. Triton, Experiments on the ﬂow 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 ﬂuid 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 ﬂows with complex
boundaries, J. Comp. Phys. 156 (1999) 209.
```
Fly UP