Proceedings of the Tenth (2000) International Offshore and Polar Engineering Conference

Seattle, USA, May 28-June 2, 2000

Copyright 2000 by The International Society of Offshore and Polar Engineers
ISBN 1-880653-46-X (Set); ISBN 1-880653-49-4 (Vol. III); ISSN 1098-6189 (Set)

Simulating Vortex Shedding at High Reynolds Numbers

P.A.B. de Sampaio
Instituto de E n g e n h a r i a N u c l e a r / C N E N , Rio de Janeiro, Brazil
A.L. G.A. Coutinho
COPPE, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil

ABSTRACT
But the sub-grid model does not need to be explicit: it can be
inherent to the numerical scheme. Indeed, Hughes (1995) has shown
the relationship between stabilized finite element formulations and the
sub-grid modeling of the unresolvable scales.

The turbulent vortex shedding arising from the cross flow past a
circular cylinder is analyzed using a Large Eddy Simulation (LES)
procedure. In particular, no explicit sub-grid stress model is
employed. Rather, the unresolved scales are dealt implicitly by a
stabilized Petrov-Galerkin finite element formulation used in
conjunction with a time-space adaptive scheme. The numerical
results are compared with available experimental data on force
coefficients and vortex shedding frequency.

Breuer (1998) has proposed the name LES without sub-grid model to
what we would tentatively call LES with a numerically implicit subgrid model. What is important to bear in mind is the distinction
between such an approach and a Direct Numerical Simulation (DNS) of
turbulence, where the time and space discretizations are fine enough to
resolve all turbulence scales. Furthermore, note that the design of a LES
with a numerically implicit sub-grid model becomes the design of the
numerical method itself. This includes not only the formulation used to
obtain the discretized equations, but also the adaptive schemes and
other algorithms that affect the way the unresolvable scales are treated
(implicitly modeled) by the computation.

KEY WORDS: vortex shedding, large eddy simulation, stabilized finite

elements, adaptive methods, circular cylinder.
INTRODUCTION
Vortex induced vibration is a major concern in Offshore
Engineering. The flow around marine structures causes vibrations that
may lead to failure by fatigue. The problem is of particular importance
in deep water oil exploitation systems, where risers and mooring lines
can be viewed as flexible slender cylinders subjected to cross flow.
Moreover, for deep water systems the variation of water currents from
the ocean surface to the seabed may be large. There is some
controversy among researchers on what is the most appropriate way to
tackle the problem (Franciss, 1999).

In this work, instead of adopting an explicit sub-grid stress model,

the effect of the unresolved scales is accounted for implicitly through
the use of a stabilized Petrov-Galerkin formulation. The simulation of
the cross flow past a circular cylinder is performed using the adaptive
parallel/vector program developed by De Sampaio and Coutinho (1999)
for the analysis of incompressible viscous flow. The discretized
equations are obtained using a finite element Petrov-Galerkin
formulation that automatically introduces streamline upwinding and
allows equal order interpolation for velocity and pressure. A remeshing
scheme, guided by the estimated error on the viscous stresses, is
combined with a local time-stepping procedure, leading to a time-space
adaptive computation of the turbulent vortex shedding problem.

We think that the development of Computational Fluid Dynamics

(CFD) will have an increasingly important role in the analysis of the
vortex-induced vibrations of these slender structures subjected to depthvarying currents. Although the coupled transient three-dimensional
simulation of flow and structure needed to address the problem directly
is still a colossal task, the time is ripe for paving the way towards such
an ambitious goal. In this paper we tackle the more basic, but
challenging problem, of simulating the turbulent vortex shedding
around a stationary cylinder at the high Reynolds numbers typical of
Offshore Engineering applications.

Turbulent analyses with Reynolds number in the range from 104 t o

106 are presented. The quality of these analyses is assessed by
comparing the numerical Drag and Lift coefficients, and the numerical
Strouhal number, with available experimental data.

A Large Eddy Simulation (LES) procedure is used to study the

turbulent flow around a circular cylinder. In a Large Eddy Simulation
the large turbulence scales are resolved by the discretization while the
small scales are taken into account through the so-called sub-grid
models.

THE CONTINUUM AND DISCRETE MODELS

We present here the finite element formulation used for the
simulation of incompressible viscous flow. The problem is defined on

461

whilst the pressure-continuity equation is

the open bounded domain ~ , with boundary F, contained in the nsddimensional Euclidean space. The governing equations are the
incompressible Navier-Stokes equations. These equations are written
using the summation convention for a =I . . . . . nsd and b=1 ..... nsd as

(:,o+

PC +t

:uo~ : r (:'a+:U, lt+:P =0

ub~')-'Ox~-LI't~x b ~ G ) J

At 8N~ 8 ~ "+1/2

! :xo :x---T- m : -rat+N'&:xo pa; :~;+x~dnI N,

(1)

+ ~ S A t :IVi +

(2)

d u "-=O
cgx a

d,-- I

r~

+x.

,o,

:xo :x,

8Xb

(4)

[//(:u:+Dtx:l]d.Q
t.:x,

In the above equations ft,

where p and /2 are the fluid density and viscosity, respectively. The
velocity and pressure dependent variables are denoted by u~ and p,
respectively.

axoj]

and ~ are the discretized fields,

interpolated using standard C o finite elements whose shape functions

are denoted by N~. The superscripts n, n+1/2 and n+l indicate the
time-level. Note that the terms explicitly written as summation of
element integrals vanish altogether for the linear shape functions used.

The model is completed introducing boundary conditions and an

initial velocity field. Velocity and traction boundary conditions are
given data g, and ?a They are prescribed on the boundary partitions

The problem is solved using a segregated solution procedure.

Pressure is computed first, then the velocity field is updated. Eqs. (3)(4) lead to symmetric positive definite matrices, allowing the use of
preconditioned conjugate-gradient solvers.

F.~ and F,~, such that F.ouF,~=F and r ~ . n F , . = ~ . Pressure and

mass flux boundary conditions are associated to the mass balance and
are given as ff and G on the boundary partitions Fp and F c , such

Discretizations similar to that presented here have been obtained by

(Zienkiewicz and Wu, 1991) following alternative approaches. The
procedure can be extended to the analysis of natural convection
problems with the inclusion of the energy balance and accounting for
buoyancy forces in the equation of momentum (De Sampaio and
Coutinho, 1999).

that F p ~ F o = F and Fpr'~Fo=O .

The continuum model is discretized using linear triangular elements
to approximate velocity and pressure. Such a choice of interpolating
spaces is not acceptable within the mixed-formulation framework, as it
violates the Babu~ka-Brezzi condition (Brezzi and Fortin, 1991).
However, the Petrov-Galerkin formulation we shall present next avoids
such a difficulty through the introduction of extra stabilizing terms
(Hughes, Franca and Ballestra, 1986), (De Sampaio, 1991). The
formulation also leads to better approximations of convection
dominated flows, for it generates streamline upwinding (Brooks and
Hughes, 1982), (De Sampaio, 1991). The reader is referred to (De
Sampaio and Coutinho, 1999) for a detailed description of the method
used in this work.

LES WITH A TIME-SPACE ADAPTIVE PROCEDURE

In this work we combine a remeshing scheme with the local timestepping algorithm introduced by De Sampaio (1993). This algorithm
sets local time-steps based on the time-scales of the convection
diffusion processes resolvable on a given mesh. These time-scales are
estimated according to local values of velocity and physical properties,
and according to the local mesh resolution. The local time-steps are
chosen according to the expression,

The following discretized model is obtained by least-squares

minimization of the time-discretized momentum balance with respect to
the velocity and pressure degrees of freedom. The pressure-continuity
equation Eq. (4) results from the combination of the least-squares
momentum minimization and the requirement of flow incompressibility
Eq.(2). The discretized momentum balance is given by

he

(5)

a,_- = Ilu"ll
where

Air

+-'2-<-~'7~

+ - - U b - -c~x
l abx z) 2

(URe) = Fcoth{URe l--Z-I

L

D.+,,Ua _fat
DN

- f:,

.(

^"

) dn

& ax,-t.ex, :x: )

E,a,..au,
LT" :x
+ Iu,

based on the local velocity modulus and on the local element size he

(3)

^"

+ F

(6)

NReJ

Note that NRe-- :llu"tlh,//2 is the local element Reynolds number

--dfI2

O Xa

t. 2 )

There are some important reasons for such a choice. First, the
weighting function used to approximate the momentum balance
becomes the SUPG weighting function (Brooks and Hughes, 1982) and
a correct amount of streamline upwinding is introduced in the
formulation. Further, note that the time-step given by Eq.(5) is
appropriate to follow the time evolution of the convection-diffusion
processes resolvable on a mesh with local size h, (De Sampaio, 1991).

:;,: +~:'~l d~

o.r

Indeed, Eq.(5)gives for pure convection A t=hJliu"ll

II
II

462

whereas for

The local time-stepping algorithm is used in conjunction with the

remeshing scheme. This permits linking spatial and time-step
refinement and naturally leads to a simultaneous time-space adaptive
procedure. Indeed, whenever the remeshing scheme creates some local
mesh refinement to better resolve a particular flow feature, the timestep distribution is also adapted accordingly, so that the corresponding
time evolution can be appropriately followed.

pure diffusion A t = p h ~ / 6 / z , Finally, Hughes (1995) has shown that

there is a relationship between the time-step given by Eq.(5), also
called the intrinsic time scale, and the modeling of the sub-grid (or
unresolvable) scales in the context of stabilized formulations.
In this work, instead of using an explicit sub-grid model, the
modeling of the unresolvable scales is inherent to the stabilized
formulation adopted and to the choice of Eq.(5) to evaluate the local
time-steps.

It is the combination of the stabilized finite element formulation, the

local time-stepping algorithm, and the remeshing scheme, that defines
the particular LES with a numerically implicit sub-grid model used in
our computations.

The role played by the stabilization terms as filters of turbulent subgrid phenomena deserves further comments. Note that if a Galerkin
method is used, the computation becomes unstable due to the negative
diffusion engendered by the Galerkin approximation of strong
convection in coarse discretizations (note that the meaning of strong
convection and coarse discretization here is relative: the Galerkin
method is adequate if the local refinement and physical data is such that
NRe<< 1 ).

CIRCULAR CYLINDER IN CROSS FLOW

The numerical techniques described above have been used to
simulate the vortex shedding past a circular cylinder in cross flow. The
problem has been non-dimensionalized using the cylinder diameter d
and the free stream velocity u 0 as reference scales for length and
velocity, respectively. The reference scale chosen for pressure was
PUo2 and time has been non-dimensionalized by d / u o . The

In this work, though, we are using a stabilized formulation, with

local time-steps defined according to Eq.(5). Thus, comparing with the
Galerkin method, we are introducing back into the approximation the
amount of diffusion required to balance the negative diffusion arising
from unresolvable scales (Kelly, Nakazawa and Zienkiewicz, 1980).
The balancing diffusion term appears from the interaction between the
non-Galerkin part of the weighting (also called the perturbation) and
the convective term in Eq.(3), i.e.

simulations are parameterized by the global Reynolds number

Re = t9 u o d/flt.
Figure 1 depicts the analysis domain. Impermeability and non-slip
velocity boundary conditions are prescribed on the cylinder surface. A
uniform velocity field with u~ = u 0 and u 2 = 0 is imposed on the face
AB. At faces AE and BD the boundary condition imposed is u2 = 0.
Free traction and a reference pressure p = 0 are the boundary
conditions prescribed at face ED. The initial condition is the velocity
field with components u~ = Uo and u 2 = 0, which is specified over all

The above term acts as a filter that damps the growth of the
numerical errors associated to sub-grid phenomena, i.e. flow features
smaller/faster than the local h, and At. In particular, note that the

analysis domain at time t = 0.

Simulations of the cross flow past a circular cylinder have been
performed for R e = 104, 5.104, 105, 5.105 and 106. The numerical

balancing diffusion term depends on the local space-time resolution,

i.e., h, and At=ct(NRe)he/llul[.

results are assessed in terms of the drag and lift forces acting on the
cylinder and of the vortex shedding frequency. The drag and lift
coefficients are respectively given by C o = 2 F o / p u ~ d
and

Of course, the discretization required to suitably model a particular

turbulent flow will be always problem dependent. Furthermore, special
boundary conditions, incorporating the universal law of the wall, may
be required to model turbulent boundary layers. Nonetheless, the
method proposed herein is capable of following the evolution of the
large/slow features resolvable in a given discretization. In this sense, it
is a LES method where no explicit sub-grid modeling is introduced.

C L = 2 F L / p u ~ d , where F o and F L are the drag and lift forces on

the cylinder (per unit span). The vortex shedding frequencyfis given in
non-dimensionalform by the Strouhal number St = f d / u o .
B

It is important to note that the time-step given by Eq.(5) varies

spatially according to local values of velocity, physical properties and
mesh size. Therefore, we need to consider a spatially varying time-step
distribution when advancing the computation in time. We employ here
an algorithm which allows each degree of freedom to advance in time
according to its own local time-step, whilst interpolated results are
periodically output at fixed times. The reader is referred to (De
Sampaio, 1993) and to (De Sampaio and Coutinho, 1999) for details
regarding the local time-stepping algorithm.

uo

The a posteriori error estimator proposed by Zienkiewicz and Zhu

(1987) is used to estimate the error on the viscous stresses and to guide
the remeshing. The scheme is designed to generate meshes containing a
controlled number of elements, in such a way that that the error on the
viscous stresses becomes evenly distributed (De Sampaio and
Coutinho, 1999). The remeshing procedure is fully automatic and
triggered during a transient analysis, whenever the relative variation of
the estimated error exceeds a preset value of 0.01.

Figure 1. Analysis domain used in the simulation of cross flow past a

circular cylinder: the length of the segments are
AM=MB=DN=NE=AF=BC=5 d and CD=FE= 10 d.

463

Everything changes with a slight increase of the Reynolds number,

beyond the critical value of 3.10L When this situation is reached, the
boundary layer undergoes a turbulent transition and reattatches to the
cylinder surface. Separation is retarded to a position behind the
cylinder, at about 120 degrees. The comparatively smaller area of low
pressure on the back of the cylinder explains the experimentally
observed abrupt drag reduction associated to the sudden change of the
point of separation.

The computations were performed on the Cray J90 machine

installed at COPPE/UFRJ. The numerical results obtained are presented
next. Figure 2 shows the time evolution of the drag and lift coefficients
during the simulated transients for Re= 5.104,Re=105 . The time
evolution of the number of elements used in these adaptive
computations is shown in Figure 3 . Figure 4 depicts the velocity field
and the adaptive mesh close to the cylinder for the Re = l0 S analysis,
at time t = 80, long after the periodic behavior has been established.
Note the adaptive refinement close to the cylinder.

This phenomenon, called as the drag crisis, could not be captured

by our computations based on coarse near wall meshes, unsuitable to
model the turbulent boundary layer and its transition. It is clear that
more refinement close to the wall, or the use of special wall boundary
conditions, is required in order to detect the drag crisis. We also
performed similar computations using the Smagorinsky model in
meshes containing about 10,000 elements, but still insufficiently
refined in the near wall region. The results were no better. Selvam
(1997) observed a moderate reduction of the drag force in his 2D LES
analyses, based on the Smagorinsky sub-grid model, using the law of
the wall boundary conditions. However, it seems that an accurate
simulation of drag crisis can only be accomplished using a full 3D
model (Selvam, 1997).

DISCUSSION OF RESULTS
The rms drag coefficients and the Strouhal numbers are presented in
Table 1. Despite the use of very coarse meshes in the near wall region
(the minimum element size is 0.01d), the results obtained in the range
from 104 to 105 are accurate enough for most engineering purposes. In
this range, the drag coefficient is overestimated by about 30% with
respect to the experimental data reported by Roshko (1961), whilst the
Strouhal number is higher than the experimental values by about 10%.

Re

Numerical

Experimental

CONCLUDING REMARKS

Cd

St

Cd

St

10,000

1.2

0.20

1.635

0.21

50,000

1.2

0.19

1.546

0.22

100,000

1.2

0.19

1.542

0.21

500,000

0.30-0.35

0.23

1.538

0.21

! ,000,000

0.4

0.32

1.540

0.22

A Large Eddy Simulation procedure, combining a stabilized PetrovGalerkin formulation and adaptive time-space discretizations, was
employed to compute the unsteady flow around a fixed circular
cylinder. The problem was simulated at the high Reynolds numbers
typical of Offshore Engineering applications, i.e. 104 to 106.
In the present method, the sub-grid modeling of the unresolvable
scales is inherent to the stabilized Petrov-Galerkin formulation used
and to the choice of suitable local time-steps. Thus, no explicit sub-grid
modeling is required. The LES procedure presented here naturally turns
into a Direct Numerical Simulation (DNS) method, provided a
sufficiently fine discretization is adopted.

Table 1. Computed drag coefficient (rms) and Stmuhal number.

But how can such coarse discretizations in the near wall region yield
reasonably accurate predictions of the force acting on the cylinder?
After all, the refinement used is insufficient to resolve the turbulent
boundary layer. Moreover, no special boundary conditions introducing
the universal law of the wall have been used.

The results obtained are accurate enough for most engineering

purposes, up to the critical Reynolds number Re=3.105, which
characterizes the beginning of the so-called drag crisis. On the other
hand, the accurate simulation of the drag crisis using the present
method will require finer meshes, especially in the near wall region,
and a full 3D model.

In order to answer that we need to consider the experimental

observations regarding the flow patterns that affect the force on the
cylinder (Panton, 1984). In particular, we must consider the points of
separation of the boundary layer. These points roughly divide the
pressure acting on the cylinder surface into high pressure, before the
points of separation, and low pressure after those points. As long as the
main contribution to the force on the cylinder is due to pressure, this
force will strongly depend on the position of the boundary layer
separation. It happens that for the sub-critical range (Reynolds numbers
up to 3.105), separation occurs in the front part of the cylinder, at about
80 degrees along the cylinder surface. This explains why the force
coefficients observed experimentally are nearly independent of
Reynolds, although the boundary layers are getting ever more thinner
with the increase of the Reynolds number.

The results obtained so far are promising, and we expect that the
CFD approach proposed in this work can be further developed to
address some important fluid-structure vibration problems of Offshore
Engineering.
ACKNOWLEDGMENTS
The authors kindly acknowledge the Center for Parallel
Computation of COPPE/UFRJ for the use of the Cray J90 computer.
REFERENCES
Breuer, M (1998). "Large Eddy Simulation of the Subcritical Flow Past
a Circular Cylinder: Numerical and Modelling Aspects", Int. J.
Numer. Meth. Fluids, Vol. 28, pp. 1281-1302.
Brezzi, F ,and Fortin M (1991). Mixed and Hybrid Finite Element
Methods, Springer, New York.
Brooks, A, and Hughes, TJR (1982). "Streamline Upwind/PetrovGalerkin Formulation for Convection Dominated Flows with
Particular Emphasis on the Incompressible Navier-Stokes

We have obtained reasonably accurate force coefficients, in the subcritical range, because our simulations predict the point of separation
reasonably well, in spite of being unable to resolve the turbulent
boundary layer. Of course, if we were solving a heat transfer problem,
the result would be useless. This serves to remind that models cannot
be judged good or bad by themselves, but only in relation to what is
expected from them.

464

Equations", Comput. Methods Appl. Mech. Engrg., Vol. 32, pp.

199-259.
De Sampaio, PAB, and Coutinho, ALGA (1999). "Simulation of Free
and Forced Convection Incompressible Flows using an Adaptive
Parallel / Vector Finite Element Procedure", Int. J. Numer. Meth.
Fluids, Vol. 29, pp. 289-309.
De Sampaio, PAB (1991). "A Petrov-Galerkin Formulation for the
Incompressible Navier-Stokes Equations using Equal Order
Interpolation for Velocity and Pressure", Int. J. Numer. Meth.
Engrg., Vol. 31, pp. 1135-1149.
De Sampalo, PAB (1993). "Transient Solutions of the Incompressible
Navier-Stokes Equations in Primitive Variables Employing
Optimal Local Time-Stepping", Proc. 8th Int. Conference on
Numer. Meth. for Laminar and Turbulent Flow, Swansea, U.K.,
pp. 1493-1504.
Franciss, R (1999). "Vortex Induced Vibration of Slender Structures of
Offshore Vessels", D.Sc. Thesis, COPPFE/UFRJ, Rio de Janeiro,
Brazil (in Portuguese).
Hughes, TJR (1995). "Multiscale Phenomena: Green's Functions, the
Dirichlet-to-Neumann Formulation, Subgrid Scale Models,
Bubbles and the Origins of Stabilized Methods", Comput. Methods
Appl. Mech. Engrg., Vol. 127, pp. 387-401.
Hughes, TJR, Franca, LP, and Ballestra, M (1986). "A New Finite
Element Formulation for Computational Fluid Dynamics: V.
Circumventing the Babu~ka-Brezzi Condition: A Stable PetrovGalerkin Formulation of the Stokes Problem Accommodating
Equal-Order Interpolation", Comput. Methods Appl. Mech. Engrg.,
Vol. 59, pp. 85-99.
Kelly, DW, Nakazawa, S, and Zienkiewicz, OC (1980). "A Note on
Anisotropic Balancing Dissipation in the Finite Element Method
Approximation to Convective Diffusion Problems", Int. J. Numer.
Meth. Engrg., Vol. 15, pp. 1705-1711.
Panton, RL (1984). Incompressible Flow, John Wiley & Sons, New
York.
Roshko, A (1961). "Experiments on the Flow Past a Circular Cylinder
at Very High Reynolds Number", J. Fluid Mech., Vol. 10, pp. 345356.
Selvam, RP (1997). "Finite Element Modelling of Flow Around a
Circular Cylinder using LES", J. Wind Eng.lnd.Aerdyn, Vol. 67 &
68, pp. 129-139.
Zienkiewicz, OC, and Wu, J (1991). "Incompressibility without Tears How to Avoid Restrictions of Mixed Formulation", Int. J. Numer.
Meth. Engrg., Vol. 32, pp. 1189-1203.
Zienkiewicz, OC, and Zhu, JZ (1987). "A Simple Error Estimator and
Adaptive Procedure for Practical Engineering Analysis", Int. J.
Numer. Meth. Engrg., Vol. 24, pp. 337-357.

CD-CL ......

,
,f' il

,li,'~ l I I !I

I'l
, I Ii t
| 'l, "E : I l,'
,,,,, I ',ii ,I
~J I
J I

-2

l0

20

30

I,

40 50
Time

,, I , I

60

, 1

70

80

Re=50,000

CD
CL ......
2

,
0

-1

-2
10

20

30

40 50
Time

60

70

80

Re=100,000

Figure 2. Drag (CD) and lift (CL) coefficients for Re=50,000 and
100,000.

465

Elements m

5500

5(XlO

4500

4000

10

20

30

40
50
'rime

60

70

80

Re=50,000

6000

Elements

5500

5000

4500

4000

10

20

30

40 50
Time

60

70

80

Re= 100,000

Figure 4. Detail of the adaptive m e s h and velocity field for Re=100,000

at time t=80.

Figure 3. Evolution o f total n u m b e r of elements in the analysis for

R e - 5 0 , 0 0 0 and 100,000.

466