Computers FluidsVol. 22, No. 2/3, pp. 229237, 1993 Printed in Great Britain. All rights reserved
COMPUTATION
00457930/93 $6.00+0.00 Copyright © 1993 Pergamon Press Ltd
OF UNSTEADY
VORTICAL
FLOWS
M . BREUER, l D . H~NEL,2t J. KLOKER 3 a n d M . MEINKE 3 tlnstitut f~ir Hydromeehanik, Universitiit Karlsruhe, Karlsruhe, Germany 2Institut fiir Verbrennung und Gasdynamik, Universit/it Duisburg, Duisburg, Germany 3Aerodynamisches Institut, RWTH Aachen, Germany
(Received 1 April 1992; in revisedform 16 September 1992) Al~ractThe NavierStokes equations are solved for unsteady, vortical flow problems. It is shown that the solution is sensitive to the formulation of the numerical algorithm, due to the appearance of small characteristic scales in time and space. Several aspects of possible influences will be discussed and demonstrated by computational results. Solutions of the flow problems shown are taken from current studies of vortex breakdown, the development of vortex streets and the simulation of jets.
INTRODUCTION
The major part of natural or technical flows is viscous and timedependent. An unsteady flow can either be generated by external excitation or develops because of a hydrodynamic instability. Examples for the latter are the formation of a Karman vortex street, the vortex breakdown or the transition to turbulent flow. The numerical simulation of such flows requires a method of solution for the NavierStokes equations, making high demands on the spatial as well as the temporal accuracy. Flows in the vicinity of instabilities are very sensitive to small disturbances. They usually exhibit different length scales in time and space, so that care and a large computational effort is necessary for the solution of such problems. As the domain of integration normally is of limited extent, artificial inflow and outflow boundary conditions have to be specified, which can have substantial influence on the solution in the interior. The spatial numerical dissipation terms, present in all solution schemes, can also have large effects on the temporal development of a flow, because the spatial and temporal accuracy of a numerical algorithm are directly coupled. The errors induced by the numerical scheme are difficult to estimate, since they are always a combination of both the spatial and temporal errors resulting from the discretization, as well as those from the boundary conditions. Due to the different length scales in time and space, an investigation of scalar model equations is not sufficient to analyze the properties of a solution algorithm for the NavierStokes equations. Numerical experiments are therefore an essential tool to clarify the different behavior of numerical schemes. In the present paper a few examples are presented which reveal numerical effects in the simulation of unsteady vortical flows. There are details of current studies in progress, which have been selected to demonstrate the diversity of problems in such flows. VORTEX
FORMATION
BEHIND
A CYLINDER
Computations of separating flow behind an airfoil at low Reynolds number and a high angle of attack have shown an influence of the numerical scheme on the instantaneous forces [1]. To study this phenomenon in more detail, the geometrically simpler, but physically similar complex flow behind a circular cylinder has been chosen for the numerical investigations. Since this fundamental vortical flow has been the object of numerous experimental and numerical studies, e.g. in Ref. [2] comparisons with experimental and other numerical data can be drawn. The flow around a cylinder exhibits a number of different flow patterns. The easiest, unsteady pattern appears immediately after an impulsive start of a cylinder. A symmetric pair of vortices is generated behind the cylinder with a characteristic extension growing up to the order of the diameter of the cylinder and with time scale tuo~/D = O (1). This flow can be considered as a basic test for the temporal accuracy of a scheme. ?Author for correspondence. 229
M. BREUER et al.
230
4
1.0
~ ' 1
0.83
0.6. _.o)J 0.4.
2
0.20.0 f oo
o;5
lo
i'.5
f.o
i.~
io
tuo= /D Fig. 1. Time history of the length o f the closed wake behind an impulsively started cylinder (Re = 3111111, Ma = 0.3, 113 x 81 grid points). N, Numerical result [4]; V , experiment [3]. NavierStokes solution with:
1, van Leer fluxvector splitting; 2, van Leer, Hgnei fluxvector splitting; 3, ~ t r d differelr~ ~heme, nodecentered; 4, cellvertex scheme; 5, central differencescheme, nodecentered solution with 225 x 225 grid points. The time history of the length of the closed wake is shown in Fig. 1 for an impulsively started cylinder. The solution with a finer grid (225 x 225 points) is in good agreement with the experimental data [3] and a numerical result [4]. Solutions on a coarser grid (113 x 81 points) show deviations for larger times, but there are no significant differences among the different methods used. For a cylinder in a steady uniform flow and with a Reynolds number larger than about 40, the flow becomes periodic due to the formation of a vortex street. The characteristics of this unsteady shear flow separating from the cylinder can be interpreted as an instability of the near wake [2]. The vortex formation starts from small vortices near the wall, which grow further downstream into the largescale major vortices. The solution for this flow reacts very sensibly to the numerical dissipation, which has a filterlike action on the small scale vortices and therefore influences further temporal development. The numerical dissipation, which is part of each solver, [e.g. 5], has therefore to be analyzed and minimized. Corresponding investigations were carried out in Refs [6, 7], using an explicit RungeKutta scheme with different spatial discretizations (nodecentered central, cellvertex central, fluxvector splitting). Figures 2(af) show the lift coefficient vs time and the instantaneous isobars of the flow around a cylinder with Ma~ = 0.3 and Redo = 3000. The solution in Figs 2(a, b), computed with a central scheme and a small value of fourthorder damping, agrees very well with results from the literature. The other two examples demonstrate the influence of the numerical dissipation, which damps the highest frequencies (smallest wavelengths) and therefore falsifies the solution. Typical damping mechanisms are the highfrequency artificial damping terms in central schemes. In Figs 2(c, d) the damping term used is 4 times higher than in Figs 2(a, b). Similar effects have approximate Riemann solvers, used for the discretization of the convective terms, if the linear eigenvalues are not represented correctly. An example is given in Figs 2(e, f) for fluxvector splitting.
(Fig. 2 Opposite) Fig. 2(a,b). Periodical flow around a circular cylinder (Re = 3000, Ma  0.3, 177 × 113 grid points); central cell vertex scheme, ~(o = 1/256: (a) lift coefficient vs time, (b) fines of constant pressure. Fig. 2(c,d). Periodical flow around a circular cylinder (Re = 3000, Ma = 0.3, 177 x i 13 grid points); central nodecentered scheme, E(4) = 1/64: (c) lift coefficient vs time, (d) lines o f constant pressure. Fig. 2(e,f). Periodical flow around a circular cylinder (Re = 3000, Ma = 0.3, 177 × 113 grid points); fluxvector splitting (van Leer): (e) lift coefficient vs time, (f) lines o f constant pressure.
8
0
b
~.01.51.0o.5
Cl 0.0
0.5
1.0
1.5 2.o
o
0
0
0
o~ o~ .m
0
o~ o~
~0
1.0 I
1.5 I
I
0.5
Cl I
0.0 I
0.5 I
1.0 I
1.5
2.0
0
f,
O
8...0
0)o
t~
b
rD~ "
2.0 I
1.5
1.0
0.5
cl 0.0
0.5
1.0
I
1,5
2.0
o
0
232
M. BREUER et al. PLANE
JETS
Distinct largescale vortex flows can be observed in free shear layers and jets. A typical application for such flows is the compressible flow in the cylinder of a piston engine. The vorticity is generated mainly by the injection during the intake stroke and it is altered considerably during the succeeding compression and expansion strokes. The physical aspects of such a flow, but also the numerical influence on the solution, are the subject of current investigation [8, 9]. In the first studies, the accuracy of different methods of solution were tested for some typical flow problems (results were presented, for example, in Refs [5, 10]). The numerical scheme, selected finally, is an explicit, secondorder Godunovtype method, based on the PLM method of Colella and Glaz [11]. This method was extended to the complete NavierStokes equations, and reformulated recently for 3D curvilinear coordinates [12]. One of the most important problems for the investigation of largescale vortex flow is the resolution of the smaller scales. An estimation of the minimum length scale ~'min in 2D and 3D flows can be made with the assumption by Kreiss et al. [13] for incompressible flow. The estimation, using typical quantities of the injected flow in piston engines, shows that the smallest scales cannot
100 x 100
200 x 200
t=3
t=5 Fig. 3. Influence of the grid size on the density contours at different times. Plane jet with ~ = 45 °, Mart = 0.5, Re D = 106; Godunovtype P L M method.
Computation of unsteady vortical flows
233
be resolved sufficiently with the available capacity of computers. Nevertheless, the major (larger) vortex structures are of great interest for the physical and technical understanding of this flow. Therefore, variations of the grid size were carried out to study the numerical influence due to the unresolved scales. Figure 3 shows an example for an oblique jet entering a quadratic chamber with three rigid adiabatic walls and with one outflow boundary (right side). The lines of constant density in Fig. 3 (left) with 100 x 100 and in Fig. 3 (right) with 200 x 200 grid points show the instantaneous vortex structure of the incoming jet. Both solutions show similar major flow structures, but much more detail is present in the case of the finer resolution. In particular, the number of vortices is larger and the gradients become steeper. The discrepancy between coarse and fine grid solutions grows with time, such that at later times the solution on the coarser grid becomes unrealistic. In a similar manner to the grid resolution, the numerical approximation of the boundary conditions may influence the solution. The most critical boundary formulations are those for the inflow or outflow boundaries. The simplest formulation for a free stream boundary, which results in stable calculations, is the extrapolation of zeroth or firstorder of the variables on the boundary. In many cases, however, vortices are approaching or even crossing such an outflow boundary. Then a simple boundary condition like this partially reflects the flow. This effect is demonstrated in Fig. 4 (top) for the jet flow through an open boundary to the right. In the literature a number of socalled nonreflecting boundary conditions are known, which, however, are not really nonreflecting, especially not for unsteady flows. In most cases the method of trial and error, and physical "understanding" is necessary to determine a suitable approximation. In the present study, a rather simple "nonreflecting" boundary condition has been found to be useful for unsteady flows. This condition is based on the relaxation of the pressure at the boundary against a prescribed exterior pressure p~, far away from the boundary [14]. This condition readspb =Pb1 "~( P o o   Pb1)"Ax/Ld, where Pb means the pressure at the boundary, Pb 1 is the pressure at the adjacent point and L d is a characteristic "decay" length, chosen here as the distance between the inlet and outflow boundary.
Fig. 4. Influence of the outflow boundary condition on the entropy contours. Plane jet with ~ = 0°, Ma~t =0.5, ReD= 10s; Godunovtypc PLM method, firstorder extrapolation (top), nonreflecting boundary condition (bottom).
234
M. B~tw~ et al.
Compared with Fig. 4 (top), this nonreflecting boundary condition, applied to the same flow, yields a nearly undisturbed jet flow through the boundary, as shown in Fig. 4 (bottom). VORTEX BREAKDOWN IN INCOMPRESSIBLE FLOW Another type of unsteady, vortical flow, studied currently in Refs [1517], is the 3D unsteady flow of a bursting vortex, called vortex breakdown. Vortex breakdown is characterized by an abrupt change in the structure of the core of a swirling flow. Its occurrence is often marked by the presence of a free stagnation point on the axis of an axial vortex flow combined with flow reversal. Basically, three types of vortex breakdown are observed, mild (double helix) breakdown, spiraltype breakdown and nearly axisymmetric bubbletype breakdown. The bursting of an isolated vortex was simulated numerically in a 3D box using a Cartesian mesh. The method of solution is a dualtimestepping scheme with artificial compressibility in the pseudo time level [16]. This method for the unsteady, 3D NavierStokes equations of an incompressible flow is secondorder accurate in time and space and has been tested for a number of different flow problems. In order to concentrate the study on pure vortical flow, and to reduce the computational work, an isolated vortex embedded in an axial flow is used as the initial condition. Major problems arise, similar to those in other unsteady, vortical flow simulations, in the formulation of the boundary conditions. The distribution of axial and circumferential velocity is prescribed in the inlet plane, the other quantities are extrapolated. The outflow plane requires a nonreflecting boundary condition to simulate undisturbed vortical flow through the boundary. In the present case a new extrapolation method was developed, which is based on the analytical continuation of the flow equations, assuming timedependent, but parabolic flow (in space) at the boundary. This boundary condition works sufficiently well. The lateral boundaries, surrounding the vortex, are free stream boundaries. They simulate, for example, the far effects of a nozzle contour or a delta wing. They can be used to control the axial and temporal development of the vortex breakdowns, e.g. by prescribing appropriate pressure gradients at the lateral boundaries. An adverse pressure gradient is usually necessary to initialize the breakdown process. With such a boundary condition the bubbletype breakdown could be provoked. But due to the selfinduced motion of the breakdown bubble, consisting of one or more ring vortices, the burst vortex moves towards the inflow boundary, which makes it impossible to carry on the numerical simulation. A new timedependent condition for the lateral boundaries was developed [18], to keep the burst vortex in the domain of integration. This condition is based on the idea that the pressure at the lateral boundaries should enforce a solution of the steadystate equations in a part near the inflow boundary, with the aim to suppress the selfinduced upstream movement. Thus, the pressure at the boundary (z, r = R) is written as pn+ l(z '
R) =p"(z, R) + or. Ap(z).
The pressure correction
\Otjz.~o.tdr+f:(~t),=o dz" is the difference between the pressure as it would result from the steady momentum equations and the pressure from the actual timedependent calculation, The weighting factor a = 0 means the fixed, prescribed pressure at the lateral boundary, whereas ~ = 1 forces the flow to become stationary. In the present calculations a small value of = ~ 0.1 was chosen to obtain a smooth change from the steady flow near the inflow boundary to the unsteady flow further downstream. With this boundary condition it was possible to avoid the selfinduced movement of the burst vortex, and thus the computations could be continued for longer times. At later times the flow mode changed from the axisymmetric breakdown to the asymmetric spiral mode. A typical sequence of streaklines are shown in F i g . 5 f o r R e = 2 0 0 a n d a n initial (dimensionless) vortex strength ~ff= 0.68. Starting from the initial straight vortex (time T = 0), the vortex swells (T = 84) and forms an axisymmetric bubblelike breakdown configuration. Without adapting the lateral boundary conditions, this bubble would move to the inflow boundary and therefore stop the calculations. But with adapted conditions the bubble stays in the domain. After a certain time (shortly after
Computation of unsteady vortical flows
>,
235
,iim~lollel.l,~llI
2.0
4.0
S.
'
"
STREAKLINES T = 84.000
(b)
o, ¢j"
";"
I
q T" o
?
X
Z
STREAKLINES T = 219.000 (c)
_____
•
~
#
°
•
>, q T"
°
o e Q e e e e e l ~
N
STREAKLINES T  270.000 Fig. 5
(ac)continued overleaf.
'
236
M. BREtmR e t
al.
(d) q
o
•
,,,4
'"
,,,
.:i.i.
q
7" q
X STREAKLINES T = 330.000 Co) q ~ ,
._______._~

~
q
• ".'.'i"."
>,
"
• . ! ' 2 . ":
O $ I IQIOeIQI~
q
7" q
X STREAKLINES T = 1818.000 Fig. 5. Development of an isolated vortex (initially) to symmetric and later to asymmetric spiraltype breakdown. Streaklines of the flow at different times (Re  200,/] = 0.68). Incompressible NavierStokes solution using a dualtimestepping method• (a) T = 8 4 ; (b) T = 2 1 9 ; (c) T =270; (d) T = 330; (e) T = 1818.
T = 219) the solution becomes asymmetrical due to small disturbances (e.g. by roundoff errors) (T = 270). Later the solution forms the spiraltype of vortex breakdown (T = 330, T = 1818), well known from the experiments. Obviously the spiral solution requires a certain, longer lifetime of the vortex, which could only be achieved with the new lateral boundary conditions. CONCLUSION As demonstrated in the above flow problems, numerical algorithms can influence an unsteady vortical solution to a much greater extent than is usually experienced in steady flow problems. Different discretization schemes, e.g. upwind or central, can yield considerably different results, although the truncation errors are formally of the same order. The numerical dissipation has been shown to be an important influential factor, which should bc minimized in the simulation of flows with smallscale structures. In addition, the resolution of the flow field is of similar importance. The deviation of two solutions with different spatial steps sizes grows with increasing time level, so that totally different solutions are obtained, when only the time of integration is chosen large enough. Boundary conditions have been shown to be critical in cases when vortices cross outflow
Computation of unsteady vortical flows
237
boundaries. Here simple techniques, often used in steadystate c o m p u t a t i o n s , are seldom satisfactory. F o r u n s t e a d y flows the outflow c o n d i t i o n s have to ensure that vortices can pass the b o u n d a r y with a m i n i m u m o f disturbances. T h e i m p o r t a n c e a n d influence o f b o u n d a r y c o n d i t i o n s was also d e m o n s t r a t e d in a different sense, when they were used to stabilize a flow structure in the interior flow field. All examples show that realistic a n d accurate simulations o f u n s t e a d y vortical flows require more care a n d sophisticated solution schemes t h a n simulations of steady flow problems. REFERENCES 1. K. Dortmann, Computation of viscous unsteady compressible flow about airfoils. In Proceedings of the llth Int. Canfi on Numerical Methods in Fluid Dynamics, 1988; Lecture Notes in Physics. SpringerVerlag, Berlin (1989). 2. M. F. Unal and D. Rockwell, On vortex formation from cylinder. J. Fluid Mech. 190, 491 (1988). 3. R. Bouard and M. Coutanceau, The early stage of development of the wake behind an impulsivelystarted cylinder. J. Fluid Mech. 101, 583 (1980). 4. Ta Phuoc Loc and R. Bourard, Numerical solution of the unsteady flow around a circular cylinder. J. Fluid Mech. 160, 93 (1985). 5. D. H~nel, Effects of numerical dissipation in solutions of the NavierStokes equations. Paper presented at the 3rd Int. Canfi on Hyperbolic Problems, Uppsala (1990). 6. M. Meinke and D. H/inel, Simulation of unsteady flows. Paper presented at the 12th Int. Confi on Numerical Methods in Fluid Mechanics, Oxford (1990). 7. M. Meinke and D. Hfinel, Time accurate multigrid solutions of the NavierStokes equations. In International Series of Numerical Mathematics, Vol. 98, pp. 289300. Birkh/iuser, Basel (1991). 8. J.J. Kl6ker, E. Krause and K. Kuwahara, Vortical structures and turbulent phenomena in a pistonenginemodel. Paper presented at the 13th Int. Canfi on Numerical Methods in Fluid Dynamics, Rome (1992). 9. J. J. Kl6ker and E. Krause, Numerical simulation of vortical and coherent structures in compressiblejet flows. Paper presented at the 4th Eur. Turbulence Conf., Delft (1992). 10. J. J. Kl6ker, Numerical influence on the calculation of compressible vortex structures. Paper presented at the Wkshp on Numerical Methods for Simulation of MultiPhase and Complex Flows, Amsterdam (1990). I 1. P. Colella and H. M. Glaz, Efficient solution algorithms for the Riemann problem for real gases. J. Comput. Phys. 59, 264 (1985). 12. J. J. Kl6ker, Numerische Simulation einer dreidimensionalen kompressiblen, wirbelbehafteten Str6mung im Zylinder eines Modellmotors. Thesis, RWTH Aachen (1992). 13. H. O. Kreiss, W. D. Henshaw and L. G. Reyna, On the smallest scale for the incompressibleNavierStokes equations. ICASE Report 888, (1988). 14. F. F. Grinstein, E. S. Oran and J. P. Boris, Numerical simulations of asymmetric mixing in planar shear flows. J. Fluid Mech. 165, 201 (1986). 15. M. Breuer and D. H~inel,Solution of the 3D, incompressible NavierStokes equations for the simulation of vortex breakdown. In Proc. 8th GAMM Conf. on Numerical Methods in Fluid Mechanics; Notes on Numerical Methods in Fluid Mechanics, Vol. 29. ViewegVerlag, Braunschweig (1989). 16. D. H~ineland M. Breuer, Numerical solution of the incompressibleNavierStokes equations for unsteady threedimensional flow. Paper presented at the 3rd Int. Cangr. on Fluid Mechanics, Cairo, Egypt (1990). 17. M. Breuer, Numerische L6sung der NavierStokes Gleichungen fiir dreidimensionale inkompressible instation/ire Str6mungen zur Simulation des Wirbelplatzens. Thesis, RWTH Aachen (1991). 18. E. Krause, The solution to the problem of vortex breakdown. In Lecture Notes in Physics, Vol. 371, pp. 3550. SpringerVerlag, Berlin (1990).
CAF ~ / ~ J