A dual time-stepping method for 3-d, viscous, incompressible vortex flows

Computers Fluids Vol.22, No. 4/5, pp. 467484, 1993 Printed in GreatBritain.All rightsreserved

A DUAL

TIME-STEPPING

0045-7930/93 $6.00+ 0.00 Copyright~7 1993PergamonPressLtd

METHOD

INCOMPRESSIBLE

VORTEX

FOR

3-D, VISCOUS,

FLOWS

M. BREUERI and D. HANEL2"J" qnstitut fiir Hydromechanik, Universit/it Karlsruhe, Karlsruhe, Germany 2Institut ffir Verbrennung und Gasdynamik, Universit/it Duisburg, 4100 Duisburg, Germany (Received 20 January 1993)

dual time-steppingmethod, combined with artificial compressibility,was developed for the solution of the 3-D Navier-Stokesequations of incompressiblefluids. The concept of artificial compressibility allows the adoption of numerical techniquesoriginallydevelopedfor compressibleflows. Thus, an incompressibleflow solveris described, in which upwind discretizations,implicit relaxationschemeswith flux-matrix splitting and a Rung~Kutta time-steppingschemewith multigrid are used. The computation of a bursting vortex, called vortex breakdown, is presented as one important application of the implicit version of the present method of solution, which of course is not restricted to vortical flows. Abstract--A

1. I N T R O D U C T I O N Viscous incompressible flows play an important role in both engineering applications of hydrodynamics and aerodynamics and in fundamental fluid dynamics. Numerical investigations of such flows are based mainly on solutions of the Navier-Stokes equations. Computational studies have been performed for several decades and various methods of solution for incompressible flows have been developed. However, the progress in computational fluid dynamics in recent years has been much greater for problems of compressible flow than for the corresponding incompressible problems, this may have been caused by the increasing interest in applications of aerospace engineering and other related branches. As a result, a large number of highly developed methods for compressible flows are in use today. However, the methodology for incompressible fluids has gained from this progress and has adopted a number of ideas and techniques originally developed for the computation of compressible flows. One class of such methods, those based on artificial compressibility, is discussed in this paper. In principle, incompressible flows can be computed by compressible flow solvers and choosing a very small Mach number. This procedure, however, becomes more and more inefficient with decreasing Mach numbers because of the limitations of the permissible time-step. Essential differences between the formulation of solution algorithms for the conservation laws of compressible and strictly incompressible fluids are caused by the lack of state equation for the pressure as a function of density. This becomes obvious with the lack of the time derivative of the density in the continuity equation, which results in a singular coefficient matrix for the time derivatives of the system of Navier-Stokes equations. The methods of solution differ mainly in the manner they ensure the missing linkage between the pressure and velocity field, and enforce conservation of mass. For 3-D incompressible flows the algorithms can generally be grouped into three types. These include method based on pressure iteration to satisfy the continuity equation, methods using the vorticity-velocity or vorticity-vector potential formulation and methods based upon artificial compressibility. The concept of artificial compressibility, originally proposed by Chorin [1] for solving steadystate problems, removes the matrix singularity by adding an artificial time derivative to the continuity equation. This renders the incompressible Euler equations totally hyperbolic, which allows the use of the highly developed algorithms of compressible flow for incompressible problems. This concept can be extended to unsteady flows by adding an artificial time derivative to all equations, this is called the dual time-stepping method [2-4]. The approach of artificial compressibility tAuthor for correspondence. 467

468

M

f:t~tv.t~

and D.

Iq,{~t~l

allows the application of all experiences and techniques for compressible flows to incompressible fluid motion. This is done in the present study. The Navier Stokes equations for 3-D incompressible flows were discretized by central as well as upwind schemes based upon approximate Riemann solvers, e.g. the flux-difference splitting of Roe [5]. Furthermore, the method of artificial compressibility permits the application of both explicit time-integration, e.g. a Runge-Kutta time-stepping scheme and unfactorized implicit relaxation schemes. Additionally convergence could be improved by multigrid techniques. The explicit code was combined with the full approximation storage (FAS) multigrid concept [6] and used for steady-state problems as an alternative to the implicit relaxation scheme, mainly used in unsteady problems. The different approaches are described in the next section. Solutions of these schemes are presented for 2-D and 3-D test problems. An important application of the present dual time-stepping method is the computation of the 3-D unsteady flow of a bursting vortex, called vortex breakdown. For this highly unsteady flow problem the implicit version was chosen, since this method satisfies the essential requirements of accuracy and numerical stability for such a problem. 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. In contrast to many breakdown investigations, the present study is concerned with the fully 3-D breakdown process. In order to concentrate on pure vortical flow, and to reduce the computational work, an isolated vortex embedded in an inviscid irrotational flow is considered. In this paper the formulation of the initial and boundary conditions is described and several computed results for different types of vortex breakdown are presented. 2. METHODS OF SOLUTION 2.1. Governing equations

Incompressible viscous fluid motion is described by the continuity and momentum equations, forming the system of Navier---Stokes equations. The divergence form of the conservation equations in general curvilinear coordinates (~, r/, ~) reads: /~* 0 , + ~ + ( ~ , + / ~ ; - R e

Q=J

'Q=J

'

,

-~ *(kv~ + Sv0 + Tv:) =0,

R=

kL'

(1)

iiooo 1

0

0

0

1

0

0

0

1

U uU + ~xP v U + ~,.p wU + ~.p

p=j-I

Rv = J

i

F g~ u~ + g20u, + g3 u¢ I g~ v~ + g2v~ + g3v: t.g t w, + g2w,7 + g31,ve, I

u V + qxP v V + qyp

~=j-i

"

w V + rhp

i2i = j

,

1 V W + ~yp I w W + ~:p J

U = u C + v~, + w~:,

~v=j

0

g2 u, + g4 Un + g5 t,t( I ] g2v~ + g4v, + gsv~ ~-g2 w, + g4 Wn+ gs W;

T v = J -I

F

o

g3u¢+gsun+g6u~

1

,

I g3v~+gsv'+g6v¢ I t. g3 we + g5 w~ + g6 w; .J

V = uqx + Vtly + wtl.,

W = u~x + v(,y + w~.

A dual time-stepping method for 3-D vortical flow

469

and g4=

g6 =

Herein Q is the vector of the primitive variables, F, (~ a n d / 1 are the inviscid (Euler) fluxes and R~, S~ and irv define the viscous terms. J is the metric Jacobian and the variables gl-g6 are functions of the metric. R corresponds to the coefficient matrix of the time derivatives and is singular, requiring a velocity field, which is free of divergence at each time level. All quantities are made dimensionless, the Reynolds number, Re, is defined with the corresponding reference quantities. 2.2. Concept of artificial compressibility

Due to the assumption of constant density, a time derivative of the pressure or density does not appear in the continuity equation. Hence the continuity equation requires a special treatment compared to the momentum equations. A common solution procedure for all the conservation equations, as is usual for compressible flows, is not applicable. The concept of artificial compressibility [1] removes this drawback by adding an artificial time derivative of pressure to the continuity equation. This yields the modified continuity equation: 1@

/3~ & + V.v = 0,

(2)

where fl is a free parameter, which corresponds to an artificial speed of sound, according to an equation of state p =//2p, similar to compressible flows. The artificial compressibility alters the type of inviscid equations from parabolic in the time-space plane, which means infinite wave speed, to totally hyperbolic with finite speed c = ~ / ~ 2 . Thus, the system of equations has similar properties to that of compressible flow. The original approach of artificial compressibility is consistent with Navier-Stokes equations at steady-state only, where the velocity field is free of divergence. The concept can be extended to unsteady flows by adding artificial time derivatives OQ/& to all equations, which is called dual time-stepping. Applying the artificial compressibility approach only to the pseudo-time z, then the physical time behavioiur is not influenced as long as the solution converges in the pseudo-time to a "steady-state" for each physical time level t. The new system of equations now reads:

with /~= diag(0,1,1,1),

/~=dlag

,1,1,1

.

The difference between equations (1) and (3) is given by the additional term / ~ , Q~, which corresponds to the derivative of the dependent variables with respect to the pseudo-time r. In contrast to the coefficient matrix/~, the matrix/~ is not singular. If the term Q~ converges to zero, the original equation (1) is solved. Setting / ~ = 0 provides a system of differential equations appropriate for solving steady-state problems. Then integration in pseudo-time can be carried out by explicit (here a Runge-Kutta time-stepping scheme) as well as implicit schemes (here an unfactorized relaxation algorithm). If, however, unsteady solutions are of interest (~ ¢ 0), only the implicit method is applied to the pseudo-time level, allowing much greater artificial time steps and leading to a better efficiency. 2.3. Spatial discretization

For spatial discretization of the conservation equations a finite-difference approach is employed to the system of equations (3) in general curvilinear coordinates. The control volume is defined as node centred, where the variables are given on the node and the control volume is arranged around it. All the derivatives of the fluxes, e.g. F¢, are approximated by p~ = Fi+ ,/2.j.k - Fi l/Zj.k A~ '

(4)

470

M BR}!U~iRand D. tt~,N~

where /e± 12,/,k are the fluxes at the cell interlaces. The other differences are written in the same way. In the following, the approximations for the Euler and viscous fluxes of the physical part are considered. The viscous fluxes are updated by means of central differences O(Ax2), which fits the elliptic nature of the viscous effects. The formulation of the inviscid mass and momentum fluxes has a strong influence on the properties of the solution, in particular for high Re flows. Three alternative approaches are chosen and tested: central differencing; the QUICK scheme (Quadratic Upstream lnterpolar for Convective Kinematics [7]); and a simplified flux-difference splitting of Roe [5], In the case of central differencing the Euler fluxes on the cell interfaces are evaluated as functions of the arithmetically averaged primitive variables. The scheme is of second-order accuracy. To suppress high-frequency error components, artificial damping terms in the form of fourth-order differences have to be added to the steady-state operator. The other alternative of projecting the variables to the cell interfaces is given by the QUICK scheme, proposed by Leonard [7]. The QUICK scheme is an upwind discretization, where upwinding is controlled by the contravariant velocities. A formulation of third-order accuracy is chosen for the convective parts, and a second-order central formulation for the remaining terms. In the case of high Re flows, a fourth-order damping term has to be added to the continuity equation, due to the use of central differences there. The third approach uses an upwind scheme, derived directly l¥om compressible solvers. Within each control volume there exists an averaged value of the variables, and hence the values at the interfaces of neighbouring cells show jumps. According to the theory of non-linear hyperbolic Euler equations, these jumps generate local Riemann problems at the cell interfaces. This Riemann problem is solved approximately using the idea of flux-difference splitting of Roe [5]. The original approach by Roe can be simplified for incompressible flow by taking into account constant density in the averaging procedure for the variables of the flux Jacobian matrices and by the modified characteristics. Second-order accuracy is achieved using the MUSCL approach [8] to interpolate the variables to the cell interfaces.

2.4. Implicit relaxation method An implicit scheme enables unrestricted time steps and is therefore advantageous for computations of both steady and unsteady solutions. The implicit scheme is derived from equation (3) with the approximations of equation (4) by defining the fluxes at the new physical time level n + 1 and the new (iteration) level v + 1 of the artificial time r. The Euler fluxes and the viscous fluxes are expanded in Taylor series and linearized with respect to the pseudo time ~:

~"~ .... '=F"~','+A"~"~AQ ' w i t h A = ? ~ , /~"+~,'+~=[{~+~,~+L"~~AQ '

withL=?~Q

(51)

Herein A is the Jacobian matrix of the Euler flux P and L is the Jacobian of the stress component *~v. The corresponding Jacobians are derived for the other terms. With the linearization of the fluxes resulting quasi-discrete delta-form of equation (3) reads:

J-~/Ar+cTd(A-Re

J L ) + =~rt -(B-Re

= _ [ ~_j

IM)

+ ~[( C - R e

,AQ''L - A i ..... + ( P ~ + d , , + / - ) : - R e

~N) AQ' '.(k~+$,

+T,~

))"*','] _"

(6)

This system of equations can be written in a short form: LHS' AQ' = RHS.

(7)

The implicit operator (LHS) represents the "numerical" part of equation (6) and the RHS of equation (7) approximates the "physical" part, the discrete Navier-Stokes equations (1) for unsteady incompressible flow.

A dual time-steppingmethod for 3-D vortical flow

471

Since the converged solution is independent of the implicit part, the derivatives of the LHS and RHS can be discretized separately and in a different way. While the discretization of the implicit operator controls the convergence, the discrete formulation of the (physical) Navier-Stokes equations has to be sufficiently accurate in space and time. The different spatial discretizations of the fluxes for the RHS of equation (6) are discussed above. Accuracy in physical time (index n) on the RHS is achieved by a second-order-accurate formulation in time for Qt: AQ,+, 3 Q , + , , , ' + , _ 4 Q , + Q .-~ Q'-* At 2At (8) The discretization of the LHS of this equation, called the numerical part, can be manipulated to achieve best convergence. For the discretization in pseudo-time a formulation of first-order accuracy is sufficient: AQ" Qn+l,v+l-Q,+ i,,. Q~ ~ - - (9) Ar Az The aim of the spatial discretization of the LHS is to achieve a diagonal-dominant solution matrix, which guarantees the stability of the relaxation method. An upwind formulation of first order combined with a matrix-splitting technique is used presently. For this reason the Jacobian matrices are diagonalized by an appropriate transformation. (10)

A = M A M -1,

where M is the eigenvector matrix of A and A is the diagonal matrix A = diag(2,, 22,23, 24). The diagonal matrix A can be split into parts with positive and negative eigenvalues, A + and A -, where the eigenvalues of A -+ are defined as: A=A++A

'

2+_L+ILI 2 '

(11)

where

2.,2=u, ,13=u+c, 2,=u-c, c = ~ [ F . The transformation back to the Jacobian splits the matrix A into two matrices A + and A - with positive and negative eigenvalues: A=MA+M-I+MA

M

t=A++A

.

(12)

All the flux Jacobians of the LHS of equation (6) are substituted by split matrices in the same way and the corresponding contributions are updated with values extrapolated from the left and right, respectively. This splitting results in the first-order upwind formulation of the implicit operator and thus in a diagonal-dominant solution matrix. The inversion of this large system of difference equations is achieved by an implicit line GaussSeidel relaxation method in each physical time-step. The solution of one physical time-step requires of the order of 10 pseudo-time steps for sufficient accuracy in time. Despite the recursive algorithm, the 3-D algorithm could be highly vectorized [9]. 2.5. Explicit Runge-Kutta method and multigrid acceleration

The explicit method, based on a Runge-Kutta time-stepping scheme, has been derived for steady-state solutions of the Navier-Stokes equations. An explicit scheme is well-suited for implementation on different vector and parallel computers because of its simpler structure, moreover with additional acceleration techniques it becomes comparable to an implicit scheme with respect to the computational time. At present a Runge-Kutta method is used, which has been successfully employed on the Euler and Navier-Stokes equations for compressible flows by Jameson [10] and other authors. A general solution scheme for an N-step Runge-Kutta method, according to the pseudo-time r (index v),

472

reads:

M~ BREUER and D. HAN'EL

Qi0~ = Q ~

"~

Q~l)=Q°-cttArJ Q(,+ 1) QN

Res(Q ~ t l ) ) ~

1 = 1 , 2 . . . . . N.

(13)

.)

=

The residual of Res(Q) is equivalent to the discretized steady-state operator, and corresponds to the RHS of equation (6), but without the time derivative. A 5-step scheme was adapted for maximum Courant number with the coefficients et = (0.25, 0.1667, 0.375, 0.5, 1.0) for central spatial differencing and e~= (0.059, 0.145, 0.273, 0.5, 1.0) for upwind schemes. Local time-steps are used to accelerate the convergence to a steady-state solution, determined by the local stability limit and constant Courant number. Further improvement of convergence was achieved by incorporating the FAS multigrid concept proposed by Brandt [6]. The present multigrid algorithm is based on an algorithm for 2-D flow problems [11, 12], which is extended here to 3-D incompressible flow. Considering a grid sequence Gk, k = 1. . . . , m, with step sizes hk = 2hk~ ~, a finite-difference approximation on the finest grid G,, may be

LmQ,,=O,

(14)

where LmQm corresponds to the discretized conservation equation (3). If the solution on the fine grid is sufficiently smooth with respect to the high-frequency solution components, equation (14) may be approximated on a coarser grid Gk ~ by a modified difference approximation Lk IQk i = ( r ~ ' l)Re~,

(15)

where (r)Res is the "fine-to-coarse defect correction" of the residual. It preserves the truncation error of the fine grid G,, on the coarser grids Gk ~ and is defined by (r~'

I )Res = ('C?)Res "Jr t k

1(I~, IQk ) -

I I ~ '(LkQk),

(16}

where the restriction operators from grid G~ to Gk L mean injection (I~ ~) and full weighting operators (Hi-~). The solution sequence used here is a V-cycle. Starting with the solution of equation (14) on the finest grid G,,, the correction equation (16) can be calculated. On the coarser grid G,, equation (15) is solved. This cycle is repeated until the coarsest grid is reached. Then the variables are interpolated back to the finest grid with some solution steps on the coarser grid levels in between. According to the FAS scheme, only the correction between the "old" fine-grid solution Q~J~ and the "new" coarse-grid solution Qk is interpolated to the fine grid to update ~+tr~°~d• Q kold +l

r~ old+ l + I~ +l(Qk ~k

- - ~rk k + l ~ kr~old +lf"

{t 7 )

Herein If +~ is the interpolation operator, where linear interpolation is used. An important influence on the multigrid convergence is given by the treatment of the boundary conditions on the coarser grids. A coarse-grid correction of the boundary conditions is therefore defined similar to that in Ref. [11], which is consistent with the FAS procedure, equation (15),

2.6. Test results The properties of the explicit and implicit method combined with the three alternative formulations for the Euler fluxes were tested by means of validated flow problems. The boundary layer flow over a flat plate, for which the exact solution of Blasius is available, represents a simple but nevertheless significant test problem to check the spatial accuracy of Navier-Stokes solvers. For high Re the numerical solution reacts very sensitively to different schemes for the inviscid fluxes and their specific damping properties. The flow field was calculated for Re = 105 using the three distinct Euler schemes. All solutions coincide well with the exact one. Comparison of the computational time required for the explicit and implicit algorithms provides a factor of about 5-10 in favour of the implicit version, depending on the applied Euler scheme. The flux-difference splitting shows the best rates of convergence followed by the QUICK scheme, indicating slightly stronger numerical dissipation. To what degree an explicit code can be

A dual time-stepping method for 3-D vortical flow

(a)

473

(b)

0

0

I O

SG

I

go

~'o.

O

2G

O

F

2G

q

O

~2

O-

I

I

Q

O

I ,

1.0

,

i

20

3.0 WU

i

4.0 "103

5.0

1.'0

2:0

3'.0 WU

4'.0 *i03

Fig. 1. Computation of the boundary layer flow over a flat plate with Re = 105 and 33 × 17 grid points. History of convergence, plotted as root-mean-square of the residual (RMS) vs work units (WU) for single grid (SG) and multigrid (2G, 3G). (a) Central discretization; (b) Q U I C K scheme.

accelerated by the FAS multigrid technique was first investigated for this test case. The histories of convergence for a computation with 33 x 17 grid points are given in Fig. 1 using central differences with a fourth-order damping term or the QUICK scheme. The single-grid method (SG) is compared to the multigrid approach with two (2G) or three (3G) grid levels. In Fig. 1 the root-mean square (RMS) of the residual is plotted as a function of work units (WU), which corresponds to the computational work of one Runge-Kutta time-step on the finest grid. The comparison between the single-grid solution and the multigrid indicates significant improvement in the rates of convergence by the multigrid technique for both schemes. The differences in the multigrid acceleration between the two different Euler schemes, however, is a matter of smoothing properties. The central-difference scheme, with a relative large damping factor of E = 2, shows better rates of convergence than the QUICK scheme. This is due to the fact that the multigrid technique requires certain smoothing properties with respect to the high-frequency solution components. Therefore, to achieve optimal speeds of convergence, the multigrid parameters have to be adapted properly. Summarizing the results, the explicit Runge-Kutta method accelerated by the FAS multigrid concept can attain rates of convergence of the same order of magnitude as the implicit version. The laminar flow through a square duct with a 90 ° bend was computed as a 3-D steady-state test case. This flow offers a good test problem for full Navier-Stokes solvers, since strong secondary flow phenomena occur both in the corner regions and by the effects of curvature in the streamwise direction. This flow problem was defined and studied experimentally by Humphrey et al. [13], and numerically by Rogers and Kwak [3]. In the present computation a grid with 65 x 33 x 33 grid points was used with up to 5 grid levels in the multigrid solution. The straight inflow section before and the outflow section behind the bend was set to a dimensionless length of 5. The curved section has a radius of curvature of the inner wall of !.8. Re = 790, based on the mean velocity at the entrance of the duct and the side length of the square cross-section. At the inflow plane the velocity is prescribed corresponding to a fully developed laminar flow in a straight square duct. The computed results are compared with the experimental results [13] and the numerical ones [3], as shown in Fig. 2. The streamwise velocity profiles through the symmetry plane ( x - y ) are plotted at two streamwise stations of the bend, 0 = 0 ° and 0 = 90 °. Good comparison were found between the two computations, whereas the agreement with the experimental data is not quite satisfactory, particularly for the latter part of the curved section. The temporal accuracy was verified for time-dependent problems, like the temporal evolution of the separation zone behind an impulsively started circular cylinder. Figure 2 shows the length of

5.0

474

M. BREUERand D. HANEL ?

(a)

0=0

o

~0

g (5

C~

c5-

0

o'.5

i.'o

2'.o U

(b)

© ,-4

?3-..

0 = 90 °

c0

°

% •

0

", 0

c5 0 ",

c5

d

e. Cb N I0

0 '{5

1 '0

115

~ ]I0

~ .5

U Fig. 2. Computation of the flow through a square duct with a 90' bend. Streamwise velocity profiles at two positions (symmetry plane) in the sections 0 = 0" and 90% with Re = 790. O, Experiment [13]; - - - , computation [3]; , computation. the separation bubble behind the cylinder as function of time. In this figure different numerical results are compared with experimental ones [14, 15]. Curves 1 and 2 are results from the present method, using the Q U I C K scheme and flux-difference splitting, respectively. The splitting slightly underpredicts the length of the bubble, which once more indicates higher dissipation. The same tendency is seen in curve 3, which is calculated with the compressible Navier-Stokes equations for M a = 0.1 [11] using a R u n g e - K u t t a scheme and flux-vector splitting. Although the deviation is small in the present example, the dissipative effects of the splitting methods [curves (2) and (3)] become more evident in highly unsteady flows like vortex shedding at higher Re [16]. A further comparison is made with results by Rosenfeld et al. [17], curve 4, which agree very well with the Q U I C K scheme, curve (1) and the experiments. 3. C O M P U T A T I O N

OF V O R T E X B R E A K D O W N

The present method of solution was used to study 3-D unsteady flow of a bursting vortex, called vortex breakdown. This phenomenon is of particular interest for the flow around delta wings,

A dual time-stepping method for 3-D vortical flow

475

3'

Length D 2

I

0

i

,

i

,

,

i

2

~

6

8

10

12

1/.

Time U~/D Fig. 3. Comparison of numerical and experimentalresults for the temporal developmentof the length of the separation bubble behind an impulsivelystarted circular cylinder with Re = 42.8. C), Coutanceau and Bouard [14]; O, Collins and Dennis [15]; 1, dual time-steppingmethod (QUICK scheme); 2, dual time-stepping method (flux-differencesplitting); 3, compressible Navie>Stokes equations (RungeKutta method, flux-vectorsplitting); 4, fractional step method, Rosenfeld et al. [17] (staggered mesh). because of the drastic influence on the aerodynamical behaviour of the flow. Therefore, prediction and control of breakdown on wings at high angle of attack plays an important role in critical flight conditions. Whereas breakdown is an undesirable effect on wings, the intense mixing taking place in the flow field of a bursting vortex can be used in combustion chambers for rapid mixing of fuel and air. Furthermore, the breakdown bubble can be employed as an aerodynamical flame holder whereby no structural part is necessary. Other examples where the breakdown phenomenon is observed are the flows in tornadoes and diffusers of water turbines. Faler and Leibovich [18] studied the breakdown phenomenon in a vortex tube. They observed 6 different types of breakdown modes, the bubble-type and spiral-type being the most important ones. The bubble-type is characterized by the formation of a breakdown bubble, which has a nearly axisymmetric outer form, whereas the flow field in the bubble is not necessarily symmetrical. The spiral-type of breakdown is marked by an abrupt kink followed by a corkscrew-shaped twising of the vortex filaments. A survey of experimental investigations is presented in Ref. [19] and the most important results are reviewed in Ref. [20] by Leibovich. The phenomenon of vortex breakdown has been studied in many theoretical and numerical investigations, but no general agreement as to the explanation of its fundamental nature has been achieved. Numerical investigations were made under the assumption of axisymmetric flow [e.g. 21-23], which seems to be too restrictive and only applicable for bubble-type vortex breakdown. 3-D investigations have been conducted by Liu and Menne [24, 25] with a Fourier decomposition, by Spall and Gatski [26] for the bubble-type and by Leibovich and Stewartson [27] for the spiral-type of breakdown. 3.1. Initial conditions

The present study is concerned with the fully 3-D breakdown process. An isolated vortex, embedded in an inviscid irrotational flow, is considered in order to concentrate the investigation on the physical details of the vortical flow and to reduce the computational work. As initial conditions the circumferential velocity distribution, independent of the axial coordinate z, is given by a polynomial in the radial coordinate r: ~ r ( 2 - r 2) vo(r) = fi * g(r),

g(r) = ~'1

forr~

1.

(18)

k--

The swirl parameter fi is set to 0.68 in all calculations. This distribution corresponds to solid body rotation in the core and to a potential vortex outside. The initial field of the axial velocity is given by . .

fw0+bz+cz 2+dz 3

w(z) = l

l+6w/az.(z-zl)

forz>~zL forz>zl,

(19)

476

M BJai~:~Rand D, Hay,Lt

where w ( z ) is assumed to be constant in the radial direction. The field of the initial pressure can be computed from the Poisson equation f\)r the pressure. The polynomial in equation (19) represents a slight acceleration to stabilize the flow near the inflow boundary. In the further downstream part, z > z~ = 4, the flow is decelerated with a constant velocity gradient 6w/6z in order to provoke breakdown by an imposed positive gradient i)1 the axial direction. 3.2. Boundar): conditions

Essential problems arise, similar to ira other unsteady vortical flow simulations, from the formulation of the boundary conditions. For the present case, boundary conditions and compatibility conditions, respectively, have to be prescribed in the inflow and outflow planes and on the lateral boundaries. In the inflow plane the distribution oF axial and circumferential velocity is prescribed, all 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 [9,28] 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 distant effects of a nozzle contour or of a delta wing. They can be used to control the axial and temporal development of the vortex breakdown, 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, 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 toward the inflow boundary, as also observed, for example, in Refs [26, 29, 30], which makes it impossible to carry on the numerical simulation. Therefore, a new time-dependent condition for the lateral boundaries was developed by Krause [28] 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 steady-state solution in a part near the inflow boundary, such that the self-induced upstream movement will be suppressed. Thus, the pressure at the boundary (z <~z~, r = R) for the region near the inflow plane is written as p"+ J(z, R ) = p"(z, R) + :¢ 'Ap(z).

The pressure correction Ap(z) =

R\ ~v!}: t )/ . . . .dr. + . . . . . )

)

= dz' 0

is an approximation for 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 c~ = 0 means fixed, prescribed pressure at the lateral boundary, whereas a = l forces the flow to become stationary. In the present calculations a small value of e ~ 0.1 was chosen to obtain a smooth change from the steady flow near the inflow boundary to the unsteady flow further downstream. How this adaptive boundary condition acts on the development of vortex flow is demonstrated in Fig. 4. This figure shows the axial variation of the axial velocity at the lateral boundary and on the axes for different times. In the first phase, where the bubble mode is formed, the bubble moves upstream to the entrance (curves 4--6); then the boundary conditions react and the main flow is accelerated (curves 7 and 8), which avoids further upstream movement. After the change from the bubble to the spiral mode, the velocity distribution approaches the initial distribution. This adaptive boundary condition is necessary, due to the fact that an appropriate pressure distribution providing a breakdown process within the domain of integration is a priori not known and therefore has to be determined iteratively. Applying this lateral boundary condition it was possible to avoid self-induced movement of the burst vortex, and thus longer observation of the flow structures in the vortex was possible.

A dual time-stepping method for 3-D vortical flow

477

1.0 0.8 0.6 0.4 0.2 ~:

0.0 -0.2 -0.4 -0.6

,

i , 2

0

l

,

l

4

~

l

6

,

l

8

, l 10

, l 12

, I 14

, 16

18

Fig. 4. Variation of the axial velocity distribution at the lateral boundary of the computational domain and on the axis due to the adaptive boundary conditions for the flow during the transition from bubble-type to spiral type vortex breakdown (as shown in Figs 7 and 8); Re = 200, fl =0.68, 6w/Sz = -0.0165. Time levels: 1, T = 84; 2, T = 132; 3, T = 165; 4, T = 201; 5, T = 219; 6, T = 240; 7, T = 300; 8, T = 330; 9, T = 390; 10, T = 510.

3.3. Numerical results of vortex breakdown F o r the c o m p u t a t i o n of the vortex flow a Cartesian grid with 41 x 41 x 60 grid points is used, the c o m p u t a t i o n a l d o m a i n is b o u n d e d by - 5 ~< x ~< 5, - 5 ~
¢i"

o. 7

,o,|,$0QI g ~

?•

.

8.0 Z

X STREAKLINES T = 60.000 Fig. 5---continued overleaf•

10.0

12.0

14.0

16.0

18.0

478

M. B t o c i ; a and D H,KNEL

i

i

°,. jslo~ In

¢" i

r

t

'

-1.0

I'(]

X "

~ ~ ' ~

a.o

......

~,9

'

zo

40 •

S,O

Z

STREAKLINES T = 72.000

(c) i

L

oi

~d

oi I

;

oi

!

I [

I

oi

~

i-+

-i~0

'/

.

a.0

-~ 5.0

z0

•

STREAKLINES T = 84.000

(d) q

i:

'!

i

oi ¢~ -L

I!

-'~ X

L0

3,0

~.0

2.0

4.o

•

Z

STREAKLINES T = 99.000 Fig. 5. Streaklines at four different time levels for the t r a n s i t i o n to b u b b l e - t y p e vortex b r e a k d o w n ; Re = 200, fi = 0.68, 6w/bz = - 0 . 0 3 . (a) T = 60; (b) T = 72; (c) T = 84; (d) T = 99,

A d u a l t i m e - s t e p p i n g m e t h o d f o r 3 - D vortical flow

479

(a) 5.0

A

Y 4.0

B C

3.0

D E

2.0

F G

1.0

N 0 P

=-0.7678 =-0.9073 =-1.0469

H I

0.0

J K

-1.0 -2.0 -3.0 -4.0 -5.0

L M

1.0469 = 0.9073 = 0.7678 = 0.6282 = 0.4886 = 0.3490 = 0.2094 = 0.0698 =-0.0698 =-0.2094 =-0.3490 =-0.4886 =-0.6282

0.0

2.0

4.0

6.0

8.0

=

I 0 . 0 12.0 14.0 16.0 16.0

Z (b) 5.0

A B

= 1.3445 = 1.1652 C = 0.9859 D = 0.6067 E = 0.6274 F = 0.4481 G = 0.2689 H = 0.0896 I =-0.0896 J =-0.2689 K =-0.4482 L =-0.6274 M =-0.8067 N =-0.9860 O =-1.1652 P =-1.3445

Y 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0

0.0

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

180

(c) 5.0

A B C D E

F

= = = = = =

1.0

G H

= 0.3044 = 0.1015

0.0

I J

=-0.1015 =-0.3044 =-0.5073 =-0.7102 =-0.9131 =-1.1160 =-1.3190 =-1.5219

Y 4.0 3.0 2.0

K

-1.0

L M N 0 P

-2.0 -3.0 -4.0 -5.0

0.0

2.0

4.0

610

8.0

10.0

1210

14.0

16.0

1.5219 1.3190 1.1160 0.9131 0.7102

0.5073

18.0

(d) 5.0

A B C D E F

Y 4.0 3.0 2.0

0

= 1.5258 = 1.3224 = 1.1190 = 0.9155 = 0.7121 = 0.5086 = 0.3052 = 0.1017 =-0.1017 =-0.3052 =-0.5086 =-0.7120 =-0.9155 =-1.1189 =-1.3224

P

=-1,5258

G H I J

1.0 0.0

K

-1.0

L

-2.0

M N

-3.0 -4-0 -5.0

0.0

2.0

4.0

6.0

6.0

10.0

12.0

14,0

16.0

18.0

Z Fig. 6. L i n e s o f c o n s t a n t v o r t i c i t y co~ in the y - z c r o s s - s e c t i o n a t x = 0 f o r the flow c o r r e s p o n d i n g to Fig. 5.

480

M. BREUERand D. HXr,;I-;I.

swelling of the vortex core takes place on the axis of the vortex, combined with the development of a free stagnation point. The breakdown bubble evolves from the swelling of the vortex core. At later times the bubble grows and moves upstream near to the inflow section (T = 99). In several experimental investigations the breakdown process was studied in a vortex tube. Sarpkaya [31] found that the axisymmetric breakdown evolves either from a double helix, or from a spiral, or directly from an axisymmetric swelling of the vortex core. The model of evolution depends on the Re and circulation number, whereas for sufficiently high values only the last form is observed. Escudier [19] documented the temporal evolution of the breakdown process by a series of photographs, showing the swelling of the core, the formation of the first small breakdown bubble. whose dimensions increase in time, and the movement of the bubble towards the inflow section. Although the conditions of the numerical and experimental investigations are not identical, the results are qualitatively in close agreement. Of particular interest is the structure of the internal breakdown region, e.g. the vorticity distribution. In Fig. 6 the distribution of the vorticity component e)~.normal to the y z plane (x = 0) is plotted for the same time levels as shown in Fig. 5 For T = 60, the breakdown regions consists of a single vortex ring. With increasing time the ring moves upstream and its strength grows. For T = 99, the formation of a second counterrotating vortex ring, which lies inside the first ring, is indicated by the lines of constant vorticity. In principle the internal structure with the two counterrotating vortex rings is comparable with the experimental observations of Leibovich [20], who found a similar double-ring structure. At later times the flow mode changed from the axisymmetric breakdown to the asymmetric spiral mode. Without adapting the lateral boundary conditions, the bubble would move to the inflow boundary and therefore would stop the calculations. A typical sequence of streaklines, for which the change from the bubble to spiral mode has been observed, is shown in Fig. 7 for a smaller initial velocity gradient 6 w / 6 z = --0.0165. Starting again from the initial straight vortex (T = 0), the vortex swells (T = 84) and forms an axisymmetric bubble-like breakdown configuration. After a certain time (shortly after T = 219) the solution becomes asymmetric (T = 270) due to small disturbances (e.g. round-off errors). Later the solution forms the spiral-type of vortex breakdown (T = 330), which corresponds to the experimentally observed spiral-type mode. The sense of the spiral is opposite to that of the basic flow, whereas the spiral itself processes with the flow. Obviously the spiral solution requires a certain, longer lifetime of the vortex, which could only be achieved with the new lateral boundary conditions. Figure 8 shows lines of constant vorticity cox in the y - x plane for the same time levels as in Fig. 7. The flow field is highly asymmetric. In the plotted cross-section it consists of a number of vortex rings, at staggered locations. The centres of the vortex ring represent the points where the corkscrew-shaped vortex pushes through the cross-section. The vortex core is preserved for several (a) 0

°

o.

i

~e.o u -",v 10

30

5.0

20

40

60

Z

STREAKLINES T = 84.000 Fig. 7--Continued opposite.

8.o

A dual time-stepping method for 3-D vortical flow

(b)

48l

o

x

•

3.0

STREAKLINES T = 219.000

• f A-

o. 7"

,°

?j

J

J

STREAKLINES T = 270.000

(d)

o. >. 7

q

fX

~.o

Z

STREAKLINES T = 330.000 Fig. 7. Streaklines at four different time levels for the transition from bubble-type to spiral-type vortex breakdown; Re = 200, ~ ~ 0.68, 6w/fz = -0.0165. (a) T = 84; (b) T = 219; (c) T = 270; (d) T = 330, C A F 22-4/~--F

482

M, BREL'ER a n d D. H~NEt~

(a) 5.0 Y

4

A B C D

= 0.6721 = 0.5825 = 0.4929 = 0.4032

1.0 i

~

= 0.2240 = 0.1344 = 0.0448

4.0

30

oo I

=-o.o448 =-0_1344

-1.0

~

-2.0

M N

--3.0 i i -4.0

0 P

=-0.2240 =-0.3136 =-0.4032 =-0.4929 =-0.5825 .... 0.6721

A B C D E F G tt I J K L M N 0 P

= 1.4024 = 1.2154 = 1.0284 = 0.8415 = 0.6545 = 0.4675 = 0.2805 = 0.0936 =-0.0934 =-0.2804 =-0.4674 =-0.6544 =-0.8413 =-1.0283 =-1.2153 =-1.4023

-5.0~

0.0

2.0

4.0

60

8.0

10.0

120

14.0

16.0

18.0

Z (b) 5.0 Y

4.0 3.0 2.0

1.0 0.0

-1.0 -2.0 -30

-40 -5.0

0 0

2.0

40

6,0

8.0

10.0

12.0

14-0

16.0

18.0

Z

(c) 5.0

Y

3.0

A B C D

20

E F

1.0

H

4.o

G

I J K L M N 0 P

0.0 -1.0 -2.0

I

-3.0 -4.0 -5.0

= = = =

1.4282 1.2365 1_0447 0.8530 = 0.6612 = 0.4694 = 0.2777 = 0.0859 =-0.1059 =-0.2976 =-0.4894 =-0.6812 =-0.8729 =-1.0647 =-1.2565 =-1.4482

i

0,0

2.0

4.0

6.0

8.0

10.0

12.0

14,0

16.0

180

(d) Y

5.0 4.0

1 / l

3.0 2.0

~ ~

A B C D E F G

1.0 0,0

H

H I J

-1.0

i

K L M N 0

-2.0 -3.0

-

P

-4.0 -5.0

0.0

2.0

4.0

6.0

8,0

10.0

120

14.0

16.0

= 0.8586 = 0.7424 = 0.8263 = 0.5101 = 0.3939 = 0.2777 = 0.1616 = 0.0453 =-0.0708 =-0.1870 =-0.3032 =-0.4194 =-0.5356 =-0.6518 =-0.7680

=-0.8841

18.0

z Fig. 8. Lines o f c o n s t a n t vorticity ~ox in the y - z cross-section at x = 0 for the flow c o r r e s p o n d i n g to Fig. 7.

A dual time-stepping method for 3-D vortical flow

483

o

• • ..t,..

o

., ",...'~'... ~..' .'

.,.~ ?'t:.~-.". ,-" '.6" -'0. '..

IOblOBOm4 1 # ~

o

T" q

X Fig. 9. Streaklines of a fully-developed spiral-type vortex breakdown; T = 1818, Re = 200, fi = 0 . 6 8 , ,~w/~z = - 0 . 0 1 8 .

twists of the spiral. If the mean flow is rotating counterclockwise a spiral with a clockwise sense is formed, which rotates with the mean flow. This causes the impression that the vortex rings in the cross-section move downstream with the flow towards the outflow section and leave the domain of integration. Also, at later times the flow field remains unsteady and asymmetric. Figure 9 shows a fully-developed spiral-type vortex breakdown at the time level T = 1818. More computational details and discussion, especially on the causes and physics involved in the breakdown process, which was not the main topic of this paper, can be found in Refs [4, 9, 28].

4. CONCLUSION The computation of 3-D vortical flows requires accurate and efficient methods of solution. The essential properties of a method of solution for the governing Navier-Stokes equations at high Re are determined by the formulation of the inviscid terms of equations. Among the different solution concepts for 3-D unsteady flows of an incompressible fluid, the method of artificial compressibility and its extension to unsteady flows by the principle of dual-time stepping has shown to be wellsuited to satisfy the requirements. In particular, this concept possesses a greater flexibility to adopt the highly developed algorithms of compressible flows than other methods. The present algorithm is written for general curvilinear coordinates and discretized on an unstaggered mesh. Second-order accuracy in space is achieved by employing three different schemes for the Euler fluxes--one central formulation and two different upwind schemes. An explicit Runge-Kutta method accelerated by the FAS multigrid technique is applied as well as an unfactorized implicit relaxation algorithm. The applicability of this method has been verified for several steady and unsteady flow problems and it proves to be a reliable and accurate solution method. For the explicit method of solution the FAS multigrid concept was implemented, which has shown a significant improvement in the convergence for steady-state solutions. The method of solution was used to study the 3-D unsteady flow of a bursting vortex. To enable the computation of this vortex flow over a long time period, this being the principal reason for the failure of many previous axisymmetric investigations, appropriate boundary conditions for the lateral boundary as well as for the outflow plane have been developed. The computations show both main types of vortex breakdown: nearly axisymmetric bubble-type as well as highly asymmetric spiral-type breakdown. The vortex structures agree in many details with the experimentally observed flow fields. The present algorithm represents a powerful tool for obtaining a detailed insight into the physics of the breakdown process as well as clarifying the causes for the transition from the bubble-type to spiral-type mode, which was not possible via prior axisymmetric studies. These are the subjects of recent investigations.

484

M. Bm.t~R and D. H),s~it REFERENCES

1. A. J. Chorim A numerical method tbr solving incompressible viscous flows. J. Compuc Phys. 2, [2 (1967). 2. C. L. Merkle and M. Athavale, Time-accurate unsteady incompressible flow algorithms based on artificmi compressibility. AIAA Paper 87 1137 (1987). 3. E. Rogers and D. Kwak, Numerical solution of the incompressible Navier-Stokes equations for steady a~,d time-dependem problems. AIAA Paper 89 1463 (1989). 4. M. Breuer and D. H/inel, Solution of the 3-D incompressible Navier Stokes equations tor tile simulation of vortex breakdown. In Proe. 8th GAMM Confi on Numerical Methods in Fluid Mechanics, Delft; Notes on Numerieal Metho& in Fluid Mechanics, Vol. 29. Vieweg Verlag, Braunschweig (1989). 5. P. L. Roe, Approximate Riemann solvers, parameters vectors and difference schemes. J. Comput. Phys. 22, 357 (1981)~ 6. A. Brandt, Guide to multigrid development. In Lecture Notes in Mathematics, Vol. 960, pp. 220 312, Springer-Verlag~ Berlin (1981). 7. B. P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Meth. Appl. Mech. Engng 19, 59 (1979). 8. B. Van Leer, Towards the ultimate conservative difference scheme. A second-order sequel u~ Godunov's method J. Comput. Phys. 32, 101 (1979). 9. M. Breuer, Numerische L6sung der Navier Stokes-Gleichungen f/ir dreidimensionale inkompressible instation/ire Str6mungen zur Simulation des Wirbelaufplatzens. Thesis, RWTH, Aachen (1991). 10. A. Jamesom Solution of the Euter equations for two-dimensional transonic flow by a multigrid method, Appl. Math. Comput. 13, 327 (1983). 11. D. H/inel, M. Meinke and W. Schrader, Application of the multigrid method in solutions of the compressible Navier- Stokes equations. Presented at the 4th Copper Mountain Conf. on Multigrid Methods, Copper Mountain. CO (1989). 12. M. Meinke and D. H/inel, Time accurate multigrid solutions of the Navier Stokes equations. In International Serie~ of Numerical Mathematics, Vol. 98, pp. 289-300. Birkhauser, Basel (1991). 13. J. A. C. Humphrey, A. M K. Taylor and J. H. Whitelaw, Laminar flow in a square duct of strong curvature. J. Huid Mech. 83, 509 (1977). 14. M. Coutanceau and R. Bouard, Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 2, unsteady flow. J. Fluid Mech. 79, 257 (1976). 15. W. M. Collins and S. C. R. Dennis, Flow past an impulsively started circular cylinder. J. Fluid Mech. 60, 105 (1973). 16. M. Meinke and D. H/inel, Simulation of unsteady flows. Presented at the 12th Int. Confi on Numerical Methods in Fluid Mechanics, Oxford (1990). 17. M. Rosenfeld, D. Kwak and M. Vinokur, A solution method for the unsteady and incompressible Navier Stokes equations in generalized coordinate systems. Presented at the A1AA 261h Aerospace Science Mtg, Reno. NV~ A1AA Paper 88 0718 (1988). 18. J. H. Faler and S. Leibovich, Disrupted states of vortex flow and vortex breakdown. Phys. Fluids' 20, 1385 (1977). 19. M. Escudier, Vortex breakdown: observation and explanations. Prog. Aerospace Sci. 25, 189 (1988). 20. S. Leibovich, The structure of vortex breakdown. A. Rev. Fluid Mech. 10, 221 (1978). 21. W. J. Grabowski and S. A. Berger, Solutions of the Navie~Stokes equations for vortex breakdown. J. Fluid ~l/lech. 75, 525 (1976). 22. M. Hafez and J. Ahmad, Vortex breakdown simulation, Part II. AIAA Paper 88-0508 (1988). 23. S. Menne, Vortex breakdown in an axisymrnetric flow. Presented at the AIAA 26th Aerospace Sciences Mtg, Reno, NV (1988). 24. S. Menne, Simulation of vortex breakdown in tubes. Presented at the Ist National Fluid Dynamics Congr., Cincinatti, OH, AIAA Paper 88-3575 (1988). 25. C. H. Liu and S. Menne, Numerical investigations of a three-dimensional vortex breakdown. Presented at the 4th Asian Congr. of Fluid Mechanics, Hong Kong (1989). 26. R. E. Spall and T. B. Gatski, A numerical simulation of vortex breakdown. Presented at the Forum on Unsteady l~Tow Separation, A S M E FlMds Engng Conf. (1987). 27. S. Leibovich and K. Stewartson, A sufficient condition for the instability of columnar vortices. J. Fluid Mech. 126, 335 (1983). 28. E. Krause, The solution to the problem of vortex breakdown. Invited lecture at the 12th Int. ConJ~ on Numerical Methods in Fluid Dynamics, Oxford. In Lecture Notes in Physics, Vol. 371, pp. 35-50. Springer-Verlag, Berlin (1990). 29. S. Menne, Rotationssymrnetrische Wirbel in achsparalleler Str6mung. Thesis, RWTH, Aachen (1986), 30. R. E. Spall, T. B. Gatski and C. E. Grosch, A criterion for vortex breakdown. Phys. Fluids 30, 11 (1987). 31. T. Sarpkaya, On stationary and travelling vortex breakdowns. J. Fluid Mech. 45, 545 (1971).