Numerical Solution of the Problem of Vortex Street Development

A method is described for the solution of time-dependent problems concerning the flow of viscous incompressible fluids in several space dimensions. The method is numerical, using a high-speed computer for the solution of a finite-difference approximation to the partial differential equations of motion. The application described here is to a study of the development of a vortex street behind a plate which has impulsively accelerated to constant speed in a channel of finite width; the Reynolds-number range investigated was 15 ≤ R ≤ 6000. Particular attention was given to those features for which comparison could be made with experiments, namely, critical Reynolds number for vortex shedding, drag coefficient, Strouhal number, vortex configuration, and channel-wall effects. The nature of the early stages of flow-pattern development was also investigated.

INTRODUCTION

UNDER certain conditions, the motion of an object through a fluid results in generation of a wake containing a double row of vortices, called a vortex street. The behavior of the flow, which has received abundant experimental and theoretical study, appears to depend upon relatively few parameters, particularly when the flow is two dimensional (i.e., is perpendicular to the axis of a very long cylindrical body). Principal among the parameters is the Reynolds number, R ≡= u₀ d/ν, where u₀ is the free-stream fluid speed, d is the height of the body (diameter, if a circular cylinder), and ν is the kinematic viscosity. For R ≤ 40, the flow is steady; for larger values of R the steady flow is unstable and vortices are shed. As R increases still higher, the vortices become less and less regular; in other words, for successively higher Reynolds numbers, the effects of turbulence become more and more pronounced.

Numerous experiments have demonstrated these and a number of related effects; in addition several authors have described theoretical studies designed to shed light on the process. References to some of these have been made in the present paper, and a nearly complete history of this extensively studied effect can be found through the associated bibliographies. A survey of the literature shows, however, that while almost every aspect has been investigated experimentally, no theoretical solution has yet been described which demonstrates in detail the complete development of a vortex street. The reason for this is easily understood; the time-wise evolution of the flow field is observed to be considerably complicated, so that the functions describing the processes must likewise be complicated, as well as difficult to extract from the basic equations. Our approach to the theoretical study has therefore been somewhat different from that ordinarily undertaken, and employs a numerical method for solving the equations. The technique is possible only through the use of a high-speed electronic computer, and all results reported were obtained with the IBM Electronic Data Processing Machine, type 7090. In addition, all flow configurations were processed directly from the computer by SC-4020 Microfilm Recorder, and are reproduced in the figures without retouching or alteration.

The value of such a theoretical study is severalfold. First, it demonstrates that the process of vortex street development depends only upon the coarse features of the flow; this is so because the method as used does not resolve many of the fine details. Second, it allows a close examination of the process through all steps of the development, thereby giving considerable insight into the manner in which the instability first starts to grow, as well as into the details of subsequent periodic vortex formation. Third, the success of the method in this particular application shows that its extensions can be made to a number of related problems with likelihood of equally valid results.

This is by no means the first application of numerical methods to incompressible-flow studies. The pioneering work of Thom [1] [2] and associates demonstrated considerable success for time-independent problems, while more recently a group at Los Alamos [3] has employed a method suggested by von Neumann for computing the time-dependent flow of an inviscid fluid. Payne [4] has done time-dependent problems in which the effects of viscosity have been included, and Abernathy and Kronauer [5] have employed a special model for the investigation of vortex street development. The technique described here, however, differs from the others in a number of respects, employing useful features from several and adding many that are new. The resulting computer method appears thereby to be both broadened in scope of applicability and increased in computer efficiency.

DESCRIPTION OF THE METHOD

in which u and v are the components of velocity along the x and y directions, respectively; p is pressure, ρ is the (constant) density; and ν is the (constant) kinematic viscosity coefficient. In addition to these variables, we introduce for convenience the usual stream and vorticity functions, ψ and ω respectively, such that

An approximate solution is accomplished through the introduction of a finite-difference mesh of cells of sides δx and δy, dividing the spatial region of interest in the manner shown in Fig 1. In this way, then, the entire continuous flow field becomes described by a finite number of quantities whose variations with time represent the changes of configuration. In addition, even the time variations proceed through finite time intervals δt.

The numerous details involved in properly solving the equations in a manner which is both stable and sufficiently accurate, have been described completely in a report by Fromm [6], so that only a brief description is appropriate here. The basic steps involved in advancing the configuration through the cycle from time t to time t + δt are as follows:

Through these steps, care must be taken to use proper finite-difference analogies for the desired boundary conditions. The manner in which the noslip conditions at walls have been incorporated is thoroughly described in reference 6; the two ends of the channel (the inflow and outflow surfaces) were handled through the convenient assumption of periodicity, so that the computational region was one period of an infinite channel with, therefore, an infinite sequence of similar objects, one in each period. The physical consequences of this periodicity assumption are discussed in the next section; the mathematical consequences are twofold. First, the end boundary conditions are easily and accurately handled; second, the drag force is simply calculated by balancing momentum changes and wall stress against object back-force.

It is apparent that the degree of accuracy of the finite-difference approximations depends critically upon the parameters δx, δy, and δt. Only a preliminary analysis has been accomplished for predicting error size as a function of these parameters, so that proof of accuracy has mainly depended upon two types of comparisons. In the first type, results from the computer have been compared with analytical solutions for several simple situations. Thus, for example, the Rayleigh problem has been studied, in which vorticity diffuses from a wall moving in its own plane; the resulting accuracy restrictions on the parameters were then applied to the choice of their values for more complicated calculations. The second type of comparison is with experiments themselves, and form the really crucial test of the computing method. Part of the vortex-street discussion is concerned with these comparisons.

Also dependent upon the artificial parameters is the stability of the finite-difference equations. A linear stability analysis has been accomplished which shows the most important of the restrictions. It will be seen, however, that at high Reynolds numbers (R > 1000), there is apparently an instability whose boundedness suggests association with the nonlinearity of the equations, but for which no analysis has been made.

DEVELOPMENT OF THE FLOW PATTERN

In all the calculations to be described, the object was a rectangular cylinder, or plate, centered between two parallel walls whose spacing, H, in all cases was 3.0 units. Two plate heights were studied, d = 0.50 units and d=0.25 units; in both cases the plate thickness was 0.125 units, while in all calculations, δx = δy = 0.125. The walls were forced to move at constant velocity u₀ = 4.0 distance units per unit time (thereby establishing the time scales used in all the discussions). All examples started, at time t = 0, with the walls and fluid impulsively accelerated to this velocity, and the first-cycle iteration procedure immediately adjusted the configuration to the nonviscous laminar-flow solution. Advancement of the configuration through subsequent time cycles resulted in a gradual transition to the viscous steady-state solution whose most prominent feature is an eddy pair just behind the plate. Since the solution procedure preserves symmetry to approximately one part in 10⁵, the steady solution persists for long times, even for large values of R. Thus we found it desirable to introduce a perturbation, accomplished by artificially increasing the value of ω by a small amount at three mesh points just in front of the plate; this was done at a time when the double eddy pattem was fairly well established. In all cases, the perturbation was small enough that no immediate change was visible in the flow pattern; nevertheless, such a small perturbation was always effective in starting the vortex shedding process within a fairly short time, provided, of course, that R was sufficiently large. For R ≤ 40 we found that the steady state flow pattern never visibly changed after introduction of the perturbation.

A typical sequence showing development of the vortex street is given in Fig 3 for R = 200, d/H = 1/6. These streamline plots show the initial nonviscous laminar flow pattem, development of the symmetric eddy pair, the subsequent commencement of vortex shedding, and the establishment of the late stage vortex street. (There are actually many time cycles between the plots shown here; the time increment per cycle was δt = 3.52 × 10^-3.) An examination of these and numerous intermediate frames leads to the following description of the onset of shedding. (Many more frames, together with a more detailed analysis, are given in the report by Fromm [6]). The initial inviscid solution (first frame of Fig 3) rapidly changes so that by a time t = 0.21, the steady state eddy region is about half its final size. At this time the flow is perturbed by the introduction of a small positive excess of vorticity just in front of the plate, but even by the time t = 0.70 (second frame) no significant effect can be detected in the streamline pattern. Shortly thereafter, however, there is a relative increase in the size of the lower eddy over that of the upper one; it becomes longer and moves slightly toward the centerline while the upper one becomes shorter and moves slightly upward away from the centerline and into the lateral stream. A quick consequence of this is a relative strengthening of the upper eddy, the upper part of which is accelerated by its contact with the stream. In addition, the upper eddy is pulled somewhat downstream and inward toward the centerline. This stage is shown in the third frame of Fig 3. The strengthening of the upper eddy, however, results in subsequent acceleration of the lower one; the upper one weakens and cannot be replenished because of its motion towards the centerline. Intermediate frames show a process which could be called streamline capture by the lower eddy from the upper one. In this way, then, the upper eddy actually disappears; it has been completely absorbed, rather than shed. At the same time, the upstream pocket left by the movement of the upper eddy is subjected to rotational acceleration by both the primary (now strengthened) lower eddy and by the lateral stream, so that a second upper eddy is formed. This is the stage shown in the fourth frame of Fig 3. Now the lower eddy is strong enough that the secondary upper eddy cannot completely absorb it. It is, however, far enough downstream that it contacts the upper part of the lateral stream, producing a high local concentration of vorticity which rapidly diffuses away. Meanwhile, the secondary upper eddy has repeated the process which launched the primary lower one into downstream motion, and by the time of the fifth frame of Fig 3, there is forming a secondary lower eddy. In this fashion the process continues, and a sample at a later time is shown in the last frame.

Fig 3 - Sequence of streamline configurations for R = 200, d/H = 1/6 at time t = 0 (upper left), t = 0.70 (upper right), t = 1.58 (middle left), t = 2.29 (middle right), t = 2.64 (lower left), and t = 3.16 (lower right).

Note that in the last frame of Fig 3 there are two manifestations of the periodicity condition imposed on the calculation. At the far left of the frame the streamlines are now slanted, showing that the incoming flow is being disturbed by the vortices shed by the upstream plate in the period to the left. (Actually, the disturbance comes from the fact that the fluid which flows out the right side of the computation region flows back in again from the left.) The second manifestation of periodicity is seen in the weakness of the newly forming vortex in the last frame; the over-all general slowing down of the flow has by this frame somewhat reduced the Reynolds number from the initial value of 200.

With this picture of the processes involved in the onset of vortex shedding, it is possible to see qualitatively why the flow is stable for low Reynolds number. The principal consideration seems to involve the relative contact that an eddy has with the lateral stream and with its adjacent partner. At low Reynolds number, the eddy region is small and the contact line between eddies is relatively long. Any perturbation which moves one of the pair into the lateral stream will still result in acceleration, but the sharing will be more immediate so that neither gains enough advantage over the other. Indeed, the perturbation will damp as a result of viscous drag of the oscillating motion on the rear of the obstacle. Thom's findings [1] tend to support this interpretation; he has shown experimentally for a circular cylinder that when a certain measure of the eddy region length is less than a certain measure of the wake thickness, then the flow is stable.

Figure 4 shows three typical filament-line plots for the same calculation as in Fig 3, obtained from the computer by the process of calculating through each time step the sequence of positions of lines of massless test particles whose velocities were given by interpolation among the four nearest velocity-point values.

Little quantitative experimental information has been found for comparison with the time dependent early stage developments. An exception is given by the results of Schwabe [7], quoted by Goldstein [8], who obtained as a function of time, the drag on a circular cylinder at R = 580. His results are shown in Fig 5, together with the drag coefficient histories from a number of our calculations at various Reynolds numbers and for d/H =1/6. (In our calculations, the reference velocity for drag coefficient was taken to be the mean fluid speed.) Since the results of Schwabe refer to a circular cylinder, it is not possible to make a direct comparison; in particular, the correspondence between his timing and ours is obscure. Schwabe's timing is based on a length scale given by the cylinder radius; if we correspondingly base it on plate half-height, then we get curve B of Fig 5, while basing it on plate half-thickness gives curve A. The latter curve agrees well with the calculational results, and an additional argument can be given for basing the time on plate half-thickness: when calculations were performed with d/H = 1/12, we found that the time of maximum drag coefficient was essentially unchanged from that obtained with d/H = 1/6.

Fig 5 - Time variations of drag coefficient for the calculations with d/H = 1/6. Schwabe curve A has timing based on plate thickness, while curve B is based on plate height. Numbers are initial values of R.

FUNCTIONALS OF THE MOTION

The drag coefficient histories shown in Fig 5 were calculated with a reference velocity which was the mean velocity of the fluid. For comparison with experiments, we have instead used as a reference velocity a value which was scaled down from an initial 4.0 by an amount proportional to that by which the velocity in front of the plate had decreased from the initial inviscid value. The results for d/H = 1/6 are included in Fig 6, which also contains the mean experimental results for a flat plate given by Flachsbart [9] and for a circular cylinder given by Thom [1]. It was noticed that this choice of reference velocity gave drag coefficient values which were consistently equal to those at a time t = 0.25 in the histories shown in Fig 5. For the larger Reynolds numbers, for which fluctuations did not allow for this choice of reference velocity, we therefore also used the values at that same time. The results give nearly constant drag coefficient for R > 500, up to R = 6000 at which C_D = 1.87.

Fig 6 - Drag coefficients for various Reynolds numbers as calculated, and as obtained experimentally for circular cylinder and plate. Datum points are from various times in various calculations; dashed line is eye drawn mean.

Also investigated was the effect of channel width on drag coefficient, at R = 300. Two channel widths were studied, d/H = 1/6 and d/H = 1/12. The reference velocity for drag coefficient was taken as the mean fluid speed. Again experimental results are available only for much larger Reynold's number, so that only the rate of variation with channel width could be compared. We thus found a common factor and scaled the two results by that same amount, as shown in Fig 7. That figure also includes the experimental results of Flachsbart [9] and of Fage and Johansen [10]. The line is a plot of C_D = 1.89 + 3(d/H). Glauert [11] has predicted the variation of drag coefficient with d/H and has compared his results with those of Fage and Johansen. Neglecting the experimental point for largest d/H as being unreliable, he obtained good agreement with his prediction, C_D = 1.818 + 4(d/H). The inclusion of more data together with the calculated results suggests, however, that his prediction, which involved numerous assumptions, gives somewhat too rapid a variation.

The Strouhal number, defined as S ≡ f d/u₀ in which f is the shedding frequency on one side, was obtained from each of the calculations. For times early enough that the reference velocity had not changed much, we obtained S = 0.119 for R = 100, increasing to S = 0.137 at R = 300. For comparison, Fage and Johansen [10] give the value S = 0.146 for R = 1.5 × 10⁵. All of these values lie slightly below those for a circular cylinder, but show the same tendency to gradually increase with R.

Also observed in the calculations was the ratio of vortex speed to wall speed, which was insensitive to Reynolds number. The average value is 0.79, which compares favorably with the theoretical value given by Heisenberg [12], 0.771.

Likewise agreeing well with experiments are velocity profiles across the wake. The calculated results are compared in Fig 8 with those obtained experimentally by Kovasznay [13] for a Reynolds number of 50 (computed) and 56 (experimental). In the figure, x/d is the number of plate heights downstream from its center, while y/d is the number of plate heights above the center line.

Fig 8 - Velocity as a function of distance from centerline for various distances behind plate. Solid lines are calculated; datum points are from experiments by Kovasznay. Zero line is shown for each curve.

Finally, the flow configurations themselves have been extensively studied. Examples have already been presented in Figs 2-4, and several additional ones are given in Figs 9-11. The filament line configurations of Figs 2 and 10 agree well with those obtained by Thom [reference 1, PI. 12, parts (c) and (f)] for nearly the same Reynolds numbers and value of d/H. Finally, the filament line configuration shown in Fig 11 was obtained from a calculation at R = 6000. While the results at this high Reynolds number showed some instability, manifested by the irregular filament lines at the left, nevertheless it has been possible to compare the configuration with several previously unpublished experimental ones for this same R and d/H, furnished us by Dr A M Lippisch, and it is seen that the agreement is good.

CONCLUSIONS AND EXTENSIONS

While the results presented here demonstrate that the method used is applicable to a variety of problems in the flow of viscous incompressible fluids, it is apparent that several directions of generalization are still open for investigation. Especially of interest is the question of how to handle free surfaces and discontinuities between fluids; some ideas on this subject have already been proposed by von Neumann and applied by his associates [3], and several other approaches are also being investigated. Also useful would be a modification allowing for the diffusion of heat from boundaries and through the fluid. This, too, is being investigated, and little difficulty is anticipated in developing the technique. It now seems likely that with these and other appropriate generalizations, the technique will be useful for a number of design and analysis problems concerned with incompressible fluid flow.

ACKNOWLEDGMENTS

We are extremely grateful to Dr A M Lippisch of the Collins Radio Company who performed numerous experiments on our behalf and permitted our use of his unpublished results. We also wish to thank Professor A Thom and Professor W A Mair for furnishing unpublished filament-line configurations which have been very useful for comparisons; several of these are included in Reference 6.

Our work was performed under the auspices of the United States Atomic Energy Commission.

REFERENCES

[2] A. Thom, C J Apelt, Field Computations in Engineering and Physics (Van Nostrand, London, 1961).

[3] A Blair, N Metropolis, A Taub, M Tsingou, A Study of a Numerical Solution to a Two-Dimensional Hydrodynamical Problem, Los Alamos Scientific Laboratory Report Number LA-1265 (1958).

[8] Modern Developments in Fluid Mechanics, edited by S Goldstein (OUP, 1938).