Numerical study of flow asymmetry and self-sustained jet oscillations in geometrically symmetric cavities

Applied Mathematical Modelling 31 (2007) 2355–2373 www.elsevier.com/locate/apm

Numerical study of flow asymmetry and self-sustained jet oscillations in geometrically symmetric cavities Tomazˇ Kolsˇek a

a,*

, Nikola Jelic´ b, Jozˇe Duhovnik

a

Department of Computer Aided Design LECAD, Faculty of Mechanical Engineering, University of Ljubljana, Asˇkercˇeva 6, 1000 Ljubljana, Slovenia b Institute for Theoretical Physics, University of Innsbruck, Techniker Strasse 25, A-6020 Innsbruck, Austria Received 1 September 2002; received in revised form 1 September 2006; accepted 9 October 2006 Available online 8 December 2006

Abstract In this paper we present the results of numerical investigation of self-sustained oscillations of a jet confined in a symmetric cavity. This work represents an attempt to reproduce empirical observations of asymmetric flows in geometrically symmetric systems and to extend the jet flow investigations to more complex possible scenarios. A well-known example of such two-dimensional flow has been reported experimentally and reproduced numerically for simple flow [E. Schreck, M. Schaefer, Numerical study or bifurcation in three-dimensional sudden channel expansions, Comput. Fluids 29 (2000) 583– 593]. It has been found that for some particular control parameter, above its critical value (bifurcation point), the jet can be deflected to either of the two sides of the cavity. In this paper we report a similar behaviour which is, however, characterized by a more complicated flow pattern. While simple flow appears only within small cavity lengths the complex flows develops with increased cavity lengths. Unlike stationary asymmetric solutions accompanied by cavity jet oscillations, as experimentally reported in e.g., [A. Maurel, P. Ern, B.J.A. Zielinska, J.E. Wesfreid, Experimental study of self-sustained oscillations in a confined jet, Phys Rev. E 54 (1996) 3643–3651], in our investigations of both simple and complex asymmetric flow we observed the slow periodical drift of the jet from one to another side of the cavity. The essential control parameters were Reynolds number Re and the ratio length to inlet width L/d. According to experiments of Maurel et al. (1996), the jet is stable and symmetric, when both L/d and Re are below certain critical values, otherwise jet oscillations appear in both experiment and our simulation (cavity oscillations regime). However, further increase of either (or both) L/d and Re leads of so called free jet type oscillations regime. This paper describes complex jet behaviour within the later oscillations regime. We believe that both simple ‘‘classical’’ and ‘‘our’’ complex stationary asymmetric solutions (as well as superimposed cavity-type and free-jet oscillations) can be explained based on physical arguments as already done in previous works. However, the origin of slow drift motion remains still to be resolved. This might be of high importance for clear distinguishing between relevant physical and numerical features in future codes developments.  2006 Elsevier Inc. All rights reserved. Keywords: Confined jet; Fluid flow; Self-sustained oscillations; CFD

*

Corresponding author. Tel.: +386 1 4771 435; fax: +386 1 4771 156. E-mail address: [email protected] (T. Kolsˇek).

0307-904X/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2006.10.010

2356

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

1. Introduction Self-sustained oscillations are frequently observed in cavity flows and flows in geometries with sudden expansion. They are related to the shear layer instability [1], Coanda effect [2–4], blind cavity effect [5] and crossflow instability [6]. It is important to understand and classify these oscillations, since they may cause either unwanted effects (application in industrial slab casting described by Lawson and Davidson [6]) or even become dangerous from the mechanical point of view (jet energy dissipation chamber described by Jelic´ et al. [7]). On the other hand these oscillations may be usefully exploited, e.g., to serve as a diagnostic tool for determining fluid properties [8]. Oscillations may appear in both laminar [8] and turbulent [9,10] flow regimes and under various geometrical and boundary conditions. Usually cavity oscillations are related to the impingement length, which depends on complex feedback mechanisms and are therefore difficult to predict. Experimental investigations give an essential information concerning flow behaviour and vortices pattern. On the other hand sources of flow behaviour may remain hidden due to technical and physical limitations. For example, long time behaviour in physical systems is very sensitive to initial conditions and may be influenced by small physical perturbations. It is therefore difficult to ensure reproducibility of experimental results. Obviously the behaviour of the system is affected by parameters, whose values range across multiple scales. Although sensitive to initial conditions, numerical simulations can be repeated many times and the results would remain the same for unchanged initial conditions. Long time memory of physical systems can be changed completely or even deleted due to random sources and perturbations, while in numerical experiments all perturbations remain under control even under conditions leading to chaotic system behaviour [11]. Finally, in many experiments it is difficult to avoid use of intrusive diagnostics methods and techniques that essentially influence the flow behaviour, especially in the cases where impingement mechanism is important. Therefore for a complete description and understanding of such systems one can use numerical method, which is less sensitive to randomly generated sources of perturbations when compared to the experiment. In this paper we investigate a geometrically simple but physically still complicated two-dimensional flow that is known for flow asymmetry and self-sustained oscillations. This system has been thoroughly investigated experimentally in the past [8], but detailed numerical investigations have not been reported. Especially the long time behaviour of impinging jets has not been investigated neither in laboratory nor in numerical investigations. As we will demonstrate, the numerical method enables one to observe a large number of solutions like stationary symmetric, stationary asymmetric, non-stationary asymmetric and alternate oscillating modes. In Section 2, we describe the physical system and briefly summarize previously reported results. The numerical method is summarized in Section 3. Our results are presented in Section 4 and finally discussed in Section 5. 2. Problem description We consider the two-dimensional geometry as shown in Fig. 1. Relevant parameters of the system are the same as in laboratory experiments of [8] as follows.

L

d

l1

l2

* * INFLOW

w

ln

*

r1

*

r2 * OUTFLOW

Fig. 1. Geometry of the computational domain.

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373 35

2357

R6 R5

30 Mode 5

R4

25 P4

P1

P2

R2

15

R1

10

P8

free jet type oscillations ty vi

R3

ca

L/d

P3

20

P7

Mode 4 P5 P6

oscilla

tions

Mode 3 Mode 2

Mode 1 stable regime

5 50

100

150

200

250

300

350

400

450

500

Re Fig. 2. The occurrence of the oscillations in the Re–L/d domain according to Maurel et al. [8].

Working fluid in water. The height of the cavity is w = 100 mm and the depth is 26 mm. The length L may take arbitrary values, but the lengths L of the same order as the height w were of interest. Therefore the control parameter ‘‘relative length’’ L/d (according to experimental results by Maurel et al. [8]) has been varied in the range 10 . . . 40 by changing the cavity length L only, whereas the inlet width d = 4 mm was held constant. These parameters correspond to the lengths in the range 40 . . . 160 mm. The height of outflow duct was 20 mm. According to experiments it has been assumed that the outflow duct height is irrelevant (unless using extreme values) for the phenomena of interest. The length of the outflow duct in numerical investigations was 120 mm in order to minimize transversal velocity component at system outlet. The second relevant parameter was Reynolds number Re = Umaxd/m (where Umax is the maximum velocity at the entrance to the cavity and m is the kinematic viscosity). Re has been varied in the range 80 . . . 600. To be able to detect the flow characteristics such as long time periodicity, we have monitored the flow variables in points l1 . . . ln and r1, r2 (see Fig. 1) and recorded them for later interpretation. It has been observed in experiments [8], that self-sustained oscillations appeared at Re > 160, L/d > 7, (Re  100)L/d < 300 (see Fig. 2). For small Re and for small L/d no oscillations were observed. Increasing either Re or L/d (or both) outside the range of self-sustaining oscillations resulted in free-jet type oscillations, whose characteristic length and frequency were dependent on the jet inlet width rather than on the cavity length. The occurrence of the oscillations can be described by the two-dimensional graph in the Re–L/d domain, as depicted by Fig. 2. Inlet velocity was of the order of 0.1 m/s. We estimated that the characteristic time of the physical relaxations processes related to the cavity length should be of an order of several seconds. However, it turned out that it is necessary to observe the system behaviour during much longer period of the order of several minutes. Ref. [8] showed, that the wavelength of the oscillations k is a fraction of the cavity length L. Empirical relation L = (N + )k has been proposed, where N is an integer corresponding to the mode selection and  is the correction factor depending on particular experimental conditions (modes of up to N = 5 have been observed in experiments). 3. Numerical method The flow in the cavity is considered laminar on the basis of estimated Re number at the inlet. It can be fully described by the equations of continuity and momentum conservation. Equations are given for the control volume CV bounded by surface S in the integral form similar for all conserved properties. They all contain local and convective rates of change on the left-hand side and diffusive and source terms on the right-hand side [12]. The continuity equation is Z Z o q dV þ qv ds ¼ 0; ð1Þ ot CV S

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2358

where q is the density and v is the fluid velocity. The momentum equation is Z Z Z Z o qv dV þ qvv ds ¼ T ds þ f b dV ; ot CV S S CV

ð2Þ

where fb is the resultant body force per unit volume and T is the stress tensor _  2 l div vI  pI: T ¼ 2lD 3

ð3Þ

_ is the rate of strain tensor I is the unit tensor, l is the dynamic viscosity, p is the pressure and D _ ¼ 1 ½grad v þ ðgrad vÞT : D 2

ð4Þ

3.1. Approximation of the governing equations A finite volume method was used to solve the above equations. The method is described in detail by Demirdzˇic´ and Muzaferija [11]. Here only the most important features implemented in the method are summarized. The constitutive equations (1) and (2) can be written in the following general form: Z Z Z Z Z o q/ dV þ q/v  ds ¼ C/ grad /  ds þ q/S  ds þ q/V dV : ð5Þ ot CV S S S CV |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Rate of change

Convection

Sources

Diffusion

The computational domain was discretized into control (finite) volumes-CVs (see Fig. 3) of polyhedral shape. For accuracy reasons and because they enable easy local grid refinement, the hexahedral shape was chosen in present investigation. All dependent variables of approximated equations were stored at the CV center (collocated variable arrangement). The following notation is used: P0 is the center of control volume, P1, . . . , Pn are centers of neighbouring control volumes, S1 is face (with normal), which is neighbour to P1, Sj, . . . , Sn are other faces with neighbours j, . . . , n, rP 0 is spatial vector of center of P0, x1, . . . , x3 are Cartesian unity vectors. For a single control volume in Fig. 3 Eq. (5) is as follows: Z Z nf Z nf Z nf Z X X X o q/ dV þ q/v  ds ¼ C/ grad /  ds þ q/S  ds þ q/V dV ; ð6Þ ot CV Sj Sj Sj CV j¼1 j¼1 j¼1 where nf is a number of faces enclosing the control volume. The approximations used for solving Eq. (6) are addressed in the following subsections.

sj

P0

sn Pn

dj

x2 x3

Pj

rP

0

x1

Fig. 3. A general three-dimensional polyhedral control volume and notation.

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2359

3.1.1. Calculation of integrals Surface integrals of convective and diffusive terms in Eq. (6) were approximated by using the value of the term w on the center of CV face by multiplying it with the face area. Volume integrals were similarly approximated by using the value w at the center of the CV multiplied by volume of the CV Z Z w  ds  wj S j ; w  dV  wP 0 V P 0 : ð7Þ CV

Sj

Surface area of faces and volume of CV were calculated using following equations: v

nj 1X sj ¼ ½ðri1  r1 Þ  ðri  r1 Þ; 2 i¼3

V P0 ¼

nf 1X rj  sj ; 3 j¼1

ð8Þ

where nvj is number of nodes of face j and ri is spatial vector of node i, rj is spatial vector of center of face j, sj is normal vector to face j. 3.1.2. Spatial variation To determine a spatial distribution of an arbitrary variable w (dependent variable /, physical properties of fluid or gradient of /), linear distribution between centers of control volumes was used wðr; tÞ ¼ wP 0 ðtÞ þ gðtÞ  ðr  rP 0 Þ;

ð9Þ

where rP 0 is spatial vector of CV center P0 and g(t) is approximation of gradient of w at point P0. Gradient of variable w at point P0 is determined by linearly interpolating w through P0 and neighbouring points. The following system of equations is obtained: dj ðgrad wÞP 0 ¼ wP j  wP 0

ðj ¼ 1; . . . ; nÞ;

ð10Þ

where dj ¼ rP j  rP 0 is vector between P0 and neighbour Pj. The least squares method was used to obtain a single value # " #1 " nf nf X X T T ðgrad wÞP 0 ¼ dj dj  dj ðwP j  wP 0 Þ : ð11Þ j¼1

j¼1

3.1.3. Integration in time The transient term of Eq. (6) has to be discretized. We have used three-time levels implicit discretization scheme Z m1 m2 3ðq/V ÞP 0  4ðq/V ÞP 0 þ ðq/V ÞP 0 o q/ dV  ; ð12Þ ot CV 2dtm where index m denotes time step. The convective, diffusive and source terms were evaluated at the new time level. 3.1.4. Assembly of algebraic equation systems Convective term of Eq. (6) was sum of convective fluxes Cj of variable / (e.g. velocity v) through a face j of CV and represented the rate at which it was transported into (or out of) the CV by the fluid motion relative to CV boundary. This term was nonlinear and needed to be linearized first. The following approximation was used: Z Cj ¼ q/v  ds  qj ðvj  sj Þ/j ; ð13Þ Sj

where * denotes values of variables at face of CV, and index j denotes face j of CV. The method of determining the convected variable /j , density qj and velocity vj had the crucial impact on accuracy and stability of the numerical method. The upwind differencing scheme – UD and central differencing scheme – CD were

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2360

alternatively used. The former is first-order accurate and defines the value of the variable /j as equal to the value of the upstream CV center /U ( /P 0 ; when flow is from P 0 to P j ; UD : ð14Þ /j ¼ /P j ; when flow is from P j to P 0 : The second one (CD) is second-order accurate and represents linear spatial distribution of variable between CV centers. Value on face j is 1 1 /CD ¼ ð/P 0 þ /P j Þ þ ½ðgrad /ÞP 0  ðrj  rP 0 Þ þ ðgrad /ÞP j  ðrj  rP 0 Þ: j 2 2

ð15Þ

CD scheme may cause unbounded solution (unphysical oscillations), which leads to poor convergence, when computational grid is coarse. When grid is sufficiently dense, the solution converges faster than UD towards the grid-independent solution. Diffusive term of Eq. (6) is sum of diffusive fluxes Dj through the face j of CV and was approximated as follows: Z  Dj ¼ l grad /  ds  lj ðgrad /Þj  sj ; ð16Þ Sj

where lj is viscosity at the center of face j, and ðgrad /Þj was determined by using Eq. (11). Source terms include surface and volume integrals. The surface integrals involving the vector q/S are in the following form:   Z nf  X 2 T l div v þ p I  sj lðgrad vÞ  ð17Þ Q/S ¼ q/S  ds  3 j¼1 and volume integrals are integrated using midpoint rule Z Q/V ¼ q/V dV  ðq/V ÞP 0 V P 0 :

ð18Þ

CV

After approximating all terms of Eq. (6) one algebraic equation per control volume was obtained. The equation contained an unknown, which linked the value of an independent variable / in the CV center with the values in the centers of neighbouring CVs a/0 /P 0 

ni X

a/j /P j ¼ b/ ;

ð19Þ

j¼1

where ni is number of CV encompassing P0, and b/ contains source terms, contributions from boundary faces and from transient term, and convective and diffusive fluxes. In order to calculate pressure field and to couple it properly to velocity field, a SIMPLE method was implemented [12,11] which uses two-step algorithm: predictor stage and correction stage. First, all variables were assigned initial values at dt = t0. Time was then advanced to t1 = t0 + dt and an iterative procedure was started to find solution to coupled nonlinear equations at the new time level. The solution of the previous time step served as initial guess for the next solution. To perform these steps, a non-commercial version of computer solver code ICCM COMET [13] was used. More details about the solution procedure may be found in [11,12]. 3.2. Boundary conditions Because of simple cavity geometry we have implemented a block-structured numerical grid of CVs, which was easy to generate with available grid generation tools. We performed numerical investigations in two dimensions only because it has been shown in cited experiments that in the regime of cavity oscillations the flow is two-dimensional. Moreover, we have checked the validity of this assumption by performing fully three-dimensional flow simulation. Fig. 4 illustrates a fair coherence of several streaklines of the simulated

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2361

Fig. 4. Illustration of the results obtained by three-dimensional flow simulation.

jet flow for typical oscillating scenario (Re = 200 and L/d = 40). A more detailed inspection of the flow at various equidistant L–w depth-planes proved the validity of two-dimensional character of the flow. Therefore we proceeded our investigation in two dimensions. It has been found that at least in the case corresponding to the cavity oscillations regime the assumption of planar jet behaviour is well confirmed. In the regime of free-jet oscillations this behaviour is, however, not well confirmed, primarily because of extremely long times needed for obtaining reliable 3D results in this regime. Symmetry boundary conditions have been applied to L–w planes. We have used a specially shaped inflow channel to obtain a fully developed parabolic velocity profile at the inlet to the cavity, as in the case of the experiment (see Fig. 1). At the end of relatively long outflow duct a zero static pressure has been applied. Initial computational with low density grid consisting of 2200 cells is shown in Fig. 5 (no inflow duct is shown for simplicity reasons). We have performed successive grid refinements (up to 11,000 cells) in order to minimize numerical errors. To ensure sufficient accuracy and convergence rate of the calculations we have selected the time step of 104 s. We have typically performed 107 iterations per case to describe the long time behaviour of the order of several minutes. The calculations of cases containing different control parameters (Re and L/d) were performed simultaneously on several PC-s of a medium performance (Pentium 4 with 2.8 GHz processors and 512 MB RAM each). To obtain a coarse picture of phenomena in Re–L/d space one has to monitor the relevant flow quantities at a spatial equidistant grid of at least 100 monitoring points. These cases should later be refined locally along the separatrices (boundaries of coherent oscillations) in the Re–L/d region of interest. In the experiment by Maurel et al. [8], the oscillation wavelength was related to the cavity length L. However, for higher lengths and higher Re the concept of wavelengths can no longer be used because of complicated flow pattern. In the absence of more advanced and reliable tools we recorded animations of particle

Fig. 5. The low density grid, 2200 cells.

2362

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

paths patterns and visually inspected them in order to establish long time periodicity. Our current resources limited the investigation to a couple of hundred cases. Concerning the initial conditions, we found appropriate to use different randomly chosen ‘‘instant’’ flow fields, obtained during the oscillations rather than the ones from steady states. We applied them further as the initial conditions for calculating flow fields under slight change of control parameters (Re and L/d). The reason for applying this procedure was our assumption that slight variations of control parameters should yield the results on the flow field not far from initial if calculated under very similar conditions (e.g., for a little longer cavity). Consequently, such initial conditions might be regarded as a kind of perturbation of the expected final states. Indeed, this assumption was useful in a wide range of control parameters and substantially reduced computational time. Near the separatrices in Re–L/d space this procedure was less efficient and seemed to lead to a kind of hysteresis phenomena (apparent sensitivity to initial conditions) which is deferred to a later study. For post processing and evaluation of the results, we used the flow field in the whole computational domain. Our method of preference for flow visualization was the streakline display because the flow pattern could better be recognized, when compared to velocity vector display. 4. Results We have performed two series of calculations, one with first-order upwind differencing scheme (UD) for convective fluxes and second with second-order central differencing (CD) scheme. UD is known to damp the numerical generated instabilities. If oscillations are detected using UD, they may be regarded as physical. However, UD may damp small amplitude physically present oscillations. Therefore, we have performed also CD calculations. Unfortunately CD may produce some unwanted numerical solution oscillations, which are superimposed on physical ones. By comparing results obtained by both methods and taking experimental results into account, one has to separate physical from numerical oscillations. To confirm the validity of our numerical method, we tried to obtain the results in the form that is most suitable for comparison with experiments by Maurel et al. [8]. For this purpose we selected a set of observation points (see examples of sampling points in Fig. 1). By monitoring simultaneously the velocity magnitude at points r1 and r2, we were able to identify the cavity side where the jet was attached to (resided at) at particular moment. To evaluate the wavelength and frequency of the jet oscillations the results were observed at the symmetry line (l1, . . . , ln according to experiments). Although in these points on-line observation (during the calculation) was enabled as well, most results at the line of symmetry and over the whole two-dimensional calculation domain were calculated in postprocessor. 4.1. First-order discretization We calculated a large number of cases of the flow field in computational domain by varying Re and relative length L/d as the essential parameters. We adjusted the inlet fluid velocity and length of the cavity to a corresponding value as in experiments. According to the definition of Re, these parameters may both alternatively be changed by varying the inlet width d. In order to simplify the calculations in this stage, a spatially flat and temporally constant velocity profile at the inlet has been chosen. We were aware, that such approximation may possibly lead to a shift of effective Re (because of the change in effective width of the jet), but it turned out, that such an assumption did not affect the common features of the flow substantially. Fig. 6 shows the velocity vector field under stationary symmetric conditions obtained for small Re and short cavity (Re = 80, L/d = 15, see P1 in Fig. 2). The same results were obtained when starting from several different initial conditions. The plotted velocity field agrees well with the one observed in experiment. For higher Re (in the range Re = 80 . . . 120), we observed an asymmetric time-independent velocity field. This situation resembles to well-known Coanda effect [2–4]. Superimposed oscillations of small amplitude have been observed, which were damped and finally vanished after sufficiently long time (several minutes). The bifurcation point cannot be precisely correlated to control parameters. We estimate, that this uncertainty is related to the assumption of two-dimensionality of the problem [4]. Illustration of the solution is shown in Fig. 7 for Re = 120 and L/d = 21 (see P3 in Fig. 2). Asymmetric steady state has not been reported in

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2363

Fig. 6. Velocity vector under stationary symmetric conditions obtained for low values of control parameters (Re = 80 and L/d = 15, P1 in Fig. 2).

Fig. 7. Velocity vector under stationary asymmetric conditions for moderate values of control parameters (Re = 120 and L/d = 21, P3 in Fig. 2).

experiments by Maurel et al. [8]. We assume that in experiments the observed phenomena can be detected in much more narrow range of Re, primarily due to random physical perturbations. In the experiment by Maurel et al. [8] self-sustained oscillations appeared above Re  160, independently of the L/d. In numerical simulations, however, the onset of self-sustained oscillations depended on the cavity length and has been observed already above Re = 110 for short cavities (L/d = 15) and became higher as Re increased (see P2 in Fig. 2). In our simulations we observed that the asymmetric steady state flow pattern can slowly drift from one to another side of the chamber. A typical illustrative instant state is shown in Fig. 8, where the results are presented for control parameters Re = 175 and L/d = 21 (see P5 in Fig. 2). The full scenario of the temporal jet behaviour can be seen in animation [14]. In general steady state solutions were always characterized by two large dominant recirculation flow regions located at the opposite sides of the jet. Oscillating solutions were characterized by additional number of smaller vortices. The largest recirculation region was always situated at the side of higher static pressure and smaller dominant one was at the opposite side of the jet. Small vortices interacted with each other and periodically merged. For high jet velocities this main recirculation region partially collapsed due to jet attraction toward the zone of depleted pressure and another additional re-circulation region appeared instead. This situation was similar to, but obviously different from crossflow instability reported by Lawson and Davidson [6].

2364

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

Fig. 8. Example of instant velocity distribution taken while self-sustained oscillations play role (Re = 175 and L/d = 21 P5 in Fig. 2).

The quantitative characterization of non-stationary behaviour was performed by observing the velocity and pressure fields at the symmetry line along different longitudinal distances (i.e., points l1, . . . , ln, see Fig. 1). The results obtained for Re = 140 and L/d = 21 (see P4 in Fig. 2) are shown in Fig. 9, where we illustrate the techniques of determining the wavelength and frequency as used in postprocessing of experimental results by Maurel et al. [8]. The first wave mode is clearly identified. The second mode determination is illustrated in Fig. 10 for Re = 500, L/d = 21, see P8 in Fig. 2. The Reynolds number Re = 500 is quite high but still bellow the onset of turbulence. The quantity (Re  100)L/d for this example is 840, which is far outside the region in which coherent oscillations are reported in experiments (i.e., 300 as we estimated above). However, main features that have been pointed out in experimental investigations are well reproduced: the wavelength of the first mode almost exactly coincides with the cavity length, while the second mode coincides with the half of the cavity length. Reported value of correction factor  = 1/4 should be considered with care. We estimate from both graphical method (see Fig. 9) as well as FFT (see Fig. 14) that the uncertainty of frequency determination is of the order of 20%. Both wavelength and frequency are to some extent variable during both measurements and simulations. On the other hand, value of  = 1/4 suggested by experimenters is of the same order. It was therefore pointless to try to obtain  from our simulations. Moreover, the empirical value 1/4 can be questioned as well. Fig. 11 shows a pair of signals which differ in amplitude and phase because the jet oscillates around permanent non-symmetrical position. Any of these signals may easily be simulated by e.g., plotting a superposition of two sinusoidal signals, such that their frequency differs nearly by a factor of two. This indicates that both first and second modes are both simultaneously present. This situation was not reported in experiments and should be described here. As shown in Fig. 12, there are two neighbouring frequencies that approximately

0.10

absolute velocity [m/s]

0.08 0.06 0.04 0.02 0.00 0

10

20 30 time [s]

40

50

Fig. 9. Example of the velocity amplitude obtained at the symmetry line (points l1, . . . , ln in Fig. 1) and quick estimation of the wave length and frequency for Re = 140 and L/d = 21 (P4 in Fig. 2).

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2365

100

distance from inlet [mm]

90 80 70 60 50 40

λ

30 20 10

T

0

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 time [s]

Fig. 10. Example of second unstable mode obtained for high Reynolds number Re = 500 and L/d = 21 (P8 in Fig. 2).

absolute velocity [m/s]

0.035 0.030 0.025 0.020

r1 r2

0.015 0.010 100

150

200

250

300

Time [s]

Fig. 11. Typical signals (velocity magnitude) observed simultaneously in monitoring points r1 and r2 for the case of moderate Reynolds number Re = 225 and L/d = 21 (P6 in Fig. 2), when first and second modes are both present.

0.005

Amplitude [m/s]

0.004 0.003 0.002 0.001 0.000 0.0

0.1

0.2 0.3 0.4 frequency [Hz]

0.5

0.6

0.7

Fig. 12. Amplitude of Fourier transform of the signal obtained in r1 from Fig. 11.

correspond to first mode ( 0.1562 Hz and 0.1367 Hz). By using these frequencies we obtained for the first mode k = 0.105 m and k = 0.120 m which were substantially above the cavity length L = 0.84 m. Higher harmonics of both of these frequencies were present as well. Finally, Fig. 13 shows an interesting situation, where an ‘‘alternately asymmetric jet’’ behaviour was observed. The system behaviour can be described by the following characteristic quantities:

2366

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

absolute velocity [m/s]

0.07 0.06

r1

0.05

r2

0.04 0.03 0.02 0.01 0.00 0

20

40

60

80 100 120 140 160 180 Time [s]

Fig. 13. Velocity magnitude oscillations in points r1 and r2 obtained for Re = 350 and L/d = 25 (P7 in Fig. 2).

• characteristic period of jet oscillations (approx. 2 s in this particular case); • characteristic time of flight of the fluid particle (measured from the entrance into the chamber to the entrance of outlet duct, for stable regime and cavity type oscillations is approx. 1 s and for long term oscillations 2–4 s); • characteristic drift period Tdrift, when the jet slowly moves from one to another quasi-stable asymmetric state (50 s); • characteristic period of long term oscillations Tlong (in our case approx. 3 min). The jet was oscillating around one or another side of the outlet edges. This situation is clearly illustrated in Fig. 14, where we present velocity magnitude oscillations observed at geometrically symmetric points r1 and r2. The results were obtained for Re = 350 and L/d = 25 (P7 in Fig. 2). The curves, showing the velocity magnitude as measured in these points during the long time, have the same shape but seem to be shifted one with respect to another for  93 s. The time period of the long term behaviour (Tlong) was thus in our case 186 s, which is much longer than any expected characteristic time of the system. This finding initiated investigations of jet behaviour in modified geometries. We performed simulations of geometries with two outlets positioned symmetrically at the inlet cavity side. This numerical domain with streaklines behaviour is illustrated in Fig. 15. The full scenario of temporal jet behaviour can be found at [15]. Initial results confirm long time periodicity of slow jet drifting from one side of the cavity to another. The streaklines shape repeated alternately in mirrored form after a well-defined time Tlong/2. In more detail, the scenario of alternately asymmetric jet behaviour is as follows:

0.010 0.009 0.008 amplitude

0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 frequency [Hz]

Fig. 14. Amplitude of Fourier transform of the signal, depicted in Fig. 13, monitored at r1.

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2367

Fig. 15. Jet and vortices sreaklines under non-stationary conditions.

• first there was a short stage appearing every Tlong/2 s (cross-points of the solid and dotted lines in Fig. 13, i.e., t = 0, t = Tlong/2, t = Tlong, . . .) in which the system was nearly symmetric and velocity field was steady, like the one shown in Fig. 6; • in the second stage the jet started to drift on either one of the two symmetric sides. The direction of the drift depended on the previous history of the jet behaviour. This stage lasted about 25 s and the jet looked similar to an asymmetric steady state (illustrated in Fig. 7). Instead of oscillations a slow drift of the jet towards one of the sides of the cavity prevailed; • during the next 25 s the oscillations (of the frequency 0.17 Hz) took place, which were characterized by small amplitude and single asymmetric re-circulation regions on both sides of the jet; • the final stage of a half-period was characterized by high frequency oscillations of 0.45 Hz and lasted up to the end of the half period. These oscillations appeared suddenly and were of large amplitude. They were gradually decaying in time while at the same time the jet slowly drifted back to the line of symmetry of the cavity. The number of re-circulation regions at the low-pressure side of the cavity dynamically changed during each oscillation period. The re-circulation vortices appeared, drifted, interacted with each other, and finally merged or collapsed. After that the same scenario repeated towards the opposite side of the cavity. Above features confirm fair qualitative reproducibility of experimental results. We performed investigations of this kind for the relative lengths L/d = 15, 21 and 25 changing carefully Re in small steps for each length and using previous flow field state as the initial condition for next. As expected, the first-order discretizations gave quantitative results, which differ from experimental ones (sometimes up to about 20%). However, we estimate that the error of experimental results is approximately the same. Moreover, we were not able to obtain fine jet structures related to free jet instability with these discretizations. Finally, the coherence of oscillations remains present even for high Re (590). Refinement of the grid in the range from 2000 to 11,000 elements did not reveal other phenomena. 4.2. Second-order discretization In order to obtain more reliable results we proceeded with calculations by applying second-order accurate central differencing (CD) scheme for convective fluxes. Using CD approximation for viscous terms yields fine

2368

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

jet oscillatory structures for both fine and coarse grids [12,16]. Henceforth CD turned out to be essential for obtaining the small-scale flow patterns, for example free-jet oscillations. Our experience shows, that • when observing the velocity magnitude signal at monitoring points r1, r2 and ln, calculated using CD, the signal became enriched by apparently random high frequencies, whereas the wavelength and main frequency remained nearly the same as obtained by applying UD, • cavity oscillations remain coherent for high Re and long cavity lengths. First-order discretizations were very useful because they served us as a filter for main physical phenomena. As demonstrated bellow, it would be a very difficult task to observe these phenomena by performing CD calculations only. The time step had to be decreased for an order of magnitude in order for solution to converge, which increased considerably the total calculation time. As a result, we needed at least a whole week calculation time for obtaining a solution, that was simulating several real minutes of fluid flow. Due to demanding computer resources we decided to proceed concentrating on an unambiguous case (according to experiments) of Re = 200, changing only the relative length of the cavity L/d. The values of 10, 15, 20, 25, 30, 35 (points

Fig. 16. Velocity magnitude and pictures of the jet profile in some characteristic times for Re = 200 and L/d = 10 (R1 in Fig. 2).

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2369

R1, . . . , R6 in Fig. 2) were investigated in detail and complete flow fields have been sampled and recorded at time intervals of 0.1 s. The main features of the flow in first case (L/d = 10, R1 in Fig. 2) are illustrated in Fig. 16. We show the time series of the signals observed at monitoring points r1 and r2. Using UD, the steady state or period of oscillations was easy to identify. However, using CD calculations the results were not easy to interpret, since the signals appeared to be stochastic and the information taken in particular locations ln and r1, r2 was not sufficient. In fact, by filtering signals one might get similar periodic signals to UD ones, but to prove the similarity the analysis had to be performed a posteriori and not during the calculations, since one had to know the range of filter pass. To visualize the jet behaviour we have chosen several instant flow fields and illustrated them in Fig. 16. The particle paths are shown at characteristic times during a cycle of the jet movement in the system. The figure describes the phases of the jet behaviour as follows: • the jet is steady for a relatively long time (see Fig. 16, t = 30 s) and attached to bottom edge of the outlet channel; • in the following 10 s weak oscillations of the jet developed with increasing amplitude. The mean jet direction is still steady; • during following 10 s (t = 40 . . . 50 s) the amplitude of the oscillations was increasing and the jet started to drift from bottom side of the outlet tube (see time t = 54.2 s); • as the jet passed the symmetry line of the cavity, the amplitude of oscillations grew (see time t = 57.3 s);

0.05

deriv of vel. magn.

0.04 0.03 0.02 0.01 0.00 94

96

98

100 102 time [s]

104

106

108

Fig. 17. The derivative of the measured velocity magnitude along the symmetry line for Re = 60 and L/d = 10 (R1 in Fig. 2).

0.07

deriv. of vel. magnit.

0.06 0.05 0.04 0.03 0.02 0.01 0.00 0

2

4

6 time [s]

8

10

12

Fig. 18. The derivative of the measured velocity magnitude along the symmetry line for Re = 200 and L/d = 15 (R2 in Fig. 2).

2370

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

• after jet has reached the upper side of outlet tube the oscillations decayed (see time t = 70.6 s); • during next 10 s (t = 70 . . . 80 s) the jet was steady and attached to the upper side of outlet tube (see time t = 80.8 s); • a similar scenario was repeated in reverse direction over the line of symmetry. Described flow behaviour is demonstrated by jet streaklines animation in [17] and by both jet and vortices streaklines in animation [18]. Careful inspection showed, that any instant state had its ‘flipped’ counterpart, hence the counterpart of streaklines at t = 53.7 s were streaklines at t = 107.9 s, counterpart of t = 70.6 s was likely to be t = 128.5 s, etc. The similarity should be regarded as a qualitative one and has been detected by observing the velocity magnitude diagram at points r1 and r2. After appropriate time shift the velocity diagram shapes fit fairly, but not exactly (as in the case with UD calculation). Similar behaviour was observed for L/d = 15 (R2 in Fig. 2). Up to the length L/d = 20 (R3 in Fig. 2) the wave-like behaviour might be identified and analysed. The long time behaviour can also be identified once the long term period has been found. We used the derivative of the velocity amplitude signal along the line of symmetry. Fig. 17 shows the signal for the case L/d = 10. The zeros of derivatives were much more reliable for determining the wavelength than the

Fig. 19. Velocity magnitude and pictures of the jet profile in some characteristic times for Re = 200 and L/d = 25 (R4 in Fig. 2).

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2371

original signal. In addition, we identified the main frequencies by Fast Fourier Transformation (FFT) analysis. The spatial and temporal acquisition of the velocity gradients yields an appropriate diagnostics of wave phenomena. We showed the most ambiguous time sequence between in the range 94 . . . 108 s, where oscillations behave stochastically. We managed to identify a wavelength, which corresponds to second unstable mode. It was obvious, that this wavelength was just embedded in a pool of higher and lower harmonics and that incommensurate frequencies were also present. On the other hand, there was a wave-like behaviour, which corresponded to first unstable mode during time sequence in the range 80 . . . 94 s. However, it seemed that during the long time observations the frequency and length of wave-like structures changed in time as the jet moved from one to another side of the outlet edges. The length L/d = 10 was very close to the onset of oscillations. This could be the reason for uncertainty of oscillation occurrence of the jet. For the lengths L/d = {15, 20, 25, 30, 40} new flow patterns appeared. For the length L/d = 15 (R2 in Fig. 2) the jet steady state was not detected. First and second unstable modes were easily identified as shown in Fig. 18. For higher cavity lengths the concept of waves cannot be applied. As we showed in Fig. 19 already for L/d = 25 (R4 in Fig. 2), the jet streaklines became complicated. Long time periodic jet behaviour and similar ‘flipped’ trajectories as for smaller lengths could still be detected. Harmonic oscillation spectra still occurred and could be identified. In addition, we observed that the velocity along jet streakline might substantially increase and become higher than the initial velocity. This might be related to the local channel width, determined by the vortices that occupy the most of the cavity space. Sometimes the jet bifurcated to two or even several branches. In general, the jet behaviour observation may be summarized as follows: • during the time which was much longer than any of expected characteristic time scales of the system the jet quasi-periodicaly followed a pattern, which became more complicated as the cavity length was increased;

Fig. 20. Velocity magnitude, pictures of the jet profiles and complete streakline pictures in two characteristic times for Re = 200 and L/d = 40.

2372

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

• both free jet instabilities and large amplitude drift of the jet i.e., deviations from a straightforward trajectory near the symmetry of the cavity appeared already for L/d = 21. As in experiments it became even more difficult to identify oscillations explicitly related to the cavity length above L/d = 35 (R6 in Fig. 2). Small scale streakline patterns along the jet trajectories remain present for long cavity lengths. The system is very rich in coherent structures and the only news is that the concept of oscillations may not be used further. In Fig. 20, we show the time series of the signal obtained in points r1 and r2, the streaklines of the jet (small pictures) and the streaklines over the whole region of the cavity in two characteristic instances, when two apparently ‘flipped’ jet trajectories were observed. This illustrates the fact that coherent structures and long time behaviour remain present for long cavity lengths. In general, there was a finite number of large vortices which were always present and a finite number of small vortices which were born and merged constantly. The behaviour of the vortices was obviously maintained by the jet kinetic energy, which steadily ‘pumped’ the energy to them, a process driven by viscosity forces. The shape of the jet trajectory and its local velocity depended on the instant distribution of the vortices. When there was enough space the jet became wide and its velocity was low, otherwise it had to go through narrow spaces between the vortices and its velocity there increased substantially above the initial jet velocity (even for several times). Some aspects of the physical mechanism of the oscillating vortex triggering and maintaining are discussed by Lawson and Davidson [6] and Gebert et al. [10]. 5. Conclusion We describe the results of investigations of the confined jet asymmetry and self-sustained oscillations during the time, which is much longer than the expected characteristic time of the system. Our results were obtained by CFD using both first and second-order differencing schemes and fit well the experimental ones reported in literature. We observed new physical phenomena at the large time-scale, such as alternating asymmetry jet flow behaviour followed by a variety of oscillations of different time and length subscale. The question to be resolved is whether this behaviour is caused by e.g., accumulation of numerical errors. In any case, this is definitely an important issue for both experimental and numerical investigations, which could help further development of the reliable and versatile codes. The coherence of long time processes and reproducibility of the results indicate that the wavelength and frequency are precisely determined by the control parameters only in a certain range. This range was reported in experiments and its validity was checked in our investigations. We have extended the range of control parameters and obtained new jet behaviour, where the concept of wavelengths cannot be used any more. Linear analysis is not an appropriate tool for describing such systems. We need new methods for quantitative description of such spatial-temporal phenomena. The first step towards an appropriate characterization of such systems is to investigate possible routs to chaos [19]. First-order discretization calculations seem to be useful for checking the basic asymmetry and cavity oscillations phenomena as well as for identifying long time behaviour. Without using it, it would be very difficult to detect long time behaviour because of additional frequencies appearing in calculations, when using secondorder discretization calculations. For more precise quantitative results and to obtain free-jet oscillations at all, second-order discretization calculations have to be used. From the technical point of view it is important to be able to influence the jet behaviour by involving some external driving forces. We showed, that random forces (whether generated physically or numerically) may influence the global system behaviour. The future work is to define a set of external forces or some synchronization, which can control the system behaviour. The solution involving external forces might be crucial for solving some scientific and industrial problems. In addition, experiments should be performed for modified geometries and in a broader range of control parameters. The two-dimensionality of the problem should also be investigated in more detail. We estimate that more realistic boundary conditions (for example using the wall b.c. instead of symmetry b.c.) would lead to discretization of the oscillation spectrum obtained in numerical experiments. Transient three-dimensional analysis for described grid density level requires a lot of computer resources, which are not available at the moment.

T. Kolsˇek et al. / Applied Mathematical Modelling 31 (2007) 2355–2373

2373

Acknowledgements These investigations were supported by the Ministry of Education, Sport and Science of the Republic of Slovenia (Contr. No. 0782-19 P0 0512). The authors thank Institute of Computational Continuum Mechanics GmbH, Hamburg, Germany for providing the computer code Comet that was used to perform the numerical calculations. We also are indebted to Prof. M. Peric for his continuous help in proper set-up of simulation domain and code discretization methods. Finally, we are obliged to Prof B. Sirok for many discussions on numerical method in chaotic systems. References [1] D. Rockwell, E. Naudasher, Self sustained oscillations of impinging free shear layers, Ann. Rep. Fluid Mech. 11 (1979) 67–94. [2] W. Cherdron, F. Durst, J.H. Whitelaw, Asymmetric flows and instabilities in symmetric ducts with sudden expansion, J. Fluid Mech. 84 (1978) 13–31. [3] M. Shapira, D. Degani, D. Weihs, Stability and existence of multiple solutions for viscous flow in suddenly enlarged channels, Comput. Fluids 18 (1990) 239–258. [4] E. Schreck, M. Schaefer, Numerical study or bifurcation in three-dimensional sudden channel expansions, Comput. Fluids 29 (2000) 583–593. [5] N.A. Molloy, Oscillatory flow of a jet into a blind cavity, Nature 224 (1969) 1192–1194. [6] N.J. Lawson, M.R. Davidson, Crossflow characteristics of an oscillating jet in a thin slab casting mould, Trans. ASME 121 (1999) 588–595. [7] N. Jelic´, T. Kolsˇek, J. Duhovnik, A. Bergant, Dissipation in a vertical needle valve induced jet in a pressure chamber, J. Mech. Eng. 46 (2000) 595–605. [8] A. Maurel, P. Ern, B.J.A. Zielinska, J.E. Wesfreid, Experimental study of self-sustained oscillations in a confined jet, Phys. Rev. E 54 (1996) 3643–3651. [9] B. Guo, A.G. Langrish, D.F. Fletcher, An assessment of turbulence models applied to the simulation of a two dimensional submerged jet, Appl. Math. Modell. 25 (2001) 635–653. [10] B.M. Gebert, M.R. Davidson, M.J. Rudman, Computed oscillations of a confined submerged liquid jet, Appl. Math. Modell. 22 (1998) 843–850. [11] I. Demirdzˇic´, S. Muzaferija, Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology, Comput. Methods Appl. Mech. Eng. 125 (1995) 235–255. [12] H. Ferziger, M. Peric´, Computational Methods for Fluid Dynamics, Springer Verlag, 1999. [13] Institute for Computational Continuum Mechanics GmbH, Comet 2.0 User manual, 2001. [14] LECAD, Faculty of Mechanical Engineering, University of Ljubljana, Animation http://www.lecad.uni-lj.si/research/supplements/ AMM06/drifting-jet.mpeg, 2006. [15] LECAD, Faculty of Mechanical Engineering, University of Ljubljana, Animation http://www.lecad.uni-lj.si/research/supplements/ AMM06/modified-chamber, m1v, 2006. [16] V. Seidl, S. Muzaferija, M. Peric´, Parallel DNS with local grid refinement, Appl. Sci. Res. 59 (1998) 379–394. [17] LECAD, Faculty of Mechanical Engineering, University of Ljubljana, Animation http://www.lecad.uni-lj.si/research/supplements/ AMM06/complex-jet.mpeg, 2006. [18] LECAD, Faculty of Mechanical Engineering, University of Ljubljana, Animation http://www.lecad.uni-lj.si/research/supplements/ AMM06/complex-full.mpeg, 2006. [19] H.G. Schuster, Deterministic Chaos: An Introduction, Physik-Verlag, 1984.