Numerical investigation of a confined jet flow in a rectangular cavity

We present the results of investigation of a low-Reynolds number uncompressible fluid flow in a rectangular two-dimensional cavity. It has been shown that the behavior of the flow may be controlled by the cavity length and width. Depending on the ratio between the length and width of the chamber the flow is either steady or non-steady. It has been found that for small chamber lengths the flow is steady and geometrically symmetric. Increased chamber lengths lead to unsteady flows with oscillating and hardly predictable jet flow pattern. Further increase of the chamber length leads again to steady but asymmetric flow. These observations are useful for possible control of the jet behavior under various conditions. This research is relevant for both basic jet-flow investigations as well as for industrial purposes. The sensitivity of the solution to perturbation of simulation parameters is identified as an important issue to be attacked in future by both experimental and numerical means.


INTRODUCTION
The jet flow confined in a cavity is a subset of numerous fluid-flow investigations because of its importance in a wide range of applications, as well as due to its relevance for understanding basics of jet behavior physics.The jet controlling is of high importance in e.g., casting processes [1], chemical engineering [2,3,4], nuclear sciences [5], as a diagnostic tool for determining fluid properties, (e.g., [6] and references given therein), etc.An inherent property of the jet confined in a cavity is self-sustained oscillations due to different kind of impingements.Oscillations appear in both laminar [6,7] and turbulent [8,9] flow regimes and under various geometrical and boundary conditions [10,11].In addition, the jet in a rectangular cavity under different conditions can be used as an excellent reference for testing and comparison of available experimental, theoretical and numerical methods and tools.
Investigation presented in this manuscript was motivated by very successful numerical study of Jelic´ et al [8] of a water jet energy dissipation chamber with a needle valve.Numerical simulations were performed in a full 3D geometry as shown in Fig. 1.Both the numerical grid and jet flow streak-lines snapshot during the oscillations of the jet in the chamber are presented.Simulations were performed for various chamber lengths and needle valve openings for optimization of the dissipation regime.The aim was to avoid possible cavitation of the valve and to minimize the damages of the chamber due to mechanical oscillations caused by the water jet.Qualitative behavior and quantitative results agreed excellently with experimental results performed in a test chamber.The computational task was performed within several days on a cluster of eight Pentium 4 processors.
Figure 1 The numerical grid and jet flow streak-lines as obtained in study of Jelic´et al. [8] These simulations encouraged Kols ˘ek et al [7] to start benchmarking investigations on an apparently simple geometry with apparently less demanding physical scenario of low Reynolds number flow following the experimental results of Maurel et al [6].
Figure 2 Streaklines of the jet as followed from the inlet to the outlet of the chamber with experimental parameters as in the study of Maurel et al. [6], in the case of periodic oscillations Fig. 2 shows the streak-lines of the jet as followed from the inlet to the outlet in the chamber with experimental parameters as in [6].In this simple geometry the inlet and outlet are at the opposite sides of a rectangular chamber.In addition to experimentally observed small amplitude oscillations that are well related to the chamber length and appear for both moderate Reynolds numbers and chamber lengths, in numerical simulationsKols ˘ek et al [7] found a slow periodical drift of the jet from one to another side of the chamber.The jet followed a complicated pattern in space and time.Full 3D simulations were performed with the chamber depth of 20 mm as in laboratory experiments.The results illustrate that the streak-lines in 3D are well aligned with each other over the depth of the chamber, confirming that the assumption of 2D flow is valid.This is a very important fact because 3D simulations were very time-consuming even with eight parallel processors.Therefore most of subsequent simulations were performed in 2D geometry.
The main control parameters in experiments of Maurel et al were the Reynolds number and the chamber length.According to these experiments the chamber width and depth were held constant in numerical simulations as well, while the chamber length and inflow velocity were varied in small steps.Until both of these parameters are moderate the jet flow path is rather simple as illustrated in Fig. 2. Increase in the parameters result in triggering of complicated jet motion, which is self-maintained (see [7] for details).Understanding and explanation of such a motion is additionally complicated due to the presence of sharp edges at the outlet, so eliminating the edge effect could help to better understand the complicated jet flow pattern.Consequently we created new computational domain which is the subject of the present document.In the present research the outlet is at the same side of the chamber as the inlet.Therefore any triggering of the oscillations due to sharp edges is expected to be far downstreams.It turned out that under such conditions harmonic oscillations related to the chamber length disappear.Instead a self interaction of the jet was observed which is related to both length and width of the chamber.While in our previous investigations we ignored the chamber width as suggested by Maurel et al [6], in current research we have taken into account the effects of the chamber width as well.It has been found that this parameter is as important as the chamber length for the oscilation phenomenon we investigate.We observed several solutions like stationary symmetric, stationary asymmetric, non-stationary asymmetric and alternate oscillating modes.Finally, we have found that while small the perturbations of the grid do not influence the global jet behavior, the detailed history of its development is very sensitive on them.As a result, this apparently "over-simplified" flow system appears to be a complicated one with a dynamics that could be possibly resolved in future only by using the chaotic behavior analysis method.
In Section 2 we describe computational domain and parameters.The numerical method is presented in section 3. The results are presented in Section 4. In Section 5 we discuss the results.

PROBLEM SET-UP
We consider the two-dimensional geometry as shown in Fig. 4. Relevant parameters of the system are the same as in laboratory experiments of [6] as follows.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 [6]) has been varied in the range 10...40 (by changing the cavity length L only), whereas the inlet width d = 4mm was held constant.These parameters correspond to the lengths in the range 40 mm ...160 mm.The width of the outflow duct was 10 mm and the length 200 mm.The second relevant parameter was Reynolds number Re = U max d/v (where U max is the maximum velocity at the entrance to the cavity and v is the kinematic viscosity).Re was held constant at Re = 200.To be able to detect the flow characteristics such as long time periodicity we have monitored the flow variables in several points.For post processing and evaluation of the results we used the flow field in the complete computational domain.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. 2 illustrates a fair coherence of several streaklines of the simulated 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 behavior is 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 [6].At the end of a relatively long outflow duct a zero static pressure has been applied.

NUMERICAL METHOD
A finite volume method was used to solve the time dependent uncompressible Navier-Stokes equations.We have used ICCM COMET CFD non-commercial flow solver code to perform numerical calculations.Second-order accurate discretization schemes were used for approximation of convective terms of momentum equations.The flow in the cavity is considered laminar on the basis of the estimated Re number at the inlet.It can be fully described by the equations of continuity and momentum conservation.In the following, the numerical method implemented in COMET is described.Equations are given for the control volume CV bounded by surfaces in the integral form simmilar for all conserved properties.The continuity and momentum equation, respectively, are -the rate of strain tensor.The constitutive equations ( 1), ( 2) can be writen in the following general form: (3) φ represents general conserved quantity (velocity, for example).The interpretation of the above generalization of the eq.( 2) is as follows: Γ is the diffusion coefficient and equals dynamic viscosity µ, q φV equals body forces per unit volume, and q φS equals µ(gradv) T - 2  3 µdiv vI -pI.A Finite volume method was used to solve the above equations.The computational domain was discretized into control (finite) volumes-CVs of hexahedral shape (see Fig 5).where n f is a number of faces enclosing the control volume.The following notation is used in the text: P 0 is the center of control volume, P 1 ...P n are centers of neighbouring control volumes, S 1 is face (with normal), which is neighbour to P 1 , S j ..S n are other faces with neighbours j...n, r P0 is spatial vector of center of P 0 , x 1 ...x 3 are Cartesian unity vectors.
Surface integrals of convective and diffusive terms were approximated by using the value of the term ψ on the center of CV face by multiplying it with the face area.Volume integrals were aproximated by using the value ψ at the center of the CV multiplied by volume of the CV: (5) Surface area of faces and volume of CV were calculated using following equations: (6) where n v j is number of nodes of face j and r i is spatial vector of node i, r j is spatial vector of center of face j, s j is normal vector to face j.To determine a spatial distribution of an arbitrary variable ψ linear distribution between centers of control volumes was used: (7) where r P 0 is spatial vector of CV center P 0 and g(t) is approximation of gradient of ψ at point P 0 .Gradient of variable ψ at point P 0 is determined by linearly interpolating through and neighboring points.The following system of equations is obtained: (8) where d j = r P j -r P 0 is vector between P 0 and neighbour P j .The least squares method was used to obtain a single value: The transient term is approximated via three-time levels implicit discretization scheme: (10) where index m denotes time step.Approximation: (11) Is used for the convective term.Furthermore, both upwind differencing scheme -UD and central differencing scheme -CD were alternatively used.Diffusive term, D j = ∫ S j µgradφ .ds ≈ µ j (gradφ) j * .s j where µ j is viscosity at the center of face j, is sum of diffusive fluxes D j through the face j of CV.Source terms include surface integrals (12) and volume integrals integrated are using midpoint rule: After approximating all terms one algebraic equation per control volume was obtained containing one unknown, which linked the value of an independent variable in the CV center with the values in the centers of neighbouring CVs: (14) where n i is number of CV encompassing P 0 , 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], [13] which uses two-step algorithm: predictor stage and correction stage.To perform these steps, a non-commercial version of computer solver code ICCM COMET [14] was used.More details about the solution procedure may be found in [12,13,14].
The above procedure is valid for both 2D and 3D space.However, since COMET considers general 3D problems, the problem dimensionality has been reduced by the implementation of symmetry boundary conditions at selected faces.
A computational grid consisting of approximately 20000 cells has been used.The time step of the numerical procedure has been selected according to the CFL (Courant, Friedrichs, and Lewy) condition.It requires that the numerical wave speed ∆x ⁄∆t be at least as fast as the physical wave speed c 1 ⁄ 2 for numerical procedure to be stable.Equivalently, the time step ∆t has to be smaller than ∆x ⁄c 1 ⁄ 2 .The maximum media speed within the cavity occured at the inlet and was of the order 0.1 m/s, whereas the distance between numerical grid points within the cavity was of the order of 2 mm.The CFL condition is fullfilled, if the time step is lower than 0.006 s.To ensure sufficient accuracy and convergence rate of the calcluations we have selected the time step of 10 -4 s.We have typically performed 10 7 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 512MB 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 localy along the separatrices (boundaries of coherent oscilations) in the Re -L/d region of interest.
In the experiment [6], 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 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 innitial 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 innitial 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 innitial 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.

RESULTS
Reynolds number Re = 200 was held constant.The attention has been paid to chamber length parameter influence to the flow regime.It has been shown that bellow L = 85mm temporaly steady solutions developed.Such a stationary solution for L = 80mm is shown in Fig. 6.However, further increase of the chamber length above L = 85 mm leads to a quasi-stationary asymmetric solution.High frequency free-jet oscillations still can still be recognized in the animation [16].As a result we found just a very narrow region of lengths of the order of 1 mm (one percent of the chamber length) in which the jet drifts from one to another side of the chamber.This situation is different from the one observed in [7], where such a drift remains inherent above certain critical chamber lengths.However, unlike investigations of [7] in the present case the Re is held at a fixed value.It can be expected that with increased Re new narrow range of length will be characterized by drifting jet.Unfortunately, with the present computational resources it would be a very time consuming task to find this range for sufficiently high number of Re values from which a generalization of the Re-L/d diagram could be assembled, as performed by [7].
In Fig. 9 we show the quasi-stationary jet behaviour (macroscopicaly stable pattern without the jet drift with superimposed wave-like structures along the jet) which is established after a small increase of the chamber length to 90 mm.Unlike non-stationary behaviour of jet shown in Fig. 7, where the direction of the jet drift is trigered each time when either of the two jet branches is impigning the jet near the inlet, in the stationary cases backward branches never touch the jet.The quasi-stationary solution above critical chamber length can be obtained at either side of the chamber, depending on the history, i.e., initial conditions, slight changes in the chamber width or the grid local density refinements.However, above the chamber length L = 100 mm, which approximately equals the chamber width, any further increase of the chamber length L has no effect on the jet pattern.There is definitely a central vortex on which the jet is settled.This situation is illustrated in Fig. 10.Based on above observations we have concluded that the diameter of the vortices might be as relevant parameter as than the chamber length concerning the stationarity of the jet.From the animation it can be seen that the self-interaction of the jet where the velocities of incoming and outcoming branches are in the opposite directions play a crucial role in instability develoment.We supposed that the width of the chamber might be a controll parameter for this process.Indeed, it has been found that slight increasing of the chamber width might stabilize the oscillations.This observation is documented in Fig. 11, where we obtained quite different solutions for the same initial conditions by slightly varying the chamber legth for 2%.For w = 100 mm we obtained a "chaotic" behaviour, while a steady solution developed by a slight increase of the chamber width.Finaly, we slightly perturbed the grid inside numerical domain without changing the boundary and initial conditions.It turned out that the history of the jet behaviour in unsteady regime is dependent on the grid structure, whereas under stationary conditions the solution is grid independent.This feature is illustrated in Fig. 12 where we present solutions for both unperturbed and perturbed grids.In both cases the global behaviour is the same whereas the detailed history is different.

CONCLUSION
Several flow patterns of a laminar jet confined in a rectangular cavity have been observed.It has been shown that the behavior of the jet might be controlled by both longitudinal and transversal cavity dimensions.In addition it has been found that the jet behavior during its oscillatory and non-stationary regimes can not be described by classical mathematical methods but resemble chaotic systems (see e.g.[17].Various solutions found by perturbing the external control parameters resulted in simmilar solutions, which differ in detailed temporal development and the dominating asymmetry direction.The reason for this is inherent for nonlinear systems where bifurcation of the solution is known for many years.However, this behavior should be carefully investigated in experiments.The origin of the assymetric drift of the jet from one to another side of the chamber has to be further investigated.It may be attributed to the physics or the numerical error accumulation.
The authors thank Institute of Computational Continuum Mechanics GmbH, Hamburg, Germany for providing the computer code Comet 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.

Figure 3 :
Figure 3: Streaklines of both jet and vortices streaklines in the case of free-jet flow (higher Reynolds number or/and longer chamber)

Figure 4
Figure4The geometry of the computational domain

( 2 )
With ρthe density, v-the fluid velocity, f b -the resultant body force per unit volume, T = 2µD .-2 -2 3 µdiv vI -pI the stress tensor, I -the unit tensor, µ -the dynamic viscosity, p -the pressure and D .

Figure 5 A
Figure 5 A hexahedral control volume and notation

Figure 6
Figure 6 Streaklines under stationary conditions obtained for L = 80 mm

Figure 7
Figure 7 Jet and vortices streak-lines under non-stationary conditions for L = 86 mm

Figure 8
Figure 8 Velocity magnitude at a central point after oscillations of the jet have developed

Figure 9
Figure 9 Jet streak-lines (quasy-stationary) for an increased length L = 90 mm, no jet oscillations

Figure 10
Figure 10 Quasi-stationary jet streak-lines for a very long chamber

Figure 11
Figure 11 Velocity magnitude at a central point after oscillations of the jet have developed

Figure 12
Figure 12 Velocity magnitude for perturbed and non-perturbed grid