Microsoft Word - MATH991101.final2.doc


Science and Technology, 5 (2000) 85-104 
© 2000 Sultan Qaboos University 
 

85 

Flow Past A Rotating Circular Cylinder 
and A Rotlet Using the Finite-Difference 
Method 

T. B. A. El Bashir 

Department of Mathematics and Statistics, College of Science, Sultan 
Qaboos University, P.O. Box 36, Al Khod 123, Muscat,  Sultanate of Oman. 

 دائرية دوارة ودوامة باستخدام طريقة الفرق المحدودة االنسياب على اسطوانة

 طيفور البشير

باستخدام بعض التحويالت   ) دوارة  (  يختص هذا البحث بدراسة االنسياب المتولد من دوامة في وجود اسطوانة دائرية              :خالصة  
هذا وقد تم التوصل إلى حل عددي يتطابق بشكل ممتاز . ودةعلـى المعادالت األساسية الحاكمة للمسألة مع طريقة الفروق المحد  

   .مع النتائج التحليلية المتاحة ، ثم تم استخدام الطريقة نفسها للحصول على نتائج جديدة لمسألة ليس لها حل تحليلي 
 
ABSTRACT : In this paper the flow generated by a rotlet in the presence of a circular cylinder is considered.  We 
introduce a transformation which simplifies the equations and boundary conditions. We use  the finite-difference 
method to obtain results in excellent agreement with all the available analytical results.  Results are presented for 
Reynolds numbers, based on the diameter of the cylinder, in the range 20Re0 ≤≤  and the rotational 
parameter, α , in the range 30 ≤≤ α  and strength of rotlet, β , in the range 30 ≤≤ β .  The results are 
found to be applicable over a wide range of values of α  and β .  The calculated values of the drag, lift and 
moment coefficients and the general nature of the streamline patterns are in good agreement with analytical 
results .  The method is then utilized to obtain new results for which no analytical solution is possible. 

      
KEYWORDS: Slow flow, cylinder, rotlet.  

    
 
                    study has been made of  the flow generated by rotating a  circular  cylinder  within a  uniform   

             stream of viscous fluid in the presence of a line rotlet. With the origin of the coordinates 
coinciding with the centre of the cylinder, the polar coordinate system ( )ϑ,r  has the boundary of the 
cylinder at r = a  and the position of the rotlet at ( )2/3, π∗c , as shown in Figure 1. There are four 
basic parameters that occur in this problem, namely the Reynolds number, defined as a2Re = vU / , 
the rotational parameter Ua /0ωα = , the non-dimensional length acc /

∗=  and the non-
dimensional strength of the rotlet Ua/Γ=β , where v  is the coefficient of kinematic viscosity of 
the fluid, U  the uperturbed main stream speed (the stream at infinity is assumed to flow parallel to 
the x -axis in the positive x direction), 0ω the angular velocity of the cylinder and Γ  the strength of 
the rotlet. At zero Reynolds number the governing equation is the biharmonic equation. In the 
absence of rotlet no solution of this equation, which satisfies both the boundary conditions on the 
cylinder and at infinity, is possible. This arises because it is necessary to maintain in the solution a 
term of the form ( ) ( )ϑsinln rr  in order for both the velocity components on the rotating cylinder to 
be satisfied. However, although a solution retaining such a term is obviously valid at points not too 
far away from the cylinder such a solution deteriorates as one moves further away. As such this 
solution fails to satisfy the boundary condition at infinity, with the exact multiple of the unwanted 
term at infinity remaining undetermined. This unknown constant can be established by treating the 



T. B. A. EL BASHIR 
 

 86

solution of the problem as the approximation to the inner flow past the circular cylinder as the 
Reynolds number tends to zero. Then matching with the solution as obtained from the outer region 
where the first approximation to the Navier Stokes equations are the Oseen equations. Full details 
regarding these expansions and matching procedure can be found in Proudman and Pearson (1957). 
Basically what this means is that the uniform flow past a circular cylinder is not a well- posed Stokes 
problem. In fact, at zero Reynolds number the problem is singular. 

 
                                                                            Y 

 

           U                               r 

                                                                                   ϑ                 X       

                                                                                                  

                                                                                                   (1,0) 

                                                                      *                  

                                                                                          ( )2/3,* πc                                                                            

Figure 1.  The geometry of the circular cylinder and the position of the rotlet 

However, the introduction of the rotlet , strength ∗=Γ Uc  so that c=β , into the flow field at 
any given distance ∗c  along the negative y -axis, allows a solution to be obtained. Analytical 
solutions to this problem have been obtained by Dorrepaal et al (1984). Their work examines the 
flows generated in a fluid by the introduction of a line singularity, such as a stokeslet or rotlet, in the 
presence of a circular cylinder and shows that a phenomenon analogous to the Stokes paradox exists 
in those flows with a uniform stream far from the cylinder. As a consequence, the uniform streaming 
flow past a circular cylinder, when a line stokeslet or rotlet of certain strength is present, is a well-
posed problem in Stokes flow. 

The solutions by Dorrepaal et al. (1984) employed an image type approach, plus a clever and 
simplistic deduction, which enabled the result to be constructed devoid of most of the analysis.  
However, the present work has established the same solution by using a Fourier Series approach and 
has confirmed this numerically by the application of a modification to the Boundary Element Method.  
The latter appears to provide an approach for the solution of the biharmonic equation, which requires 
only the position of the singularity, plus the physical values of the drag, lift and the moment on the 
circular cylinder to be known.  In addition, it seems capable of being extended to accommodate the 
presence of several different bodies as well as allowing more complex shapes for which an analytical 
solution is impossible.  It is intended to show that the presence of a rotlet in a uniform flow at non-
zero Reynolds number allows an otherwise singular problem to become well-posed as the Reynolds 
number becomes zero. 

The main aim of this paper is to solve numerically the Navier-Stokes equations for steady, two-
dimensional, incompressible viscous fluid flow past a rotating circular cylinder of radius a in the 
presence  of a rotlet of strength Γ  which is located at the point ( ) ( ) accr >∗= ∗ ,2/3,, πϑ .  At large 
distances from the cylinder it is assumed that there is a uniform flow of speed U  which is parallel to 
the negative x -axis.  Initially the strength of a rotlet is set to zero and the problem solved with 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 87

Reynolds numbers 5 and 20.  The results are in very good agreement with those obtained by Dennis 
and Chang (1970) and Fornberg (1980).  Using this as an initial estimate of the solution when a line 
singularity is present an iterative technique is developed in order to solve the problem when a rotlet, 
at ( )2/3, π∗c , of small strength is introduced into the flow.  As the drag, lift and moment on the 
circular cylinder are the most important physical quantities, as well as being easy to measure 
experimentally, see Fornberg (1980), particular attention has been paid to these quantities in this 
work. 

Having established the numerical procedure the Reynolds number is decreased towards zero 
whilst continuing to solve over a range of non-dimensional strengths of the rotlet, namely 

5.30 ≤≤ β . The question that needs to be answered is whether as Re approaches zero the value of 
the parameter β  would tend to that unique value obtained by Dorrepaal et al. (1984) in their 
analytical solution. 

Basic Equations and Boundary Conditions 

The origin of the coordinate system is fixed at the centre of the circular cylinder of radius a   
and the positive x -axis taken in the same direction as that of the uniform flow at large distances from 
the cylinder.  Polar coordinates ( )ϑ,r  are chosen such that 0=ϑ  coincides with the positive x -axis, 

 
 

( )ϑcosrx =  and  ( )ϑsinry =                                                 (2.1) 
  

A line rotlet of strength Γ  is located at the point 2/3, πϑ == ∗cr , where ac >∗ .  The steady 
flow of an incompressible fluid in a fixed two-dimensional Cartesian frame of reference can be 
described by the equations, 
 

( ) ,p1. 2 uuu ∇+∇−=∇ v
ρ

                                                  (2.2) 

    
                0=⋅∇ u ,                                                                (2.3) 

 
where p,, ρu  and v  are the velocity, density, pressure and the kinematics viscosity of the fluid, 
respectively.  Applying the curl operator to the two-dimensional equation (2.2) produces 
  

 ( ) ,2ωωu ∇=∇⋅ v                                                    (2.4) 

where Χ∇=ω u . 
In two-dimensional motion the polar resolutes of  u  can be expressed in terms of the streamfunction 
Ψ by 
 
 

r
v

r
vr ∂

Ψ∂
−=

∂
Ψ∂

= ϑϑ
,

1
                                              (2.5) 

where rv  and ϑv  are the velocity components in the r  and ϑ  directions, respectively.  By 
introducing the dimensionless variables 

,/,/,/ aa UUuuχχ Ψ=Ψ′=′=′  and ,/ Uaωω =′                         (2.6) 
 
then the governing equations in non-dimensional form become 
 



T. B. A. EL BASHIR 
 

 88

 
( )
( )yx ′′∂

′Ψ′∂
−=′∇

,
,

Re2
ω

ω                                                       (2.7) 

 ,2 ω′−=Ψ′∇                                                    (2.8) 
 

where ω′  and Ψ′  are the non-dimensional scalar vorticity and streamfunction, respectively.  For 
convenience the accent will from now on be ignored.  It is required to solve equations (2.7) and (2.8) 
subject to the no-slip conditions imposed by the circular cylinder, namely 
 

α−=
∂
Ψ∂

=Ψ
r

,0 on ,20,1 πϑ <≤=r                                       (2.9) 

and the boundary conditions at large distances from the cylinder 
 

( ) ( )ϑ
ϑ

ϑ cos
1

,sin →
∂
Ψ∂

−→
∂
Ψ∂

rr
     as ,20, πϑ <≤∞→r                           (2.10) 

 
In the presence of a line rotlet the streamfunction behaves as 

 
 

≅Ψ 1ln Rβ−        as ,01 →R                                                (2.11) 
 

where 1R  measures the distance from the rotlet  and is thus given by 
 

( )( ) ,sin2 21221 ϑrccrR ++=                                                 (2.12) 
 

where ( )2/3, πc  is the position of the rotlet.  In the above definition of the streamfunction the signs 
appearing in the expressions in (2.5) are the opposite to those given by Dorrepaal et al. (1984) but 
follow those adopted by Fornberg (1980) since it is a comparison with their results at non-zero Reynolds 
number that is to be undertaken.  For numerical convenience the perturbation streamfunction Ψ is 
introduced as 

 
 ,Vy βψ −−Ψ=                                                             (2.13) 

 

where ( )( ) 2122 sin2ln ϑrccrV ++=  with *cc = a/  and ( )Ua/Γ=β  being two non-dimensional 
parameters.  Expansion (2.13) has been taken so that 0→ψ  as ∞→r .  If the parameter 0=β  the 
problem reduces to that solved by Fornberg (1980).  However, with the Reynolds number equal to zero 
and the parameter c=β  the situation is that studied by Dorrepaal et al. (1984) except that the 
geometry in the present case corresponds to a rotation through 2/π  of their flow pattern.  Hence, their 
stream is flowing along the negative y -axis with their rotlet at ( )0,*c , whereas in the present geometry 
the stream flows along the negative x -axis with the rotlet at ( )2/3,* πc .  Using expression (2.13) in 
equations (2.7) and (2.8) gives 

 

( ) ( )
( )

( ) ( )
( ) 




















++
+

∂
∂

+





















++
+

+
∂
∂

−







∂
∂

∂
∂

−
∂
∂

∂
∂

−=
∂
∂

+
∂
∂

+
∂
∂

ϑ
ϑ

βϑ
ω

ϑ
ϑ

βϑ
ϑ
ωω

ϑ
ϕ

ϑ
ωϕ

ϑ
ωωω

sin2
cos

cos
2
Re

sin2
sin

sin
2
Re

2
Re11

22

2222

2

2

2

rccr
rc

r
rr

rccr
cr

rrrrrrrr             (2.14) 

 
 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 89

 
 
respectively.  The boundary conditions (2.9) and (2.10) are then expressed in the form 
 

( ) ( )( )2
1

sin1lnsin 2 ϑβϑψ cc ++−−=   on πϑ 20,1 <<=r ,                          (2. 16) 
 

( ) ( )
( )









++

+
−−−=

∂
∂

ϑ
ϑ

βϑα
ψ

sin21
sin1

sin
2 cc

c
r

  on πϑ 20,1 <<=r ,                   (2.17) 

 

0
1

→
∂
∂

=
∂
∂

ϑ
ϕϕ

rr
  as πϑ 20, <≤∞→r .                                            (2.18) 

 
Filon (1926) showed that in the absence of any rotlet the asymptotic form for the dimensional 
streamfunction at large distances from the cylinder and outside the wake region is given by 

 
 

ψ ~ ( ) ( ) ( )
π

πϑ
π

ϑ
22

/ln
sin

−
++ DL

CarC
r ,       as     πϑ 20, <<∞→r ,                     (2.19) 

 
where ( )aULC L 2/ ρ=  and ( )aUDC D 2/ ρ= , L  and D  being the lift and drag on the cylinder.   
Imai (1951) found higher–order terms in this streamsfunction expansion and showed how the 
coefficients relate to the moment on the cylinder.  In the case of zero Reynolds number no solution of 
the equation 04 =Ψ∇  which matches the free stream condition at infinity and satisfies the boundary 
condition on 1=r  is possible.  The solution, which satisfies the no slip condition on the cylinder and 
tends to infinity most slowly as ∞→r  is  

 
ψ ~ ( ) ( ) ( )[ ]rrrrA 2/12/lnsin +−ϑ                                              (2.20) 

 
This has been obtained by discarding the term  involving 3r .  The non-dimensional drag is directly 
related to the coefficient A, via the expression Aπ4 , but the solution suffers from the defect that it does 
not determine the value of the constant A.  The neglected inertial terms are of the order ( ) ( )( ) 22 /ln rrA  
whilst the viscous forces are of the order ( )3Re/ rA  and these terms are of comparable order when 

( )( ) ≅arrA /lnRe ( )10 . Hence, the Stokes solution should not be expected to be valid beyond  a value 
given by this expression.  That is why the Stokes solution may be an adequate representation of the 
fluid flow relatively close to the cylinder but cannot represent a uniform approximation to the total 
velocity distribution.  However, it is possible to write the Stokes solution in the form  

 
( )( ) ( )( )( ) ( )[ ] ( )ϑψ sin2/12/RelnReln rrrfrfrA +−+−= ,                             (2.21) 

 
where ( )Ref  is an arbitrary function of Re . For  ( ) 1Re <<f  and ( )Rerf  of order unity, the 
dominant term is ( )( )( ) ( )ϑsinReln rfA− .  If this is to represent the external flow, namely a uniform 
stream ( )ϑsinUr , then one must set ( )( )Reln/1 fA −= .  By substituting this value of A  and 

( )Re/1 fr = into ( )( ) arrA /lnRe , one obtains 
 

 ( )ReRe/ f ~ ( )10 .                                                       (2.22) 

,
11

2

2

22

2

ω
ϑ
ψψψ

−=
∂
∂

+
∂
∂

+
∂
∂

rrrr
(2.15)



T. B. A. EL BASHIR 
 

 90

So for ( ) ReRe =f  the Stokes solution leads to a uniform stream of order unity in that region where 
the Stokes equation ceases to be valid.  This suggests that for small Reynolds numbers the external 
uniform stream condition is reached before the Stokes approximation breaks down,  so that the Stokes 
flow represents the solution in the inner region near the cylinder but the outer flow requires the 
introduction on an Oseen variable and the procedure involves the matching of the inner and outer 
expansions.  This matching process leads to the obvious linkage between the coefficients of the terms in 
Oseen’s expansion in the outer region, closely related to Filon’s expansion with those from the Stokes 
solution in the inner region.  Hence, the appropriate value of A can be obtained, as already indicated 
from the solution in the outer region.   Any similar term in a Stokes expansion, such as ( ) ( )ϑcosln rBr  
term in the streamfunction expansion, has its coefficient similarly related to the lift.  Solutions of the 
Stokes equation represented by the form ( ) ( )ϑsinrf  or ( ) ( )ϑsinrg  fail to produce any moment 
contribution, but it can easily be seen that the solution of 04 =Ψ∇  which has a non-zero moment 
arises from the ( )ar /ln  term.  Using these results, it is possible from the Stokes expansion far from the 
cylinder to establish both the force and the moment on the cylinder. The Stokes solution produced by 
Dorrepaaal et al (1984) is able to immediately produce the drag, lift and moment on the circular 
cylinder from its expansion at large distances from the cylinder;  although the contribution from the 
singularity at the rotlet must first be removed from the  coefficient of the ( )ar /ln  term before it 
represents the moment on the cylinder. 

Since the streamfunction and vorticity equations are both elliptical in nature we should supply one 
condition for each of the variables rather than use directly conditions (2.16), (2.17) and (2.18).  It has 
been reported by Fornberg (1980, 1985), that the choice of the boundary condition for the vorticity, ω ,  
is not as sensitive as that for the streamfunction, ψ , and many authors have paid particular attention to 
the boundary condition for Ψ  at large distances from the cylinder. 

 
 

        

                                 1=ξ  /, BB ff                                                               

                                    )4,()0,( ξξ ff =                                       )0,()4,( ξξ ff =                                             
                                    )4,()0,( ξωξω =                      )0,()4,( ξωξω =  

                              0=ξ  
                                                           0=n     0=f  , 0=ω         4=η   

 

Figure 2.  The geometry of the solution domain. 

In this work a new numerical technique is introduced in order to avoid having to enforce the boundary 
condition (2.19) at large distances from the cylinder.  We introduce the transformation 
 

π
ϑ

ηξ
2

,
1

==
r

                                                         (2.23) 

and 
 ( ) ( ) rrrf /,, ϑψϑ = ,                                                      (2.24) 

where 
 ( ) ( )ηξξψηξ ,, =f .                                                      (2.25) 

 
Thus, we have 0=f  on 0=ξ  (i.e. at ∞=r ) and this requires no approximation to be made for ψ  
at large distances from the cylinder.  With the transformation (2.23) the solution domain 
( )πϑ 20,1 <≤∞<≤ r  is transformed into a finite rectangular region in the ( )ηξ,  plane 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 91

( )40,10 <≤<≤ ηξ , see Figure 2.  Substituting expressions (2.23) and (2.24) into the governing 
equations (2.14) and (2.15) one obtains 

 
( )
( )

( ) ( )

( )
( )

( )
( )

,0
4

Re
2

2/sin21
2/cos

2/sin21
2/sin12

2/sin
2

2/cos
2

,
,2

2

2

2

2
3

2

2
4

22

4

22

2

2
2

=

















∂
∂

+
∂
∂

+
∂
∂

+

∂
∂









++

+
∂
∂









++

+
+

∂
∂

+
∂
∂

+
∂
∂

+
∂
∂−

η
ω

π
ξ

ξ
ω

ξ
ξ
ω

ξ

ξ
ω

πηξξ
πηβξ

η
ω

πηξξ
πηξ

π
βξ

πη
η
ω

π
ξ

πη
ξ
ω

ξ
η
ω

π
ξ

ηξ
ω

π
ξ

cc
c

cc
c

f
f

              (2.26)                   

 
 

ωξ
ηπ

ξ
ξ

ξ
ξ

ξ −=+
∂

∂
+

∂
∂

−
∂

∂
f

fff
2

2

2
2

2

2
3 4                                     (2.27) 

 
 

and the boundary conditions (2.10) and (2.11) now become 
 

( ) ( )( ) 40,12/sin21ln2/sin 212 <≤==++−−= ηπηβπη ronfccf B                             (2.28) 

( ) ( )
( )

,40,1'
2/sin21

2/sin1
2/sin

2
<≤==








++
+

++=
∂
∂

η
πη

πη
βπηα

ξ
ronf

cc
cf

B                              (2.29) 

40,0,0,0 <≤=== ηξω onf                                             (2.30) 
 

Further, since the solution is periodic in  η  we also require that  
 

( ) ( ) ( ) ( ) ,10,0,4,,0,4, ≤≤== ξξωξωξξ forff                               (2.31) 
 

which along with all the other boundary conditions are indicated in Figure 2. 

The Solution Technique 

In order to obtain numerical solutions of the equations (2.26) and (2.27) subject to the 
boundary conditions (2.28), (2.29), (2.30) and (2.31) the region of integration 40,10 <≤≤≤ ηξ  
is covered by a rectangular mesh system of size h and k in the η  and ξ directions, respectively, and 
a modified finite-difference approximation to the differential equations (2.26) and (2.27) is 
employed.  If the standard central-difference approximation is employed, then it is not always 
possible to obtain a convergent solution especially at moderate and large values of the Reynolds 
number.   A number of schemes have been introduced with the object of improving the efficiency 
of the convergence of such schemes, e.g. upwind and downwind differencing as discussed by 
Greenspan (1968).  Gosman et al. (1969) and Runchal et al. (1969). These efficiencies arise from 
the fact that the difference equations are associated with matrices, which are diagonally dominant 
and thus amenable to iterative methods of solution.  However, these schemes suffer from a 
deficiency in that  they are only of first-order accuracy since forward or backward differences are 
employed to approximate first derivatives rather than the more accurate second-order central-
difference formulae.  The formulation in which all derivatives are approximated by central 
differences is of second-order accuracy but the matrix associated with the difference equations may 
not be diagonally dominant.  The iterative procedures may be slowly convergent or even divergent 
when this method is applied. 

There are also finite-difference schemes, which are of second-order accuracy and for which 
the associated matrices are always diagonally dominant.  These methods rely upon rather 
specialized forms of local approximations and yield  difference equations, which involve the 



T. B. A. EL BASHIR 
 

 92

exponential function.  These schemes were first introduced by Allen and Southwell (1955) in 
approximating the equation governing the vorticity during the course of finding numerical solutions 
for the steady two-dimensional flow past a circular cylinder.  Dennis and Hudson  (1978) showed 
that by a suitable adaptation of an alternative to the Allen and Southwell method suggested by 
Dennis (1960), an approximation of second-order accuracy, yielding difference equations with an 
associated matrix which is diagonally dominant, can be obtained.  These difference equations do 
not involve the exponential function and can be looked upon as a rather more complicated version 
of the central-difference formulation.  The Dennis and Hudson method contains more terms in the 
finite-difference equations than the usual central-difference approximation but the presence of these 
extra terms is very important for the associated matrix to be diagonally dominant.  Numerous 
authors have performed several numerical experiments which confirm that the method by Dennis 
and Hudson succeeds where the standard central-difference formulation fails, see for example 
Dennis (1960) and Dennis and Hudson (1978). Thus, in this paper a modified version of the finite-
difference approximation as described by Hudson and Dennis, has been used. 

It is found most convenient to set up a mesh system such that the mesh size in both the ξ  and 
η  directions are Nhk /1== , where N is a pre-assigned positive integer.  In view of the periodic 
conditions (2.31), and extra line of computation h+= 4η  for 10 ≤≤ ξ  is introduced.  Then we 
have ( ) ( )11 ++ MXN  mesh points, where 14 += NM , the mesh points 
( ) ( )MjNiii ≤≤≤≤ 0,0,ηξ  are ( )jhih, .  If subscripts 0,1,2,3 and 4 denote quantities at the 
grid points ( ) ( ) ( ),,,,,, 000000 ηξηξηξ hh +− ( )h+00 ,ηξ and ( )00 ,ηξ h−  respectively, then 
on using the Dennis and Hudson scheme, equations (2.26) and (2.27) may be written in the form 

03
0

2

4
0

**

32
0

22
0

*2*

12
0

202
0

2

2
0

2

2*

484
1

84
1

4
2

2
ω

ξξ
ββ

ξπξ
ββ

ξπξξπ
β h

f
h

ff
h

ff
h

+







+++








−+=








−+         (3.1) 

( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) 40000
2*

30000
0

2

20000

2*

100002
0

202
0

2

2*

,,
4

,,
1

,,
4

,,
12

2

ωηξηξ
β

ωηξηξ
ξπ

ωηξηξ
β

ωηξηξ
ξπ

ω
ξπ

β









−−+








+++









+++








−−=








+

EbDa

EbDa
     (3.2) 

 with 
k
h

=*β , and 

       ( ) ( )







++= ϑ
ξξ

υ
πξ

ηξ sin
1

8
Re

,
00

0
02

0
00

fh
a                                       (3.3) 

        ( ) ( )







++= ϑ

ξξπξ
βηξ cos

16
Re

88
Re

,
00

02
0

00
hh

u
h

b                                 (3.4)                                

        ( )
( )









++

+
−=

ϑξξ
ϑξ

πξ
β

ηξ
sin21

sin1
8

Re
,

0
2
0

2
0

2
0

00 cc
ch

D                                 (3.5) 

       ( ) ( )
( )









++

−=
ϑξξ

ϑβ
ηξ

sin21
cos

16
Re

,
0

2
0

200 cc
ch

E                                 (3.6) 

where ( )0,0 υu  are defined as  



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 93

( ) ( )
.

,
,

, 00
0

00
0 ξ

ηξ
υ

η
ηξ

∂
∂

−=
∂

∂
=

ff
u                                                 (3.7) 

 
The standard central-difference approximations may be obtained by setting the extra terms 
( ) ( )0000 ,,, ηξηξ ED  to be zero. 

We now briefly outline how the boundary conditions (2.28)-(2.30) can be implemented.  The boundary 
condition for the vorticity on 1=ξ  can be found by using Taylor expansion for f and ω  and inserting 
them in equation  (2.27) to get second-order accurate finite differences. 

Boundary conditions for ω : 

On ;0;0,0:40,0 0 =<≤=<≤= jMji ωηξ                              (3.8)                         

                         on ;0;0,:40,1 njjMjNi ωωηξ =<≤=<≤=                           (3.8) 

( )











+


















 





 +−+










−−−−−










−−

=
−−

hh

g
h

ghfhhhhffhh ijjnjijnijnnj

nj

1
3
2

31
26

34'
6
2

2
1

6
2

6
3

2
2

1 2
π

ω

ω      (3.9)                        

with 

( )( )
( ) 














+++
++

+
−= j

j

j
nj cccc

c
f η

πβ
η

ηβ
α

π

π

2
sin21ln

2sin21
sin1

' 2
2

2
2                             (3.10) 

( )( )
( )( ) 














++

++
−








+=
2

2
2

22
2

22

sin21

21sin
42

sin
4

j

j
jij

cc

ccc
g

η

ηβπ
η

ππ
π

π

                               (3.11) 

( ) ( ) ( )( )
( )( )322

2
23

2
232

2
sin212

cos2sin313

j

jj
j

cc

ccccc
g

η

ηηβπ
π

ππ

++

−+++
−=                          (3.12) 

on ;;,0:10,4 1iimMjNih ωωξη ==≤≤≤≤+=                              (3.13) 

on 

;;0,0:10,0 10 −==≤≤≤≤= imijNi ϖωξη                                     (3.14) 

Boundary conditions for f: 

On ;0;0,0:40,0 0 ==<≤<≤= jfiMjηξ                             (3.15)                         

                                       on ;,0:40,1 NiMj =<≤<≤= ηξ              

    















++







−= njccnjf nj 2
sin21ln

22
sin 2

πβπ
                                 (3.16) 

on ;;,0:10,4 1iim ffMjNih ==≤≤≤≤+= ξη                            (3.17) 



T. B. A. EL BASHIR 
 

 94

on ;;0,0:10,0 10 −==≤≤≤≤= imi ffjNiξη                            (3.18) 

The resulting finite-difference equations (3.1) and (3.2), subject to the boundary conditions (3.8) – 
(3.18), were solved iteratively.  From equation (2.6) the non-dimensional speed of the uniform stream at 
large distances from the cylinder is unity, whereas the non-dimensional strength of the rotlet at the 
position ( )2/3, πc , namely β , is variable.  In order to compare the results with Dorrepaal et al. (1984) 
the value of c is fixed at the magnitude used in their calculations, namely 3=c . 

Since the force components (drag and lift) and the moment are very sensitive to the method of solution, 
particular attention has been given to these quantities.  If xF  and yF  are the dimensional drag and lift 
on the cylinder, then 

( ) ( )( ) ϑϑσϑσπ ϑ rdF rrrrx ∫ =−=
2
0 1

sincos                                      (3.19) 

( ) ( )( )∫ =+−=
π

ϑ ϑϑσϑσ
2
0 1

sincos rdF rrrry                                    (3.20) 

Introducing the constitutive relations 

,2
r

p rrr ∂
∂

+−=
υ

µσ                                                      (3.21) 

,
1









∂
∂

+







∂
∂

==
ϑ
υυ

µσσ ϑϑϑ
r

rr rrr
r                                       (3.22) 

      ,
1

2 







+
∂
∂

+−=
rr

p r
υ

ϑ
υ

µσ ϑϑϑ                                            (3.23) 

into expressions (3.19) and (3.20) results in 

( ) ,sin12 2
2

22
0

ϑϑµ
ϑ

µ
ϑ

π
d

rr
p

rFx











Ψ∇+











Ψ+

∂

Ψ∂
∂
∂

−
∂
∂

= ∫                             (3.24) 

( ) ,cos12 2
2

22
0

ϑϑµ
ϑ

µ
ϑ

π
d

rr
p

rFy











Ψ∇+











Ψ+

∂

Ψ∂
∂
∂

−
∂
∂

−= ∫                           (3.25)                         

 
By using the ϑ̂  component of the Navier Stokes Equations on the surface of the cylinder where the 
component of velocity normal to the surface is zero, namely         

( )Ψ∇
∂
∂

−=






∂
∂ 21

r
p

r
µ

ϑ
                                                     (3.26) 

and substituting  ω−=Ψ∇2 , which is equation (2.8) but in its dimensional form, the equations (3.24) 
and (3.25) can be written in the form 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 95

( ) ,sin22
0 2

2
ϑϑ

ϑ
ω

ω
µ

π
d

rrr
rrFx ∫























Ψ+

∂

Ψ∂
∂
∂

−−
∂
∂

=                            (3.27) 

( )∫






















Ψ+

∂

Ψ∂
∂
∂

−−
∂
∂

−=
π

ϑϑ
ϑ

ω
ω

µ
2
0 2

2
.cos

2
d

rrr
rrFy                           (3.28) 

 
The above expressions are still dimensional, but defining the lift and the drag coefficients by 

( )aUFC xD 2/ ρ=   and ( )aUFC yL 2/ ρ=                                   (3.29) 
enables DC  and LC  to be written as 

( ) ,sin2
Re
2

2

22
0

ϑϑ
ϑ

ω
ωπ

d
rrr

rrC D






















Ψ+

∂

Ψ∂
∂
∂

−−
∂
∂

= ∫                          (3.30) 

( ) ,cos2
Re
2

2

22
0

ϑϑ
ϑ

ω
ωπ

d
rrr

rrC L






















Ψ+

∂

Ψ∂
∂
∂

−−
∂
∂

−= ∫                       (3.31) 

where Ψ,r  and ω  are all non-dimensional. 

In terms of the independent variables ζ  and DC,ϑ  and LC  become 

 

( ) ,sin2
Re
2

2

22
0

ϑϑ
ϑξ

ω
ξ
ωπ

dC D






















Ψ+

∂

Ψ∂
∂
∂

+−
∂
∂

−= ∫                            (3.32) 

( ) ,cos2
Re
2

2

22
0

ϑϑ
ϑξ

ω
ξ
ωπ

dC L






















Ψ+

∂

Ψ∂
∂
∂

+−
∂
∂

−−= ∫                         (3.33) 

 
When the boundary conditions on the cylinder are introduced, expressions (3.32) and (3.33) reduce to 
those given by Fornberg (1980), except for a reversal of the sign in LC  and the omission of a Factor 2 
in the definition of both  LC and DC  , which are both produced by the scaling we adopted.  

The moment on the cylinder 

     ϑσ ϑ
π

rdrM r∫=
2
0

                                                       (3.34) 

becomes, in dimensional form, on substituting the appropriate constitutive equation 

       ϑ
ϑ

µ
π

d
rrrr

rM












∂

Ψ∂
−

∂
Ψ∂

−
∂

Ψ∂
−= ∫ 2

2

22

22
0

2 11                                   (3.35)                         



T. B. A. EL BASHIR 
 

 96

Introducing the moment coefficient defined by  

         ( )22/ aUMCM ρ=                                                       (3.36) 
results in the non-dimensional expression 

.
22

Re
2

2

2

2
2
0

2 ϑ
ϑ

ω
π

d
rrr

rCM












∂

Ψ∂
+

∂
Ψ∂

+= ∫                                 (3.37) 

In terms of the independent variables ξ  and ϑ  the moment coefficient becomes 

,22
Re
2

2

22
0

ϑ
ϑξ

ω
π

dCM












∂

Ψ∂
+

∂
Ψ∂

−= ∫                                   (3.38) 

and the boundary conditions on the cylinder reduce this expression  to 

 [ ] .2
Re
2

0
2
0

ϑωω
π

dCM −= ∫                                          (3.39) 
 
Formulae (3.32), (3.33) and (3.39) are evaluated using Simpson’s rule.  Due to the need to evaluate the 
drag, lift and moment when 0Re =  we will work with the non-dimensional quantities 

2/Re,2/Re LD CC  and 2/Re MC  from now on, namely ( ) ( )UFUF yx µµ /,/  and ( )UaM µ/ .  
However, when the Reynolds number vanishes, that is the independence of the mainstream at infinity 
from the strength of the rotlet is no longer valid,  It is necessary to redefine the above three quantities as 
( ) ( ) ( ) ( )ΓΓ µµ /,/ aFaF yx and ( )Γµ/M  when discussing an iterative scheme in the following section. 

The Results 

Before we undertake a detailed discussion of the results it is important to stress from Dorrepaal et 
al. (1984), and from our own investigation, that at 0Re =  it is impossible to obtain a solution for an 
arbitrary value of the parameter β , which is contrary to what occurs at non-zero Reynolds numbers.  In 
the  0Re =  situation there is only one unique value of β  which will produce the uniform stream at 
infinity, and the difficulty numerically is to achieve some mechanism which will enable the numerical 
method to converge to this quantity.  Since Stokes flow with a uniform flow at infinity means that body 
must be free of any force, one can develop an iterative numerical scheme based upon this fact to enable 
the unique value of β  and the corresponding flow solution to be determined.  Exactly how this is 
achieved will be discussed in detail later within this section.  In order to compare the results obtained by 
the present numerical method with those obtained previously, the values of the Reynolds number 

5Re =  and 20Re = , in the absence of any rotlet ( )0=β  and with the rotational parameter 
1,5.0,0=α   and 2, were initially investigated.  Solutions are obtained for the mesh sizes 

20/1,15/1,10/1=h and 40/1 .  It is found that the contour lines for the streamfunction and 
vorticity distribution are in good agreement with most of the previous investigators (Dennis and Chang 
(1970); Fornberg (1980, 1985)), at both the above Reynolds numbers, repeated extrapolation of results 
derived from  2h  extrapolation of the results from 20/1,15/1,10/1=h  and 40/1   produced only a 
0.1% change in the coefficients of the lift and the drag when compared with the values derived directly 
from using 40/1=h .  As a consequence it is proposed not to implement extrapolation but instead to 
utilize a small value of h, namely 40/1=h .  Values of 2Re/  times the coefficients of the lift, drag and 
moment are presented in table 1.  The lift and drag coefficients are in good agreement with those 
produced by Fornberg (1980). 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 97

Table 1: The variation of the coefficients of lift, drag and moment with α  for 5Re = and 20 . 
 

  5Re =    20Re =   

α  2/Re LC  2/Re DC  2/Re MC  2/Re LC 2/Re DC  2/Re MC  

0.0 0.000 19.716 0.000 0.000 39.861 0.000 

0.5 6.939 19.565 0.616 25.635 3.430 2.354 

1.0 14.176 19.226 1.197 52.289 38.465 4.426 

2.0 29122 17.514 2.246 114.275 32.515 7.026 
 

Having established that the numerical procedure is producing the correct results for 0Re ≠  and 
0≠α , the condition 0=β  is relaxed.  However, since the interest is in whether the problem is well posed 

as the Reynolds number tends to zero, it is proposed to set α  to be zero in the remainder of the calculations 
to be presented in this paper. Such a restriction can easily be removed if required. 

It has already been established analytically by Dorrepaal et al. (1984) that at 0Re =  a solution in 
which there is a uniform stream of unit non-dimensional strength flowing parallel to the positive y -axis is 
possible provided that the rotlet is placed at the non-dimensional position ( )0,c , where the non-dimensional 
strength of the rotlet is c=β .  In the present situation the stream is flowing parallel to the positive x -axis 
and therefore the position of the rotlet has to be ( )2/3, πc  in order to produce a solution.  The question is 
whether the solutions at 0Re ≠  for a range of values of β , namely varying non-dimensional strengths of 
the rotlet, can be used to predict the unique analytical value of this parameter which produces the solution at 

0Re = .  The results of the solution by Dorrepaal et al. (1984) show that there is no force or moment on the 
cylinder.  Hence, in the numerical solution it is possible that the value of  β  tends to its correct value for  

0Re = if one considers the values of the parameter β  at which the quantities such as the force or moment 
on the cylinder are zero as the Reynolds number tends to zero.  To investigate this possibility values of 

LD CC Re,Re  and MCRe  are calculated as functions of the parameter β  for 5.30 ≤≤ β  at 3,5Re =  
and 1,  with 20/1,10/1=h  and 40/1 .  Figures 3, 4 and 5 show the variation of LD CC Re,Re  and 

MCRe  as a function of the strength of the rotlet, β , for 3,5Re =  and 1, respectively, when using mesh 
sizes 10x40, 20x80 and 40x160.  It is observed from figures 4 and 5 that the variations of the lift and 
moment as a function of β  are very non-linear.  However, the drag, (see Figure 3) appears to behave 
almost linearly with respect to β , over the range 5.30 ≤≤ β  for all the values of Re  considered. Further, 
all the three values of Re  considered in this paper, predict that zero drag appears to occur close to the value 
of ( ) 3== cβ . It has already been established by Dorrepaal et al. (1984) that only for a stream whose 
direction and magnitude are related to the position and strength of the rotlet will a solution be possible when 

0Re = .  In this case only for the specific value of c=β  will the problem be well-posed and produce a 
solution for which the values for the drag, lift and moment on the cylinder are all zero.  For all other values 
of β  one will have the force and moment dependent on this parameter.  It has been shown elsewhere that 
when 0Re =  non-zero values of the drag and lift implies that there are terms of the type ( ) ( )ϑsinln rr  and 

( ) ( )ϑcosln rr  at large values of r .  Similarly non-zero values of the moment require a ( )rln  term.  
Obviously the presence of such terms violates the boundary condition that 0=f  as ∞→r , and it is only 
when all such terms are absent that one will achieve the required solution.  This is why at 0Re = , it is 
impossible to obtain a convergent solution for an arbitrary value of β , unlike the situation at 0Re ≠ . 

However, for 0Re =  one can think of the lift, drag and moment as functions of β  and only when 
these quantities acquire the value of zero will β  achieve its required value, namely c=β , and produce a  



T. B. A. EL BASHIR 
 

 98

  

 
 

Figure 3. The variation of DCRe  as a function of β  ◊◊◊ .16040**,*8020,4010 ××+++×  

convergent solution.  Hence, the technique adopted was that of a Newton-Raphson  iteration on ( )βDC .  
The choice of DC  rather than LC  or MC  being based upon the linearity of this expression with respect to 
β , compared with the rather complex behavior of the latter two quantities when considered as functions of 
this variable.  The governing partial differential equation in finite-difference form, namely equation (3.1) 
and (3.2), were solved for two values of β , say 1β  and 2β , the number of iterations being terminated after 
a prescribed number since due to the boundary condition posed at 0=ξ  convergent solutions are 
impossible.  The resulting two non-zero values of ( )1βDC  and ( )2βDC are the values used in the Newton-
Raphson method in order to establish a new value of β , say 3β  at which ( )3βDC  is estimated to be zero.   
 



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 99

 

 
 

Figure 4.  The variation of LCRe  as a function of β  ◊◊◊ .16040**,*8020,4010 ××+++×  
( ) ( ) ( ) .1Re,3Re,5Re === cba  

 
Resolving for this new value of β , and using  ( )3βDC  together ( )1βDC  or ( )2βDC , whichever is closest 
to zero, to produce from the Newton-Raphson solution a new value of β . This process is continued until the 
value of ( )3βDC  reaches zero to within the required tolerance needed. The value of β  automatically 
results in the correct values of LC  and MC .  The whole iterative procedure could equally well have been 
applied to setting either the values of ( )βLC  or ( )βMC  to zero.  The resulting values of β  achieved by 
this process were 3.0004, 3.005 and 3.0010 for the mesh seizes 40x160, 20x80 and 10x40 respectively,  



T. B. A. EL BASHIR 
 

 100

 

 
 

 

Figure 5. The variation of MCRe  as a function of β  ◊◊◊ .16040**,*8020,4010 ××+++×  
( ) ( ) ( ) .1Re,3Re,5Re === cba  
 
compared with the exact analytical value of 3.  The streamline and vorticity patterns produced from the 
numerical computation and shown in Figures 6 (a) and 7 (a) are identical with those obtained from the 
analytical solution.  Figure 6 shows the streamline pattern obtained for 3=β  and 3,1,0Re =  and 5  for 
fine mesh size 40x160.  In these figures the effect of changing the Reynolds number appears to be slight, but 
if we compare this figure with Figure 10 the effect of  the Reynolds number becomes more clear.  

Figure 7 shows the vorticity pattern for 3=β  and 3,1,0Re =  and 5   for 40x160.  The effect of 
changing the Reynolds number is unclear, since the Reynolds number is small and the strength of the rotlet 
is strong. 

Figure  8 shows the streamline pattern for 5.1=β  and 3,1,0Re =  and 5 .   The region around the 
rotlet appears to be unaffected by the small change in the Reynolds number but compared with Figure 7, 
where the strength of the rotlet is double the present value, the decreasing strength of the rotlet causes a 
shrinkage of the region for the same streamline values. 

Figure 9 shows the vorticity pattern for 5.1=β  and 3,1,0Re =  and 5 .  Whilst the three figures are 
virtually indistinguishable, the decrease in the strength of the rotlet to half its earlier value, see Figure 7, has 
produced a more symmetric pattern of vorticity between the rotlet and the cylinder whilst at the same time   



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 101

                                                                  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 

Figure 6 (a,b,c, d).  The numerically obtained streamlines, Ψ , for 3=β  with a mesh size 16040 × . The 
streamlines labelled  1, 2 and 3 correspond to 7.00,1 and−=Ψ , respectively.  
( ) ( ) ( ) ( ) .5Re,3Re,1Re,0Re ==== dandcba  
the vortex lines of the same magnitude have been constrained in this region closer to the surface of the 
cylinder.  This is contrary to what appears to take place above the cylinder with the vortex lines of the 
same magnitude showing a departure from the surface. 

 
 
 

 
                                                            
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
Figure 7(a, b, c, d).  The numerically obtained vorticity, ω , for 3=β  with a mesh size 16040 × .  The  
streamlines labelled  1, 2 and 3 correspond to 4.00,4.0 and−=ω , respectively. 
( ) ( ) ( ) ( ) .5Re,3Re,1Re,0Re ==== dandcba  



T. B. A. EL BASHIR 
 

 102

 

 

 

Figure 8 (a, b).  The numerically obtained streamlines, Ψ , for 5.1=β  with a mesh size 16040 × .  The 
streamlines labelled  1, 2 and 3 correspond to 7.00,1 and−=Ψ , respectively.  

Figure 10 shows the effect of the strength of the rotlet and how the results tend to those 
obtained by Fornberg (1980). Figure 10 (f) is identical with the figure obtained by Fornberg (1980).  
Further, close to the body, the flow (when β  exceeds unity) appears to be dominated by the rotlet, 
even at a Reynolds number of 20.  The dominance of the rotlet is seen by the symmetrical nature of  

 
 

 

                
 
 

Figure 8 (c).  The numerically obtained streamlines, Ψ , for 5.1=β  with a mesh size 16040 × .  The 
streamlines labelled  1, 2 and 3 correspond to 7.00,1 and−=Ψ , respectively. 
( ) ( ) ( ) .5Re,3Re,1Re === cba  



FLOW PAST A ROTATING CIRCULAR CYLINDWER AND A ROTLET 
 

 103

                                                                          
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 

Figure 9 (a, b, c). The numerically obtained vorticity, ω , for 5.1=β  with a mesh size 16040 × .  The streamlines 
labelled  1, 2 and 3 correspond to 4.00,4.0 and−=ω , respectively.     

 

 
 
Figure 10 ( a, b, c, d, e, f).  The numerically obtained streamlines, Ψ , for 20Re =  with a mesh size 16040 × .  
The streamlines labelled  1, 2 and 3 correspond to 7.00,1 and−=Ψ , respectively. 
( ) ( ) ( ) ( ) ( ) ( ) .0,1.0,5.0,1,5.1,3 ====== ββββββ fandedcba  

 

 



T. B. A. EL BASHIR 
 

 104

the flow (when β  exceeds unity) appears to be dominated by the rotlet, even at a Reynolds number of 
20.  The dominance of the rotlet is seen by the symmetrical nature of the flow about the x-axis and y-
axis in the vicinity of the body.  As the strength of the rotlet decreases, the symmetrical flow pattern 
diminishes.  However, it is only when the value of β  falls below 0.5 that one sees the closed streamline 
behind the cylinder and even when the value of β  reaches 0.1 this closed region behind the cylinder is 
only present in the region 0<y  due to the rotation from the rotlet.  Obviously as β  tends to zero the 
closed streamline region will move to its symmetrical pattern about the x - axis. 

Conclusions 

Numerical solutions of the full, steady, two-dimensional Navier-stokes equations have been obtained 
for the fluid flow past a circular cylinder for Reynolds number 5,3,1Re =  and 20  and rotational 
parameter 20 ≤≤ α  when a rotlet of any given strength is placed at various positions outside the 
cylinder.  Although numerical solutions can be found for any value of Reynolds number greater than zero 
it is well known that when the Reynolds number is zero that a solution only exists provided there is a 
relationship between the strength of the rotlet and its location.   Further, when the Reynolds number is 
zero, the analytical solutions are such that the drag, lift and moment are all identically zero.  Therefore in 
this paper we have concentrated on obtaining numerical results with the drag, lift and moment zero and 
developed a technique for extrapolating the results obtained at small values of the Reynolds number to 
predict the solutions for zero Reynolds number.  Using this technique it is found that the numerical results 
are in excellent agreement with all the available theoretical results.  We conclude that the technique 
developed in this paper may be used with confidence to predict solutions for flows at zero Reynolds 
numbers for situations where there are no analytical solutions. 

References 

ALLEN, E., MARTIN, A. and PARIS, F., 1979. Boundary Elements in Potential and Elasticity Theory, Comp. And Stru., 10: 
351-362.   

DENNIS, S.C.R. 1960. Finite Differences Associated with Second-order Differential Equations, Quart. J. Mech Appl. Math., 13: 
487-496.  

DENNIS, S.C.R. and HUDSON J.D., 1978. A difference Method for solving the Navier-Stokes Equations, Proc, 1st Int. Conf. 
Numerical  Methods Laminar and Turbulent Flow, Pentech Press, London, 69-78. 

DORREPAAL, J.M., O’NEILL, M.E. and RANGER, K.B., 1984. Two-dimensional Stokes Flows with Cylinders and Line 
Singularities, Mathematika, 31: 65-75.           

FILON, L.N.G., 1926. The Forces on a cylinder in a Stream of Viscous Fluid, Proc. Roy. Soc., A113: 7-27. 
FORNBERG, B., 1980. A Numerical Study of Steady Viscous Flow Past as Circular Cylinder, J.Fluid Mech., 98: 819-855. 
FORNBERG, B., 1985. Steady Viscous Flow Past as Circular Cylinder, up to Reynolds number 600, J. Comput. Phys., 161: 297-

320. 
IMAI, I., 1951. On the Asymptotic Behaviour of Viscous Fluid Flow at a Great Distance from a Cylinder Body, Proc. Roy. Soc., 

A208: 487-516. 
PROUDMAN, I. And Pearson, J.R.A., 1957: Expansions at Small Reynolds Numbers for the Flow Past a Sphere and a Circular 

Cylinder, J. Fluid Mech., 2: 237-262. 
 
 
 
 
 

Received   1 November 1999  
Accepted   12 June 2000