J.
Fluid
Mech.
(1980),
vol.
98,
part
4,
pp.
819-855
Printed in Great Britain
819
A
numerical study
of
steady viscous
flow
past a circular cylinder
By
BENGT
FORNBERG
Department of Applied Mathematics, California Institute of Technology, Pasadena,
California 91125,
U.S.A.
(Received
13
June 1979)
Numerical solutions have been obtained for steady viscous flow past
a
circular c..-l;nder
at Reynolds numbers up to
300.
A
new technique is proposed for the boundwy con-
dition at large distances and an iteration scheme has been developed, based on New-
ton’s method, which circumvents the numerical difficulties previously encountered
around and beyond a Reynolds number of
100.
Some new trends are observed in the
solution shortly before a Reynolds number of
300.
As
vorticity starts
to
recirculate
back from the end of the wake region, this region becomes wider and shorter. Other
flow quantities like position
of
separation point, drag, pressure and vorticity distribu-
tions on the body surface appear to be quite unaffected by this reversal of trends.
1.
Introduction
The problem of viscous incompressible flow past a circular cylinder has for a long
time received much attention, both theoretically and numerically. In spite of many
numerical methods and calculations, the Reynolds number
Re
=
100
(based on the
diameter) appears to be the upper limit for which complete, steady flow fields have been
reliably determined. There are many reasons
for
the continuing interest in this prob-
lem and in attempts to carry numerical calculations to still higher Reynolds numbers.
One of these reasons
is
that
it
is
a good model problem for flows past other bodies of
practical importance. Steady solutions for the circular cylinder become experimen-
tally unstable around
Re
=
40.
Use of flow control methods to stabilize unstable
solutions could lead to important new classes
of
flows, which at first may be studied
more easily numerically. Many difficulties are encountered in attempts to analytically
describe the complete flow field. We believe that numerical methods can provide
further information on the limiting properties
of
the steady flow for increasing Rey-
nolds numbers. Questions such as the asymptotic development of the recirculation
region (wake bubble), drag, position of separation point, vorticity and pressure distri-
butions, etc., are all open and they are relevant to the understanding of high-Reynolds-
number flows.
Brodetsky
(1923)
suggests
a
solution
for
infinite Reynolds number in which vortex
sheets bound an infinite wake region containing stagnant flow. This solution is often
referred to as the Helmholtz-Kirchhoff free streamline model because
of
their intro-
duction
of
vortex sheets (Helmholtz
1868;
Kirchhoff
1869).
Both an infinite wake and
a finite drag
is
in agreement with experimentally and numerically observed trends for
low Reynolds numbers (up to
100
to
200).
Batchelor
(1956)
gives however arguments
against these features of the limit and suggests that the vorticity inside the wake in the
0022-1 120/80/4471-1380 $02.00
0
1980
Cambridge
University
Press
820
B.
Fornberg
limit need not be zero but can take a constant value on each side of the line of sym-
metry. Solutions of this kind (if they exist) would allow for a finite wake and therefore
no drag on the body.
The most notable phenomenon we observe at high Reynolds numbers
(Re
>
260)
is
a shortening of the wake region with vorticity being convected into
its
interior. Our
upper limit,
Re
=
300,
is
not high enough however to establish whether
or
not this new
trend will persist.
In
a
recent theoretical work, Smith
(1979)
assumes that the wake length will increase
to infinity and that the flow will tend to the Brodetsky limit with no vorticity inside
the wake. There are of course profound differences between the consequences he draws
from these assumptions and our results for
Re
>
260.
However, there are also differ-
ences for much lower Reynolds numbers as for example his figure
3
for the skin friction
shows.
Our
numerical results agree well with the comparisons he quotes and confirm
therefore the discrepency illustrated in that figure.
The flow problem is formulated mathematically in
$2.
In
$
3
we discuss the most
frequently encountered numerical difficulties. Section
4
discusses boundary conditions
at large distances. The final numerical method is presented in
$5
with the obtained
results discussed in
9
6.
This numerical method was employed for Reynolds numbers
20-300.
Some results for Reynolds numbers
2-10
are also presented, although they
were obtained by a different technique, using a direct iteration scheme based on fast
Poisson solvers. Since many methods work successfully in that range and the flow
patterns are well known, this low-Reynolds-number method will not be discussed
further.
The extensive numerical calculations for high Reynolds numbers were performed on
the Control Data Corporation STAR
100
computer, located at the CDC Service Center
in Arden Hills, Minnesota. We wish to express our gratitude to Control Data Corpora-
tion for making this computer system available to
us.
The solution of large banded
linear systems of equations was the most time-consuming part of the present calcula-
tions. These solutions ran about
200
times faster on the CDC STAR
100
than on the
Caltech IBM
3701158
computer. The IBM machine was used for Reynolds numbers
from
2
up to
10,
for preliminary tests of the high-Reynolds-number method
on
very
small grids and
for
the final data processing and graphical presentation.
2.
Mathematical formulation
In terms of stream function
Y
and vorticity
o,
satisfying
=
ay/ay,
v
=
-a~/ax,
(1)
(2)
o
=
avpx
-
aulay,
the Navier-Stokes equations can be formulated
AT+@
=
0,
(3)
where
A
=
a2/ax2
+
a2/ay2
and the Reynolds number
Re,
based on the diameter
d,
is
steady
$ow
past
a
circular cylinder
Body
82
1
*
I
I
0
I
0
Upstream
Body
Downstream
I
ti
FIGURE
1.
The
z,
y
and
E,
71
plane.
Re
=
Ud/v.
The quantity
U
is the free-stream velocity and
v
is the kinematic coefficient
of viscosity. We will from now on assume that
U
=
id
=
I.
It
is convenient to work both numerically and theoretically with the deviation from
uniform flow
instead
of
with
YP.
kcr(X,Y)
=
Y(X,Y)
-
y
(5)
A
polar co-ordinate system can be introduced by the conformal transformation
1
[+i7
=
-ln(x+iy).
n
(The variables
6
and
7
are connected to the usual polar co-ordinates
r
and
#
by
r
=
ent
and
#
=
nrry.)
The Navier-Stokes equations take in these co-ordinates the form
A$
=
-n*rZw,
(7)
here
A
=
a2pp
+
aZ/a72.
We assume symmetry and consider only the upper half-plane. The half-plane minus
the cylinder gets mapped into the semi-infinite strip
0
<
7
<
l,<
2
0
(figure
1).
On the surface of the body, the boundary conditions are
1C.(O,7)
=
-sin
(w),
(9)
corresponding to vanishing normal and tangential derivatives of
Y.
Symmetry gives
@
=
0
and
w
=
0
as boundary conditions on 7
=
0
and
’I]
=
1.
Numerically, two bound-
ary conditions must also be supplied at some outer limit
rm.
The choice of these outer
conditions will be discussed in detail later. One commonly used simple choice is to use
the ‘free stream’, i.e.
II.
=
0,
w
=
0
on this boundary.
As
we will see, this choice is very
unsatisfactory. Even when applied hundreds of radii away from the body,
it
will lead
to significant errors in the flow field (in particular in the vorticity) right up to the body
surface.
822
B.
Fornberg
3.
Numerical difficulties for increasing Reynolds numbers
Sooner
or
later numerical attempts to simulate accurately the flow at increasingly
high Reynolds numbers run into factors which limit further progress. The existence and
the nature of such factors have rarely been discussed in the numerical literature on this
problem. We believe however that the following list contains the main problem in this
respect.
1.
How to implement boundary conditions at large distances.
2.
How to implement boundary conditions at the surface of the body.
3.
How to get a reliable rate of convergence for the numerical iterations (a rate
which does not deteriorate seriously with increasing Reynolds numbers).
4.
How to approximate the vorticity transport equation in a way which (i) is stable,
(ii)
is
at
least
accurate to second order, ideally allowing Richardson extrapolation
(or
deferred correction) to fourth-order accuracy, (iii) makes the overall iteration scheme
convergent.
The item
1
above will be the subject of a detailed discussion in
94.
Here, we will only
make brief comments on the other items and then wait until
5
5,
where we will see how
the proposed method handles them.
Problem
2
arises because we have two boundary conditions for
$
and none for
w,
where one for each variable might have proved easier to work with. Many different
techniques have been tried at the boundary.
It
is
our
impression that convergence
problems have been present at this boundary in many cases.
We believe that the limiting factor in most previous work is contained in problem
3.
The physical problem becomes unstable with respect to unsymmetric disturbances
around
Re
=
40.
With symmetry imposed, stability persists much longer (for how long
is not known). Numerical methods for the steady symmetric problem have in most
cases involved iterations between the stream function and vorticity equations. Such
iterations introduce an artificial time in which instabilities can be encountered quite
early without the symmetric problem being unstable in real time. We believe this
artificial time instability is the main reason why accurate calculations have not yet
reached high Reynolds numbers.
Our
approach
to
this problem is to solve the coupled
stream function and vorticity equations rather than using one equation at a time in
some iterative manner.
The vorticity equation in steady calculations has often been solved as an elliptic
system with a method
based on the successive over-relaxation method. Local
diagonal dominance must then be assured for the difference approximation of the
vorticity transport process. This has led to such methods as upwind differencing (see
Roache
1976,
p.
64
for references). In its simplest form, using simple uncentred
approximations for
aw/ag
and
aw/av,
the local accuracy
is
only first order, which com-
putationally is extremely inefficient. Different techniques to maintain second-order
accuracy have been given. Allen
&
Southwell
(1955)
give a scheme which can be
described as a continuous transition between upwind differencing and centred
approximations. The diagonal dominance is barely maintained at all instances. In the
limit of step sizes going to zero, the approximations become centred, and the whole
method is therefore formally of second order. Dennis
(1973)
has studied this method for
finite step sizes and has found that
it
works surprisingly well and may even allow
5.
How to get adequate but not wasteful resolution of all the different scales.
Xteady
$ow
past
a
circular cylinder
823
Richardson extrapolation. Another idea to obtain second-order accuracy is to extend
the approximation width in the upstream direction. Leonard
(1979)
has described
such a method. The general experience seems to be however that upwind techniques
compare unfavourably with centred approximations with respect to accuracy. There
is no generally accepted simple scheme which satisfies all the requirements listed under
problem
4.
The flow field shows a mixture of different scales
for
high Reynolds numbers. There
is close to the body a thin boundary layer, which separates and extends downstream.
Far away from the body, a narrow wake contains a sharp perturbation from free
stream while, in all directions far out, a slow but non-trivial perturbation exists. Usual
polar co-ordinate systems which are dense enough to resolve the wake far out will be
very wasteful in other directions. Our choice of grids will exploit the fact that
w
<
1
and
@
is very smooth and governed by a simple linear equation in most of the outer
field.
It
is not always easy to tell which of these difficulties has been the most pronounced
in previous steady calculations. We have indicated some guesses in table
1,
which
summarizes a number of previous contributions to the problem.
4.
Far-field boundary conditions
Two boundary conditions must be supplied at some finite, large distance
rm
from the
cylinder
(or
at infinity if a transformation has been made which brings infinity to a
finite distance). Since both the stream function and the vorticity equations are of
elliptic nature,
it
appears natural to supply one condition for each of the variables
along each edge.
Asymptotic formulas are known for
@
and
w
as the distance
r
increases to infinity
(see for example Imai
1951).
The leading terms are
Q
=
(&Rer)hsin$ny, erfQ
=
2n-h
e-@ds;
with
C,
is the drag coefficient. In
2,
y
co-ordinates, this means that, to leading order and at
large distances, the difference
@
between stream function
Y
and free stream
y
looks
like a simple source with equal outflow in all directions balanced by an inflow in a
narrow region behind the body. At high Reynolds numbers, this simple picture is valid
first at very large distances. Figure
2
show plots of
@
for
Re
=
2
and
Re
=
200.
The
ultimate directions of some of the lines (given by
(1
1)
when
C,
has first been calculated)
have also been marked. As the Reynolds number increases,
CD
decreases and therefore,
according to
(11))
the strength of the radial outflow. The flow behind the body is
almost stagnant and the inflow in the
9
variable does not decrease correspondingly.
This causes the large circulation in
@
that we see in figure
2
(b).
SD”