PHYSICS OF FLUIDS 21, 065110 共2009兲 Macroscopic analysis of a turbulent round liquid jet impinging on an air/water interface in a confined medium Jerome Larocque,1,2,a兲 Nicolas Rivière,1,2 Stéphane Vincent,2 David Reungoat,2 Jean-Pierre Faure,1 Jean-Philippe Heliot,1 Jean-Paul Caltagirone,2 and Mathieu Moreau2,3 1 CEA-CESTA, 33114 Le Barp, France Transfert, Ecoulements, Fluides, Energétique TREFLE, UMR CNRS 8508, Site ENSCPB, Université Bordeaux 1, 16 Avenue Pey-Berland, 33607 Pessac Cedex, France 3 Institut Jean Le Rond d’Alembert, UMR CNRS 7190, 4 Place Jussieu, 75252 Paris Cedex 05, France 2 共Received 30 January 2008; accepted 6 January 2009; published online 22 June 2009兲 In this paper, a submerged vertical water jet impinging on an air/water free surface is experimentally and numerically investigated. On the one hand, a two dimensions-two components particle image velocimetry technique is carried out to perform velocity measurements from the jet injection to the unsteady free surface. The motion of this interface is recorded by a biaxial shadowgraph imaging system. On the other hand, volume of fluid and large eddy simulation methods are coupled to compute the interaction between the jet dynamics and the interface behavior. Numerical and experimental results are confronted in both the single-phase jet area and the impinging region beneath the free surface. A statistical analysis of the interface motion is led through the maximum elevation point to highlight its organized time and spatial evolution in terms of a characteristic lifetime. © 2009 American Institute of Physics. 关DOI: 10.1063/1.3085810兴 I. INTRODUCTION Hitherto, many studies have been devoted to turbulent two-phase flows which are increasingly used in industrial and environmental applications such as spray formation, coating techniques, wave breaking, or combustion. However, the complex turbulence/interface interactions are still not very well understood especially when the interface scales are of the same order of the large turbulent scales. In these turbulent two-phase flows, the scale separation between turbulence and interface is not satisfied and a strong coupling between turbulence and the interface deformations takes place. The present research is focused on the jet impingement technique in a confined medium that is widely used in industrial applications such as heating, cooling, drying, or mixing. This is one of the oldest and most attractive techniques to enhance the convective phenomenon in which heat and mass transfer occur. Within the past decade, many studies have been devoted to impinging jets on a solid plate1–3 but few investigations have been performed on jet impingement on a freely deformable interface.4–6 In this turbulent free surface flow, the challenge is to link the flow structures and their interaction with the free interface to characterize interface mass transfer and turbulent scalar transport mechanism. Gharib and Weigand7 focused on the behavior of a vortex ring approaching a free surface. They identified two mechanisms: the disconnection of the vortex ring under the free surface, followed by the connection with the interface. Many issues remain to be addressed, particularly concerning the interaction between turbulent phenomena and the interface motion. On the one hand, the interface topology is ruled by the turbulent structures which occur in each phase and, on a兲 Author to whom correspondence should be addressed. Electronic mail: jerome.larocque@cea.fr. Telephone: 共⫹33兲 05 57 04 50 62. 1070-6631/2009/21共6兲/065110/21/$25.00 the other hand, the turbulent structures of each phase are affected by the interface motion. If the interface is restricted to a free air/water interface, the dynamics of the upper fluid is often negligible as soon as the flow regime does not involve gas mixing in the liquid phase. In such a case, the shape and the motion of the interface result from a delicate balance between several forces: vertical acceleration and pressure below the free surface, gravity, and interfacial tension. The nature of the interface motion depends on which force predominates.8 In a case of a jet impingement on a deformable interface, Friedman and Katz4 also defined several distinct flow regimes which were identified primarily with the interface Richardson number related to the Froude number Ri = 1 / Fr2 = 共Di⌬␳g兲 / 共␳U2i 兲, where Ui and Di are defined as the jet velocity and the jet diameter at the interface.9 ⌬␳ and ␳ are the density jump across the interface and the fluid density, respectively. g is the gravity. In this paper, only the second regime defined by Friedman 共1.1 ⬍ Ri ⬍ 15兲 in which a strong deformable liquid-gas interface occurs 共without gas capturing bubbles in the liquid phase兲 is studied by experimental and numerical methods. To produce the air/water impingement, we consider an upward round turbulent water jet fed into a rectangular tank filled with water. Water enters in the tank from a nozzle injection pipe at Re = 1 ⫻ 104 共the Reynolds number is defined in terms of mean bulk inlet velocity Ub and inlet diameter d兲 and impinges on the free water/air surface where an unsteady moving bump is formed. The mean water level is kept constant, thanks to a water draining system, and the statistical steady turbulent state is reached. Experiments and simulations have been performed simultaneously in order to achieve first a validation of the numerical model in both turbulent “single-phase” and “interfacial” jet region, and second to provide quantitative data in the interesting case of a 21, 065110-1 © 2009 American Institute of Physics Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-2 Larocque et al. turbulent interfacial flow where the turbulent developed state strongly interacts with the free surface motion. In this region, the response of the oscillating air/water interface exciting by the turbulent jet has been examined. From an experimental viewpoint, it has been shown that surface behaviors change with jet velocity, lateral confinement, and impingement distance. Some typical behaviors have been reported. Madarame et al.10 observed self-induced oscillations caused by an upward jet, above a critical velocity value. Houard et al.11 investigated the dynamics of bidimensional jets and concluded that the oscillatory motion of surface bump below a velocity threshold is a materialization of helical instabilities which develop in the submerged laminar jet. This suggests that vortices are not only translating as a plane wave along the jet but also spiraling around the jet at the same time. Maurel et al.6 found three regimes for the free surface depending on velocities and cavity depth: stable, oscillations at a well-defined frequency, and a fountain regime with a more complex frequency spectrum. Bouchet et al.12 analyzed these oscillations in a similar plane configuration. Recently Zheng5 investigated oscillatory modes for different impingement distances. In water flow, major studies used particle image velocimetry 共PIV兲 to measure velocity.13 Whole-field measurements of velocities have also been performed in free jet flows.14–17 Nevertheless, as far as we know, the case of a free surface impinging jet has never been investigated with these whole-field techniques. In the present study, a PIV technique is used to carry out whole-field simultaneous velocity measurements in the case of a submerged axisymmetric jet impinging on a free surface. The free surface motion is also measured by a fast shadowgraph technique, thanks to two additional cameras synchronized on PIV. From a numerical viewpoint, no work is reported in literature concerning the impact of a turbulent jet on a free surface in the same configuration previously described. Nevertheless, different approaches are usually used to simulate turbulent interfacial flows. Direct numerical simulations 共DNS兲 are limited to simple configurations due to the high or excessive computational cost. However, they provide valuable information on the turbulence/interface coupling. The DNS of Bunner and Tryggvason18 reports that the turbulent transfer phenomena are modified near the interface in the case of a bubbly flow. Due to the necessity of reducing the computational cost of the numerical simulation and of characterizing the unsteady turbulent motion of the free air/water surface, the large eddy simulation 共LES兲 approach has been chosen. Contrary to the Reynolds averaged Navier–Stokes 共RANS兲 approach, the LES deterministic approach is able to control the amount of damping as it is based on a scale separation: the large scales are explicitly and completely resolved, whereas the small scales and their interactions with large scales are modeled by introducing closure assumptions. A functional point of view is adopted for the closure assumptions: the dynamics of the small structures serves only to dissipate the energy provided by the large scales of the turbulent flow. To estimate the turbulent dissipation, the mixed scale model 共MSM兲 共Ref. 19兲 is implemented in the eddyviscosity framework and is applied in each single-phase flow. Phys. Fluids 21, 065110 共2009兲 FIG. 1. Physical configuration: the confined turbulent round jet impinging on a free air-water surface. The origin of the reference frame 兵x , z , y其 is located at the center of the pipe exit. In addition, the interface capturing is ensured by a Eulerian method as the volume of fluid 共VOF兲 method which seems to be an attractive approach given its good representation of strongly deformed interfacial flow. Liovic and Lakehal20,21 tried to accurately predict and simulate turbulent interfacial flows in the framework of LES and Eulerian VOF methods. They proposed a specific law to damp and set the turbulent viscosity at the interface to zero. Their conclusions have underlined the fact that it is necessary to take into account specific subgrid terms that occur only in two-phase flows. This point has also been reported by Labourasse et al.,22 Vincent et al.,23 or Liovic Lakehal.20 In this study, the specific two-phase subgrid terms are neglected in a first approach at the interface. This assumption has to be verified and is discussed for the case of phase separation flows22,23 or turbulence bubble interaction.20,21 In Sec. II, the physical configuration, experimental, and numerical techniques are described in detail. In Sec. III, the LES numerical results are compared with the experimental measurements in the turbulent single-phase jet region. In Sec. IV, attention is turned to the turbulent air/water interface response which is analyzed. Conclusions and perspectives are finally drawn in Sec. V. II. GEOMETRY, EXPERIMENTS, AND NUMERICAL MODELING A. Physical configuration and length scales To study the interaction between turbulence and free surface, the physical configuration in Fig. 1 has been designed taking numerical and experimental constraints into account. The turbulence/free surface interaction is obtained by means of a vertical round water jet impinging on a free air-water surface. The turbulent round jet is generated at a Reynolds number Re = Ubd / ␯ = 1 ⫻ 104 in a rectangular tank which is initially filled with water on a depth H. The mean bulk inlet velocity Ub is equal to 1 m/s, while the value of the mean inlet centerline velocity U0 is 1.26 m/s. The pipe injection diameter d is equal to 0.01 m. The origin of the reference frame 兵x , z , y其 is located at the center of the pipe exit. The z axis is the axial direction of the mean upward jet. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-3 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 Before describing the experimental setup and the numerical modeling, the length scales occurring in this turbulent free surface flow should be kept in mind. It is possible to estimate the turbulent length scales such as the integral L⬁ and the Kolmogorov ␩K scales. Even if our jet cannot be considered as fully self-preserving because of a strong adverse pressure gradient due to the vertical confinement, it is interesting to get an idea of these scales in a case of a singlephase free self-similar turbulent round jet. From the experimental data of Wygnanski and Fiedler24 and Panchapakesan and Lumley,25 L⬁ ⬇ 0.17␦0.5 , 冉冊 冉 ␯3 ␩K ⬇ ⑀ 1/4 ␯3 = U3c /␦0.5 冊 共1兲 1/4 , where ⑀ is the dissipation estimated by the local mean centerline velocity Uc and the jet half-width ␦0.5. If we assumed that Uc ⬃ U0d共z − z0兲−1 and ␦0.5 ⬃ d共z − z0兲 with U0 as the inlet centerline velocity and z0 as the virtual origin of the jet 共⬃d兲, we obtain ⑀=A 冉 冊 U30 z − z0 d d −4 , 共2兲 where A = 48 is a constant determined from the experimental data of Friehe et al.26 Similarly, the respective time scales ␶⬁ and ␶K are obtained as ␶⬁ ⬇ L⬁/共0.1Ub兲, 冉冊 ␯ ␶K ⬇ ⑀ 共3兲 1/2 . It is important to note that the turbulent flow regime 共Re = 1 ⫻ 104兲 has been calibrated in such a way that large deformations occur at the interface without splitting or rupture. The unsteady deformations of the free surface are in the same order as the centimeter. Moreover, the Weber number based on the jet diameter and the inlet bulk velocity We = ␳U2bd / ␥ ⬃ 133 indicates that the effects of surface tension will be small relative to the convective effects. All characteristic length scales are shown in Table I. B. Experimental setup, techniques, and procedure 1. Setup The experiments were performed in a rectangular tank presented in Fig. 2. The tank 共440⫻ 440⫻ 300 mm3兲 is open at its top and the water level is kept constant in the test tank during the experiment thanks to an overflow device with four vertical pipes. The turbulent flow is generated by an upward submerged water jet issuing from an inlet pipe made of glass. This tube located at the bottom is 10 mm wide and 700 mm long to ensure a fully developed turbulent jet flow at the pipe exit. Indeed as previously demonstrated by Laufer,27 pipe flows become fully turbulent at 40d location. The confinement is realized, thanks to a parallelepiped cavity made of glass with a 10d ⫻ 10d cross section 共spanwise confinement L = 10d兲 TABLE I. Physical parameters. Turbulent time and length scales are estimated at z / d = 10. Parameter Symbol Value Re We Ub U0 1 ⫻ 104 133 1 m/s 1.25 m/s d H XxZxY 0.01 m 0.16 m 100⫻ 160⫻ 100 mm3 Kolmogorov length scale ␩K Integral length scale Kolmogorov time scale Integral time scale L⬁ ␶K ⬃0.04 mm ⬃1 mm ⬃1 ms ⬃10 ms Reynolds number Weber number Mean 共bulk兲 inlet flow velocity Mean centerline inlet velocity Inlet diameter Initial depth of water Confined medium ␶⬁ and a variable height, thanks to a false bottom 共streamwise confinement h = 16d兲. A horizontal overflow weir is used to keep the water level constant in the test tank during the experiment. Before any experiment, the flow pump is first connected to a closed-loop water circuit to ensure homogeneous initial conditions of temperature in the injected fluid. Then for injection, the water flow is sent to the jet tube using a control valve. Experiments are performed in water at 20 ° C in an airconditioned room. Temperature is controlled at ⫾0.1 ° C. Velocity measurements using DPIV were used to analyze turbulence in the tank and near the impingement. Visualization of the surface bump made it possible to determine the interface location. Shape and displacement DPIV measurements are carried out with two 60 mJ Nd:YAG 共yttrium aluminum garnet兲 lasers 共Quantel Twin Brio兲 with a pulse duration of 5 ns. The emitted wavelength is 532 nm. A combination of lenses focuses and expands the laser beam into a 2 mm thin vertical laser sheet 共Fig. 3兲. 2. Velocity measurements: PIV To maintain the accuracy of PIV measurements, the water is seeded with 4 ␮m nylon tracers 共rilsan A, a.k.a. PA12兲 which have a density of 1020 kg/ m3. The particles are four times smaller than the Kolmogorov scale 共between 15 and 100 ␮m in the jet兲. A charge coupled device camera 共TSI PIVCAM兲 is used to capture 2C-PIV images in a diametral plane of the flow. The camera produces 12 bit digital images with 1376⫻ 1024 pixels2 and the spatial resolution is 120 ␮m / pixel. The time delay between two laser pulses is 350 ␮s and a continuous image sequence of 150 images is obtained at a rate of 4 Hz. For each measurement location, seven sequences of images are recorded. Thus one measurement includes 1050 frame pairs during a total measurement time of 4 min 22 s, in agreement with a statistical convergence analysis. The statistical bias is ⫾0.3% in the shear layer, where fluctuations are maximum. The size of the interrogation areas is 16⫻ 16 pixels2, with overlapping interrogation windows 共typically 50%兲 in the mean two directions. A PIV cross-correlation process is completed with subpixel displacements estimated by means of a Gaussian Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-4 Larocque et al. Phys. Fluids 21, 065110 共2009兲 FIG. 2. Test tank, confinement cavity, and water loop. peak fit to the correlation peak for both horizontal and vertical displacement. PIV processing is optimized by using around ten particles in each interrogation window and images of the particles are distributed on 3–4 pixels for each particle. The common accuracy for the measured displacement is 0.05–0.10 pixel.28 The results contain between 1% and 5% spurious vectors, which could be detected by means of a median test and replaced by linear interpolation. Uncertainty on velocity measurements including measurement error and statistical bias is estimated less than 1%. It is important to note that if the camera axis is normal to the light sheet, the moving surface in front of the plane leads to an optical obstruction. The optimal way to obtain an unobstructed view of the interface is by looking upward from a slight angle below the surface. An angle of 1° from the camera provides good visualization of the interface and its turbulence beneath. Tests were performed with a calibration grid to quantify the effect of this slight angle on the velocity calculation. In our case we verify that the maximum error due to an angle of 1° is insignificant. As a result two different views are performed: a whole-field view of the jet in the cavity 共without any camera angle兲 and a close view of the interface from below 共with an angle of 1°兲, both combined for full analysis. 3. Free surface measurements: Shadowgraph imaging A shadowgraph technique 共SH兲 is used to visualize the behavior of the free surface. Two other cameras 共Vision Research Phantom V5兲 pointed at the bump are used synchronized together, ensuring two different perpendicular views with a spatial resolution of 200 ␮m / pixel. The complete experimental setup is depicted in Fig. 3. The first series are taken at a higher frequency 共200 Hz兲 to enable frequency analysis. Then both cameras are synchronized with PIV at a capture rate equal to 4 Hz. This makes it possible to analyze instantaneous quantities and to investigate more precisely the interaction between the submerged water jet and distorted free surface behavior. Visualization images of the surface bump are processed into binary images by setting an appropriate threshold value. The shape of the surface bump is extracted from these 2200 images, then providing the interface location, shape, and displacement. In this way, the maximum bump elevation point is tracked. Experimental conditions are summarized in Table II. TABLE II. Experimental conditions. Parameter FIG. 3. Experimental setup 共top view兲. Symbol Value Measured region Laser sheet depth Laser pulse delay 共PIV兲 Exposure time 共PIV兲 Recording frequency 共PIV and SH兲 Camera spatial resolution 共PIV兲 Recording frequency 共SH only兲 XxY lsd ⌬t ␦t 100⫻ 160 mm2 0.002 m 350 ␮s 5 ns f PIV RPIV f SH Camera spatial resolution 共SH兲 Total recording time 共7 ⫻ 150 frame pairs兲 RSH T 4 Hz 120 ␮m / pixel 200 Hz 200 ␮m / pixel 262.5 s Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-5 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 C. Computational methodology The numerical model is devoted to simulating turbulent two-phase flows of separated phases with a characteristic interfacial length scale of the same order of magnitude as the integral scale of turbulence. Both length scales are larger than the computational cell grid size. The two-phase numerical modeling is restricted to immiscible, incompressible and isothermal fluids. Moreover a constant surface tension is assumed. 1. Governing equations a. One-fluid approach. Let us consider two immiscible fluids, water and air referred to as k = 1 and k = 0. The variable ␾ will be indicated by ␾k in the phase k. The dynamics of each phase k is ruled by Navier–Stokes equations which control simultaneously the interface motion. Jump conditions 共mass and momentum conservation兲 at the interface are commonly added to the Navier–Stokes equation system in order to build a single set of equations valid in the whole simulation domain. As mentioned in the work of Kataoka,29 the simplification of the multifield formulation is ensured by the definition of a color function C which describes the interface evolution. This color function C is defined with the Heaviside function as follows: C共x,t兲 = 再 1 if x 苸 phase k = 1, 0 if x 苸 phase k = 0. 冎 FST共C兲 = ␥␬n␦i = ␥ ⵜ · Using the color function, local quantities such as volumic mass ␳, dynamic viscosity ␮, velocity u, and pressure p are defined by u = Cu + 共1 − C兲u , 0 冉 冉 G共x − y,t − t⬘兲⌽共y,t⬘兲dydt⬘ . 共8兲 −⬁ 冊 ⳵ ũ ¯ 关ⵜũ + ⵜTũ兴兲 + ũ · ⵜũ = − ⵜp̄ + ¯␳g + ⵜ · 共␮ ⳵t 共5兲 + 关␳1 − ␳0兴ũ␶III + ␶IV , 冊 ⳵u + u · ⵜu = − ⵜp + ␳g + ⵜ · 共␮关ⵜu + ⵜTu兴兲 + FST , ⳵t 共6a兲 ⵜ · u = 0, 共6b兲 ⳵C + u · ⵜC = 0, ⳵t 共6c兲 where FST is the surface tension force depending on the interface topology and is modeled using a continuum surface force method based on the work of Brackbill et al.,32 共9a兲 关␳1 − ␳0兴 III ␶ , ¯␳ 共9b兲 ⳵ C̄ + ũ · ⵜC̄ = − ␶III , ⳵t 共9c兲 ⵜ · ũ = 0 With respect to these definitions, the Navier–Stokes equations for two-phase flow are built in the framework of the one-fluid formalism,30,31 ␳ +⬁ + FST共C̄兲 − ⵜ · 共␶I − ␶II兲 p = Cp + 共1 − C兲p . 1 共7兲 It is assumed that the filtering operator commutes with time and spatial derivatives. Applying the filter G to the governing equations system 共6兲 leads to the filtered one-fluid equations system,22 ␮ = C␮1 + 共1 − C兲␮0 , 1 冕 冕 +⬁ ⌽̄共x,t兲 = ¯␳ ␳ = C␳1 + 共1 − C兲␳0 , ⵜC ⵜ C. 储ⵜC储 This single-field representation of the flow 共6兲 calls for certain comments. First, the so-called advection equation 共6c兲 describes the topological changes of the interface in the absence of phase change. Through this equation, the color function C is advected with the fluid velocity u and, consequently, the evolution of both phases 共air and water兲 is known. The interface is located by the color function discontinuity 关C共x , t兲 = 0.5兴. Second, the single velocity field of the interfacial flow is incompressible 共6b兲 due to the fact that each phase is incompressible. b. LES. To solve the one-fluid equation set 共6a兲–共6c兲, the computational grid cell should be fine enough to handle all turbulence length scales, i.e., the cell size should be smaller than the Kolmogorov length scale, leading to excessive computational cost in the configuration studied. To overcome this difficulty, LES formalism is used. In this approach, the largescale turbulent motions are resolved while the effects of the smaller ones are modeled 共see, e.g., Ref. 19兲. The scale separation is mathematically obtained by applying a convolution product between a large-scale-pass filter function G共x , t兲 and ¯ 共x , t兲 is defined as the quantity ␾共x , t兲. The filtered variable ␾ −⬁ 共4兲 冉 冊 where the Favre notation ⌽̃ = ␳⌽ / ¯␳ is used. Note that the Favre filtering 共.̃兲 is equal to the conventional filtering 共.兲 at locations where the fluid density is constant. The filtered ¯ are prescribed, thanks to arithmetic density ¯␳ and viscosity ␮ interpolations where the filtered color function C̄ is not constant, ¯␳ = C̄␳1 + 共1 − C̄兲␳0 , 共10兲 ¯ = C̄␮1 + 共1 − C̄兲␮0 . ␮ The filtering of Eq. 共6兲 nonlinear terms leads to the appearance of four unknown subgrid terms, Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-6 Larocque et al. ␶I = ␳u 丢 u − ¯␳ũ 丢 ũ, Phys. Fluids 21, 065110 共2009兲 共11a兲 ¯ 关ⵜũ + ⵜTũ兴, ␶ = ␮关ⵜu + ⵜ u兴 − ␮ 共11b兲 ␶III = u · ⵜC − ũ · ⵜC̄, 共11c兲 ␶IV = FST共C兲 − FST共C̄兲. 共11d兲 II T The subgrid stress tensor 共11a兲 and the subgrid viscous 共11b兲 terms are present in standard one-phase compressible LES equations, while the subgrid interfacial transport 共11c兲 and the subgrid surface tension force 共11d兲 are specific to twophase flows. 2. Subgrid scale modeling The magnitude of the different subgrid terms was a priori evaluated in the case of phase separation flows, turbulence bubble interaction,22,23 or spray atomization.33 Following this study, the subgrid stress tensor ␶I must be taken into account, while the influence of ␶II, ␶III, and ␶IV is highly dependent on the flow configurations and/or on the twophase description chosen. To date, a posteriori two-phase subgrid model evaluations20,21,34,35 are restricted to the ␶I term, while the other subgrid terms were neglected. In this work, the same hypothesis is adopted. Concerning the subgrid stress tensor, Lakehal et al.34 showed that, without modification, the standard Smagorinsky model becomes excessively dissipative near a sheared interface. This problem is overcome by using an empirical damping function34 or a variational multiscale model.35 Here, an autoadaptative viscosity model has been chosen. The subgrid stress tensor 共11a兲 is modeled by a turbulent viscosity assumption, ␶ = − ␮t共ⵜũ + ⵜTũ兲 = − 2␮tS̃. 共12兲 The MSM turbulent viscosity, proposed by Sagaut,19 accounts for the contribution of the resolved field gradients, of the kinetic energy of the highest resolved modes qsm and the ¯ = 共⌬x⌬y⌬z兲1/3. The subgrid cutoff LES filter length scale ⌬ viscosity is defined as ¯ 1+␣共2S̃ S̃ 兲␣/2关q 共x,t兲兴共1−␣兲/2 , ␮t = ¯␳CM ⌬ sm ij ij 共13兲 where ␣ is a parameter the value of which varies between 0 and 1, leading to a pure Smagorinsky model for ␣ = 1 or a pure subgrid turbulent kinetic energy model when ␣ = 0. In what follows, the value ␣ = 0.5 is chosen.36 C M is the mixed scale constant model evaluated from the standard Smagorinsky coefficient value and turbulent kinetic energy model leading to a value of C M = 0.06. qsm is computed using a filtering of the resolved velocity. The main interest of this model compared to the Smagorinsky model is that the eddy viscosity vanishes when the small resolved scale kinetic energy tends to zero without empirical input. The numerical damping of the turbulent viscosity near the wall or free surface is achieved in this way. 3. Numerical schemes The two-phase flow equations are approximated by implicit finite volumes of the second order in time and space on a fixed staggered Cartesian grid. An implicit augmented Lagrangian procedure37,38 is implemented for the treatment of the velocity-pressure coupling and the divergence free constraint. An iterative BiCGSTAB II solver,39 preconditioned under a modified and incompleted lower and upper 共LU兲 approach,40 is used to solve the linear system. The numerical approach dedicated to the treatment of interface tracking is the VOF method with total variation diminishing schemes.41,42 These schemes have been used for two reasons: 共1兲 Their implementation is simple because no reconstruction scheme is required. 共2兲 The numerical diffusion of the discontinuity is decreased in the three neighboring cells of the interface where interface tearing and stretching are predominant. 4. Geometry, boundary conditions, and grid considerations To reproduce similar water draining off conditions with respect to experiments, some mesh zones are penalized by a numerical treatment which prescribes a solid behavior in these zones. According to Khadra et al.43 and Randrianarivelo et al.,44 one possible penalty method consists in adding a Darcy term 共␮ / K兲ũ in the Navier–Stokes equation 共9兲 such that ¯␳ 冉 冊 ⳵ ũ ␮ ¯ 关ⵜũ + ⵜTũ兴兲 + ũ · ⵜũ + ũ = − ⵜp̄ + ¯␳g + ⵜ · 共␮ ⳵t K + FST共C̄兲 − ⵜ · ␶I , 共14兲 where K is the local permeability 共m2兲. It has to be noted that the filtered density and viscosity in Eq. 共14兲 are, respectively, interpolated with an arithmetic and geometric average. As mentioned by Ritz and Caltagirone,45 these averages are based on the conservation of mass and viscous flux through the interface with respect to its position in the cell. This method is called the Navier–Stokes–Brinkman penalty method. It assumes that the solid media are impermeable 共K → 0兲, whereas the fluid zones are free nonporous media 共K → ⬁兲. This approach could also be defined with ␮ → ⬁ in the solid zones and ␮ = ␮fluid in the related zones. We choose to play on the permeability to prescribe the solid behavior in the penalized zones. With this penalized technique, the Neuman conditions for water outlet in the LES simulation are more stable and have little or no effect on the motion of the turbulent interface. The penalized grid zones are represented in Fig. 4 by the 共x , y兲 points such as 5 ⱕ 兩x兩 / d ⱕ 7.5 or 5 ⱕ 兩y兩 / d ⱕ 7.5. The upper boundary condition of the calculation domain is a wall condition in order to assist the numerical stability in the corners of the upper air domain. All boundary conditions are shown in Fig. 5. The Cartesian mesh is fine and constant in the jet diameter 共12 grid points are resolved in the nozzle diameter兲 and is unrefined from the jet shear layer to the tank wall 共see Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-7 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 TABLE III. Mesh characteristics. The mesh is refined in the jet shear layer 共in the X and Y directions兲 and remains constant in the Z direction. Direction Dimensions 共mm兲 No. of points X Z Y 150 80 220 120 150 80 Grid space 共mm兲 (a) 0.75 1.8 0.75 Table III兲. The refinement is adjusted according to the integral length scale which is estimated analytically in the selfpreserving zone by Wygnanski and Fiedler24 by L⬁ = 0.17␦0.5, where ␦0.5 is the half jet width. This integral scale is about one millimeter in this zone. Moreover the grid size in the Z-direction is chosen to be constant and equal to 1.8 mm. 5. Averaging procedures (b) FIG. 4. 2D view of the calculation domain and the mesh: some zones 共5 ⱕ 兩x兩 / d ⱕ 7.5 or 5 ⱕ 兩y兩 / d ⱕ 7.5兲 are penalized in order to have an acceptable drain of water which goes out from the domain. The free air/water interface 共the dotted line兲 is initially located at z = 16d. Two time averaging procedures are used in the LES simulations. Because of the presence of two phases in this turbulent flow, some grid cells can contain air or/and water during the simulation. How can the time average in these specific “two-phase cells” be calculated? The first and classic averaging approach consists of doing the same as the one-fluid formulation which considers only one equivalent fluid for the whole two-phase domain. Whatever the fluid contained in the two-phase cells, the field of the variable ⌽共x , t兲 is continuous. Consequently, each time step contributes to the average of 具⌽共x , t兲典 as follows: 1 T 具⌽共x,t兲典 = 冕 t0+T ⌽共x,t兲dt, 共15兲 t0 where T is the integration time 共=5 s in our case兲 and t0 is the beginning of the averaging process 共=5 s in our case ensuring the steady is reached兲. The averaging time T is much longer than the timescale for the turbulent fluctuations which can be approximated by ␶⬁ 共Table I兲. The second averaging technique is based on a discontinuous view of the two-phase turbulent flow. In that way, the numerical procedure is performed in similar conditions to those of the experimental procedure which averages only the water phase PIV data. Water and air phase conditioned averages 关具⌽W共x , t兲典 and 具⌽A共x , t兲典兴 are introduced and calculated with a specific phase integration time 共TW and TA兲. This leads to the following expressions for the water quantities: 具⌽W共x,t兲典 = with TW = FIG. 5. Physical configuration: dimensions and boundary conditions. 1 TW 冕 1 TW 冕 t0+T H共C共x,t兲 − 0.5兲⌽W共x,t兲dt t0 t0+T H共C共x,t兲 − 0.5兲dt, 共16兲 t0 where H is the Heaviside function. The same procedure is used also for air quantities 具⌽A共x , t兲典 and TA. It must be noted that the continuous average is not equal to the sum of Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-8 Larocque et al. Phys. Fluids 21, 065110 共2009兲 TABLE IV. Pipe literature data. Author 共year兲 Data Re Symbol Aanen et al.a 共1999兲 Eggels et al.b 共1994兲 PIV DNS 5300 7000 䊐 䉭 Westerweel et al.c 共1993兲 Rudman and Blackburnd 共2006兲 Riviere and co-workerse 共2007兲 PIV DNS PIV 5450 11 000 1 ⫻ 104 + 〫 쎲 a Reference 14. Reference 53. Reference 77. d Reference 54. e References 49 and 50. b c (a) duce realistic inlet turbulent conditions.46–48 To impose correlated velocities at the inlet, the following procedure is applied: 共a兲 共b兲 (b) FIG. 6. Mean axial velocity field on the jet axis 共a兲 and rms axial velocity 共b兲 obtained with the continuous 共solid line兲 and phase conditioned 共dotted line兲 averaging procedures. phase conditioned averages. The continuous and phase conditioned average are related by 具⌽共x,t兲典 = TA TW 具⌽W共x,t兲典 + 具⌽A共x,t兲典. T T 共17兲 Figure 6 shows mean results for both procedures 共continuous and discontinuous兲 on the jet centerline. It is observed that each averaging procedure is equivalent until z / d ⬇ 16.5. Now the PIV measurements are only valid for z / d ⬍ 15.5; therefore the continuous averaging procedure is used for comparison between mean numerical and experimental results presented in Secs. II C 6 and II C 7. 6. Inlet turbulent conditions From a numerical point of view, the treatment of the inlet boundary conditions is very significant in terms of the accurate representation of the physical system. Synthetic turbulence generation methods have been developed to intro- Imposition of a mean velocity profile 具ui典 based on the mean experimental values.49,50 Artificial generation and superimposition of fluctuating velocity profiles based on the work of Klein et al.48 This procedure is based on the knowledge of the Reynolds tensor Rij and the generation of isotropic dimensionless fluctuations Ui. Rij is obtained from RANS simulation and Ui is generated with a random process which uses three specified integral scales L⬁1 , L⬁2 , and L⬁3 in the inlet plane. According to the work of Calmet and Magnaudet51 and Klein,52 the length scale normal to the inlet plane L⬁2 is taken to be equal to 0.2d, while the tangential length scales are approximated by 0.1d. Here are the literature references with which the results are compared 共Table IV兲. Figure 7共a兲 shows the normalized mean inlet velocity profile measured near the pipe exit, along the half pipe cross section. Thanks to the false bottom filled with water, the velocities are measured in the cylindrical glass tube with negligible optical deformations. The mean velocity profile follows the typical shape of fully turbulent pipe flows and is in good agreement with experimental data,14,49,50,53 DNS data,53,54 and semiempirical logarithmic laws.55 The maximum inlet velocity on the centerline U0 is 1.26 m/s with a bulk inlet velocity Ub of 1.02 m/s. The measured mean velocity profile is interpolated by a piecewise polynomial function and injected as a mean inlet boundary condition in the LES simulations. The Klein’s procedure is then applied and the fluctuating velocity profiles are reconstructed at each time step. Figure 7共b兲 shows the axial turbulence intensity from PIV measurements, DNS data,54 and LES data. The axial turbulence intensity is about 5% in the central region, which is consistent with the commonly admitted empirical relation uz,rms / U0 = 0.16 Re−1/8. The maximum of the PIV axial fluctuations is in the nearwall region with a peak at 15%, but is slightly shifted near the wall 共x / d ⬎ 2兲. This may be due to a lack of spatial resolution near the walls where the smallest structures are under 10 ␮m. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-9 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 FIG. 8. Numerical method for tracking the maximum bump elevation point. (a) points 共M 1 , M 2兲 and keep only the couple which has the maximum vertical coordinate zmax = Q1关z1 , z2 , C共M 1 , t兲 , C共M 2 , t兲兴. Thanks to this linear interpolation, only the coordinates xmax and y max are located on the mesh. III. RESULTS IN THE TURBULENT SINGLE-PHASE JET REGION A. Structure of the turbulent jet (b) FIG. 7. Validation of the synthetic generation of the inlet fluctuations: comparison between LES velocities 共⫺兲 and PIV measurements 共쎲兲. Mean axial velocity profile 共a兲 with the other experimental data of Aanen et al. 共Ref. 14兲 共䊏兲, Eggels et al. 共Ref. 53兲 共䉭兲, and Westerweel et al. 共Ref. 77兲 共+兲. The axial Reynolds stress 共b兲 with RANS 共Reynolds stress model兲 data 共− −兲, DNS data 共Ref. 54兲 共〫兲. 7. Tracking method for the maximum bump elevation point To follow the bump motion, the maximum bump elevation point is tracked by means of a specific numerical technique 共Fig. 8兲. This technique has to check accurately the maximum position 共xmax , zmax , y max兲 of the continuous 0.5 volume fraction isovalues predicted by the VOF model. The following procedure based on a linear interpolation in the Z direction has been used to spot the maximum bump elevation point. 共1兲 Look for and identify the pair of grid points 关M 1共x , z1 , y兲 and M 2共x , z2 , y兲兴 where C共M 1 , t兲 ⬎ 0.5 and C共M 2 , t兲 ⬍ 0.5 共empty circles in Fig. 8兲. 共2兲 Interpolate linearly the vertical coordinate of the 0.5 volume fraction isovalue with the position and the volume fraction of the points M 1 and M 2. 共3兲 Compute the interpolated elevation for every pair of Before studying the physical structure of the turbulent jet, we recall the particularities of this turbulent flow concerning the lateral and vertical confinement. Laterally, the water turbulent jet is confined in a rectangular box with a square basis of 10d. According to Djeridane’s experiments,56 the presence of such a confinement does not modify the mean jet dynamics up to 20d downstream. As expected, lateral confinement implies no effect on mean and rms quantities. This has been previously checked by our measurements without lateral confinement in the same flow conditions. Vertically, the water turbulent jet is limited by a free air-water surface initially located at z / d = 16. This latter plays a role similar to a blocking surface. That means the jet dynamics is expected to be strongly affected inside a specific zone beneath the unsteady free surface 关Fig. 9共b兲兴. The experimental and computed isocontours of mean velocity are compared in Fig. 9共a兲. The major part of the error relative to the axial inlet jet velocity U0 is mainly located in the transition zone and near the lateral walls 共⬃15%兲 where the computed mean axial velocity overpredicts the PIV measurements. These differences almost certainly originate from the large grid size, this effect will be detailed below. Usually the typical round free jet region can be subdivided into three zones: a zone of a conical jet core, a transition zone 共also called developing zone兲, and a zone of established flow 共ZEF兲 共also known as a developed zone兲. At the beginning of this latter, only the mean cross-sectional velocity profiles collapse onto a single curve and a partial selfpreserving state is reached. In our impinging turbulent round jet, the full self-preserving state 共valid in a free round jet beyond 70d兲 is not observed as all orders of the turbulent moments do not collapse. Concerning the mean value, it is widely assumed that a self-similar behavior occurs in the downstream region of a Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-10 Larocque et al. Phys. Fluids 21, 065110 共2009兲 cerning numerical60–64 and experimental13,16,17,65,66 turbulent free jets. It this way, George67 claimed that the self-similarity of a jet is not universal as the self-similar constants depend on the initial conditions. B. Mean streamwise evolution (a) (b) FIG. 9. Isocontour of numerical 共dotted lines兲 and experimental 共solid lines兲 mean axial velocity field in the plane 共x , z兲 共a兲 and simplified expected view of the jet structure 共b兲. free jet. In that case, the mean axial velocity profiles in this region Uz and the spreading rate of the jet ␦1/2 can be expressed as follows: U0 Cu = 共z − z01兲, Uz共z兲 d U z共 ␰ 兲 = exp共− K2u␰2兲, Uc ␦1/2共z兲 d = C␦ 共z − z02兲, d Figure 10共a兲 shows the numerical and experimental velocity decays along the jet axis. This curve exhibits the structure of the present impinging jet which can be subdivided into four zones: a jet core to 4d, a transition zone from z / d = 4d to 7d, then a ZEF 共horizontal stage from 7d to 13d兲 and finally the impinging region from 13d to the mean free surface location. The derivative of U0 / Uz along the jet axis highlights clearly all previous identified zones in Fig. 10共b兲. The numerical simulation overestimates the length of the jet core. According to a 95% criterion, the sharp increase of U0 / Uz takes place around z / d = 4 for PIV data and around z / d = 5 for LES results. Consequently, small discrepancies appear in the transition zone and in the ZEF zone. Figure 11 shows various mean axial velocity profiles perpendicular to the jet axis which collapse onto a single Gaussian profile represented by Eq. 共18b兲. Consequently, partial selfpreserving state is reached at around 8d. The first line in Table V shows the experimental range of these constants Cu, Ku, and C␦. Because the jet broadens linearly the cross-sectional spreading rate C␦ can easily be determined, thanks to the axial evolution of the half-width of the jet ␦0.5 关Fig. 10共c兲兴. In the experimental case C␦ = 0.081共⫾0.002兲, corresponding to a velocity conical expansion of 4.4°. The cross-sectional virtual origin is z02 = 0. The numerical jet spreading rate C␦ = 0.095 and velocity conical expansion of 5.4° are higher than the experimental values due to the overprediction of the jet core length. Whatever happens, the computed values remain in the range of those of other round jet studies. Focusing in the ZEF where theoretical analysis is available, the numerical and experimental constants related to the set of Eq. 共18兲 are compared to theoretical values in Table V. Three other theoretical algebraic relations according to the self-similar analysis of a free round jet can be built and are also in Table V, 共18a兲 C2uK2u = where ␰ = x/共z − z02兲, 4 ln共2兲 3共冑2 − 1兲 ⬃ 2.23, 共18b兲 C␦2K2u = ln共2兲 ⬃ 0.69, 共18c兲 where Cu, Ku, and C␦ are three self-similar constants which are characteristics of the jet dynamics in the ZEF. Cu is a characteristic of the axial decay rate, while Ku and C␦ represent the radial expansion of the turbulent round jet. z01 and z02 are the virtual origins of the jet. The values of these constants are strongly dependent on the inlet conditions such as the nozzle geometry,57 the mean inlet velocity profile,58 or the inlet velocity fluctuations.59 This can explain the scatter of the self-similar constant values in numerous studies con- 共19兲 C␦2/C2u = 0.75共冑2 − 1兲 ⬃ 0.31. Concerning the velocity decrease in the zone of established flow, this one can be fitted by a z−1 decay. As indicated above, determination of the self-similar velocity decay constant is more accurate with the derivative of velocity inverse U0 / Uz giving Cu = 0.14共⫾0.004兲 experimentally and 0.17 numerically. The axial virtual origins are z01 = −2d and z01 = d, respectively. All self-similar constant values are experimentally and numerically in relative good agreement with the other studies of turbulent round jets. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-11 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 (a) (c) (b) (d) FIG. 10. Inverse of mean axial velocity field on the jet axis 共a兲, identification of the different zones in the jet, thanks to the derivative of U0 / Uz 共b兲, jet half-width ␦0.5 normalized by the inlet jet diameter d 共c兲, and axial rms velocity along the centerline 共d兲. The axial rms velocity from the LES calculation is compared to the present PIV measurements 共쎲兲 and the experimental data of a turbulent round free jet configuration 关Giralt et al. 共Ref. 68兲 Re = 80 000 共䉭兲 and Quinn 共Ref. 57兲 Re = 184 000 共䊏兲兴. The position of the experimental 共− · −兲 and computed 共⫺兲 mean free surface is reported on the subplots 共a兲 and 共c兲. Axial rms streamwise velocity evolution along the axis are plotted in Fig. 10共d兲. Qualitatively, measurements and LES show the same behavior. The fluctuation intensities are approximately constant in the jet core, then suddenly increase in the developing zone, and finally decrease in the self-similar one. When approaching the free surface, the axial turbulent intensity increases again due to the impact of the upward and downward fluid just beneath the mean bump. It should be noted that only numerical results predict this increase because of a lack of measurement near the unsteady free surface in this confined configuration. This increase is also observed in the case of an impinging jet in an unconfined configuration.49 Comparison of the measurements of Giralt et al.68 and Quinn57 with the present shows an important scatter of the axial evolution of rms velocity fluctuations. It is in fact well known that the jet development is particularly sensitive to the inflow conditions which govern the primary Kelvin– Helmholtz 共KH兲 instabilities69 and secondary 共longitudinal vortex兲 instabilities70,71 from the jet core starting at the end of the transition zone. In the present configuration, the tur- bulence inflow generated in the pipe limits the KH vortices growing and leads to a smoother transition from the potential area to a fully turbulent state. Consequently, the maximum rms level reached is approximately 30% and 15% lower than in the experiments of Giralt et al. and Quinn experiments where the inflow is quasilaminar. In the LES, the synthetic inflow is nonphysical and the velocity fluctuations decrease after the injection leading to delay the transition. Consequently, the maximum fluctuation location is shifted from z / d = 7 in PIV measurement to z / d = 8 in LES results. This spurious effect is amplified by the too large synthetic fluctuation length scale to grid size ratio 共L⬁1 / ⌬x ⬇ 0.3 and L⬁2 / ⌬z ⬇ 1.1兲. To overcome this numerical difficulty, a refined mesh of a larger domain including part of the inlet pipe and/or a more realistic inflow velocity condition 共see Na and Moin,72 for example兲 are needed. Despite this limitation, it should be emphasized that most of the turbulent jet features including qualitative behavior and quantitative self-similar coefficient evaluations are correctly predicted by the LES. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-12 Larocque et al. Phys. Fluids 21, 065110 共2009兲 axis. Turbulent radial and axial intensities are intense in the jet and around 3% outside the main flow. Shear stresses feature high levels in the shear layer, resulting in two peaks at x / d = ⫾ 0.5 in the extracted slice of the flow. 1. Jet core: z / d = 2 As might be expected in the jet core, the viscous effects are negligible compared to the inertial phenomena. The computed mean velocity components fit perfectly with experimental data which are conserved until the beginning of the transition zone. Discrepancies are noteworthy on the velocity fluctuations. Indeed numerical rms and Reynolds stresses seem underestimated in this zone. Experimental centerline turbulent intensities are 5% at z / d = 2, but numerical rms velocities are weaker in the jet core. Nevertheless both experimental and numerical data show higher levels in the shear layer than on the centerline. (a) 2. Transition zone: z / d = 6 In the transition zone 共z / d = 6兲, the rise in instabilities results in a sudden increase of the turbulence level. It is observed that the computed velocity fluctuations increase faster than the experimental data, resulting in higher shear stresses and axial turbulent intensities. The transition to turbulence is, in fact, very difficult to predict accurately since complex interactions occur between the well-known varicose and helical modes.73 As far as mean velocities are concerned, good agreement is observed on radial spread and axial decrease. 3. Zone of established flow: z / d = 12 (b) FIG. 11. Self-similar behavior in the turbulent round jet: 共a兲 LES results, 共b兲 PIV measurements. Radial profiles of the mean axial velocity Uz / Uc = f共␰兲 are plotted for z / d = 8 共䊐兲, z / d = 10 共䉭兲, z / d = 12 共〫兲, and z / d = 14 共䊊兲. These profiles are approximated with a common Gaussian law f共␰兲 = exp共−K2u␰2兲 represented by the solid line. C. Mean cross-sectional evolution Figures 12 and 13 present cross-sectional mean axial velocity profiles at different distances from the jet injection, showing the jet dynamics from the inlet to the impinging region where the jet turns into radial flow beneath the free surface. The mean axial velocity follows a symmetrical bellshaped profile while the mean radial velocity profiles are relatively antisymmetrical, with a change in sign on the jet In the zone of established flow 共z / d = 12兲, the mean velocities broaden linearly with the constant spreading rate previously determined. Numerical axial profiles are slightly wider than the experiment profiles. Moreover, numerical rms and shear stresses remain 30%–40% higher. This overestimate can be connected with the jet core length also overestimated in the LES calculation. The join of the shear layers located at the end of the jet core governs the transition of the jet. For that reason, the computed transition occurs later and more suddenly leading to higher values for turbulent intensities at the beginning of the ZEF zone. 4. Zone influenced by the free surface motion: z / d = 15 In the impinging zone axial velocity profiles 共z / d = 15兲 evidence a strong decrease in the centerline and negative velocities near the vertical walls due to lateral vortices located between the jet and the free surface. The level of tur- TABLE V. Comparison of the self-similar constants. Literaturea PIV measurements LES results Cu Ku C␦ C2uK2u C2␦K2u C2␦ / C2u 关0.14, 0.18兴 关9, 10兴 关0.080, 0.093兴 2.23 0.69 0.31 0.14 0.17 10.3 8.8 0.081 0.095 2.08 2.24 0.70 0.70 0.33 0.31 a References 13, 16, 17, 65, and 66. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-13 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 (a) (d) (b) (e) (c) (f) FIG. 12. Flow characteristics in the jet core at z / d = 2 共left兲 and in the transition zone at z / d = 6 共right兲: Mean 共Expt. 쎲, Num. ⫺兲 and rms 共Expt. 䊊, Num. − · −兲 axial and radial velocity 关共a兲–共d兲兴 and shear stress 具ux⬘uz⬘典 共Expt. 䊊, Num. − · −兲 in the subplot 共c兲 and 共d兲. All quantities are normalized by the inlet jet centerline U0. bulence intensities and shear stress is correctly predicted except for the boundary layers near the lateral walls where downward jets are formed. A huge difference is also noted for the radial mean velocity at z / d = 15. The radial water draining is overpredicted and the so-called “surface currents” are wrongly computed. Surface currents 共i.e., the large transverse velocity parallel to the free surface兲 occur when a turbulent flow such as a jet or wake interacts with a free sur- face. As recently shown by the studies of Walker and Hong,74,75 these result from a combination of the Reynolds stress anisotropy and the free surface fluctuations at a high Froude number. This could indicate that the correct prediction of the sub-surface turbulence requires 共a兲 an accurate model of the convective two-phase subgrid term ␶I 共assuming that ␶II and ␶IV are negligible in the Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-14 Larocque et al. Phys. Fluids 21, 065110 共2009兲 (a) (d) (b) (e) (c) (f) FIG. 13. Flow characteristics in the zone of established flow at z / d = 12 共left兲 and in the region influenced by the free surface at z / d = 15 共right兲: Mean 共Expt. 쎲, Num. ⫺兲 and rms 共Expt. 䊊, Num. − · −兲 axial and radial velocity 关共a兲–共d兲兴 and shear stress 具ux⬘uz⬘典 共Expt. 䊊, Num. − · −兲 in the subplots 共c兲 and 共d兲. All quantities are normalized by the inlet jet centerline U0. 共b兲 case of free surface flow at a high Froude number兲 which contributes to the prediction of the Reynolds stress anisotropy and an accurate model of the advection two-phase subgrid term ␶III which contributes to the prediction of the filtered instantaneous free surface C̄共x , t兲, i.e., the free surface fluctuations. IV. ANALYSIS OF THE TURBULENT AIR/WATER INTERFACE During the stationary state, the turbulent bump moves very quickly and its motion is strongly linked to the dynamics of turbulent structures which impinge the air/water interface. It can therefore be justifiably considered that a correct representation of the turbulent bump motion cannot be ob- Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-15 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 dimensions. The symmetry of the free surface with respect to the jet axis shows and confirms that the statistics converge. Experimental shadowgraph visualization and LES results for the free surface are compared in the plane 共x , z兲 关Fig. 14共b兲兴. A hydraulic jump can be observed in this plot, located near the draining off. On the jet axis, the mean vertical coordinate of the free surface 共具zI典 = 170 mm or ␩ = 10 mm兲 is correctly predicted by the LES calculation on a 80⫻ 120⫻ 80 grid with an error of 11% compared to the experimental data 共具zI典 = 172 mm or ␩ = 12 mm兲. The value of ␩ can be linked with a phenomenological law4 in function of the interface Richardson number Ri, which is defined in terms of interface velocity Ui and spreading jet diameter Di estimated at z = H. This latter is only valid for Ri ⬎ 2. In that case of particular range of Richardson number, the water flow is considered on average to be practically radial outward and it is possible to estimate the aspect ratio ␩ / Di. Assuming the equilibrium between potential and kinetic energy 关u ⬃ 共2g␩⌬␳ / ␳兲1/2兴, it can be deduced ␩ Di = C共Ri兲−1/3 , 共20兲 where C is an appropriate empirical constant equal to 0.72 in the work of Friedman and Katz.4 Taking Ui ⬇ 0.25 m / s and Di = 4d with respect to the numerical data, we obtain Ri ⬇ 6.3 and ␩ ⬇ 15 mm for the height of the bump. The numerical and experimental free surface elevation ␩ = 10 and ␩ = 12 mm are in good agreement with this phenomenological value of 15 mm. From this, it can be thought that the mean deformation of the free surface is correctly predicted in our case. The instantaneous free surface must now be investigated to connect it with the subsurface turbulence. B. Fluctuating maximum position signals (b) FIG. 14. 3D Mean free surface elevation represented by the 0.5 isovalue of the filtered color function C̄ 共a兲 and comparison between mean experimental 共쎲兲 and numerical 共⫺兲 free surface elevation 共m兲 in the sketch 共x , z兲 共b兲. tained without a reasonable prediction of the turbulent dynamics under the free surface. In a first step, we therefore chose to focus our analysis on the nonstationary motion of the bump maximum elevation point. In this section, a time analysis is investigated and numerical LES predictions of the bump motion are compared to the experimental visualization obtained from two cameras. A. Mean free surface elevation The first investigated quantity is the mean free surface. The mean topology of the bump is directly linked to the mean flow which occurs just below the free surface. It is therefore only natural to look at the mean deformation ␩ of the interface caused by the jet, defined as ␩ = 具zI典 − H, where 具zI典 is the mean vertical coordinate of the free surface. In Fig. 14共a兲, the numerical mean free surface estimated with the 0.5 isovalue of the filtered color function C̄ is plotted in three The three-dimensional 共3D兲 position 共xmax , zmax , y max兲 of the maximum bump location is obtained numerically and experimentally with different techniques which have been detailed previously. The three time signals characterize the bump oscillations which depend on the subsurface turbulent dynamics. Two basic oscillations have been identified in literature:5,10 fluttering or horizontal oscillations and swelling or vertical oscillations. The experimental and numerical time acquisitions are, respectively, equal to 11 and 5 s with the following frequencies f = 200 Hz and f = 2000 Hz. The twodimensional 共2D兲 top view of the maximum bump location is plotted in Fig. 15 which shows that the maximum bump is located approximately inside a circle with a radius around 3d. Some important deviations from the inside of the previously defined circle are also observed: the maximum bump cartography resembles like a pointed star. This shape is observed by both simulation and experiments. 1. Statistics The maximum bump location has been first analyzed in terms of statistics. The statistical behavior of this particular point gives valuable information on the interaction between the subsurface turbulence and the interface motion. In fact, a Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-16 Larocque et al. Phys. Fluids 21, 065110 共2009兲 (a) FIG. 16. One dimensional PDF of the angle parameter ␪max 共Expt. 䊊, Num. 쎲兲 compared to a uniform distribution 共−兲 between 关−␲ , ␲兴. 冉冑 冊 ␪max = a cos 冦 (b) FIG. 15. Maximum bump location in the plane 共x , y兲. correct prediction of this statistical behavior would involve a reasonable prediction of the effects of turbulence on the free surface. The first and second order statistics of the maximum bump position and the angle parameter ␪max are summarized in Table VI. It is clear that the statistical analysis in the plane 共x , y兲 has to be performed in the radial and azimuthal direction due to the nonhomogeneous phenomena in the X and Y direction. Hence, the radius and angle parameter rmax and ␪max are defined, respectively, with the experimental and numerical center bump, ⬘ xmax rmax2 , where ⬘ = xmax − 具xmax典 xmax ⬘ = y max − 具y max典 y max ⬘ 兲2 . ⬘ 兲2 + 共y max rmax = 冑共xmax 共21兲 冧 Note that the rms radial values are calculated taking into account a zero mean value corresponding to the experimental and numerical center bump. Through these angle parameter ␪max, the influence of the boundary conditions 共the draining off兲 and confinement is investigated. In the case of a negligible confinement, the probability density function 共PDF兲 of ␪max is determined with a uniform distribution between 关−␲ , ␲兴 and can be defined as 冦 1 if ␪ 苸 关− ␲, ␲兴 PDF共␪兲 = 2␲ 0 elsewhere, 冧 共22兲 The mean and rms values corresponding to such a distribution are, respectively, equal to 0 and 冑共2␲兲2 / 12⬃ 1.813. It can be seen in Fig. 16 that the statistical experimental 共empty symbol兲 and numerical 共full symbol兲 angle distribution tend toward the previously defined uniform distribution 共solid black line兲. The first and second moments of the experimen- TABLE VI. Mean and rms values concerning the position of the maximum bump elevation point. All metric quantities are dimensionless with the jet diameter d and denoted by a star. The mean and deviation angles are calculated in radians. 具xmax典쐓 具y max典쐓 具zmax − H典쐓 具␪max典 쐓 rrms ␪rms 쐓 共zmax − H兲rms Expt. data ⫺0.404 0.362 1.254 0.138 1.025 1.837 0.215 Num. data ⫺0.025 0.089 1.525 ⫺0.087 1.330 1.761 0.348 Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-17 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 tal and numerical results are also comparable to the values of the uniform distribution 共Table VI兲. Thanks to this analysis, the confinement is assumed to be negligible as was to be expected from the previous section. Table VI also evidences that the maximum bump location is not exactly centered with the jet axis due to the experimental setup inaccuracies. The experimental mean deviation angle from the center estimated by arctan共xmax / zmax兲 is about 1.3°, while the numerical one is about 0.3°. These values are considered as negligible and do not influence the jet dynamics. Besides, both numerical and experimental data predict larger radial fluctuations rrms compared to the vertical fluctuations zrms. Two phenomena contribute to damping the axial bump oscillations. 共1兲 The gravity effects which are substantial compared to the inertial effects 共Fr = 冑1 / Ri ⬃ 0.4兲. 共2兲 The decrease in the axial energy in the subsurface region. Near the free surface, the axial energy provided by the upward turbulent jet is redistributed in the radial and azimuthal directions. Figure 17 shows one-dimensional PDF of the radial and axial position signals. Both PDFs are normalized so as to have zero mean and rms equal to 1. These are then compared to normalized Gaussian curves. It is clear that both the experimental and numerical PDF of the axial component zmax are closer to a Gaussian distribution 关Fig. 17共b兲兴. This Gaussian distribution of the free surface elevation is also observed by Hong and Walker75 at a high Froude number in the case of a turbulent round jet which evolves in parallel to a free surface. This Gaussian behavior is almost certainly originated by the random behavior of the turbulent structures which impinge and distort the free surface. Contrary to the axial direction, the PDF of the radial position of the maximum bump deviates from a Gaussian behavior 关Fig. 17共a兲兴. That means the maximum bump positions seem to be correlated in the radial direction. The radial motion of the bump maximum is not random. As the bump maximum is sometimes away from the mean center, it may be interesting to plot its altitude as a function of the distance from the center. Given that the instantaneous positions of the bump maximum do not depend on the azimuthal coordinate ␪, it can be useful to average the position of the bump maximum over the radial direction. Consequently, the mean vertical coordinate over radial direction 具Z共r兲典r with Z = zmax共r兲 − 具zmax共r兲典 is introduced as follows: 具Z典r共r兲 = 1 Ar 冕 冕 r+dr r 2␲ Z共r兲rdrd␪ , 共23兲 0 where Ar is the area between r and r + dr. Figure 18 shows that 具Z典r共r兲 is higher when the bump maximum is located near the center 共r → 0兲. However, the standard deviation Zrms共r兲 = 具共Z − 具Z典r兲2典r remains the same 共2 mm兲 wherever the maximum is even when it is located far from the center. This is certainly to be linked with the flattened profile of uz,rms approaching the free surface: fluctuations are almost identical in Fig. 13共b兲 for −2 ⬍ x / d ⬍ 2. (a) (b) FIG. 17. One dimensional normalized PDF of the maximum bump location in the radial direction 共a兲 and in the axial direction 共b兲 共Expt. 䊊, Num. 쎲兲 compared to the normalized Gaussian distribution =N共0 , 1兲 with zero mean and standard deviation equal to one 共⫺兲. The normalized variable ⌽쐓 is defined, respectively, with numerical and experimental data statistics such as 共⌽ − 具⌽典兲 / ⌽rms, where ⌽ = rmax or zmax. 2. Frequency analysis A frequency analysis was performed with each position signal 共in the radial, azimuthal, and axial directions兲 by means of fast Fourier transform 共FFT兲. Due to the topology of the free surface, the horizontal oscillations of the bump maximum in the X and Y directions are discontinuous, whereas the vertical oscillations remain continuous. No characteristic frequency was found in all signals. We remind that experimental measurements are done at a frequency of 200 Hz. FFT of the vertical position signal is plotted in Fig. 19. Only low frequencies 共⬍10 Hz兲 are observed. LES calculations and experimental measurements seem to be in agreement around this frequency range. Contrary to Zheng,5 no peak frequency has been clearly found. This is due to the Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-18 Larocque et al. Phys. Fluids 21, 065110 共2009兲 (a) FIG. 19. Power spectral density of the vertical maximum bump position zmax. The dotted line represents the measured vertical signal, while the solid line is the computed vertical coordinate of the bump maximum point. (b) FIG. 18. Mean 共a兲 and rms 共b兲 vertical coordinate of the maximum bump Z = zmax共t兲 − 具zmax典 in function of the radial coordinate normalized by the jet diameter r / d. large nozzle-interface spacing z / H which is equal to 16 in the present impinging jet experiment. In our experimental conditions, the free surface is located in the fully developed turbulent jet region. The turbulent dynamics in this region is different compared to the transitional jet region where characteristic frequencies are easily identified. In the fully developed region, it is much more difficult to identify any peak frequency as in the transitional jet region. Similarly, Zheng did not observed an apparent peak frequency for H / d ⬎ 7 using FFT algorithms. C. Characteristic time scale of the free surface motion In this section, the path lines of the maximum bump location are investigated. In fact, the trajectory of this particular point can be a history indicator of the subsurface turbulence. As previously indicated, the path lines are discontinuous in the plane 共x , y兲. Consequently, we define a cutoff length criterion ⑀ between two successive positions to iden- tify the continuous paths of the maximum bump. The value of ⑀ is expected to be greater than the spatial discretization length. In this way, several positions can be connected in the horizontal plane with the following procedure: if the distance ds between two successive positions is less than the ⑀ criterion, these positions are assumed to be continuous and linked with a solid line 共Fig. 20兲. After performing the procedure for each position and for three a priori length criteria 共⑀ = 2, 4, and 6 mm兲, continuous paths are identified and counted. These are called and defined as “filaments.” Considering the computed and measured path lines, it appears that the first value of the continuity criterion ⑀ = 2 mm is the most relevant. The other two values are too high so that discontinuous trajectories are wrongly connected. All filaments are characteristics of the unsteady free surface. Because of the statistical steady state of the turbulence dynamics, mean characteristic scales 共lifetime and length scale兲 of all filaments are estimated as follows: 冓 冔 冓 冔 N−1 具␶ f 典 = s⌬t 兺 s=1 , f 共24兲 N−1 具L f 典 = ds 兺 s=1 , f where N is the number of points for one particular filament, ⌬t is the time interval corresponding to the frequency acquisition, and 具 · 典 f is the average over all identified filaments. Table VII presents a list of all filaments in function of their characteristic lifetime ␶ f and the length criterion ⑀. Both experimentally and numerically, an increase in the length criterion involves a merging of small filaments 共corresponding to small lifetimes兲 which turn into long filaments 共corresponding to high lifetimes兲. These merging are unphysical for ⑀ = 4 and 6. Consequently only the criterion ⑀ = 2 is kept and used for the averaging over all filaments. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-19 Macroscopic analysis of a turbulent round liquid jet Phys. Fluids 21, 065110 共2009兲 TABLE VIII. Estimation 共for ⑀ = 2兲 of the mean characteristic length scale of a filament: the mean lifetime 具␶ f 典 共ms兲 and the mean length 具L f 典 共mm兲 No. of filament 具␶ f 典 具L f 典 276 152 27 23 4.1 3.3 Expt. Num. Another way to estimate the mean lifetime of the free surface is to use the time correlation function R of the bump maximum coordinates 共rmax , zmax兲, Rr = 具rmax共t兲rmax共t + ␶兲典 , 具rmax共t兲rmax共t兲典 (a) 具zmax共t兲zmax共t + ␶兲典 . Rz = 具zmax共t兲zmax共t兲典 共25兲 Figure 21 shows the time correlation function of the maximum bump radial and axial coordinates. All correlations decay very rapidly. These functions are now negative, now positive that means the convergence is not reached. The bump seems to be correlated several times during 1 s. Even if these time correlation functions are not converged, it is possible to calculate a characteristic integral lifetime. According to Piquet,76 this one is approximated with the time integration of R to the first zero of the function at time t1. ␶f = 冕 t1 R⌽dt, 共26兲 0 (b) FIG. 20. Maximum bump location in the plane 共x , y兲 and example of filament identification. Table VIII shows the results concerning the estimation of mean lifetime and length scale. A reasonable agreement is found between the numerical and experimental values of 具␶ f 典 and 具L f 典. where ⌽ = r or z. The computed and experimental lifetimes are compared in Table IX. The values in Table IX have to be compared to the averaged values over all filaments 具␶ f 典. The mean lifetime of a filament belongs to the range 关20,30兴 ms with the averaging process over all filaments, while the directional lifetimes ␶z and ␶r are estimated in the range 关20, 40兴 ms using the correlation functions. Consequently, a reasonable agreement is found between both procedures. TABLE VII. Filament counting obtained from the experimental and numerical data of the maximum bump location. ⑀ 共mm兲 and ␶ f 共ms兲 are, respectively, the distance criteria and the characteristic time of an identified filament. ␶f ⑀ 关0 ; 10关 关10; 30关 关30; 50关 2 61 4 6 31 20 2 共b兲 Experimental filament counting 74 95 66 41 4 6 25 11 69 74 共a兲 Numerical filament counting 48 26 35 20 44 28 21 19 35 25 ⬎50 17 28 34 FIG. 21. Time correlation function of the maximum elevation point in the axial 共Expt. –䊊– and Num. – –兲 and radial 共Expt. –䊐– and Num. ⫺兲 directions. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-20 Larocque et al. Phys. Fluids 21, 065110 共2009兲 TABLE IX. Comparison of the computed and experimental bump lifetime 共ms兲. Expt. Num. ␶z 共ms兲 ␶r 共ms兲 48 42 44 24 V. CONCLUDING REMARKS In the current paper, a submerged vertical turbulent jet impinging on a water/air free surface has been investigated experimentally and numerically. Simultaneous PIV and shadowgraph imaging techniques have been used to measure the velocity fields in the liquid phase and the interface motion. LES-VOF simulations have been performed to predict the turbulence/interface interaction and to asses the reliability of LES computations in a case of free surface flows. The main objective was to lead a complete comparison between numerical and experimental results and to provide new insight of interfacial motion under turbulent jet action. First the development of the turbulent round jet is characterized by mean and fluctuating velocity profiles. Although some discrepancies appear in the jet core and transition zones of the jet, the experimental and computed self-similar mean dynamics are in a relatively good agreement. The associated self-similar constant values are found in the range of other studies13,16,17,65,66 on free turbulent jets. The partial self-similar behavior of the jet is stopped by the presence of the free surface which leads to an impinging region. In this region located two jet diameters beneath the free surface, the mean axial velocity strongly decreases while turbulence intensities are significantly affected. Concerning the fluctuating quantities, a good trend is observed. The axial velocity fluctuations are reproduced with a maximum gap of 40%. Some differences appear in the jet development zone and are related to the fluctuating synthetic inlet conditions. A specific parametric study on integral length scales should be leaded in order to improve the quantitative effects of this turbulent inlet conditions. Second an original analysis of the unsteady free surface motion is provided with both numerical and experimental approaches. The tracking of the interface is achieved with the identification of the maximum bump location point. On the one hand, the computed mean free surface shape is correctly predicted compared to the measured interface shape. On the other hand, a statistical analysis of the maximum bump location shows a Gaussian distribution of the vertical oscillations, subjected to the random axial turbulent phenomena. This Gaussian behavior has also been observed by Hong and Walker75 in the case of a different turbulent free surface flow. Unlike the vertical motion, the horizontal one shows spatially and temporally organized patterns. No characteristic frequency of these oscillations is observed. That means no self-sustained oscillation appears in the present impinging turbulent jet. Despite of this, a characteristic time scale of the bump motion, connected to the interfacial patterns, is estimated. The experimental visualization and the LES calcula- tion of the free surface show a similar interface behavior. The mean elevation of the maximum bump location decreases when this one moves away from the center, while the fluctuating elevation remains constant due to the flatness of the rms velocity profile beneath the interface. Time velocity signals will be analyzed in terms of spectral density thanks to new measurements by high speed laser Doppler velocimetry and numerical probes located near the interface. In that way a finest comparison will be carried out concerning the turbulence behavior beneath the free surface. Some developments including turbulence model of the twophase subgrid term which appears in the advection equation will be proposed in the near future. 1 H. R. Gardon and J. C. Akfirat, “Heat transfer characteristics of impinging two-dimensional air jets,” Trans. ASME, Ser. C: J. Heat Transfer 88C, 101 共1965兲. 2 L. M. Huang and M. S. Elgenk, “Heat-transfer of an impinging jet on a flat surface,” ASME Trans. J. Heat Transfer 41, 1915 共1994兲. 3 D. Lacanette, A. Gosset, S. Vincent, J. M. Buchlin, and E. Arquis, “Macroscopic analysis of gas-jet wiping: Numerical simulation and experimental approach,” Phys. Fluids 18, 042103 共2006兲. 4 A. D. Friedman and J. Katz, “The flow and mixing mechanisms caused by the impingement of immiscible interface with a vertical jet,” Phys. Fluids 11, 2598 共1999兲. 5 J. Zheng, “An experimental study on a submerged water jet impinging upon a free air-water interface,” Ph.D. thesis, Kobe University, 2003. 6 A. Maurel, S. Cremer, and P. Jenffer, “Experimental study of a submerged fountain,” Europhys. Lett. 39, 503 共1997兲. 7 M. Gharib and A. Weigand, “Experimental studies of vortex disconnection and connection at the free surface,” J. Fluid Mech. 321, 59 共1996兲. 8 M. Brocchini and D. H. Peregrine, “The dynamics of strong turbulence at free surfaces. I. Description,” J. Fluid Mech. 449, 225 共2001兲. 9 S. S. Shy, “Mixing dynamics of jet interaction with a sharp density interface,” Exp. Therm. Fluid Sci. 10, 355 共1995兲. 10 H. Madarame, K. Okamoto, and M. Iida, “Self-induced sloshing caused by an upward round jet impinging on the free surface,” J. Fluids Struct. 16, 417 共2002兲. 11 S. Houard, F. Daviaud, and P. Bergé, “Phase-locking modes in a bidimensional network of coupled water jets,” Physica D 99, 318 共1996兲. 12 G. Bouchet, E. Climent, and A. Maurel, “Instability of a confined jet impinging on a water/air free surface,” Europhys. Lett. 59, 827 共2002兲. 13 T. H. Weisgraber and D. Liepmann, “Turbulent structure during transition to self-similarity in a round jet,” Exp. Fluids 24, 210 共1998兲. 14 L. Aanen, A. Telesca, and J. Westerweel, “Measurements of turbulent mixing using PIV and LIF,” Mach. Graphics Vision 8, 529 共1999兲. 15 H. Hu, T. Saga, T. Kobayashi, and N. Taniguchi, “Research on the vortical and turbulent structures in the lobed jet flow using LIF and PIV techniques,” Meas. Sci. Technol. 11, 698 共2000兲. 16 C. Fukushima, L. Aanen, and J. Westerweel, “Investigation of the mixing process in an axisymmetric turbulent jet using PIV and LIF,” Proceedings of the 10th International Symposium of Laser Techniques to Fluid Mechanics, Lisbon, Portugal, July 2000 共Springer, Berlin, 2000兲. 17 A. W. K. Law and H. Wang, “Measurements of mixing processes with combined DPIV and PLIF,” Exp. Therm. Fluid Sci. 22, 213 共2000兲. 18 A. Bunner and G. Tryggvason, “Effect of bubble deformation on the properties of bubbly flows,” J. Fluid Mech. 495, 77 共2003兲. 19 P. Sagaut, Large Eddy Simulation for Incompressible Flows—An Introduction 共Springer, Berlin, 1998兲. 20 P. Liovic and D. Lakehal, “Multi-physics treatment in the vicinity of arbitrarily deformable gas-liquid interfaces,” J. Comput. Phys. 222, 504 共2007兲. 21 P. Liovic and D. Lakehal, “Interface-turbulence interactions in large-scale bubbling processes,” Int. J. Heat Fluid Flow 28, 127 共2007兲. 22 E. Labourasse, D. Lacanette, A. Toutant, S. Vincent, P. Lubin, O. Lebaigue, J. P. Caltagirone, and P. Sagaut, “Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests, Int. J. Multiphase Flow 33, 1 共2007兲. 23 S. Vincent, J. Larocque, D. Lacanette, A. Toutant, P. Lubin, and P. Sagaut, “Numerical simulation of phase separation and a priori two-phase LES filtering,” Comput. Fluids 37, 898 共2008兲. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp 065110-21 24 Macroscopic analysis of a turbulent round liquid jet I. Wygnanski and H. E. Fiedler, “Some measurements in the selfpreserving jet,” J. Fluid Mech. 38, 577 共1969兲. 25 N. R. Panchapakesan and J. L. Lumley, “Turbulence measurements in axisymmetric jets of air and helium. I. Air jet,” J. Fluid Mech. 246, 197 共1993兲. 26 C. A. Friehe, C. W. Van Atta, and C. H. Gibson, “Jet turbulence: Dissipation rate measurements and correlations,” Proceedings in AGARD Turbulent Shear Flows, London, U.K., 1971 共AGARD, Paris, 1971兲, Current Paper No. 93. 27 J. Laufer, “The structure of turbulence in fully developed pipe flow,” National Advisory Committee on Aeronautics, Technical Note No. 2954, 1953. 28 M. Raffel, C. Willert, and J. Kompenhans, Particle Image Velocimetry 共Springer, Berlin, 1998兲. 29 I. Kataoka, “Local instant formulation of two-phase flow,” Int. J. Multiphase Flow 12, 745 共1986兲. 30 J. Magnaudet and I. Eames, “The motion of high-Reynolds number bubbles in homogeneous flow,” Annu. Rev. Fluid Mech. 32, 659 共2000兲. 31 R. Scardovelli and S. Zaleski, “Direct numerical simulation of free-surface and interfacial flow,” Annu. Rev. Fluid Mech. 31, 567 共1999兲. 32 J. U. Brackbill, B. D. Kothe, and C. Zemach, “A continuum method for modelling surface tension,” J. Comput. Phys. 100, 335 共1992兲. 33 J. Chesnel, J. Reveillon, F. X. Demoulin, and T. Menard, Subgrid modeling of liquid atomization,” Proceedings of the Sixth International Conference on Multiphase Flow 共ICMF兲, Leipzig, Germany, July 2007. 34 D. Lakehal, S. Reboux, and P. Liovic, “Subgrid scale modelling for the LES of interfacial gas-liquid flow,” Houille Blanche 6, 125 共2005兲. 35 S. Reboux, P. Sagaut, and D. Lakehal, “Large-eddy simulation of sheared interfacial flow,” Phys. Fluids 18, 105105 共2006兲. 36 P. Sagaut and R. Grohens, “Discrete filters for large eddy simulation,” Int. J. Numer. Methods Fluids 31, 1195 共1999兲. 37 M. Fortin and R. Glowinski, Méthodes de Lagrangien Augmenté. Application à la Résolution Numérique de Problèmes aux Limites 共Dunod, Paris, 1982兲. 38 S. Vincent, J. P. Caltagirone, P. Lubin, and T. N. Randrianarivelo, “An adaptative augmented Lagrangien method for three-dimensional multimaterial flows,” Comput. Fluids 33, 1273 共2004兲. 39 H. A. Van Der Vorst, “A fast and smoothly converging variant of BI-CG for the solution of non-symmetric linear systems,” SIAM 共Soc. Ind. Appl. Math.兲 J. Sci. Stat. Comput. 44, 631 共1992兲. 40 I. Gustafsson, “On first and second order symmetric factorization methods for the solution of elliptic difference equations,” Ph.D. thesis, Chalmers University of Technology, 1978. 41 S. Vincent and J. P. Caltagirone, “Efficient solving method for unsteady incompressible interfacial flow problems,” Int. J. Numer. Methods Fluids 30, 795 共1999兲. 42 S. Vincent and J. P. Caltagirone, “A one cell local multigrid method for solving unsteady incompressible multiphase flows,” J. Comput. Phys. 163, 172 共2000兲. 43 K. Khadra, S. Parneix, P. Angot, and J. P. Caltagirone, “Fictitious domain approach for numerical modelling of Navier-Stokes equation,” Int. J. Numer. Methods Fluids 34, 651 共2000兲. 44 T. N. Randrianarivelo, G. Pianet, S. Vincent, and J. P. Caltagirone, “Numerical modelling of solid particle motion using a new penalty method,” Int. J. Numer. Methods Fluids 47, 1245 共2005兲. 45 J. B. Ritz and J. P. Caltagirone, “A numerical continuous model for the hydrodynamics of fluid particles systems,” Int. J. Numer. Methods Fluids 30, 1067 共1999兲. 46 T. S. Lund, X. Wu, and K. D. Squires, “Generation of inflow data for spatially-developing boundary layer,” J. Comput. Phys. 140, 233 共1998兲. 47 A. Keating, U. Piomelli, and E. Balaras, “A priori and a posteriori tests of inflow conditions for large-eddy simulation,” Phys. Fluids 16, 4696 共2004兲. 48 M. Klein, A. Sadiki, and J. Janicka, “A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulation,” J. Comput. Phys. 186, 652 共2003兲. 49 N. Riviere, D. Reungoat, J. P. Faure, E. Biotteau, and J. P. Caltagirone, “Turbulent mixing and oscillations of a confined submerged jet impinging on a water/air free surface,” Proceedings of the 12th International Symposium on Flow Visualization, Göttingen, Germany, Sept. 2006, Paper No. ISFV12 8.5. 50 D. Reungoat, N. Riviere, and J. P. Faure, “3C PIV and PLIF measurements Phys. Fluids 21, 065110 共2009兲 in turbulent mixing: round jet impingement,” J. Visualization 10, 99 共2007兲. 51 I. Calmet and J. Magnaudet, “Large-eddy simulation of high-Schmidt number mass transfer in a turbulent channel flow,” Phys. Fluids 9, 438 共1997兲. 52 M. Klein, “Direct numerical simulation of a spatially developing water sheet at moderate Reynolds number,” Int. J. Heat Fluid Flow 26, 722 共2005兲. 53 J. G. H. Eggels, F. Unger, M. H. Weiss, J. Westerweel, R. J. Adrian, R. Friedrich, and F. T. M. Nieuwstadt, “Fully developed turbulent pipe: a comparison between direct numerical simulation and experiment,” J. Fluid Mech. 268, 175 共1994兲. 54 M. Rudman and H. M. Blackburn, “Direct numerical simulation of turbulent non-Newtonian flow using a spectral element method,” Appl. Math. Model. 30, 1229 共2006兲. 55 S. B. Pope, Turbulent Flows 共Cambridge University Press, Cambridge, 2000兲. 56 T. Djeridane, “Contribution à l’étude expérimentale de jets turbulents axisymétriques à densité variable,” Ph.D. thesis, Aix-Marseille University, 1996. 57 W. R. Quinn, “Upstream nozzle shaping effects on near field flow in round turbulent free jets,” Eur. J. Mech. B/Fluids 25, 279 共2006兲. 58 B. J. Boersma, G. Brethouwer, and F. T. M. Nieuwstadt, “A numerical investigation on the effect of the inflow conditions on the self-similar region of a round jet,” Phys. Fluids 10, 899 共1998兲. 59 M. Klein, A. Sadiki, and J. Janicka, “Investigation of the influence of the Reynolds number on a plane jet using direct numerical simulation,” Int. J. Heat Fluid Flow 24, 785 共2003兲. 60 S. Stanley, S. Sarkar, and J. P. Mellado, “A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation,” J. Fluid Mech. 450, 377 共2002兲. 61 P. C. Babu and K. Mahesh, “Upstream entrainment in numerical simulations of spatially evolving round jets,” Phys. Fluids 16, 3699 共2004兲. 62 M. Olsson and L. Fuchs, “Large eddy simulation of the proximal region of a spatially developing circular jet,” Phys. Fluids 8, 2125 共1996兲. 63 C. Bogey and C. Bailly, “Large eddy simulations of transitional round jets: Influence of Reynolds number on flow development and energy dissipation,” Phys. Fluids 18, 065101 共2006兲. 64 D. J. Bodony and S. K. Lele, “On using large-eddy simulation for the prediction of noise from cold and heated turbulent jets,” Phys. Fluids 17, 085103 共2005兲. 65 D. R. Webster, P. J. W. Roberts, and L. Ra’ad, “Simultaneous DPTV/PLIF meas of a turbulent jet,” Exp. Fluids 30, 65 共2001兲. 66 P. N. Papanicolaou and E. J. List, “Investigation of round vertical turbulent buoyant jets,” J. Fluid Mech. 195, 341 共1988兲. 67 W. K. George, in Advances in Turbulence, edited by W. K. George and R. Arndt 共Springer, New York, 1989兲. 68 F. Giralt, C. J. Chia, and O. Trass, “Characterization of the impingement region in an axisymmetric turbulent jet,” Ind. Eng. Chem. Fundam. 16, 21 共1977兲. 69 K. B. M. Q. Zaman and A. K. M. F. Hussain, “Vortex pairing in a circular jet under controlled excitation. I. General jet response,” J. Fluid Mech. 101, 449 共1980兲. 70 P. Brancher, J. M. Chomaz, and P. Huerre, “Direct numerical simulations of round jets: Vortex induction and side jets,” Phys. Fluids 6, 1768 共1994兲. 71 D. Liepmann and M. Gharib, “The role of the streamwise vorticity in the near-field entrainment of a round jet,” J. Fluid Mech. 245, 643 共1992兲. 72 Y. Na and P. Moin, “Direct numerical simulation of a separated turbulent boundary layer,” J. Fluid Mech. 374, 379 共1998兲. 73 J. Cohen and I. Wygnanski, “The evolution of instabilities in the axisymmetric jet. Part I. The linear growth of disturbances near the nozzle,” J. Fluid Mech. 176, 191 共1987兲. 74 D. T. Walker, “On the origin of the ‘surface current’ in turbulent free surface flows,” J. Fluid Mech. 339, 275 共1997兲. 75 W. L. Hong and D. T. Walker, “Reynolds-averaged equations for freesurface flows with application to high-Froude-number jet spreading,” J. Fluid Mech. 417, 183 共2000兲. 76 J. Piquet, Turbulent Flow: Models and Physics 共Springer, Berlin, 1999兲. 77 J. Westerweel, R. J. Adrian, J. G. M. Eggels, and F. T. M. Nieuwstadt, “Measurement with particle image velocimetry of fully developed turbulent pipe flow at low Reynolds numbers,” Applications Laser Techniques in Fluid Mechanics 共Springer, Berlin, 1993兲, pp. 47–99. Downloaded 31 Aug 2009 to 147.210.63.88. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp