- Email: [email protected]
Numerical simulation of the wake flow behind trapezoidal bluff bodies
Computers & Fluids 31 (2002) 99±112
www.elsevier.com/locate/comp¯uid
Numerical simulation of the wake ¯ow behind trapezoidal blu bodies R. Kahawita *, P. Wang Department of Civil Engineering, Ecole Polytechnique de Montr eal, Montr eal, Qu ebec, Canada H3C 3A7 Received 22 November 1999; received in revised form 4 October 2000; accepted 10 November 2000
Abstract Two-dimensional numerical simulations of the Benard-von Karman hydrodynamic instability behind trapezoidal blu bodies has been studied using the spline method of fractional steps. For lower Reynolds numbers, about Re=Rec < 2 (Rec is critical Reynolds number), numerical results con®rm the experimentally observed behavior reported by Goujon-Durand et al. (Phys. Rev. E 1994; 50:308), i.e. the maximum amplitude of the velocity component oscillating with the fundamental frequency follows fairly well the scaling 0:5 0:5 law Amax Re Rec and the position Xmax Re Rec . The in¯uence of the trapezoidal shape on the value of the critical Reynolds number and on the vortex shedding is brie¯y discussed. It appears that the in¯uence of the trapezoidal height H is the dominant in¯uence on the value of Strouhal number when compared with the eect of the smaller trapezoidal base width B. Ó 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction The wake ¯ow behind an obstacle has been extensively investigated numerically and experimentally, most classical work having focused on circular cylinders. A variety of vortex shedding characteristics have already been discussed in the literature in the Reynolds number range 50±200 (see, for example, the review paper of Berger and Wille [2]). In recent years, the global nature of the wake ¯ow has attracted much attention [3]. In a recent paper, Dennis et al. [17] have conducted a comprehensive investigation of the temporal development of the viscous incompressible
*
Corresponding author. Tel.: +1-514-340-4711, ext.: 4506; fax: +1-514-340-2981. E-mail addresses: [email protected], [email protected] (R. Kahawita).
0045-7930/02/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 1 5 - 9
100
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
¯ow induced by a circular cylinder. The cylinder is assumed to be started impulsively while simultaneously performing time-dependent rotational oscillations about its axis. Solutions were obtained using a series expansion for small times and a spectral ®nite dierence technique for moderate times. Good comparisons with previous experimental results were obtained. For industrial ¯owmeters that exploit the relation between the shedding frequency and the ¯owrate or velocity, a trapezoidal shape is found to be more desirable than a cylinder, since a well-de®ned vortex emission due to clean separation at its sharp edges is assured. This is the subject of the present paper. Flow of a ¯uid past a trapezoidal shape is a more complicated phenomenon than that observed around a circular cylinder. Besides the Reynolds number, the vortex shedding depends strongly on the leading dimensions of the trapezoidal body. If observations are made in the far ®eld region, it would be dicult to determine whether the vortex shedding is due to a cylinder or a trapeze. In the near ®eld region however, the ¯ow structure and pressure ®eld is signi®cantly dierent. This is somewhat analogous to the ¯ow pro®les obtained in a jet from a square ori®ce. In the near ®eld region, information on the exit geometry is still retained by the ¯ow pattern. In the far ®eld however, the jet approaches circular symmetry (``memory is lost'') and the ¯ow closely resembles what would be obtained from a round ori®ce. Goujon-Durand et al. [1] recently reported the results of an experimental investigation into the vortex shedding from a trapezoidal body. They obtained the scaling laws for the evolution of the global mode describing the envelope of the peak to peak amplitude velocity oscillation in the wake ¯ow downstream of the body. They were also able to indicate the existence of a universal curve for the renormalized amplitude of the streamwise evolution of this amplitude. More recently, Zielinska and Wesfreid [4] presented some numerical simulations of two dimensional wake ¯ow behind an equilateral triangular obstacle using the spectral ®nite element commercial code NEKTON. Their numerical results con®rmed the experiments of Goujon-Durand et al. [1] by showing that the saturation amplitude of self sustained oscillations in the wake has a well de®ned maximum Amax and corresponding downstream location Xmax which are functions of the Reynolds number. Amax diminishes and moves further downstream as the critical point is approached. Since the pioneering work of Rubin and coworkers [5±7], the cubic spline integration technique and its applications in heat transfer and ¯uid ¯ow have developed rapidly. The present investigation is devoted to the numerical simulation of the two dimensional wake ¯ow behind a trapezoidal body using the spline method of fractional steps (SMFS) [8] which in discretized form is directly applicable on a non-uniform mesh in an x±y coordinate system. The advantages over using a standard ®nite dierence procedure is that no coordinate transformation is necessary, solutions are obtained directly in the physical plane. Furthermore, any arbitrary shape may be treated since computational points may be placed directly on the surface of the body. Since the trapezoidal shape has not yet been explored numerically, it is hoped that the results will provide insight into the experimentally observed behaviour reported by Goujon-Durand et al. [1]. The in¯uence of the trapezoidal shape on the value of the critical Reynolds number and on the vortex shedding is brie¯y discussed. The numerical results indicate that the shedding frequency initially increases linearly with Reynolds number (lower Reynolds numbers). The frequency then tends to level of until a maximum is reached. Surprisingly, after attaining this maximum, the shedding frequency decreases. This behaviour is consistent with Okajima's experimental results [9] for rectangular cylinders, however in this paper it is not discussed in detail. This study may thus be
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
101
considered as contributing to the present state of knowledge on this type of ¯ow as well as contributing to the state of the art in numerical algorithms for ¯ow simulation. 2. Governing equations The non-dimensional equations in stream function and vorticity form may be written: r2 W
X
1
oX oX oX 1 2 u v rX ot ox oy Re
2
with r2
o2 o2 ox2 oy 2
and u
oW ; oy
v
oW ox
3
where x is the longitudinal coordinate in the ¯ow direction and y is the transverse coordinate. Here W and X denote the dimensionless stream function and vorticity respectively while u and v are the respective dimensionless longitudinal and transverse velocities. Re U1 D=m, D is the trapezoidal width. 2.1. Boundary conditions On the trapezoidal surface: u v W 0; X
r2 W
4
At the in¯ow region x u 1;
v0
A !
1:
and X 0
5
On the out¯ow region x A ! 1: oW oX 0; and 0 6 ox ox The boundary condition on the vorticity in the out¯ow region corresponds to the free out¯ow of a ¯uid property, (here the vorticity). Extensive numerical experimentation in fact con®rmed that no discernible eect could be observed in the results as the out¯ow boundary was displaced further downstream. On the horizontal surfaces that de®ne the lateral limits of the computational domain, y b ! 1: u 1;
v0
and X 0
7
102
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
3. Results and discussion The SMFS schemes and the boundary conditions in discretized form may be obtained in direct fashion from the procedure detailed at length in earlier articles [8,10±12]. The development of the numerical algorithm will therefore not be discussed here. It is useful to indicate however that the spline discretized form is directly applicable to a non-uniform mesh in an x±y coordinate system, thus obviating the need for a computational coordinate transformation as is usually required with ®nite dierence schemes. The coordinate system and an illustrative mesh is shown in Fig. 1a and b respectively.
Fig. 1. (a) Coordinate system and (b) illustrative mesh used.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
103
In order to avoid computations with boundary conditions at in®nity, a speci®c ®nite domain was de®ned as follows: 5 6 x 6 30 and
10:5 6 y 6 10:5
i.e. the out¯ow domain was taken to be a distance of 30 units downstream from the rear base of the trapezoidal section. The time dependent non-linear coupled partial dierential equations were solved by considering a 121 121 grid. In order to accurately describe gradients that are expected to be steep in the boundary layer regions, a non-uniform grid in both the x and y directions was used. The essential advantage of the technique for the present problem lies in the fact that boundary conditions containing derivatives as shown in Eqs. (5) and (6) may be easily incorporated into the solution procedure since values of ®rst or second derivatives may be evaluated directly and maintain the same degree of accuracy when the algorithm that represents the spline approximation to the governing equations (1)±(3) is constructed. In order to verify the consistency of the present numerical study, the ®rst cases considered were those corresponding to existing results for wake ¯ows past a circular cylinder (see, for e.g., Ref. [13]) and behind an equilateral triangle [4]. In both cases very close agreement was obtained between the results of the authors cited above and the present computations. For example, for wake ¯ows past a circular cylinder, the value of Strouhal number at Re 49:4 is St 0:1241 (experiment of Williamson [14]) and St 0:1214 at Re 48 (numerical results of Dusek [13]), which compares favourably with St 0:1256 at Re 50 for the present results. The value of the critical Reynolds number Rec , for wake ¯ows past an equilateral triangle, is Rec 38:3 [4] and Rec 35 [15], while a value of Rec 37:5 was obtained in the present study. The accuracy and the convergence of the numerical solutions were veri®ed by mesh re®nement. All computations were performed on a PC. In the present numerical simulation, the blu body of trapezoidal section is characterized by its larger base width D the smaller base being perpendicular to the incident ¯ow. The following cases with dierent values of the smaller base B and the height H have been simulated: Case
A1
A2
A3
A4
B1
B2
B3
R1
R2
H B
0.7D 0.3D
0.7D 0.5D
0.7D 0.8D
0.7D 0.9D
0.3D 0.5D
0.5D 0.5D
0.9D 0.5D
0.7D 1.0D
1.0D 1.0D
Case R1 and R2 are particular cases, corresponding to rectangles with H 0:7 and H 1:0 respectively. In order to determine the critical Reynolds number, the computational strategy employed was to start the calculation procedure at a Reynolds number greater than critical and gradually diminish it towards the critical value. The principal problem encountered with this technique was that when the Reynolds number was very close to critical, convergence towards the stationary periodic state was very slow. The present Case A4 H 0:7D, B 0:9D approaches the experimental case (H 0:7D and the angle of the lateral faces being 4°) of Goujon-Durand et al. [1] in the laboratory. By extrapolation to zero of Amax , they obtained an estimate of the value of the critical Reynolds number
104
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Rec of about 58 (for the smallest blu body). The present numerical results indicate a Rec 42 with a corresponding Strouhal number St 0:11. For Re 65, i.e. Re 1:55Rec , the present numerical results indicate that the amplitude maximum occurs at a location Xmax of about 4:25D, Goujon-Durand et al. [1] obtained Xmax 4D or 4:5D. A typical instantaneous streakline view is shown in Fig. 2. The broken line in Fig. 2 corresponds to the envelope of transverse y values for which the amplitude along the ¯ow direction is a maximum. It is encouraging that this picture is very similar to the laser illuminated ¯uorescent view of vortex emission obtained by Goujon-Durand et al. [1]. Fig. 3 indicates the periodic variation of the transverse velocity v at dierent longitudinal positions x along the line of symmetry y 0 for Re 55 for Case A3. At a point very close to the smaller base, the amplitude A of oscillation is weak. As the value of x increases, the amplitude A increases until at a downstream distance of about x 5:5D it attains its maximum value; it then decreases as shown in the ®gure. It was noted that the value of the frequency at each x is constant. This agrees with the experimental results of Goujon-Durand et al. [1]. The constant frequency implies a synchronized oscillation and is one of the major features that characterizes the existence of a global instability in the full domain of shedding as observed by Goujon-Durand et al. [1]. Fig. 4 shows the periodic variation of the longitudinal velocity u at dierent y (transverse) stations along a vertical line x 5:4D for the same parameters as those of Fig. 3. u has a maximum in amplitude at about y 0:8D. As is well known for wakes behind blu bodies, u on the line of symmetry y 0 displays a frequency which is double that observed for o-symmetric lines. As the value of y increases from y 0 to about y 0:8D, the second harmonic transforms gradually into the fundamental while the amplitude of u increases to its maximum.
Fig. 2. Typical instantaneous streakline view.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
105
Fig. 3. Transverse velocity v at y 0 for dierent x positions (Case A3, Re 55).
Fig. 4. Longitudinal velocity u at x 5:1D for dierent y positions (Case A3, Re 55).
A time sequence of the streakline view for Case A3 with Re 55 is shown in Fig. 5. The von K arm an vortex street pattern and its transport by the main ¯ow is quite clear. The peak to peak amplitude of the oscillations of u, the longitudinal (x component) of velocity, was calculated from a time history of the oscillations at points in the wake along the x direction at dierent transverse stations (values of y). As an example, for Case A3 the amplitude of the
106
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 5. Time sequence of vortex street (Case A3, Re 55).
Fig. 6. Amplitude distribution of u in vertical plane at dierent x positions (Case A3, Re 55).
oscillations in u, for Re 55 from y 0 to 4D at dierent x positions are shown in Fig. 6. The following points may be noted: At each local ®xed plane of x, as y is increased i.e. moves away from the axis of symmetry, the amplitude of the oscillations increases steadily, reaches a maximum and then decreases. For example, at x 1:5D, the local maximum amplitude is about 0.061
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
107
Fig. 7. Envelopes of the peak to peak amplitudes of u for dierent Re (Case A3).
at a transverse position of y 0:7D. If the entire ¯ow ®eld is considered, the maximum amplitude Amax is about 0.278 and occurs at a position Xmax 5:5D and at a y value of 0:8D. This means that the maximum amplitude Amax occurs at about the same position X 5:5D for both longitudinal and transverse velocities. The value of Amax and the location at which it occurs is a function of the Reynolds number. Fig. 7 presents the distribution of the maximum amplitudes of the longitudinal velocity u as a function of Reynolds number throughout the ¯ow ®eld corresponding to Case A3. With decreasing Reynolds number, Amax diminishes and the position of Xmax moves farther downstream. An interesting point is that generally speaking, for a constant value of H, as the width B is increased, the critical Reynolds number also increases. For example, for a trapezoidal body of height H 0:7, when B 0:3D, Rec 34:5; B 0:5D, Rec 36; B 0:8D, Rec 39 while when B D, i.e. a rectangle Rec 40. For a constant value of B however, as the height H is increased, the critical Reynolds number decreases. For example, for a trapezoidal body with small base width B 0:5, when H 0:5D, Rec 35; H 0:3D, Rec 31. Fig. 8a and b present the relation between the Strouhal number, St, and the Reynolds number for a constant H with dierent B and for a constant B with dierent H respectively. As shown in Fig. 8a for a ®xed trapezoidal height H, at lower Reynolds numbers, the value of Strouhal number is almost independent of the small base width B. For example, when Re < 100, as the small base width B decreases from B 0:9D to 0:3D the values of St increase by 2%. As the value of Reynolds number increases, the curves de®ning the relationship between Re and St for dierent values of B gradually separate and then level-o to reach a ceiling quite rapidly. The experimental work of Okajima [9] however, has revealed that there is a fall-o from the maximum value as the Reynolds number is increased further. Clearly, it is quite possible that Fig. 8a would indicate this
108
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 8. (a) Relationship between the Strouhal number and Reynolds number for dierent values of B (H 0:7D) and (b) relationship between the Strouhal number and Reynolds number for dierent values of H (B 0:5D).
type of behaviour. In fact, it may be observed that Case A3 begins to demonstrate this fall-o. Another useful observation is that the larger the value of width B, the more rapid is the corresponding maximum attained and the lower the value of St as shown in the same ®gure.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
109
This probably implies that since the separation points occur at the ends of the trapezoidal base with the larger width, the in¯uence of the smaller base width B on the wake ¯ow interaction appears to become signi®cant only at the lower Reynolds numbers. However, for a ®xed value of the base width B, the value of the Strouhal number evidently increases with increase in trapezoidal height H as shown in Fig. 8b. The trapezoidal height therefore plays an important role in the wake ¯ow interaction mechanism. It may also be noted that the frequency of shedding increases linearly with Reynolds number at the initial stage; this linear relation is prolonged to about Re=Rec 2. This behaviour was noted in the experiments of Goujon-Durand et al. [1]. Further increase in Reynolds number causes the frequency to increase at a slower rate until it reaches a relatively ¯at maximum. After this maximum is reached, a decline in frequency occurs as shown in Fig. 9. This interesting behaviour was observed in the experimental studies of Okajima [9]. The envelopes of the global modes corresponding to various Reynolds numbers are presented in Fig. 10. By renormalizing A and X respectively with Amax and Xmax , the points of A=Amax versus X =Xmax for dierent cases with dierent Reynolds numbers collapse on to a single universal curve con®rming that the properties of the global modes are universal. It is noted once more that these properties are almost independent of the width B for the cases H 0:7D at lower Reynolds number. However for the case B1 (H 0:3D, B 0:5D) the dierence is obvious. The amplitude at downstream locations is clearly enlarged. Two scaling laws for the maximum amplitude Amax and its location Xmax for dierent cases at various values of Reynolds numbers is presented in Fig. 11. These results are in good agreement with the scaling laws Amax Re Rec 0:5 and Xmax Re Rec 0:5 . However it is to be noted that for the Case B1 (H 0:3D, B 0:5D) at the corresponding location the Amax is increased while the Xmax is decreased as shown in Fig. 11.
Fig. 9. Variation of Strouhal number with Reynolds number for a square cylinder.
110
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 10. Renormalised curves A=Amax vs X =Xmax at various Re.
Fig. 11. Maximum amplitude Amax and position Xmax vs Re
Rec =Rec for dierent shapes.
4. Conclusions This paper has presented some results of two-dimensional numerical simulations of the Benardvon K arm an hydrodynamic instability behind trapezoidal blu bodies using the spline method of
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
111
fractional steps. Numerical results con®rm the experimentally observed behaviour reported by Goujon-Durand et al. [1]. It is noted that for a ®xed trapezoidal height H the value of Strouhal number is almost independent of the value of the small base width B at lower Reynolds number. As Reynolds number Re increases, the curves describing the relationship between Re and St for dierent values of B gradually separate. Several previous computations by the authors [10,11,16] in such diverse ®elds as natural convection heat transfer, hydrodynamics of estuaries and groundwater ¯ow have indicated that the spline method is quite versatile, of high accuracy as well as being computationally ecient. The present results constitute the ®rst attempt at applying the scheme for the computation of the wake ¯ow behind blu bodies. A signi®cant advantage of the scheme is the capability of exploiting a variable mesh that may be placed to conform to an arbitrary body shape, thus avoiding coordinate transformations. Work is currently in progress to extend the study to three dimensions.
Acknowledgements Support for this project from the Natural Sciences and Engineering Research Council of Canada under grant number RGPIN8846 is gratefully acknowledged. We are grateful to the reviewers for bringing the recent reference of Dennis et al. [17] to our attention.
References [1] Goujon-Durand S, Rener K, Wesfreid JE. Downstream evolution of the Benard-von K arm an instability. Phys Rev E 1994;50:308. [2] Berger E, Wille R. Periodic ¯ow phenomena. Ann Rev Fluid Mech 1972;4:313. [3] Huerre P, Monkewitz PA. Local and global instabilities in spatially developing ¯ows. Ann Rev Fluid Mech 1990;22:473. [4] Zieilinska BJA, Wesfreid JE. On the spatial structure of global modes in wake ¯ow. Phys Fluids 1996;7:1418. [5] Rubin SG, Graves RA. Viscous ¯ow solutions with a cubic spline approximation. Comput Fluids 1975;3:1. [6] Rubin SG, Khosla PK. Higher order numerical solutions using cubic splines. AIAA J 1976;14:851. [7] Rubin SG, Khosla PK. Polynomial interpolation methods for viscous ¯ow calculations. J Comput Phys 1977;24:217. [8] Wang P. Spline method of fractional steps in numerical models of unsteady natural convection ¯ow at high Rayleigh number. Numer Heat Transfer 1987;11:95. [9] Okajima A. Strouhal numbers of rectangular cylinders. J Fluid Mech 1982;123:379. [10] Wang P, Kahawita R. Numerical integration of partial dierential equations using cubic splines. Int J Comput Math 1983;13:271. [11] Wang P, Kahawita R, Nguyen DL. Transient natural convection with density inversion from a horizontal cylinder. Phys Fluids A 1992;4:71. [12] Wang P, Kahawita R, Nguyen DL. Numerical simulation of buoyancy Marangoni convection in two superposed immicible liquid layers with a free surface. Int J Heat Mass Transfer 1994;37:1111. [13] Dusec J. Spatial structure of the Benard von Karm an instability. Eur J Mech B 1996;15:619. [14] Williamson CHK. Oblique and parallel models of vortex shedding in the wake of a circular cylinder at low Reynolds numbers. J Fluid Mech 1989;206:579. [15] Jackson CP. A ®nite-element study of the onset of vortex shedding in ¯ow past variously shaped bodies. J Fluid Mech 1987;182:23.
112
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
[16] Wang P, Kahawita R. Oscillatory behaviour in buoyant thermocapillary convection of ¯uid layers with a free surface. Int J Heat Mass Transfer 1998;41:399. [17] Dennis SCR, Nguyen P, Kocabiyik S. The ¯ow induced by a rotationally oscillating and transalating circular cylinder. J Fluid Mech 2000;407:23.
www.elsevier.com/locate/comp¯uid
Numerical simulation of the wake ¯ow behind trapezoidal blu bodies R. Kahawita *, P. Wang Department of Civil Engineering, Ecole Polytechnique de Montr eal, Montr eal, Qu ebec, Canada H3C 3A7 Received 22 November 1999; received in revised form 4 October 2000; accepted 10 November 2000
Abstract Two-dimensional numerical simulations of the Benard-von Karman hydrodynamic instability behind trapezoidal blu bodies has been studied using the spline method of fractional steps. For lower Reynolds numbers, about Re=Rec < 2 (Rec is critical Reynolds number), numerical results con®rm the experimentally observed behavior reported by Goujon-Durand et al. (Phys. Rev. E 1994; 50:308), i.e. the maximum amplitude of the velocity component oscillating with the fundamental frequency follows fairly well the scaling 0:5 0:5 law Amax Re Rec and the position Xmax Re Rec . The in¯uence of the trapezoidal shape on the value of the critical Reynolds number and on the vortex shedding is brie¯y discussed. It appears that the in¯uence of the trapezoidal height H is the dominant in¯uence on the value of Strouhal number when compared with the eect of the smaller trapezoidal base width B. Ó 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction The wake ¯ow behind an obstacle has been extensively investigated numerically and experimentally, most classical work having focused on circular cylinders. A variety of vortex shedding characteristics have already been discussed in the literature in the Reynolds number range 50±200 (see, for example, the review paper of Berger and Wille [2]). In recent years, the global nature of the wake ¯ow has attracted much attention [3]. In a recent paper, Dennis et al. [17] have conducted a comprehensive investigation of the temporal development of the viscous incompressible
*
Corresponding author. Tel.: +1-514-340-4711, ext.: 4506; fax: +1-514-340-2981. E-mail addresses: [email protected], [email protected] (R. Kahawita).
0045-7930/02/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 1 5 - 9
100
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
¯ow induced by a circular cylinder. The cylinder is assumed to be started impulsively while simultaneously performing time-dependent rotational oscillations about its axis. Solutions were obtained using a series expansion for small times and a spectral ®nite dierence technique for moderate times. Good comparisons with previous experimental results were obtained. For industrial ¯owmeters that exploit the relation between the shedding frequency and the ¯owrate or velocity, a trapezoidal shape is found to be more desirable than a cylinder, since a well-de®ned vortex emission due to clean separation at its sharp edges is assured. This is the subject of the present paper. Flow of a ¯uid past a trapezoidal shape is a more complicated phenomenon than that observed around a circular cylinder. Besides the Reynolds number, the vortex shedding depends strongly on the leading dimensions of the trapezoidal body. If observations are made in the far ®eld region, it would be dicult to determine whether the vortex shedding is due to a cylinder or a trapeze. In the near ®eld region however, the ¯ow structure and pressure ®eld is signi®cantly dierent. This is somewhat analogous to the ¯ow pro®les obtained in a jet from a square ori®ce. In the near ®eld region, information on the exit geometry is still retained by the ¯ow pattern. In the far ®eld however, the jet approaches circular symmetry (``memory is lost'') and the ¯ow closely resembles what would be obtained from a round ori®ce. Goujon-Durand et al. [1] recently reported the results of an experimental investigation into the vortex shedding from a trapezoidal body. They obtained the scaling laws for the evolution of the global mode describing the envelope of the peak to peak amplitude velocity oscillation in the wake ¯ow downstream of the body. They were also able to indicate the existence of a universal curve for the renormalized amplitude of the streamwise evolution of this amplitude. More recently, Zielinska and Wesfreid [4] presented some numerical simulations of two dimensional wake ¯ow behind an equilateral triangular obstacle using the spectral ®nite element commercial code NEKTON. Their numerical results con®rmed the experiments of Goujon-Durand et al. [1] by showing that the saturation amplitude of self sustained oscillations in the wake has a well de®ned maximum Amax and corresponding downstream location Xmax which are functions of the Reynolds number. Amax diminishes and moves further downstream as the critical point is approached. Since the pioneering work of Rubin and coworkers [5±7], the cubic spline integration technique and its applications in heat transfer and ¯uid ¯ow have developed rapidly. The present investigation is devoted to the numerical simulation of the two dimensional wake ¯ow behind a trapezoidal body using the spline method of fractional steps (SMFS) [8] which in discretized form is directly applicable on a non-uniform mesh in an x±y coordinate system. The advantages over using a standard ®nite dierence procedure is that no coordinate transformation is necessary, solutions are obtained directly in the physical plane. Furthermore, any arbitrary shape may be treated since computational points may be placed directly on the surface of the body. Since the trapezoidal shape has not yet been explored numerically, it is hoped that the results will provide insight into the experimentally observed behaviour reported by Goujon-Durand et al. [1]. The in¯uence of the trapezoidal shape on the value of the critical Reynolds number and on the vortex shedding is brie¯y discussed. The numerical results indicate that the shedding frequency initially increases linearly with Reynolds number (lower Reynolds numbers). The frequency then tends to level of until a maximum is reached. Surprisingly, after attaining this maximum, the shedding frequency decreases. This behaviour is consistent with Okajima's experimental results [9] for rectangular cylinders, however in this paper it is not discussed in detail. This study may thus be
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
101
considered as contributing to the present state of knowledge on this type of ¯ow as well as contributing to the state of the art in numerical algorithms for ¯ow simulation. 2. Governing equations The non-dimensional equations in stream function and vorticity form may be written: r2 W
X
1
oX oX oX 1 2 u v rX ot ox oy Re
2
with r2
o2 o2 ox2 oy 2
and u
oW ; oy
v
oW ox
3
where x is the longitudinal coordinate in the ¯ow direction and y is the transverse coordinate. Here W and X denote the dimensionless stream function and vorticity respectively while u and v are the respective dimensionless longitudinal and transverse velocities. Re U1 D=m, D is the trapezoidal width. 2.1. Boundary conditions On the trapezoidal surface: u v W 0; X
r2 W
4
At the in¯ow region x u 1;
v0
A !
1:
and X 0
5
On the out¯ow region x A ! 1: oW oX 0; and 0 6 ox ox The boundary condition on the vorticity in the out¯ow region corresponds to the free out¯ow of a ¯uid property, (here the vorticity). Extensive numerical experimentation in fact con®rmed that no discernible eect could be observed in the results as the out¯ow boundary was displaced further downstream. On the horizontal surfaces that de®ne the lateral limits of the computational domain, y b ! 1: u 1;
v0
and X 0
7
102
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
3. Results and discussion The SMFS schemes and the boundary conditions in discretized form may be obtained in direct fashion from the procedure detailed at length in earlier articles [8,10±12]. The development of the numerical algorithm will therefore not be discussed here. It is useful to indicate however that the spline discretized form is directly applicable to a non-uniform mesh in an x±y coordinate system, thus obviating the need for a computational coordinate transformation as is usually required with ®nite dierence schemes. The coordinate system and an illustrative mesh is shown in Fig. 1a and b respectively.
Fig. 1. (a) Coordinate system and (b) illustrative mesh used.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
103
In order to avoid computations with boundary conditions at in®nity, a speci®c ®nite domain was de®ned as follows: 5 6 x 6 30 and
10:5 6 y 6 10:5
i.e. the out¯ow domain was taken to be a distance of 30 units downstream from the rear base of the trapezoidal section. The time dependent non-linear coupled partial dierential equations were solved by considering a 121 121 grid. In order to accurately describe gradients that are expected to be steep in the boundary layer regions, a non-uniform grid in both the x and y directions was used. The essential advantage of the technique for the present problem lies in the fact that boundary conditions containing derivatives as shown in Eqs. (5) and (6) may be easily incorporated into the solution procedure since values of ®rst or second derivatives may be evaluated directly and maintain the same degree of accuracy when the algorithm that represents the spline approximation to the governing equations (1)±(3) is constructed. In order to verify the consistency of the present numerical study, the ®rst cases considered were those corresponding to existing results for wake ¯ows past a circular cylinder (see, for e.g., Ref. [13]) and behind an equilateral triangle [4]. In both cases very close agreement was obtained between the results of the authors cited above and the present computations. For example, for wake ¯ows past a circular cylinder, the value of Strouhal number at Re 49:4 is St 0:1241 (experiment of Williamson [14]) and St 0:1214 at Re 48 (numerical results of Dusek [13]), which compares favourably with St 0:1256 at Re 50 for the present results. The value of the critical Reynolds number Rec , for wake ¯ows past an equilateral triangle, is Rec 38:3 [4] and Rec 35 [15], while a value of Rec 37:5 was obtained in the present study. The accuracy and the convergence of the numerical solutions were veri®ed by mesh re®nement. All computations were performed on a PC. In the present numerical simulation, the blu body of trapezoidal section is characterized by its larger base width D the smaller base being perpendicular to the incident ¯ow. The following cases with dierent values of the smaller base B and the height H have been simulated: Case
A1
A2
A3
A4
B1
B2
B3
R1
R2
H B
0.7D 0.3D
0.7D 0.5D
0.7D 0.8D
0.7D 0.9D
0.3D 0.5D
0.5D 0.5D
0.9D 0.5D
0.7D 1.0D
1.0D 1.0D
Case R1 and R2 are particular cases, corresponding to rectangles with H 0:7 and H 1:0 respectively. In order to determine the critical Reynolds number, the computational strategy employed was to start the calculation procedure at a Reynolds number greater than critical and gradually diminish it towards the critical value. The principal problem encountered with this technique was that when the Reynolds number was very close to critical, convergence towards the stationary periodic state was very slow. The present Case A4 H 0:7D, B 0:9D approaches the experimental case (H 0:7D and the angle of the lateral faces being 4°) of Goujon-Durand et al. [1] in the laboratory. By extrapolation to zero of Amax , they obtained an estimate of the value of the critical Reynolds number
104
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Rec of about 58 (for the smallest blu body). The present numerical results indicate a Rec 42 with a corresponding Strouhal number St 0:11. For Re 65, i.e. Re 1:55Rec , the present numerical results indicate that the amplitude maximum occurs at a location Xmax of about 4:25D, Goujon-Durand et al. [1] obtained Xmax 4D or 4:5D. A typical instantaneous streakline view is shown in Fig. 2. The broken line in Fig. 2 corresponds to the envelope of transverse y values for which the amplitude along the ¯ow direction is a maximum. It is encouraging that this picture is very similar to the laser illuminated ¯uorescent view of vortex emission obtained by Goujon-Durand et al. [1]. Fig. 3 indicates the periodic variation of the transverse velocity v at dierent longitudinal positions x along the line of symmetry y 0 for Re 55 for Case A3. At a point very close to the smaller base, the amplitude A of oscillation is weak. As the value of x increases, the amplitude A increases until at a downstream distance of about x 5:5D it attains its maximum value; it then decreases as shown in the ®gure. It was noted that the value of the frequency at each x is constant. This agrees with the experimental results of Goujon-Durand et al. [1]. The constant frequency implies a synchronized oscillation and is one of the major features that characterizes the existence of a global instability in the full domain of shedding as observed by Goujon-Durand et al. [1]. Fig. 4 shows the periodic variation of the longitudinal velocity u at dierent y (transverse) stations along a vertical line x 5:4D for the same parameters as those of Fig. 3. u has a maximum in amplitude at about y 0:8D. As is well known for wakes behind blu bodies, u on the line of symmetry y 0 displays a frequency which is double that observed for o-symmetric lines. As the value of y increases from y 0 to about y 0:8D, the second harmonic transforms gradually into the fundamental while the amplitude of u increases to its maximum.
Fig. 2. Typical instantaneous streakline view.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
105
Fig. 3. Transverse velocity v at y 0 for dierent x positions (Case A3, Re 55).
Fig. 4. Longitudinal velocity u at x 5:1D for dierent y positions (Case A3, Re 55).
A time sequence of the streakline view for Case A3 with Re 55 is shown in Fig. 5. The von K arm an vortex street pattern and its transport by the main ¯ow is quite clear. The peak to peak amplitude of the oscillations of u, the longitudinal (x component) of velocity, was calculated from a time history of the oscillations at points in the wake along the x direction at dierent transverse stations (values of y). As an example, for Case A3 the amplitude of the
106
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 5. Time sequence of vortex street (Case A3, Re 55).
Fig. 6. Amplitude distribution of u in vertical plane at dierent x positions (Case A3, Re 55).
oscillations in u, for Re 55 from y 0 to 4D at dierent x positions are shown in Fig. 6. The following points may be noted: At each local ®xed plane of x, as y is increased i.e. moves away from the axis of symmetry, the amplitude of the oscillations increases steadily, reaches a maximum and then decreases. For example, at x 1:5D, the local maximum amplitude is about 0.061
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
107
Fig. 7. Envelopes of the peak to peak amplitudes of u for dierent Re (Case A3).
at a transverse position of y 0:7D. If the entire ¯ow ®eld is considered, the maximum amplitude Amax is about 0.278 and occurs at a position Xmax 5:5D and at a y value of 0:8D. This means that the maximum amplitude Amax occurs at about the same position X 5:5D for both longitudinal and transverse velocities. The value of Amax and the location at which it occurs is a function of the Reynolds number. Fig. 7 presents the distribution of the maximum amplitudes of the longitudinal velocity u as a function of Reynolds number throughout the ¯ow ®eld corresponding to Case A3. With decreasing Reynolds number, Amax diminishes and the position of Xmax moves farther downstream. An interesting point is that generally speaking, for a constant value of H, as the width B is increased, the critical Reynolds number also increases. For example, for a trapezoidal body of height H 0:7, when B 0:3D, Rec 34:5; B 0:5D, Rec 36; B 0:8D, Rec 39 while when B D, i.e. a rectangle Rec 40. For a constant value of B however, as the height H is increased, the critical Reynolds number decreases. For example, for a trapezoidal body with small base width B 0:5, when H 0:5D, Rec 35; H 0:3D, Rec 31. Fig. 8a and b present the relation between the Strouhal number, St, and the Reynolds number for a constant H with dierent B and for a constant B with dierent H respectively. As shown in Fig. 8a for a ®xed trapezoidal height H, at lower Reynolds numbers, the value of Strouhal number is almost independent of the small base width B. For example, when Re < 100, as the small base width B decreases from B 0:9D to 0:3D the values of St increase by 2%. As the value of Reynolds number increases, the curves de®ning the relationship between Re and St for dierent values of B gradually separate and then level-o to reach a ceiling quite rapidly. The experimental work of Okajima [9] however, has revealed that there is a fall-o from the maximum value as the Reynolds number is increased further. Clearly, it is quite possible that Fig. 8a would indicate this
108
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 8. (a) Relationship between the Strouhal number and Reynolds number for dierent values of B (H 0:7D) and (b) relationship between the Strouhal number and Reynolds number for dierent values of H (B 0:5D).
type of behaviour. In fact, it may be observed that Case A3 begins to demonstrate this fall-o. Another useful observation is that the larger the value of width B, the more rapid is the corresponding maximum attained and the lower the value of St as shown in the same ®gure.
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
109
This probably implies that since the separation points occur at the ends of the trapezoidal base with the larger width, the in¯uence of the smaller base width B on the wake ¯ow interaction appears to become signi®cant only at the lower Reynolds numbers. However, for a ®xed value of the base width B, the value of the Strouhal number evidently increases with increase in trapezoidal height H as shown in Fig. 8b. The trapezoidal height therefore plays an important role in the wake ¯ow interaction mechanism. It may also be noted that the frequency of shedding increases linearly with Reynolds number at the initial stage; this linear relation is prolonged to about Re=Rec 2. This behaviour was noted in the experiments of Goujon-Durand et al. [1]. Further increase in Reynolds number causes the frequency to increase at a slower rate until it reaches a relatively ¯at maximum. After this maximum is reached, a decline in frequency occurs as shown in Fig. 9. This interesting behaviour was observed in the experimental studies of Okajima [9]. The envelopes of the global modes corresponding to various Reynolds numbers are presented in Fig. 10. By renormalizing A and X respectively with Amax and Xmax , the points of A=Amax versus X =Xmax for dierent cases with dierent Reynolds numbers collapse on to a single universal curve con®rming that the properties of the global modes are universal. It is noted once more that these properties are almost independent of the width B for the cases H 0:7D at lower Reynolds number. However for the case B1 (H 0:3D, B 0:5D) the dierence is obvious. The amplitude at downstream locations is clearly enlarged. Two scaling laws for the maximum amplitude Amax and its location Xmax for dierent cases at various values of Reynolds numbers is presented in Fig. 11. These results are in good agreement with the scaling laws Amax Re Rec 0:5 and Xmax Re Rec 0:5 . However it is to be noted that for the Case B1 (H 0:3D, B 0:5D) at the corresponding location the Amax is increased while the Xmax is decreased as shown in Fig. 11.
Fig. 9. Variation of Strouhal number with Reynolds number for a square cylinder.
110
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
Fig. 10. Renormalised curves A=Amax vs X =Xmax at various Re.
Fig. 11. Maximum amplitude Amax and position Xmax vs Re
Rec =Rec for dierent shapes.
4. Conclusions This paper has presented some results of two-dimensional numerical simulations of the Benardvon K arm an hydrodynamic instability behind trapezoidal blu bodies using the spline method of
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
111
fractional steps. Numerical results con®rm the experimentally observed behaviour reported by Goujon-Durand et al. [1]. It is noted that for a ®xed trapezoidal height H the value of Strouhal number is almost independent of the value of the small base width B at lower Reynolds number. As Reynolds number Re increases, the curves describing the relationship between Re and St for dierent values of B gradually separate. Several previous computations by the authors [10,11,16] in such diverse ®elds as natural convection heat transfer, hydrodynamics of estuaries and groundwater ¯ow have indicated that the spline method is quite versatile, of high accuracy as well as being computationally ecient. The present results constitute the ®rst attempt at applying the scheme for the computation of the wake ¯ow behind blu bodies. A signi®cant advantage of the scheme is the capability of exploiting a variable mesh that may be placed to conform to an arbitrary body shape, thus avoiding coordinate transformations. Work is currently in progress to extend the study to three dimensions.
Acknowledgements Support for this project from the Natural Sciences and Engineering Research Council of Canada under grant number RGPIN8846 is gratefully acknowledged. We are grateful to the reviewers for bringing the recent reference of Dennis et al. [17] to our attention.
References [1] Goujon-Durand S, Rener K, Wesfreid JE. Downstream evolution of the Benard-von K arm an instability. Phys Rev E 1994;50:308. [2] Berger E, Wille R. Periodic ¯ow phenomena. Ann Rev Fluid Mech 1972;4:313. [3] Huerre P, Monkewitz PA. Local and global instabilities in spatially developing ¯ows. Ann Rev Fluid Mech 1990;22:473. [4] Zieilinska BJA, Wesfreid JE. On the spatial structure of global modes in wake ¯ow. Phys Fluids 1996;7:1418. [5] Rubin SG, Graves RA. Viscous ¯ow solutions with a cubic spline approximation. Comput Fluids 1975;3:1. [6] Rubin SG, Khosla PK. Higher order numerical solutions using cubic splines. AIAA J 1976;14:851. [7] Rubin SG, Khosla PK. Polynomial interpolation methods for viscous ¯ow calculations. J Comput Phys 1977;24:217. [8] Wang P. Spline method of fractional steps in numerical models of unsteady natural convection ¯ow at high Rayleigh number. Numer Heat Transfer 1987;11:95. [9] Okajima A. Strouhal numbers of rectangular cylinders. J Fluid Mech 1982;123:379. [10] Wang P, Kahawita R. Numerical integration of partial dierential equations using cubic splines. Int J Comput Math 1983;13:271. [11] Wang P, Kahawita R, Nguyen DL. Transient natural convection with density inversion from a horizontal cylinder. Phys Fluids A 1992;4:71. [12] Wang P, Kahawita R, Nguyen DL. Numerical simulation of buoyancy Marangoni convection in two superposed immicible liquid layers with a free surface. Int J Heat Mass Transfer 1994;37:1111. [13] Dusec J. Spatial structure of the Benard von Karm an instability. Eur J Mech B 1996;15:619. [14] Williamson CHK. Oblique and parallel models of vortex shedding in the wake of a circular cylinder at low Reynolds numbers. J Fluid Mech 1989;206:579. [15] Jackson CP. A ®nite-element study of the onset of vortex shedding in ¯ow past variously shaped bodies. J Fluid Mech 1987;182:23.
112
R. Kahawita, P. Wang / Computers & Fluids 31 (2002) 99±112
[16] Wang P, Kahawita R. Oscillatory behaviour in buoyant thermocapillary convection of ¯uid layers with a free surface. Int J Heat Mass Transfer 1998;41:399. [17] Dennis SCR, Nguyen P, Kocabiyik S. The ¯ow induced by a rotationally oscillating and transalating circular cylinder. J Fluid Mech 2000;407:23.