Computation of unsteady vortical flows

Computers FluidsVol. 22, No. 2/3, pp. 229-237, 1993 Printed in Great Britain. All rights reserved

COMPUTATION

0045-7930/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~ract--The Navier-Stokes 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 time-dependent. 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 Navier-Stokes 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 Navier-Stokes 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]. Navier-Stokes solution with:

1, van Leer flux-vector splitting; 2, van Leer, Hgnei flux-vector splitting; 3, ~ t r d differelr~ ~heme, node-centered; 4, cell-vertex scheme; 5, central differencescheme, node-centered 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 large-scale major vortices. The solution for this flow reacts very sensibly to the numerical dissipation, which has a filter-like 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 Runge-Kutta scheme with different spatial discretizations (node-centered central, cell-vertex central, flux-vector splitting). Figures 2(a-f) 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 fourth-order 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 high-frequency 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 flux-vector 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 node-centered 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); flux-vector splitting (van Leer): (e) lift coefficient vs time, (f) lines o f constant pressure.

8

0

b

~.0-1.5-1.0-o.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 large-scale 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, second-order Godunov-type method, based on the PLM method of Colella and Glaz [11]. This method was extended to the complete Navier-Stokes equations, and reformulated recently for 3-D curvilinear coordinates [12]. One of the most important problems for the investigation of large-scale vortex flow is the resolution of the smaller scales. An estimation of the minimum length scale ~'min in 2-D and 3-D 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; Godunov-type 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 first-order 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 so-called non-reflecting boundary conditions are known, which, however, are not really non-reflecting, 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 "non-reflecting" 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 =Pb-1 "~-( P o o - - Pb-1)"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; Godunov-typc PLM method, first-order extrapolation (top), non-reflecting boundary condition (bottom).

234

M. B~tw~ et al.

Compared with Fig. 4 (top), this non-reflecting 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 [15-17], is the 3-D 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, spiral-type breakdown and nearly axisymmetric bubble-type breakdown. The bursting of an isolated vortex was simulated numerically in a 3-D box using a Cartesian mesh. The method of solution is a dual-time-stepping scheme with artificial compressibility in the pseudo time level [16]. This method for the unsteady, 3-D Navier-Stokes equations of an incompressible flow is second-order 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 non-reflecting 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 time-dependent, 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 bubble-type breakdown could be provoked. But due to the self-induced 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 time-dependent 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 steady-state equations in a part near the inflow boundary, with the aim to suppress the self-induced 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 time-dependent 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 self-induced 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 bubble-like 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,~ll-I

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

(a-c)--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 spiral-type breakdown. Streaklines of the flow at different times (Re -- 200,/] = 0.68). Incompressible Navier-Stokes solution using a dual-time-stepping 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 round-off errors) (T = 270). Later the solution forms the spiral-type 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 small-scale 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 steady-state 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. Springer-Verlag, 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 Navier-Stokes 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 Navier-Stokes equations. In International Series of Numerical Mathematics, Vol. 98, pp. 289-300. Birkh/iuser, Basel (1991). 8. J.J. Kl6ker, E. Krause and K. Kuwahara, Vortical structures and turbulent phenomena in a piston-enginemodel. 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 Multi-Phase 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 incompressibleNavier-Stokes equations. ICASE Report 88-8, (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 3-D, incompressible Navier-Stokes 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. Vieweg-Verlag, Braunschweig (1989). 16. D. H~ineland M. Breuer, Numerical solution of the incompressibleNavier-Stokes equations for unsteady threedimensional flow. Paper presented at the 3rd Int. Cangr. on Fluid Mechanics, Cairo, Egypt (1990). 17. M. Breuer, Numerische L6sung der Navier-Stokes 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. 35-50. Springer-Verlag, Berlin (1990).

CAF ~ / ~ J