Large-eddy simulation of a turbulent jet and wake vortex interaction

Large-eddy simulation of a turbulent jet and wake vortex interaction

Computers & Fluids 36 (2007) 772–785 www.elsevier.com/locate/compfluid Large-eddy simulation of a turbulent jet and wake vortex interaction O. Labbe´ ...

3MB Sizes 2 Downloads 64 Views

Computers & Fluids 36 (2007) 772–785 www.elsevier.com/locate/compfluid

Large-eddy simulation of a turbulent jet and wake vortex interaction O. Labbe´ *, E. Maglaras, F. Garnier ONERA, 29 Avenue de la Division Leclerc, 92322 Chaˆtillon cedex, France Received 16 August 2005; received in revised form 22 December 2005; accepted 19 June 2006 Available online 24 October 2006

Abstract In this study, temporal large-eddy simulations of the interaction between a turbulent jet flow and a trailing vortex are described. Three cases are analyzed: in the first one, the jet and the vortex axes are sufficiently well separated to not interact immediately, while in the second case, the distance between the jet and the vortex is reduced by half. In the last case the jet blows in the vortex core. In the two first cases, as the jet spreads it is progressively deflected by the continuous input of crossflow momentum. Thus it acquires azimuthal and radial components of velocity, causing the emergence of three-dimensional structures of azimuthal vorticity around it. When the jet and the vortex are superimposed, the turbulent kinetic energy does not increase, the vortex core is very buffeted. Numerical simulation results of the convection–diffusion of a passive scalar show that its distribution (initially in the jet) cannot penetrate inside the vortex core due to its solid rotation. For the cases where the jet is initially outside the vortex and when it is injected in the vortex wake, its value remains very high and cannot get out of the vortex core. This phenomenon confirms the existence of a stabilizing ‘‘dispersion buffer’’, adjacent to the core, which prevents amplification of the turbulence generated inside the core.  2006 Elsevier Ltd. All rights reserved.

1. Introduction In the near field of the aircraft wake, the exhaust jets of the turbines are entrained into two counter-rotating wingtip vortices, which at the same time roll up from the sheet of vorticity induced by the wings. To investigate the highly complex entrainment and mixing process of the turbulent jet into a trailing vortex, it is convenient to identify two overlapping regimes; the jet regime and the deflection regime. The qualitative features of jet/vortex interaction were illustrated by Miake-Lye et al. [23] who identified two distinct phases. Firstly, the jets rapidly mix with ambient air while the vorticity shed from the wings rolls up into a pair of trailing vortices. The dynamics are then dominated by the entrainment of the jets into the vortex flow in the interaction regime. In 1996, as part of a numerical investigation on the dynamics of a jet in a aircraft vortex wake, Jacquin and Garnier [18] presented a

*

Corresponding author. Tel.: +33 146 73 42 50; fax: +33 146 73 41 66. E-mail address: [email protected] (O. Labbe´).

0045-7930/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2006.06.001

smoke visualization photograph of a jet exhaust behind a representative large transport aircraft at cruise conditions. A laser lightsheet was placed in the crossflow plane at onehalf wingspan downstream. The visualization revealed a tear-drop-shaped plume, suggesting that part of the jet may have already been entrained into the wing vortex at a short distance downstream. In 1999 Brunet et al. [4] presented experimental results of a hot jet behind a generic wing planform NACA0012 equipped with two jets beneath the airfoil. Laser Doppler Velocimetry and thermocouples were used to survey the jet at four locations ranging from one-half to eight wingspans downstream. In 2002 Ferreira Gago et al. [9] used a three-dimensional temporal numerical simulation to study the interaction between a jet and a vortex wake in order to contribute to the understanding of engine jet dispersion of an aircraft in cruise condition [11,5]. After the computation of the jet development, a Lamb–Oseen vortex is imposed in the vicinity of the jet. Large-scale structures appear as counter-rotating vortex helical structures visualized by azimuthal vorticity. They rapidly lose coherence, while the vortex core keeps its axial vorticity.

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

Besides the studies mentioned previously, the computation of this complex flow has been undertaken by various authors. Gerz and Ehret [13], and Paoli et al. [26] simulated the jet/vortex interaction with large-eddy simulation (LES), while Fares et al. [6] used a Reynolds Averaged Navier– Stokes (RANS) solver. More recently Wang and Zaman [31] made an experimental investigation of jet behavior in the trailing vortex of a delta wing. Huppertz et al. [17] investigated experimentally and numerically the very near wake up to five chords behind the trailing edge of a rectangular wing and analyzed different engine positions and jet speed. Holza¨pfel et al. [15] analysed the wake vortex decay mechanisms. However, in spite of the progress made in the understanding of turbulence mechanisms in wake vortices, essential information still lacks, especially with the role of large-scale structures on a scalar field distribution. The latter helps in explaining the entrainment and mixing processes of the engine emissions in the vortex. This study is a natural continuation of the work performed by Ferreira Gago et al. [9]. However, the aim is neither to compare the results with the experimental data, nor to test different turbulence models, that has already been done in [10]. It is an attempt to give a physical insight into how the turbulent jet is wrapped around the vortex. In order to study the interaction of the jet turbulence on the vortex with respect to their initial relative locations, three different jet positions have been tested. The first one coincides with the experimental set-up [4] and in the second one, the jet is closer to the vortex; the distance between them is reduced by a factor of two. In the third case, the jet is located in the center of the vortex. The Reynolds number of flows of aeronautic interest is often very high, however many physical phenomena are also observed in flows with low Reynolds numbers. Temporal large-eddy simulations with a hybrid Mixed-Scale Model on a moderate refined mesh have been carried out for three positions of the jet. The computation is advanced in two steps; the first step is stopped when it provides an estimation of the end of the jet regime and in the second step, the interaction with a Lamb–Oseen vortex field is accounted for. The structure of the paper is as follows. In Section 2, the governing equations and subgrid modeling are given. Numerical methods and simulation parameters are given in Section 3. A comparison between the first two cases, and then the third one are described and analyzed in Section 4. Finally, in Section 5 the conclusions are presented. 2. Governing and passive scalar transport equations and subgrid modeling In the LES approach, the dimensionless formulation of compressible 3D Navier–Stokes equations are filtered spatially so that any variable f may be decomposed into a resolved (or large-scale) part f and a non resolved (or subgrid scale) part f 0 , with f ¼ f þ f 0 . This procedure may be obtained by a convolution integral of the variable

773

with a filter function depending on a filter width D. Practically, the filter width is simply given by the computational mesh cell size D = Diso, where Diso = (DxDyDz)1/3. For compressible flows, it is usual to introduce the Favre filter [8] defined as qf ¼ qf~ After some simplifications, we obtain the Favre-filtered equations of mass and momentum (reader can refer for example to [10] for details) oq o þ ðq~uj Þ ¼ 0 ot oxj o o op ob r ij o ðq~ui Þ þ ðq~ui ~uj Þ þ  ¼ sij ot oxj oxi oxj ozj

ð1Þ ð2Þ

The viscous stress tensor is given by   lð Te Þ o~ui o~uj 2 o~uk b r ij ¼ þ  dij Re oxj oxi 3 oxk where dij being the Kronecker delta and Re the Reynolds number. The right-hand sides of Eq. (2) contain the so-called subgrid stress tensor sij ¼ qðg ui uj  u~i u~j Þ, which represent the effects of unresolved scales of motion on the large-scales and need to be modeled. In this study we use an hybrid Mixed-Scale Model which is a linear combination of an eddy viscosity-type model and a similarity model. Hybrid models take advantage of both the similar and eddy viscosity-type models. They correlate well with the real subgrid scale stresses and adequately dissipate the energy in the small-scales. The subgrid scale viscosity is provided by the Mixed-Scale Model proposed by Sagaut [28]. The similarity model, formulated by Bardina et al. [1] and revised by Liu et al. [21] is not of the eddy viscosity type and reproduces the details of the stresses more accurately. After commonly adopted simplifications, the filtered equations for the energy equation has the following form: b oE o ob qj b þ pÞ~uj g  o ðb þ fð E r ij ~ui Þ þ ot oxj oxj oxj ¼ B1  B2 þ B3 þ B4 b ¼ p þ 1 qu~j u~j ; the The computable energy is given by E 2 c1 lðe TÞ oe T computable heat flux vector by b q ¼ . j

ðc1ÞRePrM 2 oxj

The Prandtl number Pr and the ratio of specific heats c = Cp/Cv are set equal to 0.7 and 1.4, respectively. Furthermore, our simulation is performed by using a low Mach number M of 0.2, which should make the calculation close to the incompressible limit. B1 ¼

1 o ðpuj  p~uj Þ; c  1 oxj

B3 ¼ skj

o ~uk ; oxj

B4 ¼ rkj

B2 ¼

o ðskj ~uk Þ oxj

o o ~uk uk  rkj oxj oxj

In this study we adopt an hybrid model for the subgrid term B1

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

774

1 o B1 ¼ 2 oxj

( 

qmsm o Te 2 ðc  1ÞPrt M oxj

1 c c ~ þ q Te ~u j  q Te qb uj ðc  1ÞcM 2

, ) b q

The SGS terms B2 and B3 depend directly upon the subgrid-scale tensor, they are thus obtained in a straightforward manner reporting the expression used to model the former tensor. Following Ghosal et al. [14] the subgrid 3=2 term B4 is modeled by B4 ¼ Cq k D . The procedure to determine the coefficient C is based on a global balance of the terms in the integrated k-equation, which leads to R ðð1  cÞB1 þ B2  B3 ÞdX m2 C¼ X with k ¼ sm2 R k3=2 D q D dX X The filtered passive scalar convection–diffusion equation (Z is initially zero in the external flow and reaches a maximum of one at the jet center) takes the following form: ( ) e o e o 1 o o Z oC 1 e~ ðq Z Þ þ ðq Z uj Þ  lð Te Þ ¼ ot oxj ReSc oxj oxj oxj

the number of discretization points, which leads to a time-developing simulation. Unfortunately, this condition is not very physical for jet simulations, because the mean flow gradients in the axial direction are neglected. The instabilities in the simulated flow are of convective nature and absolute instabilities cannot be captured. Brancher [3] showed that the Kelvin–Helmoltz instability developed in jet flows has a convective nature and the temporal approach has been validated for academic flows, such as mixing layers and jets [16,30]. 4. Computational parameters The computations are split in two phases: the jet phase (Section 4.1) and the interaction phase (Section 4.2). The equations presented in Section 2 are non-dimensionalized by scaling the velocities with the centerline velocity of the jet Vj, and the characteristic length scale is equal to the radius R. The Reynolds number is based on jet velocity and jet radius, Re = VjR/m = 5000 and the Mach number is 0.2, the Prandtl and Schmidt numbers are Pr = 0.7 and Sc = 0.7. 4.1. Jet phase

where Sc is the Schmidt number and C1 is defined by fj  Z e~ C 1 ¼ qð Zu uj Þ. The hybrid model for this equation leads to: ( , ) 1 oZ c c e~ e qb ~ qmsm C1 ¼ þ qZ uj  qZ uj b q 2 oxj where the subgrid scale viscosity is given by the MixedScale Model. 3. Numerical method and boundary conditions In convective terms, spatial derivatives are taken with sixth-order compact finite differences [20]. To minimize the aliasing error, we follow the procedure applied by Boersma and Lele [2]. The non-linear terms have been rewritten in the skew symmetric form i.e.,   oqui uj 1 oqui uj oquj oui ¼ þ ui þ quj 2 oxj oxj oxj oxj Diffusive terms are discretized by using a second-order accurate finite difference scheme. The time integration is carried out by a third-order Runge–Kutta method [22]. When the net circulation of the flow is non-zero, Pradeep and Hussain [27] demonstrated that the periodic boundary conditions typically applied in numerical simulations of an isolated vortex dynamics in transverse directions, can lead to significantly incorrect results. In the present study, an isolated vortex is simulated and non-reflexive conditions introduced by Thompson [29] are used at the lateral boundaries of the computational domain. Moreover, in the streamwise direction (y), periodic conditions are applied on boundaries in order to reduce

In the first step only the jet development is computed, not accounting for the vortex field. The computational domain is rectangular and the grid is uniform in all directions. The cross plane extends from x, z = 5.1R to x, z = 5.1R. In order to observe the development of an instability, the length of the computational domain must be an even multiple of the instability’s wavelength. Michalke and Hermann [24] give the wavelength of maximum growth of the first azimuthal disturbance. Furthermore, Gaster [12] establishes a relation between the frequency and amplification rates for a disturbance

Fig. 1. Evolution of the turbulent kinetic energy of the jet development.

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

growing with respect to time and those of a spatially growing wave having the same wave number. So by using Gaster’s transformation, we obtain the wavelength of maximal growth of the first azimuthal disturbance for a temporally evolving jet. The streamwise length of the computational domain is Ly = 6R. This corresponds to twice the wavelength of the maximum growth rate of the first azimuthal instability of a spatially evolving jet. The LES is carried out on a (69 · 61 · 69) grid corresponding to a mesh of Dy = 0.1R and Dx,z = 0.15R, which is probably a coarse grid, but enough to obtain a fully turbulent region. The Courant–Friedrichs–Lewy number is 0.6 leading to a time step Dt = 0.01. The jet axial velocity is initialized as V ¼ V 0 ð1 þ e Ve Þ where V0 is a tanh profile, r being the radial coordinate in a cross-section, and e Ve being a white noise random perturbation with a maximum amplitude 1% of V0 and k Ve k 6 1.

775

In the following, we use the  initial velocity profile: V The momentum V 0 ðrÞ ¼ 2j ð1  tanh 14 HR Rr  RrÞ . boundary layer thickness H is defined as H¼H¼ o R 1 nV 0 ðrÞon V 0 ðrÞ 1  V j dr. We restrict the present investiga0 Vj tion to one value of the jet parameter, namely R/H = 10. Among the cases studied by Michalke and Hermann [24], this value corresponds to the most unstable jet velocity profile. The initial mean temperature is provided by the Crocco–Busemann relation and the jet temperature ratio is Tj/T1 = 2, Tj, T1 being respectively the jet and ambient temperatures. The end of the jet regime is chosen when the maximum of turbulent kinetic energy is reached, around an dimensionless time t = t*/tj = 50, tj being R/Vj. This time corresponds to y/b = 0.5, b being the wingspan of the experimental model, the first cross section measured

Fig. 2. Left: transverse vorticity contour lines through the jet, solid/dashed lines indicate positive/negative vorticity Xz 2 [3, 3]. Right: passive scalar iso-surface Z = 0.4, top: t = 20, bottom t = 50.

776

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

downstream in the experiments [4]. At this time, the temperature jet has reduced by half. Fig. 1 shows the evolution of the turbulent kinetic energy and is calculated as follows. First a spectrum is acquired by taking the Fourier transform of the velocity field in the periodic direction and integrating the values in the (x, z) plane. The turbulent kinetic energy is obtained by adding the values of the one-dimensional spectrum for all wave numbers. In a first phase, this energy increases linearly, it reaches a peak around t = 50 and then decreases. The vorticity iso-contours Xz are shown in Fig. 2, where the transition to turbulence of the

jet, associated to the formation of structure with high vorticity which causes the roll-up of the jet is clearly observed between t = 20 and t = 50. On the right hand side the isosurface of passive scalar Z = 0.4 gives the state of the jet when the vortex is added. Another important jet characterization is given with the turbulent intensities plotted in Fig. 3. In the case of round jet flow simulations, profiles of turbulent stresses are formulated in cylindrical coordinates and are computed as follows. First, the instantaneous field /i (velocity) is interpolated into a polar grid. Secondly, the Fourier transform of each velocity component

Fig. 3. Profiles of mean normalized turbulent intensities hu02h i=hU c i2 ; 2 2 02 hu02 r i=hU c i ; hU y i=hU c i , at t = 50. Fig. 5. Two-dimensional sketch of computational domain for the entrainment case.

Fig. 4. One-dimensional spectra for the computational velocity field at t = 50.

Fig. 6. Lamb–Oseen vortex diffusion: comparison between the analytical solution and the numerical solution.

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

is computed. The turbulent intensities are defined as b i ðx; z; k y ; tÞ / b t ðx; z; k y ; tÞ. Here / b i; / b t and ky denote j/0i /0i j ¼ / i j respectively the Fourier transform of the velocity component, the complex conjugate of the Fourier transform velocity component and the wave number in the axial direction. Finally h/0i /0i i is computed by taking a simple average over all points in the homogeneous direction. The mean turbulent intensities (azimuthal, radial and axial) are normalized by the axial velocity in the jet center (Uc) at time t = 50. The level of turbulent intensities azimuthal and radial velocity components are of same level (7%) but the axial component is lower and more spread out. Fig. 4 presents one-dimensional energy spectrum E(ky), computed at t = 50. It is obtained by taking the Fourier transform of the velocity field in the axial direction and integrating the Fourier coefficients in the x, z plane. At high wave numbers, the energy is low and indicates that the flow is well resolved.

777

The comparison shows that profiles obtained numerically are in agreement with the analytical solution. 5. Results and discussion In order to study the interaction of the jet turbulence on the vortex with respect to the initial relative locations, three different jet positions have been tested. In the two first

4.2. Interaction phase The interaction with the vortex wake is then taken into account, and a Lamb–Oseen vortex, is fitted on the transverse part of the experimental velocity field given by 2

U h ðrÞ ¼ aU h max ðrc =rÞð1  expðbðr=rc Þ ÞÞ and oq=or ¼ qU 2h =r where a = 1/(1  eb) and b = 1.2564. The computational domain is enlarged to 169 · 169 nodes, in the cross plane. It is uniform in the domain 11.1R < x, z < 11.1R, with a mesh of 0.15R and stretched using a monotonic tanh law outside this domain. A scheme of the (x, z) flow configuration is presented in Fig. 5, with the jet position corresponding to the experimental set-up. The reference velocity Uh max corresponds to the viscous radius rc, with rc/R = .7 at the beginning of the interaction phase the initial jet to vortex velocity ratio is Vj/Uh max = 1. The dimensionless time tv is introduced as tv = t*/(R/ Uh max) and its origin is reset to zero. For computing the mean flow profiles, the simulated data, which are in Cartesian coordinates (x, y, z), are transformed to cylindrical coordinates (r, h, y) using a high order cubic spline interpolation formula. All the results presented as a function of the radius r/R are obtained in the following way; the instantaneous quantity is averaged in the axial direction y. Then the averaged quantity is interpolated into a polar grid and averaged in the azimuthal direction to obtain a profile which is a function of r. No averaging in time is performed; therefore, the mean profile is a function of time. In order to show that the simulation code diffuses correctly a laminar vortex, we have compared the numerical solution of a Lamb–Oseen vortex with the analytical one. In Fig. 6, the diffusion of the Lamb–Oseen vortex is plotted with the tangential velocity at times tv = 50 and 100. The full line corresponds to the initial Lamb–Oseen vortex.

Fig. 7. Longitudinal component of vorticity Xy, (1, 4.4; D = .2) at tv = 0 in a median plane, dashed/solid lines correspond to negative/positive values.

778

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

cases, the jet is located outside the vortex, and in the last case the jet is superimposed to the vortex core and will

be studied separately. The first one, case A, coincides with the experimental set-up described above, in which

Fig. 8. Evolution of the turbulent kinetic energy: case A (left) and case B (right).

Fig. 9. Three-dimensional view of the azimuthal vorticity Xh (green: 0.4, orange: 0.4) at tv = 55,155 and 250 for case A (left) and tv = 3, 25 and 150 for case B (right); the position of the vortex is plotted with an iso-surface of vorticity Xy (+2) in black. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

the distance between the jet and the vortex center at time tv = 0 is equal to 7.2R. In the second one, case B, the jet is closer to the vortex and the distance of experimental set-up is reduced by half i.e. 3.6R. In this case, the jet turbulence and axial velocity levels are higher than in case A. The interaction is stronger and the jet could be more disturbed. The last case C is quite different because the jet is blowing in the vortex core making the flow similar to a q-vortex (Batchelor). To illustrate the relative jet position with respect to the vortex, Fig. 7 shows the longitudinal component of vorticity Xy at time tv = 0, in the median plane normal to the axial direction. The vortex core is concentrated at the center of the computational domain (i.e. x/R = 0 and

779

z/R = 0), positive and concentric contours represent the vorticity distribution of the vortex core, while the jet patterns consist in positive and negative coherent fields. The analysis of the results is performed in terms of the turbulent kinetic energy, the azimuthal and longitudinal vorticity components, the tangential velocity for the dynamic process, the passive scalar and the scalar fluxes for the mixing characteristics. In this paragraph both cases A and B are analyzed and compared in order to evaluate the jet/vortex interaction dynamics and also the jet mixing characteristics. The evolution of the kinetic energy is shown in Fig. 8 for the two configurations. The two curves present a different behavior in the interaction process. For case A, the energy first

Fig. 10. Longitudinal component of vorticity Xy (1, 4.4; D = .2) at tv = 55, 155 and 250 for case A (left) and tv = 3, 25 and 150 for case B (right); in a median plane, dashed/solid lines correspond to negative/positive values of Xy.

780

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

Fig. 11. The tangential velocity for cases A and B at the different characteristic times.

decreases due to the jet diffusion process, such as if the jet was alone. The jet is then sufficiently close (tv = 55) to the vortex to interact with the latter and, due to amplification of the instability, the energy increases to a maximum 0.19 at tv = 155 which corresponds to an energy increase of 0.11. The energy then slowly decreases. For case B, the jet is sufficiently near the vortex to interact immediately (tv = 3) and the turbulent energy saturates. A peak (0.3) is reached at tv = 25 and the increase of energy 0.11 is the same as in case A. The energy decays then quickly. To illustrate the dynamics phenomena, for each case three characteristic times were chosen corresponding to growth, saturation and decay of the turbulent kinetic energy. At the beginning of the interaction (case A: tv = 55; case B: tv = 3), the turbulent kinetic energy grows to reach a maximum (case A: tv = 155; case B: tv = 25) followed by a period of decay called the dissipation regime (case A: tv = 250; case B: tv = 150). In Fig. 9, the helical structures surrounding the core are surfaces of azimuthal vorticity for the times mentioned

Fig. 12. Evolution of the passive scalar: case A (left) and case B (right).

Fig. 13. Turbulent flux hZ 0 u02 i: case A at tv = 155 (left) and case B at tv = 25 (right).

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

previously. These large-scale helical structures are counterrotating vortices where the green color surface indicates a positive value (+0.4) and the orange color a negative one (0.4). The position of the vortex is plotted with an isosurface of vorticity Xy(+2) in black. The dynamics of the interaction are first dominated by the entrainment and the captation of the jet by the vortex. As the jet spreads it is progressively deflected by the continuous input of crossflow momentum so that it acquires azimuthal and radial components of velocity. The deflection regime corresponds to the emergence of three-dimensional structures of azimuthal vorticity around it. The vorticity of the jet is progressively stretched and generates spiral-shaped vorticity structures. At the turbulent kinetic energy maximum, the large-scale vortical structures seem more coherent in case A than in case B. In the dissipation regime, the jet is more rapidly annihilated when it is injected closer to the vortex. In case B, the jet interacts immediately with the vortex due to the short distance between them. The jet axial velocity is higher than in case A because this velocity component decreases with time as commonly observed in the literature. As the increase of turbulent kinetic energy is the same for both cases, this augmentation seems more related to the vortex than the jet itself. As observed in the figure, the break-up of these large-scale structures takes place, however in both cases the jet does not perturb definitively the vortex core. In order to study the influence of the jet on the vortex, the longitudinal component of vorticity Xy has been plotted in Fig. 10 at time tv = 55, 155 and 250 for case A and at time tv = 3, 25 and 150 for case B in a median plane normal to the axial direction. The fully turbulent jet is deflected and entrained by the vortex velocity field and wraps around it. For case A, the vortex core does not seem very perturbed by the jet and at tv = 250, the overall vortex has almost recovered its initial state. Case B (bottom) is more noticeable because at tv = 25 the shape of the vortex is greatly altered by the jet, the negative iso-lines remaining outside the vortex core. In the diffusion zone (tv = 150) the

781

jet has been completely eliminated and the vortex core retrieves its initial shape. An enlargement of the tangential velocity profiles for the times described above, is depicted in Fig. 11. In both cases, the tangential velocity is decreasing and it is noticeable that at the same time t = 155 (case A) and t = 150 (case B) the tangential velocity is lower in case B than in case A. It is also interesting to study the influence of jet turbulence on the vortex through the passive scalar. To quantify the mixing process at different times, referenced previously, the computational passive scalar Z as a function of the radius r is plotted in Fig. 12. Results of the numerical simulation of the convection–diffusion of a passive scalar show that its distribution initially immersed in the jet, is very spread out and remains constant but weak for case A (Fig. 12 Left). In case B the scalar passive is closer to the vortex but cannot enter into the core. In addition to the investigation of the scalar field distribution, the turbulent

Fig. 15. Evolution of the turbulent kinetic energy of case C.

Fig. 14. Variance of scalar fluctuations hZ 0 Z 0 i: case A at tv = 155 (left) and case B at tv = 25 (right).

782

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

flux hZ 0 u02 i and the variance of scalar fluctuations hZ 0 Z 0 i calculated as explained previously for the jet, are displayed in the following figures. The three components (axial, azimuthal and radial) hZ 0 u02 i; hZ 0 u0h i and hZ 0 u0r i of the turbulent flux are of the same level. Moreover the structure of iso-lines of turbulent flux hZ 0 u02 i in Fig. 13 for both jet positions at their respective maximum turbulent kinetic energy, at tv = 155 for case A and at tv = 25 for case B, are correlated to the large helical structures described in Fig. 9. For case B, where the turbulence is higher, no coherent structures are visible, neither in turbulent flux, nor in azimuthal vorticity iso-surfaces. Fig. 14 presents the variance of scalar fluctuations hZ 0 Z 0 i at the same time as previously. The

Fig. 16. Case C: three-dimensional view of the azimuthal vorticity Xh (green: 0.4, orange: 0.4) at tv = 4 (top), 25 (medium) and 200 (bottom); the vortex position is plotted with an iso-surface of vorticity Xy (+2) in black. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

both cases show that the maximum hZ 0 Z 0 i is located in the helical structure which rolls up around the vortex core. The maximum hZ 0 Z 0 i in case B is ten times higher than in case A because at the maximum interaction the jet is still very turbulent. Finally, the results of both these cases where the jet is initially located outside the vortex, reveal a similar behavior. The fully turbulent jet is deflected and entrained by the vortex velocity field and starts to wrap around it. During that process the jet vorticity is progressively stretched and leads to the production of large structures in the form of

Fig. 17. Case C: longitudinal component of vorticity Xy (1, 4.4; D = .2) at tv = 4 (top), 25 (medium) and 200 (bottom), in a median plane, dashed/ solid lines correspond to negative/positive values of Xy.

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

783

helical sheets of vorticity. Due to the strong rigid-body-like flow field, the closer the jet is to the vortex, these structures break up and drive the vortex to a laminar configuration. In case C, the jet is superimposed onto the Lamb–Oseen vortex yielding an interaction strong enough to lead to a destabilization of the vortex core. As previously, the jet turbulence is already developed and injected in the vortex core. This case is similar to a q-vortex configuration, but in this configuration the jet is not reduced to an axial velocity. Fig. 15 presents the evolution of turbulent kinetic energy. After a very short increase, the curve exhibits a peak (0.23) at tv = 4, then decreases with oscillations. In order to analyze the different phenomena, three characteristic times have been chosen. The first characteristic time coincides with the maximum turbulent kinetic energy tv = 4, the second tv = 25 when the energy is reduced by half and the third tv = 200 when the energy has reached a minimum level (0.05). Similarly, as shown in Fig. 9 for cases A and B, contours of azimuthal and streamwise vorticity components are presented in Fig. 16 for the three times previously defined with the same contour levels. At time tv = 4 (top), the dynamics of the interaction are dominated by both positive and negative vorticity embedded within each other and very concentrated around the core. They progressively expand at time tv = 25 (medium), then the flow rotation suppresses the small-scale motions of turbulence, at time tv = 200 (bottom) the structures have almost disappeared. As previously, the interaction of the jet with the vortex is analyzed through the longitudinal component of vorticity Xy, see Fig. 17 at the three characteristic times in a median plane normal to the axial direction. Examining the isocontours, the shape of the vortex is very altered by the jet. At time tv = 4 (top), the longitudinal component of

vorticity Xy is concentrated. The negative vorticity (dashed lines) are already located outside of the core. At the final stage tv = 200 (bottom) the jet seems to be eliminated and the vortex core slowly retrieves its initial shape with spurious oscillations. At the intermediate time tv = 25 (medium) the iso-contours are more spread out than at tv = 4 and the vortex core starts to recover a normal shape. A time animation of the results shows that the vortex core is strongly buffeted by the jet and oscillates in the cross plane where the displacement core can reach 0.3R. This phenomenon is due to the meandering. The evolution of the longitudinal velocity in the center of the jet/vortex in

Fig. 18. Evolution of the longitudinal velocity at the center (X/R = 0, Z/R = 0).

Fig. 20. Evolution of the passive scalar at the different characteristic times in case C.

Fig. 19. Evolution of the circulation at the different characteristic times in case C.

784

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

Fig. 18 shows that two phenomena are superposed, the first one which is high frequency corresponds to the meandering of the vortex and the second one, low frequency, coincides with the oscillations of the turbulent kinetic energy of Fig. 15. Paoli et al. [26] do not observe the same behavior and find a typical vortex trajectory (descent) corresponding to a symmetry boundary condition used in the numerical simulations. Case C corresponds to a Batchelor q-vortex. The stability of this vortex is controlled by the value of the swirl parameter S, defined as S = 1.56 max(uh)/max(u2) which measures the relative tangential and axial velocity intensity. The critical value of S, given for example by Fabre and Jacquin [7] is approximately 1.5. For larger values, rotation

stabilizes all perturbations, and these latter are transformed into a regime of inertial waves. In case C, as uh is always greater than u2, the swirl parameter is in this case greater than 1.5 and therefore the instability saturates and a laminarization of the flow is reached. In Paoli et al. [26], the Swirl number is low, which leads to the deformation of the vortex structure. The circulation profiles is plotted in Fig. 19 versus r/R at the initial time and the three characteristic times. It shows the emergence of a circulation overshoot of small magnitude at the interaction peak tv = 4 also obtained from DNS of Pantano and Jacquin [25]. The propagation of the circulation overshoot is the diffusion process by which angular momentum may be finally transported outward, albeit at a small rate. In order to study the mixing process, Fig. 20 displays the computational passive scalar Z with respect to the radius r at the different referenced times. When the jet is injected in the vortex, the maximum increases with time and the passive scalar is trapped in the vortex core. The turbulent flux hZ 0 u02 i and the variance of scalar fluctuations hZ 0 Z 0 i at tv = 4 are plotted in Fig. 21, where the high values are concentrated around the vortex core. This is may be related to the existence of a ‘‘dispersion buffer’’ region discussed by Jacquin and Pantano [19], adjacent to the core, which prevents perturbations and turbulence generated inside the core from being amplified when reaching the vortex periphery. Compared to case B, the peak level is much higher. Even injected in the vortex, the jet does not destabilize the vortex. Three-dimensional structures progressively expand around the vortex and lose coherence. The vortex is greatly buffeted and undulations appear, the passive scalar seems trapped in the core. 6. Conclusion

Fig. 21. Turbulent flux (top) hZ 0 u02 i and variance of scalar fluctuations (bottom) hZ 0 Z 0 i at tv = 4 in case C.

Large-eddy simulations of the interaction between a turbulent jet and a wake vortex have been performed for the three configurations. Initially, the jet simulation is carried out until the maximum turbulent kinetic is reached. Then, the resulting jet data is injected outside of the vortex in cases A and B, while in the last case C, the jet is blown into the vortex core. According to the cases, the jet is located at a distance of 7.2R (case A), 3.6R (case B) or in the Lamb– Oseen center (case C). For cases A and B where the jet interacts with the vortex, the turbulent kinetic energy increases and due to the longitudinal jet velocity, largescale structures appear around the vortex core as counter-rotating vortex helical structures. As the downstream distance increases, the large-scale vortical structures disappear and the kinetic energy decays. The vortical structures rapidly lose coherence and the stabilizing effect of the flow rotation suppresses small-scale motions of turbulence. During all the simulations the vortex core keeps its positive axial vorticity. In case C, the jet is superimposed onto the vortex and the longitudinal vorticity shows the meandering of the vortex. For the three jet positions simulated in the present study, the vortex recovers its initial shape.

O. Labbe´ et al. / Computers & Fluids 36 (2007) 772–785

The evolution of the passive scalar clearly shows that the distribution initially included in the jet decreases rapidly and does not enter the vortex core for cases A and B. When the jet is injected in the vortex core, the passive scalar is concentrated in it and seems not be able to escape of this region. That is conform to the existence of a region adjacent to the core called the ‘‘dispersion buffer’’ discussed by Jacquin and Pantano [19]. In case C, the value of the Swirl number is sufficiently high to allow the vortex to evolve toward an equilibrium state and to become persistent. References [1] Bardina J, Ferziger JH, Reynolds WC. Improved subgrid scale models for large-eddy simulation. AIAA Paper 80-1357. [2] Boersma BJ, Lele SK. Large eddy simulation of a Mach 0.9 turbulent jet. AIAA Paper 99-1874. [3] Brancher P. Etude nume´rique des instabilite´s secondaires de jets. PhD Thesis, Ecole Polytechnique, Palaiseau, 1996. [4] Brunet S, Garnier F, Jacquin L. Numerical/experimental simulation of exhaust jet mixing in wake vortex. In: AIAA 30th fluid dynamics conference Norfolk (USA), June 28–July 01, 1999. [5] Brunet S, Jacquin L, Geffroy P. Experiment on heated jets/wake vortex interaction. ONERA RT 15/2496 DAFE/Y, 1999. [6] Fares E, Meinke M, Schroeder W. Numerical simulation of the interaction of wingtip vortices and engine jets in the near field. AIAA Paper 2000-2222. [7] Fabre D, Jacquin L. Viscous instabilities in trailing vortices at large swirl numbers. J Fluid Mech 2004;500:239–62. [8] Favre A. Equations statistiques aux fluctuations turbulentes dans les e´coulements compressibles: cas des vitesses et des tempe´ratures. CR Acad Sci Paris A 1971;273:1087–93. [9] Ferreira Gago C, Brunet S, Garnier F. Numerical investigation of turbulent mixing in a jet/wake vortex interaction. AIAA J 2002;40(2): 276–84. [10] Ferreira Gago C, Garnier F, Utheza F. Direct testing of subgrig scale models in large-eddy simulation of a non-isothermal turbulent jet. Int J Numer Meth Fluids 2003;42:999–1026. [11] Garnier F, Brunet S, Jacquin L. Modelling exhaust plume mixing in the near field of an aircraft. Ann Geophys 1997;15:1468–77. [12] Gaster M. A note on the relation between temporally-increasing and spatially-increasing disturbances in hydrodynamic stability. J Fluid Mech 1962;286:229–55.

785

[13] Gerz T, Ehret T. Wingtip vortices and exhaust jets during the jet regime of aircraft wakes. Aerospace Sci Technol 1997;14:222–4. [14] Ghosal S, Lund TS, Moin P, Akselvoll K. A dynamic localization model for large-eddy simulation of turbulent flows. J Fluid Mech 1995;286:229–55. [15] Holza¨pfel F, Hofbauer T, Darracq D, Moet H, Gamier F, FerreiraGago C. Analysis of wake vortex decay mechanisms in the atmosphere. Aerospace Sci Technol 2003;7:263–75. [16] Hu G, Sun D, Yin X. A numerical study of dynamics of a temporally evolving swirling jet. Phys Fluids 2001;13:951–65. [17] Huppertz G, Fares E, Abstiens R, Schroder W. Investigation of engine jet/wing-tip vortex interference. Aerospace Sci Technol 2004;8:175–83. [18] Jacquin L, Garnier F. On the dynamics of engine jets behind a transport aircraft, in the characterisation and modification of wakes from lifting vehicles in fluids. AGARD CP 1996;584:37.1–8. [19] Jacquin L, Pantano C. On the persistence of trailing vortices. J Fluid Mech 2002;471:159–68. [20] Lele SK. Compact finite difference scheme with spectral-like resolution. J Comput Phys 1992;103:16–42. [21] Liu S, Meneveau C, Katz J. On the properties of similarity subgridscale models as deduced from measurements in a turbulent jet. J Fluid Mech 1994;275:83–119. [22] Lowery PS, Reynolds WC. Numerical simulation of a spatially developping forced mixing layer, Rep TF-26, Stanford University, 1986. [23] Miake-Lye RC, Martinez-Sanchez M, Brown RC, Kolb CE. Plume and wake dynamics, mixing and chemistry behind a high speed civil transport. J Aircraft 1993;30:467–79. [24] Michalke A, Hermann G. On the inviscid instability of a circular jet with external flow. J Fluid Mech 1982;114:343–59. [25] Pantano C, Jacquin L. Differential rotation effects on a turbulent batchelor vortex. In: DLES-IV workshop, Enschede, The Netherlands, July 18–20, Kluwer, 2001. [26] Paoli R, Laporte F, Cuenot B. Dynamics and mixing in jet/vortex interactions. Phys Fluids 2003;15(7):1843–60. [27] Pradeep DS, Hussain F. Effects of boundary condition in numerical simulations of vortex dynamics. J Fluid Mech 2004;516:115–24. [28] Sagaut P. Large eddy simulation for incompressible flows: an introduction. 2nd edn. Springer; 2004. [29] Thompson KW. Time dependent boundary conditions for hyperbolic system II. J Comput Phys 1987;89:439–61. [30] Verzicco R, Orlandi P. Direct simulations of the transitional regime of a circular jet. Phys Fluids 1994;6:751–9. [31] Wang FY, Zaman KBMQ. Aerodynamic of a jet in the vortex-wake of a wing. AIAA J 1999;37:1270–6.