Engineering Analysis with Boundary Elements 15 (1995) 235-247 Copyright © 1995 Elsevier Science Limited Printed in Great Britain. All fights reserved 0955-7997/95/$09.50
0955-7997(95)00027-5
ELSEVIER
Vortex lock-on for a rotationaHy oscillating circular cylinder a BEM numerical study M. M. EI-Refaee Kuwait University, College of Engineering, PO Box 5969, Safat 13060, Kuwait Vortex shedding resonance or lock-on is observed when a circular cylinder is excited into rotational oscillations. Boundary element numerical simulations of this flow at subcritical Reynolds numbers of 1500 and 3000 are computed. Two lock-on or resonance flow regimes are predicted. The primary lock-on regime occurred when the reduced forcing frequency was in the vicinity of one, while the tertiary lock-on was observed at a reduced forcing frequency around three. The frequency bandwidths of the two lock-on regimes are found to be correlated with the reduced rotational amplitude of the oscillating cylinder. Behavior of the aerodynamics forces during the lock-on regimes is investigated. It was found that it is the sharp variation of the phase shift between the phase of the velocity oscillation of the body and the phase of the initial vortex formation, that causes the remarkable increase in the aerodynamic forces within the lock-on regimes. The global subcritical lock-on characteristic diagram for a rotationally oscillating cylinder is obtained. The present BEM numerical model is verified by a series of tests. The tests include mesh-size refinement studies, comparisons with experiments and other solutions for a variety of benchmark problems and for a wide range of physical parameters.
Key words: Vortex shedding resonance, rotational lock-on, circular cylinder, vortex induced vibration, rotary oscillation, boundary element method.
NOTATION P a
//
£
o3
Oscillation amplitude Reference length, c = d Drag coefficient, Cd = D / ( p~oU~ / 2 )c Cd Lift coefficient, CL = l/(Poo U2/2)c CL Fluctuating lift coefficient eL Mean drag coefficient d d Diameter Forcing frequency f Frequency o f vortex shedding fv fv for stationary flow conditions fo F, Sv, SO Stronhal numbers, F = f d / U~, Sv = fv d/ Uoo,
F/So 1 Re
&/So U
O0 = lrF0
f~
Forcing oscillation angular amplitude Air density Kinematic viscosity Angular rate Angular oscillation velocity of the forced cylinder, f~ = 27rF
Subscripts s v O
b
Separated flow Vortex shedding from oscillating cylinder Vortex shedding from stationary cylinder Body Freestream conditions
So =fod/Vo
OG
Reduced forcing frequency Sectional lift Reynolds number based on d Reduced vortex shedding frequency Horizontal velocity
1 INTRODUCTION When a fluid flows past a bluff body, a boundary layer forms along the body surface for some distance, after which it separates from the surface due to the severe adverse pressure gradient. Because o f its inherent instability, the separated shear layer eventually rolls up into vortices that shed into the body's wake at a
Separation angle on circular cylinder Peripheral (rotational) speed Reduced amplitude o f peripheral speed o f oscillations 235
236
M. M. El-Refaee
particular frequency that depends generally on the body and flow parameters such as Reynolds number. The shedding vortices affect the pressure distribution on the body surface and therefore induce additional forces of a harmonic nature. Vortex shedding resonance or lock-on is observed when an oscillating excited cylinder drives the flow vortex shedding off its basic natural (stationary) frequency.l-5 The coupling between body oscillations and the vortex shedding for a circular cylinder has been extensively investigated both theoretically and experimentally6-15 in the last few decades. Transverse oscillations are almost invariably found in air flow, round chimneys, telegraph wires, pipelines and suspension bridges while in-line motions occur far less often in air. This is because, in air, the in-line component of the periodic force is usually one order of magnitude less than the cross flow component. In water flow, marine piles and risers, submarine periscope or braced members of offshore structures can be excited into motion in both the in-line and the transverse directions. In this case, the use of lightly damped bluff structures increases the likelihood of in-line oscillations because of the larger fluid density and because of the lower velocities at which the vibrations are initiated. Jones et al. 3'4 found in their experiments at supercritical Reynolds numbers of Re > 5'5 × 106, that the maximum transverse response occurred at fifo = 0-99 when the Strouhal number (S) is equal to 0.28. A question of great practical importance is the possibility of obtaining vortex lock-on at higher harmonics. Stansby l° showed that this can indeed occur for transverse oscillations at odd superharmonics o f f = f o, 3fo, etc. It should be noted that all of the experimental data show that lock-on occurs more readily when the oscillatory amplitude is large. Vandiver6 has indicated that the lock-on frequency bandwidth that permits the synchronization is mainly governed by the oscillation amplitude. The precise limits of this range is also somewhat influenced by the Reynolds number, turbulence, and the surface roughness. 7 Griffin et al.8 have shown experimentally as well as computationally that the lock-on frequency bandwidth gets narrower as the exciting oscillation amplitude decreases for subcritical Reynolds numbers flow regimes. Ericsson9 has reported that vortex resonance (lock-on) of longitudinal oscillations occur when the forcing frequency exceeds the second harmonics ( f ~ - 2 f o ) , while transverse oscillations experience lock-on near the odd numbers superharmonics of the resonant frequency. Erricson has used the phasing of moving wall effects to explain the occurrence of the lock-on at some sub- and superharmonics frequencies. Stansby, l° in his extensive work on lock-on of subcritical flows, classified three main lock-on regimes: primary, secondary and tertiary. They refer respectively
to the situation where Sv = F, F/2, F/3 and they occur respectively when F is close to So, 2So and 3So. Although subharmonics locking-on was not detected in Ref. 10, Ongoren and Rockwell 16 have observed subharmonic synchronization when F is close to So/2. Transverse primary lock-on frequency bandwidth is found 16 between 0"85
Vortex lock-on for a rotationally oscillating circular cylinder driving frequency is equal to or three times as large as the vortex shedding frequency for a stationary cylinder, the (maximum) moving wall effect will alternatively promote and oppose the vortex shedding event, resulting in a net effect that is small. In contrast, when the driving frequency is around twice (4, 6, etc., times) as large as the vortex shedding frequency, the moving wall effect will promote all the vortex shedding events, and the net effect is very powerful. The frequency bandwidth at which the first in-line oscillations lock-on occurs, depends mainly on the reduced amplitude (a/d), while it is slightly influenced by the value of the Reynolds number in the subcritical flow regimes. It is evident in the results of Hall and Griffin1 that large amplitude in-line oscillations have a much broader lock-on frequency range than that of the low amplitude oscillations. Further evidence of the insensitivity of the flow induced vibration to Reynolds number may be found in numerous reports of the large amplitude vibration of drilling risers and pilings. 8 However, the Reynolds number has some influence on the data at relatively high oscillations amplitudes. 21 Modes of vortex formation in in-line oscillations can be categorized into two groups. 22 There is a symmetrical mode for which a pair of vortices is shed in phase from both sides of the cylinder during one oscillation cycle. Antisymmetrical modes appear in five forms. In the first antisymmetric mode I, two vortices are shed from alternate sides of the cylinder and form a couple as they subsequently move downstream. One vortex of the couple may be stronger than the other as it develops a level of greater circulation. This regime is reported to appear for reduced frequency F/So more than 2. Mode II is similar to mode I but the period for the vortex pattern is double the first mode period. This mode is evidenced for 1.76 < F/So < 2.2, at least for experiments2° at Re = 190. Modes III and IV involve the formation of counter rotating pairs and occur only for pure in-line vibrations at the boundary between symmetric and antisymmetric modes. Mode III is found for {1"5So < F < 2 } . 23 The wake therefore becomes asymmetrical in mode III as a consequence of a transition between mode IV, which was found at lower values of F/So (for 0.5 < F/So < 1.5; F/So No. 1) and the symmetric mode S that was found f o r F/So = 2. A fifth in-line antisymmetric mode (V) has been found numerically by Lecointe et al. 24 It is characterized by the fact that one vortex is shed during one cycle of the motion, always on the same rear side of the cylinder. In the above literature review, the in-line and transverse features of vortex lock-on have been addressed. As described by many investigators, synchronized vortex formation (lock-on) can occur in different regimes (primary, secondary, tertiary, etc.). In the event that the cylinder is forced into rotary motion, the issue arises as to whether there are other locking-
237
on regimes and under what conditions they occur. Comparatively few researchers25-27 have shown, in their experiments, that rotational oscillations of circular cylinders cause vortex lock-on and result in marked changes in the flow characteristics. Their measurements were carried out at wide range of Reynolds number and for a limited values of the rotational oscillation amplitude. Okajima et a l Y showed experimentally that for two ranges of Reynolds numbers (40 < Re <_ 160 and 3050 _< Re <_ 6100) and at a forced oscillation frequency around the Karman (stationary) frequency (F/So ~ 1), the vortex shedding synchronizes with the cylinder motion. Wu et al., 28 have more recently done work similar to Okajima et al. 25 at Reynolds numbers of 300 and 500 and showed that at very high oscillation frequencies and magnitudes the stagnate-fluid region behind the cylinder, as well as the vortex shedding process, could be nearly eliminated. Tokumaru and Dimotakis27 have reported that rotational oscillation at very large magnitudes can produce a significant reduction in the drag on the cylinder. Filler et al. 26 investigated experimentally the frequency response of the vortex shedding from a circular cylinder subjected to small-amplitude rotational oscillations (0"005 < vo/Uoo <_ 0.03). They showed that these small oscillations can cause primary as well as higher harmonic vortex resonance (F/So ~ 1, ~ 3) at a moderate Reynolds number range (250 < Re <_ 1200). Vortex lock-on for rotational oscillations have numerous engineering applications in addition to their importance in fundamental physics. Applications abound in offshore exploration and drilling, naval and marine hydrodynamics as well as wind engineering and electric power transmission. In recent literature reviews, investigators25'27 stressed on the need to study the lock-on characteristics of vortex shedding of a circular cylinder subjected to rotational oscillations. Unfortunately, available data in the literature are far from being sufficient for understanding the relation between the rotational oscillation parameters and the lock-on phenomena. It is the objective of this paper to take a step in that direction. In the present work, the influence of the rotational oscillation parameters (forcing frequency and oscillation amplitude) on the vortex shedding resonance (lock-on) at moderate subsonic Reynolds numbers is investigated. The computational work is conducted by using the boundary element method (BEM) developed by Wu and his co-workers, 29-31 Skerget et al. 32 and El-Reface. 33 Although the boundary element method is considered as a newcomer in computational methods compared to the prevailing space-based discretization techniques such as the FDM (finite difference method) and FEM (finite element method), it has established itself as a mainstream tool in computational fluid dynamics. Different boundary element methods have been successfully applied to solve incompressible viscous fluid flow problems in the last decade. 29- 34 Mathematically, these
M. M. EI-Refaee
238
problems are expressed as Navier-Stokes equations, which is a system of nonlinear partial differential equations. Although the method loses some of its attraction when the governing equations are nonlinear, BEM still possesses many distinguishing features, especially for relatively high Reynolds number external flow problems. 29'3°'32 In BEM, the main iterations are carried out for the unknowns on the boundary; besides that, the dimensionality of the problem is greatly reduced by BEM which in term leads to a considerable economy in the numerical work. In contrast to the prevailing computational techniques, the BEM has the unique capability of confining the computations to the small vortical region. This feature drastically improves the solution speed. Many additional attractive features of the BEM are demonstrated in the literature. 29-3~ The remainder of the present paper is organized as follows. In section 2, the BEM mathematical formulation for both the kinetics and kinematics are presented. Numerical formulation as well as solution procedure are described in section 4. A number of benchmark test cases that demonstrate the ability of the method to accurately predict the flow details involved in the simulation of vortex shedding for circular cylinders are solved and compared with available solutions and data in section 5. Grid resolution refinement results for Re = 200 are also given in section 5. In the last part of the present work, the developed BEM algorithm is used to study the vortex lock-on of a rotationally oscillating circular cylinder.
kinematics part, described by eqns (2) and (3), or eqn (4) give the relationship between the vorticity distribution at any instant and the velocity vector at the same instant of time. 2.1 Boundary and initial conditions The associated velocity boundary conditions of twodimensional flow past a finite solid body immersed in an infinite medium are given as: 17 = 17~ at infinity ~ 17 = 17s on the body J
(5)
The above are Dirichlet types of boundary conditions corresponding to the kinematics part of the problem, where Vs is the surface velocity boundary condition. The kinetic equation (1) is parabolic in time and it is necessary to prescribe the vorticity field at t = 0 as the initial conditions. In the present study, the body is assumed to be impulsively started in an infinite fluid domain. The corresponding initial conditions are given by the non-circulatory potential flow in which the vorticity is zero everywhere except on the surface of the body. In addition to the initial conditions, it is necessary to specify the vorticity values at the inner and outer boundaries at all time levels. If the outer boundary is sufficiently far away from the body, the vorticity values may be assumed to vanish at this boundary. Wu 31 has shown that the vorticity values on the solid boundary at any time may be determined from kinematics considerations alone.
2 MATHEMATICAL FORMULATION 2.2 Mathematical formulation for kinematics The time dependent laminar flow of an incompressible viscous fluid is governed by the vorticity transport equation (kinetics) 29 as follows: 0to = - ¢. + (1) Ot where t is the time, 17 is the velocity vector, and w is the magnitude of the vorticity vector. The flow kinematics are described by the conservation of mass and the vorticity definition as follows: --
v. z = o
(2)
07 = V x I7
(3)
Equations (2) and (3) can be combined to give the Poisson's equation for the velocity vector: v 217 = - v
x ,2
(4)
Equations (1), (2) and (3) are to be solved to find the values of w and the two components of velocity vector 17, respectively. It is seen that the problem has been divided conveniently into two parts. The kinetic part describes the rate of change of the vorticity at any point due to convection and diffusion, eqn (1). The
Wu and Risk 29 have shown that eqns (2) and (4), together with the boundary conditions, eqns (5), can be recast in an explicit boundary domain integral relationship for the velocity vector in terms of the vorticity vector as follows:
¢(7, t) =
1 JfR~o(Fo, t)x(Fo-F)dR ° ( 7 - 70) 2
1 {B Vs(Fb't)'E°x(Fb-7) t7_7ol 2
rbdO+ (6)
where ? is the position vector, R is the region occupied by the fluid where the vorticity is non-zero, 17oois the free stream velocity, and the subscript 'o' designates that a variable or its integration is performed in the 70 space, ~7s is the velocity at the solid boundary and Eo is the unit normal vector at the surface B. The above boundary domain integral relation (BDIR) describes the kinematics aspect of the problem and is completely equivalent to the more familiar differential continuity and vorticity definition equations specialized for two dimensional flow with finite values of velocities at the surface.
Vortex lock-on for a rotationally oscillating circular cylinder 2.2.1 Law of conservation of the total vorticity Q Upon integration of the vorticity transport eqn (1) in the domain R, and by using the standard vector identities, one obtains the following form for the law of conservation of the total vorticity (Q):36 Q + 7rR2Q = 0
(7)
where f~ is the oscillation angular velocity of the forced cylinder, and the total vorticity is given as
(8)
Q = [ todR JA
at the body surface due to the outer vorticity field and the velocity at the surface respectively. Equation (1 1) alone is not sufficient to determine the surface vortex strength % however when eqn (11) is combined with the law of conservation of the total vorticity (7), vorticity values at the surface could then be easily computed. 2.3 Mathematical formulation for the kinetics
The vorticity transport eqn (1), could be rewritten as
- v (1 ( V . w V ) -
Ow
2.2.2 Surface vorticity determination The vorticity values away from the surface are determined using eqn (1); however, in order to solve the vorticity transport equation, it is necessary to prescribe the vorticity values on the solid surface at all time levels. A generally implicit scheme for determining the surface vorticity could be derived by the application of eqn (6) to the surface of the body where the tangential velocity is Vo(t) = ~ecos (2~rFt), where F represents the reduced forcing frequency of the angular oscillation. In order to do that, the viscous region is conveniently divided into a vortex sheet of strength 3' located on the surface and an outer vorticity field, where the vorticity to is assumed to be known. The tangential component of the velocity on the surface could be obtained from the scalar product of eqn (6) with the unit tangential vector/'as follows: 1 J~ "Yox_(Fb--Fo)dBo. ~. V . 7 = vo(t ) = ~ ~b :~o[ 1 IIR ta]°x(r'b- r'°)dRo.t" "q-2--~
1 ~b +~-~
gs(r'b't)'n°(r'b--r')'t" ir - Fo[2
239
0t
V2w }
(12)
=
The corresponding boundary-domain integral formulation ofeqn (12) can be obtained by applying the Green's theorem for scalar functions and using the parabohc time-dependent fundamental solution w* as follows:
(13) where R is the region occupied by the fluid where the vorticity is non-zero, B is the solid boundary, ~7is the unit normal vector on B directed away from R, the subscript '0' designates that a variable, its differentiation or its integration are performed in the 70 space and to* is the fundamental solution of the following differential equation
vEto* 4--- - -
rbdO+
Voo'r
0 V Ot and is given by Skerget et al. a2 as
(9) Here 7b is the position vector for points on the surface where the tangential component is to be calculated. It must be noted that the region R in the second integral does not include points on the solid surface. The first integral on the right hand side of eqn (7) can be written as: f
~°X(Fb- F°) dBo. t" .
.
=
.
.
V'o[2
JB 7
1 f
,,fo.~(Fb ~-~b--_~12
dBo.r
(10)
Substituting the above expression into eqn (9), one gets
"7 1 J? "~oX(Fb- FO) dBo .i. Vod + Vss + ~oo . r 2=2--~ #~'b [Fb-Fol 2 (ll) where voa and vss are the tangential velocity components
to* =f(F, t; 7o, to) = (4~rv(t - to) ) exp
4v(t - to)J (14)
Equation (13) describes the kinetic aspect of the problem and it is equivalent to the familiar differential vorticity transport equation. The fundamental solution to* gives the vorticity distribution in an infinite unlimited stationary fluid at the time level t as a result of diffusion of a concentrated vorticity of unit strength located at ro at a preceding time level to. Therefore, the second integral in eqn (13) represents the diffusion effect of an initial vorticity distribution in a stationary fluid. The first integral in eqn (13) describes the convective part of vorticity transport, while the third integral represents the effect of continual generation of vorticity on the solid boundary. The vorticity is generated at the surface in order to satisfy both the velocity boundary condition as well as the conservation of the total vorticity.29 The
M. M. El-Refaee
240
convective term (first integral) in eqn (13) can be simplified as
-I'oa oII ,, o J' II = -
dto
o
= -
Vo" (cO*WoVo)dRo
l O P _ _u(OW~ p Os \ OnJb
dto
o
+
The total drag for the present circular cylinder, with a prescribed distribution of surface velocity, is composed of two components: pressure drag, Cop and the friction drag CDc In general when the normal velocity vanishes, the tangential components of the momentum equation on the surface are written as
R
I'J
J' JI dto
COo(fro"VoW )dRo
(18)
Integrating eqn (18) along the body surface gives
(W*WoVo"no) dBo B
o
(15)
R
P - Pref = --pu
Ii ( ) aso
(19)
ref
where Pref is the reference pressure and is taken to be the pressure of the incoming stream. The pressure coefficient is defined by
Substituted by eqn (15), eqn (13) can be written as
w(r, t) =
3 DETERMINATION OF LOADS
I I n (W'W°)'°=0dR°
P
- Pref
Cp - l pV~ +
o
dto
;oJ ;J
+ u -
R
dto
dto
o
B
Wo(go
~
2
*
(w Vo w - aJoVoW ~ *)" fro dBo
4 (W*Wovo • fro) dB o
(16)
B
wb , no)u dBo dt o
C p -- ReD
J, (003) ref
O-ff
d3o
(20)
For a circular cylinder, eqn (18) becomes
The right hand side of the above boundary integral equation for the vorticity includes the unknown velocity and vorticity as well as their gradients in the solution domain and on the solid surface (Vw.ff)B. The last integral term in the above equation vanishes when the body is oscillating only in the angular direction. Since the vorticity boundary condition is required to solve the vorticity transport equation, and in order to avoid overspecifying the problem, eqn (16) is applied on the surface of the solid to obtain a relationship between the surface vorticity and its normal derivative. The resulting boundary integral equation (BIE) has the following form:
= u
When all quantities are nondimensionalized with respect to the free stream velocity and the radius of the cylinder, eqn (17), is reduced to
4 f2~ (0a3"~ Cp - ReD Jo \O~Jb dO
(21)
Here ReD is the Reynolds number based on the cylinder diameter. Similarly, the shear stress on the body surface is given by 7- = puw. In terms of the normalized surface vorticity a3b, the skin-friction coefficient is given as Cf
7- ! pV 2 -
2 Re D
2
&b
(22)
The drag coefficients Cop, C d f a r e calculated by horizontally integrating the pressure coefficient Cp, and the friction coefficient Cf along the cylinder surface as follows:
117
CDp = ~
CpcosOdO
(23)
CDf = ~
Cf sin 0 dO
(24)
COoVow b • fro dBo
+½W(Tb't)- JIR (WbW°)'°=°dR°
- j'odto J L Wo(o. Vowa)dRo
Finally, the total drag coefficient is expressed by the addition of the two drag components as
(17)
where (Owo/Ono)b is the vorticity gradient and the subscript b denotes points on the solid boundary at which the vorticity gradient is to be computed. As will be shown later, eqn (17) is used to determine the vorticity g r a d i e n t (OWo/Ono)b at the solid surface.
CD. r = CDp "~ CDf
4 NUMERICAL F O R M U L A T I O N In order to obtain numerical solutions of the field variables, velocity and vorticity, one has to discretize
Vortex lock-on for a rotationally oscillating circular cylinder the solution boundary into E boundary elements with Are boundary nodes, and the solution domain into C cells with Arc interior nodal points. The values of the variables w, (Ow/On), w~ are assumed to vary within each element or cell according to linear function in the two coordinates x and y.31
4.1 Boundary element approximation for the kinetics A more suitable version of eqn (16) can be obtained if the timewise integration is performed from t - At to t (At is the time increment) rather than from 0 to t. The resulting equation takes the form:
II
w(7, t) =
(W*Wo)t-atdR + v
R
I
t--At
dto
X
Is (W*VoWo "fro - WoVoW* "fro) dBo
+
I'
-
t--At --At
dto
IJ
dto
Wo(~7o• VOW*)dRo
R
(W*WoV-'o"fro) dBo
(25)
For rotational oscillations, the last integral vanishes since /70 "fro = 0 on the surface. Evaluating the area and boundary integrals in the above equation for each domain and boundary elements, the resulting discretized form of the BIE is given as 34
transportation. Matrices [A] and [B] are the convective coefficient matrices of orders (Nc x ~r), while matrices [C] and [D] represent the diffusion coefficient matrices (Arc x Arc). The matrices {w}, {uw}, {vw} and { f } have order of (Ate x 1), while {Wb} and (Ow/On)b are (Arc x 1) matrices. The matrix {~o} depends only on flow variables at the preceding time level t - At. All other column matrices are to be evaluated at the new time level t. The rapid decay of the fundamental solution w* with position results in sparsely populated coefficient matrices [A], [B], [C] and [/9]. As mentioned before, the vanishing tangential velocity boundary condition provides a mechanism for vorticity generation on the boundary B. This vorticity leaves the surface mainly by diffusion and is then distributed in the flow field by both convection and diffusion. The vorticity leaving the surface can only reach a certain distance during a certain time increment At. This region of non-negligible vorticity is known as the vortical region, the size of which depends on the Reynold's number and the time increment.
4.2 Discretized boundary-domain integral equation for the surface vorticity Applying the kinematics boundary integral relation for the surface vorticity [eqn (11)] at all the boundary nodes, results in the following matrix system: [S]{w~} = [M]{wl} + {T}
-'~
I"~k'pOJP p=l
+Z
+
Hk'q -~n q=l
Uk,p (uw)t~ + E
p=l
+ Z Wk'q°')q q=l
Vk,p(vw)~
(26)
p=l
where k refers to the nodal point at which the value of w is to be computed, p refers to a non-negligible vorticity point, Nc is the total number of such nonnegligible vorticity points, q refers to point on the solid body and Ne is the total number of points on the body surface. The geometric coefficients Lk, p, Hk, q, Wk, q, Uk,p and Vk,p are obtained through temporal and spatial integration of the fundamental solution w* or its derivative (Ow*/On). The method of evaluating these coefficients are similar to the methods that are usually used in a standard finite dement technique. 34 It can be shown that the above geometric coefficients depend merely on the relative positions of the nodes k and p and on the time increment At, and the time level t. For convenience, eqn (26) is written in a matrix form as: =
241
(28)
where {wi} is the interior vorticity (Arc × 1) vector, [M] is the interior vorticity (Atc x Arc) geometric coefficient matrix, [S] is the surface vorticity (Are x Are) matrix, {wb} is the surface vorticity (Are x 1) vector and {T} being a known (Are) vector that includes the contribution of the boundary conditions as well as the total domain vorticity; 19 {Wb} could then be determined by multiplying the right hand side of eqn (28) by the matrix inverse IS] -1.
4.3 Discretized boundary-domain integral equation for the interior vorticity In order to obtain the values of the vorticity normal derivatives at the boundary, eqn (27) is particularized on the boundary to yield 0w - [Cs]{wb} -- {qo}
(29)
+ IS]{,,,.} +
The vorticity gradient (Ow/On)b is obtained by multiplying the above equation by [Ds]-l. Then, the values of {(Ow/On)b} and {wb}, are substituted into eqn (26) to
Ow
give:
{wl} = [A'][uw]+ [B']{vw}+ {A} Equation (26) describes the mechanism of vorticity
(30)
where {wi} is the vorticity values in the interior domain.
M. M. EI-Refaee
242
The size of the geometric coefficient matrices [A'] and [B'] in the above equation is proportional to the number of vortical nodes in the solution domain which is relatively small at the early time levels and is increasing as the time advances. The vector {A} includes the effect of the boundary conditions and the potential flow velocity. A successive under-relaxation iterative scheme is then used to solve for the interior vorticity {wi} in eqn (29). The values of {w}, u and v on the right hand side are lagging by one iteration such that: {&i} k+] = [A'][uw] k + [B'][vw] k + {A}
(31)
The value of above {dJi}k+l is then relaxed after each iteration as follows {03I} TM : / ~ o ( ~ I } k+l -'1- (1 -- /~o){&l} k}
(32)
where k is the iteration counter, and /3o is the under relaxation parameter. The iteration process is continued until a prescribed limit of convergence is achieved.
(3) New values of the surface normal derivative (Ow/On)b are computed explicitly using eqn (29). (b) In the kinematic step, the velocity components are explicitly computed using eqns (33) and (34). At early time levels, the velocities are computed in a confined region of non-negligible vorticity. As time advances, the kinematics computations extend to a larger region until they cover all of the entire computational domain at later time levels. In each time level, once the vorticity gradient is computed from eqn (29) and the surface vorticity is obtained from eqn (28), the value of the pressure coefficient and the skin friction coefficient are evaluated from eqns (21) and (22). The load coefficients are then calculated from eqns (23) and (24).
4.4 Numerical treatment of the kinematics
5 VALIDATION TESTS
The velocity components are obtained from the explicit integral relationship between the velocity field and the vorticity field, namely, eqn (6). The area integral is evaluated by dividing the flow field into triangular elements where the vorticity is assumed to vary linearly. The total integral is replaced by the sum of integrals over elements associated with non-zero vorticity. For computational purposes, the velocity components u and v are given in matrix form as: 3°
In order to evaluate the accuracy of the present method in predicting the details of vortex lock-on, a number of benchmark test problems were solved and compared with the available data and solutions. The first test case is the natural or unforced vortex shedding from a circular cylinder at a Reynolds number of 200. For this condition the inflow boundary condition is a uniform flow, and the vortex shedding pattern is developed without forcing. The power spectral density corresponding to the x-component of the velocity Ux at a point (x, y) = (2, 2) is shown in Fig. 1, where the units are scaled by the cylinder diameter. The highest peak in the spectrum is located at the natural shedding frequency So of 0.188. Other peaks can be seen at higher harmonics of this shedding frequency. It is evident in the figure that the present spectrum distribution is qualitatively comparable with the numerical results of Hall and Griffin. 1 The sensitivity of the results to the grid resolution has been examined for the flow of Re = 200. Table 1
{u} = [U]{w} + {ub}
(33)
{v} = [V]{w} + {%}
(34)
where {w} refers to the non-zero vorticity vector, and [U] and [V] are the geometric coefficient matrices. Exact expressions for these functions can be found in Refs 30 and 34. The contributions of the far stream velocity as well as the velocities at the solid surface are included in the two vectors {Ub} and {%}. 4.5 Solution procedure
10 3
The overall numerical procedure to advance the solution by one time step At, consists of a kinetic step followed by a kinematic step: (a) In the kinetic step the following iteration scheme is used to solve for the vorticity at the new time level, n + 11: (1) Starting with an old set of known interior vorticity values {w~} at an earlier time level n, a new set of vorticity values at a new time level, n + 1, {w~'+l}, is iterated for convergence by using eqns (29), (30) and (31). (2) New values of the surface vorticity, {Wb}"+1, are obtained explicitly from eqn (28).
~
10z
--- Hall & Griffin [1] Present
10 I
s
10Q
I
~,~ 10-1 /
i
'I
I,
I I I!
[I ii iI i I
I
~
i
10.2 r~ 10-3 10--4
~'~ J
I
I
I
I
1 'w
I
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Frequency (S)
Fig. 1. Spectrum corresponding to the horizontal component Ux at a point in the wake (x, y) = 2, 2.
Vortex lock-on for a rotationally oscillating circular cylinder
243
T a b l e 1. ( R e = 2 0 0 )
Parameter Mean dra~ ceeff., Cd Fluctuating drag amp_, Cd Lift amp., CL Strouhal number, So Base pressure
Low resolution
Medium resolution
High resolution
grid (50 x 40) a
grid (60 X 60)
grid (80 X 80)
1.24
1.27
1.32
0-08 (So) 0"11 (So) 0-43 0.192
0-07 (So) 0"09 (2So) 0.51 0.189
0-04 (So) 0"06 (2So) 0-58 0.188
-0-93
-0.91
Experiment35
Other N u m . 24
-0"89
1.30
1.29
0.184
0.039 (So) 0"53 (2So) 0"569 0.195
-0"9
Cpb aGrid (50 × 40) means 50 equal spaced mesh in 0 direction, and 40 exponentially spaced mesh in the radial direction: 1 < r < 80. shows the computed values o f the force coefficients and Strouhal numbers for different resolutions. Other numerical solutions and data are also listed in the table. One can conclude from Table 1 that the Strouhal frequency decreases as the mesh resolution increases and the force coefficients are more influenced by the grid resolution than the shedding frequency. Moreover, as indicated in the table, the high resolution solution (80 × 80) is in good agreement with experiments. In this regard, the present author has showed earlier3° that the solution becomes weakly influenced by the grid resolution beyond the above high resolution grid (80 x 80). In the second validation test, a uniform flow around stationary circular cylinder at a higher Reynolds number o f Re = 1500 was solved and compared with experiments. Figure 2 displays the good agreement between the present predicted surface pressure distribution and the experimental data of Norberg. 35 Table 2 presents the values o f the mean drag coefficient Cd, the Strouhal frequency So and the mean base pressure coefficient C'pb for both the present solution and experiments. A close inspection o f Table 2 reveals a good comparison between the present solution and the experiments. 35 In order to demonstrate the ability of the present technique to capture the details of the lock-on regimes, a third numerical experiment was conducted to determine the lock-on characteristics o f an in-line oscillating circular cylinder at a Reynolds numbers of 4000. Figure 3
Table 2. (Re = 1500) Parameter
Present
Mean pressure drag coeff., (~dp Strouhal frequency, SO Mean base pressure
~
0.890
0.216
0.212
-0.765
O""T~
I
OI
I
I
I
I
I
l
, 0.00 0.05
0.10
0.15
0.20 0.25
0.30 0.35 0.40 0.4.5 0-50
Forcedf~quucy O~
Fig. 3. Lock-on for in-line oscillations at Re = 4000.
0.6 Re = 80-1OO(5.19]
--
Re = 2OO,*oIIiom(ll
~1 b -
0.S
-
-0.75
2.0!- _o
•
o,I- \
0.925
cocFF., Cpb
o Exp[35]
I.O~
Experiments35
A
0.4
O0
!1201
- - k .. 2OO.imwm,oimio, - - lie = 4000, Immem ,olmJom
+ Re = 4OO0[19]
Q
0.2
-I.$ 0.1
-2.0
o
I
2o
I
4o
6o
*o Angk5
too
12o
140
ioo
Fig. 2. Surface pressure distribution over an unforced circular cylinder at Re = 1500.
"~:Ig on
~s0
(e)
qon lock-on
0.0
I 1.2
1 1.4
I 1.6
I !.8
I 2.0
I 2.2
I 2.4
I 2.6
lock.on
I 2.8
3.0
F/So Fig. 4. Lock-on characteristic diagram for in-line oscillations
244
M. M. EI-Refaee
shows the variation of the frequency ratio (F/Sv), fluctuating lift coefficient CL and the mean drag coefficient Cd as a function of the in-line forcing frequency F. Comparison with experiments 19 shows a good agreement between the present results and experiments. A strong relation between the flow force coefficients and the lock-on frequency is also evident in Fig. 3. The lift and drag coefficients considerably increase within the lock-on regime. Further evidence of the suitability of the present scheme to accurately predict the lock-on characteristics is achieved by reproducing the primary lock-on frequency bandwidth of the in-line oscillating circular cylinder. 1 The present predictions are compared with the available numerical solution I and measurements 5A9'2° in Fig. 4, where the results are plotted as function of the oscillation amplitude (2a/d) and the reduced forcing frequency (F/So) for various subcritical Reynolds numbers (100 < Re <_4000). The amplitude (2a/d) represents the peak-to peak displacement in the longitudinal direction. The remarkable conclusions that one can draw from this figure are: (i) the present predictions (at Re = 200, 4000) are in good agreement with the experiments; 5'19 (ii) the present solution at Re = 200 is in excellent agreement with the numerical solution of Hall and Griffin;1 (iii) the lock-on frequency bandwidth is getting narrower as the oscillation amplitude decreases; (iv) at low and moderate oscillation amplitudes (2a/d) < 0.2, the lock-on frequency bandwidth is mainly governed by the oscillation amplitude and not by the value of the Reynolds number; (v) for higher oscillation amplitudes, the lock-on bandwidth is slightly increased with an increase in the Reynolds number.
6 ROTATIONAL LOCK-ON (RESULTS AND DISCUSSION)
It is evident in the above validation tests that the present BEM scheme is accurately capable of capturing the features of the in-line lock-on characteristics. It is, therefore, used to study the rotational lock-on characteristics of an oscillating circular cylinder, which is the main objective of the present work. In order to determine the rotational lock-on regimes, the frequency bandwidth of each r_egime, and the behaviour of the force coefficients (Cl, Cd) during the lockingon regimes, the present algorithm is used to obtain solutions for different values of the lock-on controlling parameters: the reduced forcing frequency, F/S o, and the reduced rotational oscillation amplitude, ~0Figure 5 gives the reduced vortex shedding (Sv/So), the phase shift (A
,,
/
"u14
",,
o%jC, I •
t
0.0
++I~+~I"~P ~ ........
1
1
I
I
"
"
Sv/So = I
....... 1.0
I
~'l~::::::l /"
l
.
0.5
~I -
l#
1.5
~ 0 5 I- ~ • F ~
q~
1.5
.... 2.0
I
I
I
S,JSo = 113
" Present n~dution - - Prinmry l~k-on re@tin 2.5
3.0
3.5
4.0
4.5
5.0
Reduced forcing frequency (F/So)
Fig, 5. Variation of the reduced shedding frequency, the average phase shift, and the force coefficientsas functions of the reduced forcing frequency. It is seen in the figure that there is a substantial departure of the reduced vortex shedding frequency (Sv/So) from unity, not only in the vicinity of F/So = 1.0, but also at higher superharmonic frequency, F/So = 3. In fact, the variation of (Sv/So) with (F/So) clearly exhibits a resonant-type behaviour in two frequency ranges: 0.6 < F/So <_ 1.3 and 2.4 _< F/So <_3.3. The first range is the primary lockon regime and the second one represents the tertiary lock-on regime. The reason that the primary lock-on as well as the tertiary lock-on occur in the vicinity of F/So = 1 and 3 is best explained by using the phase shift theory27 and the moving wall effect (ejection of circulation). 28 In Fig. 9, the phase of the near-wake vortex formation ~v relative to the phase of the forcing velocity oscillation
Vortex lock-on for a rotationally oscillating circular cylinder Phase o Phase
of velocity oscillation of vortex shedding
1 I I I (b) F/S0 -- 3.0
I
I
I
I
i
I
t
I I I I (a) F/S 0 = 1.0
I
I
t
I
I
t
I
(c) F/S0 = 2.0
.VYx/X/V2/
-Ir/2
~- lV/2 -Ir/2
~,w/2
0.0
I I I I I I I I I I I 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Reduced
time
(F~)
Fig. 6. P h a s e o f t h e f o r c i n g velocity oscillation ~b, relative to p h a s e o f v o r t e x s h e d d i n g ~v as a f u n c t i o n o f t h e r e d u c e d time at t h r e e values o f t h e r e d u c e d f o r c i n g f r e q u e n c y .
the case when the forcing frequency is twice the vortex shedding frequency (F/So = 2), the vortex shedding (Fig. 6(c)) is alternatively in-phase ( A ~ ) = 0), then out-of-phase (A~ =7 0 with the body oscillating velocity. The moving wall effect first promotes separation, and then opposes separation every other vortex shedding event, causing the net effect (ACL, ACd) to be small. When the reduced frequency (F/So) departs from 1
Time = 21.863
245
and 3, the phase shift A~, remarkably, retains the same steady behaviour as in case (a) (Fig. 6) but with a phase shift A ~ > 0. The magnitude A~, however, mainly depends on how far the reduced frequency is departed from 1 or 3. Except for (F/So) equal to 1 or 3, the vortex shedding experienced an almost (not fully) in-phase synchronization within the two lock-on regimes. The relations between the time-averaged phase shift , ~ and the force coefficients (CL, Cd) are also illustrated in Fig. 5. The rapid variation in the phase shift during the lock-on regimes results in an increase in the force coefficients.25 In this regard, it should be noted that the increase in the lift amplitude AGE during the rotational locking-on is somewhat smaller than the corresponding increase for the in-line or transverse locking-on. This, indeed, is attributed to the generation of a secondary forcing circulation that opposes the initial Karman circulation which in turn, damps the adverse effects of the synchronized vertical force. The increase in the mean drag Add, however, is more pronounced during the rotational locking-on as compared to that of the in-line or transverse lockingon. Generally, the separation point moves further upstream as the forcing peripheral velocity component increases. In the case of in-line and transverse oscillations the peripheral velocity component represents only a fracture of the oscillating velocity, whereas it is exactly equal to the total oscillating velocity for the rotational
Time = 21.863
Time = 21.863
Time = 23.863
Time = 23.863
Time = 23.863
Time = 31.863
Time = 31.863
Time = 31.863
(a) F/So = 0.0
(b) F/So = 1.0
(c)
F/So = 0.45
Fig. 7. T i m e e v o l u t i o n o f flow p a t t e r n s at t h r e e settings o f the r e d u c e d f o r c i n g f r e q u e n c y .
246
M . M . EI-Refaee
oscillation. It is therefore anticipated that separation will be further promoted when the cylinder is forced into rotational oscillation which in turn increases the resultant horizontal force (drag). Figure 7 compares the vortex structure in the near wake among three cases of reduced forcing frequency settings (F/So) and at three different time levels. The first case represents the no forcing (stationary) case (a), while the second case (b) gives the flow structures when the reduced frequency lies in the primary lockon regime (F/So ~ 1.0). The third case is for a value of the reduced frequency outside the lock-on regime (F/So = 0-45), case (c). As indicated in the figure, the vortex structures in the near wake for the primary lock-on case, are qualitatively similar to those of the no forcing case. However, the lock-on flow exhibits a considerable increase in the upper and lower vortex strengths as compared to the no forcing flow which indeed is a property of synchronization. On the other hand, it is clear in Fig. 7(c) that there will be no appreciable changes in the initial vortex structures as the reduced frequency takes values outside the lock-on regimes (F/So = 0-45), which is a sign of non-synchronization. So far, the rotational lock-on characteristics have been predicted at a specified value of the reduced oscillation amplitude, ~0 = 0.35. In order, however, to construct a global lock-on characteristic diagram for a rotationally oscillating circular cylinder at subcritical Reynolds numbers, the present scheme is further used to obtain predictions at different values of reduced amplitudes and Reynolds numbers. The lock-on characteristic diagram, Fig. 8, which collects all the predictions, marks the boundaries between the lock-on and non-lock on regions. A close look at Fig. 8 reveals the fact that the lock-on frequency bandwidth depends mainly on the value of the reduced oscillating velocity amplitude and weakly on the Reynolds number. The boundaries for the lock-on and non-lock on regions
o Present, b = 3000 m P r ~ t , 1~ = 1500
..= !.5
8
"6
l~rm~ lock-on regime
1.0
~
o Eat~rim~ts[261
Tertairy lock-on re~r~ o
~
)
o
0.5
-g ogJ "o
0.0
0.5
1.0
1,5
2.0
2.
3.0
3.J
4.0
4.5
Reduced forcing frequency (]P/So)
Fig. 8. Rotational lock-on characteristic diagram.
5.0
compares favourably with the only small amplitude data of Filler et al. 27
7 CONCLUSION Based on the results and discussion presented in the last section, the following conclusions can be drawn for the parameter ranges covered in the present investigation. (i) The efficacy and validity of the present BEM to handle massive separated flows are examined in a series of numerical experiments. The tests include a grid refinement test as well as a number of benchmark problems. (ii) It was concluded that it is the ejection of the circulation into the near wake, rather than the wall jet effect, that initiates the rotational lock-on. Essentially, rotational lock-on occurs when the circulation ejected by the oscillating rotational motion of the body is synchronized with the initial vortex shedding from the stationary cylinder. (iii) The near-wake vortex structure is phase-locked with the body oscillating velocity during two regimes of synchronization, the primary lock-on regime (F/So ~ 1-0), and the tertiary lock-on regime (F/So ~ 3.0). (iv) The frequency bandwidth for each regime depends mainly on the reduced rotational velocity amplitude ~0. It gets narrower as the oscillation amplitude decreases. (v) Moreover, it was found that the above frequency bandwidth is weakly related to the value of the Reynolds number. (vi) The average value of the phase shift /Xq, between the initial vortex shedding and the body oscillating velocity changes sharply during the two lock-on regimes. The sharp change in /X,ep is accompanied by a considerable increase in the force coefficients (CL, Cd). (vii) The increase in lift (A(~]) during the rotational lock-on regimes is smaller than that of the other corresponding forcing oscillating mechanisms (in-line and transverse motions). (viii) The increase-in drag (ACd), however, is much stronger than that of the corresponding in-line and transverse oscillations. (ix) The synchronization characteristics defined in the lock-on characteristic diagram (Fig. 8) hold only for the considered oscillation amplitudes (0-02 < ~0-< 1"6). At higher amplitudes, however, the bandwidth of synchronization is weakly correlated with the amplitude; at the other extreme, i.e. at a sufficiently small amplitude, the synchronization range may cease to exist altogether. (x) It is evident in Fig. 8 that the present predictions compare favourably with the only experimental small amplitudes data of Filler et al. 26 (xi) While the above results were obtained for moderate subcritical Reynolds numbers (Re = 1 500300), preliminary evidence (Fig. 8) suggests that the description of the flow phenomena presented here is qualitatively the same over a larger range of subcritical Reynolds numbers.
Vortex lock-on for a rotationally oscillating circular cylinder REFERENCES 1. Hall, M. S. & Griffin, O. M. Vortex shedding and lock-on in a perturbed flow. Journal of Fluids Engineering, 1993, 115, June, 283-90. 2. Parkinson, (3. V. & Ferguson, N. Amplitude and surface pressure measurements for a circular cylinder in vortexexcited oscillation at subcritical Reynolds numbers. Paper 18, Meeting on Ground Wind Load Problems in Relation to Launch Vehicles, Langley Research Center, 1966. 3. Cincotta, J. C., Jones, (3. W. & Walker, R. W. Experimental investigation of wind-induced oscillation effects on cylinders in two-dimensional flow at high Reynolds numbers. Paper 20, Meeting on Ground Wind Load Problems in Relation to Launch Vehicles, Langley Research Center, 1966. 4. Jones, (3. W. Jr, Cincotta, J. C. & Walker, R. W. Aerodynamics forces on a stationary and oscillating circular cylinder at high Reynolds numbers. NASA-TR, 1969. 5. Tatsuno, M. Vortex street behind circular cylinder oscillating in the direction of flow. Bulletin of Research Institute for Applied Mechanics of Kyosho University, 1972, 36, 25-37. 6. Vandiver, J. K. Dimensionless parameters important to the prediction of vortex-induced vibrations of long, flexible cylinders in ocean currents. Journal of Fluids and Structures, 1993, 7, 423-55. 7. Ericsson, L. E. Karman vortex shedding and the effect of body motion. AIAA Journal, 1980, 18 (8), 935-44. 8. Griffin, O. M. & Ramberg, S. E. Some recent studies of vortex shedding with applications to marine turbulars and
risers. ASME Journal of Energy Resources Technology, 1982, 104, 2-13. 9. Ericsson, L. E. Steady and unsteady vortex-induced asymmetric loads, review and further analysis. AIAA Paper 79-1531, Williamsburg, VA, 1979. 10. Stansby, P. K. The locking-on of vortex shedding due to cross-stream vibration of circular cylinders in uniform and shear flows. Journal of Fluid Mechanics, 1976, 74 (4), 641-65. 11. Nuzzi, F., Magness, C. & Rockwell, D. Three-dimensional vortex formation from an oscillating, non-uniform cylinder. Journal of Fluid Mechanics, 1992, 238, 31-54. 12. Armstrong, B. J., Barnes, F. H. & Grant, I. A comparison of the structure of the wake behind a circular cylinder in a steady flow with that in a perturbed flow. Physics of Fluids, 1987, 30, 19-26. 13. Belvins, R. D. The effect of sound on vortex shedding from cylinders. Journal of Fluid Mechanics, 1985, 161, 217-37. 14. Williamson, C. H. K. & Roshko, A. Vortex formation in the wake of an oscillating cylinder. Journal of Fluid Structures, 1988, 2, 355-81. 15. Rockwell, D. Active control of globally-unstable separated flows. ASME International Symposium of Nonsteady Fluid Dynamics (Proceedings), 1990, FED-Vol. 92, 379-94. 16. Ongoren, A. & Rockwell, D. Flow structure from an oscillating cylinder, Part I: Mechanisms of phase shift in near wake. Journal of Fluid Mechs., 1988, 191, 197-224. 17. Ericsson, L. E. & Reding, J. P. Dynamic stall analysis in light of recent numerical and experimental results. Journal of Aircraft, 1976, 13, April, 248-55.
247
18. Swanson, W. M. The Magnus effect: a summary of investigation to date. Journal of Basic Engineering, 1961, 83, Sept., 461-70. 19. Tanida, Y., Okajima, A. & Watanbe, Y. Stability of circular cylinder oscillating in a uniform flow or in a wake. Journal of Fluid Mechanics, 1973, 61 (4), 769-84. 20. Griffin, O. M. & Ramberg, S. (3. Vortex shedding from a cylinder vibrating in line with an incident uniform flow. Journal of Fluid Mechanics, 1976, 75, 257-71. 21. Barbi, C., Favier, D. P., Maresca, C. A. & Telionis, D. P. Vortex shedding and lock-on of a circular cylinder in oscillatory flow. Journal of Fluid Mechanics, 1986, 107, 527-44. 22. Ongoren, A. & Rockwell, D. Flow structure from an oscillating cylinder, Part II: Mode competition in the near wake. Journal of Fluid Mechanics, 1988, 197, 225-46. 23. Couder, Y. & Basdevant, C. Experimental and numerical study of vortex couples in two dimensional flows. Journal of Fluid Mechanics, 1986, 173, 225-52. 24. Lecointe, Y. & Piquet, J. Periodic and multiple periodic behavior of locked-in vortex shedding. Proceedings ofl6th ONR Symposium on Naval Hydrodynamics, Ed. National Academy of Sciences, Washington, 1986, pp. 470-89. 25. Okajima, A., Takata, H. & Asanuma, T. Viscous flow around a rotationally oscillating circular cylinder. Institute of Space & Aeronautical Science, University of Tokyo, Report 532, 1975. 26. Filler, J. R., Marston, P. L. & Mih, W. C. Response of shear layers separating from a circular cylinder to small amplitude rotational oscillations. Journal of Fluid Mechanics, 1991, 231, 461-89. 27. Tokumaru, P. T. & Dimotakis, P. E. Rotary oscillation control of cylinder wake. Journal of Fluid Mechanics, 1991, 224, 77-90. 28. Wu, J., Mo, J. & Vakili, A. On the wake of a cylinder with rotational oscillations. AIAA 2nd Shear Flow Conference, Tempe, Florida, 1989. 29. Wu, J. C. & Risk, Y. M. Integral - - representation approach for time dependent viscous flow. Lecture Notes in Physics. Springer-Verlag, New York, 1978, pp. 558-64. 30. E-Refaee, M. M. Boundary layer control of separated flow over circular over cylinders - - a BEM parametric study. International Journal of Boundary of Element Methods, (in press). 31. Wu, J. C. Velocity and extraneous boundary conditions of viscous flow problems. AIAA Paper No. 75-47, AIAA 13th Aerospace Science Meeting, 1975. 32. Skerget, P., Alujevic, I., Zagar, Z. & Rek, Z. Boundary element method for laminar and turbulent flow of incompressible viscous fluid. ZANN. Z. angew. Math. Mech. 1988, 68, T413-T416. 33. EI-Refaee, M. M. The total pressure boundary element method (TPBEM) for steady flow problems. Journal of Engineering Analysis with Boundary Elements, (in press). 34. Risk, Y. An integral - - representation approach for timedependent viscous flows. PhD thesis, Georgia Institute of Technology, Atlanta, 1980. 35. Norberg, C. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio. Journal of Fluid Mechanics, 1994, 258, 287-316. 36. Sankar, N. L. & Wu, J. C. Viscous flow around oscillating airfoils - - a numerical study. AIAA llth Fluid and Plasma Dynamics Conference, Seattle, 1978.
