- Email: [email protected]
A numerical analysis of two-dimensional flow past a rectangular prism by a discrete vortex model
Compeers and RuMs Vol. 10, No. 4, pp, 243-259, 1982
0045-T930~2/04~A3-17503.00/0
PrintedinGreatBritain.
© Ig~2PersamonPressLtd.
A NUMERICAL ANALYSIS OF TWO-DIMENSIONAL FLOW PAST A RECTANGULAR PRISM BY A DISCRETE VORTEX MODEL S. NAGANO College of General Education, University of Tokyo, Komaba, Tokyo, Japan and M. NAITOI" and H. TAKATA Department of Aeronautics, University of Tokyo, Hongo, Tokyo, Japan (Received 27 May 1981; in revisedform 18 March 1982) Abstr=ct--A discrete vortex model is developed to analyse the two-dimensionalfully separated unsteady flow past a rectangular prism. The effects of viscous diffusion of vortices and the loss in vorticity &fterthe shedding of vortices are incorporated into the model. Numerical computations are carried out for a stationary prism with a thickness ratio ranging from 0.5 to 2.0. The formation of Karman vortices, the vortex shedding frequency, and the fluid forces acting on the body are favourably compared with the experimental results by various research workers. The method of analysis is also shown to be applicable to the flow past a prism that is in forced vibration. 1. INTRODUCTION
There are two main procedures for the numerical analysis of unsteady separated flow past bluff bodies. One is to solve the Navier-Stokes equations using a finite difference or a finite element method. This method is limited to flow at low Reynolds number of less than about a few thousand at the present time. The other procedure that seems to be promising for high Reynolds number flow is the use of the so-called discrete vortex model. The discrete vortex model, however, cannot yet be said to have been fully established for various types of bluff bodies. Examples that have been analysed with some success by the method are limited to flows past a few simple bodies [1-8]. It may be worth noting that the point vortices, shed from the separation points in these flows are usually swept away from the body down to the rear of the body. The ability of the discrete vortices to represent continuous vortex sheets and vortex clouds in the wakes can be said to have been verified in these studies. In the present study, we try to analyse the two-dimensional flow past a rectangular prism. One of the difficult points in this flow arises from the fact that the vortices shed from the front edges have to go a long way near the surface of the body. They not only have to represent the vortex sheets and vortex clouds in the downstream wake, but also simulate the recircnlating flow near the surface of the body, where interactions between the vortices and the body will play an important role. In this sense, the present study may serve as a test case of the applicability of a discrete vortex model to flow past more complicated bluff bodies. 2. FORMULATION
OF THE
MODEL
2.1 Con/ormal mapping In a discretevortexmodel, the separatedflowpasta bluffbody isrepresentedby an inviscid potentialflow plus a large number of vortices.The boundary conditionimposed along the surfaceof the body is thatthe normal velocitycomponent should vanishthere.This condition can be satisfied by introducing images of the potential vortices. We, therefore, transform the physical plane to a plane where image vortices can be found easily. The transformation consists of two steps in the present study. In the first step the exterior of a rectangle in the physical plane (z-plane) is mapped to the interior of a rectangle in an intermediate plane (~-plane). If we denote the lengths of the lateral and longitudinal sides of the rectangle in the physical plane by tPresent address: CivilAviation Bureau, Ministry of Transportation,Chiyoda, Tokyo, Japan. 243
S. NAGAN0 et aL
244
2b and 2c, respectively, the mapping function is given by
z=-~A {zn~ + (E-k'2)~l + c
(I)
with -
ibk
where k is a constant which varies with the thickness ratio c
c/b of the
rectangle as (Kober[9])
'
E'-keK
-b = E - k ' - ~ '
K and E being the complete elliptic integrals of the first and second kind respectively. They are defined by
K = K 2' k
=
(1
-- k 2
sin 2 ~)-1/2 d ~
and
-
Jo
(I
ke sin e ~)i/2 d~.
The quantities with a prime are defined as k ' = ( l - k ~ ) I/2,
ql" ~ t K'=K(-~,k')andE'=E(-~,k).
zn in eqn (1) and sn and dn that will appear in the following are the Jacobian elliptic functions. In the second step the interior of the rectangle in the £-plane is mapped to the exterior of a circle of radius R in the transformed plane (A-plane). The mapping function is given by
= R(dn~,- iksn~).
(2)
Though the radius R of the circle can be chosen arbitrarily, the scale factor at the far field can be made unity if we choose R = b l 2 ( E - k'ZK). The outline of the mappings from a square to a circle is shown schematically in Fig. 1. It is seen that the coordinate meshes in the X-plane are quite similar to those in the z-plane at a small distance from the body.
2.2 Complex potential and induced velocities The complex potential function w in the circle plane may be written as
w= U(A + ~ ) - i(), - ~ ) 21rhf.sin(21rft) +~iFm{Iog(A-Am)
log(A
R---~2)I
(3)
where the first term on the r.h.s, describes a non-separated potential flow past a circle, the second term an oscillatory flow which is necessary when the circle is vibrating transversely (at a frequency f and an amplitude h) in a plane normal to the approaching flow, and the third term describes the flow due to M potential vortices and their images in the circle. In eqn (3) r,,
A numerical analysis of two-dimensionalflow past
't
245
, D z= ~,b)
z=(-c,b)
G'
z-Diane z = (c, - b)
z = ( - c , -b)
A
H
E
(2K')
H'
F
/ ~ ( K ' )
r-Diane G' (-K)
GO
c__..
(K)
u '
;'
G
~,-DIone
-
'
Fig. 1. Conformalmappingof a rectangleontoa circlein two steps.
denotes the circulation (positive when rotating clockwise) of a vortex at 7,m and an over bar denotes the complex conjugate. The origin of the z-plane is placed at the center of the rectangle which is in transverse vibration. The velocity at an arbitrary point in the physical plane is given, if the point does not coincide with the center of a vortex, by dw dw d?, dw ?2 u - iv = d--z"= d--~-'dz= d--~-'{?,"+ 2R2(2k 2- 1)7,2+ R4}'n'
(4)
where dk/dz is deduced from eqns (1) and (2). In the evaluation of the velocity at a vortex center in the physical plane, one has to subtract from the total potential the complex potential corresponding to the vortex whose velocity is to be determined. It may be noted that the modified potential is not the one in the ?,,plane but in
246
S. NAGANOel' al.
the z-plane• This procedure gives as the velocity of the vortex of strength F. at position zn R2 _ R2 u. - iv. = { U ( I - A,-~. ) - 2i~'hf ( l + -~, ) sin(2wft )
+~irm
~ irs(
1
1
(re#n) ~n 2
x {A4 + 2R~(2k2_ 1)A2+ R4}~/2 iF. An{R2(2k 2 - 1);t.2 + R 4} ~-~ {An4 + 2R2(2k ~- 1)Anz + R4} 3t~"
(5)
2.3 Forces acting on a body The fluid forces acting on a body can be calculated from the generalized Blasius theorem •
2
c~
(6)
where p denotes the density of the fluid, X and Y represent the component of force in the x and y directions, respectively, and the integrals are to be carried out along a closed path around the body• There should be no vortices in the domain between the surface of the body and the integration path. Substituting the complex potential into eqn (6) and carrying out the indicated integrals, we finally have M
X + i Y = p ~_, Fm(vm- ium) + 16ip~r3hf2RZ(k 2- 1) cos(2crft) M
+ ip X
Fm[R2{)tm4 +
2R2(2k2- 1)Am2+ R4}-v2 (/~m + iVm)
m=l
+ (L2{~2 + 2R:(Ek 2- 1)L 2 + Re}-!/2 - 1)(urn - iv=)] •
g
OF.,
R2
(7)
where u, and vm represent the x and y components of velocity at the center of the ruth vortex in the physical plane, and again the second term on the r.h.s, of eqn (7) is necessary only when the body is vibrating laterally. The last term in eqn (7) is retained for the reason stated later in Section 3.4, and z* in this last term denotes the point of the image of the ruth vortex inside the integration path. In the present study, X and Y correspond to the drag and lift acting on the rectangle, respectively, which may be expressed in terms of the drag and lift coefficients defined as CD=9--~U
and
CL=9--~U.
(8)
3. METHOD OF CALCULATION 3.1 Introduction of vortices One of the most important points in a discrete vortex model is how to determine the strengths of the nascent vortices and where to introduce them in the flow field. We have made a systematic study of various models proposed so far including those of Clements[2], Kuwahara[5] and Sarpkaya[7]. The final method which is believed to be physically the most reasonable results from the following considerations. The vortices are shed from the edges into the shear layers as a result of the separation of the boundary layers along the faces of the prism. Therefore, the rate of vorticity shedding should
A
numericalanalysisof two-dimensionalflowpast
247
be determined through the relationship dr
1. 2 = ~ u,,
(9)
where U, is the velocity at the outer edge of the boundary layer. It is true, as the flow is unsteady, that the outer edge of the boundary layers will also change with the time. Here, however, a suitable fixed point is chosen as a substitute for the time varying outer edge of the boundary layer. At the front edges, for instance, the vortices are shed as a result of the separation of the boundary layer along the face AB. Therefore, the point where the rate of vorticity shedding from edge A (Fig. 1) is calculated is fixed at z = (- c - 8, b); 6 is the distance from the face to the fixed point. Note that the vorticity shedding at the rear edges may result from the separation of either the boundary layer along the longitudinal faces or one along the rear face. The fixed point for edge C is, therefore, either at z = (c, - b - 8) or at z = (c + 8 - b) depending upon the situation of the flow near the rear edge. The points for edge B and D are chosen in a similar way. Though the distance 8 could be different at the front and rear edges, it is assumed to be the same in order to make the number of disposable parameters as small as possible. There is no rigorous criterion to determine the distance 3. But, since 8 is a substitute for the distance from the surface to the outer edge of the boundary layer, one can estimate the order of its magnitude resorting to a boundary layer theory. If one regards the flow along the front face as the approaching stagnation flow, the thickness of the boundary layer at the front edges may be calculated as 8 - 3.5b1~/~ (e.g. Schlichting[10]). Hence, if Re is in the range 104 ~ 10~, then 812b is in the range 0.02-0.005. Other information may be provided from the velocity distribution near the edges in the discrete vortex representation of the flow. According to our preliminary calculations past a square prism, the distance from edge A to the position where the velocity component parallel to the front face took a maximum value varied between 812b = 0.005 and 812b = 0.015. Taking these results into consideration, the distance 8 is fixed at 812b = 0.01 in all the calculations discussed in this paper. The positions, where the nascent vortices are introduced, are determined so as to satisfy the Kutta condition (dwldA = 0) at all edges. Since the Kutta condition alone is not sufficient to determine the positions of the nascent vortices, they are assumed to appear on hypothetical streamlines of the separated flow. For instance, the vortices shed from the front edges are assumed to appear on the extension of AB in view of the Kutta condition that the flow separates tangentially at the edges. Note that the distances from the four edges to the positions of the respective nascent vortices are not fixed but are determined at every time interval so as to satisfy the Kutta conditions for the four edges simultaneously. The positions of introduction of the nascent vortices are thus allowed to oscillate. It is again noted that, although positive vortices always appear at edge A and negative ones at edge B, the sense of rotation of the nascent vortices cannot be determined apriori at the rear edges. The vortices shed from edge C, for instance, should be allowed to appear either on the extension of BC (when the boundary layer along the lower face separates and the vortices are negative) or on the extension of DC (when the boundary layer along the rear face separates and the vortices are positive). 3.2 Convection of vortices The velocities induced at the center of the vortices in the physical plane can be calculated by eqn (5). Therefore, the movement of the nth vortex may be calculated by
d,. (Ow). which is approximated by the following first order difference scheme
(dw),
z , ( t + A t ) = z , ( t ) + 3-i
.At.
(10)
248
S. NAOANOet al.
We must note, however, that the path of a vortex near the separation point is usually highly curved. Therefore, if we calculate the movement of the general vortices at a time step At, then the nascent vortices should be calculated at a much smaller time interval. Actually, the movement of the nascent vortices is pursued at At~5 (Sarpkaya[7]). The selection of the time step At is very important in a practical calculation. The accuracy of the calculation becomes higher, while the computer time required becomes greater, when a smaller time step is specified. The same argument also applies to the time interval (AT) between the introduction of the nascent vortices. After a large number of preliminary calculations the following time step and interval are employed throughout the calculations discussed in this paper.
At=O.OS(2b/U) and AT=2At. 3.3 DiMusion of vorticity As has been pointed out at the beginning of this paper, the vortices shed from the front edges of the rectangle have to go a long way along the surface of the body. Therefore, the interactions between the vortices and the body will play an important role in the present flow. Special attention must be paid to the treatment of vortices when they come close to the body. The usual assumption of removing vortices that have approached closer than a pre-fixed distance to the surface of the body does not seem to give a reasonable model in the present flow. In this connection, it is noted that a potential vortex induces the velocity whose magnitude is in inverse proportion to the distance from its center. While, the velocity q induced by a vortex in a viscous fluid at rest is given by
where v is the kinematic viscosity of the fluid, t the time, and r the distance from the center of the vortex (Lamb [l l]). q increases with r from zero to a maximum (at r = r,) and then falls asymptotically to zero. The distance r, is a function of time and viscosity, and ff the vortex is concentrated in a point at t = 0, the viscous core radius r, is given by r, = 2.242X/-~. The difference in the velocities induced by a potential vortex and a viscous vortex is slight when r is larger than r,, but it becomes serious as r gets smaller than r,. In other words the potential vortices tend to induce unrealistic velocities when they come too close to the surface of the body or too close to each other. It seems reasonable, therefore, to let each point vortex bear a characteristic radius ro which varies with the kinematic viscosity of the fluid and the time that has elapsed after the shedding of the vortex (i.e. ro =/3r,,/3 being a constant), and to remove from the calculation the vortices that have approached closer than the radius r0 to the surface of the body. Physically the removal may be regarded to correspond to the destruction of vorticity in the shear layer by interaction with the boundary layers along the body. Also, a second characteristic radius rfi (r~ =/3'r.,/3' being another constant) is introduced, and any two vortices are coalesced to an equivalent single vortex when their separation is smaller than the sum of the second characteristic radii in order to avoid the excessive induced velocities to each other, as well as to reduce the total number of vortices in a smooth way. Thus, the distances for the destruction and the coalescence of vortices are assumed to be proportional to the viscous core radius of each vortex, /3 and /3' being the constants of proportionality. It may be worth noting that the diffusion of vorticity in the separated flow is due not only to the kinematic viscosity of the fluid but also to the apparent viscosity due to turbulence. However, the correction to the viscosity may be absorbed in the selection of the values of/3 and/3'. Unfortunately, there is no rigorous way of determining the proper values of /3 and/3' on a theoretical basis. Hence, after carrying out a number of preliminary calculations on flow past a square prism, the constants/3 and ~' in the present study are assumed to be unity as the first approximation, that is,/3 =/3' = 1 and hence ro = rfi = r,. In order to see the potential ability of the present diffusion model of vorticity, separate numerical calculations have been carded out on flow near the leading edge of a semi-infinite flat
A numerical analysis of two-dimensional flow past
249
Re-lO
Re-103]
Re-t0
S"
°" f: oo
°0%°°0
~ % ° ooo- o~
o
0
0
0
0
0
""
o
%
° 0~
Fig. 2. Leading edge separation of a semi-infinite thick flat plate. All the results are those at T = 14. The recirculating flow region is still in the course of development at this time at Re = 106.
plate. Some of the results are shown in Fig. 2, where the Reynolds number is based on the thickness of the plate. Figure 2 is not quoted for quantitative discussions, but is shown just as an illustrative example. For instance, at low Re, the diffusion model as it stands will only poorly represent the actual diffuse vorticity field that would exist in the real fluid flow. Nevertheless, qualitative characteristics of the flow may be seen: When Re is small (Re ~ 102), the flow does not separate. At a little higher Re (Re ~ 103) only a short separation bubble is formed. With further increase in Re the vortices begin to make a recirculating motion in the separated flow region and this becomes greater as Re is increased. It may, therefore, be said that, qualitatively, at least, the diffusion model has some ability to simulate the interactions of the shear layer with the surface of the body. 3.4 Loss of vorticity It is widely known for various bluff bodies that the amount of vorticity found in the fully formed street vortices is some 60% or less of the amount generated in each shear layer before separation (e.g. Fage and Johansen[12]). There seem to be three main causes of the loss in vorticity. They are (i) the destruction of vorticity near the wall due to the interaction with the boundary layers along the body, (ii) the cancellation of opposite-signed vortices due to mixing and (iii) the dissipation of vorticity due to three-dimensional turbulence. In the discrete vortex representations of the flows tried by various research workers, some mechanisms for the loss in vorticity due to (i) and (ii) have usually been included, but the loss in vorticity in these studies has been much less than that observed in the real fluid. Our preliminary calculations past a square prism, where the vorticity loss due to (i) and (ii) were included in the manner stated in the previous section, revealed that the vorticity clusters formed just behind the body retained so strong a vorticity as to let many of the upstream elementary vortices hit the longitudinal faces of the prism, thus preventing the periodic shedding of vortices. It seems necessary, therefore, to include some extra mechanism for the loss in vorticity in our discrete vortex representation.t
tThe outline of this paper was presented at the Euromech Colloquium held at Imperial College in July 1979 (Bearman and Graham[18]). After the Colloquium, Prof. Sarpkaya kindly sent us copies of his two papers. We found in one of his papers that he and his colleague adopted a vorticity loss model which would work in nearly the same manner as ours (Sarpkaya and Shoaff[8]). Also we found in the other paper by him that, parallel to our application of the discrete vortex model to an oscillating rectangular cylinder (see Section 4), he had been applying his model to an oscillating circular ¢yfinder (Sarpkaya[19]).
250
S. NAOANOet al.
Presumably, the extra loss in vorticity in the real fluid is due mainly to (ih'), whose precise mechanism in the separated flow, however, is far from being known as yet. What we can do at the moment is to incorporate into the model an artificial loss of vorticity in a phenomenologlcal way. We, therefore, carried out an experiment in a low turbulence wind tunnel and measurements were made of the change in the strength of vortices shed from a square prism. The vorticity flux was evaluated from the time-mean velocity measurement at various stations from the front to the rear edges of the square, while the strength of downstream vortices was estimated in a manner analogous to Fage and Johansen[12]. Figure 3 shows the estimated variation in the strength of vortices with the streamwise distance. The strength of vortices decays only minimally just after the front edge separation (-0.5 _-
u
~'0.5 x
I
5,0
J
IO,O
I
15,0
x/2b
Fig. 3. Space-wise loss in vorticity obtained in the experiment for a square prism at Re = 2.5 x 104.
1.0
8
~-o.s
T)
T2
T
Fig. 4. TL,ne-wise loss in vorticity incorporated into the numerical model. F(TO/F(T =0)=0.99 and F(T2)/F(T = 0) = 0.61 where T1 = 2 d U and T2 = T1 + 8blU.
A numerical analysis of two-dimensional flow past
251
T ffi4 (half a period), all quantities are nondimensionalized hereafter by the far upstream uniform velocity U and the length of the lateral side of the rectangle 2b. Unless otherwise stated, the Reynolds number (Re ffi 2bU/u) of the flow is set at 2 x 104 in the following calculations. It must be noted that the Reynolds number of the flow appears in the model only through the diffusion process of the vortices. Although the Reynolds number, through which some qualitative effects of viscosity are taken into consideration, will not necessarily correspond rigorously to that of the real flow, the effects of changing the Reynolds number of the flow in the model will be examined later. 4. RESULTS AND DISCUSSIONS
The numerical calculations are carried out for four types of stationary rectangular prisms with a c/b ratio ranging from 0.5 to 2.0, and for a square prism that is in forced vibration at various frequencies. The development of the wake behind a stationary square prism after the impulsive start from rest is shown in Fig. 5; The small circles denote the positive vortices shed from the front upper edge, the x's the negative ones from the front lower, the triangles the positive and the crosses the negative ones from the rear upper, the y's the negative and the diamonds the positive ones from the rear lower edges, respectively. These marks do not represent the strength of the elementary vortices but just show their positions. The small amplitude vibration of the prism for half a period at the early stage (T < 4.0) is seen to be large enough to amplify asymmetry in the flow. One of the interesting points is that some of the vortices shed from the front edges get close to the longitudinal faces of the prism and then go upstream back into the dead flow region. Recirculating flow regions are formed between the shear layers and the longitudinal faces. Many of the vortices shed from the rear edges are also swept back into the recirculating flow regions.
T- 1 ~
I
I::?::: :~:;:,,
~
~`:
"
T" 2 ~'"~:
T" 3 ,~,,,¢J~"
T" 9 ~.~7;, ~
C.",', ...~,=
i.
T=11
l"l .7.,,~ ~
b
T-6
=x
Fig. 5. Growth of wakes behind a square prism afterthe impulsive start from rest. Numbers in the figure denote the non-dimensional time after the start.
252
S. NAGANOet al.
Figure 6 shows the fully developed quasi-steady periodic states of the evolution of the wakes during nearly a full cycle of computation for the four types of stationary prism. The double circles and the thick x's denote the centers of positive and negative vortex clusters that are replaced to equivalent vortices. The periodic shedding of vortices and the formation of the vortex clusters are quite clear when c/b is smaller than or equal to unity• The periodicity is a little obscured at c/b = 2.0. This point will be discussed again later. Figure 7 compares the enlarged picture of the typical vortex arrangement near the longitudinal face of the square prism with the separation streamline estimated from the mean velocity measurement in our experiment. The curve (A) in Fig. 7 shows the outer edge of the shear layer that was obtained in the experiment. The curve connects the positions at which the maximum time-mean velocity was detected in the traverse of a hot-wire at various stations from the front to the rear edges of a square prism. Curve (B) connects the positions where the time variation of the instantaneous velocity was the largest. Therefore, curve (B) may be regarded as to represent the time-mean position of the shear layer in the experiment. The vortex sheets formed by the elementary vortices in the numerical calculation are seen to well simulate the separation streamline in the real flow. The rate of shedding of vorticity from the four edges is shown in Fig. 8 as a function of time for the four stationary rectangular prisms. After the initial period of flow est~iblishment, the rate (a)
(b) ~,,..++
T-32
~+
~.'~.,
~ .;~
%'.+',~+
";t'.
.i'
,+
XM~ ~;.,: ...... ~P.a ~, .y,r.."..,
+,: ..... :..:t
;,,
x
~ _ + p ~ . . , ?;;
{4"~ •
-,,
~" + -
%
,v.+
Fig. 6. Vortexarrangementsovernearlya-fullcycleof steadyperiodicflow.(a), c/b = 0.5; Co),clb = 0.62.
A numericalanalysisof two-dimensionalflowpast • (C)
253
(d)
• • ~,
• ~,.,;,~
% ~
":~"
x
~-32
....~, ":.
:
.~,~,:~
.'~,' . ,
."
•
T-36 ~ ' k ~ , ~ . :
~ %":~o.
.~.. .%
.::'~:,
o
,
X
2"
."."~"';~i~,~.~...... T-38
"'~d~'"
",,.
f
Fig.6(Contd).Vortexarrangementsovernearlya fullcycleofsteadyperiodicflow.(c),clb = 1.0;(d),clb = 2.0.
is seen to vary periodically except for the case with cJb = 2.0. It is also quite clear that the rate of shedding of vorticity from the rear edges is considerably smaller than that from the front edges. This is not surprising since the rear edges are in the dead flow region and the flow velocity near them is usually small. It may also be noted that the periodic shedding is established in the shortest time when db is equal to 0.62. In Fig. 9 the time variation of the fluctuating lift coefficient is plotted. Here again the periodicity is quite clear (except for the case with cJb = 2.0) and these plots are used to determine the Strouhal number of the vortex shedding in each case. Since CL is plotted against the nondimensional time, the Strouhal numbers for (a)-(c) can be calculated as an inverse of the average period of the last few full cycles of the fluctuating lift. They are tabulated in Table 1. In accordance with the fact that the periodic shedding of vortices and the formation of vortex clusters are somewhat obscured in Fig. 6(d), where clb -- 2.0, the time variation of CL in Fig. 9(d) cannot be said to be fully periodic either. Therefore, a Fourier component analysis has been made of the unsteady lift to find the Strouhal number for this particular case. And it has been found that there are two dominant Strouhal numbers as shown in Table 1. In this connection it is noted that Otsuki et al. [13] also found in their experimental study two dominant CAF
10:4
- C
254
S.
T-26
NAGANO
a
(A) ,~ ,, x
~
,
~ j ~~ ~ ( B )
~'-°
,
T 28
(A ~o(Ao o
T-30
~.~----¢--- (B) ~o~-~ ~ o OoO ° ,
aL
°~
:
~ "
(A) o
(A)
o o o ~o~
o,
'
(B)
~
~
(
B
)
-'C;-.'~"~o- #
o~ ~
-
Fig. 7. Enlarged picture of vortex arrangements near the longitudinal faces of a square prism at various instants of time. ©, positive vortex near the upper face; x, negative vortices near the lower face. The positions of the negative vortices are shown with their y coordinatenegated.
X
2.0
00
o xx
(a)
O0 ~'0000 0 Xx
"o
"o 1 , 0
+
oo o ~ O"'vvo^
Or~"-~xx x X X x
xxXXx
0
xxx 0
xo
x ~^
o x0 ,-,x
xxxln,l~Sr.,n x,XOx 'x'0 OXx x v -~ -x ~00000,., x Xxxx x 0o00 xx x 0oO x x x
A
10
20 x x o x ~
(b) ~:'
+-'
x
~
0 x x x 00o00 O
"~o ] . , 0 + ,
~O Oo Vx
xXx O x xO
XxxxX
j~÷+_
0O°o
xy O 0x
0O0o0
.
.--
,
xx O
^o
--
o
x~ O
O^
x
~" "~ 1 . 0 --
o6~^
, ~,,_,
- ~ --OUO O
¥Ooo~oOXxx XxxX A
& + A& & + _ +~+ 4- ~ A~A÷ ~&~A.~;
0
~';+ ~"
+A~
-+
,
o,.,x
X
X XxX
^
0
O
°o°
,-,
xx ~)
x 0
XxxxX
o
x
~
~
~,A +÷ +4.
A 4-
x
°° o
x^O
Ox
Q^ x ~,.,,
o~X
O
x
xxX
vo00
"
,t
+++
T X
_
0
0
~o
O"
& 4.+ 4. A A 4. ~ A & ++& + -A A
0°o°
O
oOOo
O~ X __O UO ~ "
~
x
~J xx x O
30 x~ X
,
x
+ ++ " + - +
x
~÷+ez~.9.+$+,_a.+.+"^,,a,~
20
x x~'
0 x
x~O^O~XX x O O x x O ,.OOoO + Xxx x Oooo xxx uO
.&A . ~÷÷-14¢ A &
x_*~,~
oo
XX
x
T
~x
xO0
+÷++++~+
,~ X ~
x,.,
A
i0 2,0
xO x
20 ooo
x x,.,oU
++ ~A
0
30
Ox
Xxxx x "0Oo0 XXx
(C)
X ~ X ! WoX..xxXXXX
"~ "-.
0 0oO
^O°oo Xo v Ox
4.
o x
x ^~
10 2.0
xx
"
2.0 --
0 xo
A &
"0
Xx
o x
X
x
OXX O
Ooo+
OX
x O xxxX
O
,,
+ + 4-
A A .A&
*
*++++~,&~M.a.q..e+*a;÷+, T
30
(~)
0o 0
4.0
XXxOO x
x~ 1,0
o
X X
÷+÷ A A A + A A A4. ÷ ÷ & + +A
"
+ +
A
t+-+~=-~-~++M~-~+++~,'N~***~"~+*~* ~--~~*~-~t*~;+'~;* * 10 20 30 T
Fig. 8. Rate of shedding of vorticityfrom the four edges. (0 front upper, x front lower, A rear upper, and + rear lower edges)
Strouhal numbers for a stationary prismatic bar with rectangular cross section when the ratio c/b was in the range from 2.0 to 2.8. One was on the branch that was connected to the lower clb and the other on the branch connected to the higher clb. Another point which should be noted in Fig. 9 is the amplitude of the unsteady lift. It increases gradually as clb is increased from 0.5 to 1.0, will" e it tends to decrease when c/b is increased to 2.0. The root mean square amplitudes of the unsteady lift coefficient are tabulated in Table 1.
255
A numericalanalysisof two-dimensionalflowpast
2.0
(o)
0
o% _ o °~°°°°'~-
i0
20 ,.,°°°°o
_ ~°°°o°°°
,-°°°%., o-.°° Oo o 3() Ooooooo T
%ooo°
oo
-2,0 2.0
(b) oo
J 0
~0o 0000
n
°°°%;
°ooooO~
oo° 0
=
i0 o
0000
0
°° '3
o
o
0
20,o O00
0o
on
oOn -o
oI n 30000°0
°OoooO
o
o°°%
o T
-2.0
.01
(c)
ooo
o
oO°°O
i°0io
o
o°O o
o
°°oOoooO
300
0000 o
oo 0
-2.0
0
o
o
o
o
o° T
oo o
o
o
0000
(d)
0
o
o
oooo OooOo
2.0
o
o
10,9 v o o O0 0 Ooo
O0o
000 o, o 20 oo
o
o 0
0 O0
o
0%
0
0
oi
oo
o
0030
Oo ° ° ° o o
TOo
Oo
0
o° ° o
Fig. 9. Time variation of fluctuatinglift coefficient.(a), c/b = 0.5; (b), clb = 0.62; (c), c/b = 1.0; (d), clb = 2.0. Table 1. Resultsof the computationsfor the four types of stationaryrectangularprism. St denotesthe Stroulminumberand an overbardenotesthe timeaveragedvalue. C(x,out0t ion
c/b
St
~/(CL)
Co
(0)
0.5
0.161
0.584
2.2q7
(b)
0.62
0.153
0.797
2.253
(c)
1.0
0.142
1.173
2.059
(d)
2.0
0.102 0.20tl
0.916
1.181
The time variation of the drag coefficient is shown in Fig. 10. The mean level of CD is seen to be the highest when clb is around 0.5-0.62. It decreases gradually as clb is increased further. The time mean drag coemcient is also tabulated in Table 1. The period of the cyclic time variation of Co is clearly halved compared to that of CL when clb is equal to 0.5 and 0.62, but the behavior is not obvious when c/b is equal to or greater than unity. If the Fourier component analysis were made of the unsteady drag, the half period variation in CD would be detected for the latter cases too. Figure 11 compares the results obtained in the present numerical analysis with those of the experiments by various research workers. The top figure compares the Strouhal number. Although the Strouhal numbers obtained in the numerical calculations are generally about lO~ higher than the experimental ones, the general trend of the experimental results is well predicted. In the middle figure the root mean square values of the fluctuating lift coefficient are shown. Though not many experimental results are available, the agreement of the analytical
256
S. NAGANOet al.
3.0 (O)
2.0
o
o
0
0
o
o oo oOOO o°° o o°°Ooo° °oo °o o 000OoOooooooOooo%00oO0ooo o OoO oo
1.0
o
o° °o
°°°oO
Oo oO o
0
ooOOC~ooOOooOoo° 0
o
oOo
OoO
o
oOoooo °
Oo°°° °Oo
o
oooOOO o
o
oo
°o
o
ooo
o °°
oo° °o°
(b)
oo°°Oo
0
00 0
0~0 0
o=oO 0 0°°°000000°0°°=°°°°%°0 °O°°o O 0 0o~
O0
o
°o°o0 0 0
(C)
0
0 O
0
0%0oo0o°00o° _ ..,.,..oo _~ ^o°° 0 0 0 vv ~r~ 0 ou Oo ^o O 0 o o ~ o ~ °°oOOO °°o o °°o°°°°°° (cl) o" o oo -',.,o Oo i
1
i
i0
20
30
Fig. 10. Timevariationof drag coefficient.(a) c/b = 0.5; (b), db= 0.62; (c), db = 1.0;(d), c/b = 2.0.
results with the experiment is quite good. It is interesting to note that the root mean square value of the fluctuating lift coefficient seems to take a maximum at c/b between 1.0 and 2.0. As is seen
in the bottom figure, the mean drag coefficient obtained in the present study predicts the general trend of the experimental result quite well; although, the sharp peak in CD near clb = 0.6 has not been detected in the limited number of numerical calculations. Therefore, it may be said that the present model as a whole, is able to represent the separated flow past a rectangular prism. A further supplementary study has been made to see the effects of changing the Reynolds number of the flow in the present model. In the computation for the stationary square prism at Re = 2 x 104, the Reynolds number of the flow was suddenly changed to a different value at T = 35, the computation was continued further, and the results were compared. The result at Re = 2 x 105 was almost completely coincident with that at Re = 2 x 10~. When Re was reduced to 104, a slight difference was seen in the time variation of CL but the periodic pattern as a whole was similar to that at Re = 2 x 104. No obvious difference was found in the downstream vortex arrangements either. However, the result at Re = 2 x 103 was quite different from the other three cases. The Karman vortices were not formed and the periodic vortex shedding was not sustained. The radius of an elementary vortex grew fast at this Reynolds number and too frequent combination and destruction of vortices took place just after leaving the separation points. The consequence was the failure of the model to represent the continuous shear layers by a sufficiently large number of elementary vortices. The situation might have been improved by employing a larger time interval AT between the introduction of the nascent vortices to avoid the frequent combination of nearby vortices. We did not, however, try a detailed examination of this point. It may be worth noting that the time step At has to be made smaller if Re is greater than, say, 106. The reason is that, when Re is higher, the core radius of a vortex remains smaller. The vortices can then approach closer to the body or to each other. Therefore, the
A numerical analysis of two-dimensional flow past
257
0.2
0,1
I
,
I
1,0
c/b
2.0
c/b
2,0
1.0
,
0
!
,
I
1,0
0
3.o I
.e/J~."
\.
2.0
1,0
T o
. i
,
1,0
!
c/b 2.0
Fig. 11. Comparison of the present results (shown by O) with experiments by Bearman a al.[14] (A, Re = 2 - 7 x 104), Nakaguchi et al. [15] ( x, Re = 2 - 6 x 104), Nakamura et al. [16] (V, Re = 3.3 - 13 x 104), and Okajima et all17] (12],Re = 6 - 1 2 x 104).
time step At should be made smaller, as At ~ l/~/R-e to be rigorous, if the accuracy of the calculation is to be kept constant. To sum up the results of the supplementary study, the present model with the parameters given in the paper is suited for flow at Re in the order of 104- 10s. The order of Re can be raised easily by adopting smaller At. The present model is also applicable to the analysis of flow past a vibrating rectangular prism.'r To show this by examples, the numerical calculations are carried out for the square prism that is subjected to a sinusoidal transverse vibration in a plane normal to the uniform flow. The amplitude of the vibration h is fixed at 0.1, while the nondimensional frequency of the vibration So (= 2bf/U) is varied in the range from 0.125 to 0.155. An example of the vortex pattern formed behind a vibrating square prism is shown in Fig. 12. The lateral distance between the positive and negative vortex streets is seen to be much smaller when compared to the vortex patterns behind a stationary prism. Figure 13 shows the time variation of the fluctuating lift obtained for the square prism in forced vibration at various frequencies. The displacement of the prism Y is also shown by the thin broken line. It is known that, when a bluff body in a stream is forced to vibrate, the vortex shedding frequency is locked to the frequency of the forced vibration if So is close to the natural vortex shedding Strouhal number St of the non-vibrating body. According to the experimental work by Otsuki et a1.[13], the locking range for a square prism in forced tSee footnote on p. 249.
258
S. NAOAN0 et al.
0
x
Fig. 12. Typical vortex pattern
2.0
behind a
vibrating square prism,
(e)
d
o
°o
O0
o
o°°
o'---.. . . . - "
o
000
0
0
2.0 I
(f)
-o O0o
o°°o
o
o
o
01=OOoO;,oo,0o:,o
_.y
~-- "-''-... /-" "-'"",..o
"" ' L
|
000
.... "
0
o
2~,,.~ "-. o
0
oo
0000
,oo
:o ..j "'--
o .... "
30
0000
0
°°°
-2, O I" 2,0
(g)
J
Oo
o
"P°o o ° / : ~ ' % o 10 . / ; , ~ ' , . - ~ --o.~"" ~ "-~ / °o° "oo
0 o
°°°°
---',
0
Y
...-."d°o
(h)
2,0 o
_
_..,,,_,,,,~
..o
0
~'--JY'o
o
O0
-, 30
._.
./
o
"'R
o°"J'-'l;
T" O0000 °
0O
o
000
00
oOoO o o
o
o
o
00
0 O0
-2.0
0
0 00
o / ....
""--% ° ' ' / -ZU " o.,
0
.... " T
°°°
00 0
T
"----
"
0
00000
000 /
°°o°
30,,........~
"'-- . . . . .
00
00
o
°°
0
-2.0
(J
S~ = 0.135.
o°°° °°°°° °°°° Y o °° ,-~-.....0"-..,i 0 o ....o..... .. 2,0 ,,,.I.......~o
o "-..o
0
0
o Y
o o
.
"~-~"
°ooO
°°-"'~'
20
°°°o
o o
o
" .... "
o
o o Ooo
o'°-J30
o
",z_ oo
oo
-2.0 Fig. 13. Time variation of fluctuating lift on a vibrating square prism. (e), S~ =0.125; (f), S~ =0.135; (g), So = 0.147; (h), S, = 0.155; ---, displacement of the square prism.
transverse vibration with a nondimensional amplitude 0.1 was 0.90_< SJ St -<1.08. Figure 13 shows that the phase difference between C,. and Y varies with the time at (e) SJSt = 0.88 and (h) SJSt = 1.09, i.e. the locking-in does not take place in these cases. While at (f) SJSt = 0.95 and (g) SiS, = 1.04, the phase differences approach asymptotically to constant values, and the vortex shedding is locked to the v~ration of the prism. F~lhermore, it is clear that the fluctuating lift lags behind the displacement of the prism when S~ is larger than St as in (g), while CL is ahead when $, is smaller than St (f). These phase differences between CL and Y in the locking-in regime are in good qualitative agreement with the experimental observations. It will be noted that the vorticity loss model shown in Fig. 4, which was assumed based upon the measurement of flow around a stationary square prism, has been conventionally approximated in the calculations of the flow past the vibrating square prism. There is no obvious reason, however, why the vorticity loss of vortices shed from a vibrating body may be assumed to be the same as that for a stationary one. As the vorticity loss has large influence on the magnitude of fluid forces acting on a body, a further study on the mechanism of loss in vorticity behind a vibrating body seems to be necessary in order to make more detailed quantitative discussions on the fluid forces acting on a vibrating body. This point is left for the future study.
A numerical analysis of two-dimensional flow past
259
5. CONCLUDING REMARKS
A discrete vortex model has been established to analyse the two-dimensional separated flow past a rectangular prism. The flow field around a rectangle has been transformed to the exterior of a circle via elliptic functions. The effects of the viscous diffusion of vortices and the loss in vorficity have been incorporated into the model. The diffusion model seems to have some ability for simulating interactions of the shear layers with the body. It is also effective in enabling the formation of the recirculating flow region along the surface of the body as well as in avoiding excessive induced velocities. Perhaps the most critical part of the model lies in the incorporation of the artificial loss of vorticity. It would be much better if the loss of vorticity could be incorporated into the model in a form of a general turbulent dissipation concept. However, with a careful choice of the parameters stated in the paper, the present model as a whole has been shown to have a si~ificant ability for representing separated flow past bluff bodies even when the interactions between the shear layers and the bodies play an important role in determining the characteristics of the flow. The numerical computations have been carried out at the Computer Center at the University of Tokyo. The authors would like to express their gratitude to Mr. Y. Machida who carried out the experimental work.
REFERENCES I. A. J. Chorin, Numerical study of slightly viscous flow. J. Fluid Mech. $7, 785 (1973). 2. R. R. Ciements, An inviscid model of two-dimensional vortex shedding. J. Fluid Mech. $7, 321 (1973). 3. R. R. Clements and D. J. Maull, The representation of sheets of vorticity by discrete vortices. Prog. Aero. ScL 16, 129 (1975). 4. M. Kiya and M. Arie, A contribution to an inviscid vortex-shedding model for an inclined fiat plate in uniform flow../. Fluid Mech. 82, 223 (1977). 5. K. Kuwahara, Numerical study of flow past an inclined fiat plate by an inviscid model. J. Phys. SOc. Japan 35, 1545 (1973). 6. K. Kuwahara, Study of flow past a circular cylinder by an inviscid model. J. Phys. Soc. Japan 45, 292 (1978). 7. T. Sarpkaya, An inviscid model of two-dimensional vortex shedding for transient and asymptotically steady separated flow over an inclined plate. J. Fluid Mech. 68, 109 (1975). 8. T. Sarpkaya and R. L. Shoaff, An inviscid model on two-dimensional vortex shedding for transient and asymptotically steady separated flow over a cylinder. A/AA Y. 17.11, 1193 (1979), 9. H. Kober, Dictionary of Conformal Representations. Dover, New York (1957). 10. H. Schlichtlng, Boundary Layer Theory. McGraw-Hill, New York (1960). 11. H. Lamb, Hydrodynamics. Dover, New York (1945). 12. A. Fage and F. C. Johansen, On the flow of air behind an inclined plate of infinite span. Proc. Roy. Soc. A, 116,170 (1927). 13. Y. Otsuki, K. Washizu, H. Tomizawa and A. Ohya, A note on the aeroelastic instability of a prismatic bar with square section. J. Sound Vib. M, 233 (1974). 14. P. W. Bearman and D. M. Trueman, An investigation of the flow around rectan._gulat cylinders. Aero. Quart. 23, 229 0972). 15. H. Nakaguchi, K. Hashimoto and S. Muto, An experimental study on aerodynamic drag of rectangular cylinders. J. Japan Soc. Aero. Space SOl. 16, 1 (1968). 16. Y. Nakamura and T. Mizota, Unsteady lifts and wakes of oscillating rectangular prisms..A. SOc. Civil Engng 151, EM6, 855 (1975). 17. A. Okajima, T. Mizota and H. Takata, An experimental study on the pressure distn~oution along the surface of transversely oscillating rectangular cylinders in a uniform flow. Preprint lor JSME Meeting No. 774-13 (1977). 18. P. W. Bearman and J. M. R. Graham, Vortex shedding from bluff bodies in oscillating flow: A report on Euromech 119. J. Fluid Mech. 99, 225 (1980). 19. T. Sarpkaya, Vortex induced oscillatlons--a selective review. J. Appl. Mech. 48-2, 241 (1979).
0045-T930~2/04~A3-17503.00/0
PrintedinGreatBritain.
© Ig~2PersamonPressLtd.
A NUMERICAL ANALYSIS OF TWO-DIMENSIONAL FLOW PAST A RECTANGULAR PRISM BY A DISCRETE VORTEX MODEL S. NAGANO College of General Education, University of Tokyo, Komaba, Tokyo, Japan and M. NAITOI" and H. TAKATA Department of Aeronautics, University of Tokyo, Hongo, Tokyo, Japan (Received 27 May 1981; in revisedform 18 March 1982) Abstr=ct--A discrete vortex model is developed to analyse the two-dimensionalfully separated unsteady flow past a rectangular prism. The effects of viscous diffusion of vortices and the loss in vorticity &fterthe shedding of vortices are incorporated into the model. Numerical computations are carried out for a stationary prism with a thickness ratio ranging from 0.5 to 2.0. The formation of Karman vortices, the vortex shedding frequency, and the fluid forces acting on the body are favourably compared with the experimental results by various research workers. The method of analysis is also shown to be applicable to the flow past a prism that is in forced vibration. 1. INTRODUCTION
There are two main procedures for the numerical analysis of unsteady separated flow past bluff bodies. One is to solve the Navier-Stokes equations using a finite difference or a finite element method. This method is limited to flow at low Reynolds number of less than about a few thousand at the present time. The other procedure that seems to be promising for high Reynolds number flow is the use of the so-called discrete vortex model. The discrete vortex model, however, cannot yet be said to have been fully established for various types of bluff bodies. Examples that have been analysed with some success by the method are limited to flows past a few simple bodies [1-8]. It may be worth noting that the point vortices, shed from the separation points in these flows are usually swept away from the body down to the rear of the body. The ability of the discrete vortices to represent continuous vortex sheets and vortex clouds in the wakes can be said to have been verified in these studies. In the present study, we try to analyse the two-dimensional flow past a rectangular prism. One of the difficult points in this flow arises from the fact that the vortices shed from the front edges have to go a long way near the surface of the body. They not only have to represent the vortex sheets and vortex clouds in the downstream wake, but also simulate the recircnlating flow near the surface of the body, where interactions between the vortices and the body will play an important role. In this sense, the present study may serve as a test case of the applicability of a discrete vortex model to flow past more complicated bluff bodies. 2. FORMULATION
OF THE
MODEL
2.1 Con/ormal mapping In a discretevortexmodel, the separatedflowpasta bluffbody isrepresentedby an inviscid potentialflow plus a large number of vortices.The boundary conditionimposed along the surfaceof the body is thatthe normal velocitycomponent should vanishthere.This condition can be satisfied by introducing images of the potential vortices. We, therefore, transform the physical plane to a plane where image vortices can be found easily. The transformation consists of two steps in the present study. In the first step the exterior of a rectangle in the physical plane (z-plane) is mapped to the interior of a rectangle in an intermediate plane (~-plane). If we denote the lengths of the lateral and longitudinal sides of the rectangle in the physical plane by tPresent address: CivilAviation Bureau, Ministry of Transportation,Chiyoda, Tokyo, Japan. 243
S. NAGAN0 et aL
244
2b and 2c, respectively, the mapping function is given by
z=-~A {zn~ + (E-k'2)~l + c
(I)
with -
ibk
where k is a constant which varies with the thickness ratio c
c/b of the
rectangle as (Kober[9])
'
E'-keK
-b = E - k ' - ~ '
K and E being the complete elliptic integrals of the first and second kind respectively. They are defined by
K = K 2' k
=
(1
-- k 2
sin 2 ~)-1/2 d ~
and
-
Jo
(I
ke sin e ~)i/2 d~.
The quantities with a prime are defined as k ' = ( l - k ~ ) I/2,
ql" ~ t K'=K(-~,k')andE'=E(-~,k).
zn in eqn (1) and sn and dn that will appear in the following are the Jacobian elliptic functions. In the second step the interior of the rectangle in the £-plane is mapped to the exterior of a circle of radius R in the transformed plane (A-plane). The mapping function is given by
= R(dn~,- iksn~).
(2)
Though the radius R of the circle can be chosen arbitrarily, the scale factor at the far field can be made unity if we choose R = b l 2 ( E - k'ZK). The outline of the mappings from a square to a circle is shown schematically in Fig. 1. It is seen that the coordinate meshes in the X-plane are quite similar to those in the z-plane at a small distance from the body.
2.2 Complex potential and induced velocities The complex potential function w in the circle plane may be written as
w= U(A + ~ ) - i(), - ~ ) 21rhf.sin(21rft) +~iFm{Iog(A-Am)
log(A
R---~2)I
(3)
where the first term on the r.h.s, describes a non-separated potential flow past a circle, the second term an oscillatory flow which is necessary when the circle is vibrating transversely (at a frequency f and an amplitude h) in a plane normal to the approaching flow, and the third term describes the flow due to M potential vortices and their images in the circle. In eqn (3) r,,
A numerical analysis of two-dimensionalflow past
't
245
, D z= ~,b)
z=(-c,b)
G'
z-Diane z = (c, - b)
z = ( - c , -b)
A
H
E
(2K')
H'
F
/ ~ ( K ' )
r-Diane G' (-K)
GO
c__..
(K)
u '
;'
G
~,-DIone
-
'
Fig. 1. Conformalmappingof a rectangleontoa circlein two steps.
denotes the circulation (positive when rotating clockwise) of a vortex at 7,m and an over bar denotes the complex conjugate. The origin of the z-plane is placed at the center of the rectangle which is in transverse vibration. The velocity at an arbitrary point in the physical plane is given, if the point does not coincide with the center of a vortex, by dw dw d?, dw ?2 u - iv = d--z"= d--~-'dz= d--~-'{?,"+ 2R2(2k 2- 1)7,2+ R4}'n'
(4)
where dk/dz is deduced from eqns (1) and (2). In the evaluation of the velocity at a vortex center in the physical plane, one has to subtract from the total potential the complex potential corresponding to the vortex whose velocity is to be determined. It may be noted that the modified potential is not the one in the ?,,plane but in
246
S. NAGANOel' al.
the z-plane• This procedure gives as the velocity of the vortex of strength F. at position zn R2 _ R2 u. - iv. = { U ( I - A,-~. ) - 2i~'hf ( l + -~, ) sin(2wft )
+~irm
~ irs(
1
1
(re#n) ~n 2
x {A4 + 2R~(2k2_ 1)A2+ R4}~/2 iF. An{R2(2k 2 - 1);t.2 + R 4} ~-~ {An4 + 2R2(2k ~- 1)Anz + R4} 3t~"
(5)
2.3 Forces acting on a body The fluid forces acting on a body can be calculated from the generalized Blasius theorem •
2
c~
(6)
where p denotes the density of the fluid, X and Y represent the component of force in the x and y directions, respectively, and the integrals are to be carried out along a closed path around the body• There should be no vortices in the domain between the surface of the body and the integration path. Substituting the complex potential into eqn (6) and carrying out the indicated integrals, we finally have M
X + i Y = p ~_, Fm(vm- ium) + 16ip~r3hf2RZ(k 2- 1) cos(2crft) M
+ ip X
Fm[R2{)tm4 +
2R2(2k2- 1)Am2+ R4}-v2 (/~m + iVm)
m=l
+ (L2{~2 + 2R:(Ek 2- 1)L 2 + Re}-!/2 - 1)(urn - iv=)] •
g
OF.,
R2
(7)
where u, and vm represent the x and y components of velocity at the center of the ruth vortex in the physical plane, and again the second term on the r.h.s, of eqn (7) is necessary only when the body is vibrating laterally. The last term in eqn (7) is retained for the reason stated later in Section 3.4, and z* in this last term denotes the point of the image of the ruth vortex inside the integration path. In the present study, X and Y correspond to the drag and lift acting on the rectangle, respectively, which may be expressed in terms of the drag and lift coefficients defined as CD=9--~U
and
CL=9--~U.
(8)
3. METHOD OF CALCULATION 3.1 Introduction of vortices One of the most important points in a discrete vortex model is how to determine the strengths of the nascent vortices and where to introduce them in the flow field. We have made a systematic study of various models proposed so far including those of Clements[2], Kuwahara[5] and Sarpkaya[7]. The final method which is believed to be physically the most reasonable results from the following considerations. The vortices are shed from the edges into the shear layers as a result of the separation of the boundary layers along the faces of the prism. Therefore, the rate of vorticity shedding should
A
numericalanalysisof two-dimensionalflowpast
247
be determined through the relationship dr
1. 2 = ~ u,,
(9)
where U, is the velocity at the outer edge of the boundary layer. It is true, as the flow is unsteady, that the outer edge of the boundary layers will also change with the time. Here, however, a suitable fixed point is chosen as a substitute for the time varying outer edge of the boundary layer. At the front edges, for instance, the vortices are shed as a result of the separation of the boundary layer along the face AB. Therefore, the point where the rate of vorticity shedding from edge A (Fig. 1) is calculated is fixed at z = (- c - 8, b); 6 is the distance from the face to the fixed point. Note that the vorticity shedding at the rear edges may result from the separation of either the boundary layer along the longitudinal faces or one along the rear face. The fixed point for edge C is, therefore, either at z = (c, - b - 8) or at z = (c + 8 - b) depending upon the situation of the flow near the rear edge. The points for edge B and D are chosen in a similar way. Though the distance 8 could be different at the front and rear edges, it is assumed to be the same in order to make the number of disposable parameters as small as possible. There is no rigorous criterion to determine the distance 3. But, since 8 is a substitute for the distance from the surface to the outer edge of the boundary layer, one can estimate the order of its magnitude resorting to a boundary layer theory. If one regards the flow along the front face as the approaching stagnation flow, the thickness of the boundary layer at the front edges may be calculated as 8 - 3.5b1~/~ (e.g. Schlichting[10]). Hence, if Re is in the range 104 ~ 10~, then 812b is in the range 0.02-0.005. Other information may be provided from the velocity distribution near the edges in the discrete vortex representation of the flow. According to our preliminary calculations past a square prism, the distance from edge A to the position where the velocity component parallel to the front face took a maximum value varied between 812b = 0.005 and 812b = 0.015. Taking these results into consideration, the distance 8 is fixed at 812b = 0.01 in all the calculations discussed in this paper. The positions, where the nascent vortices are introduced, are determined so as to satisfy the Kutta condition (dwldA = 0) at all edges. Since the Kutta condition alone is not sufficient to determine the positions of the nascent vortices, they are assumed to appear on hypothetical streamlines of the separated flow. For instance, the vortices shed from the front edges are assumed to appear on the extension of AB in view of the Kutta condition that the flow separates tangentially at the edges. Note that the distances from the four edges to the positions of the respective nascent vortices are not fixed but are determined at every time interval so as to satisfy the Kutta conditions for the four edges simultaneously. The positions of introduction of the nascent vortices are thus allowed to oscillate. It is again noted that, although positive vortices always appear at edge A and negative ones at edge B, the sense of rotation of the nascent vortices cannot be determined apriori at the rear edges. The vortices shed from edge C, for instance, should be allowed to appear either on the extension of BC (when the boundary layer along the lower face separates and the vortices are negative) or on the extension of DC (when the boundary layer along the rear face separates and the vortices are positive). 3.2 Convection of vortices The velocities induced at the center of the vortices in the physical plane can be calculated by eqn (5). Therefore, the movement of the nth vortex may be calculated by
d,. (Ow). which is approximated by the following first order difference scheme
(dw),
z , ( t + A t ) = z , ( t ) + 3-i
.At.
(10)
248
S. NAOANOet al.
We must note, however, that the path of a vortex near the separation point is usually highly curved. Therefore, if we calculate the movement of the general vortices at a time step At, then the nascent vortices should be calculated at a much smaller time interval. Actually, the movement of the nascent vortices is pursued at At~5 (Sarpkaya[7]). The selection of the time step At is very important in a practical calculation. The accuracy of the calculation becomes higher, while the computer time required becomes greater, when a smaller time step is specified. The same argument also applies to the time interval (AT) between the introduction of the nascent vortices. After a large number of preliminary calculations the following time step and interval are employed throughout the calculations discussed in this paper.
At=O.OS(2b/U) and AT=2At. 3.3 DiMusion of vorticity As has been pointed out at the beginning of this paper, the vortices shed from the front edges of the rectangle have to go a long way along the surface of the body. Therefore, the interactions between the vortices and the body will play an important role in the present flow. Special attention must be paid to the treatment of vortices when they come close to the body. The usual assumption of removing vortices that have approached closer than a pre-fixed distance to the surface of the body does not seem to give a reasonable model in the present flow. In this connection, it is noted that a potential vortex induces the velocity whose magnitude is in inverse proportion to the distance from its center. While, the velocity q induced by a vortex in a viscous fluid at rest is given by
where v is the kinematic viscosity of the fluid, t the time, and r the distance from the center of the vortex (Lamb [l l]). q increases with r from zero to a maximum (at r = r,) and then falls asymptotically to zero. The distance r, is a function of time and viscosity, and ff the vortex is concentrated in a point at t = 0, the viscous core radius r, is given by r, = 2.242X/-~. The difference in the velocities induced by a potential vortex and a viscous vortex is slight when r is larger than r,, but it becomes serious as r gets smaller than r,. In other words the potential vortices tend to induce unrealistic velocities when they come too close to the surface of the body or too close to each other. It seems reasonable, therefore, to let each point vortex bear a characteristic radius ro which varies with the kinematic viscosity of the fluid and the time that has elapsed after the shedding of the vortex (i.e. ro =/3r,,/3 being a constant), and to remove from the calculation the vortices that have approached closer than the radius r0 to the surface of the body. Physically the removal may be regarded to correspond to the destruction of vorticity in the shear layer by interaction with the boundary layers along the body. Also, a second characteristic radius rfi (r~ =/3'r.,/3' being another constant) is introduced, and any two vortices are coalesced to an equivalent single vortex when their separation is smaller than the sum of the second characteristic radii in order to avoid the excessive induced velocities to each other, as well as to reduce the total number of vortices in a smooth way. Thus, the distances for the destruction and the coalescence of vortices are assumed to be proportional to the viscous core radius of each vortex, /3 and /3' being the constants of proportionality. It may be worth noting that the diffusion of vorticity in the separated flow is due not only to the kinematic viscosity of the fluid but also to the apparent viscosity due to turbulence. However, the correction to the viscosity may be absorbed in the selection of the values of/3 and/3'. Unfortunately, there is no rigorous way of determining the proper values of /3 and/3' on a theoretical basis. Hence, after carrying out a number of preliminary calculations on flow past a square prism, the constants/3 and ~' in the present study are assumed to be unity as the first approximation, that is,/3 =/3' = 1 and hence ro = rfi = r,. In order to see the potential ability of the present diffusion model of vorticity, separate numerical calculations have been carded out on flow near the leading edge of a semi-infinite flat
A numerical analysis of two-dimensional flow past
249
Re-lO
Re-103]
Re-t0
S"
°" f: oo
°0%°°0
~ % ° ooo- o~
o
0
0
0
0
0
""
o
%
° 0~
Fig. 2. Leading edge separation of a semi-infinite thick flat plate. All the results are those at T = 14. The recirculating flow region is still in the course of development at this time at Re = 106.
plate. Some of the results are shown in Fig. 2, where the Reynolds number is based on the thickness of the plate. Figure 2 is not quoted for quantitative discussions, but is shown just as an illustrative example. For instance, at low Re, the diffusion model as it stands will only poorly represent the actual diffuse vorticity field that would exist in the real fluid flow. Nevertheless, qualitative characteristics of the flow may be seen: When Re is small (Re ~ 102), the flow does not separate. At a little higher Re (Re ~ 103) only a short separation bubble is formed. With further increase in Re the vortices begin to make a recirculating motion in the separated flow region and this becomes greater as Re is increased. It may, therefore, be said that, qualitatively, at least, the diffusion model has some ability to simulate the interactions of the shear layer with the surface of the body. 3.4 Loss of vorticity It is widely known for various bluff bodies that the amount of vorticity found in the fully formed street vortices is some 60% or less of the amount generated in each shear layer before separation (e.g. Fage and Johansen[12]). There seem to be three main causes of the loss in vorticity. They are (i) the destruction of vorticity near the wall due to the interaction with the boundary layers along the body, (ii) the cancellation of opposite-signed vortices due to mixing and (iii) the dissipation of vorticity due to three-dimensional turbulence. In the discrete vortex representations of the flows tried by various research workers, some mechanisms for the loss in vorticity due to (i) and (ii) have usually been included, but the loss in vorticity in these studies has been much less than that observed in the real fluid. Our preliminary calculations past a square prism, where the vorticity loss due to (i) and (ii) were included in the manner stated in the previous section, revealed that the vorticity clusters formed just behind the body retained so strong a vorticity as to let many of the upstream elementary vortices hit the longitudinal faces of the prism, thus preventing the periodic shedding of vortices. It seems necessary, therefore, to include some extra mechanism for the loss in vorticity in our discrete vortex representation.t
tThe outline of this paper was presented at the Euromech Colloquium held at Imperial College in July 1979 (Bearman and Graham[18]). After the Colloquium, Prof. Sarpkaya kindly sent us copies of his two papers. We found in one of his papers that he and his colleague adopted a vorticity loss model which would work in nearly the same manner as ours (Sarpkaya and Shoaff[8]). Also we found in the other paper by him that, parallel to our application of the discrete vortex model to an oscillating rectangular cylinder (see Section 4), he had been applying his model to an oscillating circular ¢yfinder (Sarpkaya[19]).
250
S. NAOANOet al.
Presumably, the extra loss in vorticity in the real fluid is due mainly to (ih'), whose precise mechanism in the separated flow, however, is far from being known as yet. What we can do at the moment is to incorporate into the model an artificial loss of vorticity in a phenomenologlcal way. We, therefore, carried out an experiment in a low turbulence wind tunnel and measurements were made of the change in the strength of vortices shed from a square prism. The vorticity flux was evaluated from the time-mean velocity measurement at various stations from the front to the rear edges of the square, while the strength of downstream vortices was estimated in a manner analogous to Fage and Johansen[12]. Figure 3 shows the estimated variation in the strength of vortices with the streamwise distance. The strength of vortices decays only minimally just after the front edge separation (-0.5 _-
u
~'0.5 x
I
5,0
J
IO,O
I
15,0
x/2b
Fig. 3. Space-wise loss in vorticity obtained in the experiment for a square prism at Re = 2.5 x 104.
1.0
8
~-o.s
T)
T2
T
Fig. 4. TL,ne-wise loss in vorticity incorporated into the numerical model. F(TO/F(T =0)=0.99 and F(T2)/F(T = 0) = 0.61 where T1 = 2 d U and T2 = T1 + 8blU.
A numerical analysis of two-dimensional flow past
251
T ffi4 (half a period), all quantities are nondimensionalized hereafter by the far upstream uniform velocity U and the length of the lateral side of the rectangle 2b. Unless otherwise stated, the Reynolds number (Re ffi 2bU/u) of the flow is set at 2 x 104 in the following calculations. It must be noted that the Reynolds number of the flow appears in the model only through the diffusion process of the vortices. Although the Reynolds number, through which some qualitative effects of viscosity are taken into consideration, will not necessarily correspond rigorously to that of the real flow, the effects of changing the Reynolds number of the flow in the model will be examined later. 4. RESULTS AND DISCUSSIONS
The numerical calculations are carried out for four types of stationary rectangular prisms with a c/b ratio ranging from 0.5 to 2.0, and for a square prism that is in forced vibration at various frequencies. The development of the wake behind a stationary square prism after the impulsive start from rest is shown in Fig. 5; The small circles denote the positive vortices shed from the front upper edge, the x's the negative ones from the front lower, the triangles the positive and the crosses the negative ones from the rear upper, the y's the negative and the diamonds the positive ones from the rear lower edges, respectively. These marks do not represent the strength of the elementary vortices but just show their positions. The small amplitude vibration of the prism for half a period at the early stage (T < 4.0) is seen to be large enough to amplify asymmetry in the flow. One of the interesting points is that some of the vortices shed from the front edges get close to the longitudinal faces of the prism and then go upstream back into the dead flow region. Recirculating flow regions are formed between the shear layers and the longitudinal faces. Many of the vortices shed from the rear edges are also swept back into the recirculating flow regions.
T- 1 ~
I
I::?::: :~:;:,,
~
~`:
"
T" 2 ~'"~:
T" 3 ,~,,,¢J~"
T" 9 ~.~7;, ~
C.",', ...~,=
i.
T=11
l"l .7.,,~ ~
b
T-6
=x
Fig. 5. Growth of wakes behind a square prism afterthe impulsive start from rest. Numbers in the figure denote the non-dimensional time after the start.
252
S. NAGANOet al.
Figure 6 shows the fully developed quasi-steady periodic states of the evolution of the wakes during nearly a full cycle of computation for the four types of stationary prism. The double circles and the thick x's denote the centers of positive and negative vortex clusters that are replaced to equivalent vortices. The periodic shedding of vortices and the formation of the vortex clusters are quite clear when c/b is smaller than or equal to unity• The periodicity is a little obscured at c/b = 2.0. This point will be discussed again later. Figure 7 compares the enlarged picture of the typical vortex arrangement near the longitudinal face of the square prism with the separation streamline estimated from the mean velocity measurement in our experiment. The curve (A) in Fig. 7 shows the outer edge of the shear layer that was obtained in the experiment. The curve connects the positions at which the maximum time-mean velocity was detected in the traverse of a hot-wire at various stations from the front to the rear edges of a square prism. Curve (B) connects the positions where the time variation of the instantaneous velocity was the largest. Therefore, curve (B) may be regarded as to represent the time-mean position of the shear layer in the experiment. The vortex sheets formed by the elementary vortices in the numerical calculation are seen to well simulate the separation streamline in the real flow. The rate of shedding of vorticity from the four edges is shown in Fig. 8 as a function of time for the four stationary rectangular prisms. After the initial period of flow est~iblishment, the rate (a)
(b) ~,,..++
T-32
~+
~.'~.,
~ .;~
%'.+',~+
";t'.
.i'
,+
XM~ ~;.,: ...... ~P.a ~, .y,r.."..,
+,: ..... :..:t
;,,
x
~ _ + p ~ . . , ?;;
{4"~ •
-,,
~" + -
%
,v.+
Fig. 6. Vortexarrangementsovernearlya-fullcycleof steadyperiodicflow.(a), c/b = 0.5; Co),clb = 0.62.
A numericalanalysisof two-dimensionalflowpast • (C)
253
(d)
• • ~,
• ~,.,;,~
% ~
":~"
x
~-32
....~, ":.
:
.~,~,:~
.'~,' . ,
."
•
T-36 ~ ' k ~ , ~ . :
~ %":~o.
.~.. .%
.::'~:,
o
,
X
2"
."."~"';~i~,~.~...... T-38
"'~d~'"
",,.
f
Fig.6(Contd).Vortexarrangementsovernearlya fullcycleofsteadyperiodicflow.(c),clb = 1.0;(d),clb = 2.0.
is seen to vary periodically except for the case with cJb = 2.0. It is also quite clear that the rate of shedding of vorticity from the rear edges is considerably smaller than that from the front edges. This is not surprising since the rear edges are in the dead flow region and the flow velocity near them is usually small. It may also be noted that the periodic shedding is established in the shortest time when db is equal to 0.62. In Fig. 9 the time variation of the fluctuating lift coefficient is plotted. Here again the periodicity is quite clear (except for the case with cJb = 2.0) and these plots are used to determine the Strouhal number of the vortex shedding in each case. Since CL is plotted against the nondimensional time, the Strouhal numbers for (a)-(c) can be calculated as an inverse of the average period of the last few full cycles of the fluctuating lift. They are tabulated in Table 1. In accordance with the fact that the periodic shedding of vortices and the formation of vortex clusters are somewhat obscured in Fig. 6(d), where clb -- 2.0, the time variation of CL in Fig. 9(d) cannot be said to be fully periodic either. Therefore, a Fourier component analysis has been made of the unsteady lift to find the Strouhal number for this particular case. And it has been found that there are two dominant Strouhal numbers as shown in Table 1. In this connection it is noted that Otsuki et al. [13] also found in their experimental study two dominant CAF
10:4
- C
254
S.
T-26
NAGANO
a
(A) ,~ ,, x
~
,
~ j ~~ ~ ( B )
~'-°
,
T 28
(A ~o(Ao o
T-30
~.~----¢--- (B) ~o~-~ ~ o OoO ° ,
aL
°~
:
~ "
(A) o
(A)
o o o ~o~
o,
'
(B)
~
~
(
B
)
-'C;-.'~"~o- #
o~ ~
-
Fig. 7. Enlarged picture of vortex arrangements near the longitudinal faces of a square prism at various instants of time. ©, positive vortex near the upper face; x, negative vortices near the lower face. The positions of the negative vortices are shown with their y coordinatenegated.
X
2.0
00
o xx
(a)
O0 ~'0000 0 Xx
"o
"o 1 , 0
+
oo o ~ O"'vvo^
Or~"-~xx x X X x
xxXXx
0
xxx 0
xo
x ~^
o x0 ,-,x
xxxln,l~Sr.,n x,XOx 'x'0 OXx x v -~ -x ~00000,., x Xxxx x 0o00 xx x 0oO x x x
A
10
20 x x o x ~
(b) ~:'
+-'
x
~
0 x x x 00o00 O
"~o ] . , 0 + ,
~O Oo Vx
xXx O x xO
XxxxX
j~÷+_
0O°o
xy O 0x
0O0o0
.
.--
,
xx O
^o
--
o
x~ O
O^
x
~" "~ 1 . 0 --
o6~^
, ~,,_,
- ~ --OUO O
¥Ooo~oOXxx XxxX A
& + A& & + _ +~+ 4- ~ A~A÷ ~&~A.~;
0
~';+ ~"
+A~
-+
,
o,.,x
X
X XxX
^
0
O
°o°
,-,
xx ~)
x 0
XxxxX
o
x
~
~
~,A +÷ +4.
A 4-
x
°° o
x^O
Ox
Q^ x ~,.,,
o~X
O
x
xxX
vo00
"
,t
+++
T X
_
0
0
~o
O"
& 4.+ 4. A A 4. ~ A & ++& + -A A
0°o°
O
oOOo
O~ X __O UO ~ "
~
x
~J xx x O
30 x~ X
,
x
+ ++ " + - +
x
~÷+ez~.9.+$+,_a.+.+"^,,a,~
20
x x~'
0 x
x~O^O~XX x O O x x O ,.OOoO + Xxx x Oooo xxx uO
.&A . ~÷÷-14¢ A &
x_*~,~
oo
XX
x
T
~x
xO0
+÷++++~+
,~ X ~
x,.,
A
i0 2,0
xO x
20 ooo
x x,.,oU
++ ~A
0
30
Ox
Xxxx x "0Oo0 XXx
(C)
X ~ X ! WoX..xxXXXX
"~ "-.
0 0oO
^O°oo Xo v Ox
4.
o x
x ^~
10 2.0
xx
"
2.0 --
0 xo
A &
"0
Xx
o x
X
x
OXX O
Ooo+
OX
x O xxxX
O
,,
+ + 4-
A A .A&
*
*++++~,&~M.a.q..e+*a;÷+, T
30
(~)
0o 0
4.0
XXxOO x
x~ 1,0
o
X X
÷+÷ A A A + A A A4. ÷ ÷ & + +A
"
+ +
A
t+-+~=-~-~++M~-~+++~,'N~***~"~+*~* ~--~~*~-~t*~;+'~;* * 10 20 30 T
Fig. 8. Rate of shedding of vorticityfrom the four edges. (0 front upper, x front lower, A rear upper, and + rear lower edges)
Strouhal numbers for a stationary prismatic bar with rectangular cross section when the ratio c/b was in the range from 2.0 to 2.8. One was on the branch that was connected to the lower clb and the other on the branch connected to the higher clb. Another point which should be noted in Fig. 9 is the amplitude of the unsteady lift. It increases gradually as clb is increased from 0.5 to 1.0, will" e it tends to decrease when c/b is increased to 2.0. The root mean square amplitudes of the unsteady lift coefficient are tabulated in Table 1.
255
A numericalanalysisof two-dimensionalflowpast
2.0
(o)
0
o% _ o °~°°°°'~-
i0
20 ,.,°°°°o
_ ~°°°o°°°
,-°°°%., o-.°° Oo o 3() Ooooooo T
%ooo°
oo
-2,0 2.0
(b) oo
J 0
~0o 0000
n
°°°%;
°ooooO~
oo° 0
=
i0 o
0000
0
°° '3
o
o
0
20,o O00
0o
on
oOn -o
oI n 30000°0
°OoooO
o
o°°%
o T
-2.0
.01
(c)
ooo
o
oO°°O
i°0io
o
o°O o
o
°°oOoooO
300
0000 o
oo 0
-2.0
0
o
o
o
o
o° T
oo o
o
o
0000
(d)
0
o
o
oooo OooOo
2.0
o
o
10,9 v o o O0 0 Ooo
O0o
000 o, o 20 oo
o
o 0
0 O0
o
0%
0
0
oi
oo
o
0030
Oo ° ° ° o o
TOo
Oo
0
o° ° o
Fig. 9. Time variation of fluctuatinglift coefficient.(a), c/b = 0.5; (b), clb = 0.62; (c), c/b = 1.0; (d), clb = 2.0. Table 1. Resultsof the computationsfor the four types of stationaryrectangularprism. St denotesthe Stroulminumberand an overbardenotesthe timeaveragedvalue. C(x,out0t ion
c/b
St
~/(CL)
Co
(0)
0.5
0.161
0.584
2.2q7
(b)
0.62
0.153
0.797
2.253
(c)
1.0
0.142
1.173
2.059
(d)
2.0
0.102 0.20tl
0.916
1.181
The time variation of the drag coefficient is shown in Fig. 10. The mean level of CD is seen to be the highest when clb is around 0.5-0.62. It decreases gradually as clb is increased further. The time mean drag coemcient is also tabulated in Table 1. The period of the cyclic time variation of Co is clearly halved compared to that of CL when clb is equal to 0.5 and 0.62, but the behavior is not obvious when c/b is equal to or greater than unity. If the Fourier component analysis were made of the unsteady drag, the half period variation in CD would be detected for the latter cases too. Figure 11 compares the results obtained in the present numerical analysis with those of the experiments by various research workers. The top figure compares the Strouhal number. Although the Strouhal numbers obtained in the numerical calculations are generally about lO~ higher than the experimental ones, the general trend of the experimental results is well predicted. In the middle figure the root mean square values of the fluctuating lift coefficient are shown. Though not many experimental results are available, the agreement of the analytical
256
S. NAGANOet al.
3.0 (O)
2.0
o
o
0
0
o
o oo oOOO o°° o o°°Ooo° °oo °o o 000OoOooooooOooo%00oO0ooo o OoO oo
1.0
o
o° °o
°°°oO
Oo oO o
0
ooOOC~ooOOooOoo° 0
o
oOo
OoO
o
oOoooo °
Oo°°° °Oo
o
oooOOO o
o
oo
°o
o
ooo
o °°
oo° °o°
(b)
oo°°Oo
0
00 0
0~0 0
o=oO 0 0°°°000000°0°°=°°°°%°0 °O°°o O 0 0o~
O0
o
°o°o0 0 0
(C)
0
0 O
0
0%0oo0o°00o° _ ..,.,..oo _~ ^o°° 0 0 0 vv ~r~ 0 ou Oo ^o O 0 o o ~ o ~ °°oOOO °°o o °°o°°°°°° (cl) o" o oo -',.,o Oo i
1
i
i0
20
30
Fig. 10. Timevariationof drag coefficient.(a) c/b = 0.5; (b), db= 0.62; (c), db = 1.0;(d), c/b = 2.0.
results with the experiment is quite good. It is interesting to note that the root mean square value of the fluctuating lift coefficient seems to take a maximum at c/b between 1.0 and 2.0. As is seen
in the bottom figure, the mean drag coefficient obtained in the present study predicts the general trend of the experimental result quite well; although, the sharp peak in CD near clb = 0.6 has not been detected in the limited number of numerical calculations. Therefore, it may be said that the present model as a whole, is able to represent the separated flow past a rectangular prism. A further supplementary study has been made to see the effects of changing the Reynolds number of the flow in the present model. In the computation for the stationary square prism at Re = 2 x 104, the Reynolds number of the flow was suddenly changed to a different value at T = 35, the computation was continued further, and the results were compared. The result at Re = 2 x 105 was almost completely coincident with that at Re = 2 x 10~. When Re was reduced to 104, a slight difference was seen in the time variation of CL but the periodic pattern as a whole was similar to that at Re = 2 x 104. No obvious difference was found in the downstream vortex arrangements either. However, the result at Re = 2 x 103 was quite different from the other three cases. The Karman vortices were not formed and the periodic vortex shedding was not sustained. The radius of an elementary vortex grew fast at this Reynolds number and too frequent combination and destruction of vortices took place just after leaving the separation points. The consequence was the failure of the model to represent the continuous shear layers by a sufficiently large number of elementary vortices. The situation might have been improved by employing a larger time interval AT between the introduction of the nascent vortices to avoid the frequent combination of nearby vortices. We did not, however, try a detailed examination of this point. It may be worth noting that the time step At has to be made smaller if Re is greater than, say, 106. The reason is that, when Re is higher, the core radius of a vortex remains smaller. The vortices can then approach closer to the body or to each other. Therefore, the
A numerical analysis of two-dimensional flow past
257
0.2
0,1
I
,
I
1,0
c/b
2.0
c/b
2,0
1.0
,
0
!
,
I
1,0
0
3.o I
.e/J~."
\.
2.0
1,0
T o
. i
,
1,0
!
c/b 2.0
Fig. 11. Comparison of the present results (shown by O) with experiments by Bearman a al.[14] (A, Re = 2 - 7 x 104), Nakaguchi et al. [15] ( x, Re = 2 - 6 x 104), Nakamura et al. [16] (V, Re = 3.3 - 13 x 104), and Okajima et all17] (12],Re = 6 - 1 2 x 104).
time step At should be made smaller, as At ~ l/~/R-e to be rigorous, if the accuracy of the calculation is to be kept constant. To sum up the results of the supplementary study, the present model with the parameters given in the paper is suited for flow at Re in the order of 104- 10s. The order of Re can be raised easily by adopting smaller At. The present model is also applicable to the analysis of flow past a vibrating rectangular prism.'r To show this by examples, the numerical calculations are carried out for the square prism that is subjected to a sinusoidal transverse vibration in a plane normal to the uniform flow. The amplitude of the vibration h is fixed at 0.1, while the nondimensional frequency of the vibration So (= 2bf/U) is varied in the range from 0.125 to 0.155. An example of the vortex pattern formed behind a vibrating square prism is shown in Fig. 12. The lateral distance between the positive and negative vortex streets is seen to be much smaller when compared to the vortex patterns behind a stationary prism. Figure 13 shows the time variation of the fluctuating lift obtained for the square prism in forced vibration at various frequencies. The displacement of the prism Y is also shown by the thin broken line. It is known that, when a bluff body in a stream is forced to vibrate, the vortex shedding frequency is locked to the frequency of the forced vibration if So is close to the natural vortex shedding Strouhal number St of the non-vibrating body. According to the experimental work by Otsuki et a1.[13], the locking range for a square prism in forced tSee footnote on p. 249.
258
S. NAOAN0 et al.
0
x
Fig. 12. Typical vortex pattern
2.0
behind a
vibrating square prism,
(e)
d
o
°o
O0
o
o°°
o'---.. . . . - "
o
000
0
0
2.0 I
(f)
-o O0o
o°°o
o
o
o
01=OOoO;,oo,0o:,o
_.y
~-- "-''-... /-" "-'"",..o
"" ' L
|
000
.... "
0
o
2~,,.~ "-. o
0
oo
0000
,oo
:o ..j "'--
o .... "
30
0000
0
°°°
-2, O I" 2,0
(g)
J
Oo
o
"P°o o ° / : ~ ' % o 10 . / ; , ~ ' , . - ~ --o.~"" ~ "-~ / °o° "oo
0 o
°°°°
---',
0
Y
...-."d°o
(h)
2,0 o
_
_..,,,_,,,,~
..o
0
~'--JY'o
o
O0
-, 30
._.
./
o
"'R
o°"J'-'l;
T" O0000 °
0O
o
000
00
oOoO o o
o
o
o
00
0 O0
-2.0
0
0 00
o / ....
""--% ° ' ' / -ZU " o.,
0
.... " T
°°°
00 0
T
"----
"
0
00000
000 /
°°o°
30,,........~
"'-- . . . . .
00
00
o
°°
0
-2.0
(J
S~ = 0.135.
o°°° °°°°° °°°° Y o °° ,-~-.....0"-..,i 0 o ....o..... .. 2,0 ,,,.I.......~o
o "-..o
0
0
o Y
o o
.
"~-~"
°ooO
°°-"'~'
20
°°°o
o o
o
" .... "
o
o o Ooo
o'°-J30
o
",z_ oo
oo
-2.0 Fig. 13. Time variation of fluctuating lift on a vibrating square prism. (e), S~ =0.125; (f), S~ =0.135; (g), So = 0.147; (h), S, = 0.155; ---, displacement of the square prism.
transverse vibration with a nondimensional amplitude 0.1 was 0.90_< SJ St -<1.08. Figure 13 shows that the phase difference between C,. and Y varies with the time at (e) SJSt = 0.88 and (h) SJSt = 1.09, i.e. the locking-in does not take place in these cases. While at (f) SJSt = 0.95 and (g) SiS, = 1.04, the phase differences approach asymptotically to constant values, and the vortex shedding is locked to the v~ration of the prism. F~lhermore, it is clear that the fluctuating lift lags behind the displacement of the prism when S~ is larger than St as in (g), while CL is ahead when $, is smaller than St (f). These phase differences between CL and Y in the locking-in regime are in good qualitative agreement with the experimental observations. It will be noted that the vorticity loss model shown in Fig. 4, which was assumed based upon the measurement of flow around a stationary square prism, has been conventionally approximated in the calculations of the flow past the vibrating square prism. There is no obvious reason, however, why the vorticity loss of vortices shed from a vibrating body may be assumed to be the same as that for a stationary one. As the vorticity loss has large influence on the magnitude of fluid forces acting on a body, a further study on the mechanism of loss in vorticity behind a vibrating body seems to be necessary in order to make more detailed quantitative discussions on the fluid forces acting on a vibrating body. This point is left for the future study.
A numerical analysis of two-dimensional flow past
259
5. CONCLUDING REMARKS
A discrete vortex model has been established to analyse the two-dimensional separated flow past a rectangular prism. The flow field around a rectangle has been transformed to the exterior of a circle via elliptic functions. The effects of the viscous diffusion of vortices and the loss in vorficity have been incorporated into the model. The diffusion model seems to have some ability for simulating interactions of the shear layers with the body. It is also effective in enabling the formation of the recirculating flow region along the surface of the body as well as in avoiding excessive induced velocities. Perhaps the most critical part of the model lies in the incorporation of the artificial loss of vorticity. It would be much better if the loss of vorticity could be incorporated into the model in a form of a general turbulent dissipation concept. However, with a careful choice of the parameters stated in the paper, the present model as a whole has been shown to have a si~ificant ability for representing separated flow past bluff bodies even when the interactions between the shear layers and the bodies play an important role in determining the characteristics of the flow. The numerical computations have been carried out at the Computer Center at the University of Tokyo. The authors would like to express their gratitude to Mr. Y. Machida who carried out the experimental work.
REFERENCES I. A. J. Chorin, Numerical study of slightly viscous flow. J. Fluid Mech. $7, 785 (1973). 2. R. R. Ciements, An inviscid model of two-dimensional vortex shedding. J. Fluid Mech. $7, 321 (1973). 3. R. R. Clements and D. J. Maull, The representation of sheets of vorticity by discrete vortices. Prog. Aero. ScL 16, 129 (1975). 4. M. Kiya and M. Arie, A contribution to an inviscid vortex-shedding model for an inclined fiat plate in uniform flow../. Fluid Mech. 82, 223 (1977). 5. K. Kuwahara, Numerical study of flow past an inclined fiat plate by an inviscid model. J. Phys. SOc. Japan 35, 1545 (1973). 6. K. Kuwahara, Study of flow past a circular cylinder by an inviscid model. J. Phys. Soc. Japan 45, 292 (1978). 7. T. Sarpkaya, An inviscid model of two-dimensional vortex shedding for transient and asymptotically steady separated flow over an inclined plate. J. Fluid Mech. 68, 109 (1975). 8. T. Sarpkaya and R. L. Shoaff, An inviscid model on two-dimensional vortex shedding for transient and asymptotically steady separated flow over a cylinder. A/AA Y. 17.11, 1193 (1979), 9. H. Kober, Dictionary of Conformal Representations. Dover, New York (1957). 10. H. Schlichtlng, Boundary Layer Theory. McGraw-Hill, New York (1960). 11. H. Lamb, Hydrodynamics. Dover, New York (1945). 12. A. Fage and F. C. Johansen, On the flow of air behind an inclined plate of infinite span. Proc. Roy. Soc. A, 116,170 (1927). 13. Y. Otsuki, K. Washizu, H. Tomizawa and A. Ohya, A note on the aeroelastic instability of a prismatic bar with square section. J. Sound Vib. M, 233 (1974). 14. P. W. Bearman and D. M. Trueman, An investigation of the flow around rectan._gulat cylinders. Aero. Quart. 23, 229 0972). 15. H. Nakaguchi, K. Hashimoto and S. Muto, An experimental study on aerodynamic drag of rectangular cylinders. J. Japan Soc. Aero. Space SOl. 16, 1 (1968). 16. Y. Nakamura and T. Mizota, Unsteady lifts and wakes of oscillating rectangular prisms..A. SOc. Civil Engng 151, EM6, 855 (1975). 17. A. Okajima, T. Mizota and H. Takata, An experimental study on the pressure distn~oution along the surface of transversely oscillating rectangular cylinders in a uniform flow. Preprint lor JSME Meeting No. 774-13 (1977). 18. P. W. Bearman and J. M. R. Graham, Vortex shedding from bluff bodies in oscillating flow: A report on Euromech 119. J. Fluid Mech. 99, 225 (1980). 19. T. Sarpkaya, Vortex induced oscillatlons--a selective review. J. Appl. Mech. 48-2, 241 (1979).