J. Fluid Mea.. (1980), vol. 98, parl4, pp.

819-855 819
Printed in (keat Britain

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,']inder
at Reynolds numbers up to 300. A new technique is proposed for the boundary 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 limitmg properties of the steady flow for increaaing 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-Reyno Ids-
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

oo2~1120/80/4471-1380 802.00 <rJ 1980 Cambridge University Preas

820 B. Fornherg
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) aBBumes 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 aBBumptions 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 discuBBes boundary conditions
at large distances. The final numerical method is presented in §5 with the obtained
results discussed in §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 discuBBed
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 CDe STAR 100 than on the
Caltech IBM 370/158 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 'I' and vorticity w, satisfying
U = 0'1'/fJy, v= - 0'1'lox, (1)
w = av/ox-8u/fJy, (2)
the Navier-Stokes equations can be formulated
6.'1'+w = 0, (3)

6.w+
_
Re o'Y ow 0'1' OW
)= 0' (4)
2 ( ox fJy fJy ox

where 6. = o'/ox'+ o'/fJy' and the Reynolds number Re, based on the diameter d, is
 Stemly flow past a circular cylinder 821

y Body

--- -',
,......... ,,
, "
,,
,, ,,
,, ,,
,,,
, \
,,
\
I
Upstream Body Downstream
- x

FIGURE 1. The x, y and g, 1Jplane.

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 = !d = 1.
It is convenient to work both numerically and theoretically with the deviation from
uniform flow
1jr(x,y) = 'Y(x,y)-y (5)
instead of with 'Y.
A polar co-ordinate system can be introduced by the conformal transformation

;+i'l = !.In(x+iy).
1T
(6)

(The variables; and are connected to the usual polar co-ordinates rand () by r = e"'
'1
and () ~ 1T'I.)The Navier-Stokes equations take in these co-ordinates the form

t.1jr = -1T W, (7)

Re ow. OW o1jrow o1jrOW

t.w ="2 1Tr(cos() -sm () ; (8)
( 0; a"
- 0; +
a" 0'1 0; )
)
here
t. = o'/o;'+o'/a,,'.
We assume symmetry and consider only the upper half-plane. The half-plane minus
the cylinder gets mapped into the semi-infinite strip 0 ,,;; ,,;;1, ; 0 (figure 1).
'1 '"
On the surface of the body, the boundary conditions are
1jr(0,') = -sin (1T'I), (9)

o1jr~i'l)= -1TSin(1T'I): (10)

corresponding to vanishing normal and tangential derivatives of 'Y. Symmetry gives

1jr 0 and w = 0 as boundary conditions on
'1 = 0 and '1 = 1. Numerically, two bound-
~

ary conditions must also be supplied at some outer limit r~. The choice of these outer
conditions will be discussed in detail later . One commonly used simple choice is to use
the 'free stream', Le.1jr ~ 0, w = 0 on this boundary. As we will see, this choice is very
unsatisfactory. Even when applied hundreds ofradii away from the body, it will lead
to significanterrors in the flowfield(inparticular in the vorticity) right up to the body
surface.
822 B. Fornherg

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 num~ricalliterature 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.
5. How to get adequate but not wasteful resolution of all the different scales.
The item 1 above will be the subject of a detailed discussion in §4. Here, we will only
make brief comments on the other items and then wait until §5, where we will see how
the proposed method handles them.
Problem 2 arises because we have two boundary conditions for >frand 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 ow/a/[,and ow/a", 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
 Steady flow 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 q 1
and 1fris 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 roofrom 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 1frand w as the distance r increases to infinity
(see for example Imai 1951). The leading terms are

,1. =-('l-erfQ)
CD (11)
2

w ~
_ CDRe 51e-Q' (12)
4../11 r '

with Q ~ (tRer)'sint11'l, erfQ = 211-' foQ e-"ds;

CDis the drag coefficient. In x, y co-ordinates, this means that, to leading order and at
large distances, the difference 1frbetween 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 1frfor Re 2 and Re
~ ~ 200. The
ultimate directions of some of the lines (given by (11) when CDhas first been calculated)
have also been marked. As the Reynolds number increases, CDdecreases and therefore,
according to (11), the strength of the radial outflow. The flow behind the body is
almost stagnant and the inflow in the 1frvariable does not decrease correspondingly.
This causes the large circulation in 1frthat we see in figure 2 (b).
824 B. Fornberg

;::
S
..!J
..c
o
o""-
~
.
~o
. 0
§> 'ED
o e

~'"
~oo 000
~ ~"

00 00 0IN 00
00
.."
00 000
oo~
000
""00 "'~ O..,jl
__1>':1

oou:.ooooo
__ 00
00_0000<:0
_lQ~1.O o~
o
oo
'"

a;
00
~
'"
 Steady floW past a circular cylinder 825

(0)

-3,0
0.5 ......-
-0,25 3,0

-0'25 2.0
1.0

-0'25/
(b)
-0,5

0.5
0.25
FIGURE 2. The 1jf field for (a) R. = 2 and (h) Be ~ 200.

The vorticity far from the body is concentrated in a thin streak downstream with
maximum strength inversely proportional to the distance. Since the question of
boundary condition for w turns out to be mnch simpler than the condition for ifF,we
consider that case first.

Boundary condition for w at roo

Vorticity is transported along streamlines and dissipated by viscosity. It is only
because of the dissipative term that a boundary condition at the outflow side is needed.
If wrong boundary values are supplied, the introduced errors will decrease exponen-
tially upwind from this boundary. The only concern we need have is purely numerical.
The computed values must not fall apart into staggered mesh oscillations far into the
region. That would normally happen if values for w are given on the boundaries and all
826 B. FornlJerg
approximations are centred and accurate to second order. Upwind approximations of
the vorticity transport would solve the problem, but with a loss of accuracy. Our
choice is to use centred approximations, but impose fJw/fJ; = 0 on the boundary (by
equating the values on the last two mesh lines). As we noted, the problem did not need
any information from the outer boundary and this method prevents the possibility of
staggered mesh oscillations.

Boundary condition for >jroJ,roo

The question of boundary conditions for >jris more involved. The quantity >jrsatisfies in
1/ co-ordinates
f"
!!.>jr+ 1T""W = O. (7)
Actual physical information has to be supplied through the outer boundary, and it will
propagate immediately to the interior. A distance r = 100 corresponds to; = 1.47 and
r = 10' to f, = 4.40. The use of a large rootogether with inaccurate boundary values
does not give a reasonable accuracy. Also, as we will see, we want to avoid numerical
calculations for large values of r because of increasing differences between the scales
in the problem.
We will consider four different choices of boundary condition for >jr.
(a) Free stream (using >jr= 0 on the outer edge).
(b) One term of the Oseen approximation (equation (11»).
(c) Normal derivative zero (o1fr/of, = 0).
(d) 'Mixed condition' connecting >jrand o1fr Iof,on the boundary.
We will briefly discuss the justification for each of these choices and then use two tests
to assess their accuracy and usefulness. These tests are:
(i) to plot >jrand w when we move rooin and out and see how much the flow picture
changes;
(ii) to display the maximum vorticity on the body surface as a function of roowith
the different choices of boundary condition.

(a) Free stream

We know that >jrforlarge f,will tend to a linear function in 1/ (with a jump at 1/ = 0) and
that the values of>jr satisfy a Poisson equation in f, and 1/. The use of the free stream,
i.e. imposing >jr= 0 at some small ;00' seems very crude. Nevertheless, as table 1 showed,
this has been used in most cases even for very small values of roo.Therefore the method
has to be studied for accuracy. Figure 3 (a) shows the result of test 1 at Reynolds
number Re = 2 using roo = 23.1 and roo= 91.5. We can easily conclude that large
errors have been introduced. In particular, the vorticity field has changed even up to
the body surface. Figure 4(a) shows that the maximum vorticity on the body surface
is as much as 20 % in error with respect to roo= 23.1. Figures 4(b,c) show that large
errors also persist for high Reynolds numbers. (Note that larger roohave been used in
figures 4 b, c.)

(b) One term of the Oseen approximaJ,ion

The use of this boundary condition requires knowledge of the drag coefficient CD' This
can be evaluated from the current approximation of the flow field at each iteration by
means of a line integral around the body. Figures 2 (a, b) indicate that this boundary
 Steady flow pa8t a circular cylinder 827

, ,
~ 9
'"

0-
-
".. -,- N, M ..~
I , (.
, , ,
~ 9
'" '"

0-
Ot-
..
828 B. Fornberg

,,
].40 (aJ

.
~iii }-35
Free
stream
X,
,,
,,
>-
.., \ ,
~ I-3D \
c \
c
One-term \
~ 1.25 \
~ Oseen
'X ... X...
~ ...
:.
::;
1.20
a~
-=O-X---i(-
a! _
...
-X~
---------
].15 '0
2 5 10 500 ,
1.10
o 0.5 1.0 1.5 2.0

\
.u 5-70
Free
stream ~
\
(b)

<'!
Jj \
,
..,>-
0
"'c0 5.65
.
:0
\\
\
"'.
"'.
X
C \
." -:i
. .
...
...
"€ Mixed ........
-_
.
0
~0
> conditions
-x
n___::~::8 _ K-
" 5.60
.., :z ___nnhn
c
::; ,-
"" ,, iJ1/1 ;!- -:t 00
-=0'. &, 0
! a~ ~ ~ ~
5.55
2 5 10 20 50 100 200 500 1000 2000 5000 ,
0 0.5 ].0 1.5 2.0 2.5

Figure 4 (at b). For legend see page 829.

condition can be expected to be accurate at most for very low Reynolds numbers.
Figure 3(b) and 4(a) show however quite significant errors already for Re = 2. We did
not try this method for higher Reynolds numbers.

(c) Normal derivative zero

This condition is as easy to apply as free stream (for instance, it requires no evaluation
of CD). It is a 'softer' condition than >fr= 0 and is consistent with the picture of
equilines of>fr eventually going in the ~ direction (radially out from the body every-
where apart from in the wake where it is radially in towards it). Test number 1 in
figure 3 (c) shows this method to be superior to the previous methods at a low Reynolds
number. Both the justification and the actual performance of this boundary condition
decrease however with increasing Reynolds numbers as figures 2 and 4(a, b) show.
 Steady flow pa8t a circular cylinder 829

1
12.3 vi (c)
:gl
OJ
12.2 ~I
\
1
X
Free stream \
u 12.1 1
" ~I \
'" .,
" 12'0
>-
~I \
\
"0 Mi~ed\
conditions X,
~" ........
0 11.9 I
x- :;; 1_ _ ______
."
'" : ,x"-
'f0 11.8 '01 /
1 /
."
:E
>

11.7 ~=J
3~
1 X

", v
0'"
11.6 "I <io 0-
:! N
"', v
2 5 10 20, '"
50.100 200 500 1000 2000 5000 r
11.5
0 0,5 ',0 1.5 2.0 2.5 3.0
FIGURE 4. Comparison between different outer boundary conditions for 1jrat (a) Re = 2,
(b) Re ~ 40 and (0) Re
= 200.

We feel that, for practical purposes, this condition completely solves the problem of an
easy and accurate boundary condition for >frat Reynolds numbers up to about 40.

(d) 'Mixed condition' metlwd

Let us start the description of this technique with some observations. In g, 7Jco-ordi-
nates, the vorticity decays like O(e-") along the g axis. This vorticity enters as right-
hand side in the equation (7) for >fr.These small values of ware scaled up by the factor
r2 which is large at large distances. We know that the net effect of this is that >frcon-
verges to a linear function of 7J(with a jump at 7J= 0) as r tends to infinity. From
figure 2 (b) however we know that, for reasonable distances and high Reynolds numbers,
>frbehaves in a non-trivial way for all values of 7J,not just for 7J O. It is desirable to
apply a boundary condition so close to the body that the thin '" transition around
7J 0 and the slow variation for all 7Jcan both be resolved by the same grid. The
~

further out the boundary condition is applied, the more severe will the mixing of scales
become.
Figure 5 shows what wand >frlook like in g, 7Jco-ordinates at He = 200. A dotted line
is drawn where we would like to implement the boundary condition. According to our
previous discussion of boundary conditions for w, we can expect to have fairly accurate
values of w available close to this line. There are two reasons why >fr,determined from
the linear equation (7), takes non-zero values below this line. Let us write >fr= >fr, + >fr2
with the two components >fr,and >fr2 corresponding to the two reasons below.
(i) There is a thin streak of vorticity along the g axis. It decays and gets thinner but
the factor r2 in (7) enhances the influence of vorticity at large distances. The net effect
of this streak of vorticity is that 1/1for large values of Econverges to a function of 7Jwith
a jump of 7J= 0 (as given by equation (11)). In the thin region where vorticity is

)
830 B. Fornberg

w
0.6 ~=l' o r~=o ~=l
I

Path of extremum of
,I vorticity for different T 10
,
,,, 20
- End of recirculation region-
40

_ - - Distance from body at which-

we wish to implement outer
100
boundary conditions

200

400
2

1000

2000

4000

Limit Limit
.-0-42 .-0.25
FIGURE 5. Contours of w and Vrat Re = 200.

present, 821r,/8~2 is very small compared to 821r,/8'12 so 1r, satisfies approximately

(13)

for constant ;.
(ii) The values of 1r close to the body influence the values of 1rfar out since they are
coupled by equation (7). Far out, we can think of this influence as a potential flow with
1r2 satisfying
821r2 821r2
0.
_ (14)
8f'+8'12-

(The right-hand side in (7) is already taken care of in equation (13) for 1r,.) The bound-
ary conditions for 1r2are 1r2 = 0 at the three sides 'I = 0, 'I = 1,; = 00 and unknown at
~ = ~oo. As figures 2(b) and 5 clearly show, this term 1r2is important along the full
 Steady flow past a circular cylinder 831
extent of the outer boundary. The Oseen approximation assumes the flow in the wake
to depend on the drag only and fails to describe this part >fr,of the solution for >fro
These observations lead us to the following method:
(a) We suppose a guess >frofor the variable >fris available at I; = I;~.
(b) With w known at two mesh lines close to I;~, we can use (13) to solve for >fr,
(two-point boundary value problems) on these lines and thus get an approximation for
o>fr,lol; at I; = I;~.
(c) We evaluate >fro - >fr,at I;= I;~.Since >fr,satisfies (14), we can easily (Fourier
>fr, ~

decomposition) find o>fr,lol;at I; = I;~. (Note that this does not require any numerics to
be performed beyond I;~.)
(d) Given the guess >fro >fr,+>fr,we have thus found o>fr,lol;and o>fr,lol;,i.e. given
~

>frwe know o>fr101;. This is exactly the right form of a mixed condition to be imposed on
(7) at I; = I;~. We can put this relation in the form: Given >fron the next to the last mesh
line, we can find it on the last line. That formulation is ideally suited for iterative
methods to solve (7).
Figures 4 (b, c) show how this condition works for He 40, 200 respectively. We did
~

not try it for Reynolds numbers lower than 40 since the simpler condition o>frlol; = 0
worked sufficiently well there.
As a summary, it appears that the use of the simple o>fr101; = 0 for Reynolds numbers
up to about 40 and the mixed condition for higher Reynolds numbers solves the
question of boundary condition for >frat large distances.
Dennis (1976) also points out the insufficiency of the free stream as the outer bound-
ary condition and proposes a technique based on matching w, owlol;, >fr and o>frlol;
between an inner region governed by the Navier-Stokes equations and an outer region
governed by a simplified version of the equations. The approach suggested for the
outer region differs in many respects from our observations of the flow pattern and OUf
conclusions from it. In particular computations are performed throughout this outer
region (to find both >frand wi, and the simplified equations have terms in 0'101;'
dropped. The term o'>fr101;' is assumed to be small compared to o'>fr lo~' at large
distances. This leads to parabolic equations which can be marched for both >frand W.
Both the economy and the justification of such an approach seem unclear to us in
the light of our previous discussion.

5. Numerical method for steady solutions

The main problem that a method for steady solutions at high Reynolds numbers
must solve is the one of assured fast rate of convergence for the discrete approximating
equations. Since the real physical problem may be unstable, any iteration technique
which can resemble artificial time may have convergence problems. This danger is
entirely circumvented with the use of Newton's method for the discrete nonlinear
system of equations. With an initial guess close enough to a simple solution, quadratic
convergence is assured. Many iterative methods fail to converge as moving wave
patterns appear in the artificial time. Because of the quadratic convergence, such
waves cannot appear in this method. The price to be paid for this guaranteed perform-
ance is a comparatively high computational cost. This is no longer a major concern
with the introduction of extremely powerful array processors.
We have already remarked that the problem involves different scales. Vorticity is
832 B. For1Iberg

,
(.) \
lOi jI
-~-
IOj /"
-"'~--
:
I

51 \
-10 -5 0 5 10 20 30 40 SO 60 70 80 90

(b) (,)
IOi
2i

(d)

200i

,
,
-400 -200 o 200

FIGURE 6. Steps in the mapping of the inner region, c = 0.2. (a) x plane;
(b) P plane; (c). plane; (d) 'plane.

only present at the body and in a streak downstream from it. This region is the only
one where the problem is trnly nonlinear. Outside that region, we have only to solve

!lift = O. (15)

The solution to this equation only changes on a large scale. This clearly suggests the
use of two different grids with different numerical methods on them.
The dotted region in figure 6(a) is a suitable computational region for the nonlinear
part of the flow field. A conformal mapping of a region of such shape to a rectangle can
be achieved explicitly. For instance, setting p = x' maps the upper half-plane minus
the unit circle (figure 6a) into the region in figure 6(b), and z = c(p-l/p) maps the
segment of the unit circle onto the imaginary axis and leads to region in the figure 6 (c).
In these figures and in the following computations, the constant c was given the value
0.2. Finally, we want to obtain a mapping which resembles the original x, y system far
downstream and the z co-ordinates close to the surface. This suggests a final step
I:= z(1 +Z2)C-3. The real and imaginary axis are kept from the z plane but the cubic
term cancels the cube root in the initial step (figure 6d).
We put a rectangular grid in the I: plane as is shown dotted in figure 6(d). Its images
in the other planes are shown in figures 6(a-<:). (Changing the constant c affects the
width of the mapped wake region in the x plane.) The complete mapping from x to I:
can be written as

(16)
 Steady flow past a circular cylinder 833
and it can be inverted in closed form (choosing principal branches) by
t = c31;,

1 (17)
p= a-3a'

With the notation x ~ [; + in and 1;= a + ip equations (7) and (8) take, in the 1;plane,
the forms
02V 02V
0, (18)
oa2 +OP2+w/M =
02W 02W He
oa2+Op2+"2 oa 'op-op'oa(OV oW oV ow
(OW oa oW OP
)/
- oa' o[;+op'o[; M ) = 0, (19)

where M = (oa/o[;)2+ (OP/O[;)2,The functions oa/o[; and oP/o[; can be evaluated as the
real and imaginary parts of the complex function
01; _ 01; 82 op_ (1+3z2)(P+l/p) (20)
ox-oz'op'ox- 3xc2
In a further transformation step, the grid was compressed close to the body and
stretched far out by the change
a = c,(&" -1), (21)
where c, and C2are constants (0,008 and 0,118 respectively in our work), A rectangular
region in the y, p plane is now discretized equidistantly, On this rectangular grid in the
y, p plane, we approximate the transformed Navier-Stokes equations in the straight-
forward way with centred, second-order approximations, On this inner grid, V was
given on all the boundaries and also 0V / Oy on the body surface, For w, no condition was
given on the surface, w 0 on top and bottom boundaries and ow/Oy = 0 on the right
~

side boundary (large distance), When we used our fine grid (65 x 114 points in the p and
y directions respectively, corresponding to an outer boundary approximately 600 radii
from the body), the condition involving the normal derivative OV/Oy at the surface
was approximated using a centred approximation with an artificial point inside the
body, Use of (18) at the surface allowed this inside P?int to be eliminated, giving a
relation between w at the surface and V at the mesh line nearest the surface, The less
accurate method of approximating OV/Oy by a one-sided stencil out from the body
(and no artificial point inside tI,e body) is also satisfactory in this problem, and was
used in our calculations on the coarse inner grids (33 x 44, 33 x 52 and 33 x 86),
The discretizations we have described lead to a large coupled nonlinear system of
equations with the same number of equations as unknowns. Newton's method is now
used to solve this system, To describe the structure of the Jacobian matrix for this
system, we need some additional notation, Assuming the grid in the y, p plane is of
size N x M, we let

37 PLK98
834 B. F<Yr11herg

,
I
~ ~
! I
~ £.)2 W3W4
i I I
WN_I

I
'I

,, .
,,
, .. .
... .
RWN_l
~

- RBC~,
Rw,
Rw,
w,
W,
Rw,
.
W, .
.. ..
RWN_l
~ ~
RBCWN
I I U-I ~BCI/I;

FIGURE 7. Structure of the linear system in Newton's method.

be vectors ofIength M - 2 containing the values of w along the M - 2 interior points in

the fJ direction on the N different grid lines. (For example WI contains the values of the
vorticity on the body surface.) Similarly. we use

for the function ifr. Further. let

Rifr.. Rifr3' RifrN_1

be vectors containing the residuals in the vorticity and stream function equations
centred along lines 2 to N -1. Finally.

are the residual vectors for the boundary conditions on lines 1 and N. (Since ifrl is
known. we do not enter ifrI as an unknown. RBCifr I refers to the second condition,
iJifr/iJy given.) .
The Jacobian coefficient matrix takes a shape which directly reflects the shapes of
the approximating stencils. The particular ordering of the equations and the unknowns
in the linear system shown in figure 7 may at first look irregular. This ordering, however.
leads to a coefficient matrix which easily simplifies to a form well suited to numerical
solution. As a first step in solving the system, we superpose multiples of the equations
in the first N - 2 (block) rows to eliminate all entries in the bottom right major block.
The bottom left major block then assumes a 13-diagonal form, which is completely
confined inside the dotted region in figure 7. Figure 8 shows the structure of this band
 Steady flow plUJt a circular cylinder 835

r T ,
,,
/
/
r ,
,, /
/
, ,
,, /
,, /

, , /
,, / ,,

:,/ ,
/' I
, / / /'
I
I"
..."1
I I'
,,'I
I
t' (" (" I
I
:
I
I
I
I
I
1
I
I
I
I
I I I I
I

I ,'I
:
I
:
I r
I
-1
:
t ,.. I ,/ I ,,/ I
1,,/ I ,/ 1,../ "I I
1/ I" f
, " I"
, ,
, ,
, /
, /

,
, /
, ,/ /
/

/
/
/
,
/
i
,
/
/ .L

FIGURE 8. Structure of the l3-dia.gonal banded

" coefficient matrix.

General grid size Coarse grid

.
Grids used in present work

Fine grid
,
MxN 33 x 48 65 x 114
Size of reduced Bandwidth 4xM-7 125 253
banded system to Number of
solve equations (M-2)x(N+l) 1519 7245
Computer time for CDC STAR-l00,
nwnerical solution LV -decomposition 2.558 29.33 s
Back substitution 0.22 s 1.25 s
IBM 370/158'
(estimate)
L U -decomposition 8 min lh40min
TABLE 2. Sizes and computer times for banded linear systems.

matrix (with the lines displaced sideways so that the diagonals go vertically down in
the picture). The approximate cost of solving this banded system is shown in table 2.
The L U -decomposition for the 65 x 114 grid requires approximately 470 x 10" additions
and multiplications. Using FORTRANwith no special optimization, this required 29.33
seconds (including all non-arithmetic operations like loop control, partial pivoting
27-2
836 B. Fornherg

FIGURE 9. An example of the superposition of outer and inner grids.

etc.) on the ODO STAR 100, giving an average speed of about 16 Mflops (Million
floating point operations per second, 64 bit wordlength).
The 13-diagonal system we obtained after the initial elimination could have been
generated more directly first by analytically eliminating'w between (18) and (19) and
then by approximating the resulting fourth-order biharmonic-type equation for yrwith
standard centred approximations. It is not clear which is the better method by which
to generate the final banded system.
As we have observed earlier, yr behaves non-trivially in all directions from the body
even at large distances. If w is known inside the wake, ifr can be found everywhere
from the simple equation (7). To solve this equation, we use a modified polar grid
g, r/ with g, r/ given by
(22)

On this grid, we discretize g and .,' equidistantly. Figure 9 shows a typical case of inner
and outer grids superposed onto each other.
We can now describe the complete method. Assuming we have an initial guess for
yr and w, we carry out the following steps:
1. Perform one Newton iteration on the inner grid. This gives an improvement in w
(and in yr).
2. Interpolate w to the outer grid using two-dimensional cubic splines (fourth-order
accurate).
3. Solve the Poisson equation for yr on the outer grid (with the' mixed' boundary
condition at the outer boundary). This solution is obtained by 'black-red' ordered
successive over-relaxation (SOR).
4. Interpolate these new yr values back to the edges of the inner grid. Return to
step 1.
When the grid sizes become large, step 1 is the most expensive one. It also converges
quadratically while the iterations between the grids only converge linearly. For these
 Steady flow past a circular cylinder 837
reasons, we perform only one Newton iteration for each outer iteration. In the step 4,
over- and under-relaxation can be used when the values of If on the edge from the
previous iteration are corrected with better approximations. The use of relaxation
factors between 0,5 and 1.5 (rather than 1.0 corresponding to using only the new
approximation) led in all cases to an overall convergence with a factor of about 10
during each full cycle over steps 1-4. Six to eight digits of accuracy in solving the non-
linear system (truncation errors are probably larger) thus required about eight itera-
tions. Complete computer times, including all overhead, ranged from between 1 and
2 minutes for a 33 x 48 inner grid and between 6 and 8 minutes for a 65 x 114 inner
grid. These times were independent of the Reynolds number.
Direct jumps in Reynolds numbers 20-40--100--130--160--200--230--260 were used
with the 65 x 114 inner grid. In the proximity of Reynolds number 300 smaller con-
tinuation steps were needed: 260--275-285-290--295-297'5-300. If the initial guess was
not sufficiently close to the solution, the Newton iterations diverged, first in the areas
where the guess (the solution at a previous Reynolds number) was least accurate. This
was in the region where vorticity recirculation was observed. The changes were there
very rapid as a function of the Reynolds number. Since Reynolds number 300 also is
close to the upper limit at which results are reliable due to truncation errors (calcula-
tions on a grid twice coarse began to deviate significantly shortly before Reynolds
number 200) and since further mesh refinement would have been too costly, no effort
was made to reach still higher Reynolds numbers. We found no indication of any kind
that any principal problem was present, which might have prevented further increases
in Reynolds number, had a finer grid been used.
We will now discuss how the proposed method handles the difficulties that were
listed in §3.
1. Boundary conditions at large distances were discussed in detail in §4. We applied
the' mixed' boundary condition to the outer grid only. A Dirichlet condition, obtained
by interpolation from the outer grid, was used for the inner grid. Since the outer grid
was handled by SOR, the' mixed' condition was easily implemented.
2. The problem with boundary conditions at the surface of the body was that we
had two conditions for If and none for w, when one for each variable might have been
easier. With Newton's method, all that is required to obtain quadratic convergence is
the correct number of equations, an isolated solution, and a sufficiently good initial
guess. The fact that there was not one condition for each of the variables If and Whas
become completely irrelevant.
3. Quadratic convergence is assured in Newton's method for isolated solutions. This
excludes every possibility that waves or similar features may move in artificial time
during the iterations. (The linearly convergent outer-inner iteration to find boundary
values for If for the inner grid could conceivably have failed to converge.)
4. Since convergence on the inner grid has been assured, there is no longer any need
to introduce special stencils or other methods to stabilize the vorticity transport
approximation at the expense of the accuracy. Second-order centred approximations
satisfy all the three requirements. Since the interpolation between the grids was
fourth-order accurate, overall fourth-order accuracy could have been obtained by
Richardson extrapolation.
5. The chosen grids, with the inner one refined at the body surface, seemed to give
suitable resolution in all regions. It was possible to use a much cheaper numerical
838 B. Farnberg

,---,,--

-~.._-----_.

25 30 35 40 45 50

----
~-----

25 30 35 40 45 50

25 30 35 40 45 50

.
,~

01 2 5 10 15 20 25 30 35 40 45 50

~
(e)~'
~.
'.--~r.\~
012 5 10 15 20 25 30 35 40 45 50

if)

012 5 10 15 20 25 30 35 40 45 50
FIGURE 10 (a-f). For legend see page 839.

technique iu those regions where the equations were without special nonlinear
difficulties.

6. Results
There are many questions on the limiting flow fields at high Reynolds numbers that
remaiu open. Although this investigation does not solve any of them, distinct trends
for iucreasiug Reynolds numbers can be seen iu several cases. The most striking feature
we notice is a recirculation of vorticity from the end of the wake region, starting around
a Reynolds number of 260. This leads to quite sudden changes iu some flow quantities
like the length and the width of the wake region while other quantities, such as the
vorticity and pressure distributions on the body surface and the drag, seem quite
unaffected. The rest of this section contains a discussion of the results, mainly focusing
on the numerical aspects rather than on the physical consequences of the observations.
Aswas mentionediu theiutroduction, the results for Re = 2to 10 were not obtaiued by
the present method, and are of comparatively low accuracy and graphical
resolution.
Figures 10 and 11 show the streamliues and the vorticity fields for Re = 2 to 300. We
notice a recirculation of vorticity, clearly developed beyond Re = 290 but suggested
by the vorticity fields already at Re ~ 200, 230. This widening of the region with
vorticity causes the separating streamline to close nearer the body. (The small' island'
of low vorticity in figure 11 (I) should be connect..d with the thin streak approaching
it from the left. The vertical mesh resolution was not sufficient for representing this
streak continuously.)
To ensure as far as is possible that the observed features of the flow fields are real
and not peculiarities of the discret.. approximations, .we checked the independence of
these solutions to changes in the grid sizes. Halving the densities of both outer and
inner grids from 129x 132, 65 x 114, respectively, to 65 x 52, 33 x 48 (also, corres-
ponding to a smaller total region) only offered sufficient resolution power up to a
Reynolds number of around 150. Figures 12 (a)-(h) compare the flow fields using these
different grids for Re~ 20 and 200. The discrepancies are quite large at Re ~ 200 but
the figures show that the beginning of a trend towards a recirculation of vorticity, is
present on both grids. At Re = 290, only the fine grid could be used, but we have a
parameter c available in the conformal mapping for the inner grid. This constant c was
changed from 0' 2 to 0, 25 and 0.3. This corresponds to halving the vertical extent of the
irmer grid (i.e. doubling the vertical grid density). Figures 13 (a)-(f) show the effect
 840 B. Fornherg

(aJ

30 35 40 45 50
10 15 20 25
012 5

(bJ

20 25 30 35 40 45 50
012 5 10 15

(c)

-- 25 30 35 40 45 50
10 15 20

(If)

20 25 30 35 40 45 50
10 15

(e)

5 10 15 20 25 30 35 40 45 50

if)
+ 50
01 2 5 10 15 20 25 30 35 40 45

FmURE 11 (a-J). For legend see pa.ge 841.

these changes had on the computed flow fields. We note a slight movement of the line
for equal vorticity. The vorticity in this outer area is almost constant and the point-
wise change is much smaller than the movement would seem to suggest.
Figures 14 (a)-(d) show details of the flow pictures (with data taken from the inner
grid) close to the body for Reynolds numbers 20 and 300 and figure 15 the distribution
of vorticity on the body surface. The structure and position of the separation point is
not completely understood. Theory based on the Helmholtz-Kirchhoff model pre-
dicts that, as the Reynolds number goes to infinity (Bordetsky 1923) the separation
point may move forward to an angle of 550 (measured from the front). Our results (and
this possible limit) are shown in figure 16. The width of the wake bubble and the
position of its end (measured from the centre of the circle) are shown in figures
17 (a, b). Table 3 compares the length with earlier results.
There are several ways to calculate the pressure in the flow field and the pressure
on the body surface in particular. In g, "Ico-ordinates, the pressure satisfies (r = e")

(23)
P, = - ~e w. + {'I', '1'" - '1'. '1',. + 1T('l'1+ 'Y:nj(1Tr)',

(24)
 Steady flow past a circular cylinder 841

(g)
+
.01 2 5 10 15 20 25 30 35 40 45 50

(h) ~
:-----------....!._---
-=
-----------
- -
------ ----
:1 r
5 10 15 ~25 T
20 30 35 40 45 50
-~~ --~---
--- - -- - ---~--,- -------
~
--- --- =-----
(I) --------
=-------
-- ------- -----
+ I -
I --
=r--=--=
=-=r-- ------------.-------
~----- ===---=-
0'] 2 5 10 15 20 25 30 35 40 45 50

, -----

+ I-,-~=-=-r= =-=~ ----

012 5 10 15 20 25 30 35 40 45 50

(k) ~~---
~=
+ ~-- ,---
012 5 10 15 20 25 30 35 40 45 50

~- -
(/) =
-~
-------
(: ~

+
0;1:2 5 10 15 20 25 30 35 40 45 50

FIGURE 11. Vorticity fields at different Reynolds numbers. For (a}-(l) see figure 10.

and the Poisson equation

- (Ptt+P..) = M(;(UY,)' - (;II_'Yt) ('1';,+'Yt)) , (25)

with boundary conditions

Pt = - ~e w, on g = 0,

P, = ° on 'T/=
°
and 'T/= 1, (26)

P = ° on g = lioo (infinity). )
In equation (25), '1' can be replaced by ifr (with no further changes in the equation).
With ifr and therefore '1' known from our calculations, we employed three different
methods to find the pressure at the front stagnation point:
1. We solved the Poisson equation (25) with boundary conditions (26).
2. We evaluated the line integral (23) from far upstream, g with 'T/= 1.
= goo,to g =
3. We evaluated the line integral (23) from g = liooto g = g, with 'T/=° I, with (24)
along a half-circle to g = Ii"'T/ = 0, with (23) to g 0, 'T/= 0, and finally with (24) along
=
the surface to g= 0, 'T/ = 1. Figure 18 shows this path.
Table 4 gives the results obtained with these three methods.
Method 1seemsquite accurate for low Reynolds numbers only. A loss of accuracy for
842 B. Fornberg

----------
_____
(aJ
- =~ r-; ~-
.
-~-_._-
-~
- --
----

on
oI 2 5 10 15 20 25 30 35 40 45 50

012 5 10 15 20 25 30 35 40 45 50

(e) -~~-----------
~= ) --- ~ --_/~~
+
012 5 10 15 20 25 30 35 40 45 50

(d)

+
012 5 10 15 20 25 30 35 40 45 50
~-;==~ - - -'=. ==~=' !I~~~~-=c~~-~-==-=-:
(e) -_
:- ---=====-~~~~=~:~:=~~ --
~c-
.
oI2 5 10 15 20 25 30 35 40
~

45
-=

50

~~~.::;:;;:;:::--=~~-== =~;.;-=~ =-~---:==~~~=~

(f)
~~~_~=':~~=:~_~~~~~=~~=~~_'
01 2 5 10 15 20 25 30 35 40 45 50

(g)

012
. 5 10 15 20 25 30 35 40 45 50

(h) -===-~-===--co_~_____ ___ u__

- - -- - --- -

.
012
.

5
-----

10 15 20
-----

25 30
-
--

35 40 45 50
FIGURE 12. Compa.rison between results on coarse and fine computational grids. Be = 20. (a)
coarse grid, (b) fine grid, (e) coarse grid, (d) fine grid. Re ~ 200, (e) coarse grid, (j) fine grid,
(g) coarse grid, (h) fine grid.

high Reynolds numbers is caused by an increasingly singular behaviour of the right-

hand side of (25), in particular towards the end of the wake bubble.
Method 2 works well for high Reynolds numbers, but cannot handle low Reynolds
numbers well. This is because vorticity is then present quite far in front of the body.
This integration can only use values from the outer grid, since the inner grid does not
reach the inflow axis. Vorticity was assumed to be absent in this outer region. Although
interpolation of w from the inner grid was performed, the term - 2w,/ Re was not avail-
able with high accuracy.
Method 3 was chosen in order to avoid the problem in the two previous methods. It
probably gave the best overall accuracy. A consistency check was obtained by chang-
 Steady flow pa8t a circular cylinder 843

(a)
<- c:::: <" -------:3-S
+
012 5 10 15 20 25 30 35 40 45 50
~~" =c-
(b)
~.

+
~-
c::: c = > '-'S
012 5 10 15 20 25 30 35 40 45 50

(c)
'"
~=
5 10 15 20 25 30 35 40 45 50

(d)
~~.
012 5 10 15 20 25 30 35 40 45 50

(e)

+
~012 5 10 15 20 25 30 35 40 45 50

if) ~~----::--=
C ---
+
012 5 10 15 20 25 30 35 40 45 50
FIGURE 13. Numerical solution at Be = 290 for different values of the parameterc in the mapping
for the inner grid. (a) 0 = 0.2, (b) 0 = 0.25, (0) 0 ~0.3, (d) 0 = 0.2, (e) 0 = 0,25, (f) 0 = 0.3.

ing the radius of the outer circular part of the path. In no case did a change within the
range r = 2 to r = 16 lead to more than a 0.8 % variation in the calculated pressure.
The integration from the back to the front stagnation point used data from the inner
grid, and should be of high accuracy. For this part of the path equation (24) simplifies
to 2
~-&~ ~
since '1", and '1". both are zero on the surface.
Table 4 gives the values we have accepted as most likely for the front and back
stagnation point pressures. They are also displayed in figure 19. We note that the
trend for the rear stagnation pressure appears to be to pass zero towards positive
values. This is consistent with Batchelor's limiting solution but is not consistent with
the Kirchhoff-Helmholtz solution, where the preBBure is zero everywhere inside the
wake. (A Reynolds number of 300 is much too low however to give reliable evidence on
the final limit.) The pressure distributions over the body surface is shown in figure 20.
Table 5 compares the front stagnation preBSures we obtained with previously published
results.
844 B. Fornherg

0'5

0.2
0.1

~O~I~

(a)

(b)
:\\
0.\
-0'2_

0 -O'L
I I
'I 2
FIGURE 14 (a, b). For legend see page 845.

The drag D on the cylinder is given by

D = pU"aCD, (28)
where p is the fluid density, U the free stream velocity, a the radius of the cylinder and
CD the drag coefficient. This non-dimensional coefficient CD can be evaluated by a line
integral around the body, choosing for example a circle with radius r. The formula for
CDis

(29)
 Steady flow past a circUlar cylinder 845

- -

(d)

FIGURE 14. Details of the flow pictures for Be = 20 and Be = 200. Be = 20, (a) stream function,
(b) vorticity. Re
= 200, (c) stream function, (d) vorticity.

Evaluation of Cn using different radii r gives to some extent a consistency check on the
solution. This integral is however subject to severe numerical cancellations, especially
around the end of the wake region at high Reynolds numbers. This can be seen by direct
inspection of the integrands. It may also be concluded from the fact that the trape-
zoidal rule and Simpson's rule on the same data give quite different results (with the
first ones more nearly independent of the radius). In figure 21 we show Cn as a function
of r for Re 20 to 300. With the exception of the last case, there are three curves corre-
sponding to the coarse and fine grids (outer grids 65 x 52 and 129 x 132 respectively)
and the Richardson extrapolated result. The differences between these curves and
their independence of r for small radii suggest that
at Re = 20, Cn ~ 2.0001 :t 0.0002,
at Re = 40, Cn ~ 1'4980:t 0'0005,
at Re = 100, Cn = 1.058:t 0'001,
=
at Re 200, Cn = 0'829:t 0'002,
at Re = 300, Cn = 0.722:t 0.003.
These values of Cn are displayed in figure 22.
846 B. Fornherg

-15

-10

-5

-I o 8
o 1.0
+1
0,5 o ~

FIGURE 15. The vorticity distribution over the body surface.

.3

.. ----------------------
.
"i!: 2

..
.~
e
0.
...x'
..x.X
'" ,X'
, , X'
X

o
2 4
10 20 40 100200 400
Re
FIGURE16. The position of the separation point.

The error estimates are, however, hy no means completely certain. For example, they
do not take into account possible errors caused by a too restrictive computational
region in front of the body for low Reynolds numbers. We do h..wever believe that
error to be small.
The irregularity close to the surface in figure 21 comes from the fact that we used
values oflfr and w from the outer grid only. That grid was designed for evaluation oflfr
at large distances and is too coarse for drag calculations close to the surface. Table 6
compares our estimates for CD against others published previously. We believe the
discrepancies in almost all cases are mainly due to the inaccurate ways in which the
outer boundary conditions have been handled. For example calculations at Re = 2
 Steady floW past a circular cylinder 847

4 {aJ

:0
.0 3
.
.
~.0

"~.
-0 2
.0
CO

o 20 40 100 200 230 260 290 300 Re

40 (b)

. /
-, ~--"",
:0
~.0

.
.0

"
" 30
" ""
,,'

/ ~/
-0 /
.
.
'" 20 /
/

.
-0

:2 /'
/
/ "
/

/
~10 /
/
/
,,' "
,,

0 20 40 100 200 230 260 290 300 Re

FIGURE 17. The width (a) and length (b) of the wake bubble. The length is measured
from the centre of the cylinder.

showed a change ill Cn from 7.87 to 6.92 (eITors about 19 % and 5 %) when a free-
stream boundary condition was applied at r = 23.1 and r = 91.5 respectively. The
results ill closest agreement with ours are those obtained by Dennis (1976) for Re = 20
and 40 with use of his improved boundary conditions (Cn 1.998 and Cn = 1.494 ~

respectively).
There are several questions pertinent to the flow field which we did not specifically
address. For example. there is a possibility that non-unique. steady, perhaps un-
symmetric, solutions for high Reynolds numbers may exist. We did not detect any
evidence of simple bifurcations or tumillg poillts on the main solution branch up to
848 B. Fornherg

8
t 1:6'"
'" 1
- 1

8
'" t 1 t 1 1:6
'"
]'
00
8 ~1
'" t 1% 1% 'tI
§
- '" Q

~8 - J
- I~ I '" - -0"''''
--
~~~j
1
.,g

...'?00
0'"
1 t I~
0
~I
I
P:<

J i
~~00
,. -
",,:,
I: 1I
""
'. ~.g
"
~0
,. '" ,Q
0
O-oeDlQ
"ItI~~CO?
,. 00
1':''1 fi
"''''''' "'''' "" ~...
I
0

... _ "8
Q
00
0'"'" I~~ § ...0
"'':'
" "" I t .&'
f .s
.
.;;:
~00
cqr--QOOO
... ,.0000 1 ".
.;
00 00 ,.
~cN~ '" ~~tJ
l~Q Eo<

+-
 Steady flow past a circular cylinder 849

~=~_ ~=~, ~=o ~=o ~=~,- -- - -- - -

n=1 n=1 n=O n=O

FIGURE 18. The integration path for the pressure calculation.

1-4 (a)
x,
50 1-2
,
\
0. ,
0 ,
~1-0 'x
~., \
\
~0-8 \
,

"0
.::
.. 0-6
X,
'X,
-X
--X--X-X -limit 0,5-
.
0:
~~0-4

0-2

0
2 4 10 20 40 100 200 400 1000
Re

-1-4 (b)

-1-2
"0.
-0
0
:8 -1-0
~0
~~-0-8
"
~~..
-0-6

.
0:
~~-0-4

-0-2

0
2 4 10 20 40 100 200 400 1000
Re
FIGURE19. The pressure at (a) the front stagnation point and (b) the rear stagnation point.
850 B. F<m1herg

Pressure at front stagnation point

Pressure at
Method Probable rear
Reynolde correct stagnation
nwnber 1 2 3 value point
20 0-649 0-604 0-640 0-64 -0-27
40 0-586 0-562 0-571 0-57 -0-23
100 0-543 0-528 0-523 0-53 -0-17
200 0-490 0-514 0-506 0-51 -0-12
300 0-415 0-508 0-499 0-51 -0-09

TABLE 4. Pressure at stagnation points calculated by different methods.

Pressure

1-4
Re = 2

1-2

1.0 Re=4

0-8
Re=10
Re = 20
0-6 Re=40
Re = 100
Re =200.
300
0-4

0-2

1-5 1-0 0-5 o 0

o
0-5 0-4 0-3 0-2 0-1 o Re"
= 300
Re = 200
-0-2 Re = 100
'Re= 40
Re = 20
Re=10
-0-4

Re =4
-0-6

-0-8 Re= 2

-1-0
FIGURE 20. The distribution of pressure over the body surface.
Steady flow past a circular cylinder 851

0 0
'"
0
'"
I I I I I
'"
'"

0 -
0
I I I I I '"
'"
'"

...
8 - -
I I I
'" '"
'"
'" '"

0
0 '"
'" '" '"
~'" '" '"
'" '" '"

0
0 ..
.. '"
I I I I I
'"

~'"
...
[
'" ..
0 '" 0
'"
'"
'" I I '"
I I
.9
'" '"
~~.. ~m
......'" -
~0 0'"'" ...... 0
..
666
"''''''' "''''
66 8
I'<
.0

~0
~-
~'"
"''''
"''" <
'" Eo<
'" '" 66
"'''' I I

'"

..
0 0.....
...... '"
"'..
'" "'''''''
666 00
1"''''
 852 B.Fornherg

2.2
...... ' Re = 20
.................. .......
:.:... __ _ ___ _ _ _ __~...,
2.0 '-'''''-':-''-,,"'\''\..A,~~,

1-8

1.6 ........ ......... ............ Re =40

........... , ,,--,,--,-,,---,-"',.....-
- - - -
'.' .'
1-4

1.2
........ Re = 100
............. .................
rJ ....... ..'
.h.................. - -- - ------
1.0
, ,..: Re = 200
0.8
.-- .
':-_ _/~I"': - ", .:" - - - -- - - --
Re= 300 ~- _
,," -- - - - --- - --, /...~'~\~~~ / {/"-- -
0.6 '," I I
,
I! ,
0.4 "

0.2

1.21.5 3 4 5
~ 1,0 2,0 30 40 ~O 100 200 300'400500600
1.0 0.1 0.2 0.3 0,4 0,5 0,6 0,7 0,8 0.9 1.0 H 1.2 1.3 1.4 1.5 l-6 1.7 1.8 1.9 2.0 2.1 t
FIGURE 21. The computed drag coefficients On as functions of Reynolda number and radius of
integration path. coarse grid; ---, fine grid; -, extrapolation.

Be = 300. Any such singularities would most likely have led to a breakdown or slowing
up of the inner-<mter iteration or a singularity of the Jacobian matrix in the Newton
method. We monitored the determinant of this matrix and figure 23 shows the log-
arithm of it. Although we did not find any evidence of non-uniqueness we cannot
exclude the possibility of entirely different classes of solutions, even at Reynolds
numbers lower than Be = 300.

This research was supported by Control Data Corporation and by D.O.E. (Office of
Basic Energy Sciences).
Steady flow past a circular cylinder 853
854 B. Fornherg

I
,
2
\
\
\
\,
,
" "-
"-
"-
G'l '- - -- ---- x__x__
x---JtX1{_

o 20 40 100 200 230 260 290 300 Re

FIGURE 22.The drag coefficient 8B a. function of Reynolds number.

24000
23800
23600
23400
~23200
~23000
'"
-.!.-22800
0
;;;
.2 22600

22400

22200
22000
21800
0

FIGURE 23. The determinant of the Jacobian matrix for different Reynolds numbers.

REFERENCES
ALLEN, D. N. DE G. & SOUTHWELL, R. V. 1955 Relaxation methods applied to determine the
motion, in two dimensions, of a viscous fluid pasta. fixed cylinder. Quart. J. Meeh. Appl.
Math. 8, 129.
BATCHELOR, G. K. 1956 A proposal concerning laminar wakes behind bluff bodies at large
Reynolds number. J. Fluid Meck. 1, 388.
BRODETSKY. S. 1923 Discontinuous fluid motion past circular and elliptic cylinders. Proc. Roy.
Soc. A 102, 542.
 Steady flow past a circular cylinder 855
DENNIS,S. C. R. 1973 The numerical solution of the vorticity transport equation. Proc. 3rd
Int. Conf. on Numerical Method8 in Fluid Meek., vol. 2 (ed. H. Cabannes & R. Tewam),
Lecture notes in Physics, vol. 19, p. 120. Springer.
DENNIS, S. C. R. 1976 A numerical method for calculating steady flow past a cylinder. Proc.
5th Int. Conf. on Numerical Metluxl8 in Fluid Dynamics (ed. A.1. van de Vooren & P. J.
Zandbergen), Lecture notes in Physics, vol. 59, p. 165. Springer.
DENNIS, S. C. R. & CHANG,G. Z. 1970 Numerical solutions for steady flow past a circular cylin.
der at Reynolds numbers up to 100. J. Fluid Meek. 42,471.
GUSHCHIN,V. A. & SCHENNIKOV,V. V. 1974 A numerical method of solving the Navier-Stok.8S
equations. Zh. vychist. Mat. mat. Fiz. 14, 512.
HAMIELEc,A. E. & RAAL.J. D. 1969 Numerical studies of viscous flow aroundcircular cylin-
ders. PhY8. Fluids 12, 11.
HELMHOLTZ,
H. v. 1868 Dber discontinuirliche FliissigkeitB-Bewegungen, Phil. Mag. 36 (4),
337.
lMAI, I. 1951 On the asymptotic behaviour of viscous fluid flow at 8 great distance from a
cylindrical body, with special reference to Filon's paradox. Proc. Roy. Soc. A 208.487.
INGHAM,D. B. 1968 Note on the numerical solution for WlSteady viscous flow past a circular
cylinder. J. Fluid Meek. 31, 815.
JAIN, P. C. & KAWAGUTI,M. 1966 Numerical study of a viscous fluid past a circular cylinder.
J. PhY8. Soc. Japan 21, 2055.
JAIN, P. C. & SANKARARAO, K. 1969 Numerical solution ofunstea.dy viscous incompressible
fluid flow paat a circular cylinder. Phys. Fluids Suppl. II 57.
KAWAGUTI,M. 1953 Numerical solution of the Navier-Stokes equations for the flow around a
circular cylinder at Reynolds number 40. J. Phya. Soc. Japan 8,747.
KELLER. H. B. & TAKAMI, H.. 1966 Numerical studies of viscous flow about cylinders. In
Numerical Solutions of Nonlinear Differential Equations (ed. D. Greenspan), p. 115. Wiley.
KmCHHOFF, G. 1869 Zur Theorie freier Fhissigkeits-Strahlen. J. reine angew. Math. 70, 289.
LEONARD,B. P. 1979 A stable and accurate convective modelling procedure based on quadratic
npstream interpolation. Compo Meth. Appl. Meek. Engng 19, 59.
NIEUWSTADT,F. & KELLER, H. B. 1973 Viscous flow past circular cylinders. Compo Fluids. I,
59.
PATEL, V. A. 1976 Time-dependent solutions of the viscous incompressible flow past a circular
cylinder by the method of series truncation. Compo Fluids 4, 13.
PAYNE, R. B. 1958 Calculations of unsteady viscous flow past a circular cylinder. J. Fluid
Meeh. 4, 81.
ROACHE,P. J. 1976 Oomputational Fluid Dynamics. Albuquerque: Hennosa.
SMITH, F. T. 1979 Laminar flow of an incompressible fluid past a bluff body: the separation,
reattachment, eddy properties and drag. J. Fluid Meek. 92, 171.
SON, J. S. & HANRATTY,T. J. 1969 Nwnerica1 solution for the flow aroWld a cylinder at Rey-
nolds numbers of 40, 200 and 500. J. Fluid Meek. 35, 369.
TA, P. L. 1975 Etude nmmkique de l'ecoulement d'Wl fluide visqueux incompressible autour
d'un cylindre fixe ou en rotation. Effet Magnus. J. Moo. 14, 109.
TAKAMI,H. & KELLER, H. B. 1969 Steady two-dimensio.nal viscous flow of an incompressible
fluid past a circular cylinder. Pkys. Fluids Suppl. II51. .
TROM, A. 1933 The flow past circular cylinders at low speeds. Proc. Roy. Soc. A 141, 651.
THOMAN,D. C. & SZEWCZYK,A. A. 1969 Time dependent viscous flow over 8 circular cylinder.
PhY8. Fluids Suppl. II76.
TuANN, S.-Y. & OLSON,M. D. 1978 Numerical studies of the flow aroWld a circular cylinder by
a finite element method. Oomp. Fluids 6, 219.
UNDERWOOD,R. L. 1969 Calculations of incompressible flow past a circular cylinder at mod.
erate Reynolds numbers. J. Fluid Meek. 37,95.